计算化学公社
标题: 使用molclus程序做团簇构型搜索和分子构象搜索 [打印本页]
作者Author: sobereva 时间: 2015-1-5 18:29
标题: 使用molclus程序做团簇构型搜索和分子构象搜索
注：本文更新频繁，每次出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使用的流程框图：
(, 下载次数 Times of downloads: 497)
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）和哪些程序的组合。绿色代表支持，棕色代表不支持。绿色里如果有文字的话，说明需要在模板文件里体现相应的关键词
(, 下载次数 Times of downloads: 232)
----以下是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模拟部分所有涉及到的文件都可以在这里下载，如果有的地方不清楚直接看里面的文件就明白了。
(, 下载次数 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的相关的文件在此：
(, 下载次数 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的相关的文件在此：
(, 下载次数 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个构型。它们的结构如下
(, 下载次数 Times of downloads: 334)
这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方法绘制的图像，氢键（蓝色等值面）、范德华作用（绿色）和弱位阻作用（棕色）都生动直观地显示了出来，对于把握结构复杂的分子团簇中的相互作用很有帮助。
(, 下载次数 Times of downloads: 314)
4 分子构象搜索实例：麻黄素（Ephedrine）
molclus的使用方法通过第3节的例子已经讲得很透彻了，将molclus用于分子构象搜索的过程是大同小异的。这里我们对麻黄素（Ephedrine）进行构象搜索。分子结构如下，可见它有许多可旋转的键，不做构象搜索很难确定哪种构象能量最低。
(, 下载次数 Times of downloads: 335)
此例Amber模拟部分所有涉及到的文件都可以在这里下载
(, 下载次数 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帧组合成了动画，如下所示
(, 下载次数 Times of downloads: 323)
在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
笔者把其中部分结构做成了动画，如下所示
(, 下载次数 Times of downloads: 346)
我们用VMD看看cluster.xyz中的第2帧，思考一下为什么它能量几乎最低。从结构上可以发现，它的柔性链的结构很舒展，没什么位阻，而且羟基氢和氮的孤对电子离得比较近，可能有内氢键的特征。为了确认这一点，我们用Multiwfn+VMD作RDG图，如下所示
(, 下载次数 Times of downloads: 327)
可见，羟基氢和氮的孤对电子区域之间的等值面上明显有一块蓝色，很明确地证明了存在内氢键，这是这个结构能量低的重要因素之一。
我们可以看看这24个结构中苯环的取代基的构象是怎么分布的。用4.2节的做法令苯部分进行align，然后在VMD的Graphics - Representations-Trajectory里的Draw Multiple Frames里输入0:24，就可以将这24个结构同时显示出来，如下所示
(, 下载次数 Times of downloads: 343)
可见，这些构象中苯环部分结构没有任何差异，只是取代基部分构象发生很大变化，分布很广，几乎把所有可能性都覆盖了。这个例子证明利用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就能自动对其中的结构通过计算程序进行优化、能量计算和筛选，比起本文例子这样做基于经典力场的分子动力学产生初始结构更方便、更好用得多得多，本文的读者一定别忘了看本文第二节末尾提到的相应的帖子。
作者Author: diaolanxinyu 时间: 2015-1-5 18:44
感谢sobereva老师分享！
作者Author: chrinide 时间: 2015-1-5 20:25
赞一个
作者Author: xmseet 时间: 2015-1-5 21:51
顶，绝对的支持~
作者Author: fhh2626 时间: 2015-1-6 11:56
原理略有些简单啊。。用MD方法搜索团簇的效率还是比MC低了不少，这样很难得到团簇的global minima吧
还有一个问题就是MD所得到的结构用于量化的OPT的话很容易不收敛的，因为距离极小值的距离比较远
如果要进一步开发的话，建议放弃MD软件，用basin hopping进行数据采样，这样得到的结构是一个分子力学上的极小值，既可以节约采样的时间，也可以节约量化结构优化的时间
作者Author: sobereva 时间: 2015-1-6 12:28
我开发过MC搜索团簇的程序，对于搜索原子团簇适合，但这种方法明显不适合分子团簇（本文强调的是分子团簇），或者说效率比MD明显要低，而且收敛难度反倒远比基于MD大得多，毕竟MC过程根本没有很好地考虑原子间的作用，别提几何优化收敛了，由于随机移动原子造成的构型不合理，甚至一开始就可能SCF不收敛，甚至分子团簇产生了错误的结构（比如氢原子干脆优化后发生转移了）。
分子动力学得到分子团簇全局极小点没有任何问题，是分子团簇构型搜索有效的方法，更是众多文献讨论分子团簇（包括很大团簇）的标准方法。basin hopping、metadynamics之类采样方未必效率更高，而且有局限性。
“原理简单”不是缺点反倒是优点，用简单的原理，很好地解决实际需要解决的问题，这才是好的、适合大众使用的程序。
反之原理复杂，甚至引入一堆参数，需要使用者需要大量经验积累，颇为麻烦地才能使用来解决看似不复杂的问题，甚至还得改反复改代码，那样的程序我不推崇，不是大众向的，而只是少数专家们发paper用的。
作者Author: tjuchan 时间: 2015-1-6 14:45
那个settings.ini中的ipause参数=3表示总不暂停？，但是文件中写的是=0表示不暂停。好像是，，
作者Author: sobereva 时间: 2015-1-6 15:00
帖子里写错了，已更正。
其实无所谓，只要不是1和2就行，所以0和3效果都一样
作者Author: tjuchan 时间: 2015-1-6 16:44
综合性有点强，很难把几类软件串在一起。
作者Author: ruanyang 时间: 2015-1-6 18:08
Sob 老师强，赞一个！
作者Author: wei 时间: 2015-1-6 20:35
1）该软件的思路简单但巧妙，确能有效解决团簇构型和构象搜索问题。
2）MD我用过gromacs的比较老的版本，它的安装和上手还是有些麻烦；对于体系包含1：1两种分子（或三种）的比较复杂的体系，有时不太好建模。能否MD软件也能支持MS做的MD轨迹文件呢，它还是比较好上手和建模的；我记得作者曾经也写过个MS轨迹文件转换的程序，这应当不算太困难。
3）QM能支持gaussian、orca、mopac，这个基本满足各种精度的需要了。
作者Author: sobereva 时间: 2015-1-6 21:24
gromacs还好。而且也有第三方编译的windows版，上来直接就能用。至于上手，其实效仿文章的例子，很快就能弄会，参数文件也就改改模拟时间、保存频率、温度这些，没多少需要变动的。
建模的话强烈推荐packmol（http://www.ime.unicamp.br/~martinez/packmol/），想建立什么样的都能建立，十分方便灵活。
MS模拟的轨迹可以通过perl脚本导出成xyz（http://sobereva.com/143），直接就能给molclus使用，不需要对molclus做任何修改。但我觉得与其学MS，还不如就用gmx或amber，MS价格不菲。
作者Author: panger 时间: 2015-1-7 10:01
强就一个字
作者Author: 花非花 时间: 2015-1-7 18:08
SOB威武！这个软件对阴离子、阳离子构成的化合物是否适用？
作者Author: sobereva 时间: 2015-1-7 18:48
同样适用。其实这和molclus本身没直接关系。对于一个体系，只要动力学程序能模拟，量化程序能算，那么就都可以用molclus做搜索。
作者Author: tjuchan 时间: 2015-1-8 22:02
那个不是说/6-31G*不能描述弱相互作用吗！至少得加弥散函数吗。
作者Author: sobereva 时间: 2015-1-8 22:09
对于几何优化，并且是氢键这样静电作用为主的，在优化的时候为了省时间（也避免收敛困难），不加弥散可以接受。几何结构对于弥散函数的敏感度也并不高。但是算能量的时候就一定要加弥散了。
作者Author: tjuchan 时间: 2015-1-9 08:45
哦，那对于我要计算的类似于苯二聚体，这种色散主导的弱相互作用呢？正寻思这用MC替代例子中的MD，因为没接触过MD，什么top文件都不会写》。
作者Author: sobereva 时间: 2015-1-9 12:28
假设你用你的MC的程序也能产生出.xyz轨迹，那么也可以直接结合molclus用。
但是MC不比MD更合适，MC在移动原子后，可能分子就严重变形、散架了，或者有明显不合理的接触。如果人为添加一些约束、筛选条件避免这种情况，不仅麻烦，得到的结构以及产生结构的效率也并不比MD有优势。
不管什么体系，色散还是静电主导，只要是分子团簇，文中的过程都是适用的（但也不排除个别情况需要人为干预一下）。
作者Author: Aminus 时间: 2015-1-11 23:34
Sob出手必属精品
作者Author: tjuchan 时间: 2015-1-14 20:02
chk文件能保存吗？
作者Author: sobereva 时间: 2015-1-14 21:33
没有保存和备份chk文件的功能，因为molclus考查的只是能量。
对于完整版购买者，如果确实需要这个功能或者需要其它功能/改进，都可以根据需要加上。
作者Author: tjuchan 时间: 2015-1-15 10:11
请问一下，算能量使用你推荐的orca，如何添加BSSE呢？因为我相求的是两分子之间的相互作用能而不是说总体能量？谢谢
作者Author: sobereva 时间: 2015-1-15 14:47
ORCA没法像Gaussian那样自动划分片段做BSSE，而必须按照counterpoise方法的实现步骤，自行定义鬼原子来手动做。定义坐标时在原子名后面写上冒号则它就成为鬼原子，如O : 7.439917 6.726792 7.762120
顺便注意，若像文中那样用了D3校正且在def2-QZVP级别来计算，根本不需要考虑BSSE。因为此时BSSE不仅已经非常小了，而且原本的DFT-D3校正的参数在拟合的时候也是在这个级别下没考虑BSSE的情况下进行的。
作者Author: ChemiAndy 时间: 2015-1-16 07:54
fhh2626 got a point，不是原理过于简单，而是有限温度MD无法确保搜索到所有的low lying minima，因为有些minima之间可能有很高的能垒，低温穿不过去，而提高温度的话容易跑散，evaporation问题。
MC做构象搜索更有效率。实际上不必担心跑出完全不合理的构象。因为“不合理构象”必然是高能的，必然在筛选中去除。它们代表了能垒顶部的一些构象，去除掉完全没有问题。其次，MC可以灵活定义允许的运动方式，比如，只允许水分子以整体方式转动和平动，这可以进一步加速构象搜索的效率。希望在下一个版本中加以考虑。
作者Author: sobereva 时间: 2015-1-16 15:54
势垒很高翻不过去的话只需要多做几圈退火即可，跑散了再收回来，只需要写个简单脚本自动来做即可。或者直接改改动力学程序，把参考温度变化设成正弦曲线，每到最低点的时候自动取一次结构，照样是基于MD的，对于分子团簇就算有多大的势垒也能容易地解决。也不用怕蒸发飞了，毕竟有周期边界条件。而且也可以人为给各个分子在体系质心上加入个弱限制势来避免。
MC对于分子团簇有极为麻烦的问题，其一是难以保持分子结构。如果随机移动原子，那么分子可能就散了，优化后成键关系也都变了。如果随机移动、旋转分子（比如相对于分子质心），那么分子总会出现不合理的接触，如果说通过检测条件去掉这些不合理的碰撞，反复无数次地生成猜测的结构，那么最终得到的这些分子之间也会非常散，优化起来效率极低，或者根本就得不到能通过检测的结构（比如10个小分子的团簇，通过随机生成初始结构，想得到一个又不很凌乱、空隙很大，同时又没有不合理的接触这是极为困难的）。而只有分子动力学可以“自然而然”地生成既具有随机性，而且结构又相对紧凑的初猜。所以，对于分子团簇，MC通常只会比MD效率低。除非真是有那种明确知道有很高势垒的情况，比如聚合物链内旋转（有限温度下不容易转动），那么可以考虑MC和MD相结合的方法，只把涉及势垒最高的自由度上以一定间隔用MC方法变动一次。但对于目前考虑的小分子团簇而言，基本没有这个问题。
原子团簇，用MC是绝对没问题的，而且效果很好，但不适合分子团簇。
如果非不愿意用MD，那么可以考虑用packmol批量生成初始结构，既紧凑，而且又有随机性（每次运行时用不同的seed）。
另外，molclus的输入是.xyz文件，并没有局限于只能用MD。不管什么方法、什么程序、现成的还是自编的，只要生成的初猜能写成.xyz轨迹形式，molclus都可以处理。
作者Author: mei 时间: 2015-1-19 04:00
请问以上介绍的方法可以用于研究不同分子组成的分子簇的构型吗？如果可以，如何将不同的分子加入盒子里，以及如何得到拓扑文件？
作者Author: sobereva 时间: 2015-1-19 04:38
完全可以。
文中用genbox加入了4个水分子，比如再加入3个NH3，就再运行一次genbox，把插入的分子设成NH3分子结构文件并且输入结构选为上一步genbox得到的结构即可。
用packmol可以更确切地、精准地构建混合体系。
程序内建的力场没自带要研究的分子的拓扑文件的情况，对于amber就用antechamber就行了，gromacs用此文的方法：
几种生成有机分子Gromacs拓扑文件的工具
http://sobereva.com/266
作者Author: kexueshi 时间: 2015-1-22 11:46
正在学习，谢谢sob！确实只需要注册就能下载，没有级别限制啊！太感谢了。
作者Author: tjuchan 时间: 2015-1-24 17:24
首先，我想请问一下，文献(关于两分子作用力的文献，在优化阶段)，作者经常说"全优化"、“对所有可能的构型进行优化”(这里都是谈两分子)，这里的全优化是什么意思啊？所有的可能的构型是不是有点主观性(所有是不是有点绝对)？
其次，我这里采用了帖子的做法，发现经过MOPAC这一步，筛选下来的构型就一两个(尽管能量差选0.01kcal/mol)，相比文献得到的构型少很多，但是我得到构型基本上是最稳定的(和文献最稳定的构型相似)，文献中虽然有些构型能量差很小，但是确实是完全不同的构型，所以我对帖子中这句“则如果两个结构能量相差>**KJ/mol就会被认为是两个不同的结构”是不是有点欠考虑，因为就像很多文献看到的有些构型虽然完全不一样，但是能量确实差很小。
作者Author: sobereva 时间: 2015-1-24 21:29
“全优化”并没有确切意思，很含糊。可能是指优化时没有外加约束。“所有可能的”也很含糊，只能是说初猜结构非常多，已经几乎穷尽所有情况。
可以把能量差判断阈值设更小。实际上我也打算在isostat里加入另一种判断方式，即几何结构偏差的判断，但最近还没时间写，过一段时间会加进去。
作者Author: gkobebryant 时间: 2015-1-26 15:03
good
作者Author: holgrove 时间: 2015-1-31 10:09
学习了
作者Author: 北纬18° 时间: 2015-2-3 15:50
可惜 molclus收费啊
作者Author: sobereva 时间: 2015-2-3 16:10
收费和CASTEP卖几十万比，真是九牛一毛
何况免费版能处理原子上限已经不少了
作者Author: sobereva 时间: 2015-2-9 14:08
为了解决你提到的这个问题，刚刚对isostat进行了更新并相应地发布了molclus 1.1版，见http://bbs.keinsci.com/forum.php?mod=viewthread&tid=765
这个改进很重要
作者Author: 北纬18° 时间: 2015-2-12 05:47
isostat 处理的应该是mopac 的输出的out文件，为什么给出的能量是 a.u. ？
作者Author: sobereva 时间: 2015-2-12 11:25
程序会自动转换
作者Author: 北纬18° 时间: 2015-2-13 01:07
EPH 的out文件中，TOTAL ENERGY 都是 -1899 ev 左右, 除以 27.2114 ev 也就 - 70 au 。
而例子给的能量 E 在 -0.057 au 左右， 不太明白是怎么算的。
作者Author: sobereva 时间: 2015-2-13 01:21
应当取FINAL HEAT OF FORMATION，这才是半经验方法真正有意义的输出。不要管TOTAL ENERGY。
作者Author: 北纬18° 时间: 2015-2-13 02:30
原来是以FINAL HEAT OF FORMATION为准。
作者Author: 北纬18° 时间: 2015-2-13 03:08
http://www.keinsci.com/research/molclus.html 主页还是 1.0 的，1.1的在哪？
作者Author: sobereva 时间: 2015-2-13 10:45
忘改链接了，现已更新
作者Author: therotyonth 时间: 2015-2-13 20:44
感谢sobereva老师分享！
作者Author: hzfish 时间: 2015-3-6 23:03
金属氧化物团簇吸附气体分子或者其它分子是否适用呢？
作者Author: sobereva 时间: 2015-3-7 07:23
molclus仅是个接口，能不能用就看两个问题：
1 是否有合适的程序能将此问题动过动力学（或其它方法，如蒙特卡罗）模拟出来得到轨迹，或者说产生一批分布不同的初始构型
2 MOPAC/Gaussian/ORCA能否计算
用MS或其它一些程序实现步骤1是没问题的。步骤2也肯定能做到，对这些量化程序来说属于常规计算。所以molclus可以用。
作者Author: hzfish 时间: 2015-3-10 09:42
谢谢sob的解答！！
作者Author: xxj9618 时间: 2015-4-1 21:23
sob老师，现在我已经用gromacs跑REMD得到了一系列轨迹文件，但是它是dt= 0.002；nsteps = 50000000 应该是 100 ns。由于我的分子是由28个氨基酸构成，并且轨迹文件相当大。我想请问一下，这个如果使用molclus的话，用Gaussian优化的可能性。
作者Author: sobereva 时间: 2015-4-2 10:27
你可以用MOPAC先算一个结构，看看耗时多少。然后从轨迹中均匀提取合适数目的帧，利用molclus调用MOPAC批量优化。
MOPAC是算得动的，计算量和取的帧数称正比。
作者Author: 北极天99 时间: 2015-4-23 17:35
看了好几遍，里面好多东西都看不懂。所以想问老师，我平时全部都是用高斯计算，如果我要使用molclus搜索构象，是不是至少还得再安装VMD和gromacs两个软件，还有其他的吗？
作者Author: sobereva 时间: 2015-4-23 18:48
其它的不用。需要装VMD和其它一种主流的动力学程序（amber、gmx、NAMD等等都可以。建议gmx，免费，照着上面的例子，即便没有使用基础，也很快就能搞自己的体系）
作者Author: 北极天99 时间: 2015-4-23 21:56
好的 谢谢sob老师
作者Author: 小范范1989 时间: 2015-4-25 11:00
[img]sob老师好，针对您的这个软件，我想咨询一下这么几个问题：
1：如果购买，是不是能开发票？我现在先想练练会不会用，熟悉之后购买，毕竟sob给出的价格很优惠。
2：如果说购买了您的那个能算更多分子的molclus。假如说我以后换了电脑，是不是还能安装？
3：我是高斯用户，经常来优化比较大的有机分子。是不是可以通过高斯view画出分子，然后使用molclus先预先优化，得到比较合理的结构文件，再放到linux的高斯优化？
4：我试了试您的3.5：molclus+高斯。遇到下面的问题
4.1：不明白这个您说的使用流程成的第一步，对我来说，我用高斯view画出分子，还需要先经过一个动力学的软件先跑？我这只有高斯。
4.2：我运行您的例子，包含了traj.xyz，还有template.gjf。就会出现下面图像的样子，自动打开了我的高斯。
路径我写为了gaussian_path= "D:\gauss-gaussview\install\g09w.exe" 这是我自己win电脑上安装高斯的路径。
现在我就卡在这，不知道再怎么办。谢谢sob的指点。
作者Author: sobereva 时间: 2015-4-25 21:34
1 可以开，但需要额外加收快递费。
2 能。直接拷过去即可
3 单纯的优化分子并不需要用moilclus。自行写一批优化任务的输入文件，批量执行即可。
4.1 需要先做动力学产生一大批初始构象，然后可以利用molclus进行反复筛选、排序
4.2 g09w.exe改为g09.exe。真正的高斯是g09.exe，g09w.exe只不过是它的一个图形外壳而已。
作者Author: zidu113 时间: 2015-4-27 19:50
sobereva老师：
我在win 7下，使用g09下运行你发给我的molclus程序（非免费win版）,我发现还是不能使用，生成的isomers.xyz是空文件。提示错误如下：
*** Configuration 50 ***
Loading geometry...
Checking geometry...
Generating gau.gjf...
Running Gaussian: "c:\g09w\g09.exe" gau.gjf gau.out
No executable for file l1.exe.
Search path GAUSS_EXEDIR is : " "
:no error
Error: Optimization did not converge!
Wall clock time elapsed in this cycle: 0s
麻烦您能指点一下，我可能是那里没还没有做对，谢谢。
作者Author: sobereva 时间: 2015-4-27 21:03
你得配置Gaussian的环境变量。参见Multiwfn手册附录1（直接贴在下面了）
(, 下载次数 Times of downloads: 188)
作者Author: zidu113 时间: 2015-4-27 21:41
，太谢谢sobereva老师了。问题圆满解决。明天再看看在ubuntu系统上有没有问题，还望老师能多多指点。
作者Author: zidu113 时间: 2015-4-29 19:09
sobereva老师：
你好.不好意思,又要麻烦你指点我,如何在ubuntu系统上运行你发给我的molclus程序（非免费win版).我按照网上的操作,试图安装g09,设置环境变量及文件权限，一直没有成功。恳请老师能否给你个较详细的g09安装步骤教程。
作者Author: sobereva 时间: 2015-4-29 20:01
解压到/sob下
建立/sob/g09/scratch
在/root/.bashrc里加入
export g09root=/sob
export GAUSS_SCRDIR=/sob/g09/scratch
source /sob/g09/bsd/g09.profile
在g09/Default.Route里面设定默认参数，如
-M- 3000MB
-P- 4
如果运行时提示
files in the gaussian directory are world accessible. this must be fixed
则在Gaussian目录下运行chmod 750 -R *
作者Author: zidu113 时间: 2015-4-29 21:42
谢谢sob老师每次都能这么及时回复.明天早上去实验室按您的指导再设置一下,相信这次应该没问题了.
作者Author: nkallwar 时间: 2015-6-12 15:46
此处发现失效链节一枚，辛苦sob老师更改一下哈
作者Author: sobereva 时间: 2015-6-12 17:05
谢谢，已改
作者Author: fzmn 时间: 2015-7-28 11:02
赞一个，感谢SOB老师的分享！
作者Author: spearous 时间: 2015-8-14 00:15
个人感觉，还是要看问题的势能面形状。如果是势能面是不少深度类似的势阱（如团簇），MC比较好；如果是由深浅相差较大的势阱组成，MD多次退火比较好。不好一概而论。其他算法没用过，不太清楚。也许有普适算法？
作者Author: sobereva 时间: 2015-8-15 02:44
原子团簇和分子团簇的特征有相当大的差异，这已经在势能面结构的基本形状上有决定性的差别了。分子团簇搜索首先必须保持分子的完整性，且不让分子之间不合理接触严重，这在MD中是自发的。而至于MD、MC在探索势能面的效率上，则又是另一个问题。
作者Author: liyuanhe211 时间: 2016-3-1 21:50
通过动力学方式获得分子构象搜所的初始结构时，直接（类似随机的）截取动力学中的一些点做初猜是否合适？
可否用Amber做到以这些随机点为初始结构，缓慢降温至0K后再用半经验、DFT优化？若研究对象是较复杂的分子（如近100个原子的有机体系），是否有必要？
作者Author: sobereva 时间: 2016-3-2 03:20
随机截取并不比均匀截取更好，因为没有均匀截取的系统性，以及保证涵盖的全面性。
amber用你那样的做法虽然可以，但不很理想。直接取结构的话，其实并不能说温度，除非把粒子动能也读进去作为初始条件。这样降为0K，每个点都需要进行MD模拟，花时间，不比直接优化好多少，如果是先用分子力场预优化一下倒是可以。100个原子这种情况按照文中那样普通MD就可以，如果怀疑有一些自由能垒很高妨碍转化，可以在gromacs里做周期性退火，然后molclus取其中对应恰降到0K的结构，会很方便也很有效。
作者Author: liyuanhe211 时间: 2016-3-2 10:29
我明白了，那就直接Amber了。十分感谢！
作者Author: 止水 时间: 2016-4-4 21:31
赞一个！！1
作者Author: lidanhoney 时间: 2016-5-30 19:48
老师，用高斯来优化3.5中的例子，是不是要很长时间，一开始运行，就没反应了。
作者Author: sobereva 时间: 2016-5-31 00:38
自行监控gau.out看看算到什么情况了
作者Author: yaochuang 时间: 2016-6-13 09:17
Sob老师，
如果想用molclus计算的结果发表文章的话，应该怎样描述程序或是引用什么文献呢？
作者Author: sobereva 时间: 2016-6-13 09:35
引Tian Lu, Molclus program, website: http://www.keinsci.com/research/molclus.html
作者Author: yaochuang 时间: 2016-6-13 10:05
好的，谢谢Sob老师！
作者Author: 鸢翔flybird 时间: 2016-6-13 16:27
本帖最后由 鸢翔flybird 于 2016-6-13 17:06 编辑
Sob老师您好，有几个问题需要请教：
1. 动力学搜索构象容易漏掉极小点吗？2. 看到您在《gentor：扫描方式做分子构象搜索的便捷工具》写道“对于体系较大的，可以用Gromacs周期性退火，每100 ps内从0 K升温到500 K再返回到0 K， 跑50 ns，得到500个处于0 K的构象，用molclus对它们优化”，请问这样只是得到能量最低结构的概率变大了，还是各个极小点不遗漏的概率也变大了？
3. 本帖中关于分子构象的教程是在非溶剂化的条件下跑动力学的，对于极个别时候气相和溶剂下分子结构差异大的情况，如果最终目的是想要精确地计算这种分子在水中的自由能，那么在气相下跑动力学+仅仅在量化优化时加smd溶剂化模型，这样得到的结果是不是不可靠？而应该在水盒子里跑动力学，而后溶剂化下量化优化？
4. 如果分子本身带电荷，跑动力学时是否一定要加离子抗衡？如果是的话，量化优化及计算单点能时是否需要将加入的抗衡离子去除？
5. 自己设计的新有机分子（100个原子以内）用Gaff力场是否够精确？还是要重新拟合参数？
6. 用Linux版molclus，在我们组的集群上运行高斯优化那一步时，出现了图1所示错误，好像环境变量有问题，高斯程序属性如图2，程序所在文件夹里l1.exe也在，请问如何解决呢？
作者Author: sobereva 时间: 2016-6-13 17:39
1 看你跑的动力学时间长短，以及温度等因素。除了势能面扫描以外没有任何构象搜索方法能确保找到所有极小点。
2 主要是增加找到最小点的几率，虽然极小点找全几率也会变大，但这不是主要目的，高能的极小点还是容易漏掉。
3 即便是溶剂下的情况，气相跑动力学也没问题，毕竟主要目的只是比较随机地产生初猜。动力学时候就要考虑溶剂效应的话建议用GB隐式溶剂模型。
4 不用PBC的话没必要加抗衡离子。用PBC的话应该加上，但量化计算时候可以去掉。
5 gaff一般可以。想更准确可以用MM3、MMFF94之类。
6 你先确认能否直接在命令行里用g09命令调用高斯，能的话，settings.ini里面高斯的路径就直接填g09就完了。如果直接用g09调用不了，说明还没配置好高斯。
作者Author: 鸢翔flybird 时间: 2016-6-13 20:12
谢谢Sob老师解答。关于第6个问题，我不太理解“直接在命令行里用g09命令调用高斯”的意思，是在putty等软件的命令窗口里“module load gaussian/09-d01”吗？我们组做高斯任务都是用脚本pbs提交到计算节点的，只有做类似于chk转为fchk这种小任务时才会module load一下，以便在提交任务的节点上直接完成。这样的话怎么用molclus调用高斯呢？
作者Author: sobereva 时间: 2016-6-13 23:12
pbs提交方式运行我不清楚
作者Author: therotyonth 时间: 2016-6-21 20:10
感谢sobereva老师分享！
作者Author: bianwenbo 时间: 2016-7-1 06:51
在Material Studio中有一个 Conform Analysis 同样是进行构象搜索的，不过那个似乎是基于分子力场的计算，可以对比较大的分子进行初步筛选，排除掉不合理的构象。对于包含非常多可自由旋转的长链烷烃类可能是较好的选择。不知道理解的对不对
作者Author: liyuanhe211 时间: 2016-7-18 19:38
用Amber做多个升温-降温循环的MD、产生柔性较高的复杂分子的初猜。在做MD之前需要用Amber力场做极小化吗（事先用MMFF94优化过）？
作者Author: 清岚 时间: 2016-7-20 10:34
老师，我在用您给的离子跑molclus程序是出现一下问题：在settings.ini中我设置的是igaucontinue是1，但是在服务器中运行molclus后的结果是attach://5270.jpg，好像并没有执行，用第一个模板优化没有收敛就直接结束了，请问老师这跟我写的脚本有关系吗？您能不能帮忙看一下脚本有没有问题attach://5271.jpg？
作者Author: sobereva 时间: 2016-7-20 13:09
不用，MD跑起来之后，一开始优不优化都一样。但如果MD一开始就崩了，还是得先优化一下。
作者Author: sobereva 时间: 2016-7-20 13:15
settings.ini传上来我看一眼
要确保改过那个参数后settings.ini已经保存了
作者Author: liyuanhe211 时间: 2016-7-20 14:18
我明白了，谢谢解答！
作者Author: 清岚 时间: 2016-7-20 17:07
本帖最后由 清岚 于 2016-7-20 17:13 编辑
恩恩 谢谢老师，我又检查了一遍，确定保存了，然后又试了几遍发现可以了。
我还有一个问题想问，我用您给的template.gjf模板跑的6个Li的团簇，traj.xyz文件给出的初始构型数是15种，但是没有收敛的，都是error，这样的话能不能试着改一下template.gjf中的关键词什么的，我看了”量子化学计算中帮助几何优化收敛的常用方法“这个帖子，可是不知道该怎么修改。这是我的settings.ini文件：attach://5274.jpg
作者Author: sobereva 时间: 2016-7-21 00:14
怎么修改你看molclus自动备份出来的每个优化任务产生的.out文件，看轨迹、能量怎么变化的，以此确定如何解决，比如能量还在下降可以把maxcyc加大。这种Li6团簇还是比较容易优化的
作者Author: 清岚 时间: 2016-7-21 08:55
好的，非常感谢老师，我会好好看的
作者Author: hopedream 时间: 2016-9-29 15:50
写的真好
作者Author: whtu 时间: 2016-10-11 11:00
sob老师，在服务器上执行molclus后只能根据traj.xyz文件逐个地进行计算，而且是没有提交到节点上进行计算的。molclus能不能实现提交到服务器上计算，在多个节点的情况下就可以实现多个任务同时计算，充分利用资源以及节约等待时间。
作者Author: sobereva 时间: 2016-10-11 13:10
提交方式运行应该可以，我没试过，但应该肯定有办法
作者Author: whtu 时间: 2016-10-17 11:03
1、sob老师，可以实现提交到节点计算，但是我还是没有办法实现在节点上同时执行多个结构优化计算，如4.4 molclus+Gaussian中gau00001.gjf、gau00002.gjf……同时优化，不知道sob老师有什么方法？
2、在分子构象搜索中，如果做的是真空下的搜索，那么MD部分对力场的选择是不是基本没什么要求了，选什么样的力场也就无所谓了？
3、如果做的是离子液体的阳离子或阴离子的构象搜索，如：1-(2-hydroxyethyl)-3-methylimidazolium 阳离子，在MD部分和molclus+MOPAC部分怎么体现或设置其电荷为+1？
作者Author: sobereva 时间: 2016-10-18 04:27
1 把molclus目录多拷几份，用不同文件夹命名试试
2 看你具体以什么方式做构象搜索，如果是分子力场做MD产生初猜结构，分子力场显然还是有可能影响搜索结果的，至少力场得让MD过程中能采样到你感兴趣的结构附近。
3 MD时候就通过原子电荷来体现。MOPAC时候你只需要设整体的电荷，量化计算SCF收敛后的波函数自动对应阴阳离子状态。
作者Author: jiameiye 时间: 2017-4-25 23:31
请问老师有提交molclus任务的PBS脚本模板吗？
作者Author: jiameiye 时间: 2017-4-25 23:50
#!/bin/bash
#PBS -l nodes=1:ppn=20
#PBS -l walltime=50:00:00
#PBS -N Molclus
#PBS -o $input.stdout
#PBS -q qdef
#PBS -A lp_leuven_mg
. /user/leuven/319/vsc31921/molclus_1.3.5_Linux_free/gau09e.vars
if [ ! -d $VSC_SCRATCH_NODE/$input ]
then
mkdir -p $VSC_SCRATCH_NODE/$input
fi
export GAUSS_SCRDIR=$VSC_SCRATCH_NODE/$input
#Specify molclus directory
MOLCLUS_DIR=/user/leuven/319/vsc31921/molclus_1.3.5_Linux_free
export workdir=$PBS_O_WORKDIR
#cd workdir and execute molclus
cd $workdir
./molclus >
rm -r $GAUSS_SCRDIR/*
请问老师，用这个PBS脚本运行molclus的语句是怎么写？
作者Author: sobereva 时间: 2017-4-26 05:53
PBS的事我不大清楚
作者Author: 冰释之川 时间: 2017-4-26 09:42
本帖最后由 冰释之川 于 2017-4-26 09:44 编辑
traj.xyz拷到molclus所在目录下：如果是Linux版，进入molclus文件夹，输入./molclus。然后molclus就开始运行了。
作者Author: jiameiye 时间: 2017-4-26 15:16
谢谢你，直接当前节点运行的我会，我说的是PBS脚本提交任务到任务排队系统，不是直接运行的。
作者Author: jiameiye 时间: 2017-4-26 15:17
好的，谢谢老师。
作者Author: jiameiye 时间: 2017-4-26 15:33
我这里只能进入登录节点，任务需要通过作业调度系统排队提交到计算节点，所以我需要一个PBS作业脚本提交运行molclus任务，脚本里会说明调用高斯来做构型搜索和优化。现在的问题是不知道settings.ini 怎么调用，最终输出文件是高斯的.log， 还是分析总结的结果文件。
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |