适用于QM/MM模拟的一些技巧

引言

多尺度方法通过采用结合不同尺度的方法描述体系的各个部分使得模拟复杂体系中的物理化学过程成为可能，最具代表性的就是量子力学/分子力学(Quantum Mechanics/Molecular Mechanics, QM/MM)方法。对发生化学反应、光激发过程的核心区域（如酶活性中心，溶液中的溶质）采用量子力学方法，确保能准确描述电子结构的重排；对环境区域，则用分子力场，描述原子级别的运动。这类方法大大降低了计算耗时（不需要整个体系用量子力学方法），然而引入了建模的复杂性。

过去几年本人一直在与QM/MM的建模斗争。最初由于QM/MM建模的繁琐，曾用簇模型代替。然而酶催化反应总是动态的，需考虑采样问题，此时需要模拟完整的溶液中的酶体系；另一方面，截断的簇模型也难以考虑反应中心与环境间中长程的相互作用（如静电稳定化作用）。尽管后续学习到传承下来的方法，也常为手工操作之繁琐而困扰。单个体系的模拟尚且如此复杂，何况对大批量体系的高通量模拟。

本人也试用过多种软件来辅助这一过程，但仍未有十分满意的工作流。尽管功能全面的软件、平台不在少数，如Amber、ASH、PyChemshsell等，但每学习一个就相当于学习一套新的语言规范，况且还要处理安装依赖、环境之类的问题，这些平台自身也不能完成所有事情。本文并不讨论如何建立一个自洽的多尺度模拟工作流，也许将来会尝试做这件事。仅仅简单记录QM/MM模拟过程的一些技巧，为自动化多尺度模拟铺路。

用PyMOL选择QM区域

进行QM/MM计算，常常先进行分子动力学模拟进行构型采样，然后选取活性中心几何结构较有利于反应的构型进行后续QM/MM能量最小化和势能面扫描；或将轨迹聚类分析得到的每个类别代表性帧作为后续自由能曲线计算的初始构型。无论如何，都需要先定义QM区域。除了参与反应的底物必须包含在QM区域，和底物形成盐桥或强氢键的残基也应包含于其中；QM-MM边界的化学键最好不是双键，一般是非极性共价单键；QM区域尽可能保持电中性（可能会有电荷不为整数的情况，这个时候程序会调整边界MM原子或所有MM原子的电荷，使得QM/MM哈密顿量下整个体系仍处于电中性，Amber24手册 p159）。

目前QM区域选择仍然是十分需要人类直觉，配合可视化会较为方便。为此写了一个qm_model.py，作为PyMOL的辅助脚本：

from __future__ import print_function
import pymol
from pymol import cmd
from parmed.amber import AmberParm

def qmmm_help():
    text = """
    useful commend for building a QM/MM model with pymol
    
    hide everything, resname WAT+Na\+
    sele ligand, resname ROH+4GB+0GB
    center ligand
    zoom ligand
    sele pocket, byres (ligand around 5)
    show stick, pocket
    sele center_res, resname ROH+4GB+0GB
    select qm_waters, resn HOH+WAT and (center_res around 3)
    select nearby_res, byres (center_res around 3) and not center_res
    sele center_atom, resname 4GB and name O4
    select nearby_res, byres (center_atom around 4) and not center_res

    select qm_region, none
    select qm_region, qm_region or center_res
    select qm_region, qm_region or (resid 168+308+366 and sidechain)

	qm2amber qm_region,qm_region.dat
	qm2chemshell qm_region,qm_region.dat
    qm_charges qm_region,complex_solv.prmtop
    """

    print(text)

cmd.extend("qmmm_help", qmmm_help)

def qm2amber(selection, qm_region_file=None):
    qm_atom_ids = cmd.identify(selection, mode=0)
    qm_atom_mask = "@"+",".join(str(n) for n in qm_atom_ids)
    print(qm_atom_mask)

    if qm_region_file:
        file = open(qm_region_file,"w")
        file.write("@")
        for item in qm_atom_ids: 
            file.write(f"{item},")
        file.close()


cmd.extend("qm2amber", qm2amber)

def qm2chemshell(selection, qm_region_file=None):
    qm_atom_ids = cmd.identify(selection, mode=0)
    print(qm_atom_ids)
    
    if qm_region_file:
        file = open(qm_region_file,"w")
        file.write("set qm_atoms { ")
        for item in qm_atom_ids: 
            file.write(f"{item[1]} ")
        file.write("}")
        file.close()

cmd.extend("qm2chemshell", qm2chemshell)



def qm_charges(selection,top_file):
    qm_atom_ids = cmd.identify(selection, mode=0)
    parm = AmberParm(top_file)
    charges = [parm.atoms[i-1].charge for i in qm_atom_ids]  # serials 是 1-based
    total_charge = sum(charges)
    print(total_charge)

    return total_charge

cmd.extend("qm_charges", qm_charges)

