求助：使用文献参考力场模拟液体发现密度远低于标准值

寻觅者

in文件如下
#################### 基础设置 ####################

units               real                       # 真实单位
atom_style          full                       
boundary            p p p                      # 周期性边界
pair_style          lj/cut/coul/long 14.0       # long range cut off 14.0  
pair_modify         mix arithmetic tail yes    # mix arithmetic=Lorentz-Berthelot,tail yes=ATC
kspace_style        pppm 0.0001                # long range force pppm
bond_style            harmonic
angle_style            harmonic
dihedral_style      harmonic
improper_style      cvff

#################### 原子坐标及力场文件读取 ####################

read_data           PureNMP.data
include             forcefield_PureNMP.data

#################### 变量及组设定 ####################

variable            Temp equal 298                            # Temperature in K
variable            Pres equal 1                       # Pressure in atm
variable            RandomSeed equal 114514            # The random seed for velocity

variable            NiniNVT equal 100000               # Initialize the NVT ensemble
variable            NequNPT equal 1000000              # Equilibrium and product in the NPT ensemble

variable            Nthermol equal 100000              # 100ps看一次热力学数据

#################### NVT系综平衡 ####################

neighbor            2.0 bin
neigh_modify            every 1 delay 0 check no

run_style           verlet

min_style           cg                                 #共轭梯度法
min_modify          dmax 0.05
minimize            1.0e-6 1.0e-6 1000 10000

reset_timestep 0

velocity            all create 50.0 ${RandomSeed}

thermo              ${Nthermol}
thermo_style        custom step temp etotal density press vol
thermo_modify       flush yes

#################### NVT系综初始化 ####################
fix                 NvtEns all nvt temp 50.0 ${Temp} 100.0

run                 ${NiniNVT}
reset_timestep 0

#################### NPT系综平衡 ####################

unfix               NvtEns

fix                 NptEns all npt temp ${Temp} ${Temp} 100.0 iso ${Pres} ${Pres} 1000.0

run                 ${NequNPT}


力场选用文献中推荐的GAFF力场，参数来自AMBERtools21中的gaff.dat文件，但是文献中算得密度是1.018g/cm3，而我算得密度只有0.90g/cm3，不知道错误原因在何处，还请各位大神指点一二

ChemiAndy

请问你的初始态是如何构造的？如何确保模拟已经收敛？是否进行了多样本统计平均处理？

液态体系的密度的分子模拟估算其实并不简单，真实液体的密度是体系经过各种可能构型之后的一个统计平均值，而且体系呆在某种构型的概率符合玻尔兹曼分布，即与此构型的能量呈反指数关系。进行分子模拟时，由于温度通常较低，很难通过自然演进在模拟的可及范围内找到模拟体系真正的最低能构象。绝大多数模拟其实从一开始就呆在初始态附近的一个低能构型附近晃悠，你得到的结果自然也就与真实值相去甚远，更遑论力场质量对结果的扭曲了。

在不谈论力场质量前提下，一般模拟液体的密度（实质是寻找最大概然构型），要从多个完全不同的初始构型出发，得到多个收敛的模拟的密度结果然后取平均。一般来说最少要构建3个不同的初始构型。构造初始构型的方法，最好的是蒙特卡洛方法；其次是分子动力学淬火方法，让体系升高到很高，彻底远离初始构型，再祈祷它随机落到一个理想的、与初始构型完全无关的构型上去。然后再以该构型出发去做一次密度的估算，依次往复。。。

hdhxx123

ChemiAndy 发表于 2021-7-1 11:57
请问你的初始态是如何构造的？如何确保模拟已经收敛？是否进行了多样本统计平均处理？

液态体系的密度的 ...

“祈祷”这个用词好好玩

Kangtor

ChemiAndy 发表于 2021-7-1 11:57
请问你的初始态是如何构造的？如何确保模拟已经收敛？是否进行了多样本统计平均处理？

液态体系的密度的 ...

哈哈哈

Kangtor

你怎么跑的动力学？

寻觅者

ChemiAndy 发表于 2021-7-1 11:57
请问你的初始态是如何构造的？如何确保模拟已经收敛？是否进行了多样本统计平均处理？

液态体系的密度的 ...

老师您好
我的初始态是直接用packmol构建的，分子在15nm盒子里平均分布，最后平衡了有1ns，在这1ns里我的体系总能量几乎稳定了，关于多样本统计，我只改变了初始温度和压力，而这些样本最终的密度平衡结果不变，均为0.9左右
问题在于，我用msi2lmp得到的初始模型，使用cvff力场参数，在相同的情况下能够得到密度约为0.98，非常接近实验值，这是非常令我不解的，因为cvff在模拟有机分子的精度没有gaff高，而且gaff力场是有文献算得的1.016
按照老师您的方法，那我是不是应该重新构建初始态，多次进行模拟，然后取结果的平均值？

寻觅者

Kangtor 发表于 2021-7-1 14:49
你怎么跑的动力学？

我的in文件已经贴在上面呢，就是让体系驰豫到标准态，用NPT系综控制

ChemiAndy

你这个模拟有两个问题：

1. 首先Packmol是基于纯几何方式在盒子里面堆积分子。存在相当大的概率导致盒子内的分子排列在能量上是极不合理的。因此使用Packmol搭建的盒子，最好进行NVT高温淬火处理，使得体系能够有较大概率跑到附近的低能构型上去。

2. 初始密度问题。你用较大的盒子往里塞分子，导致初始盒子的密度低于实验值，然后期望NPT在1ns内将盒子shrink到理想的值附近，太naive了。构型弛豫的效率是极低的，特别是在298K这样的温度下，1ns听上去很长，对液体体系来讲却只是转瞬而已。而且你NVT松弛的不充分，体系仍然在你一开始搭建的那个不合理构型附近辗转徘徊。

当然，以上分析都是general的，我并不了解你的体系的特点。有的液体体系粘度很低，比如大部分非极性小分子液体体系，构型弛豫相对容易；有的体系粘度很大，比如极性分子液体体系和高分子体系，这样的体系的模拟是极大的挑战。要读文献看看人家是怎么克服采样困难的。

你的模拟，建议从实验密度值开始搭建盒子，之后NVT高温淬火（1000K/数百个PS），然后在298K下NVT能量稳定后，开始NPT。至少500ps且确定能量/体积稳定后，开始Production run，1ns。重复三次取平均值。 注意每次运行Packmol都要更改随机数种子，否则得到的是一样的构型。

寻觅者

ChemiAndy 发表于 2021-7-2 04:40
你这个模拟有两个问题：

1. 首先Packmol是基于纯几何方式在盒子里面堆积分子。存在相当大的概率导致盒子 ...

好的，谢谢老师！