ligprep

conda activate ambertools
# 高斯优化
# 注意检查体系的电荷和自选多重度(eg.将T3的gjf电荷0 1改为1 1)
antechamber -i lig.sdf -fi sdf -o lig.gjf -fo gcrt -at gaff2 -gn "%nproc=4" -gm "%mem=4GB" -gk "#B3LYP/6-31G* em=gd3bj pop=MK iop(6/33=2,6/42=6) opt" -rn MOL  -nc 1  
antechamber -i C12_lig.mol2 -fi mol2 -o lig.gjf -fo gcrt -at gaff2 -gn "%nproc=8" -gm "%mem=8GB" -gk "#B3LYP/6-31G* em=gd3bj opt scrf=solvent=Water" -rn MOL  -nc -4
-nc后面加电荷数，电荷为1就是1，电荷为-1就是-1，在gif文件里有体现
自旋多重度计算：（未配对的正电子数-未配对的负电子数）*2+1=1 一般都为1，遇到金属离子需要特别注意
gif文件如下

%nproc=8
%chk=molecule
%mem=8GB
#B3LYP/6-31G* em=gd3bj opt scrf=solvent=Water


remark line goes here


-4   1
第一个-4指的是电荷数，电荷为-4就是-4
第二个1指的是自旋多重度

antechamber \
    -i PA-coA2.mol2 \
    -fi mol2 \
    -o PA-coA_gas.gjf \
    -fo gcrt \
    -at gaff2 \
    -gn "%nprocshared=8" \
    -gm "%mem=8GB" \
    -gk "#B3LYP/6-31G* em=gd3bj pop=MK iop(6/33=2,6/42=6) SCF=(XQC,VShift=500,MaxCycles=512)" \
    -rn MOL \
    -nc -4


echo "退出码: $?"
head -8 PA-coA_gas.gjf

%nprocshared=8
%chk=molecule_noMg.chk
%mem=8GB
#B3LYP/6-31G* em=gd3bj pop=MK iop(6/33=2,6/42=6) opt


remark line goes here


-4   1

DRAK2蛋白-配体MD（noMg & Mg体系）

环境准备(755)

conda activate ambertools
source ~/.gmx2024.4.sh
export PATH=/home/dddc/rtluo/software/g16:$PATH

一、配体准备（两个体系共用）

1. 高斯计算RESP电荷

%nprocshared=16
%mem=16GB
%chk=molecule.chk
#B3LYP/6-31G* em=gd3bj pop=MK iop(6/33=2,6/42=6) opt scf=(maxcycles=512,xqc,noincfock)
# 提交高斯计算（约21小时）
nohup g16 lig.gjf > lig.log &


# 确认正常结束
tail -3 lig.log
# 应显示：Normal termination of Gaussian 16

2. ligprep目录准备

system/ligprep/
├── lig.log      ← 高斯计算结果（真空RESP）
├── lig.gjf      ← 高斯输入文件
└── lig.mol2     ← 配体坐标文件

二、蛋白准备

noMg体系

conda activate pdb2pqr
cd /home/databank/ydn/DRAK2/MD/noMg/system1


pdb2pqr30 noMg_protein.pdb noMg_protein.pqr \
    --ff AMBER --ffout AMBER \
    --with-ph 7.4 \
    --pdb-output noMg_protein_pH74.pdb

Mg体系

conda activate pdb2pqr
cd /home/databank/ydn/DRAK2/MD/Mg/system2


pdb2pqr30 Mg_protein_re.pdb Mg_protein_re.pqr \
    --ff AMBER --ffout AMBER \
    --with-ph 7.4 \
    --pdb-output Mg_protein_re_pH74.pdb

检查蛋白

# 检查HIS质子化状态（HIE/HID/HIP）
grep "HIS\|HIE\|HID\|HIP" protein_pH74.pdb | grep "^ATOM" | \
    awk '{print $4, $5, $6}' | sort -u


# 检查二硫键（有CYS的SG原子，若两个SG距离<2.5Å则有二硫键）
grep " SG " protein.pdb


# 检查缺失残基（看残基编号是否连续）
grep "^ATOM" protein_pH74.pdb | awk '{print $6}' | sort -n | uniq


