技术背景
当我们使用分子力场进行分子动力学模拟时，通常包括成键相互作用（键长项、键角项和二面角项）和非成键相互作用（范德华力）。而其中成键相互作用，可以根据作用的范围，进一步定义成1-2相互作用、1-3相互作用和1-4相互作用，也就是键长键角和二面角。在FF19SB力场中，为了进一步提升力场模拟的精确度，引入了一个CMAP参数，用于计算(C_{alpha})位置上的二面角耦合相互作用，也可以理解为1-5相互作用，作为力场的一个修正项。这里我们介绍具体的双二面角耦合势能项的计算方法。
输入与输出
具体的双二面角耦合相互作用，可以参考如下图示（图片截取自参考文献1）：
力场项的输入是(C_{alpha})周边的两个二面角：(phi)和(psi)，而我们需要计算的是这两个参数对应的一个势能修正项：(E_{mod}=E(phi,psi))。根据实验数据和量化计算出来的数据，Amber中给出了每一个氨基酸种类的二面角格点势能：
CMAP
%FLAG CMAP_COUNT 1
%FLAG CMAP_TITLE
Gly CMAP
%FLAG CMAP_RESLIST 1
GLY
%FLAG CMAP_RESOLUTION 24
%FLAG CMAP_PARAMETER
0.82366 1.09817 1.13106 1.38658 2.11377 3.63194 5.20835 4.52792
4.24857 3.42100 2.57125 2.21561 4.26659 2.21561 2.57125 3.42110
4.24858 4.52792 4.29107 3.63184 2.03731 1.32036 1.12797 1.09814
1.01976 0.89767 0.89484 1.21525 1.97095 3.70570 4.38855 4.65256
4.43948 3.89754 4.21954 2.92219 1.83300 1.61528 2.13407 3.13273
4.03304 4.37836 4.21500 3.54162 2.14871 1.44646 1.20387 1.08264
0.94294 0.79652 0.83328 1.24136 2.06053 3.43883 4.17385 4.53101
4.43115 4.03641 3.71837 2.02965 1.28863 1.21833 1.89263 2.85593
4.80291 4.23497 4.08718 3.49858 2.17811 1.53114 1.22608 1.06909
1.06124 0.90425 0.95129 1.35925 1.80900 3.18916 3.85055 4.27470
4.17515 3.68355 2.71193 1.96823 0.88051 1.02084 1.59189 2.52883
3.47767 4.03244 3.86779 3.39508 2.23969 1.66969 1.33019 1.16679
1.08157 1.02005 1.10493 1.29605 1.79270 2.88632 3.63402 4.05349
3.97657 3.35415 2.24207 1.18576 0.74267 0.67495 1.31122 2.24278
3.21905 3.69186 4.76949 4.00571 2.34764 1.77622 1.34849 1.12849
1.09360 0.76295 0.64220 1.00733 1.73708 2.50787 3.41571 4.45847
3.84910 3.11193 1.97839 1.18177 0.33803 0.43176 1.23577 2.17947
2.92428 3.38497 3.39138 3.10177 3.35228 2.15282 1.45287 1.28536
0.61764 0.23640 0.15650 0.65857 1.56522 2.44967 3.17669 3.78043
4.22216 2.79953 1.66579 0.59084 0.15579 0.43594 1.52040 2.27449
2.76979 3.16911 3.12455 2.92674 2.45824 1.88927 1.36443 0.95121
-0.00243 -0.14438 0.06931 0.84132 1.89921 2.91978 3.46170 3.60722
3.77366 2.18724 1.52870 0.20777 0.22187 1.10049 2.08960 2.88893
3.46272 3.66327 3.40280 3.02172 2.23074 1.32476 0.62376 0.25154
-0.24862 0.05886 0.73935 1.98217 3.25167 4.06285 4.21447 3.90130
2.90300 1.70374 0.74057 0.52755 1.12943 2.10802 2.78056 3.65981
4.30377 4.58476 4.35318 3.65296 2.33204 0.81283 0.00000 -0.29835
0.44379 1.32981 2.39709 3.82578 5.02459 5.42450 4.99639 4.00610
2.64886 1.63847 1.37610 1.71448 2.28503 2.64949 3.24747 4.17266
5.17265 5.75041 5.58476 4.21912 2.20540 0.64843 0.02065 -0.08593
2.15530 3.40478 4.51993 5.80241 6.12627 6.00746 5.12631 4.58197
2.80154 2.55737 2.87599 2.95502 2.71629 2.60279 3.64721 4.53964
5.90407 6.82108 6.15044 4.18011 2.15201 0.96807 0.81255 1.19704
4.21011 5.19583 5.45356 5.26216 5.22957 4.91772 4.44698 3.93537
3.74038 4.11322 3.92595 2.94210 2.11295 2.42266 3.12536 4.83315
6.43292 6.19670 5.13463 3.65052 2.40663 2.07125 2.29840 3.03351
5.33877 4.55338 3.87255 3.87071 3.62149 3.87069 4.31250 4.77933
5.37137 5.20874 3.65026 2.18763 2.20027 2.18700 3.65026 5.21405
5.24065 4.77943 4.31250 3.87553 3.63325 5.19982 3.88530 4.53465
4.20747 3.02450 2.29840 2.07135 2.40664 3.65062 5.13204 6.19680
6.43292 4.83351 3.12587 2.96667 2.11029 2.93241 3.92594 4.11322
3.74082 3.93551 4.44698 4.92068 5.53058 5.24416 5.35209 5.19605
2.15668 1.19704 0.81255 0.96806 2.15211 4.17062 6.15034 6.82108
5.90407 4.54177 3.63728 2.60267 2.71853 2.95502 2.87589 2.55737
2.80717 3.89294 5.12631 6.00746 6.12372 5.80230 5.01898 3.40679
0.45450 -0.08700 0.02065 0.64843 2.20540 4.21912 5.58476 5.75041
5.17265 4.16788 3.24747 2.64460 2.28503 1.72593 1.37610 1.63847
2.64538 4.01006 5.00425 5.42091 5.02459 3.80593 2.39681 1.32981
-0.24862 -0.29835 0.00000 0.81283 2.33214 3.64847 4.35318 4.59008
4.30377 3.65991 2.78056 2.10802 1.84494 0.52754 0.74067 1.61716
2.90300 3.89516 4.21447 4.06285 3.25166 1.98217 0.73946 0.07427
-0.00243 0.25154 0.64043 1.82987 2.23143 3.02038 3.40280 3.66327
3.46272 2.88893 2.08960 1.10049 0.22171 0.20862 1.52897 2.18595
3.24376 3.60722 3.46170 2.91968 1.89590 0.84888 0.07066 -0.14460
0.61764 0.94007 1.36444 1.88927 2.45691 2.92674 3.12469 3.16911
2.76979 2.27377 1.52040 0.44049 0.15579 0.59084 1.66579 2.79477
4.25589 3.78043 3.17669 2.45043 1.56512 0.65855 0.15650 0.23719
1.09360 1.28549 1.45287 2.15281 2.67683 3.10177 3.39138 3.38497
2.91731 2.17947 1.23582 0.43170 0.33793 1.00844 1.97383 3.11193
3.84910 4.46392 3.41571 2.50787 1.73708 1.00733 0.64220 0.76295
0.97094 1.11640 1.33640 1.77623 2.30465 3.21763 3.71128 3.69186
3.21905 2.24277 1.31122 0.67238 0.74272 1.18578 2.24207 3.35415
3.97657 4.05349 3.63266 2.88632 1.79270 1.23680 1.00391 1.05657
0.98044 1.16708 1.33038 1.66969 2.23969 3.39508 3.86779 4.03243
3.47590 2.52883 1.59193服务器托管网 1.02084 0.88051 1.46518 2.71203 3.68355
4.17470 4.27471 3.85203 3.18916 1.80900 1.35925 0.95109 0.90998
0.92229 1.06909 1.22608 1.53114 2.17810 3.49858 4.08708 4.23358
3.81032 2.85593 1.89254 1.21833 1.28863 2.03059 3.37664 4.03661
4.43079 4.53101 4.17385 3.43883 2.06053 1.24136 0.83440 0.79652
1.01976 1.08264 1.20394 1.44647 2.14870 3.54162 4.22527 4.37836
4.03304 3.13273 2.13407 1.61528 1.83300 2.92219 4.10899 3.89754
4.43948 4.65256 4.38856 3.70570 1.97095 1.21515 0.89442 0.89767
这里是一个24*24大小的Resolution表格，对应的是在空间中将两个二面角的大小划分成24份，每份为15度，而这里的值给出的是每一个离散的二面角组合的势能值。拿到这个表之后，我们需要对它做一个连续化处理，也就是使用所谓的插值法，这里选用的是三次样条插值。
三次样条插值
三次样条插值算法的原理还是比较简单粗暴的，直接根据通用公式(f(x)=ax^3+bx^2+cx+d)进行计算，那么就要求至少有四个输入输出值(x_1,y_1,x_2,y_2,x_3,y_3,x_4,y_4)用于插值：
[left(begin{matrix}
x_1^3&&x_1^2&&x_1&&1
x_2^3&&x_2^2&&x_2&&1
x_3^3&&x_3^2&&x_3&&1
x_4^3&&x_4^2&&x_4&&1
end{matrix}right)left(
begin{matrix}
abcd
end{matrix}
right)=left(
begin{matrix}
y_1y_2y_3y_4
end{matrix}
right)
]
x_1^3&&x_1^2&&x_1&&1
x_2^3&&x_2^2&&x_2&&1
x_3^3&&x_3^2&&x_3&&1
x_4^3&&x_4^2&&x_4&&1
end{matrix}right)left(
begin{matrix}
abcd
end{matrix}
right)=left(
begin{matrix}
y_1y_2y_3y_4
end{matrix}
right)
]
这样就可以通过计算矩阵的逆来得到一组插值参数：
[left(
begin{matrix}
abcd
end{matrix}
right)=
left(begin{matrix}
x_1^3&&x_1^2&&x_1&&1
x_2^3&&x_2^2&&x_2&&1
x_3^3&&x_3^2&&x_3&&1
x_4^3&&x_4^2&&x_4&&1
end{matrix}right)^{-1}left(
begin{matrix}
y_1y_2y_3y_4
end{matrix}
right)
]
begin{matrix}
abcd
end{matrix}
right)=
left(begin{matrix}
x_1^3&&x_1^2&&x_1&&1
x_2^3&&x_2^2&&x_2&&1
x_3^3&&x_3^2&&x_3&&1
x_4^3&&x_4^2&&x_4&&1
end{matrix}right)^{-1}left(
begin{matrix}
y_1y_2y_3y_4
end{matrix}
right)
]
如此一来我们就把一组离散的数据插值成一个可以求导的连续函数：
[f(x)=ax^3+bx^2+cx+d
frac{df}{dx}=3ax^2+2bx+c
]
frac{df}{dx}=3ax^2+2bx+c
]
双三次插值
回到双二面角耦合相互租用的场景，我们按照15度一个点，对整个二面角的自变量空间做了一个离散化，可以先用可视化的形式看一下对应的能量分布。
import numpy as np
import matplotlib.pyplot as plt
string = """
0.82366 1.09817 1.13106 1.38658 2.11377 3.63194 5.20835 4.52792
4.24857 3.42100 2.57125 2.21561 4.26659 2.21561 2.57125 3.42110
4.24858 4.52792 4.29107 3.63184 2.03731 1.32036 1.12797 1.09814
1.01976 0.89767 0.89484 1.21525 1.97095 3.70570 4.38855 4.65256
4.43948 3.89754 4.21954 2.92219 1.83300 1.61528 2.13407 3.13273
4.03304 4.37836 4.21500 3.54162 2.14871 1.44646 1.20387 1.08264
0.94294 0.79652 0.83328 1.24136 2.06053 3.43883 4.17385 4.53101
4.43115 4.03641 3.71837 2.02965 1.28863 1.21833 1.89263 2.85593
4.80291 4.23497 4.08718 3.49858 2.17811 1.53114 1.22608 1.06909
1.06124 0.90425 0.95129 1.35925 1.80900 3.18916 3.85055 4.27470
4.17515 3.68355 2.71193 1.96823 0.88051 1.02084 1.59189 2.52883
3.47767 4.03244 3.86779 3.39508 2.23969 1.66969 1.33019 1.16679
1.08157 1.02005 1.10493 1.29605 1.79270 2.88632 3.63402 服务器托管网 4.05349
3.97657 3.35415 2.24207 1.18576 0.74267 0.67495 1.31122 2.24278
3.21905 3.69186 4.76949 4.00571 2.34764 1.77622 1.34849 1.12849
1.09360 0.76295 0.64220 1.00733 1.73708 2.50787 3.41571 4.45847
3.84910 3.11193 1.97839 1.18177 0.33803 0.43176 1.23577 2.17947
2.92428 3.38497 3.39138 3.10177 3.35228 2.15282 1.45287 1.28536
0.61764 0.23640 0.15650 0.65857 1.56522 2.44967 3.17669 3.78043
4.22216 2.79953 1.66579 0.59084 0.15579 0.43594 1.52040 2.27449
2.76979 3.16911 3.12455 2.92674 2.45824 1.88927 1.36443 0.95121
-0.00243 -0.14438 0.06931 0.84132 1.89921 2.91978 3.46170 3.60722
3.77366 2.18724 1.52870 0.20777 0.22187 1.10049 2.08960 2.88893
3.46272 3.66327 3.40280 3.02172 2.23074 1.32476 0.62376 0.25154
-0.24862 0.05886 0.73935 1.98217 3.25167 4.06285 4.21447 3.90130
2.90300 1.70374 0.74057 0.52755 1.12943 2.10802 2.78056 3.65981
4.30377 4.58476 4.35318 3.65296 2.33204 0.81283 0.00000 -0.29835
0.44379 1.32981 2.39709 3.82578 5.02459 5.42450 4.99639 4.00610
2.64886 1.63847 1.37610 1.71448 2.28503 2.64949 3.24747 4.17266
5.17265 5.75041 5.58476 4.21912 2.20540 0.64843 0.02065 -0.08593
2.15530 3.40478 4.51993 5.80241 6.12627 6.00746 5.12631 4.58197
2.80154 2.55737 2.87599 2.95502 2.71629 2.60279 3.64721 4.53964
5.90407 6.82108 6.15044 4.18011 2.15201 0.96807 0.81255 1.19704
4.21011 5.19583 5.45356 5.26216 5.22957 4.91772 4.44698 3.93537
3.74038 4.11322 3.92595 2.94210 2.11295 2.42266 3.12536 4.83315
6.43292 6.19670 5.13463 3.65052 2.40663 2.07125 2.29840 3.03351
5.33877 4.55338 3.87255 3.87071 3.62149 3.87069 4.31250 4.77933
5.37137 5.20874 3.65026 2.18763 2.20027 2.18700 3.65026 5.21405
5.24065 4.77943 4.31250 3.87553 3.63325 5.19982 3.88530 4.53465
4.20747 3.02450 2.29840 2.07135 2.40664 3.65062 5.13204 6.19680
6.43292 4.83351 3.12587 2.96667 2.11029 2.93241 3.92594 4.11322
3.74082 3.93551 4.44698 4.92068 5.53058 5.24416 5.35209 5.19605
2.15668 1.19704 0.81255 0.96806 2.15211 4.17062 6.15034 6.82108
5.90407 4.54177 3.63728 2.60267 2.71853 2.95502 2.87589 2.55737
2.80717 3.89294 5.12631 6.00746 6.12372 5.80230 5.01898 3.40679
0.45450 -0.08700 0.02065 0.64843 2.20540 4.21912 5.58476 5.75041
5.17265 4.16788 3.24747 2.64460 2.28503 1.72593 1.37610 1.63847
2.64538 4.01006 5.00425 5.42091 5.02459 3.80593 2.39681 1.32981
-0.24862 -0.29835 0.00000 0.81283 2.33214 3.64847 4.35318 4.59008
4.30377 3.65991 2.78056 2.10802 1.84494 0.52754 0.74067 1.61716
2.90300 3.89516 4.21447 4.06285 3.25166 1.98217 0.73946 0.07427
-0.00243 0.25154 0.64043 1.82987 2.23143 3.02038 3.40280 3.66327
3.46272 2.88893 2.08960 1.10049 0.22171 0.20862 1.52897 2.18595
3.24376 3.60722 3.46170 2.91968 1.89590 0.84888 0.07066 -0.14460
0.61764 0.94007 1.36444 1.88927 2.45691 2.92674 3.12469 3.16911
2.76979 2.27377 1.52040 0.44049 0.15579 0.59084 1.66579 2.79477
4.25589 3.78043 3.17669 2.45043 1.56512 0.65855 0.15650 0.23719
1.09360 1.28549 1.45287 2.15281 2.67683 3.10177 3.39138 3.38497
2.91731 2.17947 1.23582 0.43170 0.33793 1.00844 1.97383 3.11193
3.84910 4.46392 3.41571 2.50787 1.73708 1.00733 0.64220 0.76295
0.97094 1.11640 1.33640 1.77623 2.30465 3.21763 3.71128 3.69186
3.21905 2.24277 1.31122 0.67238 0.74272 1.18578 2.24207 3.35415
3.97657 4.05349 3.63266 2.88632 1.79270 1.23680 1.00391 1.05657
0.98044 1.16708 1.33038 1.66969 2.23969 3.39508 3.86779 4.03243
3.47590 2.52883 1.59193 1.02084 0.88051 1.46518 2.71203 3.68355
4.17470 4.27471 3.85203 3.18916 1.80900 1.35925 0.95109 0.90998
0.92229 1.06909 1.22608 1.53114 2.17810 3.49858 4.08708 4.23358
3.81032 2.85593 1.89254 1.21833 1.28863 2.03059 3.37664 4.03661
4.43079 4.53101 4.17385 3.43883 2.06053 1.24136 0.83440 0.79652
1.01976 1.08264 1.20394 1.44647 2.14870 3.54162 4.22527 4.37836
4.03304 3.13273 2.13407 1.61528 1.83300 2.92219 4.10899 3.89754
4.43948 4.65256 4.38856 3.70570 1.97095 1.21515 0.89442 0.89767
"""
phi = np.arange(24) * 15
psi = np.arange(24) * 15
resolutions = np.array([float(s) for s in string.split()], np.float32).reshape((24, 24))
plt.figure(figsize=(10, 10))
plt.title('GLY Resolutions')
plt.contourf(phi, psi, resolutions, cmap=plt.cm.hot)
plt.xlabel(r'$phi$')
plt.ylabel(r'$psi$')
plt.xticks(phi)
plt.yticks(psi)
plt.grid()
plt.show()
得到的结果如下图所示：
此时需要注意的是，我们手头仅有每一个格点上的对应能量值，但是我们需要的是整个自变量空间的能量值，而且需要是连续的，这样才能够进行自动微分，进而求解对应的力。这时候需要用到双三次样条插值法。我们假设需要插值的点处在某一个方格子中间，注意不是格点了，是在某一个方格子的中间。那么此时我们可以以这个方格的位置为中心，构造出一个九宫格。虽然是九宫格，但其实是4横4纵一共16个能量点。这里给出的插值目标函数为：
[E_{interp}=sum_{i,j=1}^4c_{i,j}left(frac{phi-phi_L}{Deltaphi}right)^{i-1}left(frac{psi-psi_L}{Deltapsi}right)^{j-1}
]
]
由于系数我们是当做一维来计算的，所以需要做一个展开令(a_{J}=a_{4i+j-4}=c_{i,j})，同时，我们可以将角度的变换关系全都放到系数项里面，这样表示(phi)和(psi)时在形式上就会更加简单。经过变换后会得到一个矩阵方程：
[sum_{J}M_{I,J}a_J=E_I
]
]
其中(M_{I,J})可以通过格点的(phi,psi)与中心方格的中心点的(phi_L,psi_L)构造出来，这里不展开写公式，详细内容参见代码实现。那么最后得到系数矩阵(M_{I,J})和能量表之后，就可以进一步得到需要插值的系数：
[a=M^{-1}E
]
]
具体的代码实现如下：
M = np.zeros((16, 16))
for i in range(16):
for j in range(16):
ii = i // 4
ij = i % 4
ji = j // 4
jj = j % 4
M[i][j] = (-1.5 + ij) ** ji * (-1.5 + ii) ** jj
M = np.matrix(M, np.float32)
c_0 = np.einsum('ij,j->i', M.I, E_sub.reshape(-1))
这样就可以得到一组参数：
[ 0.80195737 -0.09009341 0.089855 0.00695395 0.00651729 0.04415806
0.00911433 -0.01563221 0.1316725 0.02541727 -0.02819001 0.01033084
0.03029101 -0.00083557 0.00954252 -0.00721776]
根据这个参数，可以代入到上面计算势能的插值函数(
E_{interp}=sum_{i,j=1}^4c_{i,j}left(frac{phi-phi_L}{Deltaphi}right)^{i-1}left(frac{psi-psi_L}{Deltapsi}right)^{j-1}
)中，计算出格点内给定位置的势能值。完整的(phi-psi)空间参数计算Python代码如下所示：
import numpy as np
from itertools import product
import matplotlib.pyplot as plt
string = """
0.82366 1.09817 1.13106 1.38658 2.11377 3.63194 5.20835 4.52792
4.24857 3.42100 2.57125 2.21561 4.26659 2.21561 2.57125 3.42110
4.24858 4.52792 4.29107 3.63184 2.03731 1.32036 1.12797 1.09814
1.01976 0.89767 0.89484 1.21525 1.97095 3.70570 4.38855 4.65256
4.43948 3.89754 4.21954 2.92219 1.83300 1.61528 2.13407 3.13273
4.03304 4.37836 4.21500 3.54162 2.14871 1.44646 1.20387 1.08264
0.94294 0.79652 0.83328 1.24136 2.06053 3.43883 4.17385 4.53101
4.43115 4.03641 3.71837 2.02965 1.28863 1.21833 1.89263 2.85593
4.80291 4.23497 4.08718 3.49858 2.17811 1.53114 1.22608 1.06909
1.06124 0.90425 0.95129 1.35925 1.80900 3.18916 3.85055 4.27470
4.17515 3.68355 2.71193 1.96823 0.88051 1.02084 1.59189 2.52883
3.47767 4.03244 3.86779 3.39508 2.23969 1.66969 1.33019 1.16679
1.08157 1.02005 1.10493 1.29605 1.79270 2.88632 3.63402 4.05349
3.97657 3.35415 2.24207 1.18576 0.74267 0.67495 1.31122 2.24278
3.21905 3.69186 4.76949 4.00571 2.34764 1.77622 1.34849 1.12849
1.09360 0.76295 0.64220 1.00733 1.73708 2.50787 3.41571 4.45847
3.84910 3.11193 1.97839 1.18177 0.33803 0.43176 1.23577 2.17947
2.92428 3.38497 3.39138 3.10177 3.35228 2.15282 1.45287 1.28536
0.61764 0.23640 0.15650 0.65857 1.56522 2.44967 3.17669 3.78043
4.22216 2.79953 1.66579 0.59084 0.15579 0.43594 1.52040 2.27449
2.76979 3.16911 3.12455 2.92674 2.45824 1.88927 1.36443 0.95121
-0.00243 -0.14438 0.06931 0.84132 1.89921 2.91978 3.46170 3.60722
3.77366 2.18724 1.52870 0.20777 0.22187 1.10049 2.08960 2.88893
3.46272 3.66327 3.40280 3.02172 2.23074 1.32476 0.62376 0.25154
-0.24862 0.05886 0.73935 1.98217 3.25167 4.06285 4.21447 3.90130
2.90300 1.70374 0.74057 0.52755 1.12943 2.10802 2.78056 3.65981
4.30377 4.58476 4.35318 3.65296 2.33204 0.81283 0.00000 -0.29835
0.44379 1.32981 2.39709 3.82578 5.02459 5.42450 4.99639 4.00610
2.64886 1.63847 1.37610 1.71448 2.28503 2.64949 3.24747 4.17266
5.17265 5.75041 5.58476 4.21912 2.20540 0.64843 0.02065 -0.08593
2.15530 3.40478 4.51993 5.80241 6.12627 6.00746 5.12631 4.58197
2.80154 2.55737 2.87599 2.95502 2.71629 2.60279 3.64721 4.53964
5.90407 6.82108 6.15044 4.18011 2.15201 0.96807 0.81255 1.19704
4.21011 5.19583 5.45356 5.26216 5.22957 4.91772 4.44698 3.93537
3.74038 4.11322 3.92595 2.94210 2.11295 2.42266 3.12536 4.83315
6.43292 6.19670 5.13463 3.65052 2.40663 2.07125 2.29840 3.03351
5.33877 4.55338 3.87255 3.87071 3.62149 3.87069 4.31250 4.77933
5.37137 5.20874 3.65026 2.18763 2.20027 2.18700 3.65026 5.21405
5.24065 4.77943 4.31250 3.87553 3.63325 5.19982 3.88530 4.53465
4.20747 3.02450 2.29840 2.07135 2.40664 3.65062 5.13204 6.19680
6.43292 4.83351 3.12587 2.96667 2.11029 2.93241 3.92594 4.11322
3.74082 3.93551 4.44698 4.92068 5.53058 5.24416 5.35209 5.19605
2.15668 1.19704 0.81255 0.96806 2.15211 4.17062 6.15034 6.82108
5.90407 4.54177 3.63728 2.60267 2.71853 2.95502 2.87589 2.55737
2.80717 3.89294 5.12631 6.00746 6.12372 5.80230 5.01898 3.40679
0.45450 -0.08700 0.02065 0.64843 2.20540 4.21912 5.58476 5.75041
5.17265 4.16788 3.24747 2.64460 2.28503 1.72593 1.37610 1.63847
2.64538 4.01006 5.00425 5.42091 5.02459 3.80593 2.39681 1.32981
-0.24862 -0.29835 0.00000 0.81283 2.33214 3.64847 4.35318 4.59008
4.30377 3.65991 2.78056 2.10802 1.84494 0.52754 0.74067 1.61716
2.90300 3.89516 4.21447 4.06285 3.25166 1.98217 0.73946 0.07427
-0.00243 0.25154 0.64043 1.82987 2.23143 3.02038 3.40280 3.66327
3.46272 2.88893 2.08960 1.10049 0.22171 0.20862 1.52897 2.18595
3.24376 3.60722 3.46170 2.91968 1.89590 0.84888 0.07066 -0.14460
0.61764 0.94007 1.36444 1.88927 2.45691 2.92674 3.12469 3.16911
2.76979 2.27377 1.52040 0.44049 0.15579 0.59084 1.66579 2.79477
4.25589 3.78043 3.17669 2.45043 1.56512 0.65855 0.15650 0.23719
1.09360 1.28549 1.45287 2.15281 2.67683 3.10177 3.39138 3.38497
2.91731 2.17947 1.23582 0.43170 0.33793 1.00844 1.97383 3.11193
3.84910 4.46392 3.41571 2.50787 1.73708 1.00733 0.64220 0.76295
0.97094 1.11640 1.33640 1.77623 2.30465 3.21763 3.71128 3.69186
3.21905 2.24277 1.31122 0.67238 0.74272 1.18578 2.24207 3.35415
3.97657 4.05349 3.63266 2.88632 1.79270 1.23680 1.00391 1.05657
0.98044 1.16708 1.33038 1.66969 2.23969 3.39508 3.86779 4.03243
3.47590 2.52883 1.59193 1.02084 0.88051 1.46518 2.71203 3.68355
4.17470 4.27471 3.85203 3.18916 1.80900 1.35925 0.95109 0.90998
0.92229 1.06909 1.22608 1.53114 2.17810 3.49858 4.08708 4.23358
3.81032 2.85593 1.89254 1.21833 1.28863 2.03059 3.37664 4.03661
4.43079 4.53101 4.17385 3.43883 2.06053 1.24136 0.83440 0.79652
1.01976 1.08264 1.20394 1.44647 2.14870 3.54162 4.22527 4.37836
4.03304 3.13273 2.13407 1.61528 1.83300 2.92219 4.10899 3.89754
4.43948 4.65256 4.38856 3.70570 1.97095 1.21515 0.89442 0.89767
"""
def E_interp(c, phi, psi):
# 根据计算好的参数以及输入的phi和psi值，计算势能值
rc = c.reshape((4, 4))
E = 0
for i in range(4):
for j in range(4):
E += rc[i][j] * phi ** i * psi ** j
return E
def get_c0(E_sub):
# 计算插值参数
M = np.zeros((16, 16))
for i in range(16):
for j in range(16):
ii = i // 4
ij = i % 4
ji = j // 4
jj = j % 4
M[i][j] = (-1.5 + ii) ** ji * (-1.5 + ij) ** jj
M = np.matrix(M, np.float32)
c_0 = np.einsum('ij,j->i', M.I, E_sub.reshape(-1))
return c_0
resolutions = np.array([float(s) for s in string.split()], np.float32).reshape((24, 24))
def get_resolution_0(resolutions, phi_L=0., psi_L=0.):
# 遍历空间中每个网格的中心，计算势能值
Einterp = np.zeros((23, 23))
for i in range(23):
for j in range(23):
idx = list(product([i-1, i, i+1, (i+2) % 24], [j-1, j, i+1, (j+2) % 24]))
E_sub = np.array([resolutions[x] for x in idx]).reshape((4, 4))
c0 = get_c0(E_sub)
Einterp[i][j] = E_interp(c0, phi_L, psi_L)
return Einterp
phi = np.arange(23) * 15
psi = np.arange(23) * 15
hi, si = 0., 0.
# 网格中心的势能分布
resolutions_0 = get_resolution_0(resolutions, hi, si)
hi, si = -0.5, -0.5
# 网格格点的势能分布
resolutions_1 =get_resolution_0(resolutions, hi, si)
# 绘图
plt.figure(figsize=(20, 10))
plt.subplot(1,2,1)
plt.title('GLY Resolutions Center')
plt.contourf(phi+7.5, psi+7.5, resolutions_0, cmap=plt.cm.hot)
plt.xlabel(r'$phi$')
plt.ylabel(r'$psi$')
plt.xticks(phi+7.5, rotation=60)
plt.yticks(psi+7.5)
plt.grid()
plt.subplot(1,2,2)
plt.title('GLY Resolutions Grid')
plt.contourf(phi, psi, resolutions_1, cmap=plt.cm.hot)
plt.xlabel(r'$phi$')
plt.ylabel(r'$psi$')
plt.xticks(phi, rotation=60)
plt.yticks(psi)
plt.grid()
plt.show()
输出的结果如下所示：
将这个结果对比原始的Resolution，我们发现格点上的势能分布是完全一致的。但是如果使用的是网格中心的势能值用来统计分布，势能面就会出现较大的差异，但其实我们做插值的时候，也要考虑是否存在“过拟合”的问题。
Amber算法复现
在Amber里面，就是直接使用双三次样条插值进行计算，这里不对其算法进行展开介绍，但是可以贴一部分用于复现的代码：
def get_c1(E_sub):
dM = np.matrix([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[1., 0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0.],
[1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 1., 0., 0., 0., 2., 0., 0., 0., 3., 0., 0., 0.],
[0., 0., 0., 0., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 1., 1., 1., 1., 2., 2., 2., 2., 3., 3., 3., 3.],
[0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 1., 0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1., 0., 0.],
[0., 1., 2., 3., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 1., 2., 3., 0., 1., 2., 3., 0., 1., 2., 3., 0., 1., 2., 3.],
[0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 0., 0., 0., 2., 0., 0., 0., 3., 0., 0.],
[0., 0., 0., 0., 0., 1., 2., 3., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 2., 3., 0., 2., 4., 6., 0., 3., 6., 9.]],
np.float32)
dE = np.array([E_sub[1][1], E_sub[2][1], E_sub[1][2], E_sub[2][2],
0.5 * (E_sub[2][1] - E_sub[0][1]), 0.5 * (E_sub[3][1] - E_sub[1][1]),
0.5 * (E_sub[2][2] - E_sub[0][2]), 0.5 * (E_sub[3][2] - E_sub[1][2]),
0.5 * (E_sub[1][2] - E_sub[1][0]), 0.5 * (E_sub[2][2] - E_sub[2][0]),
0.5 * (E_sub[1][3] - E_sub[1][1]), 0.5 * (E_sub[2][3] - E_sub[2][1]),
0.25 * (E_sub[2][2] + E_sub[0][0] - E_sub[2][0] - E_sub[0][2]),
0.25 * (E_sub[3][2] + E_sub[1][0] - E_sub[3][0] - E_sub[1][2]),
0.25 * (E_sub[2][2] + E_sub[0][2] - E_sub[2][1] - E_sub[0][3]),
0.25 * (E_sub[3][3] + E_sub[1][1] - E_sub[3][1] - E_sub[1][3])],
np.float32)
c_1 = np.einsum('ij,j->i', dM.I, dE)
return c_1
这个双三次样条插值算法的基本原理，不是像之前我们提到的一样使用格点中心点作为中心来进行插值，而是选取了中心网格的左上角格点作为中心(0,0)，然后对中间网格的四个格点、以及边的一阶差分、二阶差分进行插值。这样做的好处是，在网格的边缘处，不仅仅是势能面光滑，且一阶导数可以连续。
总结概要
本文介绍了最新的一些分子力场中有可能使用到的1-5相互作用——双二面角耦合项的计算。而常规的计算方式是，通过量化的计算得到每一个Residue的碳对应的两个二面角的数值在空间中的离散化数值。然后在分子模拟的过程中使用插值的方案，对相关的条目进行计算，例如使用双三次样条插值。但其实这种插值算法的使用有可能导致一些其他的问题，比如深度学习中可能会经常提到的“过拟合”问题。由于这个条目本身就只是一个修正项，从数值大小上来说并不算大，因此这些细节是否需要优化还有待商榷。
版权声明
本文首发链接为：https://www.cnblogs.com/dechinphy/p/cmap.html
作者ID：DechinPhy
更多原著文章：https://www.cnblogs.com/dechinphy/
请博主喝咖啡：https://www.cnblogs.com/dechinphy/gallery/image/379634.html
参考文献
- Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-Phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations.
服务器托管，北京服务器托管，服务器租用 http://www.fwqtg.net
机房租用，北京机房租用，IDC机房托管， http://www.fwqtg.net
相关推荐: [数据校验/数据质量] 数据校验框架(Java)：hibernate-validation
0 前言 其一，项目中普遍遇到了此问题，故近两天深入地研究了一下。 其二，能够自信地说，仔细看完本篇，就无需再看其他的Java数据校验框架的文章了。 1 数据校验框架概述 1.0 数据校验框架的产生背景 以Web项目为例，用户需要填写表单信息保存提交。 页面输…
