Towhee-经验分享

以下是zhang_jaj 的问题:
你好，我想用towhee做NVT下argon气在surface面上nucleation过程。首先我想把potential选为LJ，surface与particles也是LJ作用，

1）问题是这个surface是需要用particles来make吗，如果是的话那么构成它的particles该怎么固定

2 ）或者，直接用towhee里面的LJ9-3 wall吗，想知道这种wall是与particles怎么作用的，它不是用小粒子构成的吗？

3）模拟完成后，液体与固体的结构判断该怎么做，有什么程序代码吗？

4）这个nucleation的rate该怎么算（是单位时间单位体积内产生临街晶核的个数），我想知道程序中怎么实现，我没发现有那个参数相关的，或者是需要自己来加东西吗？

5）还有这个nucleation barrier该怎么来实现？

6）towhee不能并行，那么其实一个任务就是在一个process上处理的；如果我的个人电脑比服务器的主频高，那么是不是用台式机比服务器更快呢（因为台式机是4核的，它是不是能自己把towhee的任务来用4核一起处理？）？

7）希望您能给点关于towhee代码学习的经验，就是对于我要做的这个方向的。

这是我的回答:
实在不敢当, 我现在也是学习阶段, 因为没做过类似的体系, 对于你的问题只能有一些建议.

1）问题是这个surface是需要用particles来make吗，如果是的话那么构成它的particles该怎么固定

可以使用particle模型来模拟, 固定比较简单, 在move设置中不让这些粒子移动就行了.

2 ）或者，直接用towhee里面的LJ9-3 wall吗，想知道这种wall是与particles怎么作用的，它不是用小粒子构成的吗？

这个我以前没有用过, 不过文献当中确实有人这样做. 这种势能的使用情况可以参考例子Fris_Walls, 另外上网搜了一下, 这几个帖子和文章可能对你有用.
http://jcp.aip.org/resource/1/jcpsa6/v67/i5/p2384_s1
http://jcp.aip.org/resource/1/jcpsa6/v73/i8/p4050_s1
http://lammps.sandia.gov/threads/msg24356.html
http://www.sklogwiki.org/SklogWi ... ard-Jones_potential
这种势是粒子势能转化而来的, 第四个参考文献有推导过程, 主要的优势是可以大大提高速度. lammps 的帖子跟你的情况有点像, 使用MD做的, 可以借鉴.

3）模拟完成后，液体与固体的结构判断该怎么做，有什么程序代码吗？

这个主要是分析RDF吧, 代码在 utils 中有, 可以根据你的情况进行修改.

4）这个nucleation的rate该怎么算（是单位时间单位体积内产生临街晶核的个数），我想知道程序中怎么实现，我没发现有那个参数相关的，或者是需要自己来加东西吗？

这个是需要后处理才能得到的吧. towhee本身不是计算成核过程的. 应该需要编写代码实现. 不过MC 如何计算时间呢? 如果用MC步数是不是会有问题? 是不是用MD更合适一点?

5）还有这个nucleation barrier该怎么来实现？

这个应该是进行不同温度点的模拟, 作速率-温度关系图, 然后通过拟合热力学公式得到吧.

6）towhee不能并行，那么其实一个任务就是在一个process上处理的；如果我的个人电脑比服务器的主频高，那么是不是用台式机比服务器更快呢（因为台式机是4核的，它是不是能自己把towhee的任务来用4核一起处理？）？

towhee的不能用多个核并行计算一个任务, 而事实上也不需要. 一般来说, 模拟多个任务可以等价于模拟一个大任务, MC的本质决定它是这样的性质, MD可不行. 所以尽量还是使用towhee提供的并行模式. 注意选择不同的随机种子.

7）希望您能给点关于towhee代码学习的经验，就是对于我要做的这个方向的。

