gromacs_docker

在容器中使用gpu需要安装：nvidia-container-toolkit

运行命令：

docker run -it --net=host --gpus all --name 容器名 -e NVIDIA_DRIVER_CAPABILITIES=compute,utility -e NVIDIA_VISIBLE_DEVICES=all 镜像名

首次使用报错：

docker: Error response from daemon: could not select device driver “” with capabilities: [[gpu]].

解决办法：

sudo curl -s -L https://nvidia.github.io/nvidia-container-runtime/gpgkey | \
  sudo apt-key add -
distribution=$(. /etc/os-release;echo $ID$VERSION_ID)
sudo curl -s -L https://nvidia.github.io/nvidia-container-runtime/$distribution/nvidia-container-runtime.list | \
  sudo tee /etc/apt/sources.list.d/nvidia-container-runtime.list
sudo apt-get update

sudo apt-get install nvidia-container-runtime

02 验证执行下列命令：

which nvidia-container-runtime

输出 /usr/bin/nvidia-container-runtime，表示安装成功。 03 docker 使用:

docker run -it --gpus all **

此时，设置使用设备上全部的显卡。

目前尚不支持debian12

Dockerfile.lower

为一些老旧的cpu一些加速指令集用

Intel(R) Xeon(R) CPU E5-2650 v2 @ 2.60GHz

模拟教程

http://www.mdtutorials.com/gmx/

给小分子生成top学文件的工具

Sobtop（http://sobereva.com/soft/Sobtop） 这是我开发的GROMACS拓扑文件产生工具，主要产生GAFF、AMBER力场的拓扑文件，但由于其力场库可以自行非常方便地修改和扩充，因此Sobtop本质上是完全普适、通用的。Sobtop可谓是最理想、最灵活、最易用的产生GROMACS拓扑文件的工具。此程序用起来超级简单，什么额外的程序以及特殊的运行环境都不需要装，解压即用。Sobtop使用极其方便，照着屏幕上的提示敲几下键盘，itp、top和gro文件就产生了，另外也可以要求产生rtp文件。Sobtop的主页有非常详细的产生各类体系拓扑文件的例子，并给出了详细的相关要点的说明。从例子中你会体会到Sobtop的设计特别注重兼顾便利和灵活，初级用户会体会到它极其便利，而高级用户则会体会到通过Sobtop构建复杂体系拓扑文件特别灵活好用。

下载（访问码：npg8）

• acpype：这是一个Python脚本，可以在http://svn.code.sf.net/p/ccpn/code/branches/stable/ccpn/python/acpype/下载。使用前必须在机子里安装AmberTools（免费），acpype会调用其中的Antechamber先产生Amber格式的拓扑文件然后再转成GROMACS的。acpype用法很简单，要处理xxx.mol2就执行./acpype.py -i xxx.mol2，算完后会新产生一个xxx目录，里头有_GMX后缀的.gro、.itp、.top，直接在GROMACS里用即可。默认情况下，产生的拓扑文件是基于GAFF力场的，另外也会输出_OPLS后缀的基于OPLS力场的文件，但属于实验性质不建议用。.mol2文件用常用的GaussView就可以产生（但必须确保在gview里看到的分子结构中没有诡异的成键方式），也可以通过OpenBabel或Antechamber将其它格式转成.mol2。默认情况下acpype分配的原子电荷是Antechamber产生的AM1-BCC，虽然能用，但明显不如RESP/RESP2电荷理想。建议大家按前述做法用Multiwfn计算出RESP或RESP2电荷，自行写入分子拓扑信息的[atoms]的原子电荷那一列。

小分子top学文件生成

生成拓扑（topology）文件的流程涉及多个步骤，通常包括分配原子类型、确定分子内的相互作用参数以及生成最终的拓扑文件。以下是一个通用的流程总结：