# 确认Mg在文件中（Mg体系）
grep "MG" Mg_protein_re_pH74.pdb
tail -5 Mg_protein_re_pH74.pdb

移走多余文件，确保system目录下只有一个pdb

mkdir -p backup
mv noMg_protein.pdb noMg_protein.pqr backup/
ls *.pdb
# 应只剩 noMg_protein_pH74.pdb

三、目录结构

MD/
├── noMg/
│   ├── total_control-v2.py
│   ├── lig_resp_cal.py
│   ├── atom_name_check.py
│   ├── md_parm_gen.py
│   ├── amber_to_gmx-add_restraint.py
│   ├── pre_equ.py
│   └── system1/
│       ├── noMg_protein_pH74.pdb  ← 蛋白文件（唯一pdb）
│       ├── noMg_lig.mol2          ← 配体mol2
│       ├── ligprep/               ← 高斯计算文件
│       ├── parameters/            ← 自动生成
│       ├── mdp/                   ← mdp参数文件
│       └── pre-equ/               ← 自动生成
└── Mg/
    ├── total_control-v2.py
    ├── ...
    └── system2/
        ├── Mg_protein_re_pH74.pdb ← 蛋白文件（唯一pdb）
        ├── Mg_lig.mol2            ← 配体mol2
        └── ...

四、脚本修改（运行前必做）

noMg体系

# 修改GPU编号
sed -i "s/'GPU_DEVICES': \[0, 1, 2\]/'GPU_DEVICES': [5]/" \
    /home/databank/ydn/DRAK2/MD/noMg/pre_equ.py


grep GPU_DEVICES /home/databank/ydn/DRAK2/MD/noMg/pre_equ.py

Mg体系

Atom index (2) in position_restraints out of bounds (1-1)

Fatal error: 1 atoms are not part of any of the T-Coupling groups

# 修复1：添加MG到METAL_IONS
sed -i "s/METAL_IONS = {'MN', 'ZN', 'SO4'}/METAL_IONS = {'MN', 'ZN', 'SO4', 'MG'}/" \
    /home/databank/ydn/DRAK2/MD/Mg/amber_to_gmx-add_restraint.py


# 修复2：添加MG到溶质组
sed -i "s/solute_keywords = \['Protein', 'MOL', 'SO4', 'MN'\]/solute_keywords = ['Protein', 'MOL', 'SO4', 'MN', 'MG']/" \
    /home/databank/ydn/DRAK2/MD/Mg/amber_to_gmx-add_restraint.py


# 修改GPU编号
sed -i "s/'GPU_DEVICES': \[0, 1, 2\]/'GPU_DEVICES': [5]/" \
    /home/databank/ydn/DRAK2/MD/Mg/pre_equ.py


# 确认修改
grep "METAL_IONS" /home/databank/ydn/DRAK2/MD/Mg/amber_to_gmx-add_restraint.py
grep "solute_keywords" /home/databank/ydn/DRAK2/MD/Mg/amber_to_gmx-add_restraint.py
grep "GPU_DEVICES" /home/databank/ydn/DRAK2/MD/Mg/pre_equ.py

五、运行自动化前处理脚本

# noMg
cd /home/databank/ydn/DRAK2/MD/noMg
CUDA_VISIBLE_DEVICES=5 nohup python total_control-v2.py 1 > script.log 2>&1 &
tail -f script.log


# Mg
cd /home/databank/ydn/DRAK2/MD/Mg
CUDA_VISIBLE_DEVICES=5 nohup python total_control-v2.py 1 > script.log 2>&1 &
tail -f script.log

lig_resp_cal.py       ← 读取高斯log，生成配体参数（lig.prep/lig.frcmod）
手动
cd /home/databank/ydn/DRAK2/MD/C18:1/Mg/system/ligprep


cd system/ligprep


antechamber \
  -i lig.log \
  -fi gout \
  -o lig.prep \
  -fo prepi \
  -c resp \
  -nc -4 \
  -rn LIG \
  -s 2


parmchk2 \
  -i lig.prep \
  -f prepi \
  -o lig.frcmod \
  -a y
  
