在GROMACS中使用OPC水模型

GROMACS中OPC水的拓扑文件

image-20250210152116532

与常用的三位点水模型相比，OPC多了一个虚拟电荷位点，总共有4个电荷位点和在氧原子上的单个Lennard-Jones位点参与相互作用。

首先确认OPC水的势函数类型与需要配合使用的力场是否兼容，（力场）拓扑文件的开头一般会有[ defaults ]字段，定义非键相互作用的计算规则：

例如：

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

GROMACS中非键相互作用为对势加和，范德华作用可由Lennard-Jones或Buckingham势计算，其中Lennard-Jones写成： \[ V_{LJ}(r_{ij})=\frac{C_{ij}^{(12)}}{r_{ij}^{12}}-\frac{C_{ij}^{(6)}}{r_{ij}^{6}} \tag{1} \] 或 \[ V_{LJ}(r_{ij})=4\epsilon_{ij}\left(\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^6\right) \tag{2} \]

[ defaults ] :

  • nbfunc Use 1 (Lennard-Jones) or 2 (Buckingham).
  • comb-rule
  • gen-pairs is for pair generation. The default is ‘no’, i.e. get 1-4 parameters from the pairtypes list. When parameters are not present in the list, stop with a fatal error. Setting ‘yes’ generates 1-4 parameters that are not present in the pair list from normal Lennard-Jones parameters using fudgeLJ
  • fudgeLJ is the factor by which to multiply Lennard-Jones 1-4 interactions, default 1
  • fudgeQQ is the factor by which to multiply electrostatic 1-4 interactions, default 1
  • NN is the power for the repulsion term in a 6-NN potential (with nonbonded-type Lennard-Jones only), starting with GROMACS version 4.5, grompp also reads and applies NN, for values not equal to 12 tabulated interaction functions are used (in older version you would have to use user tabulated interactions).

Note that gen-pairs, fudgeLJ, fudgeQQ, and NN are optional. fudgeLJ is only used when generate pairs is set to ‘yes’, and fudgeQQ is always used. However, if you want to specify NN you need to give a value for the other parameters as well.

首先需要将OPC水的原子类型加入需要配合使用的力场文件中：

1
2
3
;To use this opc.itp, you must also add below lines into [ atomtypes ] of ffnonbonded.itp under AMBER forcefield
;OW_opc      8      15.9994  0.0000  A   0.316655  0.89036
;HW_opc      1       1.0080  0.0000  A   0.0       0.0

如果是使用tleap建模的流程，则将原子类型加入complex_GMX.top的[ atomtypes ]中

这里给出范德华参数C6和C12：

1
2
3
OW_opc   OW_opc      15.9994  0.0000    A     0.316655      0.89036
HW_opc   HW_opc       1.0080  0.0000    A     0.0           0.0
MW       MW           0.0     0.0000    A     0.0           0.0

SOB给出的opc.itp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
[ moleculetype ]
; molname	nrexcl
SOL		2

[ atoms ]
; id  at type     res nr  res name  at name  cg nr  charge    mass
  1   OW_opc  1       SOL       OW       1       0        16.00000
  2   HW_opc  1       SOL       HW1      1       0.6791   1.00800
  3   HW_opc  1       SOL       HW2      1       0.6791   1.00800
  4   MW      1       SOL       MW       1      -1.3582   0.00000

#ifndef FLEXIBLE

[ settles ]
; i	funct	doh	dhh
1	1	0.08724	0.137116

#else
;OPC paper didn't define force constant for flexible mode, so the FC values are set to TIP4P-Ew ones
[ bonds ]
; i     j       funct   length  force.c.
1       2       1       0.08724 502416.0 0.08724        502416.0 
1       3       1       0.08724 502416.0 0.08724        502416.0 
        
[ angles ]
; i     j       k       funct   angle   force.c.
2       1       3       1       103.6  628.02  103.6  628.02  

#endif


[ virtual_sites3 ]
; Vsite from                    funct   a               b
4       1       2       3       1       0.14772952     0.14772952
;
; a = b = distance(O-MW) / [ cos (angle(MW-O-H)) * distance(OH) ] /2
;	      = 0.01594/(cos(103.6/2/180*3.141592653589793)*0.08724)/2

