计算化学公社

标题: 快速将RESP电荷从.chg文件迁移至.mol2文件 [打印本页]
作者
Author: chemzhang    时间: 2022-3-3 15:30
标题: 快速将RESP电荷从.chg文件迁移至.mol2文件
本帖最后由 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文件，得到下图的一系列文件
(, 下载次数 Times of downloads: 49)

最初得到的BF4-.mol2文件如下所示。文件中的电荷是来自高斯优化过程的MULLIKEN电荷
  1. @<TRIPOS>MOLECULE
  2. BF4-.out
  3. 5 4 0 0 0
  4. SMALL
  5. MULLIKEN_CHARGES
  6. ****
  7. Gaussian 16 # opt hf/6-31g geom=connectivity
  8. @<TRIPOS>ATOM
  9.       1 B          -0.0001    0.0004    0.0003 B       1  UNL1        1.2546
  10.       2 F          -0.8199   -1.0524   -0.4986 F       1  UNL1       -0.5639
  11.       3 F           0.3731   -0.2852    1.3441 F       1  UNL1       -0.5634
  12.       4 F           1.1707    0.1117   -0.8035 F       1  UNL1       -0.5639
  13.       5 F          -0.7238    1.2257   -0.0421 F       1  UNL1       -0.5634
  14. @<TRIPOS>BOND
  15.      1     4     1    1
  16.      2     2     1    1
  17.      3     5     1    1
  18.      4     1     3    1
复制代码


然后我希望将mulliken电荷替换为chg文件中储存的resp电荷，运行chg2mol2.sh即可
  1. [chemzhang ~]$ ./chg2mol2.sh
  2. please enter your .mol2 file's path (eg. ~/xxx.mol2): BF4-.mol2
  3. please enter your .chg file's path (eg. ~/xxx.chg): BF4-.chg
  4. MULLIKEN charges from BF4-.mol2:
  5. 1.2546 -0.5639 -0.5634 -0.5639 -0.5634


  8. RESP charges from BF4-.chg:
  9. 1.3753509763 -0.5940223436 -0.5936477442 -0.5939881453 -0.5936927433

  11. please enter the scaled factor (eg. 0.8): 0.8
  12. scaled RESP charges:
  13. 1.1002807810 -0.4752178749 -0.4749181954 -0.4751905162 -0.4749541946

  15. current molecule's name in BF4-_RESP.mol2: UNL1
  16. DO you wish to update it? [y/n]: y
  17. please enter the new name: BF4
  18. 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文件要方便太多
  1. @<TRIPOS>MOLECULE
  2. BF4-.out
  3. 5 4 0 0 0
  4. SMALL
  5. SCALED RESP_CHARGES WITH SCALED FACTOR 0.8
  6. ****
  7. Gaussian 16 # opt hf/6-31g geom=connectivity
  8. @<TRIPOS>ATOM
  9.       1 B          -0.0001    0.0004    0.0003 B       1  BF4        1.1002807810
  10.       2 F          -0.8199   -1.0524   -0.4986 F       1  BF4       -0.4752178749
  11.       3 F           0.3731   -0.2852    1.3441 F       1  BF4       -0.4749181954
  12.       4 F           1.1707    0.1117   -0.8035 F       1  BF4       -0.4752178749
  13.       5 F          -0.7238    1.2257   -0.0421 F       1  BF4       -0.4749181954
  14. @<TRIPOS>BOND
  15.      1     4     1    1
  16.      2     2     1    1
  17.      3     5     1    1
  18.      4     1     3    1
复制代码