关于代码, 其实我也没有太多经验. 我主要是对结果进行后处理. towhee 本身是一个general Monte Carlo code, 不是专门为某些体系设计的, 但是对于你的体系应该是可以胜任的. 主要针对体系的代码, 差别也主要是建模和分析结果两个过程, 核心代码都是类似的. 我个人认为如果时间有限, 不要陷到源代码中去.当然能够多了解一些代码对于软件的使用也是非常有好处的.
towhee 的代码由于需要我简单的读过一些, 还是比较清楚的. 可以尝试从总体上去学习下.
对于结果的后处理的代码, 建议从修改towhee 提供的代码开始. 使用语言以自己最熟悉的为准, 如果都凑合, 推荐 fortran90.
多学习 linux bash 脚本也是很有好处的, 我现在简单的数据处理基本都用它.
最后强调, 不要陷到代码中去, 能够完成自己的用途就可以了, 越简单越直接越好. 多多思考跟你体系相关的问题, 多查文献.

///希望能多多向你学习towhee知识，导师不管，自己一个人摸索。担心毕业十大问题啊，谢谢了

大家基本都差不多, 坚定信心, 可以做的很好. 其实模拟计算这个方向比起做实验还是比较透明的, 重复文献一般没有问题, 技术问题一般论坛上都会回答.
我也是在学习的过程, 大家多多讨论, 相互学习.

————————————————————————————————————————————————————————————————————

你好，想再请教下您。关于towhee参数设置问题：

1）.pm*的值是不是按照从小到大进行排列的呢。example中有个是

pmtraat 1.0 
pmtracm 0.67  
pmrotate 1.0
这种情况，当要判断是进行哪种move时，是不是只进行pmtraat呢？

2）.给体系加上一个LJ-9-3 wall时，里面有一个参数是ljfrho（the number density of atoms in the integrateed wall potential）是wall上与体系内分子作用的点密度吗？还有ljfsig与ljfeps都是人为的指定大小吗？非常感谢！

1. 是的，尝试概率逐渐递增。例如
pm1 0.1
pm2 0.5
pm3 1.0
pm4 1.0
那么pm1～4的概率分别为0.1，0.4，0.5，0. 只要出现1.0，那么后面的pm概率就都为0了
2.不是，参数的具体意义和设置看我发给你的几个链接。

————————————————————————————————————————

你好，还想向您请教个问题关于hmatrix参数的。mannul上说这个坐标矩阵可以不正交，但是必须是right hand rule的。我在想如果用其他非正交的坐标矩阵时，程序在运行过程中是不是速度就低了很多？？？（因为个人认为正交坐标里各个原子间进行距离计算时要方便的多（如我在一个面上固定一个单层石墨烯））

还有个问题是，对于NVT系宗，如果我先让N=1000，T=50K后，该怎么选取box的大小呢（气体），这时不同的box大小应该就对应不同的P吧。是利用PV=NRT吗？这个box的大小是依据我想要的P用那个公式来调节吗？，还是normal atmosphere下算个大约的box大小？非常感谢！！！

没有用过非正交的盒子，不过这个对速度影响应该不大吧。你可以试试。

对于气体的话可以使用PV=NRT来估算，注意单位。压力是根据你的体系要求指定的吧。

————————————————————————————————————————

还有关于LJ-9-3wall里面的ljfrho参数，这个我实在是理解不了是怎么算的，用哪些原子，除以多大的空间（和作用范围有关吗）？因为英语很菜，对那句话理解不清楚，非常感谢！

我看了一下LJ93wall的参数手册，lifrho是指每立方埃内的分子数。比如你的wall是Si，每个晶胞中分子数除以晶胞体积。如果你的wall是SiO2，可能就要建立位置重叠的两个wall，Si和O分别建立wall，这样的数密度分别是SiO2晶胞中Si的数目除以晶胞体积，O的数目除以晶胞体积。

——————————————————————————————————————————————

这个看了还是不太明白。按照mannul应该是把wall当做一种一般性的surface，而这个ljfrho就是这个面上的点分布密度。还有下面的那个sigma和epislon都是这个wall与体系中粒子的作用参数。（这个参数该怎么定义呢），还有如果把这个wall当做Si之类的来计算ljfrho，那么这些Si原子的排列该怎么确定(或者说这是一个无序的surface吗)？你之前给的那个文章链接我一直下不了，不知道你有吗，可以给我一份吗（zhjjaj@163.com）？关于LJ9-3的文章非常感谢！http://jcp.aip.org/resource/1/jc ... _s1?isAuthorized=no