atom_name_check.py    ← 对比原子名和坐标，生成LIG.PDB
md_parm_gen.py        ← tleap组装复合物，生成Amber参数
常见报错
问题已经定位：蛋白 PDB 中存在与 Amber 模板不兼容的氢原子名称：


MET HB1
HIE HD1
CSER HXT


最稳妥的方法是删除蛋白现有氢原子，让 tleap 按 ff14SB 模板重新补氢。保留 HIE 等残基名称，因此组氨酸质子化状态仍由残基名控制。


cd /home/databank/ydn/DRAK2/MD/C14:0/Mg


cp system/14C_protein.pdb 14C_protein_withH.pdb.bak


pdb4amber \
  -i system/14C_protein.pdb \
  -o system/14C_protein_noH.pdb \
  --nohyd


mv system/14C_protein_noH.pdb system/14C_protein.pdb


确认错误氢原子已消失：


grep -E " HB1 | HD1 | HXT " system/14C_protein.pdb


正常应无输出。确认 Mg 没有丢失：


grep -c " MG " 14C_protein_withH.pdb.bak
grep -c " MG " system/14C_protein.pdb


两个数字应相同。然后删除失败参数并重建：


rm -rf system/parameters
python md_parm_gen.py


amber_to_gmx.py/python amber_to_gmx-add_restraint.py 1       ← 转换GROMACS格式，插入位置约束，生成index.ndx
pre_equ.py            ← EM→NVT×3→NPT×3预平衡


改成cpu
cd /home/databank/ydn/DRAK2/MD/C14:0/Mg


sed -i \
's/gmx_mpi mdrun -v -deffnm/gmx_mpi mdrun -v -nb cpu -pme cpu -bonded cpu -update cpu -deffnm/g' \
pre_equ.py


grep -n "mdrun" pre_equ.py




tleap 报错


cd /home/databank/ydn/DRAK2/MD-DRAK2/DCA-16C/Mg/system
cp pro.pdb pro_withH.pdb.bak
pdb4amber \
    -i pro.pdb \
    -o pro_noH.pdb \
    --nohyd

六、手动修复index.ndx

Fatal error:
Group Solute referenced in the .mdp file was not found

gmx_mpi make_ndx -f parameters/gmx.gro -n parameters/index.ndx << EOF
q
EOF

cd /home/databank/ydn/DRAK2/MD/noMg/system1/parameters


gmx_mpi make_ndx -f gmx.gro -n index.ndx -o index.ndx << EOF
name 18 solute
name 19 solvent
q
EOF


grep -i "solute\|solvent" index.ndx

cd /home/databank/ydn/DRAK2/MD/Mg/system2/parameters


gmx_mpi make_ndx -f gmx.gro -n index.ndx -o index.ndx << EOF
21 | 14
name 23 solute
22
name 24 solvent
q
EOF


grep -i "solute\|solvent" index.ndx

七、正式MD（100ns）

生成md.tpr

# noMg
cd /home/databank/ydn/DRAK2/MD/noMg/system1
mkdir -p md


gmx_mpi grompp \
    -f mdp/md_1um.mdp \
    -c pre-equ/npt3/npt3.gro \
    -r pre-equ/npt3/npt3.gro \
    -p parameters/gmx.top \
    -n parameters/index.ndx \
    -o md/md.tpr


# Mg（同上，路径替换为Mg/system2）

运行100ns

# noMg
cd /home/databank/ydn/DRAK2/MD/noMg/system1/md
CUDA_VISIBLE_DEVICES=5 nohup gmx_mpi mdrun -v -deffnm md > md_run.log 2>&1 &
tail -f md_run.log


# Mg

cd /home/databank/ydn/DRAK2/MD/Mg/system2/md
CUDA_VISIBLE_DEVICES=4 nohup gmx_mpi mdrun -v -deffnm md > md_run.log 2>&1 &
tail -f md_run.log

确认完成

tail -5 md_run.log
# 应显示：Performance: xxx ns/day

八、续跑200ns/300ns

# 1. 修改mdp步数
# 200ns: nsteps=100000000
# 300ns: nsteps=150000000
sed -i 's/nsteps          = 50000000/nsteps          = 100000000/' \
    mdp/md_1um.mdp