1. 准备分子结构文件 输入文件: 提供分子结构文件，常见格式包括 .mol2、.pdb、.xyz 等。这些文件通常包含分子中各个原子的坐标和连接信息。
2. 分配原子类型 选择力场: 决定使用哪种力场（如 AMBER、GAFF、UFF 等）来分配原子类型。 分配原子类型: 使用自动工具根据分子的化学环境为每个原子分配适当的原子类型。 检查并手动修正未能自动分配的原子类型。 可能会使用如 assign_AT.dat 文件的自定义规则来辅助分配。 确认原子类型: 确保所有原子类型已正确分配。
3. 分配 Lennard-Jones (LJ) 参数 LJ_param.dat 文件: 检查并确保所使用的 LJ_param.dat 文件中包含所有已分配的原子类型的 LJ 参数。 检查参数: 如果工具提示找不到某些原子的 LJ 参数，手动添加或修改 LJ_param.dat 文件以包含缺失的参数。
4. 生成键合参数 选择生成方法: 直接使用预构建参数: 如果可能，使用预构建的键合参数（如键长、键角、二面角等）。 使用 mSeminario 方法: 如果需要更高精度的键合参数，可以选择基于 Hessian 矩阵的 mSeminario 方法。此方法需要外部量子化学软件生成的 Hessian 矩阵文件（如 Gaussian 的 .fchk 文件）。 混合方法: 结合预构建参数和 mSeminario 方法，缺失的参数可以通过猜测或其他方法补充。
5. 处理报错和修正 常见问题: 处理工具报错，如未找到特定原子类型的参数、文件格式错误等。 解决方案: 根据错误提示，修改相应的参数文件或重新分配原子类型，确保工具能够顺利生成拓扑文件。
6. 生成拓扑文件 确认生成文件: 当所有参数和原子类型都正确分配后，工具将生成最终的拓扑文件（通常是 .top 文件）。 输出文件检查: 检查生成的拓扑文件，确保所有参数和连接信息准确无误。
7. 后续处理 进一步优化: 如果有需要，可以使用其他工具或手动编辑拓扑文件以进行进一步优化。 使用拓扑文件: 将生成的拓扑文件与力场参数、坐标文件结合使用，开始分子动力学模拟或其他计算。

amber 立场原子查看：https://emleddin.github.io/comp-chem-website/AMBERguide-AMBER-atom-types.html

文件介绍

.itp 文件 内容: .itp 文件通常用于存储分子类型的详细拓扑信息（如配体、小分子、辅因子等）。它包括分子中的原子类型、键、角度、二面角、非键相互作用等，但不包含坐标信息。 用途: .itp 文件可以被包含在主拓扑文件（.top 文件）中，用于为特定的分子类型定义拓扑。这种方法有助于模块化和重用拓扑定义，尤其是在多个体系中需要使用相同的小分子拓扑时。 .top 文件 内容: .top 文件是整个系统的主拓扑文件，通常包含：

力场信息: 包括使用的力场类型和相关参数。 分子类型的定义: 引用或包含.itp文件以定义特定分子的拓扑。 体系中的分子数量: [ molecules ] 部分列出整个系统中包含的分子类型及其数量。 其他全局参数: 如组分的排除规则（[ exclusions ]）、组合规则（[ defaults ]）等。 用途: .top 文件是系统的主拓扑定义文件，GROMACS 在模拟过程中会根据该文件来了解系统中的所有分子及其相互作用。

小分子复合物合并

1.预处理受体蛋白、配体分子 //将受体蛋白另存为 protein.pdb //将配体分子进行合适的质子化后另存为 mol.mol2

2.为小分子生成top文件和gro文件 打开sobtop，把mol.mol2拖入sobtpo //选择1，生成top文件 //选择3，尽可能使用GAFF力场 //选择0，进入下一步 //选择4，如果可能，预先构建成键参数，任意猜测缺少的选项 //回车，生成的top文件生成在sobtop软件根目录下 //回车，生成的itp位置限制文件在sobtop软件根目录下 //选择2，生成gro文件 //回车，生成的gro文件在sobtop软件根目录下 //回车，退出sobtop软件

//将sobtop软件中的 mol.gro mol.itp mol.top //三个文件剪切到实验工作环境目录下

3.产生蛋白质连带一个离子的拓扑文件 gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top //选择AMBER99SB-ILDN力场和TIP3P水(注意，如果静电荷不为零: Total charge in system x.000 e 则需要添加抗衡离子) //得到的topol——Protein_chain_ A.itp对应蛋白，topol_Ion_chain_A2.itp对应离子,posre开头的文件对应限制势itp

4.合并gro文件，另存为complex.gro文件 //将mol.gro加入到pdb2gmx产生的protein.gro的末尾，并将第二行的原子数改为（蛋白质原子数+小分子原子数）,建议作图确认结构合理性，注意文件中的空格和回车问题 //注意，protein.gro文件末尾的三个小数是晶格的坐标，不要删除或大幅修改