————————————————————————————————————————————————

我从头讲好了。我们一般模拟时使用的是全粒子势能比如UFF，amber等等，这种势能每个原子作为独立的粒子，有自己的力场参数。这样在模拟时需要计算每一对粒子的相互作用能。如果我们的体系是气体或者液体分子和一个表面的相互作用，当然也可以使用全粒子势能参数去模拟，但是我们也可以把表面所有粒子的力场参数积分使其变成一个整体。9-3 wall 势能是从 12-6 粒子势能进行空间积分得到的。积分推导的过程见：http://www.sklogwiki.org/SklogWi ... ard-Jones_potential , 建议自己可以简单推导一下。这样做的好处是大大减少了粒子的数目（表面当成了一个具有特定形状的粒子），从而减少了计算量。推导过程中，原来12-6中的参数sigma ，epsilon意义没有变化，直接继承下来，另外新加入了一个参数ljfrhho，和体系的组成相关，即体相单位体积内的粒子数目。
我只是以Si和SiO2 为例，你要根据你的体系去得到sigma epsilon和ljfrhho。请问你的体系中wall的组成是什么物质？费了这么大劲，不知讲清楚了没有。那个文献我也没有，你可以求助一下

——————————————————————————————————————————————————————

感觉你讲了后，确实对这个势有了很多了解（之前特别迷糊）。我想用单层graphene和双层graphene试着模拟Ar的nucleation过程。这个2D的如单层的a为2.46（Angstrom），那我算的时候是不是应该就是2/（2.46*[sqrt（3）/2]*2.46）=0.38....

再比如一个fcc的（111）面，加入他的a=2，那么计算出来的是不是应该就是4/（2^3）=0.5？这样的话就跟选fcc的哪个切面无关了呀（或者我还是理解错了）？

还有我以为加上一个LJ9-3wall（如graphene）就是一个虚拟的wall；但如果我再在其他位置加一个同样组成的如graphene-wall时，并且给出这个graphene中各原子的坐标，而且这个wall不是LJ9-3的势。或者说在towhee中模拟的只要是加了一个用来吸附的表面，那么它就应该是当做LJ9-3wall来处理吗。这样他的原子排列还是不好确定（如模拟外延生长时要选不同的面）？
所以加入一个自己引入的面是不是应该和加一个LJ9-3是完全两个概念吧？

还有，以你的经验，你觉得像graphene有近1000作用的碳原子，这个得当做一个独立的大分子来建立初始文件，你觉得应该用哪种方式建呢，这样就得在input_style下面有1000个作用的vibration手动参数输入，感觉太恐怖了。

最近在这几个问题上实在搞不定，就只能想到麻烦你了

——————————————————————————————————————

还有一个问题是，我用ms建了graphene的模型，导出pdb文件格式，但是他是以ATOM位头的，pdb2towhee转换好像只能是HETATM打头的，这个格式转换是个问题啊。是不是应该自己编一个转换的脚本呢。

还有就是如Ar气在graphene上nucleation的问题，我是把Ar放在towhee_coords里，然后把graphene放在towhee_tempalte里面，这个石墨烯用的是trappe-ua的力场，问题是每个碳原子都要输入一个vibration成键的参数，感觉工作量特别大。不知道石墨烯有没有比较好的输入方式。谢谢！

单层和双层石墨烯使用9-3 Wall 势能模拟可能不太合适，因为积分推导过程中是认为表面厚度无限大。1000个C原子使用全原子势能也是可以接受的吧。
你理解密度的计算还是不太正确。9-3 wall 积分时认为表面厚度无限大，那么它的密度就等于体相的密度，跟切的表面没有关系。
没有明白你第二段的意思。你能再简单地讲一下吗？towhee中大分子的建模确实比较头疼。不过，如果你的graphene在模拟过程中是固定不动的，那就简单了。你把1000个原子的石墨烯当作1000个独立的分子，这样就不需要输入vibration 参数了。