# 2. 生成新tpr（从md.gro续跑）
gmx_mpi grompp \
    -f mdp/md_1um.mdp \
    -c md/md.gro \
    -r md/md.gro \
    -p parameters/gmx.top \
    -n parameters/index.ndx \
    -o md/md_200ns.tpr


# 3. 续跑
# 注意：-deffnm必须用md，与原来一致，否则checkpoint不匹配
cd md
CUDA_VISIBLE_DEVICES=5 nohup gmx_mpi mdrun \
    -s md_200ns.tpr \
    -deffnm md \
    -cpi md.cpt \
    -append \
    -nb gpu -pme gpu \
    > md_200ns.out 2>&1 &


tail -f md_200ns.out
# 应显示：continuing from step 50000000, 100000.0 ps


续跑：
gmx_mpi grompp -f ../mdp/md_1um.mdp -c ../md_400ns/md.part0004.gro -r ../md_400ns/md.part0004.gro -p ../parameters/gmx.top -n ../parameters/index.ndx -o md_500ns.tpr


CUDA_VISIBLE_DEVICES=3 nohup gmx_mpi mdrun     -s md_500ns.tpr     -deffnm md_500ns     -cpi ../md_400ns/md.cpt     -noappend     -nb gpu -pme gpu     > md_500ns.out 2>&1 &

九、基本分析

cd md/


# 能量（Potential、Temperature、Pressure稳定说明体系正常）
gmx_mpi energy -f md.edr -o energy.xvg << EOF
Potential
Temperature
Pressure
EOF


# RMSD（蛋白骨架相对初始结构的偏差）
gmx_mpi rms -s md.tpr -f md.xtc -o rmsd.xvg << EOF
Backbone
Backbone
EOF


# RMSF（各残基的柔性）
gmx_mpi rmsf -s md.tpr -f md.xtc -o rmsf.xvg << EOF
Backbone
EOF

十、MMGBSA计算

含Na拓扑
parmed ../../../parameters/complex.prmtop << EOF
strip :WAT
outparm native-Na.prmtop
quit
EOF


mkdir -p md/analysis/MMGBSA
cd md/analysis/MMGBSA


# 1. 生成蛋白-配体index
# noMg体系（Protein组1 + MOL组13）
echo -e "1 | 13\nname 18 Protein_MOL\nq" | \
    gmx_mpi make_ndx -f ../../md.gro -o protein_lig.ndx


# Mg体系（Protein组1 + MG组14 + MOL组15）
gmx_mpi make_ndx -f ../../md.gro -o protein_lig.ndx
# 交互输入：
# > 1|14|15 蛋白+mg+配体
# > q


若续跑可以合并轨迹（在没叠加的情况）
gmx_mpi trjcat -f md.xtc md.part0002.xtc  md.part0003.xtc -o md_combined.xtc


# 2. 提取轨迹（去水去离子，修正PBC）
gmx_mpi trjconv \
    -f ../../md.xtc \
    -s ../../md.tpr \
    -n protein_lig.ndx \
    -o protein_lig.xtc或者dry.xtc \
    -pbc mol -ur compact
# noMg选18（Protein_MOL）
# Mg选21（Protein_MG_MOL）


可以用循环：
cd '/home/databank/ydn/DRAK2/MD2/C18:1/Mg/system'


for i in 1 2 3 4
do
    echo "========== 开始处理 md${i} =========="


    if [[ ! -f "md${i}/md.gro" || ! -f "md${i}/md.xtc" || ! -f "md${i}/md.tpr" ]]; then
        echo "md${i} 缺少 md.gro、md.xtc 或 md.tpr，跳过"
        continue
    fi


    mkdir -p "md${i}/analysis/RMSD"
    cd "md${i}/analysis/RMSD" || exit 1


    # 创建 Protein + MG + MOL 索引组
    printf "1|14|15\nq\n" | gmx_mpi make_ndx \
        -f ../../md.gro \
        -o protein_lig.ndx


    # 提取 Protein + MG + MOL 无水轨迹
    printf "21\n" | gmx_mpi trjconv \
        -f ../../md.xtc \
        -s ../../md.tpr \
        -n protein_lig.ndx \
        -o dry.xtc \
        -pbc mol \
        -ur compact


    echo "========== md${i} 处理完成 =========="
    echo


    cd ../../../ || exit 1
