君，已阅读到文档的结尾了呢~~
扫扫二维码，随身浏览文档
手机或平板扫扫即可继续访问
gromacs文件介绍and一些杂知识
举报该文档为侵权文档。
举报该文档含有违规或不良信息。
反馈该文档无法正常浏览。
举报该文档为重复文档。
推荐理由：
将文档分享至：
分享完整地址
文档地址：
粘贴到BBS或博客
flash地址：
支持嵌入FLASH地址的网站使用
html代码：
&embed src='/DocinViewer-4.swf' width='100%' height='600' type=application/x-shockwave-flash ALLOWFULLSCREEN='true' ALLOWSCRIPTACCESS='always'&&/embed&
450px*300px480px*400px650px*490px
支持嵌入HTML代码的网站使用
您的内容已经提交成功
您所提交的内容需要审核后才能发布，请您等待！
3秒自动关闭窗口查看: 11594|回复: 87
使用molclus程序做团簇构型搜索和分子构象搜索
在线时间1584 小时
帖子威望eV
重要补充说明1：molclus从1.2版开始加入了genmer工具用来产生团簇的初始构型，详见《genmer：生成团簇初始构型的超便捷工具》（），比起本文例子中通过分子动力学程序产生分子团簇初猜结构方便得多得多。
重要补充说明2：molclus从1.3版开始加入了gentor工具用来通过扫描方式产生分子构象搜索用的初始构型，详见《gentor：扫描方式做分子构象搜索的便捷工具》（），因此对小分子做构象搜索时就不是必须得像本文那样通过动力学来产生初始构象了，省事不少。
使用molclus程序做团簇构型搜索和分子构象搜索
文/Sobereva(3)
First release：2015-Jan-5& &Last update：2016-Jun-18
前言：molclus是北京科音自然科学研究中心开发的一款用于分子团簇构型搜索和分子构象搜索的程序，需要结合第三方分子动力学和量子化学程序使用，本文将对此程序进行介绍。第1节介绍基于分子动力学做构型/构象搜索的基本原理，第2节介绍molclus程序的基本使用，3、4节通过实例进行演示。
1 构型/构象搜索原理
分子团簇的构型搜索和分子的构象搜索是计算化学经常涉及的问题。前者用来寻找分子团簇的各种构型，比如可以回答这样的问题：5个某分子组成的团簇，总共有多少种结构？能量各是多少，哪个构型最稳定？后者用来寻找特定分子的各种可及的构象，比如一个柔性分子，有很多键都可以旋转从而构成不同构象，我们想知道哪个构象能量最低，各种构象能量差是多少，由此我们还可以进一步根据Boltzmann分布估计出特定温度下不同构象出现的比率是多少。
构型/构象搜索的算法很多，比如势能面扫描、分子动力学、蒙特卡罗、基因算法、LMOD等等，适用场合不同。对于分子团簇构型搜索和分子构象搜索而言，分子动力学方法是比较好的选择。其基本过程是：在合适的方法下做一定长度的分子动力学模拟，由于热运动，结构能够反复穿越势垒到达势能面不同区域。每隔一定步数取一次结构，将这些分布在势能面的不同位置上的结构作为初猜结构分别进行优化，最终就会得到势能面上各个极小点位置的构型和能量。
对于分子团簇的构型搜索和分子的构象搜索而言，并不会涉及到键的形成和断裂，因此不需要用极耗时的从头算分子动力学，而只需要计算量极低的分子力场来跑分子动力学就够了。这样的程序很多，比如免费的gromacs、tinker、NAMD、lammps，现在amber实际上也算基本免费了（用免费的ambertools里的sander来跑就行了，不需要用收费的PMEMD部分），以及收费的Material Studio。至于结构优化和能量计算部分，可以用分子力场、半经验、DFT、高精度后HF方法等来做，这就取决于需要的精度和运算能力了。
2 molclus程序
molclus有Windows、Linux和MacOSX版，由Fortran编写。molclus主页为：。程序会不断改进、更新。molclus的免费版可处理的原子数目上限为30个。如果需要处理30个以上原子，需购买完整版，价格96元。源代码有偿提供，价格399元。付费购买2年内可免费升级，并提供技术支持，免费版不提供技术支持。购买和下载事宜请访问molclus的主页。
molclus程序可以视为是分子动力学程序与量子化学程序之间的一个接口。几乎所有动力学程序都可以用，包括自编的，只要能输出.xyz格式的轨迹或者输出的轨迹能被免费的VMD程序（）识别即可。量化程序目前支持Gaussian、ORCA、MOPAC2012。其中Gaussian是大家最常用的，主要适合小体系。ORCA由于RI技术，在纯泛函下比Gaussian快得多得多，因此可以用于中等体系，当然对于小体系也没问题。MOPAC2012专门做半经验计算，它支持的PM6-DH+通常足矣得到定性合理的结果，而速度比DFT快两个数量级，故很适合用于中、大体系的优化。可以说，molclus支持的Gaussian、ORCA和MOPAC2012把不同尺度、不同精度范围都覆盖了。如果不熟悉ORCA、RI技术和MOPAC2012，请务必参看《大体系弱相互作用计算的解决之道》（）当中的介绍。分子团簇和弱相互作用关系极为密切，必须选择合适的计算级别，如果不熟悉弱相互作用计算的话建议看看《乱谈DFT-D》（）。
2.2 使用流程
molclus直接解压后就能用。基本使用流程如下：
(1)用分子动力学程序模拟分子团簇或单个分子，每隔一定步数将结构写入一次轨迹。力场并不要求很准确。
(2)用免费的VMD程序将动力学程序输出的私有轨迹格式转化为标准的.xyz格式轨迹文件，命名为traj.xyz。
(3)设定好量化程序的几何优化模板文件，模板文件需放到当前目录下，以template命名。对于Gaussian程序，后缀为.gjf；对于ORCA，后缀为.inp；对于MOPAC2012，后缀为.mop。molclus会将模板文件里面的[GEOMETRY]字段替换为traj.xyz里的结构并传递给量化程序执行，因此应该在模板文件设定好基组、方法和相关关键词（比如内存分配量、辅助收敛的设定等）。
(4)设定好molclus的参数文件settings.ini，详见2.3节。
(5)启动molclus。molclus会载入traj.xyz、settings.ini和模板文件，然后自动调用量子化学程序对traj.xyz的每一帧结构进行优化，并将优化好坐标和能量写入到当前目录下的isomers.xyz文件中。
(6)启动molclus中附带的isostat工具，它会载入isomers.xyz，对其中的结构根据能量进行排序并进行统计、去除能量和几何结构相似度过高的重复结构。isostat在当前目录下输出的cluster.xyz文件就是根据能量从低到高排列并且去除了重复结构后的文件，可以载入VMD来方便地观看每个结构。之后还可以输出指定温度下各个构型的波尔兹曼分布比例（原理见，理应用自由能进行计算，isostat用的是单点能作为其近似）。
如果想使用比结构优化更高级别的方法精确计算每个结构的能量，可以把settings.ini里的ieneonly改为1，把模板文件里的和优化有关的关键词都去掉，再把cluster.xyz改名为traj.xyz然后启动molclus，molclus就会调用量化程序对其中每个结构再做一次单点计算。
每次程序运行前，会对当前目录下.out、.arc等量化程序之前产生的临时文件进行清理，所以如果有重要的临时文件在当前目录下，执行molclus前应当先将它们挪至它处。
2.3 参数文件
程序目录下的settings.ini在molclus启动时会被载入，其中记录了运行参数，在下面进行介绍。大部分参数都不用改。
isys：1=Windows；2=Linux。在什么系统中运行molclus就设为什么
iprog：调用的量化程序。1=Gaussian；2=MOPAC2012；3=ORCA
ngeom：0=处理traj.xyz里每一帧；n=处理前n帧；n,m=处理n~m帧（比如2,5即处理2，3，4，5帧）
ieneonly：0=对traj.xyz里每一帧进行优化；1=对每一帧只计算单点能
distmax：n=如果结构中任意两个原子间的距离超过n埃则跳过这一帧。做团簇动力学的时候，有可能某个分子跑飞了，对这样的结构优化没有意义，利用这个参数就可以忽略掉此结构
ipause：1=对于Gaussian和ORCA，计算出错（如不收敛）则暂停，以便用户对输出文件进行检查，然后可以按回车继续处理下一帧结构；2=每处理完一个结构后都暂停；0=总是不暂停
ibkout：1=每处理完一帧结构，把程序的输出文件备份，以便之后检查或根据自己的需要提取其它数据。比如对于第5帧结构，Gaussian的输出文件会备份到gau00005.out，ORCA的会备份到orca00005.out，MOPAC的会备份到MOPAC00005.arc；0=不做备份，因此molclus运行完毕后将找不到量化程序中间输出的文件
----以下是Gaussian专用的参数----
gaussian_path：Gaussian的可执行文件的路径，比如&/sob/g09/g09&、&d:\study\g09w\g09.exe&。也可以只写&g09&，前提是程序所在目录已经加入到了PATH环境变量里。如果路径中有空格，则必须有双引号，没空格的话无所谓，后同
igaucontinue：1=如果优化是因为超过了步数上限而引起的，则会基于另一个模板文件template2.gjf对最后一帧结构进行继续优化；0=不做这个处理。igaucontinue这个选项主要有两个主要用处。当初始结构偏离平衡结构太远时，无论Gaussian输入文件里opt=maxcyc=n设多大，实际步数都不允许超过程序内定上限（和坐标变量数成正比），因此可能因为步数上限不够而导致没收敛，因此可以基于template2.gjf继续优化。另一种情况是基于template.gjf里的关键词在优化时发生了震荡而没有收敛，若在template2.gjf里写上备用方案的关键词，诸如calcfc、GDIIS、opt(maxstep=x,notrust)等等，可能在这些关键词下继续优化的时候就能收敛了。关于辅助收敛的方法，在此贴有汇总：《量子化学计算中帮助几何优化收敛的常用方法》（）
energyterm：molclus是从Gaussian输出文件中末尾的archive字段部分读取能量的，这个选项用于设定读取哪项。比如是HF/DFT计算，这个参数应设为HF=，即读取HF=这个字符串后面紧跟着的那个数（HF/DFT能量）。如果是MP2或双杂化泛函计算，就应当设MP2=。以此类推
ibkchk：和ibkout的用法一样，但这个控制的是如何备份.chk文件
----以下是ORCA专用的参数----
orca_path：ORCA可执行文件的路径，比如&D:\study\orca\orca.exe&
ibktrj：1=把每次优化后产生的.trj改名为orca[帧号].xyz。这是优化过程的轨迹文件，因此可以用VMD直接观看优化过程
----以下是MOPAC2012专用的参数----
mopac_path：MOPAC2012可执行文件的路径，比如&C:\Program Files\MOPAC\MOPAC2012.exe&。
3 分子团簇构型搜索实例：(H2O)4
此例我们通过molclus结合Gromacs 4.6.7分子动力学程序以及MOPAC2012、Gaussian09 D.01量子化学程序搜索(H2O)4团簇构型，这个团簇已经被很多文章研究过了。这里用这种很小的团簇进行示例仅是为了节约时间，对无论多大的团簇计算过程基本都是一样的。
Gromacs是免费、最主流、速度快且很好用的小分子和生物大分子的分子动力学模拟程序，所以我们这里用Gromacs，编译方法参见《Gromacs 5.0与4.6.7编译方法》（）。大家也可以用其它自己擅长的动力学程序，只要能产生VMD能支持的轨迹就行，因为我们要靠VMD来把轨迹转为molclus能读取的.xyz轨迹格式（当然啦，如果自己有其它办法能转成.xyz轨迹格式也行）。
动力学我们用NVT系综。对于当前体系实际上不用周期边界条件也行。不过对于较大团簇的话，当模拟温度较高时，可能模拟过程中外侧的某个分子跑着跑着就飞了，越飞越远，之后的轨迹就没意义了。而如果用了周期边界条件，它飞过边界后还会跑回来。
Gromacs模拟部分所有涉及到的文件都可以在这里下载，如果有的地方不清楚直接看里面的文件就明白了。
(55.96 KB, 下载次数: 178)
18:18 上传
点击文件名下载附件
3.1 准备模拟体系
H2O我们就用SPCE水模型，它的拓扑文件是Gromacs自带的，就是gromacs/top/gromos54a7.ff/spce.itp。我们用gview之类工具画一个水分子，保存成H2O.pdb。由于spce.itp里氧原子名字为OW，两个氢原子名分别为OW1和OW2，所以我们也把pdb里的原子名改为与之对应的。
然后建立个文件夹用于存放模拟相关的文件。把H2O.pdb放进去，并建立一个空盒子文件emptybox.gro，盒子各边长都为1nm，内容如下：
emptybox
0
& &1.0 1.0 1.0复制代码
运行此命令往这个小盒子里加入4个H2O
genbox -cp emptybox.gro -ci H2O.pdb -o H2O_box.gro -nmol 4
editconf -f H2O_box.gro -o H2O_box2.gro -d 2
此命令重新建立了盒子，是围绕这四个水在各个方向延展2nm得到的，这使得这4个水处在这个大盒子中间区域，而且盒子尺寸也大于非键作用cutoff距离的两倍。（也有其它建立初始结构的方法，比如用packmol生成初始结构，可控性更高）
编写当前体系的拓扑文件system.top，内容如下
#include &gromos54a7.ff/forcefield.itp&
#include &gromos54a7.ff/spce.itp&
[ system ]
4*H2O
[ molecules ]
SOL 4复制代码之所以[ molecules ]部分这么写，是因为spce.itp当中水的分子名是SOL。
3.2 动力学模拟
我们用1fs步长，模拟1000000步（1ns），每隔10000步（10ps）保存一次坐标，共会得到101帧。保存间隔太大则需要更多的模拟时间，而太短的话则相邻两帧结构往往差异较小，优化后会容易收敛到相同的结构。对于这种小团簇来说，动力学模拟的耗时是可以忽略不计的，所以间隔可以设大些，但太大也没意义。关键是获得的帧数要足够多。
特别需要注意的一点是，像水团簇这样以静电作用为主导的团簇构型（氢键的本质主要是静电作用），如果直接模拟的话，由于氢键作用太强，团簇结构很稳固，总是在最稳定的那几个结构中变化，难以跑到高能结构上。虽说提高模拟温度可以让不同结构之间变化得更容易，但是实际做过就会发现，温度设很高后虽然团簇没这么牢固了，结构变化更自由了，却很容易跑散，最终还是得不到能量较高的构型。因此，我们要人为把静电作用削弱，这里我们用epsilon-r=4关键词来把相对介电常数由默认下真空的1改为4，即静电相互作用被削弱到了1/4。此时，氢键维持的高度有序的结构就容易被热运动破坏了，缺乏取向性的范德华作用对于保持团簇状态的重要性也大为提升了。
由于氢键被我们削弱了，范德华作用的强度又远远低于氢键，因此我们模拟所用的温度也必须比一般化学直觉要低得多，不然肯定跑散了。应该对模拟温度进行几次测试，摸索出尽可能高但又保证团簇不散架的温度。对于此例，我们发现60K比较理想。注意不宜一上来就控温到60K，否则由于初始动能太大靠前的一些帧容易散架。因此我们需要利用gromacs的退火设定，设起始温度为0K，在前60ps内线性升温至60K（即1ps升温1K），然后一直保持在60K。
非键相互作用用最简单直接的cut-off方式计算就行了，1.4nm截断距离完全绰绰有余了。
以下是我们模拟用的参数文件md.mdp
integrator = md
dt = 0.001 ; ps
nsteps = 1000000
nstxtcout = 10000
;
pbc = xyz
rlist = 1.4
coulombtype = cut-off
rcoulomb = 1.4
vdwtype = cut-off
rvdw = 1.4
;
epsilon-r=4
Pcoupl = no
Tcoupl = v-rescale
tau_t = 0.5
ref_t = 100
tc_grps = system
annealing = single
annealing_npoints = 2
annealing_time = 0 60
annealing_temp = 0 60
constraints = h-bonds复制代码（ref_t这个参数此时无意义，因为退火设定已经设了温度了，但不写的话程序会报错）
执行下列命令开始模拟
grompp -f md.mdp -c H2O_box2.gro -p system.top -o md1.tpr
mdrun -v -deffnm md1
因为体系很小，几秒钟就模拟完了。
3.3 观看和转换轨迹
把H2O_box2.gro载入VMD，然后在VMD main里选择它，在上面点右键，选delete frames，点delete按钮，这样目前的这一帧，即初始结构就被删掉了。再右键选Load Data into Molecule，选md1.xtc。当前共有101帧，编号为0~100。从轨迹上看，体系没有跑散，而且结构波动还是比较大的，因此适合作为molclus的输入结构。
像当前例子这么小的体系，虽然极小点结构最多超不过10个，但是由于轨迹中的很多结构优化后都会收敛到相同结构，而且还有可能碰上优化不收敛之类的问题，因此我们要从轨迹中取至少几十帧才有可能把各种构型都找全，特别是找到其中高能结构。如果目的只是为了找到能量最低构型，难度显然会低不少，用的帧数也可以少一些。
我们这里把轨迹中所有帧都保存成.xyz轨迹格式。在VMD Main窗口里点击当前项目，点右键选Save coordinate，file type设为xyz，然后点save，输入traj.xyz，这个文件就可以用作molclus的输入了。如果觉得帧数太多算不动，或者想排除掉对应升温过程的帧，也可以在保存轨迹时使用First和Last设定保存范围，以及stride设定每隔几帧保存一次。
3.4 molclus+MOPAC2012做初步优化
这一节涉及到的molclus的相关的文件在此：
(37.5 KB, 下载次数: 89)
18:20 上传
点击文件名下载附件
虽说这个体系比较小，DFT优化一个结构不会很耗时，但是我们这里要优化101帧结构，直接用Gaussian做DFT优化还是颇耗时的。而且轨迹中的结构由于热运动，偏离极小点往往比较远，优化不仅需要的步数很多，还经常容易发生震荡，或者坐标系出现问题（Gaussian的一个老毛病）导致优化了半天最终却失败。一个比较机智的方法是先用MOPAC2012在半经验方法下预优化一下，瞬间就能把这么多结构都优化完，而且几乎不会报错。但是优化出的结构的质量还只能算是作为精确计算的初猜的水平，没法直接用，因此之后还得再用DFT来refine一下。虽然Gaussian09也支持很多半经验方法，但比起MOPAC2012还是慢很多，而且没有很适合弱相互作用计算的。MOPAC2012里支持的PM6-DH+用于氢键计算在半经验的层次上算是最理想的选择之一（个人感觉比PM7要好）。
在molclus文件夹下编写一个MOPAC输入文件的模板文件，文件名应当为template.mop，内容为
PM6-DH+ precise
molecule
All coordinates are Cartesian
[GEOMETRY]
复制代码其中precise代表用比默认更高的精度进行计算。[GEOMETRY]在molclus执行过程中将被自动替换为traj.xyz的每一帧的结构并产生MOPAC输入文件MOPAC.mop。
把molclus目录下的settings.ini设定好，大部分不用改，只需要确认isys、iprog和mopac_path是否设对，请根据2.3节的说明合理地设定。
把上一节得到的traj.xyz拷到molclus所在目录下。如果是Windows版，直接双击molclus.exe。如果是Linux版，进入molclus文件夹，输入./molclus。然后molclus就开始运行了。
运行过程中会提示处理到了第几帧，优化是否顺利完毕，最终能量为多少。每个顺利优化完成的结构都会被记录到当前目录下的isomers.xyz里。本例的MOPAC优化转眼间就都顺利完成了，isomers.xyz最后共有101帧。用文本编辑器打开isomers.xyz就会看到每个优化后的结构坐标和能量。
isomers.xyz里这101帧里绝大多数都是重复的结构，我们用isostat工具按能量对之进行排序并去除重复结构。Windows版直接双击isostat.exe，Linux版运行./isostat就启动了isostat。启动后如果直接敲回车，它就会载入当前目录下的isomers.xyz，然后问你区分重复结构的能量阈值和几何阈值各是多少。由于几何优化的收敛精度所限，即便实际上两个结构对应于同一个极小点，结构也会稍有偏离，能量也会有所差异。在isostat中，如果两个结构间能量偏差和几何偏差都小于阈值的话，则这两个结构会被认为是相同的而去除其一，或者说，它们被归到了一个簇里。这里我们只使用能量阈值而不用几何偏差作为判断（这么做没特殊含义，这里仅作为示例），输入能量阈值时直接敲回车用默认的0.25kcal/mol，然后输入几何阈值时输入一个很大的数比如99。我们然后屏幕上就会看到isomers.xyz中每个构型的能量，以及是被当成了新的结构还是被认为和之前的结构是重复的：
Energy of structure& &&&1 :& &&&-0. a.u.,&&new cluster
Energy of structure& &&&2 :& &&&-0. a.u.,&&new cluster
Energy of structure& &&&3 :& &&&-0. a.u.,&&new cluster
Energy of structure& &&&4 :& &&&-0. a.u.,&&new cluster
Energy of structure& &&&5 :& &&&-0. a.u.,&&new cluster
Energy of structure& &&&6 :& &&&-0. a.u.,&&attributed to cluster& &&&3
Energy of structure& &&&7 :& &&&-0. a.u.,&&attributed to cluster& &&&4
Energy of structure& &&&8 :& &&&-0. a.u.,&&attributed to cluster& &&&5
然后看到排除了重复结构后，能量由低到高排序的结果。E是体系能量，DE是相对于能量最低的构型的相对能量，DGmin是当前构型与与它构型最相近的构型的几何偏差（这里没意义，因为刚才没启用几何阈值来区分构型）。Count代表isomers.xyz中根据能量和距离阈值判断与这个结构相同的结构共有多少。可见能量最低的3个结构在轨迹中出现几率占了主导。
The lowest energy is& &&&-0. a.u.
Sorting clusters according to energy...
#& & 1 Count:& &21&&E=& &&&-0.399205 a.u.&&DGmin=& &0.17&&DE=& & 0.00 kcal/mol
#& & 2 Count:& &23&&E=& &&&-0.396947 a.u.&&DGmin=& &0.20&&DE=& & 1.42 kcal/mol
#& & 3 Count:& &37&&E=& &&&-0.396480 a.u.&&DGmin=& &0.20&&DE=& & 1.71 kcal/mol
#& & 4 Count:& & 1&&E=& &&&-0.395613 a.u.&&DGmin=& &0.33&&DE=& & 2.25 kcal/mol
#& & 5 Count:& & 1&&E=& &&&-0.394521 a.u.&&DGmin=& &0.38&&DE=& & 2.94 kcal/mol
#& & 6 Count:& &10&&E=& &&&-0.392905 a.u.&&DGmin=& &0.17&&DE=& & 3.95 kcal/mol
#& & 7 Count:& & 2&&E=& &&&-0.391901 a.u.&&DGmin=& &0.50&&DE=& & 4.58 kcal/mol
#& & 8 Count:& & 2&&E=& &&&-0.388161 a.u.&&DGmin=& &0.90&&DE=& & 6.93 kcal/mol
#& & 9 Count:& & 2&&E=& &&&-0.387285 a.u.&&DGmin=& &0.47&&DE=& & 7.48 kcal/mol
#& &10 Count:& & 1&&E=& &&&-0.385500 a.u.&&DGmin=& &0.36&&DE=& & 8.60 kcal/mol
#& &11 Count:& & 1&&E=& &&&-0.384237 a.u.&&DGmin=& &0.90&&DE=& & 9.39 kcal/mol
这11个结构连同能量都被记录到了当前目录下的cluster.xyz中，格式与traj.xyz一致。之后再按回车可以退出isostat，如果输入温度值可以得到此温度下各个构型的波尔兹曼分布百分比。
PM6-DH+的计算结果，主要就是作为个初筛，接下来我们用DFT对这11个结构进行二次优化和筛选。
3.5 molclus+Gaussian做优化
这一节涉及到的molclus的相关的文件在此：
(8.49 KB, 下载次数: 123)
11:27 上传
点击文件名下载附件
为了节省时间，我们用便宜但对于当前体系来讲几何优化精度可以接受的M062X/6-31G*来做优化。编写Gaussian输入文件的模板文件template.gjf，内容如下（请自行合理地添加%nprocs、%mem等关键词）
# M062X/6-31G* int=ultrafine opt(tight,maxcyc=40) nosymm
Template file
0 1
[GEOMETRY]
同时，我们也写一个备用方案的模板文件template2.gjf。当基于template.gjf优化达到步数上限不收敛时，如果settings.ini里igaucontinue被设为了1，就会取最后的结构基于template2.gjf继续优化。此例template2.gjf的内容为
# M062X/6-31G* int=ultrafine opt(tight,gdiis,maxstep=2,notrust,maxcyc=100) nosymm
Template file
0 1
[GEOMETRY]
以这种方式写这两个模板文件是有一些考虑的。经过MOPAC2012预优化，初始结构偏离准确结构已经不算太远了，当基于template.gjf的默认优化设定进行优化时，如果40步还没收敛，很可能已经发生了震荡，继续优化下去可能白费功夫，因此用了maxcyc=40来设定步数上限（默认值比这个高）。在template2.gjf中，我们尝试改用其它优化方法GDIIS，并且还把步长上限设定在了很小的值，这对于解决优化末期的震荡往往很有效。nosymm关键词其实不需要写，但写了的好处在于可以令Gaussian输出文件中的结构在观看时不会突跃（这是因为默认情况下Gaussian总把结构调整到标准朝向），而且写了也不会多花时间，因为初始结构偏离对称性还是不小的，Gaussian自动判断不出对称性。
模板文件里用了int=ultrafine和opt=tight，这是因为M062X这类明尼苏达系列泛函对积分格点要求高，优化时用较高质量格点比较放心。用较严的收敛限是因为分子团簇势能面有时会很缓（尤其是范德华作用主导的情况），用默认收敛限的话构型和能量都可能不很精确，可能有零点几kcal/mol的误差，甚至会收敛到带虚频的结构上。
检查settings.ini中的各项设置。现在iprog应该从上一节用的2改为1，确认gaussian_path是否设对，igaucontinue是否已被设为了1。
删掉原先的traj.xyz，把上一节得到的cluster.xyz改名为traj.xyz作为这次的输入结构。执行molclus。根据molclus的输出我们发现，第5、6、8、9、10、11号结构第一次优化都没有收敛，但自动切换为template2.gjf继续优化后都顺利收敛了，说明利用双模板文件的方案是有效的（如果不用tight收敛限的话，第一次优化就收敛的可能性明显大得多）
算完之后，运行isostat，结果为
The lowest energy is& &-305. a.u.
Sorting clusters according to energy...
#& & 1 Count:& & 3&&E=& &-305.560998 a.u.&&DGmin=& &0.49&&DE=& & 0.00 kcal/mol
#& & 2 Count:& & 2&&E=& &-305.560941 a.u.&&DGmin=& &0.27&&DE=& & 0.04 kcal/mol
#& & 3 Count:& & 2&&E=& &-305.558619 a.u.&&DGmin=& &0.25&&DE=& & 1.49 kcal/mol
#& & 4 Count:& & 3&&E=& &-305.558223 a.u.&&DGmin=& &0.25&&DE=& & 1.74 kcal/mol
#& & 5 Count:& & 1&&E=& &-305.557159 a.u.&&DGmin=& &0.49&&DE=& & 2.41 kcal/mol
可见，之前MOPAC2012优化出的结构经过二次优化，又去掉了些重复的，目前有5个构型。它们的结构如下
1.jpg (78.66 KB, 下载次数: 4)
10:12 上传
这5个结构确实都不相同。1和5比较相似，都是非平面的，主要差异在于1比5多了一条氢键。2~4比较相似，四个氧都基本处在同一个平面上，4条氢键组成一个环，区别仅在于这四个水的氢的朝向。结构2当中，1-3号、2-4号水的氢都在平面同侧；结构3当中，1-2号、3-4号水的氢都在平面同侧；结构4当中，1-3-4号水的氢都在同侧。其中结构2的能量最低容易理解，因为对称度最高，氢键形成得最自在，而结构4的氧原子就不在一个平面上了。
为了稳妥起见，大家可以对它们都做振动分析看看是否确实无虚频，这也顺便得到了各种热力学校正量，可用于随后获得吉布斯自由能等数据。上面得到的5个结构经检验确实都无虚频。
PS：如果懒得对每个结构手动一一写输入文件做振动分析，也可以直接利用molclus。也就是把模板文件里的opt有关的关键词去掉，写上freq，然后把settings.ini里的ieneonly改为1，并确认ibkout设为了1以使得程序备份每个结构的输出文件。然后运行molclus，产生出的比如gau00003.out就是第3个结构的振动分析的输出文件。我们可以用诸如
grep &Thermal correction to Gibbs Free Energy& *.out
这样的命令把每个结构下的吉布斯自由能校正量都一次性输出出来，将它加到下一节得到的高精度单点能上就获得了高精度的自由能，然后可以根据《根据Boltzmann分布计算分子不同构象所占比例》（）准确计算出这些构型在特定温度下的分布比例。
3.6 高精度单点能计算
现在只有5个结构了，直接手工一一处理也比较方便了。M062X/6-31G*的结构还说得过去，但算的能量只是定性合理，准确比较不同异构体的能量显然还是得用更高精度再计算一下单点能。
利用molclus也可以方便地对这5个结构做单点能计算。一般认为MP2算氢键比较不错，我们用Gaussian在MP2/aug-cc-pVTZ下做单点，template.gjf内容写为# MP2/aug-cc-pVTZ
Template file
0 1
[GEOMETRY]
把settings.ini里的ieneonly改为1，说明只做单点不做优化。把energyterm由HF=改为MP2=，说明molclus将从输出文件末尾读取MP2能量而非HF/DFT能量。对于其它计算级别，如果不知道energyterm该怎么写，找个小体系试试，看一下末尾部分输出信息就知道了。
删掉原先的traj.xyz，把上一节得到的cluster.xyz改名为traj.xyz，运行molclus对其中的结构进行计算。然后运行isostat，用默认的判断阈值，我们就看到了这些结构的MP2/aug-cc-pVTZ下的相对能量：
#& & 1 Count:& & 1&&E=& &-305.360246 a.u.&&DGmin=& &0.27&&DE=& & 0.00 kcal/mol
#& & 2 Count:& & 1&&E=& &-305.358657 a.u.&&DGmin=& &0.25&&DE=& & 1.00 kcal/mol
#& & 3 Count:& & 1&&E=& &-305.358541 a.u.&&DGmin=& &0.25&&DE=& & 1.07 kcal/mol
#& & 4 Count:& & 1&&E=& &-305.353735 a.u.&&DGmin=& &0.49&&DE=& & 4.09 kcal/mol
#& & 5 Count:& & 1&&E=& &-305.352220 a.u.&&DGmin=& &0.49&&DE=& & 5.04 kcal/mol
用VMD再看一下刚得到的cluster.xyz，会发现在MP2/aug-cc-pVTZ下，居然形成平面构型的三个水团簇变成了能量最低的，而原先有6条氢键的那个现在竟成了能量最高的。
在《大体系弱相互作用计算的解决之道》（）一文中笔者大力推荐在ORCA里使用双杂化泛函PWPB95-D3(BJ)结合质量很高的def2-QZVPP基组做计算，精度很好，而且由于RI技术使得计算速度很快。这里演示一下molclus+ORCA在这个级别下对这些构型做单点能计算。
我们编写ORCA的单点计算模板文件template.inp，内容如下（假设用4核，每核分配约3GB内存）
!RI-PWPB95 d3 def2-QZVPP def2-QZVPP/JK def2-QZVPP/C rijk nopop pal4 grid4 tightscf
%maxcore 3000
* xyz 0 1
[GEOMETRY]
*
复制代码
把settings.ini里的iprog设为3。然后运行molclus即开始计算。结果如下
#& & 1 Count:& & 1&&E=& &-305.771954 a.u.&&DGmin=& &0.27&&DE=& & 0.00 kcal/mol
#& & 2 Count:& & 1&&E=& &-305.770464 a.u.&&DGmin=& &0.25&&DE=& & 0.93 kcal/mol
#& & 3 Count:& & 1&&E=& &-305.770325 a.u.&&DGmin=& &0.25&&DE=& & 1.02 kcal/mol
#& & 4 Count:& & 1&&E=& &-305.766111 a.u.&&DGmin=& &0.49&&DE=& & 3.67 kcal/mol
#& & 5 Count:& & 1&&E=& &-305.765025 a.u.&&DGmin=& &0.49&&DE=& & 4.35 kcal/mol
这个级别下，用Intel 4核机子计算每个构型才花了约2分钟。虽然定量结果和MP2/aug-cc-pVTZ的有差异，但是能量顺序一致。可见M062X/6-31G*错误地预测了能量顺序，因此在优化后做一下高精度计算是很有必要的。之所以看似有六个氢键的结构能量反倒最高，应该解释为氢键键角都太大，因此强度都较弱。
在J.Phys.Chem.A,115,12034一文中，作者获得了三种水团簇在金标准CCSD(T)/CBS级别下的能量，相对能量为0.00, 0.85, 3.55kcal/mol，对相应结构PWPB95-D3(BJ)/def2-QZVPP给出的能量是0.00, 0.93, 3.67kcal/mol。可见，PWPB95-D3(BJ)/def2-QZVPP的误差非常小。因此，对于分子团簇单点能计算，笔者强烈推荐用ORCA做PWPB95-D3(BJ)/def2-QZVPP，100个原子以内都能搞定，性价比远胜其它方法。
PS：如果想用molclus结合ORCA做优化，template.inp里写上opt关键词，并确认ieneonly设为了0、iprog设为了3即可。
3.7 分子团簇搜索小结
通过(H2O)4这个例子，基本上把分子团簇构型搜索需要注意的各种问题都讨论了，其处理过程对于多种分子构成的混合团簇也完全适用。一定要注意不同分子团簇特点不同，构型搜索需要的参数、具体策略往往差异很大，没有唯一通用的处理办法，必须根据实际情况反复摸索、灵活变通、合理巧妙运用molclus才可能得到想要的结果。对于缺乏经验的人，一次就做成功的可能性比较有限。
在不同计算级别下优化，极小点构型、数目往往会有差异（特别是那些高能、不太稳定的构型）。而且优化参数、筛选策略也影响最终是否能找全极小点。建议以不同优化方式多试试，否则很容易漏掉高能构型。
顺带一提，只要随便逮几个分子，按照此节方法用molclus搜索一下它们组成的团簇构型，算算能量差，做个振动分析算下boltzmann分布比率，然后再根据此帖《Multiwfn支持的弱相互作用的分析方法概览》（）的方法做一通相互作用本质的分析，就很容易能出一篇文章，但切勿用此方法大量灌水！下图就是用帖子里提到的RDG方法绘制的图像，氢键（蓝色等值面）、范德华作用（绿色）和弱位阻作用（棕色）都生动直观地显示了出来，对于把握结构复杂的分子团簇中的相互作用很有帮助。
2.png (43.67 KB, 下载次数: 6)
18:23 上传
4 分子构象搜索实例：麻黄素（Ephedrine）
molclus的使用方法通过第3节的例子已经讲得很透彻了，将molclus用于分子构象搜索的过程是大同小异的。这里我们对麻黄素（Ephedrine）进行构象搜索。分子结构如下，可见它有许多可旋转的键，不做构象搜索很难确定哪种构象能量最低。
3.jpg (33.79 KB, 下载次数: 3)
18:24 上传
此例Amber模拟部分所有涉及到的文件都可以在这里下载
(608.63 KB, 下载次数: 95)
18:24 上传
点击文件名下载附件
4.1 动力学模拟
虽然也可以用Gromacs进行模拟，但是产生gromacs拓扑文件常用的prodrg、ATB工具都只支持gromos力场，这是联合原子力场，对于当前体系这样含有非极性氢的情况会自动去掉这些氢，而量化计算是需要全原子的，因此之后还得再加氢，略麻烦。虽然在《几种生成有机分子Gromacs拓扑文件的工具》（）当中提到用MKTOP、OBGMX也能生成诸如全原子力场的拓扑文件，但比较寨。我比较推荐直接用Amber来做，我们这里用的是amber12，Ambertools里的antechamber工具可以方便地生成小分子在amber中模拟所需要的文件。如果非要用gromacs来模拟，也可用文中提到的acpype基于antechamber产生Gromacs格式的GAFF全原子力场拓扑文件。
PS1：免费的Ambertools从14版开始已经包含了amber的主要的动力学模块sander了，可以不必花钱买amber了。
PS2：如果不会装Amber/Ambertools，可参见《Amber14安装方法》（）。
首先建立个文件夹，在里面再建立antech文件夹，把麻黄素结构文件Ephedrine.pdb拷进去，执行
antechamber -i Ephedrine.pdb -fi pdb -o Ephedrine.prepin -fo prepi -c bcc -s 2
parmchk -i Ephedrine.prepin -f prepi -o Ephedrine.frcmod
这样就得到了麻黄素的库文件Ephedrine.prepin和补充的参数文件Ephedrine.frcmod。若用文本编辑器打开Ephedrine.prepin，会看到当前分子名叫INT。（PS：应确认输入的.pdb文件里的分子结构合理、没有原子缺失。pdb文件里是否有记录原子连接关系的CONNECT字段无所谓。）
把Ephedrine.prepin和Ephedrine.frcmod拷到上一级目录，然后执行下面的命令启动tleap
tleap -s -f leaprc.ff10
在tleap里输入
source leaprc.gaff
loadamberprep Ephedrine.prepin
loadamberparams Ephedrine.frcmod
saveamberparm INT EPH.prmtop EPH.inpcrd
就得到了拓扑文件EPH.prmtop和结构文件EPH.inpcrd。
为了让动力学能够充分采样麻黄素的势能面各个区域，而且由于不必像分子团簇那样担心跑散，我们用很高的温度模拟（1000K）。总共模拟1ns，步长2fs，每5000帧写入一次轨迹，即间隔10ps。周期边界条件对此体系是多余的因此这里不用。此例是做真空下的构象搜索，所以igb=0。相应的参数文件md1.in如下
1ns simulation at 1000K
&cntrl
imin=0,nstlim=500000,dt=0.002,ntpr=50,ntwr=100,ntwx=5000,ntc=2,
tempi=1000,temp0=1000,ntt=3,ntb=0,cut=12.0,gamma_ln=2.0,igb=0
/
复制代码
执行此命令进行模拟，很快就完成了
sander -O -i md1.in -o md1.out -p EPH.prmtop -c EPH.inpcrd -r md1.rst -x md1.mdcrd
4.2 转换轨迹
将EPH.prmtop拖入VMD，然后在相应项目上面点右键选Load Data into Molecule，然后选md1.mdcrd，总共有100帧。播放轨迹的话，会看到分子在不断旋转，很难比较不同帧之间结构的差异。笔者建议把分子的苯环部分做一下align，这样苯环部分就基本固定了，只会看到柔性链的构象在不断变动。具体做法是在VMD里选择Extensions-Analysis-RMSD Trajectory Tool，在左上角把内容改为index 16 to 25（对应于苯环上原子的编号），去掉&noh&复选框，然后点ALIGN按钮。再次播放轨迹，就可以清楚地看到由于模拟温度足够高，导致柔性链的构象不断显著变化，达到了我们的目的。笔者把这100帧组合成了动画，如下所示
4.gif (299.14 KB, 下载次数: 10)
18:25 上传
在VMD main窗口里点当前项目，点右键选Save coordinates，格式选xyz，然后保存成traj.xyz。
4.3 molclus+MOPAC2012做初步优化
像第3节的例子一样，我们先用半经验方法进行初筛，可以将要考虑的结构数目降低近一个数量级。我们编写MOPAC的模板文件template.mop：
PM6-DH+ precise
molecule
All coordinates are Cartesian
[GEOMETRY]复制代码
确认settings.ini中的iprog设为了2，ieneonly设为了0，并且已经设好了mopac_path。然后把traj.xyz拷到molclus目录中，运行molclus对其中的结构进行优化。然后运行isostat对结果进行统计，直接敲三下回车用默认的能量和几何判断阈值，结果如下
#& & 1 Count:& & 8&&E=& &&&-0.059452 a.u.&&DGmin=& &0.24&&DE=& & 0.00 kcal/mol
#& & 2 Count:& & 1&&E=& &&&-0.057176 a.u.&&DGmin=& &0.29&&DE=& & 1.43 kcal/mol
#& & 3 Count:& & 9&&E=& &&&-0.056859 a.u.&&DGmin=& &0.18&&DE=& & 1.63 kcal/mol
#& & 4 Count:& & 4&&E=& &&&-0.056548 a.u.&&DGmin=& &0.29&&DE=& & 1.82 kcal/mol
#& & 5 Count:& & 1&&E=& &&&-0.056523 a.u.&&DGmin=& &0.12&&DE=& & 1.84 kcal/mol
#& & 6 Count:& & 1&&E=& &&&-0.056513 a.u.&&DGmin=& &0.12&&DE=& & 1.84 kcal/mol
#& & 7 Count:& & 2&&E=& &&&-0.055854 a.u.&&DGmin=& &0.16&&DE=& & 2.26 kcal/mol
#& & 8 Count:& & 4&&E=& &&&-0.055561 a.u.&&DGmin=& &0.23&&DE=& & 2.44 kcal/mol
#& & 9 Count:& & 4&&E=& &&&-0.055465 a.u.&&DGmin=& &0.20&&DE=& & 2.50 kcal/mol
#& &10 Count:& & 2&&E=& &&&-0.054854 a.u.&&DGmin=& &0.26&&DE=& & 2.88 kcal/mol
#& &11 Count:& & 1&&E=& &&&-0.054807 a.u.&&DGmin=& &0.28&&DE=& & 2.91 kcal/mol
#& &12 Count:& & 5&&E=& &&&-0.054777 a.u.&&DGmin=& &0.24&&DE=& & 2.93 kcal/mol
#& &13 Count:& & 4&&E=& &&&-0.054723 a.u.&&DGmin=& &0.19&&DE=& & 2.97 kcal/mol
#& &14 Count:& & 1&&E=& &&&-0.054690 a.u.&&DGmin=& &0.55&&DE=& & 2.99 kcal/mol
#& &15 Count:& & 1&&E=& &&&-0.054521 a.u.&&DGmin=& &0.24&&DE=& & 3.09 kcal/mol
#& &16 Count:& &12&&E=& &&&-0.054478 a.u.&&DGmin=& &0.21&&DE=& & 3.12 kcal/mol
#& &17 Count:& & 7&&E=& &&&-0.054461 a.u.&&DGmin=& &0.18&&DE=& & 3.13 kcal/mol
#& &18 Count:& & 1&&E=& &&&-0.054459 a.u.&&DGmin=& &0.16&&DE=& & 3.13 kcal/mol
#& &19 Count:& & 1&&E=& &&&-0.054366 a.u.&&DGmin=& &0.23&&DE=& & 3.19 kcal/mol
#& &20 Count:& & 5&&E=& &&&-0.054328 a.u.&&DGmin=& &0.27&&DE=& & 3.21 kcal/mol
#& &21 Count:& & 1&&E=& &&&-0.054317 a.u.&&DGmin=& &0.25&&DE=& & 3.22 kcal/mol
#& &22 Count:& & 1&&E=& &&&-0.054276 a.u.&&DGmin=& &0.16&&DE=& & 3.25 kcal/mol
#& &23 Count:& & 1&&E=& &&&-0.054258 a.u.&&DGmin=& &0.16&&DE=& & 3.26 kcal/mol
#& &24 Count:& & 5&&E=& &&&-0.052824 a.u.&&DGmin=& &0.28&&DE=& & 4.16 kcal/mol
#& &25 Count:& & 2&&E=& &&&-0.052341 a.u.&&DGmin=& &0.17&&DE=& & 4.46 kcal/mol
#& &26 Count:& & 5&&E=& &&&-0.052125 a.u.&&DGmin=& &0.41&&DE=& & 4.60 kcal/mol
#& &27 Count:& & 5&&E=& &&&-0.051873 a.u.&&DGmin=& &0.16&&DE=& & 4.75 kcal/mol
#& &28 Count:& & 4&&E=& &&&-0.051050 a.u.&&DGmin=& &0.21&&DE=& & 5.27 kcal/mol
#& &29 Count:& & 2&&E=& &&&-0.050232 a.u.&&DGmin=& &0.16&&DE=& & 5.78 kcal/mol
现在只有29个结构了，可以把得到的cluster.xyz放到VMD里进行观看了解它们的结构。
如果感兴趣的只是能量最低的一些构象，可以把cluster.xyz靠后的高能结构删掉，以降低下一节量化计算的耗时。
4.4 molclus+Gaussian做优化
我们接下来用B3LYP-D3(BJ)/6-31G*对MOPAC2012优化后的结构进一步优化。虽然对于当前体系DFT-D3色散校正不是必须，但由于是“免费”的，所以不如加上。我们编写几何优化的模板文件template.gjf，内容如下
# B3LYP/6-31G* empiricaldispersion=GD3BJ opt=maxcyc=40
Template file
0 1
[GEOMETRY]
另编辑一个备用方案的模板文件template2.gjf，把关键词部分替换为# B3LYP/6-31G* empiricaldispersion=GD3BJ opt(gdiis,maxstep=2,notrust,maxcyc=80)即可。其实，当前体系凭经验就能知道优化不会难收敛，远比优化团簇容易得多，所以其实用不用备用方案也无所谓。
把settings.ini中的iprog改为1，并且确认energyterm目前为HF=、ieneonly为0。删掉之前的traj.xyz，把上一节得到的cluster.xyz改名为traj.xyz。然后运行molclus，再运行isostat下默认阈值下进行处理，结果减少到了24个：
#& & 1 Count:& & 3&&E=& &-520.107194 a.u.&&DGmin=& &0.22&&DE=& & 0.00 kcal/mol
#& & 2 Count:& & 1&&E=& &-520.107187 a.u.&&DGmin=& &0.21&&DE=& & 0.00 kcal/mol
#& & 3 Count:& & 1&&E=& &-520.105306 a.u.&&DGmin=& &0.17&&DE=& & 1.18 kcal/mol
#& & 4 Count:& & 1&&E=& &-520.102523 a.u.&&DGmin=& &0.28&&DE=& & 2.93 kcal/mol
#& & 5 Count:& & 1&&E=& &-520.102399 a.u.&&DGmin=& &0.19&&DE=& & 3.01 kcal/mol
#& &24 Count:& & 3&&E=& &-520.095625 a.u.&&DGmin=& &0.23&&DE=& & 7.26 kcal/mol
笔者把其中部分结构做成了动画，如下所示
5.gif (19.03 KB, 下载次数: 6)
18:27 上传
我们用VMD看看cluster.xyz中的第2帧，思考一下为什么它能量几乎最低。从结构上可以发现，它的柔性链的结构很舒展，没什么位阻，而且羟基氢和氮的孤对电子离得比较近，可能有内氢键的特征。为了确认这一点，我们用Multiwfn作RDG图，如下所示
6.jpg (40.98 KB, 下载次数: 5)
18:27 上传
可见，羟基氢和氮的孤对电子区域之间的等值面上明显有一块蓝色，很明确地证明了存在内氢键，这是这个结构能量低的重要因素之一。
我们可以看看这24个结构中苯环的取代基的构象是怎么分布的。用4.2节的做法令苯部分进行align，然后在VMD的Graphics-Representations-Trajectory里的Draw Multiple Frames里输入0:24，就可以将这24个结构同时显示出来，如下所示
7.png (31.27 KB, 下载次数: 6)
13:44 上传
可见，这些构象中苯环部分结构没有任何差异，只是取代基部分构象发生很大变化，分布很广，几乎把所有可能性都覆盖了。
这个例子证明利用molclus结合分子动力学和量子化学程序可以很好地也比较方便地搜索小分子构象。如果需要更准确的能量，可以再像(H2O)4那个例子一样用高精度的方法对各个结构计算单点。
本文介绍了molclus程序的用法，通过(H2O)4分子团簇的构型搜索和麻黄素的构象搜索这两个例子充分证明了molclus的重要实用性。可能有人觉得本文介绍的过程不够傻瓜化，需要学一些动力学程序和量化程序的基本使用、需要有一些分子模拟和量化的基础才行，没法像一些商业化的程序那样点点按钮就出结果。但实际上，那样看似傻瓜的程序反倒不好用，可控程度差，不够灵活，精度也不高，教学用或者粗浅地研究可能比较适合，而做为深入研究目的，则往往派不上用场。而本文的molclus+动力学+量化程序的构型/构象搜索方法，只要稍加揣摩，多加练习，就能在日后的研究中得心应手，起到重要价值。
另外再提及一点，虽然本文一直提到的都是通过动力学来产生初始构型/构象，但实际上，molclus本身并没有局限于只能用动力学程序产生初始构型/构象。对于特殊问题，如果动力学不适合或不方便，大家也可以用其它算法，比如势能面扫描、蒙特卡罗、基因算法等等，以及其它各种类型程序来产生初始构型/构象，只要能把得到的初始结构都写进一个.xyz轨迹文件就行了，然后molclus就能自动对其中的结构通过量化程序进行优化、能量计算和筛选。
价格相当公道:-)
赞！简直不能更好
感谢分享！
思想家公社的门口Blog：
北京科音自然科学研究中心：
Multiwfn主页：
计算化学公社：
思想家公社QQ群1号：00人已满），思想家公社QQ群2号：（讨论计算化学为主，加入时必须注明研究方向）
Money and papers are rubbish, get a real life!
PS: 此账号为诸Sob共用
在线时间228 小时
帖子威望eV
Level 4 (黑子), 积分 1137, 距离下一级还需 363 积分
感谢sobereva老师分享！
在线时间96 小时
帖子威望eV
Level 3 能力者, 积分 478, 距离下一级还需 22 积分
在线时间220 小时
帖子威望eV
Level 4 (黑子), 积分 1303, 距离下一级还需 197 积分
顶，绝对的支持~
在线时间195 小时
帖子威望eV
Level 4 (黑子), 积分 990, 距离下一级还需 510 积分
原理略有些简单啊。。用MD方法搜索团簇的效率还是比MC低了不少，这样很难得到团簇的global minima吧
还有一个问题就是MD所得到的结构用于量化的OPT的话很容易不收敛的，因为距离极小值的距离比较远
如果要进一步开发的话，建议放弃MD软件，用basin hopping进行数据采样，这样得到的结构是一个分子力学上的极小值，既可以节约采样的时间，也可以节约量化结构优化的时间
在线时间1584 小时
帖子威望eV
原理略有些简单啊。。用MD方法搜索团簇的效率还是比MC低了不少，这样很难得到团簇的global minima吧
我开发过MC搜索团簇的程序，对于搜索原子团簇适合，但这种方法明显不适合分子团簇（本文强调的是分子团簇），或者说效率比MD明显要低，而且收敛难度反倒远比基于MD大得多，毕竟MC过程根本没有很好地考虑原子间的作用，别提几何优化收敛了，由于随机移动原子造成的构型不合理，甚至一开始就可能SCF不收敛，甚至分子团簇产生了错误的结构（比如氢原子干脆优化后发生转移了）。
分子动力学得到分子团簇全局极小点没有任何问题，是分子团簇构型搜索有效的方法，更是众多文献讨论分子团簇（包括很大团簇）的标准方法。basin hopping、metadynamics之类采样方未必效率更高，而且有局限性。
“原理简单”不是缺点反倒是优点，用简单的原理，很好地解决实际需要解决的问题，这才是好的、适合大众使用的程序。
反之原理复杂，甚至引入一堆参数，需要使用者需要大量经验积累，颇为麻烦地才能使用来解决看似不复杂的问题，甚至还得改反复改代码，那样的程序我不推崇，不是大众向的，而只是少数专家们发paper用的。
思想家公社的门口Blog：
北京科音自然科学研究中心：
Multiwfn主页：
计算化学公社：
思想家公社QQ群1号：00人已满），思想家公社QQ群2号：（讨论计算化学为主，加入时必须注明研究方向）
Money and papers are rubbish, get a real life!
PS: 此账号为诸Sob共用
在线时间267 小时
帖子威望eV
Level 4 (黑子), 积分 964, 距离下一级还需 536 积分
那个settings.ini中的ipause参数=3表示总不暂停？，但是文件中写的是=0表示不暂停。好像是，，
在线时间1584 小时
帖子威望eV
那个settings.ini中的ipause参数=3表示总不暂停？，但是文件中写的是=0表示不暂停。好像是，，
帖子里写错了，已更正。
其实无所谓，只要不是1和2就行，所以0和3效果都一样
思想家公社的门口Blog：
北京科音自然科学研究中心：
Multiwfn主页：
计算化学公社：
思想家公社QQ群1号：00人已满），思想家公社QQ群2号：（讨论计算化学为主，加入时必须注明研究方向）
Money and papers are rubbish, get a real life!
PS: 此账号为诸Sob共用
在线时间267 小时
帖子威望eV
Level 4 (黑子), 积分 964, 距离下一级还需 536 积分
综合性有点强，很难把几类软件串在一起。
在线时间406 小时
帖子威望eV
Level 5 (御坂), 积分 1600, 距离下一级还需 2400 积分
Sob 老师强，赞一个！
在线时间45 小时
帖子威望eV
Level 3 能力者, 积分 295, 距离下一级还需 205 积分
1）该软件的思路简单但巧妙，确能有效解决团簇构型和构象搜索问题。
2）MD我用过gromacs的比较老的版本，它的安装和上手还是有些麻烦；对于体系包含1：1两种分子（或三种）的比较复杂的体系，有时不太好建模。能否MD软件也能支持MS做的MD轨迹文件呢，它还是比较好上手和建模的；我记得作者曾经也写过个MS轨迹文件转换的程序，这应当不算太困难。
3）QM能支持gaussian、orca、mopac，这个基本满足各种精度的需要了。
在线时间1584 小时
帖子威望eV
1）该软件的思路简单但巧妙，确能有效解决团簇构型和构象搜索问题。
2）MD我用过gromacs的比较老的版本， ...
gromacs还好。而且也有第三方编译的windows版，上来直接就能用。至于上手，其实效仿文章的例子，很快就能弄会，参数文件也就改改模拟时间、保存频率、温度这些，没多少需要变动的。
建模的话强烈推荐packmol（），想建立什么样的都能建立，十分方便灵活。
MS模拟的轨迹可以通过perl脚本导出成xyz（），直接就能给molclus使用，不需要对molclus做任何修改。但我觉得与其学MS，还不如就用gmx或amber，MS价格不菲。
思想家公社的门口Blog：
北京科音自然科学研究中心：
Multiwfn主页：
计算化学公社：
思想家公社QQ群1号：00人已满），思想家公社QQ群2号：（讨论计算化学为主，加入时必须注明研究方向）
Money and papers are rubbish, get a real life!
PS: 此账号为诸Sob共用
在线时间183 小时
帖子威望eV
Level 4 (黑子), 积分 1251, 距离下一级还需 249 积分
强就一个字
在线时间1 小时
帖子威望eV
Level 2 能力者, 积分 30, 距离下一级还需 120 积分
SOB威武！这个软件对阴离子、阳离子构成的化合物是否适用？
在线时间1584 小时
帖子威望eV
SOB威武！这个软件对阴离子、阳离子构成的化合物是否适用？
同样适用。其实这和molclus本身没直接关系。对于一个体系，只要动力学程序能模拟，量化程序能算，那么就都可以用molclus做搜索。
思想家公社的门口Blog：
北京科音自然科学研究中心：
Multiwfn主页：
计算化学公社：
思想家公社QQ群1号：00人已满），思想家公社QQ群2号：（讨论计算化学为主，加入时必须注明研究方向）
Money and papers are rubbish, get a real life!
PS: 此账号为诸Sob共用
Powered by}