pdb 转化有问题可以修改pdb2towhee.f90 ，也可以编脚本实现

——————————————————————————————————

非常感谢啊，我今天下午去把这个大分子上的每个原子当做独立的分子，程序成功运行了。我想了下，其实每个c原子在力场都是一个独立体对待的，所以固定不动时把它当做很多歌分子也是可行的（下午给我乐的，呵呵）。其实这个LJ-9-3我发现它就是一个弥散的surface吧，我那段话说的意思主要是为了确定这个surface表面的原子排布不好确定，现在我想通了，它应该没有这个原子排列的概念吧。希望你能给我举个例子，比如我分别用bcc与fcc的铁的来做这个LJ-9-3 wall，加入他们的晶格常数都是2（Angtrosm），那么这个ljfrho分别是怎么求的，值时多少？而ljfrho下面那个ljgsig与ljfeps的值，是不是需要重新建个模型用Fe的势与体系的势进行一个cycle，然后从towhee-out中看Fe与体系中各原子的sigma和epsilon值？

还有towhee mannul中关于lj-9-3介绍下面那个公式里的sigma与epsilon应该就是ljfntypes下面的一系列ljfsig与ljfeps值吧，这样的话，如果体系有n种与surface作用的原子，那么在积分过程中就应该有n中lj-9-3的积分式（参数不同）吧

关于那个pdb2towhee的我也搞定了，用":1, $s/ATOM  /HETATM  /g"，然后就可以了。

有进展值得庆贺。比如bcc的铁，每个晶胞的体积为2^3=8, 包含2个Fe原子，那么它的ljrho 应该是2/8=0.25。ljgsig 和ljfeps 你说的是对的。学习软件最主要的还是靠手册学习，towhee的手册写得还是挺详细的，我当时也是一条条自己看的，然后结合例子进行学习。

——————————————————————————————————————————————————

关于涉及到石墨烯的建模，因为石墨烯是固定的，所以我按照你说的那种方法把每个原子当做一个独立的分子来建模。这样存在着问题：
1.由于突然分子数多了所以对于inix， iniy， iniz 那三个值就需要增加，但是不知道这三个值到底是什么意思，把盒子划分开有意义吗（因为原子坐标都是从towhee_coords中读取的）？

2.当我用c60试验的时候，这种方法一切正常，但是当我用涉及到1000多碳原子的石墨烯用这方法时，运行什么的都正常，但是在输出文件中有问题：信息只输出到关于各原子间势那块就没了，也就是显示了LJ公式，然后是sigma   epsilon这一列，下面没东西了。。。不知道怎么回事，但是程序确实还在运行中，而且还不停有box_*.pdb文件产生，好生纳闷啊

希望你能帮我排除下故障，谢谢啊

我下午去试了下，好像是因为没有算完，所以到那里就不显示了。。。

——————————————————————————————————————————————

楼主你好，我也是刚刚接触到towhee软件的，我是做三元萃取体系气液平衡的，现在需要模拟纯组分噻吩(thiophene)的气液平衡情况，在模拟之前看了一下要用的力场文件，发现里面缺少一部分我所需要的键长键角文件，请问那个力场文件是不是没有显示全部可以处理的键长键角数据呢？（因为之前读过一篇文献说towhee钟的Trappe-EH力场可以处理N原子的问题，但是我进入到towhee 7.0.2版本的那个力场文件后发现里面根本没有N原子的数据。所以才有这种疑问##）。如果真的没有我该自己添加力场的数据么？需要的话是简单地从别的联合原子力场中找到N原子数据，然后粘贴我现在用的这个力场中吗？小弟初学，麻烦你了！

对这个力场我不是很熟悉，最好的办法是找文献中Trappe的参数，然后按照towhee的格式添加到力场文件中。如果Trappe力场确实没有你所需要的参数，可以从别的力场得到参数，一定要两者的力场形式相同（比如都是LJ 12-6）。