done


要对齐
cpptraj
> parm native.prmtop
> trajin dry.xtc
> autoimage
> rms（对齐）
> trajout dry-align.xtc
> run
cd /home/databank/ydn/DRAK2/MD-DRAK2/DCA-16C/Mg/system


for i in 2 3 4
do
    echo "========== md${i}: 生成 dry-align.xtc =========="


    cd "md${i}/analysis/RMSD" || exit 1


    cpptraj > cpptraj_align.log 2>&1 << EOF
parm ../../../parameters/native.prmtop
trajin dry.xtc
autoimage
rms
trajout dry-align.xtc
run
quit
EOF


    echo "md${i} 对齐完成"
    ls -lh dry-align.xtc


    cd ../../../ || exit 1
done
cpptraj << 'EOF'
parm native_na.prmtop
trajin dry-na.xtc
autoimage
rms first :1-291@CA
trajout dry-na-align.xtc
run
quit
EOF


查看组别 
grep "^\[" protein_lig.ndx
# 3. 准备igb5_0-200ns.in
cat igb5_0-200ns.in
# startframe=0, endframe=50000, interval=100（100ns轨迹）


1000帧取1帧
cpptraj
> parm native.prmtop
> trajin dry-align.xtc 1 10000000000 1000
> trajout check.pdb
> run
cd /home/databank/ydn/DRAK2/MD-DRAK2/DCA-16C/Mg/system


for i in 1 2 3 4
do
    echo "========== md${i}: 生成 check.pdb =========="


    cd "md${i}/analysis/RMSD" || exit 1


    cpptraj > cpptraj_check.log 2>&1 << EOF
parm ../../../parameters/native.prmtop
trajin dry-align.xtc 1 last 1000
trajout check.pdb
run
quit
EOF


    echo "md${i} 检查轨迹生成完成"
    ls -lh check.pdb


    cd ../../../ || exit 1
done




去Mg：
parmed native.prmtop << EOF
strip :MG
outparm native_noMg.prmtop
quit
EOF
parmed pro.prmtop << EOF
strip :MG
outparm pro_noMg.prmtop
quit
EOF
# 4. 运行MMGBSA
nohup mpirun -np 20 MMPBSA.py.MPI \
    -i igb5_0-200ns.in \
    -o final.dat \
    -cp ../../../parameters/native.prmtop \
    -rp ../../../parameters/pro.prmtop \
    -lp ../../../parameters/lig.prmtop \
    -y protein_lig.xtc \
    -eo mmgbsa.csv > mmpbsa.log 2>&1 &
    
    base="/home/databank/ydn/DRAK2/MD-DRAK2/C12:0/Mg/system"
for i in 3 4; do
    workdir="${base}/md${i}/analysis/MMGBSA"


    mkdir -p "$workdir"
    cp "${base}/md1/analysis/MMGBSA/igb5_0-500ns.in" "$workdir/"


    (
        cd "$workdir" || exit 1


        nohup mpirun -np 20 MMPBSA.py.MPI \
            -i igb5_0-500ns.in \
            -o final.dat \
            -cp ../../../parameters/native.prmtop \
            -rp ../../../parameters/pro.prmtop \
            -lp ../../../parameters/lig.prmtop \
            -y ../RMSD/dry-align.xtc \
            -eo mmgbsa.csv \
            > mmpbsa.log 2>&1 &


        echo "md${i} 已启动，PID=$!"
    )


    sleep 2done
    nohup mpirun -np 20 MMPBSA.py.MPI     -i igb5_300-500ns.in     -o final.dat     -cp ../../../../parameters/native.prmtop     -rp ../../../../parameters/pro.prmtop     -lp ../../../../parameters/lig.prmtop     -y ../../RMSD/dry-align.xtc     -eo mmgbsa.csv > mmpbsa.log 2>&1 & 
nohup mpirun -np 20 MMPBSA.py.MPI \
    -i ../igb5_300-500ns.in \
    -o final.dat \
    -cp ../../../../../parameters/native_noMg.prmtop \
    -rp ../../../../../parameters/pro_noMg.prmtop\
    -lp ../../../../../parameters/lig.prmtop \
    -y ../../../RMSD/dry-align_noMg.xtc  \
    -eo mmgbsa-nomg.csv > mmpbsa-nomg.log 2>&1 &


