首页 >> 仿真软件说明 >>gromacs >>mdp >> integrator
详细内容

integrator

integrator

翻译虽然叫积分器,但是该方法是用于处理计算运动的方法


steep: 能量最小化的最速下降算法,根据最陡路径下到了最低处再找最陡路径,垂直等高线的方向,一开始速度快,后面的难收敛

配合:

emtol:容差,单位为 [kJ mol-1 nm-1],设置为emtol=100

emstep:最大步长,单位为 [nm],设置为emstep=0.01


cg: 能量最小化的共轭梯度算法,结合当前的梯度与上一步结合的方向进行,两个梯度进行结合

配合:

emtol:容差,单位为 [kJ mol-1 nm-1],设置为emtol=100

emstep:最大步长,单位为 [nm],设置为emstep=0.01


image.png

md: 蛙跳法,对符合牛顿公式的运动进行积分。

md-vv: 速度Verlet算法, 用于积分运动方程,当需要非常精确的积分方法并采用温度和/或压力耦合时, 使用速度Verlet积分方法可能更好。注意,对于几乎所有的生产模拟,md 积分器都足够精确。

md-vv-avek: 与 md-vv 相同的 速度Verlet算法,不同之处在于动能被确定为 md 积分器中两个半步动能的平均值,因此更精确。使用 Nose-Hoover 和/或 Parrinello-Rahman 耦合,计算成本略有增加。


对于给定的步长,半步平均的动能和温度稍微准确一些;比起整步动能(使用md-vv),采用半步平均的动能(md和md-vv-avek)得到的平均动能的差异更接近于小步长极限情况下获得的动能。

对NVE模拟, 这种差别通常不明显, 因为粒子的位置与速度仍然完全相同;

差别在于模拟温度含义的解释,而不在于所产生的轨迹,虽然半步平均方法得到的动能更精确,意味着当时间步长变大时它受的影响小些,但这种方法的噪声更大。

与整步动能方法相比,半步平均动能方法得到的体系总能量(动能和势能的总和)的RMS偏差更大(在大多数情况下约高一倍)。漂移仍然是相同的,并且轨迹也完全相同。


对NVT模拟,两种方法将会有差异, 正如控温部分所讨论的,因为会调节粒子的速度以使得模拟的动能,以任何一种方式计算,都能达到与设定温度相应的分布。

在这种情况下,三种方法不会给出完全相同的结果。由于速度和位置都定义于同一时刻 , 速度Verlet积分方法可用于一些方法中, 

特别是那些严格正确的控压方法,实际上这些控压方法不能使用蛙跳式积分方法,相比蛙跳式积分方法,速度Verlet方法多耗费的时间可以忽略不计,但目前需要的通讯调用多一倍。

在大多数情况下, 尤其是对于大的体系,其通讯速度对并行化非常重要,热力学系综之间的差异以极限消失。

当只需要NVT系综时, 蛙跳式可能是首选的积分方法。

对控压模拟,热力学细节非常重要,只有速度Verlet方法能给出真正的系综。在任一情况下, 都可能需要采用双精度进行模拟以得到正确的热力学校正细节。


image.png


image.png

Verlet算法基础,但是没有依赖速度,所以没法得到速度,但是我们必须要有速度,该方法的速度误差较大,因此蛙跳法和Velocity Verlet进行修正

对于所有的方法步长不能太大,否则会因为误差导致体系的能量会炸裂


sd: 一种精确有效的蛙跳随机动力学积分器。使用约束时,每个积分步骤需要约束坐标两次。根据力计算的成本,该法可能会耗费大量模拟时间。

配合:

tc-grps:一组或多组原子的温度,单位 ref-t [K] ,

tau-t:每组的逆摩擦常数,单位tau-t [ps] 设定,忽略参数 tcoull,使用ld-seed种子初始化随机生成器。

当用作恒温器时, tau-t 最好设置为 2 ps,因为这样摩擦力较水的内摩擦低,而恰足以去除多余的热量。

注意:该法温度偏差的衰减速度是具有相同 tau-t 的 Berendsen 恒温器的两倍,也即温度回原更快。


sd2: 这曾经是默认的sd方法,但现在已被弃用。该法需要每一步每个坐标分配四个高斯随机数。施加约束后温度会稍高。

(2022.2版本已经不再支持该选项) 有约束时, 每个积分步坐标需要约束两次。取决于力的计算代价,这可能会占用显著一部分模拟时间。

精确地对随机动力学进行继续运行是不可能的,因为在模拟过程中没有存储随机数产生器的状态。


