CO2超临界模拟

超临界CO₂性质模拟

模拟设置

mdp设置

控温技术：V-rescale

Berendsen热浴会导致体系的粒子速度分布偏离Maxwell分布，导致粒子间的速度差异被拉大
V-rescale 解决了粒子速度不满足Maxwell分布的问题
控压技术：Parrinello-Rahman

Berendsen压浴无法产生精确的NPT系综，对大体系不明显，但是Parrinello-Rahman可以产生正确的NPT系综，其缺点是容易产生类似Nose-Hoover热浴的震荡行为，因此不适合温度变化的过程

topol构建

从该网站TraPPE: Transferable Potentials for Phase Equilibria Force Field (umn.edu)http://trappe.oit.umn.edu/获取CO₂的LJ参数

carbon dioxide
#	(pseudo)atom	type	epsilon/kB [K]	sigma [Ang.]	charge [e]	epsilon/kJ/mol
1	C	O=[C]=O	27	2.7918	0.7	2.24E-01
2	O	[O]=C=O	79	3	-0.35	6.56E-01
3	O	[O]=C=O	79	3	-0.35	6.56E-01
#	stretch	type	length [Ang.]
1	‘1 - 2’	O=(C=O)	1.16		K_b(J/K)	NA
2	‘2 - 3’	O=(C=O)	1.16		1.38E-23	6.02E+23
#	bend	type	theta [degrees]	k_theta/kB [K/rad^2]
1	‘2 - 1 - 3’	O=©=O	180	0

topol.top

TraPPE力场的CO₂的LJ参数需要转换一下：

其中： epsilon($\epsilon$)=值*K_b*NA*0.001

topol中sigma($\sigma$)的单位为nm

对于键长使用constraints进行限制

其他部分见下节

