以下是 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 吗 ) ？你之前给的那个文章链接我一直下不了，不知道你有吗，可以给我一份吗（ [email protected] ）？关于 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 师兄