脚本虽然比较简陋，但我本身没有系统接受过有关计算机的教育，bash语言也是自学而来，我的目标是能用就行。脚本中有几个主要的功能，我都加了相应的注释说明。只要有linux基础都可以看懂，并根据自己的需求做出修改。
  1. #! /usr/bin/env bash
  2. # transplant the RESP charges from .chg file to .mol2 file
  3. # chemzhang 2022-2-28

  5. # assign input files
  6. read -r -p "please enter your .mol2 file's path (eg. ~/xxx.mol2): " mol2file
  7. read -r -p "please enter your .chg file's path (eg. ~/xxx.chg): " chgfile

  9. # locate to the row where mol2 file stores charge value
  10. mol2NRATOM=$(grep -n ATOM "$mol2file" | awk 'BEGIN {FS=":"} {print $1}')    # locate to the row containing "ATOM" flag in mol2 file
  11. mol2NRBOND=$(grep -n BOND "$mol2file" | awk 'BEGIN {FS=":"} {print $1}')    ## locate to the row containing "BOND" flag in mole2 file
  12. mol2NRfirst=$(("$mol2NRATOM"+1))  # equivalent to mol2NRfirst=$(echo "$mol2NRATOM+1" | bc)
  13. mol2NRlast=$(echo "$mol2NRBOND-1" | bc) # equivalent to mol2NRlast=$(("$mol2NRATOM"-1))
  14. num=$(("${mol2NRlast}"-"${mol2NRfirst}"+1)) # the number of to-be-updated rows

  16. # store MULLIKEN charge from mol2 file as a shell array
  17. k=0
  18. MULLIKEN=()
  19. for i in $(seq $(("$mol2NRfirst")) 1 "$mol2NRlast");
  20. do MULLIKEN[$k]=$(head -n "${i}" "${mol2file}" | tail -n 1 | awk '{print $NF}') ;
  21. #   echo "$k   ${MULLIKEN[$k]}";
  22.    ((k++));
  23. done
  24. echo "MULLIKEN charge from ${mol2file}:"
  25. echo "${MULLIKEN[@]}"
  26. echo -e "\n"

  28. # store MULLIKEN charges from mol2 file as a shell array
  29. k=0
  30. RESP=()
  31. for i in $(seq 1 1 ${num});
  32. do RESP[$k]=$(head -n "${i}" "${chgfile}" | tail -n 1 | awk '{print $NF}') ;
  33. #   echo "$k    ${RESP[$k]}";
  34.    ((k++));
  35. done
  36. echo "RESP charge from ${chgfile}:"
  37. echo "${RESP[@]}"
  38. echo -e "\n"

  40. # scale the RESP charges
  41. read -r -p "please enter the scaled factor (eg. 0.8): " sf
  42. scaled_RESP=("${RESP[@]}")
  43. for i in $(seq 0 1 $((${#scaled_RESP[@]}-1)));
  44. do scaled_RESP[$i]=`awk -v x=${scaled_RESP[$i]} -v y="$sf" 'BEGIN{printf "%.10f", x*y}'`;
  45. #   echo "$k    ${scaled_RESP[$k]}";
  46. done
  47. echo "scaled RESP charges:"
  48. echo "${scaled_RESP[@]}"
  49. echo -e "\n"

  51. # replace the MULLIKEN charges in mol2 file with scaled_RESP charges
  52. cp "${mol2file}" "${mol2file}".copy
  53. for i in $(seq 0 1 $(("${num}"-1)));
  54. do sed -i "s/${MULLIKEN[$i]}/${scaled_RESP[$i]}/g" "${mol2file}";
  55. done

  57. # post-processing. edit this part depending on your own commands.
  58. sed -i "s/MULLIKEN_CHARGES/SCALED RESP_CHARGES WITH SCALED FACTOR ${sf}/g" "${mol2file}"
  59. mv "${mol2file}" "${mol2file%.*}"_RESP.mol2
  60. mv "${mol2file}".copy "${mol2file}"
  61.    #re-name the molecule in _RESP.mol2 file
  62. name=$(head -n "${mol2NRfirst}" "${mol2file%.*}_RESP.mol2" | tail -n 1 | awk '{print $((NF-1))}')
  63. echo "current molecule's name in ${mol2file%.*}_RESP.mol2: ${name}"
  64. read -r -p "DO you wish to update it? [y/n]: " flag
  65. if [ "$flag" != 'y' ];
  66. then
  67.     echo "original name "${name}" will be kept"
  68. else
  69.     read -r -p "please enter the new name: " newname
  70.     sed -i "s/ ${name} / ${newname} /g" "${mol2file%.*}"_RESP.mol2
  71. fi
  72.    #generate frcmod file for amber
  73. #parmchk2 -i "${mol2file%.*}"_RESP.mol2 -f mol2 -o "${mol2file%.*}".frcmod
  74. echo "done"
复制代码



作者
Author: zeus12138    时间: 2025-7-12 09:13
本帖最后由 zeus12138 于 2025-7-12 09:56 编辑

chemzhang先生您好，请问我使用了您的脚本，如何进行引用或者如何标注出您呢？chemzhang先生，我在使用您的脚本的时候发现只有第一个原子电荷写入到了mol2，但是后面的原子和第一个原子的原子电荷都是一样的，想请教下您我要怎么解决嘞
作者
Author: zeus12138    时间: 2025-7-13 08:40
zeus12138 发表于 2025-7-12 09:13
chemzhang先生您好，请问我使用了您的脚本，如何进行引用或者如何标注出您呢？chemzhang先生，我在使用您的 ...

感谢老师，我解决了第二个的电荷重复的问题，是因为没有重新仔细看sob老师在sobtop上对于例2的讲解，以及没有使用ORCA对于小分子优化。



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3