力场拟合工具xff开发杂记

类别:    标签: gmx js   阅读次数:   版权: (CC) BY-NC-SA

引介

xff是一个用于拟合力场参数的在线工具, 可以生成GROMACS格式的拓扑文件. 目前它的功能尚简陋, 只能基于力常数矩阵(力常数为体系能量对坐标的二阶导数, 力常数矩阵常称为Hessian, 中文也有称海森矩阵, 或海森, 海塞的, 本文中称为海森, 以避免与力场参数中的力常数混淆)生成力场的简谐势参数. 后续可能会添加更多的功能, 也可能不会, 凭天意.

这篇文本是我开发xff工具过程中遇到或想到东西的杂乱记录, 算不得正式说明. 我暂时也没有时间将它爬梳成连贯有序, 逻辑自洽的正式文档, 那就采取这种所谓的论语体吧, 想到哪是哪, 想到哪说哪. 既然是杂记, 难免存在思维跳跃之处, 希望那只是因为我头脑活跃, 而不是因为痴呆健忘.

我不想把它梳理成比较正式的文档还有一个私心, 是我希望开发过程透明, 尽量展示而不是隐藏其中的一些细节, 以便有需要的人可以更清楚地了解实际开发过程是如何进行的, 并能从中学到些东西. 在正式发表的论文中, 无论公式推导还是实验步骤, 多是大混乱后清理出的道路, 看起来逻辑合理, 走起来干净通畅, 但实际过程往往远非如此. 数学王子高斯, 后人评价他像只狐狸, 会清理掉自己走过的痕迹, 单留下漂亮的结果让人赞叹而无法学习. 要学习, 只得去研究他未烧掉的草稿, 甚是可惜. 我当然没有狂妄到无知的地步, 只是在这里留下草稿, 凭需要者自取.

杂记

代码

gmxtools.

界面说明

测试

c6h6.fchk来自JCTC作者所给示例, 频率缩放因子为0.957.

苯键长项, 键长r单位Å, 力常数K/2单位kcal/mol/A²  
原子i-j 类型 r(py) r(js) K/2(py) K/2(js-QR) K/2(js-Jacobi)
1 2 C-C 1.390 1.390394 365.36 365.363991 352.426152
1 3 C-C 1.390 1.390238 366.18 366.180366 353.483247
1 4 C-H 1.084 1.084071 350.91 350.909252 350.909252
2 5 C-C 1.390 1.390238 366.18 366.180366 353.483247
2 6 C-H 1.084 1.084071 350.91 350.909252 350.909252
3 9 C-C 1.390 1.390234 366.19 366.188070 353.491432
3 12 C-H 1.084 1.084136 350.98 350.976623 350.976619
5 7 C-C 1.390 1.390234 366.19 366.188070 353.491432
5 8 C-H 1.084 1.084136 350.98 350.976623 350.976619
7 9 C-C 1.390 1.390398 365.35 365.354387 352.415876
7 10 C-H 1.084 1.084071 350.91 350.909266 350.909265
9 11 C-H 1.084 1.084071 350.91 350.909266 350.909265
平均 C-C 1.390 1.390289 365.91 365.909209 353.131898
C-H 1.084 1.084093 350.93 350.931714 350.931712
js-QR对角化程序可以得到与py脚本一致的结果
对称化方法得到的值稍有差距
苯键角项, 键角θ单位°, 力常数K/2单位kcal/mol/rad²  
原子i-j-k 类型 θ(py) θ(js) K/2(py) K/2(js-QR) K/2(js-Jacobi)
2 1 3 C-C-C 120.0 119.998987 96.96 96.960206 66.029786
2 1 4 C-C-H 120.0 119.971551 30.87 30.869577 28.686460
3 1 4 C-C-H 120.0 120.029462 30.87 30.867021 28.770805
1 2 5 C-C-C 120.0 119.998987 96.96 96.960206 66.029786
1 2 6 C-C-H 120.0 119.971551 30.87 30.869577 28.686460
5 2 6 C-C-H 120.0 120.029462 30.87 30.867021 28.770805
1 3 9 C-C-C 120.0 120.002046 97.94 97.940921 66.478172
1 3 12 C-C-H 120.0 119.998832 30.16 30.156745 28.108635
9 3 12 C-C-H 120.0 119.999123 30.16 30.156688 28.108598
2 5 7 C-C-C 120.0 120.002046 97.94 97.940921 66.478172
2 5 8 C-C-H 120.0 119.998832 30.16 30.156745 28.108635
7 5 8 C-C-H 120.0 119.999123 30.16 30.156688 28.108598
5 7 9 C-C-C 120.0 119.998968 96.96 96.960406 66.029859
5 7 10 C-C-H 120.0 120.030060 30.87 30.867122 28.770812
9 7 10 C-C-H 120.0 119.970972 30.87 30.869815 28.686559
3 9 7 C-C-C 120.0 119.998968 96.96 96.960406 66.029859
3 9 11 C-C-H 120.0 120.030060 30.87 30.867122 28.770812
7 9 11 C-C-H 120.0 119.970972 30.87 30.869815 28.686559
平均 C-C-C 120.0 120.000000 97.29 97.287178 66.179272
C-C-H 120.0 120.000000 30.63 30.631161 28.521978
对称化对键角项影响更大.
苯二面角项, 二面角φ单位°, 力常数K/2单位kcal/mol/rad²  
原子i-j-k-l 类型 φ(js) K/2(js) K/2(js-Jacobi)
3 1 2 5 C-C-C-C 0.000000 48.833360 48.833373
3 1 2 6 C-C-C-H 180.000000 25.439312 25.439316
4 1 2 5 C-C-C-H 180.000000 25.439312 25.439316
4 1 2 6 H-C-C-H 0.000000 17.199663 17.199665
2 1 3 9 C-C-C-C 0.000000 48.798538 48.798536
2 1 3 12 C-C-C-H 180.000000 25.169302 25.169301
4 1 3 9 C-C-C-H 180.000000 25.416863 25.416866
4 1 3 12 H-C-C-H 0.000000 17.069950 17.069951
1 2 5 7 C-C-C-C 0.000000 48.798538 48.798536
1 2 5 8 C-C-C-H 180.000000 25.169302 25.169301
6 2 5 7 C-C-C-H 180.000000 25.416863 25.416866
6 2 5 8 H-C-C-H 0.000000 17.069950 17.069951
1 3 9 7 C-C-C-C 0.000000 48.799333 48.799324
1 3 9 11 C-C-C-H 180.000000 25.416703 25.416697
12 3 9 7 C-C-C-H 180.000000 25.169312 25.169312
12 3 9 11 H-C-C-H 0.000000 17.069785 17.069783
2 5 7 9 C-C-C-C 0.000000 48.799333 48.799324
2 5 7 10 C-C-C-H 180.000000 25.416703 25.416697
8 5 7 9 C-C-C-H 180.000000 25.169312 25.169312
8 5 7 10 H-C-C-H 0.000000 17.069785 17.069783
5 7 9 3 C-C-C-C 0.000000 48.832735 48.832722
5 7 9 11 C-C-C-H 180.000000 25.439419 25.439413
10 7 9 3 C-C-C-H 180.000000 25.439419 25.439413
10 7 9 11 H-C-C-H 0.000000 17.199838 17.199835
平均 C-C-C-C 0.000000 48.810306 48.810302
C-C-C-H 180.000000 25.341818 25.341817
H-C-C-H 0.000000 17.113162 17.113161
py脚本暂时无法计算二面角项, 只给出xff的结果
由于分子对称性很高, 对称化对二面角结果影响很小
苯反常二面角项, 二面角ξ单位°, 力常数K/2单位kcal/mol/rad²  
原子i-j-k-l 类型 ξ(js) K/2(js) K/2(js-Jacobi)
1 2 3 4 C-C-C-H_imp 0.000000 120.715427 120.715415
2 1 5 6 C-C-C-H_imp 0.000000 120.715427 120.715415
3 1 9 12 C-C-C-H_imp 0.000000 120.412621 120.412619
5 2 7 8 C-C-C-H_imp 0.000000 120.412621 120.412619
7 5 9 10 C-C-C-H_imp 0.000000 120.715598 120.715586
9 3 7 11 C-C-C-H_imp 0.000000 120.715598 120.715586
平均 C-C-C-H_imp 0.000000 120.614549 120.614540
py脚本暂时无法计算反常二面角项, 只给出xff的结果
由于分子对称性很高, 对称化对二面角结果影响很小

