构建及优化[小大]分子力场参数脚本ztop.py

ggdh

I.脚本特色
1. 调用Multiwfn，amber等产生top文件(目前支持产生gaff2力场，输出格式支持amber，gromacs, charmm，dl_poly格式）（具体用法在本贴）
2. 根据量化结果，优化力场参数（目前支持平衡键长键角以及对应的力常数，电荷），以及补充缺失参数（具体用法在本贴）
3. 通过片段组合法，从小分子top文件制作大分子top文件（具体用法在这里）
4. 通过结构匹配，从小分子top文件中将部分参数转移到大分子top文件中（具体用法在这里）
5. 自动化二面角参数拟合（具体用法在这里）
6. 以上所有操作既可以单步独立完成，又可以联合多步，通常一行命令即可

待实装功能
下面是在计划中一些主要待实装的功能，都做出来，不知道还要鸽多久。就先把目前版本放出来再慢慢升级。
1. 调用libpargen产生opls-aa力场，调用openff支持SMIRNOFF力场
2. 调用intermol支持lammps和Desmond格式
3. 调用chargemol支持DDEC电荷
4. 支持生成全元素的LJ参数，从而支持全元素的top文件生成
5. 产生聚合物时增加波动和随机性
6. 产生周期性结构的top文件, 周期性体系的扩胞操作
7. 产生整个体系的top文件

反馈：兄弟们有好的功能建议或者是急需某待开发功能的也欢迎回帖或直接联系我（QQ：32589927）
功能一多bug就多，特别是搞拓扑文件会遇到各种奇奇怪怪的情况，作者无法都考虑到，遇到bug记得反馈！

脚本核心价值观：
解决各种top文件的制作问题，以及，超懒

II.安装依赖环境
安装依赖环境是难点
操作系统：windows 和 linux 都可以运行本脚本，但是需要调用Linux程序（AmberTools, Libpargen)的情况下，只能在Linux下运行
需要python3库：networkx(图论), openmm(动力学), parmed(力场文件io和转换)，pandas和numpy(数据处理), matplotlib(绘图)
需要的第三方软件：Multiwfn，AmberTools，Gaussian，Libpargen(暂时不用）
下面是具体的安装步骤：
1.首先安装最新版的anaconda3，设置好环境变量后，输入

1. which python

复制代码

确认当前python来自anaconda

2. 安装Multiwfn参考：Multiwfn的安装

3. 安装ambertools（需要使用其中的antechamber，parmchk2，tleap）参考：Amber14安装方法
在集群上普通用户使用gnu编译器编译AmberTools19方法, 或者直接使用下面conda安装

4. 用conda安装的时候，容易出现不同包之间的版本冲突问题，建议创建一个虚拟环境专门运行本脚本，下面的命令建立了一个名为ztop的虚拟环境，并在其中安装必要的包

1. conda update --all
   
2. conda create -n ztop -c conda-forge  networkx pandas numpy matplotlib parmed ambertools openmm

复制代码

之后要运行本脚本之前先运行下面的命令激活该环境

1. conda activate ztop

复制代码

运行下面的命令取消激活

1. conda deactivate

复制代码

安装中可能遇到的问题以及解决方法：
1. conda速度太慢，尝试使用conda和conda-forge的清华源(或者其他国内源）
2. networkx版本太低，导致。可以尝试运行时提示"No module named 'networkx.algorithms.graph_hashing'", 尝试运行

1. conda install networkx=2.5.1

复制代码

或者

1. pip install networkx

复制代码

3. --checkenv检测到openmm部分时提示"undefined symbol: omp_get_num_procs"， 尝试运行：

1. conda install --channel conda-forge llvm-openmp
   
2. conda update intel-openmp
   

复制代码

4. openmm的版本不够，openmm7.6和之前的版本不兼容，本程序目前支持7.6版本，所以如果你的openmm版本不够，ztop会认为你没有安装openmm，用下面的命令升级openmm到最新版

1. conda update -c conda-forge openmm

复制代码

安装本脚本
使用git的方式下载，因为本脚本会经常修复一些bug，会立马上传到git里面

1. git clone https://gitee.com/coordmagic/coordmagic.git

复制代码

得到coordmagic文件夹，如果以后有更新，cd到该文件夹下运行

1. git pull origin master

复制代码

进行更新
其中可以看到3个文件夹，分别是coordmagic，ztop，ztop_test，以及一个文件ztop.py
给ztop.py赋予可执行权限，然后将其所在文件夹加到PATH路径中，最后运行：

1. ztop.py --checkenv

复制代码

查看依赖的程序和python包的状态。

III.上手指南
下面按模块进行上手介绍，采用的例子文件在ztop_test目录下的PI_poly_init下，除了这里的介绍，大家还可以用：

1. ztop.py --help

复制代码

来查看选项的用法，下面的介绍先展示命令，然后对命令进行解释。
首先，切换到练习目录下

1. cd ztop_test

复制代码

下面的操作都将在这个目录下进行。

一、快速产生top文件(paramgen模块)
这个模块相当于一个懒人脚本，通过调用Multiwfn， AmberTools等程序，快速生成top文件。

例子1，直接从Gaussian结果文件制作gaff2-resp的gromacs的top文件
resp电荷的介绍见：RESP拟合静电势电荷的原理以及在Multiwfn中的计算

1. cp PI_poly_init/DON.log ./
   
2. cp PI_poly_init/DON.fchk ./

复制代码

DON.log 和DON.fchk 是Gaussian优化结构获得的log和fchk文件, 对应的gjf文件可查看PI_poly_init/DON.gjf

1. ztop.py -g "DON.log" -o DON.top,DON.gro

复制代码

-g选项开启产生top模式，-o指定输出文件的格式，.top和.gro为gromacs的格式，运行成功的话，当前目录下会产生gaff2-resp的DON.top DON.gro。

例子2，从mol2文件开始，制作gaff-am1bcc的amber的prmtop文件
这个就是模仿acpype.py了

1. cp PI_poly_init/DON.mol2 ./

复制代码

先准备好mol2文件, 这里用pdb，xyz以及Gaussian的log文件都可以。

1. ztop.py -g "DON.mol2;q=am1-bcc;ff=gaff" -o DON.prmtop,DON.inpcrd

复制代码

注意这里-g选项的用法，第一个是输入文件，然后用分号";"隔开后面的各种选项，q=am1-bcc指定电荷类型，ff=gaff指定力场类型。这里antechamber调用sqm算am1-bcc电荷需要等几分钟。

例子3，从mol2文件开始，让程序产生gjf文件，最终制作gaff2-resp2的charmm的psf文件

1. cp PI_poly_init/DON.mol2 ./
   
2. ztop.py -g "DON.mol2;q=resp2;ff=gaff2"

复制代码

因为调用Multiwfn去制作resp2电荷需要气相和液相的fchk文件，而目前没有，因此程序会报错退出，并且在当前目录下产生了DON_gas.gjf 和DON_sol.gjf两个文件, 用户需要打开这两个文件，修改其中的一些参数，比如溶剂，nproc，泛函，基组等。然后提交给Gaussian计算。Gaussian计算完成后，会把log和fchk文件放到当前目录下，继续运行

1. ztop.py -g "DON.mol2;q=resp2;ff=gaff2" -o DON.psf,DON.crd

复制代码

就可以顺利完成任务，注意这里如果让程序去产生Gaussian输入文件，会直接在log文件中输出静电势，从而省掉Multiwfn算esp时间。另外在生成gjf文件的时候可以用--gjf 选项控制gjf文件中的部分参数，比如：

1. ztop.py -g "DON.mol2;q=resp2;ff=gaff2" --gjf "nproc:20;mem:20GB;opt:0;solvent:chloroform;method:PBE1PBE;basis:def2TZVP

复制代码

这个选项大部分的意思熟悉Gaussian可以直接看出来，其中opt:0的意思是不要进行优化。
关于resp2电荷的介绍参考RESP2原子电荷的思想以及在Multiwfn中的计算，resp2默认是resp2-0.5，如果要产生resp2-0.6电荷，用下面的命令：

1. ztop.py -g "DON_gas.log;q=resp2-0.6;ff=gaff2"

复制代码

该命令虽然没有指定输出，但是仍然会在当前目录下产生 DON_resp-0.6.chg和DON_resp-0.6.mol2文件，这两个文件中都含有电荷信息。
如果只想调用Multiwfn产生电荷，而不用AmberTools产生力场（比如在window下），那么可以输入

1. ztop.py -g "DON.log;q=cm5;charge_only"

复制代码

这个会在当前目录下产生DON_cm5.chg，其中是cm5*1.2电荷。

二、从量化结果优化力场参数(paramrefine模块)
背景知识
小分子的力场参数通常包括键长，键角，二面角和相应的1,4作用，VDW参数，原子电荷。产生力场参数的方法见：几种生成有机分子GROMACS拓扑文件的工具。 这些工具的思路基本上是根据分子的拓扑结构确定原子类型。然后根据原子类型，在相应的数据库中找到键长，键角，二面角，VDW参数。然后通过量子化学软件计算得到RESP或DDEC电荷。此外SMIRNOFF力场（一种伏特加）则是绕过的原子类型的指定，直接根据分子拓扑结构直接从数据库文件中获得参数。无论是哪种方法，都是从预先储存好的数据文件中获得参数。这就容易带来误差：首先，原子类型种类有限，难以覆盖各种千变万化的实际结构；其次，原子类型的判断可能出现失误。我想既然已经用量化软件计算了电荷，为什么不顺便通过量化的结构，把力场参数也优化一番。因为最终的top文件，比如Gromacs的top，amber的prmtop中的力场参数，都是不依赖原子类型的（比如每一个键的平衡键长和力常数都可以不一样，即使他们对应的原子类型一样）。因此，本模块的制作的初衷就是直接把量化计算的结果，转换为力场参数，而不依赖原子类型的数量以及判断正确与否。有了此功能后，我们还能把激发态对应的键长键角电荷参数赋予分子，让MD能够在激发态势能面上跑一跑；也能够对一些没有原子类型的情况产生力场（目前由于产生LJ参数的功能尚未完成，支持全元素的计划还没实现）

例子4，优化top文件中的平衡键长和键角
首先运行例子1，得到DON.top和DON.gro文件，然后运行

1. ztop.py -p DON.top -x DON.gro -r e -o DON_rf.top

复制代码

-p 选项指定top文件， -x选项指定坐标文件，-r选项开启优化力场参数模式，-r e的意思是优化平衡键长和键角
此时程序会在当前目录找一个名为DON.log的文件（例子1中，已经把这个目录拷贝到当前目录下了），如果该log中含有内坐标信息（opt任务，或者单点任务中加上关键词Geom(redundant))，那么就会把该log中的键长键角设为平衡键长键角。
如果log文件名不是DON.log，而是比如XYZ.log，那么命令应该是

1. ztop.py -p DON.top -x DON.gro -r "e=xyz.log" -o DON_rf.top

复制代码

此外，还可以和例子1中的命令合并，从fchk和log一步得到优化键角的top文件

1. ztop.py -g "DON.log" -r e -o DON.top,DON.gro

复制代码

然后程序会给出下面的输出，另外，平衡键长键角的新旧数据对比保存在当前目录下的PARAM_REFINE下。


例子5，优化top文件中的平衡键长和键角和相应的力常数
力常数就是把Hessian矩阵投影到bond和angle上，采用的是modified Seminario方法，详见https://pubs.acs.org/doi/10.1021/acs.jctc.7b00785。
从log和fchk文件开始：

1. ztop.py -g DON.log -r "e;k" -o DON.top,DON.gro

复制代码

-r "e;k"意思就是同时优化键长键角和力常数，默认读取当前文件夹下的DON.log，和DON.fchk，注意这里的任务必须是含有freq的任务。
从top和坐标文件开始

1. ztop.py -p DON.top -x DON.log -r "e=DON.log;k=DON.fchk" -o DON.prmtop,DON.inpcrd

复制代码

其中-r “e=DON.log;k=DON.fchk” 指定平衡键长和力常数对应的log和fchk，其实这里用默认的-r "e;k"就可以，因为文件名和输入的top是一样的。然后程序会给出下面的输出，另外，新旧力场参数对比保存在当前目录下的PARAM_REFINE下。

讲一下modified Seminario(ms)和普通的Seminario的区别，这两个做键长是一样的，做键角时不同，以氨气为例，假如其中一个H-N在摆动，会同时造成2个H-N-H角的变化，那也就是说这两个键角力常数会同时阻止这个N-H的摆动，普通的Seminnario方法没有考虑这一点，他会把N-H的摆动分别投影到两个键角上，导致键角力常数被高估，于是modified Seminario方法搞出了一个scale factor去减少这个力常数，它们通过计算频率对少量分子做了验证，认为ms方法得到的频率更符合量化结果。最终的结果就是modified Seminario得到的键角力常数更小。但是那篇文章也不一定靠谱。所有如果大家想用原版的Seminario方法的话(原版的分子刚性更大，更不容易变形，加上关键词origin就好：

1. ztop.py -g DON.log -r "e;k;origin" -o DON.top,DON.gro

复制代码

三、从小分子top拼接产生大分子top（residuecomposer模块）
用ztop构建大分子的拓扑文件

四、二面角参数拟合（torsionfit模块）
用ztop拟合二面角参数

一些类似软件
https://github.com/alanwilter/acpype
yyds，一个单文件5000行的python代码，本软件的目的就是把他干掉，目前尚未成功。
QUBEKit
本贴程序的很多思路和功能(比如modseminario和待开发的vdw参数)就是借鉴的这个程序，不过这个不支持制作大分子（因为大分子量化算不动哈哈），而且最后的输出力场格式不太清楚（有知道的兄弟告知一声）
Force Field Toolkit Plugin
VMD的插件，产生charmm力场的charmm格式的拓扑。
QMDFF
产生的力场形式好像只有QMDFF本身能用，不支持gromacs，amber，openmm等高速MD，15年之后好像就没有继续开发了。
forcebalance
它的原理是首先设定目标值（包含能量，力，相互作用能，结合能，频率，密度，蒸发焓等), 然后给一系列结构，优化力场参数，使得力场获得的值和目标值尽可能接近。能够考虑到团簇或者bulk相关性质，相当于是拟合了所有的参数，这就要求样本量够大，这对于大分子来说计算量很大。其次如果一次拟合的参数过多，会让人担心因为参数之间依赖产生的过度拟合问题。

致谢：感谢sob提供的Multiwfn和技术支持，感觉NSFC NO. 51873160提供的经费支持。

lyj714

大佬nb

心向暖阳

千呼万唤始出来

牧生

不明觉厉，但又感觉可能以后用的着

wangyueda

期待二面角参数尽快上线！

sun666

请问大佬能否把amber的mcpb.py功能搬进去，以便对带金属的分子生成拓扑文件。这样功能更加强大

k64_cc

要不要考虑一下挪到Github上去，提issue应该比加QQ方便。

还有要不要考虑增加手动指定virtual site的功能？

ggdh

k64_cc 发表于 2021-8-3 10:55
要不要考虑一下挪到Github上去，提issue应该比加QQ方便。

还有要不要考虑增加手动指定virtual site的功 ...

gitee上也能提issue把，过两天全元素的功能实现了就有点特色了 再弄到github上，每次上github就要科学上网才有速度。。麻烦。。。
virtual site 我先看看怎么弄，看看parmed是否支持

k64_cc

ggdh 发表于 2021-8-3 11:16
gitee上也能提issue把，过两天全元素的功能实现了就有点特色了 再弄到github上，每次上github就要科学上 ...

能是能，主要是我没看着正文里有项目链接……

parmed肯定是不支持，不过GMX是支持的——但是OpenMM读GMX top又不一定支持，所以这事其实挺麻烦的。就目前结果而言，virtual site看起来好像是能拯救nonpolarizable forcefield的样子，尤其是在halogen bond和chalcogen bond上。

coolrainbow

赞啊，天下苦此久矣。每次有个有机分子要生成力场，就有无尽的工作要做。。

sun666

求问楼主Ni(NH3)62+ 这种离子配合物的力场该如何生成呐

心向暖阳

sun666 发表于 2021-8-13 18:05
求问楼主Ni(NH3)62+ 这种离子配合物的力场该如何生成呐

你把例子操作一遍 有了ztop 还要啥mcpb

sun666

心向暖阳 发表于 2021-8-13 19:25
你把例子操作一遍 有了ztop 还要啥mcpb

谢谢，我是在想能不能直接产生金属配合物的力场。 能的话造福人类了，最近愁这个呐

Zhangpc

请问出现python package simtk not available如何解决（应该是import simtk.testInstallation异常，我试了pip install simtk无果），谢谢！
ztop.py --checkenv
Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.
check parmed...parmed version 3.4.3 detected
check networkx...networkx version 2.6.2 detected
check Multiwfn...success!
check AmberTools...success!
check openmm...Error!!! python package simtk not available

飒飒萨达

谢谢老师，收藏了