————————————————————————————————————————————————————

楼主师兄你好，我又有几个问题想问一下你。
1.两个力场的混合规则不同，可不可以用什么方法给其中一个导入一种混合规则呢？如果可以怎么做？
2.对于手动输入电荷时，如果关于力场的描述文件中的Coulombic interactions项下有电荷值，我是不是就该输入此电荷值，如果Coulombic interactions没有值，我是不是就不输入电荷值呢，即为0.0d0？入Dubb2004势，用来建立SiO2分子时？
3,在建立一个分子时，如果有多个units的atom组成，那么对于towhee_coords文件里面是不是应该按照input_style下的unit值对应的把一种atom输入完再输入另一种呢？
4，您知道这个程序的结果文件pdb或者movie文件里怎么判断材料是fcc还是bcc之类的吗，或者您知道那个软件有这个功能吗？
5.（补充的）对于不同的力场文件，如argon的epsilon值都不相同，是不是意味着argon在固体情况下，对于不同的力场他的晶格常数不同呢？
6，对于涉及同位素的，就拿H和D来说吧，如果我又H的力场文件，可不可以把H的mass修改一下当做D的立场文件来用呢？（对于LJ势）
非常感谢！

1. 混合规则不同一般不推荐联合使用，可以测试一下看看影响是否很明显。
2. 不考虑库仑作用可以把电荷写为0，如果全部都不考虑可以把库仑作用关闭
3. 不明白？？举个例子？
4. 这个需要你自己去判断吧？我没有发现有这样的功能。
5. 应该是这样的，注意也比较一下力场的形式
6. 可以的。

————————————————————————————————————

你好：
1.我不知道混合规则不同的能不能用，所以就没用，下午可以去试一试
2.也就是说，只要我的体系不想考了库伦作用，我就可以都不写电荷值吧
3.问题解决，因为当时我导出来的sio2的pdb文件里si和o是间隔的写的，对于towhee需要把si的坐标写完，再写o的坐标。
4.因为我做的是气体低温结晶过程，所以需要判断到底是冻结成了什么结构，这个很郁闷啊
6.关于这个我知道对于vasp是可以这么修改的，但是对于LJ势的力场文件，因为不知道参数和原子的质量有没有关系，所以不知道这么修改得到的同位素力场可不可以用。

问题：还想问一下，不知道您对形核过程了解吗，我不知道该怎么在模拟过程中计算体系的形核势垒。如果说是用成核时的体系总能量与初始的总能量相减的话，那么这个势垒是不是特别依赖与初始模型的原子排布（因为我是看您的教程用packmol生成的文件，所以不知道这个生成的排列对于计算势垒时是否具有一般性）
towhee里面有umbrella sampling和flux forwar sampling两种抽样吗？如果没有想要加的话是在编译后的程序里加呢，还是加了之后重新编译呢？
非常感谢！！！！

6. 对于同位素：量子力学都可以使用相同的参数，分子力学当然也可以了对于成核不太了解。形核势垒可以通过多次模拟把初始结构的影响降低吧。umbrella sampling 在towhee 中是有的，后面那个不太清楚，你可以看一下手册。如果要修改程序，当然需要重新编译了！

——————————————————————————————————————

如果可以直接把力场文件修改当做同位素来用的话，那假如我需要涉及H2和D2两种气体的混合，那我是不是需要把复制一个H2的力场修改后重命名或者放到其它目录来调用呢？

可以这么做，也可以在原来的力场中加入D原子类型

——————————————————————————————————————————————————————

你好，关于同位素力场文件，改了mass这个值就能用我有点疑问。如果说关于H和D的，质量不同，然后其他什么都相同（力场），那么他们两个单独进去的模拟应该表现的情况一样是吗？因为mass好像对LJ中的势能没有影响吧。这样的话，其实同位素就不可分了吧？

