关于towhee 回答zhang_jaj的问题-转载小木虫

Posted 2020-11-01 Simyang

以下是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下面得输太多太多行了。。。
关于这个问题，学长，你怎么看？

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