快速将RESP电荷从.chg文件迁移至.mol2文件

chemzhang · 发表于 Post on 2022-3-3 15:30:22

本帖最后由 chemzhang 于 2022-3-4 10:15 编辑

在使用Amber做分子动力学时常常需要含有RESP电荷的mol2文件。而multiwfn计算的RESP电荷储存在chg文件之中。尽管notepad++可以方便地将chg文件中的电荷数据迁移至原有的mol2文件，但在Linux系统中我希望能够直接通过代码完成迁移，从而舍弃掉鼠标操作，实现对Multiwfn更好的利用。因此我写了一段shell程序分享给大家。

使用脚本的时候只需在linux端直接运行该程序即可，然后根据提示逐步完成操作。比较方便的做法是把该脚本、mol2文件与chg文件放在同一目录下，免去输入路径容易出错的情形。
以amber官网中计算离子液体为例,Simulations of a room-temperature ionic liquid (ambermd.org),对BF4-离子的mol2文件，教程里记录了复杂的修改过程。而我直接用g16优化BF4-离子后得到高斯输出文件，然后利用open babel将gout文件转为mol2格式，再利用multiwfn根据chk文件获取chg文件，得到下图的一系列文件
202203031508238053..png

最初得到的BF4-.mol2文件如下所示。文件中的电荷是来自高斯优化过程的MULLIKEN电荷

@<TRIPOS>MOLECULE
BF4-.out
5 4 0 0 0
SMALL
MULLIKEN_CHARGES
****
Gaussian 16 # opt hf/6-31g geom=connectivity
@<TRIPOS>ATOM
1 B -0.0001 0.0004 0.0003 B 1 UNL1 1.2546
2 F -0.8199 -1.0524 -0.4986 F 1 UNL1 -0.5639
3 F 0.3731 -0.2852 1.3441 F 1 UNL1 -0.5634
4 F 1.1707 0.1117 -0.8035 F 1 UNL1 -0.5639
5 F -0.7238 1.2257 -0.0421 F 1 UNL1 -0.5634
@<TRIPOS>BOND
1 4 1 1
2 2 1 1
3 5 1 1
4 1 3 1

复制代码

然后我希望将mulliken电荷替换为chg文件中储存的resp电荷，运行chg2mol2.sh即可

[chemzhang ~]$ ./chg2mol2.sh
please enter your .mol2 file's path (eg. ~/xxx.mol2): BF4-.mol2
please enter your .chg file's path (eg. ~/xxx.chg): BF4-.chg
MULLIKEN charges from BF4-.mol2:
1.2546 -0.5639 -0.5634 -0.5639 -0.5634
RESP charges from BF4-.chg:
1.3753509763 -0.5940223436 -0.5936477442 -0.5939881453 -0.5936927433
please enter the scaled factor (eg. 0.8): 0.8
scaled RESP charges:
1.1002807810 -0.4752178749 -0.4749181954 -0.4751905162 -0.4749541946
current molecule's name in BF4-_RESP.mol2: UNL1
DO you wish to update it? [y/n]: y
please enter the new name: BF4
done

复制代码

运行脚本时按照提示逐步操作。首先输入mol2文件与chg文件地址（如上所述，将该脚本、mol2文件与chg文件放在同一目录下在这一步会方便不少）。然后脚本会自动读取原mol2文件中的MULLIKEN电荷与chg文件中的RESP电荷并显示出来，然后将chg文件中的RESP电荷写入_RESP.mol2文件。
下一步是询问是否要约化电荷。其可以直接根据你输入的scaled factor计算出约化电荷并写入_RESP.mol2文件，免去了手动或使用excel计算等麻烦。如果不想约化，将相关的代码注释掉或者在使用时输入1。这里我设置约化系数为0.8。
由于AMBER中使用packmol建立模拟盒子需要对residue有命名（且不超过三个字母/数字）。因此脚本还会获取当前的residue的名称并显示出来，询问你是否需要改名。我这里将分子由obabel自动生成的UNL1重新命名为BF4。
最终得到的BF4-_RESP.mol2文件，里面的电荷已经根据chg中的内容替换为经过约化的RESP电荷，名称也已成为BF4,可直接用于后续amber的frcmod文件生成，比教程里人为指定成键情况修改mol2文件要方便太多

@<TRIPOS>MOLECULE
BF4-.out
5 4 0 0 0
SMALL
SCALED RESP_CHARGES WITH SCALED FACTOR 0.8
****
Gaussian 16 # opt hf/6-31g geom=connectivity
@<TRIPOS>ATOM
1 B -0.0001 0.0004 0.0003 B 1 BF4 1.1002807810
2 F -0.8199 -1.0524 -0.4986 F 1 BF4 -0.4752178749
3 F 0.3731 -0.2852 1.3441 F 1 BF4 -0.4749181954
4 F 1.1707 0.1117 -0.8035 F 1 BF4 -0.4752178749
5 F -0.7238 1.2257 -0.0421 F 1 BF4 -0.4749181954
@<TRIPOS>BOND
1 4 1 1
2 2 1 1
3 5 1 1
4 1 3 1

