Gromacs-初识
记录首次成功的动力学模拟，用的是随便找的BODIPY。
程序安装
暂时conda install -c conda-forge gromacs=2024.5，等编译了再说
从sobtop开始
多数教程都是用的gromacs本身的工具，但由于BODIPY里面氟和硼的存在，力场参数大寄，所以笔者首先琢磨sobtop产生输入的方式。
Step0：准备工作
如果模拟的是CHNO等常见原子，可以直接跳过这一步。而如果像笔者一样模拟的体系带诸如硼一类的罕见元素，你中奖了！
去文献里检索所需的原子的LJ势参数，把参数填进sobtop的LJ_param.dat里。例如：
1
2
;Atom type ;vdW ;ε(kcal/mol)
B_NNFF 1.772 0.095 ;Ref:10.1049/mnl.2009.0112
除此之外，还要修改assign_AT.dat以便指认原子类型：
1
2
3
4
$B_NNFF
nbond 4
element B
bond 1 N 1 N 1 F 1 F
Step1：拟合RESP电荷
首先用DFT进行结构优化和频率分析。示例:
1
# opt freq b3lyp/def2svp em=gd3bj
计算完毕后，进行formchk，然后把输出文件用gaussview保存一份mol2的出来。将multiwfn目录下examples/RESP/RESP_noopt.sh复制到当前目录，运行：
1
./RESP_noopt.sh IPR.mol2
计算完毕后得到IPR.chg，备用。
Step2：产生拓扑文件
运行sobtop：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
IPR.mol2 // 载入mol2文件
7 // 设置电荷
10 // 使用multiwfn导出的chg文件
IPR.chg // 上一步计算出的RESP电荷
0 // 返回
1 // 生成拓扑文件
3 // 指认GAFF原子类型，也可以用4，无所谓了，反正大头是自定义的
5 // 很明显，GAFF力场不认识硼。使用5指认自定义原子类型
0 // 下一步
3 // 尽量用GAFF参数，没有的基于hessian拟合
IPR.fchk // 输入hessian文件
[ENTER] // 在当前目录输出top文件
[ENTER] // 在当前目录输出itp文件
2 // 生成gro文件
[ENTER] // 在当前目录输出
完成后，当前目录下有三个文件，分别为：IPR.gro、IPR.itp、IPR.top。将这些文件转移到进行gromacs计算的目录，运行：
1
gmx editconf -f IPR.gro -o boxed_bodipy.gro -c -d 1.0 -bt cubic
产生的boxed_bodipy.gro就是最终的设置过盒子的gro文件了。
Step3：能量极小化
进行MD之前，需要消除局域高能量构象。此例用的是DFT优化过的分子，当然没必要，但可以用以测试。编写参数文件minim.mdp。这里提供一个示例，来自[GROMACS] 求助：能量最小化的mdp文件怎么写。：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
;define = -DFLEXIBLE
integrator = cg
nsteps = 10000
emtol = 100.0
emstep = 0.01
;
nstxout = 100
nstlog = 50
nstenergy = 50
;
pbc = xyz
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
DispCorr = EnerPres
;
constraints = none
运行：
1
2
gmx grompp -f minim.mdp -c boxed_bodipy.gro -p IPR.top -o em.tpr
gmx mdrun -v -deffnm em
由于就40个原子，光速结束。能量最小构象保存在em.gro中，可以使用VMD打开，理论上看着应当非常合理，绝对不会是像笔者第一次跑那样把吡咯环干碎了。
em.edr记录了能量。导出势能：
1
2
3
gmx energy -f em.edr -o energy
12 // 选择Potential导出势能
0 // 结束输入
输出文件为energy.xvg，忽略前面的注释，最后有势能随步数变化的数据。此例体系小，就俩点，点多可以绘图观看。
em.trr记录坐标。可以导出轨迹pdb文件：
1
gmx trjconv -f em.trr -s em.tpr -o em.pdb
也可以VMD直接观看：
1
vmd -f em.trr -s topol.tpr
Step4：平衡
有时需要平衡NVT与NPT
与能量极小化一样，这都是不必要的，纯属测试。
待写