求助: 如何使用gromacs对tip4p水模型进行normal mode分析

Lacrimosa

使用3点水模型做出来的结果还算正常，同样的方法用在四点水模型上后出现了大量虚频。恳请各位老师指点。
体系很简单，就是一个900个分子的水盒子，力场用的tip4p/2005，npt跑平衡后做了能量最小化，结果如下：
Potential Energy  = -5.41678855181009e+04
Maximum force     =  9.08945716315534e-03 on atom 1865
Norm of force     =  4.01353831285450e-04

normal mode使用的mdp如下：
define = -DFLEXIBLE
integrator = nm
dt                       = 0.002     ;
nsteps                   = 10000   ;
nstenergy                = 5000
nstlog                   = 5000
nstxout-compressed       = 2000
continuation             = yes
constraints              = none
cutoff-scheme            = Verlet
coulombtype              = PME
vdwtype                  = Cut-off
DispCorr                 = EnerPres
emtol = 0.01
emstep = 0.01
nstlist = 5
rlist = 1.0
rvdw = 1.0
rcoulomb = 1.0


计算结果大致如下：
@    title "Eigenvalues"
@    xaxis  label "Eigenvalue index"
@    yaxis  label "Eigenvalue [Gromacs units]"
@TYPE xy
@ subtitle "mass weighted"
     7    -1.57368e+11
     8     -1.3606e+11
     9    -1.20261e+11
    10    -1.06565e+11
    11    -9.77809e+10
    12     -8.9158e+10

sobereva

尝试用双精度版本LBFGS方法做em以及正则模式分析。如果也有问题，而所有三点水模型都没问题，估计是gmx在振动分析时对虚拟点考虑不周。也可以试试关掉色散校正并且库仑作用用cutoff算

Lacrimosa

sobereva 发表于 2023-8-23 07:34
尝试用双精度版本LBFGS方法做em以及正则模式分析。如果也有问题，而所有三点水模型都没问题，估计是gmx在振 ...

感谢老师回复，根据您的建议做了一些尝试还是解决不了，可能确实是虚位点的问题了。请问有什么其他的程序可以替代gromacs实现这个分析么？

Lacrimosa

自答一下:
除正则分析之外，还可以通过对质量权重的速度自相关函数(mass-weighted velocity autocorrelation function, 下文简写为mw-vacf)进行傅里叶变换来获得VDOS，小体系低温下的测试结果与正则分析基本吻合
通过此文 J. Chem. Phys. 150, 194111 (2019)中的eq(22)以及eq(23)可以计算体系的自由度，其数量应为3N，其中N为体系中的总原子数。实际计算结果会与3N略有差异。
需要注意以下几点：
#1. gmx velacc只能计算vacf，使用eq(22)需要对每个原子进行计算然后乘以其相应的质量后求和，按照此流程算出的mw-vacf结果与用gmx dos算出的不同（此处使用的是gmx 2021.3），暂时还无法确定gmx dos是如何计算mw-vacf的。
#2. 该方法基于harmonic approximation，体系应在尽可能低的温度下进行模拟。
#3. 体积对计算结果有着明显的影响，如要将正则分析的结果与该方法的结果进行对比，应在相同体积下进行模拟。

对于#1，也可以对每一种原子计算vacf, 然后乘以其质量以及数量后对所有种类求和