注：本文更新频繁，每次出molclus新版本都会更新。文中评论区的很多讨论都是对于已过时的老版本的，不要作为最新版本的参考。
使用molclus程序做团簇构型搜索和分子构象搜索 Using molclus program to perform configuration search of clusters and conformational search of molecules
文/Sobereva @北京科音
First release：2015-Jan-5 Last update：2023-Jul-17
前言：molclus（主页：http://www.keinsci.com/research/molclus.html）是北京科音自然科学研究中心开发的一款免费的用于原子/分子团簇构型搜索和分子构象搜索的程序包，目前已被大量知名期刊的高水平研究文章所使用。本文将对此程序包的基本特征和其中的molclus组件的用法进行介绍。第1节介绍基于分子动力学做构型/构象搜索的基本原理，第2节介绍molclus程序的基本使用，第3、4节通过实例进行演示。
1 构型/构象搜索原理
团簇的构型搜索和分子的构象搜索是计算化学经常涉及的问题。团簇构型搜索的目的是找到原子团簇（如Li6）或者分子团簇（如水的五聚体）的最稳定的一批构型，包括坐标和相对能量。构象搜索是对于柔性分子而言的，柔性分子有许多可旋转的键，因此势能面上有很多极小点，构象搜索目的是找到最稳定的一批构象，得知哪个构象能量最低、构象的能量彼此间相差多少。在得到最稳定的一批团簇构型或者分子构象之后，我们还可以根据它们的自由能，进一步根据Boltzmann分布估计出特定温度下它们的分布比率，见《根据Boltzmann分布计算分子不同构象所占比例》（http://sobereva.com/165）。
构型/构象搜索的算法很多，比如势能面扫描、分子动力学、蒙特卡罗、基因算法、LMOD等等，适用场合不同。团簇构型搜索和分子构象搜索都可以基于分子动力学实现，基本过程是：在合适的方法下做一定长度的分子动力学模拟，由于热运动，体系能够反复穿越势垒到达势能面不同区域。每隔一定步数取一次结构，将这些分布在势能面的不同位置上的结构作为初猜结构分别进行几何优化，最终就会得到势能面上各个极小点位置的构型和能量，之后可以对能量进行排序取最低的一批。还可以通过分子动力学做周期性退火来实现构象搜索，效果往往更好。比如每200 ps期间让温度从0 K上升至500 K然后再平缓降到0 K，如果跑2000 ps，就会产生10个对应0 K的结构，对这些结构再进行依次优化，即可得到一批极小点结构和能量。退火过程中每次当体系温度较高时，体系就有较大动能，从而能够较容易地穿越势垒而跑到其它势阱区域，然后温度降低时，体系就会陷入势阱的极小点附近。对于分子团簇的构型搜索和分子的构象搜索而言，并不会涉及到化学键的形成和断裂，因此不需要用极耗时的从头算分子动力学，而只需要计算量极低的分子力场来跑分子动力学就够了。这样的程序很多，比如免费的GROMACS、Tinker、NAMD；Amber实际也算基本免费了（用免费的Ambertools里的Sander模块来跑就行了，不需要用收费的PMEMD部分）。用xtb基于半经验DFT层面的GFN-xTB方法跑动力学耗时也并不很高，而且没有折腾参数和弄拓扑文件的麻烦，还可以在动力学过程中描述键的形成和断裂，故也可以用于原子团簇的搜索。至于结构优化和能量计算部分，可以用分子力场、半经验、DFT、高精度后HF方法等来做，具体用什么就取决于体系特征、需要的精度和运算能力了。
另外，一种简单粗暴但很实用的做分子构象搜索的方法是“系统式搜索法”，即让分子的某些二面角按照指定规则旋转并将所产生的结构做为优化的初始结构。比如设定4个键可以旋转，每个键每120度旋转一次，那么总共就可以产生3*3*3*3=81个初始构象，对这些构象都做优化，就可以得到一批极小点结构，都进行优化后，肯定能找到能量最低构型。molclus程序包中的gentor组件就可以实现这个目的。也有一些专门的构象生成程序如Confab、Frog2，比gentor考虑的因素更多，也更傻瓜化。
对于团簇构型搜索，一种也是看起来很暴力但实际效果不错的做法是把指定类型和指定数目的分子或原子随机挨在一起，形成球形团簇，将这样产生的几十乃至几百个团簇当做初猜结构来依次进行优化，就可以得到一批团簇极小点构型，只要生成的初始团簇数目较多，其中总会有能量最低结构。molclus里的genmer组件就可以实现产生团簇初猜结构的目的。
简单来说，基于molclus做构型/构象搜索所要干的主要事情就是：用户先基于第三方分子动力学程序，或者基于molclus自带的gentor/genmer组件，或者用Confab等构象生成程序，先产生一批初始结构，然后molclus就可以调用量子化学程序或分子力学程序对这些结构依次进行优化（还可以做单点、振动分析等任务），然后用molclus自带的isostat组件对结果进行归簇、排序，从而最终得到能量最低结构和能量较低的一批结构。这个过程看起来没什么复杂的，但molclus将整个过程变得极尽方便省事，而且使用者可以对整个流程进行十分灵活的控制，是想要较低的耗时还是要较高的精度，都可以完全由使用者决定。目前市面上有不少收费的构象搜索程序，比如Spartan、GMMX、HyperChem、Sybyl、MOE、M$里的Conformer、MacroModel之类，我相信，只要你稍微花点时间用明白了molclus，之后一定就再也不想用那些程序了，又不灵活又贵又没法达到像样的精度。
2 molclus程序
2.1 简介
molclus主页为：http://www.keinsci.com/research/molclus.html，程序的可执行文件可在此免费下载。molclus会不断改进、更新。
如果molclus在你的文章中被使用，必须引用，基本格式为（可根据期刊要求修改）：
Tian Lu, molclus program, Version x.x, http://www.keinsci.com/research/molclus.html (accessed 月 日, 年)
这里的月、日、年就是你下载Molclus程序的日期。笔者在未来会专门撰写介绍Molclus的文章发表，届时请引用论文（发表后会在molclus官网和此帖中提及）。
在通过分子动力学程序，或者molclus自带的gentor/genmer组件，或者Confab等程序产生记录了大量初始结构的.xyz格式的轨迹文件之后，molclus就可以调用计算程序对其中的结构进行优化或者能量计算了。molclus目前支持调用以下程序：
• Gaussian：支持G09及之后版本。这是最最常用的做量子化学计算的程序，收费。Gaussian也支持UFF/AMBER/Dreiding分子力场计算，比半经验方法又快出两数量级乃至更多
• ORCA：网址https://orcaforum.cec.mpg.de/。免费。这是第二流行的量子化学程序。由于其RI技术利用得很好，在纯泛函、双杂化泛函下比Gaussian快得多得多，还支持又非常精确又能算得动中等大小体系的DLPNO-CCSD(T)。
• MOPAC：支持2012及之后的版本，网址http://openmopac.net。免费。专门做半经验计算，它支持的PM6-DH+或PM7对有机体系通常足矣得到定性合理的结果，而速度比DFT快两个数量级，故很适合用于中、大体系的优化
• xtb：做GFN-xTB方法计算的程序，免费。此方法相当于半经验的DFT，精度整体强于各种半经验方法，特别是普适性更好，而且xtb速度比MOPAC更快，在此文中有安装方法和基本用法介绍：《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》（http://sobereva.com/421）
• Open Babel：网址http://openbabel.org/wiki/。免费。此程序主要用于转换化学领域的文件格式，但也可以实现MMFF94高精度有机分子力场以及其它一些力场的计算
可以说，molclus支持的Gaussian、ORCA、MOPAC、xtb和Open Babel程序把不同尺度、不同精度范围都覆盖了，精度&耗时的基本关系是：高级别后HF > DFT > 半经验/GFN-xTB > 分子力场，这些方法都可以由上述程序实现。如果读者不熟悉ORCA和RI，建议参看《详谈Multiwfn产生ORCA量子化学程序的输入文件的功能》（http://sobereva.com/490）当中的介绍，如果不熟悉MOPAC，可参看本文中的相关信息：《大体系弱相互作用计算的解决之道》（http://sobereva.com/214）。分子团簇和弱相互作用关系极为密切，必须选择合适的计算级别，否则团簇搜索将毫无意义。不熟悉弱相互作用计算的话建议看看《乱谈DFT-D》（http://sobereva.com/83）、《谈谈“计算时是否需要加DFT-D3色散校正？》。（http://sobereva.com/413）。
2.2 使用流程
以下是molclus使用的流程框图：
Molclus通常有以下使用方式：
• genmer结合molclus：用于团簇分子或原子团簇构型搜索。使用简便
• gentor结合molclus：用于分子构象搜索。使用简便，但不支持环状区域构象搜索（环状区域是指诸如环己烷这样有多种构象的环状区域）
• Confab或Frog2构象生成程序结合molclus：用于分子构象搜索。使用最为傻瓜化，但不支持环状区域构象搜索、不支持普通有机分子以外的情况
• 经典力场分子动力学程序结合molclus：用于分子团簇搜索和分子构象搜索。使用稍繁琐，被动力学模拟的体系必须有恰当的力场
• xtb程序跑动力学结合molclus：普适性极强，用于构象搜索、原子/分子团簇构型搜索皆可，使用容易（不受力场约束，不需要准备拓扑文件）
molclus的基本使用流程如下：
(1) 产生记录所有初始结构的多帧traj.xyz文件。通常有这些做法：
● 用分子动力学程序模拟团簇或单个分子，每隔一定步数将结构写入一次轨迹；也可以做周期性退火，将对应0 K的帧写入轨迹。令轨迹名为traj.xyz（如果程序不支持输出xyz轨迹格式，可用免费的VMD程序将其私有轨迹格式转化为.xyz格式的轨迹）。
● 用molclus程序自带的gentor产生含有一批分子初始构象的traj.xyz文件，或者用genmer工具产生含有一批团簇初始构型的traj.xyz文件。
● 用Confab直接产生含有分子的一批随机的构象的traj.xyz文件。
(2) 对于让molclus调用Gaussian、ORCA或MOPAC的情况，需提供这些量子化学程序的模板文件并放到当前目录下。对于Gaussian程序，名为template.gjf；对于ORCA，名为template.inp；对于MOPAC，名为template.mop。模板文件里要写好调用这些程序计算时要用的关键词，如理论方法、基组、相关关键词（内存分配量、辅助SCF收敛的设定、优化的关键词和选项等），但是坐标部分写成[GEOMETRY]。molclus计算时会自动将[GEOMETRY]字段替换为traj.xyz里的结构，然后传递给这些量子化学程序执行。对于调用xtb和Open Babel的情况，不需要设定模板文件，因为这俩本身就是通过命令行调用的。
(3) 设定好molclus的参数文件settings.ini指定任务类型和运行设置，详见2.3节。
(4) 启动molclus。molclus会载入traj.xyz、settings.ini和模板文件，然后自动调用量子化学或分子力学程序对traj.xyz的每一帧结构进行优化，并将优化好坐标和能量写入到当前目录下的isomers.xyz文件中（molclus很灵活，也可以对每帧结构做比如优化+振动分析+高精度单点能计算，这种情况molclus会把坐标、自由能、自由能的热校正量和虚频数目写入到isomers.xyz中，也同时输出在屏幕上）
(5) 启动molclus中附带的isostat工具，载入上一步得到的isomers.xyz，此工具会对其中的结构根据能量进行排序并进行统计、去除能量和几何结构相似度过高的重复结构。isostat在当前目录下输出的cluster.xyz文件就是根据能量从低到高排列并且去除了重复结构后的文件，可以载入VMD来方便地观看每个结构。之后还可以根据提取的各结构的能量，输出指定温度下各个构型的波尔兹曼分布比例（原理见《根据Boltzmann分布计算分子不同构象所占比例》http://sobereva.com/165）。isostat支持并行计算来大大降低耗时，在启动时可以设置并行核数。
如果想使用比结构优化更高级别的方法精确计算每个结构的能量，可以把settings.ini里的itask改为1，把模板文件改成单点任务所用的关键词，再把cluster.xyz改名为traj.xyz，然后启动molclus。此时molclus就会调用量子化学程序对其中每个结构再做一次单点计算，可以再次用isostat进行统计筛选。
程序的安装与运行：Molclus本身不需要安装就可以直接启动。在Windows下启动molclus就是双击molclus.exe图标。在Linux下启动就是进入molclus所在目录，输入./molclus（有时需要先用chmod +x ./molclus加上可执行权限）。每次程序运行前，会对当前目录下.out、.arc等量子化学或分子力学程序之前产生的临时文件进行清理，所以如果有重要的临时文件在当前目录下，执行molclus前应当先手动将它们挪至别处。
Windows版molclus包里有isostat.exe和isostat_32bit.exe。如果你的机子是64位系统，一定要用前面的。后面那个纯粹是给仍然在用32bit系统的人用的，对内存可使用量限制较大，对巨大的isomers.xyz进行筛选可能会崩溃，而且不支持并行计算。
如果你想中断molclus的运行，可以直接在molclus运行窗口里按Ctrl+C终止。还有一种做法是在当前目录下创建名为STOP的空文件，molclus每开始循环traj.xyz里的一个结构时，一旦发现当前目录下有名叫STOP的文件就会自动结束任务。
如果你想让molclus在远程服务器上一直运行，不会因为断开与服务器的连接后自动中断，应当运行比如nohup ./molclus > out.txt &，然后要断开终端时输入exit以优雅的方式断开，这样molclus就会在服务器上一直跑直到结束。molclus输出信息都会被输出到out.txt里，可以随时打开查阅，也可以运行tail -f out.txt来实时监控新写入到此文件里的内容便于跟踪计算进度。
对于在超算等特殊情况下使用，必须需要通过提交的方式执行上述量子化学程序，有可能molclus没法自动调用那些程序。为此molclus程序包里还提供了xyz2QC工具来解决这个问题，请参看《将多帧xyz文件转化成量子化学输入文件的工具：xyz2QC》（http://bbs.keinsci.com/thread-12468-1-1.html）。如果你之前没接触过xyz文件格式，建议你看一看《谈谈记录化学体系结构的xyz文件》（http://sobereva.com/477）。另外，有molclus用户针对在集群上通过PBS作业提交方式使用molclus的方法做了介绍，见：http://bbs.keinsci.com/thread-15598-1-1.html。
如果你不会装Gaussian，看http://sobereva.com/439；如果不会装ORCA，看http://sobereva.com/451；如果不会装MOPAC，看http://sobereva.com/262。如果不会装xtb，看http://sobereva.com/421。Open Babel在Windows下安装很简单，直接启动安装程序然后一直下一步即可，对于Redhat派系的Linux系统（如CentOS），用yum install epel-release添加EPEL源后用yum install openbabel即可安装。
2.3 参数文件settings.ini
molclus目录下的settings.ini在molclus启动时会被载入，其中记录了运行参数，在下面进行介绍（此文件里的注释也有描述）。实际情况需要改什么就在原文件基础上修改，注意格式别改乱了，用不到的设置也不要自己删掉。通常做一个新任务只需要改两、三个参数就够了。
iprog：molclus调用的量子化学或分子力学程序。1=Gaussian；2=MOPAC；3=ORCA；4=xtb; 5=Open Babel
ngeom：0=处理traj.xyz里每一帧；n=处理前n帧；也可以直接指定考虑的帧号，如2,5-8,12-14（如果只处理比如第6帧需要写6-6）。利用这个结合下面说的iappend=1可以实现续算等目的
itask：对traj.xyz里每一帧所做的任务类型。0=几何优化；1=计算单点能；2=振动分析；3=优化+振动分析；-1=热力学组合方法计算（仅对Gaussian有效）
注：对于Gaussian、ORCA、MOPAC这样需要设置template模板文件的情况，itask仅决定按照什么任务的形式来从输出文件中提取数据，而所做的任务的关键词需要照常写（如itask=0的时候，模板文件里得照常写优化的关键词opt及相关选项，例如对于Gaussian可以是opt(gdiis,maxstep=5,notrust)）。而对于xtb和Open Babel，并不需要设模板文件，itask直接就决定了在调用时做什么任务。
ibkout：1=每处理完一帧结构，把程序的输出文件备份，以便之后检查或根据自己的需要提取其它数据。比如对于第5帧结构，Gaussian的输出文件会备份到gau00005.out，ORCA的会备份到orca00005.out，MOPAC的会备份到MOPAC00005.arc；0=不做备份，因此molclus运行完毕后将找不到计算程序中间输出的文件；2=只备份成功的任务；3=只备份失败的任务；-1=对于调用Gaussian或ORCA的情况，如果发现当前目录下有以前备份好的文件（如gau00002.out、orca00005.out），则Molclus会直接读里面的数据而不再重新计算
distmax：n=如果结构中任意两个原子间的距离超过n埃则跳过这一帧。做团簇动力学的时候，有可能某个分子跑飞了，对这样的结构优化没有意义，利用这个参数就可以忽略掉此结构。如果设0的话，可以从一个文本文件里读取更精细的设置，能设置某一批和另一批原子间必须满足的最近和最远距离，违背要求的结构将被跳过，详见http://bbs.keinsci.com/thread-28128-1-1.html
ipause：1=对于Gaussian和ORCA，计算出错（如不收敛）则暂停，以便用户对输出文件进行检查，然后可以按回车继续处理下一帧结构；2=每处理完一个结构后都暂停；0=总是不暂停
iappend：0=新产生的结构会覆盖原有的isomers.xyz文件；1=新产生的结构会续写到当前目录下之前的isomers.xyz文件（如果有的话）末尾，若没有则会新建之
freeze：设置含有优化过程的任务中原子的冻结情况。0代表不做冻结，设成诸如1,4-6,8-10代表对1,4,5,6,8,9,10原子的笛卡尔坐标进行冻结。对Gaussian/ORCA/MOPAC/xtb有效
为了避免有些人糊涂，在这里提供一个表格，明确说明molclus支持哪些任务（itask）和哪些程序的组合。绿色代表支持，棕色代表不支持。绿色里如果有文字的话，说明需要在模板文件里体现相应的关键词
----以下是Gaussian专用的参数----
gaussian_path：调用Gaussian的命令，比如"/sob/g09/g09"、"D:\study\g16w\g16.exe"。在linux下如果Gaussian已正常安装了，可以只写诸如g09或g16。如果路径中有空格，路径两边必须加上双引号，后同
[注意：如果是Windows版Gaussian，gaussian_path必须指向诸如g09.exe而不能是g09w.exe，因为后者只是前者的一个图形界面而已。另外，运行前还必须进入“控制面板”-“高级系统设置”-“高级”-“环境变量”，添加GAUSS_EXEDIR环境变量，使之指向Gaussian目录（不是指向可执行文件名），如D:\study\g09w，否则Gaussian没法以命令行方式运行，也因此molclus将没法调用Gaussian。另外，不要用太老的G09版本，比如年代久远的G09 A.01、A.02版可能没法被molclus支持]
igaucontinue：1=如果优化失败，则会基于另一个模板文件template2.gjf对之前最后一帧结构继续优化，如果还失败了会再基于template3.gjf继续优化；0=不做这个处理。igaucontinue=1的主要用处在于给用户更多备选方案。比如可以在template.gjf里写常规的优化关键词（并且步数上限设得比较有限，如80），而在template2.gjf、template3.gjf里使用更有助于收敛的关键词，诸如calcfc、GDIIS、opt(maxstep=x,notrust)等等。这样当第一次优化不收敛的时候，就会自动利用其它模板文件里的关键词进一步优化。关于辅助收敛的方法，在此贴有汇总：《量子化学计算中帮助几何优化收敛的常用方法》（http://sobereva.com/164）。另外，有时候优化过程中途出现三个原子排成直线而报错，靠igaucontinue=1可以实现基于之前的结构继续优化，通常也能进一步优化下去。
energyterm：molclus是从Gaussian输出文件中末尾的archive字段部分读取电子能量的，这个选项用于设定读取哪项。对于HF/DFT计算，这个参数应设为HF=，即读取HF=这个字符串后面紧跟着的那个数（HF/DFT能量）。如果是MP2或双杂化泛函计算，就应当设MP2=；如果是CCSD(T)，就应当设CCSD(T)=，以此类推，不懂的话仔细看http://sobereva.com/488。如果模板文件设的是TDDFT激发态计算，且这个参数设的是TD=，则molclus会读取激发态能量（TD里root选项指定的态），从而可以做激发态构型/构象搜索。
molclus也可以从Gaussian的热力学组合方法输出中读取指定的热力学量，当itask被设为了-1，template.gjf里的关键词比如设成了# G4，如果将energyterm定义为"G4 Free Energy="（此处的双引号必须写），则molclus就会从输出文件里面找这个字段来读取G4方法算的自由能。
对于freq (itask=2)和opt freq (itask=3)任务不需要设这个，因为总是读取的直接输出的自由能（但如果用提供了template_SP.gjf让molclus在这两种任务后自动算高精度单点来得到较准确自由能，则需要设这个来指定从单点任务中读取什么能量，见后文）。
ibkchk：和ibkout的用法一样，但这个控制的是如何备份.chk文件
----以下是ORCA专用的参数----
orca_path：ORCA可执行文件的路径，比如"D:\study\ORCA4\orca.exe"、"/sob/ORCA4/orca"
ibkgbw：和ibkout的用法一样，但这个控制的是如何备份.gbw文件
ibktrj：1=把每次优化后产生的后缀为_trj.xyz的文件改名为orca[帧号].xyz备份出来。这是优化过程的轨迹文件，因此可以用VMD直接观看优化过程
ibkhess：1=把每次振动分析或优化+振动分析后产生的.hess的文件改名为orca[帧号].hess备份出来。这是Hessian文件，以便用户之后以某些方式利用之
mpioption： ORCA并行运算时传递给MPI的额外选项。默认为none，即什么也不设。此选项有特殊的用处，例如在root下通过较新版本OpenMPI运行ORCA，不带--allow-run-as-root选项就运行不了，此时可以设定mpioption= --allow-run-as-root来解决此问题。如果参数里面有空格，应当写成比如mpioption= " '--bind-to core' "
----以下是MOPAC专用的参数----
mopac_path：MOPAC可执行文件的路径，比如"C:\Program Files\MOPAC\MOPAC2016.exe"、"/sob/MOPAC2016/MOPAC2016.exe"。
----以下是xtb专用的参数----（molclus总是通过"xtb"命令调用xtb，因此不用指定其可执行文件路径）
xtb_arg：xtb程序运行时额外接的参数，要用双引号扩住，比如"--chrg 1 --uhf 1"代表体系电荷为1，alpha电子比beta电子多1个。（虽然自行使用xtb做优化时应加上--opt，但不要写在此处，因为molclus会根据任务类型自动加上）
----以下是Open Babel专用的参数----（以正常方式安装好Open Babel后molclus即可直接调用，不用指定其可执行文件obabel的路径）
obabel_FF：Open Babel中使用的力场。建议用MMFF94，是最准确的有机分子力场，但支持的元素也仅限有机体系所涉及的。目前的Open Babel里还可以用UFF、ghemical、GAFF力场
obabel_param：Open Babel运行时额外加入的参数，比如--steps 2500代表几何优化做2500步。（注意，2.4版之前的Open Babel使用默认的共轭梯度法优化时有bug，应加上--sd改为最陡下降法优化）
----以下是Shermo专用的参数----
Shermo_path：Shermo程序（http://sobereva.com/soft/shermo）的可执行文件的路径，必须用Shermo >=2.3版。如果此路径设置正确，则molclus调用Gaussian或ORCA做振动分析或优化+振动分析后，会自动调用Shermo来计算自由能而不直接从Gaussian和ORCA输出文件里读取。Shermo在这方面计算上比Gaussian和ORCA自己的代码更为强大、灵活、稳健
T：Shermo程序计算自由能用的温度(K)
P：Shermo程序计算自由能用的压力(atm)
sclZPE：Shermo程序计算自由能时ZPE部分用的频率校正因子
ilowfreq：低频模式的处理方式。0=RRHO模型（和Gaussian相同，不推荐），1=将低于100波数的频率提升到100波数，2=用Grimme的准RRHO方法模型（强烈推荐！对柔性大体系比RRHO模型精确得多，详见Shermo程序原文里的介绍）
2.4 技巧&细节
xyz是一种格式非常简单的可用记录多帧结构的格式，可以直接用文本编辑器手动编辑，因此用户可以自由修改traj.xyz、cluster.xyz和isomers.xyz的内容，从而控制molclus处理哪些结构，或者手动提取特定结构、删除多余的结构等。molclus完全基于.xyz文件，这使得molclus非常灵活。不懂xyz格式的话可以参看《谈谈记录化学体系结构的xyz文件》（http://sobereva.com/477）。
在便宜级别下做优化和振动分析，而在高精度级别下算电子能量，二者相加来得到较准确的自由能是司空见惯的做法，看比如《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（http://sobereva.com/327）。Molclus的振动分析（itask=2）和优化+振动分析（itask=3）任务直接支持这种做法。具体做法是在当前目录下提供template_SP.gjf（template_SP.inp）文件，设置Gaussian（ORCA）高精度单点计算的关键词，格式同template.gjf（template_SP.inp）。之后，当优化或优化+振动分析任务结束后，程序发现当前目录下有template_SP.gjf（template_SP.inp）时就会用其中的信息构成单点任务文件并调用Gaussian（ORCA）计算，然后自动将新算出来的单点能与前面算的自由能热校正量相加作为最终能量输出。几个细节：
(a)如果template.gjf和template_SP.inp同时存在，优先考虑template_SP.inp
(b)opt、opt freq任务用的程序与自动算高精度单点的可以不同，比如iprog=1调用Gaussian时也可以通过在当前目录下提供template_SP.inp来令程序自动调用ORCA算高精度单点
(c)利用template_SP.gjf基于Gaussian自动算高精度单点时，Gaussian用户需要在settings.ini里对energyterm参数恰当设置，用于指定如何读取这个级别的单点任务输出文件里的能量
molclus运行意外原因中断的续算方式：比如traj.xyz里总共有1000个结构要算，算到第180个结构的时候没算完就意外中断了，此时就把settings.ini里的ngeom设为180-1000，把iappend设为1，然后重新启动molclus计算，程序就会把之前没算的第180~1000个结构算完，并且把新产生的结构续写到之前产生的isomers.xyz里。这样，最后得到的isomers.xyz里就有1000帧了，和一次性全跑完得到的精确一样。
molclus可以同时运行多个以试图节约时间。大家都知道支持并行的计算程序的计算耗时不会和调用的核数线性相关，随着核数增加并行效率必然会愈发降低。比如你的机子核数很多，有56核，traj.xyz里有1000个结构需要处理，你希望降低molclus处理这个traj.xyz的总耗时，那么你可以比如把molclus文件夹复制成四份，比如分别叫molc1、molc2、molc3、molc4，分别将这四个目录下的settings.ini里的ngeom设为1,250、251,500、501,750、751,1000。然后开启四个终端，分别进入这四个目录并且运行molclus。等都运行完了，使用这条命令把四个目录下产生的isomers.xyz合并在一起：
cat molc1\isomers.xyz molc2\isomers.xyz molc3\isomers.xyz molc4\isomers.xyz > isomers.xyz
之后用isostat对这样的isomers.xyz产生的统计结果，和只跑一个molclus但是一次性考虑所有1000个结构的情况是完全一样的，而总耗时能降低不少。类似地，如果你有很多台机子可用，也可以这样把一个任务分开在不同机子上计算，然后把不同机子上的isomers.xyz合并到一起。
注：对于分割运行的数目非常多的情况，手动操作麻烦，可以考虑使用ggdh开发的辅助工具《分解合并Molclus任务的小脚本splitMC》（http://bbs.keinsci.com/thread-14361-1-1.html）。
为了灵活，molclus也可以直接指定设置文件和初始结构文件的路径。例如molclus /RAS/Layer.ini代表将/RAS/Layer.ini作为设置文件。再比如molclus lock.ini /RAS/chu2.xyz代表用当前目录下的lock.ini的设定对/RAS/chu2.xyz里的结构进行计算。因此设置文件和初始结构文件并非只能是当前目录下的settings.ini和traj.xyz。
isostat运行时可以结合各种参数。比如执行isostat Rize.xyz -Nout 8 -Eout 2.4 -Edis 0.5 -Gdis 0.2 -T 350 -nt 32命令，代表将当前目录下的Rize.xyz作为输入文件，往cluster.xyz里导出的簇是能量最低的8个且相对能量还必须在2.4 kcal/mol以内。使用0.5 kcal/mol和0.2埃分别作为归簇的能量和几何偏差判据。计算Boltzmann分布使用350K。使用32个线程并行计算。凡是通过命令行指定的参数就不需要在交互式界面里再输入。忘记这些选项怎么写的话可以随时运行isostat -h查看。
有用户对isostat程序的运行机制感兴趣，这里具体交代一下。isostat使用能量和几何结构偏差作为判断两个结构是否相似的标准，同时小于两个阈值才算相似。几何结构在判断时是基于距离矩阵来比较的（因此对比前isostat并不会对两个结构做对齐(align)，也不要求原子间序号顺序一致）。isostat会循环输入的.xyz里的每个isomer。对于输入文件里第一个isomer，或者当前isomer和之前已有的所有簇都不相似，那么这个isomer将被作为一个新的簇，此簇的能量、结构也等同于这个isomer。若当前isomer与之前某个簇相似（能量和结构差异都同时小于自设的阈值），那么这个isomer就被认为归入了这个簇，因此这个簇的容量会+1；如果与此同时这个isomer的能量比这个簇的能量更低，那么这个isomer将被作为这个簇的代表，即使用这个isomer的能量和结构作为这个簇的能量和结构。等所有isomer都循环完之后，程序就会输出各个簇的绝对能量，相对于能量最低的簇的能量差(DE)，以及这个簇与与它结构最相近的另一个簇的原子间距偏差最大值(DGmin)。显然，最终得到的簇的数目<=isomer的数目。
两个isomer的结构差异在isostat中是这样算：根据isomer的原子坐标计算距离矩阵，将其非对角元转化为一维数组形式，并从小到大进行排序，这称为原子间距数组。两个isomer间的结构差异用二者的距离间距数组的差值数组的绝对值的最大值来衡量。
traj.xyz和settings.ini也不是必须叫这俩名字、必须放到当前目录下。如果molclus启动时发现当前目录下没有这俩文件，会直接提示用户输入，因此也可以命名为便于分辨的名字。
![]() ![]() ![]()
下面，本文将示例如何将GROMACS或AMBER分子动力学程序、molclus程序和量子化学程序相结合进行团簇构象搜索和分子构型搜索。如果你不打算利用这俩分子动力学程序产生初始结构（其实一般也没必要，因为过程较麻烦），那么可以把下面例子中的动力学部分直接略过，直接参考基于molclus做搜索的部分，之后再看以下文章分别了解怎么用gentor、genmer、xtb、Confab明显更便利地产生初始结构并结合molclus做搜索：
《gentor：扫描方式做分子构象搜索的便捷工具》（http://bbs.keinsci.com/thread-2388-1-1.html）
《genmer：生成团簇初始构型结合molclus做团簇结构搜索的超便捷工具》（http://bbs.keinsci.com/thread-2369-1-1.html）
《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）
《将Confab或Frog2与Molclus联用对有机体系做构象搜索》（http://bbs.keinsci.com/thread-20063-1-1.html）
如果你虽然想基于动力学程序产生初始结构，但是又不太会用GROMACS、AMBER等基于经典力场的分子动力学程序，那么建议用上面帖子里提到的通过xtb程序跑动力学，方便得多得多。
3 分子团簇构型搜索实例：(H2O)4
此例我们通过molclus结合GROMACS分子动力学程序以及MOPAC、Gaussian 09 D.01量子化学程序搜索(H2O)4团簇构型，这个团簇已经被很多文章研究过了。这里用这种很小的团簇进行示例仅是为了节约时间，对无论多大的团簇计算过程基本都是一样的。
GROMACS是免费、最主流、速度快且很好用的小分子溶液和生物大分子的分子动力学模拟程序，所以我们这里用GROMACS 4.6.7，编译方法参见《Gromacs 5.0与4.6.7编译方法》（http://sobereva.com/247），较新版本的见《GROMACS的安装方法（含全程视频演示）》（http://sobereva.com/457）。大家也可以用其它自己擅长的动力学程序，只要能产生VMD能支持的轨迹就行，因为我们要靠VMD来把轨迹转为molclus能读取的.xyz轨迹格式（当然啦，如果自己有其它办法能转成.xyz轨迹格式也行）。
动力学我们用NVT系综。对于当前体系实际上不用周期边界条件也行。不过对于较大团簇的话，当模拟温度较高时，可能模拟过程中外侧的某个分子跑着跑着就飞了，越飞越远，之后的轨迹就没意义了。而如果用了周期边界条件，它飞过边界后还会跑回来。
GROMACS模拟部分所有涉及到的文件都可以在这里下载，如果有的地方不清楚直接看里面的文件就明白了。
H20_4_MD.rar.rar
(55.96 KB, 下载次数 Times of downloads: 2412)
3.1 准备模拟体系
对H2O我们就用SPC/E水模型，它的拓扑文件是Gromacs自带的，就是gromacs/top/gromos54a7.ff/spce.itp。我们用gview之类工具画一个水分子，保存成H2O.pdb。由于spce.itp里氧原子名字为OW，两个氢原子名分别为HW1和HW2，所以我们也把pdb里的原子名改为与之对应的。
然后建立个文件夹用于存放模拟相关的文件。把H2O.pdb放进去，并建立一个空盒子文件emptybox.gro，盒子各边长都为1nm，内容如下：
运行此命令往这个小盒子里加入4个H2O
genbox -cp emptybox.gro -ci H2O.pdb -o H2O_box.gro -nmol 4
注：本文的Gromacs的操作是对于4.x版而言的。如果你是gromacs >= 5.0版用户，不用按上面步骤建立emptybox.gro，直接写gmx insert-molecules -box 1 1 1 -ci H2O.pdb -o H2O_box.gro -nmol 4就行了。而且后文的mdrun应改为gmx mdrun，grompp应该为gmx grompp。
然后运行
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
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
执行下列命令开始模拟
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设定每隔几帧保存一次。
特别注意：标准的.xyz文件对各个原子记录的是元素名，而不应当是原子名。但是，按照上述方式载入gro和轨迹文件，再转换出的traj.xyz里记录的却是原子名，如果你用文本编辑器打开traj.xyz就会看到里面每个水包含OW、HW1、HW2，这都是GROMACS模拟过程中水分子里的氧原子和两个氢原子的原子名。一般情况下，用molclus读入traj.xyz之前，应当手动把这些原子名替换为对应的元素名，这样才能100%确保molclus能正常工作、结果有意义。但是，由于周期表里没有叫做OW、HW1、HW2的元素，因此molclus在读取的时候会进一步试图通过首字母判断元素，因此molclus是可以恰当地将这三种原子判断为氧原子和氢原子的，因此接下来的分析也没有问题。但为了稳妥，建议最好还是先把traj.xyz里的原子名手动替换为元素名，否则元素判断错误时结果将明显和实际不符（比如如果有个碳原子的原子名叫CA，那么就会被molclus视为钙原子，那之后的任务肯定就没意义了）。
3.4 molclus+MOPAC做初步优化
这一节涉及到的molclus的相关的文件在此：
molclus_MOPAC2012.rar
(37.5 KB, 下载次数 Times of downloads: 1284)
虽说这个体系比较小，DFT优化一个结构不会很耗时，但是我们这里要优化101帧结构，直接用Gaussian做DFT优化还是颇耗时的。而且轨迹中的结构由于热运动，偏离极小点往往比较远，优化不仅需要的步数很多，还经常容易发生震荡，或者坐标系出现问题（Gaussian的一个老毛病）导致优化了半天最终却失败。一个比较机智的方法是先用MOPAC在半经验方法下预优化一下（或者用xtb程序预优化），瞬间就能把这么多结构都优化完，而且几乎不会报错。但是优化出的结构的质量还只能算是作为精确计算的初猜的水平，没法直接用，因此之后还得再用DFT来精炼一下。虽然Gaussian也支持很多半经验方法，其中PM6-D3和PM7描述弱相互作用也不赖，但速度比起MOPAC还是慢很多。本例用的PM6-DH+用于氢键计算在半经验的层面上算是最理想的选择之一。
在molclus文件夹下编写一个MOPAC输入文件的模板文件，文件名应当为template.mop，内容为
- PM6-DH+ precise
- molecule
- All coordinates are Cartesian
- [GEOMETRY]
把molclus目录下的settings.ini设定好，大部分不用改，只需要确认iprog和mopac_path是否设对，请根据2.3节的说明合理地设定。
把上一节得到的traj.xyz拷到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.38549952 a.u., new cluster
Energy of structure 2 : -0.39920546 a.u., new cluster
Energy of structure 3 : -0.39694736 a.u., new cluster
Energy of structure 4 : -0.39290484 a.u., new cluster
Energy of structure 5 : -0.39648044 a.u., new cluster
Energy of structure 6 : -0.39694738 a.u., attributed to cluster 3
Energy of structure 7 : -0.39290484 a.u., attributed to cluster 4
Energy of structure 8 : -0.39637520 a.u., attributed to cluster 5
...
然后看到排除了重复结构后，能量由低到高排序的结果。E是体系能量，DE是相对于能量最低的构型的相对能量，DGmin是当前构型与与它构型最相近的构型的几何偏差（此例这里没意义，因为刚才没启用几何阈值来区分构型）。Count代表isomers.xyz中根据能量和距离阈值判断与这个结构相同的结构共有多少。可见能量最低的3个结构在轨迹中出现几率占了主导。
The lowest energy is -0.39920546 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的相关的文件在此：
molclus_G09_M062X_6-31Gx.rar
(8.49 KB, 下载次数 Times of downloads: 1918)
为了节省时间，我们用便宜但对于当前体系来讲几何优化精度可以接受的M06-2X/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]
以这种方式写这两个模板文件是有一些考虑的。经过MOPAC预优化，初始结构偏离准确结构已经不算太远了，当基于template.gjf的默认优化设定进行优化时，如果40步还没收敛，很可能已经发生了震荡，继续优化下去可能白费功夫，因此用了maxcyc=40来设定步数上限（对此体系默认值比这个高）。在template2.gjf中，我们尝试改用其它优化方法GDIIS，并且还把步长上限设定在了很小的值，这对于解决优化末期的震荡往往很有效。nosymm关键词其实不需要写，但写了的好处在于可以令Gaussian输出文件中的结构在观看时不会突跃（这是因为默认情况下Gaussian总把结构调整到标准朝向），而且写了也不会多花时间，因为初始结构偏离对称性还是不小的，Gaussian自动判断不出对称性。
模板文件里用了int=ultrafine和opt=tight，这是因为M06-2X这类明尼苏达系列泛函对积分格点要求高，优化时用较高质量格点比较放心（从G16开始int=ultrafine已成默认，不需要再手动写）。用较严的收敛限是因为分子团簇势能面有时会很缓（尤其是范德华作用主导的情况），用默认收敛限的话构型和能量都可能不很精确，可能有零点几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.56099830 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
可见，之前MOPAC优化出的结构经过二次优化后又去掉了些重复的，目前有5个构型。它们的结构如下
这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个结构经检验确实都无虚频。
注：为了方便，可以直接利用molclus对每个结构做振动分析。也就是把模板文件里的opt有关的关键词去掉，写上freq，然后把settings.ini里的itask改为2。运行molclus时，每计算完一个结构后就会提示这个结构的自由能、自由能热校正量，以及虚频数目，并且写入isomers.xyz。最终的isomers.xyz中，在每个结构的注释行里，G是自由能，Gcorr是自由能热校正量，Nimag是虚频数。如果你希望得到更高精度的自由能，可以比如在Linux下使用grep Gcorr isomers.xyz |cut -c 29-38命令提取对应Gcorr数据的列，将它加到下一节得到的高精度单点能上就获得了高精度的自由能，然后可以根据《根据Boltzmann分布计算分子不同构象所占比例》（ http://sobereva.com/165）准确计算出这些构型在特定温度下的分布比例。
如果想直接让molclus调用Gaussian时直接做opt和freq，就把模板文件里写上opt和freq关键词，把itask设为3。这样得到的isomers.xyz里存的是优化后的结构，优化后结构下的G、Gcorr、Nimag也都体现在里面。
3.6 高精度单点能计算
现在只有5个结构了，直接手工一一处理也比较方便了。M06-2X/6-31G*的结构还说得过去，但算的能量只是定性合理，准确比较不同异构体的能量显然还是得用更高精度再计算一下单点能。
利用molclus也可以方便地对这5个结构做单点能计算。一般认为MP2算氢键比较理想（虽然MP2计算其它类型相互作用很不划算因此强烈不建议使用），故我们用Gaussian在MP2/aug-cc-pVTZ下做单点，template.gjf内容写为 - # MP2/aug-cc-pVTZ
- Template file
- 0 1
- [GEOMETRY]
把settings.ini里的itask改为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条氢键的那个现在竟成了能量最高的。
在《大体系弱相互作用计算的解决之道》（ http://sobereva.com/214）一文中笔者大力推荐在ORCA里使用双杂化泛函PWPB95-D3(BJ)结合质量很高的def2-QZVPP基组做计算，精度很好，而且由于RI技术使得计算速度很快。这里演示一下molclus+ORCA在这个级别下对这些构型做单点能计算。
我们编写ORCA的单点计算模板文件template.inp，内容如下（假设用4核，每核分配约3GB内存） - ! RI-PWPB95 D3 def2-QZVPP def2/JK def2-QZVPP/C rijk nopop pal4 grid4 tightscf noautostart
- %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关键词，并确认itask设为了0、iprog设为了3即可。
3.7 分子团簇搜索小结
通过(H2O)4这个例子，基本上把分子团簇构型搜索需要注意的各种问题都讨论了，其处理过程对于多种分子构成的混合团簇也完全适用。一定要注意不同分子团簇特点不同，构型搜索需要的参数、具体策略往往差异很大，没有唯一通用的处理办法，必须根据实际情况反复摸索、灵活变通、合理巧妙运用molclus才可能得到想要的结果。对于缺乏经验的人，一次就做成功的可能性比较有限。
在不同计算级别下优化，极小点构型、数目往往会有差异（特别是那些高能、不太稳定的构型）。而且优化参数、筛选策略也影响最终是否能找全极小点。建议以不同优化方式多试试，否则很容易漏掉高能构型。
顺带一提，只要随便逮几个分子，按照此节方法用molclus搜索一下它们组成的团簇构型，算算能量差，做个振动分析算下Boltzmann分布比率，然后再根据此帖《Multiwfn支持的弱相互作用的分析方法概览》（ http://sobereva.com/252）的方法做一通相互作用本质的分析，就很容易能出一篇文章，但切勿用此方法大量灌水！下图就是用帖子里提到的RDG方法绘制的图像，氢键（蓝色等值面）、范德华作用（绿色）和弱位阻作用（棕色）都生动直观地显示了出来，对于把握结构复杂的分子团簇中的相互作用很有帮助。
4 分子构象搜索实例：麻黄素（Ephedrine）
molclus的使用方法通过第3节的例子已经讲得很透彻了，将molclus用于分子构象搜索的过程是大同小异的。这里我们对麻黄素（Ephedrine）进行构象搜索。分子结构如下，可见它有许多可旋转的键，不做构象搜索很难确定哪种构象能量最低。
此例Amber模拟部分所有涉及到的文件都可以在这里下载
Ephedrine_MD.rar
(608.63 KB, 下载次数 Times of downloads: 1070)
4.1 动力学模拟
本例用Amber动力学程序来做动力学部分，用的是Amber12，Ambertools里的antechamber工具可以方便地生成基于GAFF力场的小分子在Amber中模拟所需要的文件。如果想用GROMACS来模拟也完全可以，用《几种生成有机分子GROMACS拓扑文件的工具》（ http://sobereva.com/266）提到的acpype基于antechamber产生GROMACS格式的GAFF全原子力场拓扑文件跑动力学即可。
PS：免费的Ambertools从14版开始已经包含了amber的主要的动力学模块sander了，可以不必花钱买amber了。如果不会装Amber或Ambertools，可参见《Amber14安装方法》（ http://sobereva.com/263）。
首先建立个文件夹，在里面再建立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
quit
就得到了拓扑文件EPH.prmtop和结构文件EPH.inpcrd。
为了让动力学能够充分采样麻黄素的势能面各个区域，而且由于不必像分子团簇那样担心跑散，我们用很高的温度模拟（1000K）。总共模拟1ns，步长2fs，每5000帧写入一次轨迹，即间隔10ps。周期边界条件对此体系是多余的因此这里不用。此例是做真空下的构象搜索，所以igb=0。相应的参数文件md1.in如下
- 1ns si<font>mul</font>ation 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帧组合成了动画，如下所示
在VMD main窗口里点当前项目，点右键选Save coordinates，格式选xyz，然后保存成traj.xyz。
4.3 molclus+MOPAC做初步优化
像第3节的例子一样，我们先用半经验方法进行初筛，可以将要考虑的结构数目降低近一个数量级。我们编写MOPAC的模板文件template.mop：
- PM6-DH+ precise
- molecule
- All coordinates are Cartesian
- [GEOMETRY]
确认settings.ini中的iprog设为了2，itask设为了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*对MOPAC优化后的结构进一步优化。虽然对于当前体系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=、itask为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
笔者把其中部分结构做成了动画，如下所示
我们用VMD看看cluster.xyz中的第2帧，思考一下为什么它能量几乎最低。从结构上可以发现，它的柔性链的结构很舒展，没什么位阻，而且羟基氢和氮的孤对电子离得比较近，可能有内氢键的特征。为了确认这一点，我们用Multiwfn+VMD作RDG图，如下所示
可见，羟基氢和氮的孤对电子区域之间的等值面上明显有一块蓝色，很明确地证明了存在内氢键，这是这个结构能量低的重要因素之一。
我们可以看看这24个结构中苯环的取代基的构象是怎么分布的。用4.2节的做法令苯部分进行align，然后在VMD的Graphics - Representations-Trajectory里的Draw Multiple Frames里输入0:24，就可以将这24个结构同时显示出来，如下所示
可见，这些构象中苯环部分结构没有任何差异，只是取代基部分构象发生很大变化，分布很广，几乎把所有可能性都覆盖了。这个例子证明利用molclus结合分子动力学和量子化学程序可以很好地也比较方便地搜索小分子构象。
如果需要更准确的能量从而更可靠地判断哪个构象是能量最低、获得更准确的能量差以获得可靠的构象分布比例，可以像(H2O)4那个例子一样用更高精度的方法对各个结构计算单点，比如M06-2X/ma-TZVP、B2PLYP-D3(BJ)/may-cc-pVTZ。此例中，在MOPAC粗略优化后，其实没必把所有结构都用B3LYP-D3(BJ)/6-31G*优化，只要把能量低的一批（比如10个）优化即可，此时只需要把settings.ini里的ngeom设为10，接下来molclus调用Gaussian对traj.xyz里的结构优化就只考虑前10个了。
5 总结
本文介绍了molclus程序的用法，通过(H2O)4分子团簇的构型搜索和麻黄素的构象搜索这两个例子充分证明了molclus的重要实用性。可能有人觉得本文介绍的过程不够傻瓜化，得具备基本的计算化学知识，没法像一些商业化的程序那样点点按钮就出结果。但实际上，那样看似傻瓜的程序其实非常难用，可控程度差，很不灵活，精度也烂，教学用或者粗浅地研究可能比较适合，而做为正经的研究目的则基本派不上用场。而本文示例的分子动力学+molclus+量子化学程序联用的构型/构象搜索方法，只要稍加揣摩，就能在日后的研究中得心应手，起到重要价值。当然，千万别忘了，molclus本身绝对没有局限于只能用经典力场的动力学程序产生初始构型/构象，因为还可以用molclus自带的gentor、genmer组件，或者Confab、Frog2构象生成工具，或者xtb跑动力学等其它别的方式产生traj.xyz，然后molclus就能自动对其中的结构通过计算程序进行优化、能量计算和筛选，比起本文例子这样做基于经典力场的分子动力学产生初始结构更方便、更好用得多得多，本文的读者一定别忘了看本文第二节末尾提到的相应的帖子。
| 
[复制链接 Copy URL]
发表于 Post on 2015-1-5 18:29:02
收藏 Add to favorites
发表于 Post on 2015-1-5 18:44:02
楼主 Author