SF6

sf6.fchk来自sobtop所带示例, 频率缩放因子为1.

SF6键长项, 键长r单位Å, 力常数K/2单位kcal/mol/A²  
原子i-j 类型 r(py) r(js) K/2(py) K/2(js-QR) K/2(js-Jacobi)
1 2 S-F 1.578 1.578071 233.64 232.849562 232.849562
1 3 S-F 1.578 1.578071 233.64 232.849562 232.849562
1 4 S-F 1.578 1.578071 233.64 232.849562 232.849562
1 5 S-F 1.578 1.578071 233.64 232.849563 232.849563
1 6 S-F 1.578 1.578071 233.64 232.849562 232.849562
1 7 S-F 1.578 1.578071 233.64 232.849563 232.849563
平均 S-F 1.578 1.578071 233.64 232.849562 232.849562
此分子对称性高, 对称化影响很小, 所得力常数与py稍有差距
SF6键角项, 键角θ单位°, 力常数K/2单位kcal/mol/rad²  
原子i-j-k 类型 θ(py) θ(js) K/2(py) K/2(js-QR) K/2(js-Jacobi)
2 1 3 F-S-F 133.1 90.000000 57.30 14.113457 92.917974
2 1 4 F-S-F 127.0 90.000000 57.30 87.565933 92.955537
2 1 5 F-S-F 133.2 90.000000 57.30 14.109553 92.875212
2 1 7 F-S-F 127.0 90.000000 57.30 88.486213 93.013775
3 1 4 F-S-F 134.8 90.000000 57.30 80.156500 79.481116
3 1 6 F-S-F 134.9 90.000000 57.30 55.510993 92.917974
3 1 7 F-S-F 134.8 90.000000 57.30 80.069979 79.523690
4 1 5 F-S-F 134.9 90.000000 57.30 79.565591 79.449826
4 1 6 F-S-F 128.9 90.000000 57.30 72.386283 92.955537
5 1 6 F-S-F 134.9 90.000000 57.30 55.450648 92.875212
5 1 7 F-S-F 134.9 90.000000 57.30 79.480339 79.492366
6 1 7 F-S-F 128.9 90.000000 57.30 73.014011 93.013775
2 1 6 F-S-F 120.8 180.000000 30.96 1969.376713 1969.377275
3 1 5 F-S-F 141.2 180.000000 30.96 1912.596606 1912.595899
4 1 7 F-S-F 129.5 180.000000 30.96 1942.410168 1942.410442
平均 F-S-F N/A 90.000000 57.30 64.992458 88.456000
F-S-F_Linear N/A 180.000000 30.96 1941.461162 1941.461205
py脚本所用随机向量方法与xff不同
线性键角势的力常数与常规力常数无法直接比较
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 圆球堆积接触数的估算
后一篇: 使用gnuplot绘制山体阴影图

访问人次(2015年7月 9日起): | 最后更新: 2022-12-17 03:04:57 UTC | 版权所有 © 2008 - 2022 Jerkwin