随机或速度Langevin动力学在牛顿运动方程中增加了摩擦项和噪声项

image.png

其中γi为摩擦常数[1/ps], ri为噪声过程


当与体系的时间尺度相比, 1/γi较大时, 随机动力学可视为具有随机温度耦合的分子动力学。与使用Berendsen温度耦合的MD相比, SD的优势在于它产生的系综是已知的。

模拟真空中的体系时, SD还有另外一个优点: 对整个平动转动自由度不存在误差累积.。

当 1/γi小于体系的时间尺度时, 动力学与MD完全不同, 但取样仍然正确。


sd与sd2这两种积分方法的精度等价于常规MD的蛙跳式积分方法与速度Verlet积分方法的精度, 

不同在于存在约束时复杂SD积分方法的采样温度稍微过高(尽管误差小于使用整步速度动能项的速度Verlet积分方法)。简单积分方法几乎与等同于常用的离散Langevin方程方法。

但以脉冲方式施加摩擦项与速度项,这种方案的优点在于依赖速度的项在整个时间步中都起作用, 由于对力的积分取决于坐标与速度, 这样就使得方案对力进行正确积分很简单。


bd: 作为布朗或位置郎之万动力学(position Langevin dynamics)的 Euler 积分器,

配合:

bd-fric:速度是力除以摩擦系数,单位 [amu ps-1]

ref-t:随机热噪声


当 bd fric = 0 时,每个粒子的摩擦系数计算为 mass/tau-t,与sd一样使用ld种子初始化随机生成器。

在高摩擦的极限情况下, 随机动力学退化为Brown动力学, 也被称为位置Langevin动力学. 它施加于过阻尼体系, 即惯性效应可忽略的体系. 方程为

image.png

image.png

其中riG为服从高斯分布的噪声, 其μ=0,σ=1

摩擦系数 γi可以取为对所有粒子相同, 或者取为  γi=mi γi, 其中的摩擦系数 γi对不同组的原子可取不同值。

由于假定体系处于过阻尼状态, 因此可以使用大的时间步长,应使用LINCS方法进行约束, 因为SHAKE方法对大的原子位移不收敛。


nm: Normal mode analysis。在tpr文件中的结构上执行简正模式分析,仅双精度编译GROMACS可用。

GROMACS可以进行简正模式分析, 计算是通过对角化质量加权的Hess矩阵H进行的:

image.png

其中M包含了原子质量, 

R的列为本征向量, 

λi为本征值, 

ω为相应的频率.

对通常的简正模式计算, 在计算Hess矩阵前需要进行彻底的能量最小化, 需要的能量容差取决于体系的类型, 粗略的指示值为0.001 kJ mol-1. 应当使用共轭梯度或L-BFGS方法进行能量最小化, 并使用双精度版本.

在这些计算中会涉及多个GROMACS程序. 首先应该使用mdrun进行能量最小化, 其次再使用mdrun计算Hess矩阵. 注意, 创建运行输入文件时, 应使用来自全精度轨迹文件中的最小化的构型, 因为结构文件不够精确. g_nmeig程序进行对角化计算, 并将简正模式按其频率进行排序. mdrun和g_nmeig都应当使用双精度版本. 可以使用g_anaeig程序分析简正模式. 要创建任意温度下的结构系综中的任意一组简正模式, 可使用g_nmens程序. 


tpi: Test particle insertion,插入测试粒子,主要计算过量化学势使用

拓扑结构中的最后一个分子是测试粒子。

使用该法时应当在 mdrun 处借助 -rerun 选项提供一个轨迹文件,该轨迹不应包含要插入的分子。

分子的插入将在在每帧中的随机位置以随机指向执行 nsteps 次。

当 nstlist 大于1时,使用相同的邻域列表(以及相同的长程能量(当 rvdw 或 rcoulomb > rlist 时),仅用于单原子分子)在半径为 rtpi 的球体中围绕相同的随机位置执行 nstlist 插入。

由于 neighborlist 构建的成本很高,因此这样对同一个列表执行多个额外的插入几乎不再消耗额外算力。

随机种子用 ld 种子设置。

玻尔兹曼加权的温度用 ref-t 设置,这应该与原始轨迹模拟的温度相匹配。

tpi 正确实现了色散校正。

所有相关量都会写入 mdrun -tpi 选项指定的文件。插入能量的分布被写入 mdrun -tpid 选项指定的文件。不会输出轨迹或能量文件。并行 tpi 给出了与单节点 tpi 相同的结果。