[ exclusions ]
1	2	3	4
2	1	3	4
3	1	2	4
4	1	2	3

OPC离子参数：

image-20250210190045066

使用LB combined rule：

image-20250210190122263

参考https://github.com/intbio/gromacs_ff/blob/master/amber99bsc1.ff/ions.itp写opc_ionsitp文件，重复：

1
2
3
4
5
6
7
[ moleculetype ]
; molname       nrexcl
IB+             1       ; big positive ion

[ atoms ]
; id    at type         res nr  residu name     at name  cg nr  charge
   1    IB                  1       IB+            IB      1    1.00000

范德华参数在力场的ffnonbonded.itp中或top文件加入：

参考此处的建议，采用Kenneth M. Merz, Jr等人拟合的单价离子的IOD参数集：

image-20250210193856226

比较Amber和GROMACS中C6和C12 \[ \begin{align} C_{12}^{ij}&=\epsilon^a_{ij}R^{12}_{min,ij}\\ &=4\epsilon^g_{ij}\sigma^{12}_{ij} \end{align} \]

\[ \begin{align} C_{6}^{ij}&=2\epsilon^a_{ij}R^{6}_{min,ij}\\ &=4\epsilon^g_{ij}\sigma^{6}_{ij} \end{align} \]

因而 \[ \epsilon^g_{ij}(kJ/mol)=\epsilon^a_{ij}(kcal/mol)*4.182 \]

\[ \sigma_{ij}=R_{min,ij}(Angstrom)/2^{\frac{1}{6} }/10 \]

Rmin/2 (Å) ε (kcal/mol) σ (nm) ε (kJ/mol)
Li+ 1.321 0.00642 0.11769 0.02683
Na+ 1.47 0.03038 0.13096 0.12706
K+ 1.743 0.16869 0.15528 0.70548
Rb+ 1.81 0.22132 0.16125 0.92558
Cs+ 1.99 0.38035 0.17729 1.59063
Tl+ 1.866 0.26895 0.16624 1.12474
Cu+ 1.213 0.00137 0.10807 0.00573
Ag+ 1.502 0.03963 0.13381 0.16572
F– 1.737 0.16427 0.15475 0.68697
Cl– 2.16 0.52989 0.19243 2.21598
Br– 2.315 0.64855 0.20624 2.71224
I– 2.59 0.80294 0.23074 3.35789

在complex_GMX.top中加入

1
2
3
4
5
6
7
8
9
10
11
12
Li       Li           6.94    0.00000   A     0.11769       0.02683 
Na       Na          22.99    0.00000   A     0.13096       0.12706 
K        K           39.10    0.00000   A     0.15528       0.70548 
Rb       Rb          85.47    0.00000   A     0.16125       0.92558 
Cs       Cs         132.91    0.00000   A     0.17729       1.59063 
Tl       Tl         204.38    0.00000   A     0.16624       1.12474 
Cu       Cu          63.546   0.00000   A     0.10807       0.00573 
Ag       Ag         107.87    0.00000   A     0.13381       0.16572 
F        F           19.00    0.00000   A     0.15475       0.68697 
Cl       Cl          35.45    0.00000   A     0.19243       2.21598 
Br       Br          79.90    0.00000   A     0.20624       2.71224 
I        I          126.9     0.00000   A     0.23074       3.35789

GROMACS 蛋白质水溶液建模

我们已通过AmberTools23中的tleap工具和acpype获得蛋白质-配体复合物的拓扑文件，tleap输入文件如下

1
2
3
4
5
6
7
8
9
10
11
source leaprc.protein.ff19SB
source leaprc.gaff2
loadamberparams ../subA_resp/subA.frcmod
loadoff ../subA_resp/NTA.lib
loadamberparams ../subB_resp/subB.frcmod
loadoff ../subB_resp/NTB.lib
complex = loadpdb complex_input.pdb #注意：要使得pdb中配体原子名称与mol2中匹配，蛋白质端基要TER
check complex
savepdb complex complex_leap.pdb
saveamberparm complex complex.prmtop complex.inpcrd
quit