复制代码

脚本虽然比较简陋，但我本身没有系统接受过有关计算机的教育，bash语言也是自学而来，我的目标是能用就行。脚本中有几个主要的功能，我都加了相应的注释说明。只要有linux基础都可以看懂，并根据自己的需求做出修改。

#! /usr/bin/env bash
# transplant the RESP charges from .chg file to .mol2 file
# chemzhang 2022-2-28
# assign input files
read -r -p "please enter your .mol2 file's path (eg. ~/xxx.mol2): " mol2file
read -r -p "please enter your .chg file's path (eg. ~/xxx.chg): " chgfile
# locate to the row where mol2 file stores charge value
mol2NRATOM=$(grep -n ATOM "$mol2file" | awk 'BEGIN {FS=":"} {print $1}') # locate to the row containing "ATOM" flag in mol2 file
mol2NRBOND=$(grep -n BOND "$mol2file" | awk 'BEGIN {FS=":"} {print $1}') ## locate to the row containing "BOND" flag in mole2 file
mol2NRfirst=$(("$mol2NRATOM"+1)) # equivalent to mol2NRfirst=$(echo "$mol2NRATOM+1" | bc)
mol2NRlast=$(echo "$mol2NRBOND-1" | bc) # equivalent to mol2NRlast=$(("$mol2NRATOM"-1))
num=$(("${mol2NRlast}"-"${mol2NRfirst}"+1)) # the number of to-be-updated rows
# store MULLIKEN charge from mol2 file as a shell array
k=0
MULLIKEN=()
for i in $(seq $(("$mol2NRfirst")) 1 "$mol2NRlast");
do MULLIKEN[$k]=$(head -n "${i}" "${mol2file}" | tail -n 1 | awk '{print $NF}') ;
# echo "$k ${MULLIKEN[$k]}";
((k++));
done
echo "MULLIKEN charge from ${mol2file}:"
echo "${MULLIKEN[@]}"
echo -e "\n"
# store MULLIKEN charges from mol2 file as a shell array
k=0
RESP=()
for i in $(seq 1 1 ${num});
do RESP[$k]=$(head -n "${i}" "${chgfile}" | tail -n 1 | awk '{print $NF}') ;
# echo "$k ${RESP[$k]}";
((k++));
done
echo "RESP charge from ${chgfile}:"
echo "${RESP[@]}"
echo -e "\n"
# scale the RESP charges
read -r -p "please enter the scaled factor (eg. 0.8): " sf
scaled_RESP=("${RESP[@]}")
for i in $(seq 0 1 $((${#scaled_RESP[@]}-1)));
do scaled_RESP[$i]=`awk -v x=${scaled_RESP[$i]} -v y="$sf" 'BEGIN{printf "%.10f", x*y}'`;
# echo "$k ${scaled_RESP[$k]}";
done
echo "scaled RESP charges:"
echo "${scaled_RESP[@]}"
echo -e "\n"
# replace the MULLIKEN charges in mol2 file with scaled_RESP charges
cp "${mol2file}" "${mol2file}".copy
for i in $(seq 0 1 $(("${num}"-1)));
do sed -i "s/${MULLIKEN[$i]}/${scaled_RESP[$i]}/g" "${mol2file}";
done
# post-processing. edit this part depending on your own commands.
sed -i "s/MULLIKEN_CHARGES/SCALED RESP_CHARGES WITH SCALED FACTOR ${sf}/g" "${mol2file}"
mv "${mol2file}" "${mol2file%.*}"_RESP.mol2
mv "${mol2file}".copy "${mol2file}"
#re-name the molecule in _RESP.mol2 file
name=$(head -n "${mol2NRfirst}" "${mol2file%.*}_RESP.mol2" | tail -n 1 | awk '{print $((NF-1))}')
echo "current molecule's name in ${mol2file%.*}_RESP.mol2: ${name}"
read -r -p "DO you wish to update it? [y/n]: " flag
if [ "$flag" != 'y' ];
then
echo "original name "${name}" will be kept"
else
read -r -p "please enter the new name: " newname
sed -i "s/ ${name} / ${newname} /g" "${mol2file%.*}"_RESP.mol2
fi
#generate frcmod file for amber
#parmchk2 -i "${mol2file%.*}"_RESP.mol2 -f mol2 -o "${mol2file%.*}".frcmod
echo "done"

复制代码

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[Multiwfn资源与经验] 快速将RESP电荷从.chg文件迁移至.mol2文件

评分 Rate