LJ 势能确实跟质量没有关系。不过模拟的整个过程是这个样子的：
初始速度和坐标---> 根据势能函数求作用力 --> 根据牛顿第二定律求加速度（！！！，跟质量有关） --> 由加速度求速度 ---> 由速度求位置 --> ......循环

————————————————————————————————————————————————————————————

这个是分子动力学的吧啊，如果是蒙卡的话，那就没法区分H2与D2了吗，是吧？

不好意思弄混了。我想应该是没有办法区分的

——————————————————————————————————————————————

那就傻了，难道只能换lammps了吗？我在towhee里面的力场没有看到有D2和T2的。

不知你研究的是同位素的动力学效应还是热力学效应，如果是动力学效应，MC 本身是不能做的。如果是热力学效应，说实话，我目前还理解不了它们之间有什么差别

——————————————————————————————————————————————————————————————————————————

我是想做H2和D2在低温结晶过程的模拟，看他们的nucleation过程等，但是据我所知，他们的三相点是明显不同的，所以气体液化的温度也不同。但是用同一个力场的话感觉就是全同的了。。。体现不出什么差别来？

——————————————————————————————————————————————————————————————————————————

你好，看来只能又麻烦你了啊。之前你给我说towhee里面有umbrella sampling这个算法，但是我不会用它来确定我的nucleation过程的势垒，不知道这个具体该怎么做，麻烦能给我说个大概的过程吗？非常感谢，还有那个势函数我不用去纠结了，目前先不让做那个了，只做氩气的。

——————————————————————————————————————————————————————————————————————————

你好，又要麻烦您了！之前和您讨论过建立石墨烯的模型，当时是把所有的C原子当做独立的molecular来对待，即势函数只考虑了非键参数，没有考虑键参数。这个的前提是让石墨烯不参与move（这样的话把C-C作用搞成0更好，还不知道该怎么来做，难道要单独由几何法得到C-Ar参数做一个力场？）。但是最近一个师兄说我的模型太理想了，应该让C原子也move，如果这样的话，就需要考虑石墨烯中的非键作用。如果这样的话不知道这个石墨烯的模型该怎么建？
如果把所有的键都算上的话，在towhee_input下面得输太多太多行了。。。关于这个问题，学长，你怎么看？

———————————————————————————————————————————————————————————————————————————————

 

首先在这里我要衷心的感谢043114076师兄，在学习towhee过程中他给予了我极大的帮助和指导。他不但具有丰富的学识，而且非常友善有耐心。

关于input 
首先确保格式正确无误，然后参数设置方面，参数都是有物理意义的，重要的是对monte carlo方法本身的理解，不要局限于towhee本身的manual（不是说不重要）。当然最初接触要多加调试，会加深理解，最后总结的时候要与参数物理意义结合，归结于理论。

参数的设置对于模拟时间和结果的准确性影响很大：

模拟参数的设置，如controlstyle只对输出参数的形式有影响，可通过自己想要怎么样处理数据进行设置。equili 和product过程的差别主要在于是否去做结果分析。equili 过程是不能做分析的，主要是为了得到热力学平衡的结构；product 的结果才能做RDF 、RMSD等等的分析。仔细比较一下不同controlstyle的差别，就会发现两者只是生成output的差别。一般是采用manual，按照自己想要的输出结果进行设置。

如stepstyle的选择
'cycles' 跑nstep步mc循环。每次循环包含分子数次mc尝试移动，所以速度慢点。在一个cycles中每个移动的概率是不一样的。
'moves' 跑nstep次mc尝试移动。（速度非常快，每次概率一样）

要充分理解两者的本质区别，从而进行取舍。

初始构型的设置主要是影响平衡时间，对最终结果影响不大，但能够快速达到平衡对计算来说非常具有优势。若没有通过其他途径得到初始构型，那就直接从input里设置；自己设置时一定要注意hmatrix，盒子大小的设置，对于平衡影响很大，一般液相盒子通过分子数和液相密度计算得到，气相盒子通过理想气体状态方程计算得到。

若有从MC或MD方法得到的初始构型可以转化为initial或coords文件进行计算，从而能够快速进行计算。