//蛋白质的限制势itp文件在pdb2gmx的时候已经产生，但小分子的没有，genrestr是对输入的结构产生坐标或距离限制势itp文件的工具，接下来运行命令，进行限制势的产生 gmx genrestr -f mol.gro -o posre_mol.itp //选择组的时候选择0，system默认的位置限制势常数是1000kJ/mol/nm2，已经足够大 //将下列语句插入到mol.itp文件的末尾，注意，复制时连“#”井号一同复制，最好在末尾添加之前空一行，方便检查文件错误

#ifdef POSRES #include "posre_mol.itp" #endif

//这样当mdp中使用define = -DPOSRES的时候配体的位置也会被限制了

//把配体的itp文件引入整体的拓扑文件topol.top，把分子数也设置好 //注意，在引入的时候需要将小分子的mol.itp文件引入到蛋白质链之前，因为mol.itp最开头定义了[atomtypes]因此，这个itp要最优先被引入 //将下列语句插入到引入蛋白质的itp文件引入之前； #include "mol.itp" 并在末尾的[molecules]中引入mol 1，将topol.top的格式与complex.gro中分子出现的顺序对应 //即topol.top中应该是类似这样的顺序：

; Include forcefield parameters #include "amber99sb-ildn.ff/forcefield.itp"

; Include chain topologies #include "mol.itp" #include "topol_Protein_chain_A.itp" #include "topol_Ion_chain_A2.itp"

; Include water topology #include "amber99sb-ildn.ff/tip3p.itp"

#ifdef POSRES_WATER ; Position restraint for each water oxygen [ position_restraints ] ; i funct fcx fcy fcz 1 1 1000 1000 1000 #endif

; Include topology for ions #include "amber99sb-ildn.ff/ions.itp"

[ system ] ; Name CATIONIC TRYPSIN in water

[ molecules ] ; Compound #mols Protein_chain_A 1 Ion_chain_A2 1 mol 1

//注意，末尾的mol和1之间有一个空格，且下一步添加溶剂后有可能会缺失回车，需要纠正

5.设定盒子，加水，加离子，能量极小化 //设定盒子 gmx editconf -f complex.gro -o complex_box.gro -d 0.8 -bt cubic //设定的盒子是立方盒子，可能会增加计算量，但是不容易产生边界相互作用的问题

//加水 gmx solvate -cp complex_box.gro -o complex_SOL.gro -p topol.top //注意这一步加水后有可能topol.top文件最后一行的SOL可能会串行，需要手动添加回车，避免其与mol 1连在同一行，容易在后续处理中报错

//将mdp模板中的em.mdp,pr.mdp,md.mdp拷贝到当前目录 //产生临时tpr文件 gmx grompp -f em.mdp -c complex_SOL.gro -p topol.top -o em.tpr -maxwarn 1 //这里如果警告较多可以将maxwarn的数值改大一些

//加离子，使得整个体系变为电中性（与第三步中的电荷数相对应） gmx genion -s em.tpr -p topol.top -o system.gro -neutral //这里选择分组时选择SOL对应的分组 //产生的带有离子且电中性的体系为system.gro

//能量极小化 gmx grompp -f em.mdp -c system.gro -p topol.top -o em.tpr

gmx mdrun -v -deffnm em

6.限制性动力学100ps gmx grompp -f pr.mdp -c em.gro -p topol.top -r em.gro -o pr.tpr

gmx mdrun -v -deffnm pr //注意，对于复杂体系如果用常规2fs步长可能模拟刚开始就会崩溃，因此此处做限制性动力学期间，步长用较小的1fs以求稳妥。

7.产生索引文件 gmx make_ndx -f pr.gro 依次输入： x1|x2|x3（定义“蛋白-离子-配体”组，新的组号接在原有的组号之前，组号为R）（lijilo注释：|为右大括号右侧的按键，用“shift”+“、”便可输入该符号） !R(把其他部分定义为一个组，新的组号为R+1) name R protein_lig(改名为蛋白_配体) name R+1 envir（改名为环境） q //得到的index。ndx里的组名就和模板文件里面的md.mdp中的组对应了

8.常规动力学1ns gmx grompp -f md.mdp -c pr.gro -p topol.top -o md.tpr -n index.ndx -maxwarn 10 //如果这里警告多，可以将maxwarn后的数值改大一些

gmx mdrun -v -deffnm md //由于配体、离子与蛋白的相互作用较为明显，所以这里任务之中将他们一起作为一个控温组和同一个消除屏东转动的组 //如果需要跑较长时间的动力学，需要修改md.mdp文件