我们得到了complex_GMX.top、complex_GMX.gro作为进一步建模的起点。

要使用水模型，我们需要在complex_GMX.top文件中加入以下字段。首先include水模型的拓扑文件

1
2
3
4
5
6
7
8
9
10
; Include water topology
#include "opc.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif
; Include topology for ions
#include "opc_ions.itp"

首先创建盒子和加水

1
2
gmx editconf -f complex_GMX.gro -o newbox.gro -bt cubic -d 1.2 -c
gmx solvate -cp newbox.gro -cs tip4p.gro -p complex_GMX.top -o solv.gro

editconf创建立方盒子(-bt cubic)，分子边缘距离盒子边缘至少1.2 nm(-d 1.2)，将分子置于中心(-c)

solvate用tip4p.gro中的溶剂分子填充盒子，它会在complex_GMX.top中写入溶剂分子数目：

1
2
3
4
5
6
7
[ system ]
complex in water

[ molecules ]
; Compound        nmols
 complex          1     
SOL             25863

然后添加离子

1
gmx grompp -f ions.mdp -c solv.gro -p complex_GMX.top -o ions.tpr

这里ions.mdp是一个基本的参数文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
; LINES STARTING WITH ';' ARE COMMENTS
title		    = Minimization	; Title of run

; Parameters describing what to do, when to stop and what to save
integrator	    = steep		; Algorithm (steep = steepest descent minimization)
emtol		    = 1000.0  	; Stop minimization when the maximum force < 10.0 kJ/mol
emstep          = 0.01      ; Energy step size
nsteps		    = 50000	  	; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist		    = 1		    ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet
ns_type		    = grid		; Method to determine neighbor list (simple, grid)
rlist		    = 1.0		; Cut-off for making neighbor list (short range forces)
coulombtype	    = cutoff	; Treatment of long range electrostatic interactions
rcoulomb	    = 1.0		; long range electrostatic cut-off
rvdw		    = 1.0		; long range Van der Waals cut-off
pbc             = xyz 		; Periodic Boundary Conditions

之后用Monte Carlo方法在合适的位置添加离子

1
gmx genion -s ions.tpr -o complex_solv_GMX.gro -p complex_GMX.top -pname NA -nname CL -neutral

加入了6个Na+：

1
2
3
4
5
[ molecules ]
; Compound        nmols
 complex          1
SOL         25857
NA               6

习惯用拓扑文件名为 complex_solv_GMX.top

在拓扑文件中复合物拓扑之后加入，以便使用对复合物重原子的位置限制

1
2
3
4
; Include Position restraint file
#ifdef POSRES
#include "posre_complex.itp"
#endif

模拟工作流

运行前，这些文件需复制到工作文件夹

1
2
3
4
5
6
opc_ions.itp
opc.itp
em.mdp
complex_solv_GMX.gro
posre_complex.itp
complex_solv_GMX.top

能量最小化

1
2
gmx grompp -f em.mdp -c complex_solv_GMX.gro -p complex_solv_GMX.top -o em.tpr
gmx mdrun -deffnm em

其中参数文件为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
; minim.mdp - used as input into grompp to generate em.tpr
title       = Minimization
integrator  = steep     ; Algorithm (steep = steepest descent minimization)
emtol       = 500.0    ; Stop minimization when the maximum force < 50.0 kJ/mol/nm
emstep      = 0.01      ; Energy step size
nsteps      = 50000     ; Maximum number of (minimization) steps to perform
comm-mode   = Linear    ; remove translation of center of mass
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
cutoff-scheme   = Verlet
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = PME       ; Treatment of long range electrostatic interaction
rcoulomb        = 1.2       ; Short-range electrostatic cut-off
rvdw            = 1.2       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions (yes/no)

之后便是温度和盒子的预平衡，使得系统充分驰豫，以便在NPT系综下进行生产相模拟。该部分涉及热浴和压浴。