使用initial时：
linit(logical)       
        .TRUE. 生成初始构型。
        .FALSE.不生成初始构型，从towhee_initial中读入。

计算过程停电终止，若想接着前面的计算继续下去，按照帮助文档上说是把backup里的内容复制到initial里，然后再运行就可以。但尝试一次直接计算100000的，再尝试一次先计算50000再利用initial计算50000，这两次运行的结果并不太一样，最后的结构可能是不一样的，但是热力学信息应该是类似的。要么可能是运行的步数不够多。因为monte carlo是随机过程，每次随机数都是不同的。

使用coords时：
可以把pdb 文件转换为coords 文件，utils 中有这个小程序。 
coords文件的格式比较简单，手册中的说明也比较清楚。注意在input中初始结构的参数（initstyle 等）也要相应的修改。
模拟移动类型比率设定，概率为累加概率，即实际概率为当前数值减前一项数值，必有一项概率为1.后续项则均为0.
pmvol(double precision) 盒子体积变化，对于非常消耗模拟时间，比例通常设置非常小，而且在NVT系综时此项不起作用，可设为0.
pm2boxrbswap
pm2boxcbswap
这两项都是粒子交换，根据自己所计算体系进行选择，在高温时候可适当明显减小比例会大大优化模拟结果。
各种移动的接受概率不能太高也不能太低，否则效率都是不高的。要多试试，不同的体系是需要经验的。有文献说0.3-0.5之间模拟效果最好。
电荷使用手动指定电荷的方法（manual），在input file 指定电荷。这样就可以使用文献当中的电荷值。电荷的分布对于模拟结果影响很大，一般不设置的话都是默认采用软件自带的从烷烃和烯烃等体系迁移过来的电荷值。若要修正电荷，可通过量化方法进行拟合，或者与实验数据拟合进行不断调整。
例如我运行中出现运行towhee ,出现assemble：unknown match_style for nonbond，仔细核查了输入文件都没有发现问题，总以为是力场里非键参数的问题，结果只是格式出错了，解决办法如下：

在linux下运行：
dos2unix towhee_input
dos2unix towhee_ff_TraPPE-UAn
之后我试了这个方法，没有转化成功，不知道是什么原因，但是师兄电脑运行很好。然后师兄又给我一个建议：如果没有dos2unix 这个命令，可以使用linux的vi 打开你的上述两个文件，把所有^M 去掉。
果然用vi打开input 文件以后，发现了很多^M，删除以后运行正常。
这是属于格式问题，原因是在windows 编辑的文件与linux不兼容，有两个解决的办法：
1. 一般情况下，使用dos2unix 命令是可以进行转换的。
2. 即使不能转换也没有关系，可以在linux 下使用vi 打开你的文件，把那些不应该出现的特殊的符号删除掉，最常见的是"^M"
另外编辑towhee_input 文件，推荐使用unicode类型的编辑器，我使用的是geany, 一个免费开源的软件，win和linux都可以使用。
每次编写完input，检查完各参数的设置后一定要检查格式，一个小小的标点符合，一个空格都会使程序运行不下去。
运行过程中，一旦出现问题就会中断并提示错误，仔细琢磨提示的错误，再进行修改，如我遇到过：出现 'command not found',有时候是因为命令打错了，大小写没有注意，有时候是由于命令与输入文件不在同一个路径内；出现'已杀死'，有时候是因为终端运行任务太多，超过电脑负荷，有时候是因为运行程序所需要的输出文件不完全，或者终端写代码时候与input文件不负荷；出现 'expected string',有时候是因为输入文件里的参数漏了某些前后对应的，有时候是因为格式问题。。。。。。