对于带电分子,使用具有精细网格的PME是最准确也是最有效的,因为系统中的电势每帧只需计算一次。


tpic: Test particle insertion into a predefined cavity location,测试粒子插入预定义空腔位置。该程序与 tpi 相同,只是从轨迹中读取一个额外的坐标,作为插入位置。

要插入的分子应该居中于 (0,0,0) 。GROMACS 不会替你完成居中,因为对于各个模拟,使用不同的居中方式可能会更好。同样 rtpi 设置该位置周围球体的半径。

每帧仅执行一次邻近搜索,不使用 nstlist。并行 tpic 给出了与单节点 tpic 相同的结果。


插入粒子,参考

#include "oplsaa.ff/forcefield.itp"

#include "oplsaa.ff/tip4pew.itp"


[ moleculetype ]

; Name          nrexcl

Methane         3


[ atoms ]

;   nr       type      resnr residue  atom   cgnr     charge       mass

     1       opls_066  1     CH4      C      1          0     16.043

 ;这个直接代表甲烷就行


[ System ]

Methane in water


[ Molecules ]

SOL               395

Methane           1 ;插入1个甲烷



您还需要将测试粒子添加到gro文件中。只需编辑conf.gro(或任何其他.gro文件使用),并在末尾添加一行包含测试粒子的位置(就在框坐标之前)。我添加的行看起来像这样:

396CH4      C 1581   0.000   0.000   0.000

只需要在不插入前的gro文件插入粒子坐标就好

实际位置并不重要;GROMACS只需要一个测试粒子的占位符。此外,您需要在.gro文件的第二行将系统中的粒子总数加1。

我们只需要一个TPI参数文件。只需从您的散装水模拟中复制prd.mdp,并将积分器更改为tpi。您应该将nsteps更改为要尝试的每帧插入次数。我选择了100000步进行模拟。您还需要将截止方案更改为组,因为Verlet尚未在TPI中实现。

$ gmx grompp -f mdp/tpi.mdp -o tpi.tpr -po tpi.mdp -pp tpi.top -c conf.gro

$ gmx mdrun -s tpi.tpr -o tpi.trr -x tpi.xtc -c tpi.gro -e tpi.edr -g tpi.log -rerun prd.xtc


在本例中,日志文件名为tpi.log,其中包含一行平均体积和平均过量化学势。我的两行看起来像这样:

<V> =1.18704e+01  nm^3          <mu>=8.81230e+00  kJ/mol

<mu>的输出单位为kJ/mol,但如果我们将其转换为kcal/mol,我们得到2.106 kcal/mol。


在热力学和物理化学中,超额化学势(excess chemical potential)是一个用来描述混合物中某一组分的化学势与其在纯物质状态下的化学势之间差异的概念。

它是超额函数(excess function)的一个方面,用于表征混合物或溶液的非理想行为。


具体来说,对于混合物中的组分i,超额化学势定义为该组分在混合物中的化学势与其在理想溶液中的化学势之差。数学上可以表示为:


\[ \mu_i^E = \mu_i - \mu_i^{"} \]

其中:

- \(\mu_i^E\) 是组分i在混合物中的超额化学势。

- \(\mu_i\) 是组分i在混合物中的化学势。

- \(\mu_i^{"}\) 是组分i在纯物质状态或理想溶液状态下的化学势。


超额化学势的概念有助于我们理解混合物中组分之间相互作用所导致的非理想行为,如非理想混合、溶液中的活性系数、相分离等现象。

在实际应用中,超额化学势可以通过实验测定,例如通过测量混合物的蒸汽压、凝固点下降、沸点升高或者渗透压等热力学性质来得到。

此外,它也是许多热力学模型和活度系数模型的计算基础,如享利定律、范特霍夫定律、Raoult定律以及更为复杂的电解质溶液理论。

在化学工程和工业应用中,了解超额化学势对于过程设计和优化至关重要,因为它影响了混合物组分在分离过程中的行为,比如蒸馏、吸收和萃取等单元操作。


mimic: 启用 MiMiC QM/MM 以运行混合分子动力学模拟。记住使用该选项还需要启动用 MiMiC 编译的CPMD。在这种模式下,所有关于积分的选项(T耦合、P耦合、步长和步数等)都被忽略,因为相应积分将CPMD完成。与力计算相关的选项(截止值、PME参数等)正常工作。QMMM-grps 读取对 QM 原子的选择(该选项为2022.2版本新增,不适用于部分老旧版本)。



最新评论
请先登录才能进行回复登录
技术支持: CLOUD | 管理登录
seo seo