计算化学公社

 找回密码 Forget password
 注册 Register
Views: 2395|回复 Reply: 0

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

[复制链接 Copy URL]

6

帖子

1

威望

101

eV
积分
127

Level 2 能力者

发表于 Post on 2022-3-3 15:30:22 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 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电荷
  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


  6. RESP charges from BF4-.chg:
  7. 1.3753509763 -0.5940223436 -0.5936477442 -0.5939881453 -0.5936927433

  8. please enter the scaled factor (eg. 0.8): 0.8
  9. scaled RESP charges:
  10. 1.1002807810 -0.4752178749 -0.4749181954 -0.4751905162 -0.4749541946

  11. current molecule's name in BF4-_RESP.mol2: UNL1
  12. DO you wish to update it? [y/n]: y
  13. please enter the new name: BF4
  14. 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

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

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

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

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

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

  45. # replace the MULLIKEN charges in mol2 file with scaled_RESP charges
  46. cp "${mol2file}" "${mol2file}".copy
  47. for i in $(seq 0 1 $(("${num}"-1)));
  48. do sed -i "s/${MULLIKEN[$i]}/${scaled_RESP[$i]}/g" "${mol2file}";
  49. done

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



chg2mol2.sh

2.81 KB, 下载次数 Times of downloads: 25

评分 Rate

参与人数
Participants 4
威望 +1 eV +6 收起 理由
Reason
亚切ss + 2 好物!
hdhxx123 + 3 好物!
丁越 + 1 赞!
sobereva + 1

查看全部评分 View all ratings

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 03:51 , Processed in 0.655196 second(s), 26 queries .

快速回复 返回顶部 返回列表 Return to list