tail -f mmpbsa.log


mpirun -np 20 MMPBSA.py.MPI     -i igb5_300-500ns.in     -o final.dat     -cp ../../../../parameters/native.prmtop     -rp ../../../../parameters/pro.prmtop    -lp ../../../../parameters/lig.prmtop     -y ../../RMSD/dry.xtc      -eo mmgbsa.csv > mmpbsa.log 2>&1 &

# 1. 计算所有Na+到配体的距离（300-500ns，每1000帧取一帧）
cpptraj << 'EOF'
parm native_na.prmtop
trajin dry-na-align.xtc 150001 250001 1000
distance na293 :293 :292 out dist_na.dat
distance na294 :294 :292 out dist_na.dat
distance na295 :295 :292 out dist_na.dat
distance na296 :296 :292 out dist_na.dat
distance na297 :297 :292 out dist_na.dat
distance na298 :298 :292 out dist_na.dat
distance na299 :299 :292 out dist_na.dat
distance na300 :300 :292 out dist_na.dat
distance na301 :301 :292 out dist_na.dat
distance na302 :302 :292 out dist_na.dat
distance na303 :303 :292 out dist_na.dat
distance na304 :304 :292 out dist_na.dat
distance na305 :305 :292 out dist_na.dat
distance na306 :306 :292 out dist_na.dat
distance na307 :307 :292 out dist_na.dat
distance na308 :308 :292 out dist_na.dat
distance na309 :309 :292 out dist_na.dat
distance na310 :310 :292 out dist_na.dat
distance na311 :311 :292 out dist_na.dat
run
quit
EOF


# 2. 计算平均距离并排序（最近的在最上面）
awk 'NR==1{next} {for(i=2;i<=NF;i++) sum[i]+=$i; count++} \
    END{for(i=2;i<=NF;i++) printf "Na%d\t%.2f Å\n", i+291, sum[i]/count}' \
    dist_na.dat | sort -k2 -n
    
建pro_na302.prmtopcpptraj << 'EOF'parm ../../../parameters/complex.prmtopparmstrip :WATparmstrip :MOLparmstrip :293-301,303-311parmwrite out pro_na302.prmtopquitEOF
show sticks, resn MOL
select ion_test, inorganic
show spheres, ion_test
set sphere_scale, 0.5, ion_test
显示某原子：pymol
select na302, resi 302
show spheres, na302
color orange, na302


nohup mpirun -np 20 MMPBSA.py.MPI \    -i igb5_0-500ns-250.in \    -o final_na302.dat \    -cp ../RMSD-NA/native_na302.prmtop \    -rp ../RMSD-NA/pro_na302.prmtop \    -lp ../../../parameters/lig.prmtop \    -y ../RMSD-NA/dry-na302-align.xtc \    -eo mmgbsa_na302.csv > mmpbsa_na302.log 2>&1 &echo "PID: $!"

cd /home/databank/ydn/DRAK2/MD2/Mg/system2/md1
gmx_mpi mdrun -v -deffnm md -cpi md.cpt


cd /home/databank/ydn/DRAK2/MD2/Mg/system2/md2
gmx_mpi mdrun -v -deffnm md -cpi md.cpt


cd /home/databank/ydn/DRAK2/MD2/Mg/system2/md3
gmx_mpi mdrun -v -deffnm md -cpi md.cpt


nohup bash -c 'export CUDA_VISIBLE_DEVICES=5; gmx_mpi mdrun -v -deffnm md -cpi md.cpt -gpu_id 0' > md_run.log 2>&1 &

常见报错及修复

FileNotFoundError: tleap

conda activate ambertools

FATAL: atom does not have a type

pdb4amber -i protein.pdb -o protein_fix.pdb --nohyd

posre out of bounds (1-1)

1 atoms not in T-Coupling groups

Group Solute not found

Inconsistency in checkpoint

-deffnm与checkpoint记录的文件名不一致

OXT no type

grep -v "OXT" protein.pdb > tmp.pdb