气液平衡确实比较难以模拟，跟目前理论方法有关，加深对monte carlo方法的理解会有很大的帮助，也会对选择参数有指导。这些参数都是有物理意义的，重要的是对方法本身的理解，不要局限于towhee本身的manual（不是说不重要）。多多调试也会加深理解，使误差缩小，但要明白这些参数的调整只是为了更快达到平衡或者更好地收敛并不改变方法的本质。
关于output文件
计算结果output最后的block averages,里的 block energy, density, virial press, mol fracs,这四项有时候会突然全部变成 0.  例如最后输出10个block averages,从第5开始变成0 ，只有一直如此；还有时候是从第五个变成0到第七个，然后第八个又恢复正常值。
为什么会出现这种不正常的现象呢？
后来我仔细检查了初始文件、结构、中间的结构、output 等等，并和例子文件对照，然后调整了盒子大小，改变一下粒子交换的比例，目标接收概率等解决了这个问题。
在系统温度接近临界温度时，模拟体系的波动显著增强，容易出现某模拟盒的粒子数为零的情况，导致模拟进程中断。因此，必须小心调整模拟条件。高温（接近临界温度）时候，将粒子交换步所占的比率显著减小，能够获得较好的模拟效果。

关于力场
若计算体系里分子片段在软件里都有定义，那就直接采用自带力场；若有些分子片段没有定义，要么从其他分子片段里迁移，要么就自己按照所要使用力场的格式要求编写。
力场修正，有(1) 完全拟合实验数据。(2) 部分源于量化计算，部分源于对实验数据的拟合。(3) 完全拟合量化方法计算得到的数据。
一般的分子内相互作用（如键长、二面角等）可完全通过量化计算获得准确的结果，但是分子间的相互作用（如库伦相互作用、范德华相互作用等）则很难完全通过量化计算得到准确通用的结果，需要不断拟合，最好是参考实验数据或者状态方程计算的数据。要注意无论是经验的还是量子力学方法的。这些参数之间是compromise的结果，所以自己调整时要尽量小心。除了参数，模拟体系的模型选择也很重要，好的模型可能使得更快达到平衡。

关于Utils
towhee 的utils 目录下有一些非常有用的小程序，也可以写一些简单的脚本。如pdb2towhee可以将pdb文件转化成其他形式的文件从而在输入文件中使用。Fitcoex可以得到临界温度和临界密度。rdf2pmfair可以得到RDF图（radial distribution funtion)等。
举一例：
我按照手册上的步骤
先    cd/Utils
      make pdb2towhee
然后我也发现Utils里有了这个可执行文件。
接下又到带有pdb文件的计算结果文件夹里执行这个   pdb2towhee命令
可是总是说 command not found。
后来我把可执行命令pdb2towhee复制到结果文件夹里就可以执行
./ pdb2towhee 这个命令了，按照提示一步一步输入代码。都是一些简单的小问题。
出现这个问题是由于执行命令和结果文件不在同一个路径。也可以尝试把执行命令复制到bin文件中。


关于模拟方法的改进以及与其他方法结合

气液相平衡的时候先用MC或者MD得到平衡的气相和液相，再用GEMC。

如果用MC的话，towhee就可以，关键是运行时只设置一个盒子，分别得到气相和液相的平衡结构。control style 设置我一般是使用manual 自己指定每个参数。（这个关系不大，只是输出的参数有差别）观察能量和压力等参数达到平衡以后，再运行两个盒子的任务。一般气相平衡很容易，液相要多运行很多步。平衡的判断最重要的是查看MSD。如果会MD的话，也可以试试，平衡的速度更快。

结语
气液相平衡的模拟采用的蒙特卡洛方法，首先要理解这个方法的原理，才能加深对于软件学习以及参数设置的理解，从而可以避免过多无谓的调试，尽快地达到平衡进行计算缩小误差。学习towhee 首先是input 文件，参照manual 理解每一个参数的意义，各参数之间是compromise的关系，注意前后对应。能够快速达到平衡，减小误差

然后是力场，一般体系都是使用软件自带力场，有些时候只需要修改一些参数就行了。但是遇到软件没有定义的分子基团，就要自己进行编写力场了，格式严格按照手册上的要求，注意对齐，不要有多余的空格。目前我刚接触towhee一个月，很多东西仍然在学习，把自己的学习经历贴出来跟大家分享，希望得到大家的指导。最后还要再次感谢043114076师兄