如果不担心崩溃，可以直接在NPT系综中驰豫：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
title       = NPT equilibration 
define                   = -DPOSRES ; restraint the complex 
; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 500000     ; 2 * 500000 =  1 ns
dt          = 0.002     ; 2 fs
; Output control
nstxout     = 0         ; save coordinates every 1.0 ps
nstvout     = 0         ; save velocities every 1.0 ps
nstcalcenergy    = 500  ; calculate energies every 1.0 ps
nstenergy   = 500       ; save energies every 1.0 ps
nstlog      = 500       ; update log file every 1.0 ps
energygrps  = System
nstxout-compressed  = 500      ; write .xtc trajectory every 1.0 ps
compressed-x-grps   = System 
; Bond parameters
continuation    = no           ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints 
constraints     = h-bonds     ; all bonds (even heavy atom-H bonds) constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
cutoff-scheme   = Verlet
ns_type         = grid      ; search neighboring grid cells
nstlist         = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb        = 1.0      ; short-range electrostatic cutoff (in nm)
rvdw            = 1.0       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype     = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order       = 4         ; cubic interpolation
fourierspacing  = 0.16      ; grid spacing for FFT
; vdW 
;vdwtype          = PME
;vdw-modifier     = Potential-Shift
;lj-pme-comb-rule = Lorentz-Berthelot
; Temperature coupling
tcoupl      = V-rescale                     ; modified Berendsen thermostat
tc-grps     = System                         ; coupling group
tau_t       = 0.1                           ; time constant, in ps
ref_t       = 300                           ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl      = C-rescale             ; pressure coupling is on for NPT
pcoupltype  = isotropic                     ; uniform scaling of box vectors
tau_p       = 0.3                           ; time constant, in ps
ref_p       = 1.0                           ; reference pressure, in bar
compressibility = 4.5e-5                    ; isothermal compressibility of water, bar^-1
refcoord_scaling    = com
; Periodic boundary conditions
pbc         = xyz       ; 3-D PBC
; Dispersion correction
;DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = yes        ; velocity generation off after NPT

一般先限制复合物驰豫，然后放开复合物平衡。

含虚拟位点的体系在gpu计算上有一些限制，不能在gpu上update，这里实使用gmx2024：

1
2
3
4
5
6
gmx grompp -f ../npt2.mdp -c em.gro -p complex_solv_GMX.top  -o npt2.tpr -r em.gro 
gmx mdrun -deffnm npt2 -ntmpi $NGPUS -ntomp $OMPTHREADS -gpu_id 0 -nb gpu  -pme gpu -bonded gpu 
gmx grompp -c npt2.gro -t npt2.cpt  -p complex_solv_GMX.top  -f ../run.mdp -o run.tpr
gmx mdrun -deffnm run -ntmpi $NGPUS -ntomp $OMPTHREADS -gpu_id 0 -nb gpu  -pme gpu -bonded gpu 
gmx grompp -c run.gro -t run.cpt  -p complex_solv_GMX.top  -f ../run2.mdp -o run2.tpr
gmx mdrun -deffnm run2 -ntmpi $NGPUS -ntomp $OMPTHREADS -gpu_id 0 -nb gpu  -pme gpu -bonded gpu

Reference

  1. Izadi, S.; Anandakrishnan, R.; Onufriev, A. V. Building Water Models: A Different Approach. J. Phys. Chem. Lett. 2014, 5 (21), 3863–3871. https://doi.org/10.1021/jz501780a.
  2. Sengupta, A.; Li, Z.; Song, L. F.; Li, P.; Merz, K. M. Parameterization of Monovalent Ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB Water Models. J. Chem. Inf. Model. 2021, 61 (2), 869–880. https://doi.org/10.1021/acs.jcim.0c01390.
Science > Theoretical and Computational Chemistry > Molecular Modeling
#GROMACS #Molecular Modeling #Classical Force Field
在GROMACS中使用OPC水模型
http://example.com/2025/02/10/2025-02-10-在GROMACS中使用OPC水模型/
Author
Graminilune
Posted on
February 10, 2025
Licensed under