9.分析RMSD gmx rms -f md.xtc -s md.tpr -o rmsd_protein.xvg //两次都选骨架部分（protein） //注意观测在模拟过程中骨架的RMSD波动是否较大

gmx rms -f md.xtc -s md.tpr -o rmsd_lig.xvg //第一次选择骨架部分（protein）第二次选择配体（MOL） //可以消除蛋白质整体的运动，观察小分子配体相对于蛋白质的运动

由于sobtop有点问题

改用方法：参考博客

RESP 原子电荷生成

利用 ORCA+Multiwfn 生成小分子的 RESP 原子电荷（假设要计算FAD.mol2）：

默认axv2指令集，运行核心数目12

cd /root/workdir
./RESP_ORCA.sh mole.mol2 

ps: RESP_ORCA.sh 我进行了调整官方用的版本是界面版本，我用的是 Multiwfn_noGUI 可以参考源码仓库的 RESP_ORCA.sh 文件。

运行结束就可以得到一个 mole.chg 文件

小分子键的优化

利用 ORCA+Multiwfn，来生成小分子键的优化

这一部分参照了这个视频的流程和这篇文章，主要是为了得到 .hess 文件）：

(1)将 FAD.mol2 文件转成 ORCA 输入文件 .inp 格式， 具体操作为：打开 Multiwfn 并导入第一步中用过的 FAD.mol2 文件，依次输入：

我这里使用 Multiwfn_noGUI 版本

容器映射关系： - ./script/RESP_ORCA.sh:/root/script/RESP_ORCA.sh

Multiwfn_noGUI mole.mol2 # 交互式运行
Multiwfn_noGUI mole.mol2 > /dev/null << EOF
100
2
12
\n
-11
1
0
4
1
0
q
EOF

操作注释

100 # #选择其他功能的part 1
2 #文件转换功能
12 #生成 ORCA 输入文件
\n #输入你想保存的路径和名称，如果按回车就会保存在当前文件夹
-11 #选择 ORCA 版本
1 # 选择 ORCA 5.0 更高版本
0 #更改任务类型，因为我们要进行键的优化嘛
4 #优化
1 # B97-3c
0 # exit
q # exit

最后会得到一个 mole.inp 文件

(2)使用 ORCA 得到 .hess 文件；在.inp 文件（假设它叫 mole.inp）所在位置打开 cmd，输入：

需要使用绝对路径计算

/root/orca_6_0_0_shared_openmpi416_avx2/orca mole.inp > mole.out
# 不行就试试ORCAPATH\orca mole.inp > mole.out，其中ORCAPATH是你的orca的绝对路径）等待计算完成就可以啦。

(3) 使用sobtop生成topo

容器中已经安装在/root/sobtop

使用前将生成的.chg 文件和.hess 文件和.mol2 放到同一个文件夹中，执行：

操作前准备：

cp -r /root/sobtop/* ./
chmod +x ./sobtop
chmod +x ./atomtype
/root/workdir/sobtop mole.mol2

操作

/root/workdir/sobtop mole.mol2 > /dev/null << EOF
7
10
/root/workdir/mole.chg
0
2
/root/workdir/mole.gro
-1
4
1
2
2
/root/workdir/mole.hess
/root/workdir/mole.top
/root/workdir/mole.itp
0
EOF

操作注释

'.mol2 path' #键入原本的 .mol2 文件的绝对路径
7 #添加电荷
10 #添加由 Multiwfn 生成的 .chg 文件
'.chg path' #键入第一步中得到的 .chg 文件的绝对路径
0 #返回
2 #产生 .gro 文件
'your path' #你希望储存的路径和名称，按回车就会生成在默认（应该是sobtop.exe所在的）文件夹里
-1 #设置力场的产生方法
4 #选用 DRIH 方法
1 #产生 GROMACS 拓扑文件
2 #使用 GAFF 原子类型，没法识别的自动用 UFF 原子类型
2 #通过DRIH方法得到力常数
'.hess path' #键入第二步中得到的 .hess 文件的绝对路径
'' #.top的储存路径和名称，按回车就会生成在默认文件夹里（和sobtop.exe在同一个文件夹内）
'' #.itp的储存路径和名称，同上
#输出结束且无报错之后
0

然后，就得到了相应的.gro、.top和.itp文件了，接下来可以进行gromacs的分子动力学运算了。