; CO2_GMX.itp created by acpype (v: 2022.6.6) on Fri Dec  2 12:12:40 2022
[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               yes             0 0

[ atomtypes ]
;name   bond_type     mass     charge   ptype   sigma         epsilon       
 o        o           0.00000  0.00000   A     3.05e-01   6.56e-01 ; 
 c1       c1          0.00000  0.00000   A     2.8e-01   2.24e-01 ; 

[ moleculetype ]
;name            nrexcl
 COO              3

[ atoms ]
;   nr  type  resi  res  atom  cgnr     charge      mass       ; 
     1    o     1   COO     O    1    -0.35     16.00000 ; 
     2    o     1   COO    O1    2    -0.35     16.00000 ; 
     3   c1     1   COO     C    3      0.7     12.01000 ;

; [ bonds ]
; ;   ai     aj funct   r             k
;      1      3   1    1.16e-01       6.6526e+05 ;      O - C     
;      2      3   1    1.16e-01       6.6526e+05 ;     O1 - C     

[ constraints ]
    1      3   1    1.16e-01      
    2      3   1    1.16e-01     

[ angles ]
;   ai     aj     ak    funct   theta         cth
     1      3      2      1    1.8000e+02    8.7864e+02 ;      O - C      - O1    

[ system ]
 CO2

[ molecules ]
; Compound        nmols
 COO              323

mol.gro

CO2
3
    1COO      O    1   3.020   2.554   3.060  0.3762 -0.1350  0.1096
    1COO     O1    2   3.051   2.380   2.909  0.9207 -0.1455  0.2320
    1COO      C    3   3.037   2.467   2.985  0.5628 -0.4036  0.4654

体系使用packmol构建，5*5*5 nm 盒子， 323 个分子

# packmol.inp
tolerance 2.0
output structure.pdb
structure CO2.pdb
number 323
inside cube 0. 0. 0. 50.
end structure

1	packmol < packmol.inp

模拟

检验模拟结果与实验的符合程度。

将模拟数据与实验数据相比，压力较低时，密度偏小。压力较大时，密度偏大。

拓扑文件组成


相互作用类型	指令	原子数	函数类型	参数
必需	`defaults`		非键函数类型; 组合规则cr; 生成对(no/yes); fudge LJ(); fudge QQ()
必需	`atomtypes`		原子类型; 质量m(u); 电荷q(e); 粒子类型; Vcr; Wcr
	`bondtypes`			`bonds`
	`pairtypes`			`pairs`
	`angletypes`			`angles`
	`dihedraltypes`*			`dihedrals`
	`constrainttypes`			`constraints`
LJ	`nonbond_params`	2	1	V^cr; W^cr
Buckingham	`nonbond_params`	2	2	a(kJ mol^-1); b(nm^-1); c6(kJ mol^-1 nm⁶)
必需	`moleculetype`			分子名称; n_ex^nrexcl
必需	`atoms`	1		原子类型; 残基编号; 残基名称; 原子名称; 电荷组编号; 电荷q(e); 质量m(u)
必需	`system`			系统名称
必需	`molecules`			分子名称; 分子数

原子数：对于每个指令需要制定的原子类型
函数类型：选择的函数类型
nrexcl: 计算非键相互作用时，需要排除相邻的n个键

原子类型性质

Proterty	Symbol	Unit
Type
Mass	m	a.m.u
Charge	q	electron
epsilon	$\epsilon$	kJ/mol
sigma	$\sigma$	nm

非键参数

非键参数由范德华参数V（c6 或$\sigma$，取决于结合规则）和W（c12 或 $\epsilon$）组成，ptype为粒子类型，见下表

Particle	Symbol
atom	A
shell	S
virtual site	V（或 D）

在topol[ defaults ] block中定义了结合规则

结合规则1：

$V_{ii}=C_{ij}^{(6)} =4\epsilon _{i} \sigma _{i}^{6}\ [kJ\ mol^{-1}\ nm^{6}]$

$W_{ii}=C_{ij}^{(12)} =4\epsilon _{i} \sigma _{i}^{12}\ [kJ\ mol^{-1}\ nm^{12}]$

结合规则2和3：

$V_{ii}=\sigma_{i}\ [\ nm\ ]$

$W_{ii}=\epsilon_{i}\ [\ kJ\ mol^{-1}\ ]$

[ defaults ]

1
2
3

[ defaults ]
; nbfunc    comb-rule    gen-pairs    fudgeLJ    fudgeQQ
1           2            yes          0.5        0.8333

[ defaults ]:

nbfunc为非键函数类型, 取1(Lennard-Jones)或2(Buckingham)
comb-rule为结合规则的编号.
gen-pairs用于对的生成. 默认为no, 即得到对类型列表的1–4的参数. 当列表中不存在参数时, 会导致致命错误并停止运行. 设置为yes会根据使用fudgeLJ的正常Lennard-Jones参数, 生成不存在于配对列表中的1–4参数
fudgeLJ为Lennard-Jones 1–4相互作用的修正因子, 默认为1
fudgeQQ为静电1–4相互作用的修正因子, 默认为1
N为6-N势中排斥项的次数(只用于非键Lennard-Jones), 自GROMACS 4.5版本开始, mdrun会读取和使用N, 当N不等于12时, 使用表格相互作用函数(在旧版本中, 你必须使用用户自定义的表格相互作用).

总结

CO2的top参数需要自己查询文献，同时需要将获得的参数改为gromacs格式
topol的组成还是很值得玩味的，我这只粗略地列出了一部分，详细可见手册以及参考[1]。
第一次做非生物体系，可能有很多设置错误，此文仅为一个记录，仅供参考。

引用

[1] GROMACS中文手册：第五章　拓扑文件|Jerkwin,http://jerkwin.github.io/GMX/GMXman-5/#55-约束算法

[2] CO2的EPM2模型参数如何加到itp文件中？ - 分子模拟 (Molecular Modeling) - 计算化学公社 (keinsci.com)(http://bbs.keinsci.com/forum.php?mod=viewthread&tid=18748&highlight=EPM)