用MOPAC做结构优化

MOPAC (Molecular Orbital PACkage) 是一款专注于半经验方法的量子化学程序。尽管高斯等主流的量子化学程序都支持半经验方法，但在计算速度上都无法与MOPAC相比。除了速度快之外，MOPAC支持的半经验方法也很全面，包括最新的PM6、PM6的各种弱相互作用校正版本和PM7等。MOPAC2016的安装教程见《MOPAC2016安装教程》一文。本文简单介绍一下用MOPAC做大分子体系的结构优化的方法。

输入文件

输入文件后缀为.mop，先看例子：

PM6-DH+ MOZYME
ALA
test
N     -0.73442210    0.95158026   -1.92049098
H     -0.00239630    1.61342378   -1.75893173
H     -1.49770492    1.40459869   -2.38111407
H     -0.39139292    0.20893021   -2.49564477
C     -1.19272716    0.40442879   -0.63539109
H     -1.55976838    1.19906435   -0.01997654
C     -2.32004689   -0.61481022   -0.88419234
C     -0.01727161   -0.29321959    0.07396847
H     -1.95300567   -1.40944578   -1.49960689
H     -2.65364309   -1.01307694    0.05122050
H     -3.13675951   -0.13008050   -1.37705905
O      1.15958265    0.11183605   -0.11181243
O     -0.26311155   -1.40132646    0.94377413

输入文件的第一行是关键词，指明方法、计算任务和参数等。本例中没有指定任务类型，MOPAC默认做结构优化，而不是单点计算。第二行和第三行内容任意，也可是空行。下面则是分子的坐标。MOPAC中支持使用笛卡尔坐标，也支持使用内坐标，还支持两者混写，具体可参阅在线手册。一般来说，我们用笛卡尔坐标较多。

输出文件

使用MOPAC2016.exe XXX.mop命令，即可运行该任务，输出文件的后缀为out。下面是该例的输出文件的一些重要部分：

首先是结构信息：

该部分给出了体系结构的一些信息，包括化学式、点群等。下面是程序根据一般的成键规则判断的体系中原子的带电情况，对本例的丙氨酸体系，笔者故意使用了内盐形式：

程序也将其判断出来了。这个小功能在我们平时的建模中其实挺实用，比如PDB文件或CIF文件中记录的结构常常有缺失、重复或氢的位置不合理，可尝试用MOPAC检查一下。笔者最近从一篇文献中拿到一个BN纳米管的结构，原先我并没有在意文献中的结构是否合理，直接拿过来计算了。后来用MOPAC做结构优化发现程序识别出结构带16个单位正荷，于是仔细检查了结构，发现原文的作者没有对端部的B用H原子饱和：

若未设定体系的电荷和多重度，MOPAC会自动判断。若要指定电荷，使用charge=n，而指定多重度则直接写singlet、doublet、triplet等或通过MS=n来设定。

接下来便是结构优化信息：

本例使用EF (EigenFollowing)算法优化结构，这是对原子数较少的体系的默认算法。而对变量数超过100的体系，默认使用L-BFGS算法。下面的CYCLE是指结构优化的一圈，而非SCF计算的一圈。后面会输出梯度模 (gradient norm) 和生成焓，默认当梯度模小于1时，结构优化收敛。

之后是收敛结构的一些信息：

半经验方法一般使用标准状况下的生成焓来拟合参数，因此输出的能量也是生成焓，这与一般的电子结构方法不同。

最后，CARTESIAN COORDINATES部分便是我们最关心的优化后的结构：

可以存成xyz文件用VMD等可视化程序查看。MOPAC官方没有比较好的可视化程序，卢天曾开发过一个小工具可用于读取每一步的坐标，并转化为xyz文件，用VMD查看轨迹。在使用时需要注意加上prnt=2关键词，这样MOPAC才会在每一步输出坐标。而使用L-BFGS时，即使使用prnt=2，也不会每一步都输出结构，此时可以加上EF关键词，强制使用EF算法优化。具体使用方法参见：

http://sobereva.com/212

结构优化的设置

MOPAC中可冻结部分原子，此时只要在该坐标后面写0即可，需要优化的坐标后写1（可省略），如

PM6-DH+ MOZYME
ALA
test
N     -0.73442210 0   0.95158026 0  -1.92049098 0
H     -0.00239630 1   1.61342378 1  -1.75893173 1
H     -1.49770492 1   1.40459869 1  -2.38111407 1
H     -0.39139292 1   0.20893021 1  -2.49564477 1
C     -1.19272716 0   0.40442879 0  -0.63539109 0
H     -1.55976838 1   1.19906435 1  -0.01997654 1
C     -2.32004689 0  -0.61481022 0  -0.88419234 0
C     -0.01727161 0  -0.29321959 0   0.07396847 0
H     -1.95300567 1  -1.40944578 1  -1.49960689 1
H     -2.65364309 1  -1.01307694 1   0.05122050 1
H     -3.13675951 1  -0.13008050 1  -1.37705905 1
O      1.15958265 0   0.11183605 0  -0.11181243 0
O     -0.26311155 0  -1.40132646 0   0.94377413 0

此时就只优化体系中的H原子。当然，若为了实现这个目的，在MOPAC中有比较简单的写法，即同时写上noopt opt-H两个关键词便可，表示不做优化但优化H原子。

方法的选择

对于大体系，尤其是弱相互作用体系，推荐使用PM7和PM6的带弱相互作用版本，如PM6-DH+、PM6-D3H4等。笔者曾经优化过一个约900原子的DNA螺旋结构，发现PM6-DH+得到的结构最好。当然，凭一个例子无法判断好差，反正MOPAC算起来飞快，几个方法同时试一下也无妨。

传统半经验方法的标度为O(N3)，MOZYME是一种使半经验方法实现线性标度的方法，一般都推荐用上。MOPAC官网有一个使用和不使用MOZYME的对比：

可以看出，使用MOZYME技术，能极大地提高速度和降低内存的使用，处理1000原子的体系都是轻轻松松的。

Comments

• 本文只介绍了MOPAC最基本的功能。事实上程序还支持很多功能，比如溶剂化效应、过渡态优化、IRC扫描、振动分析、激发态计算、周期性计算，还有很多性质的计算，读者可以参考官方的在线手册。
• MOPAC也支持优化过渡态，但笔者使用后发现体验不佳，结果也不一定可靠，本文也不做介绍了。
• MOPAC单核的运行速度已经非常快了。理论上也支持并行，且还支持GPU加速。笔者没有试过GPU加速，但是对于CPU并行，我实在是看不出并行和串行的快慢区别。
• 半经验方法归根到底不是一种高精度的方法，优化出的结构的可信度也不高。但作为DFT计算的预处理还是可以的。