可以把该脚本放置于pdb的同一个目录，载入pdb后。cd到该路径，load qm_model.py就可以使用其中功能（或者放到一个每次打开PyMOL，程序都会优先扫描的目录）。

qmmm_help输出一系列有助于QM/MM建模的命令
qm2amber和qm2chemshell会输出按程序输入格式定义的qm区域
qm_charges能够计算当前选择的QM区域电荷
有待继续更新

1	`<iframe src="https://github.com/chem1912/blog_picture/tree/main/videos/qm_region.mp4" width="800" height="600"></iframe>`

作业非正常终止的情形下自动保存模拟文件

这不仅仅适用于QM/MM MD模拟，也适用于任何计算耗时较长、重新计算成本高昂的模拟。在集群上进行量化、分子动力学模拟，通常建议将相关文件复制到计算节点的临时目录，在该目录进行相应计算，避免在用户目录中频繁读写文件，从而加快模拟。然而这种方式不便于对模拟过程的监测，且一旦作业被取消或集群断电导致的作业终止，模拟的数据就可能会丢失。通常的作业脚本中，只有程序正常结束，运行到脚本末尾才会把临时目录下的文件复制回用户目录。如果作业临时终止，脚本后面的部分就不会被执行。

模拟文件的保存之重要性，对于高级别的QM/MM增强采样模拟尤为重要。这类模拟以周为时间单位，耗费数月时间只为得到单步反应的自由能曲线。本人从一次炎夏的B3LYP QM/MM Metadynamics得到了该教训——之后的模拟，即使会损失部分效率，也将模拟文件实时输出到用户目录中。如今采用以下的方法就可保证效率的同时及时保存数据：

#!/bin/bash
#SBATCH --job-name=md_relax
#SBATCH --output=MD.o-%j
#SBATCH --error=MD.e-%j
#SBATCH --partition=cpu
#SBATCH --nodelist=b001              
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=8                     


JOB_ID=${SLURM_JOB_ID}
JOB_NAME=${SLURM_JOB_NAME}

ulimit -c 0
ulimit -s unlimited
ulimit -v unlimited

# ---> Define the working directory at computing node <---
export STARTDIR=$(pwd)
export RUNDIR=/tmp/$JOB_ID.tmpdir

mkdir -p $RUNDIR
chmod 700 ${RUNDIR}
cp -r $STARTDIR/* $RUNDIR

# ----> enter working directory and execute your job <----
cd $RUNDIR

source /etc/profile.d/modules.sh
module purge 
module load amber/20
module load orca/5.0.1
export AMBERVER=20
export AMBERHOME=/share/apps/amber/20
export BINDIR=$AMBERHOME/bin
export LIBDIR=$AMBERHOME/lib
export INCDIR=$AMBERHOME/include
export DATDIR=$AMBERHOME/dat
export LD_LIBRARY_PATH=/share/intel/oneapi-2022.0.1/mpi/2021.5.0/lib:/mlatom/software/miniconda3/envs/amber/lib:$LD_LIBRARY_PATH

ulimit -s unlimited
ulimit -c 0


cleanup() {
    echo "the task was canceled, save the results ..."
    rm $RUNDIR/MD*
    cp -r $RUNDIR/* $STARTDIR/ 
    echo "$(date): copy done, exit"
  
    echo "$(date): SIGNAL CAUGHT" 
    echo "$(date): STARTDIR=$STARTDIR" 
    echo "$(date): RUNDIR=$RUNDIR" 
    echo "$(date): RUNDIR contents:" 
    ls -la $RUNDIR 
    
    exit
}

trap 'cleanup' SIGTERM SIGINT EXIT

echo "This job begins at: $(date)" 

sander -O -i cmd_qmmm.in -o cmd_qmmm.log -c rep.c0_with_wat.rst -p complex_solv.prmtop -r cmd_qmmm.rst -x cmd_qmmm.nc  &
SANDER_PID=$!
wait $SANDER_PID

# ----> copy files back <----
rm $RUNDIR/MD*
cp -r $RUNDIR/* $STARTDIR

cd $STARTDIR
rm -rf $RUNDIR

echo $(date)

以上的脚本适用于slurm集群，用scancel取消作业，shell捕捉到scancel发送的信号，就会执行cleanup中的命令。该脚本的关键在于，sander需要在后台运行，并用wait等待进程，以保持shell活跃，有足够的时间执行文件复制。

这种方式同样适用于LSF集群，使得在用bkill终止作业或其他原因导致的作业中断情形下，临时目录的文件能够被复制回用户目录。

Science > Theoretical and Computational Chemistry > QM/MM

#QM/MM #PyMol

适用于QM/MM模拟的一些技巧

http://example.com/2026/06/22/适用于QMMM模拟的一些技巧/

Author

Graminilune

Posted on

June 22, 2026

Licensed under

支持QM/MM-Ewald的Amber24编译 Next