求助，GAFF力场中二面角计算时无法重现官网势能曲线

liwl926

关于二面角的计算问题，详细学习了Amber力场二面角参数的解释http://sobereva.com/1，并查看了sobtop里bonded_param.data,但是仍无法重现https://ambermd.org/antechamber/gaff.html#test1中二面角的函数中Figure 4. An example: rotation profiles of X-c-c-X的曲线。sobtop中公式：E_dihedral = sum[ (PK/IDIVF) * (1 + cos(PN*phi - phase)) ]；参数：X -c -c -X    4    1.200       180.000           2.000，当phi=180时，E_dihedral=0，而当phi=0时，E_dihedral值依然为0，这与amber主页二面角的函数描述的Figure 4不符，更与实际不符，一定是什么地方理解有误. 也尝试了对NP=0，1，2求和，依然不对。求教，非常感谢！

sobereva

看GAFF力场原文图1可知，这个是对于文中图4里1号分子而言的，旋转二面角的时候牵扯到分子里不止一个二面角项。图中的应当是总的效果，而不仅仅是一项

liwl926

sobereva 发表于 2024-5-16 21:08
看GAFF力场原文图1可知，这个是对于文中图4里1号分子而言的，旋转二面角的时候牵扯到分子里不止一个二面角 ...

谢谢sob老师提示，今天又学习了一天。但按照公式E_dihedral = sum[ (PK/IDIVF) * (1 + cos(PN*phi - phase)) ]和参数，即便考虑全部四组二面角，仍无法获得GAFF力场原文图3（即官网图4）曲线，而且如果考虑全部二面角都扭转180°，在角度为0°时应该是另一个极小点，而不是图中的能量最高点（约5.4 kcal/mol）。求指教！持续感谢！！

sobereva

liwl926 发表于 2024-5-17 15:19
谢谢sob老师提示，今天又学习了一天。但按照公式E_dihedral = sum[ (PK/IDIVF) * (1 + cos(PN*phi - phas ...

你可以先用Gaussian做势能面扫描，尝试复现一下文中的量子化学的势能面扫描结果。
如果成功了，再去用gromacs计算旋转到不同二面角的结构的势能，如果也和图吻合了，再一个个删去二面角项，反复对比，早晚能搞清楚原因。