GROMACS version: 2022.4
GROMACS modification: Yes (Colvars module)
Hi all,
I have been trying to model an sII THF hydrate in GROMACS 2022.4 for about a month now. My starting structure is a hydrate determined from x-ray crystallography, with THF molecules added manually to large cages found in the hydrate.
When I perform an energy minimization on my starting structure, I consistently find that my THF molecules move to one side of the box (I’ve included the gro file). However, the actual water hydrate structure is mostly maintained. This should not be the case - this system is found experimentally and should be stable.
Here are the things that I have checked/tried:
- I’ve run tried using different water models. These include TIP3P, TIP4P, and TIP4P/ICE. When I simulated the hydrate alone (no THF), it was only stable with the TIP4P models. However, it did not fix my issue.
- I checked my THF forcefield. This is what we believed to be the issue, since it was made by LigParGen, and then updated by me with accurate charges and epsilon/sigma values from literature. However, after simulating a THF-only system we achieved densities close to experimental values, which supports our force field being relatively accurate.
- I recently also tried restraining the THF in their positions in the center of the cages. I believe I implemented these restraints correctly, and with the default value of 1000 I was unable to keep the THF within their cages.
My PI and I have the suspicion that there might be something wrong with the system setup, or the mdp for energy minimization. The only way I can get energy minimization to run is to not use any constraints (in my case h-bond LINCS). When I do use LINCS, my simulation is stopped by warnings, and when I examine the pdb I find the THF twisted. So something is happening to the THF which is unrealistic, but I’m very stuck on how to solve it.
So my main question is, am I doing something wrong with running my simulations like in my mdp? Or is there simply something wrong with my system/forcefield? Any help is appreciated, like I said this has been a headache for a month and I really need to get it figured out.
I’ve attached image of the ‘successful’ em w/o LINCS. I also have pasted the topol.top, em.mdp, and thf.itp files below. I can send a zip of all supporting files to people individually. Thank you for your help, and let me know what other information you need!
Image
thf.itp
[ atomtypes ]
opls_800 C00 12.0110 0.000 A 3.50000E-01 2.76144E-01
opls_801 C01 12.0110 0.000 A 3.50000E-01 2.76144E-01
opls_802 C02 12.0110 0.000 A 3.50000E-01 2.76144E-01
opls_803 C03 12.0110 0.000 A 3.50000E-01 2.76144E-01
opls_804 O04 15.9990 0.000 A 2.90000E-01 5.50000E-01
opls_805 H05 1.0080 0.000 A 2.45000E-01 1.25520E-01
opls_806 H06 1.0080 0.000 A 2.45000E-01 1.25520E-01
opls_807 H07 1.0080 0.000 A 2.45000E-01 1.25520E-01
opls_808 H08 1.0080 0.000 A 2.45000E-01 1.25520E-01
opls_809 H09 1.0080 0.000 A 2.45000E-01 1.25520E-01
opls_810 H0A 1.0080 0.000 A 2.45000E-01 1.25520E-01
opls_811 H0B 1.0080 0.000 A 2.45000E-01 1.25520E-01
opls_812 H0C 1.0080 0.000 A 2.45000E-01 1.25520E-01
; Non-bonded LJ Potentials from ‘Molecular Dynamics Study on the Growth Mechanism of Methane plus Tetrahydrofuran Mixed Hydrates’ Wu et. al.
;[ nonbond_params ]
; i j func V(c6)(nm) W(c12)
;opls_800 opls_113 1 .3430 0.2547
;opls_801 opls_113 1 .3430 0.2547
;opls_802 opls_113 1 .3430 0.2547
;opls_803 opls_113 1 .3430 0.2547
;opls_804 opls_113 1 .3163 0.1509
; Map of atoms in THF
;
; O04
; H05 / \ H0C
; \ / \ /
; C00 C03
; / \ /
; H06 \ / H0B
; H07 - C01 ----- C02 - H0A
; | |
; H08 H09
;
[ moleculetype ]
; Name nrexcl
THF 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 opls_800 1 THF C00 1 0.2840 12.0110
2 opls_801 1 THF C01 1 0.0000 12.0110
3 opls_802 1 THF C02 1 0.0000 12.0110
4 opls_803 1 THF C03 1 0.2840 12.0110
5 opls_804 1 THF O04 1 -0.5680 15.9990
6 opls_805 1 THF H05 1 0.0000 1.0080
7 opls_806 1 THF H06 1 0.0000 1.0080
8 opls_807 1 THF H07 1 0.0000 1.0080
9 opls_808 1 THF H08 1 0.0000 1.0080
10 opls_809 1 THF H09 1 0.0000 1.0080
11 opls_810 1 THF H0A 1 0.0000 1.0080
12 opls_811 1 THF H0B 1 0.0000 1.0080
13 opls_812 1 THF H0C 1 0.0000 1.0080
[ bonds ]
2 1 1 0.1529 224262.400
3 2 1 0.1529 224262.400
4 3 1 0.1529 224262.400
5 1 1 0.1410 267776.000
6 1 1 0.1090 284512.000
7 1 1 0.1090 284512.000
8 2 1 0.1090 284512.000
9 2 1 0.1090 284512.000
10 3 1 0.1090 284512.000
11 3 1 0.1090 284512.000
12 4 1 0.1090 284512.000
13 4 1 0.1090 284512.000
5 4 1 0.1410 267776.000
[ angles ]
; ai aj ak funct c0 c1 c2 c3
1 2 3 1 112.700 488.273
2 3 4 1 112.700 488.273
2 1 5 1 109.500 418.400
2 1 6 1 110.700 313.800
2 1 7 1 110.700 313.800
1 2 8 1 110.700 313.800
1 2 9 1 110.700 313.800
2 3 10 1 110.700 313.800
2 3 11 1 110.700 313.800
3 4 12 1 110.700 313.800
3 4 13 1 110.700 313.800
4 3 11 1 110.700 313.800
5 4 12 1 109.500 292.880
10 3 11 1 107.800 276.144
5 1 7 1 109.500 292.880
1 5 4 1 109.500 502.080
5 4 13 1 109.500 292.880
4 3 10 1 110.700 313.800
3 2 8 1 110.700 313.800
8 2 9 1 107.800 276.144
5 1 6 1 109.500 292.880
12 4 13 1 107.800 276.144
3 2 9 1 110.700 313.800
6 1 7 1 107.800 276.144
3 4 5 1 109.500 418.400
[ dihedrals ]
; IMPROPER DIHEDRAL ANGLES
; ai aj ak al funct c0 c1 c2 c3 c4 c5
[ dihedrals ]
; PROPER DIHEDRAL ANGLES
; ai aj ak al funct c0 c1 c2 c3 c4 c5
4 3 2 1 3 2.301 -1.464 0.837 -1.674 -0.000 0.000
3 4 5 1 3 1.715 2.845 1.046 -5.607 -0.000 0.000
4 5 1 2 3 1.715 2.845 1.046 -5.607 -0.000 0.000
9 2 3 4 3 0.628 1.883 0.000 -2.510 -0.000 0.000
11 3 2 1 3 0.628 1.883 0.000 -2.510 -0.000 0.000
8 2 3 4 3 0.628 1.883 0.000 -2.510 -0.000 0.000
7 1 2 3 3 0.628 1.883 0.000 -2.510 -0.000 0.000
6 1 2 3 3 0.628 1.883 0.000 -2.510 -0.000 0.000
10 3 2 1 3 0.628 1.883 0.000 -2.510 -0.000 0.000
13 4 3 2 3 0.628 1.883 0.000 -2.510 -0.000 0.000
12 4 3 2 3 0.628 1.883 0.000 -2.510 -0.000 0.000
8 2 1 7 3 0.628 1.883 0.000 -2.510 -0.000 0.000
9 2 1 6 3 0.628 1.883 0.000 -2.510 -0.000 0.000
9 2 1 7 3 0.628 1.883 0.000 -2.510 -0.000 0.000
12 4 3 11 3 0.628 1.883 0.000 -2.510 -0.000 0.000
13 4 3 11 3 0.628 1.883 0.000 -2.510 -0.000 0.000
10 3 2 8 3 0.628 1.883 0.000 -2.510 -0.000 0.000
13 4 3 10 3 0.628 1.883 0.000 -2.510 -0.000 0.000
10 3 2 9 3 0.628 1.883 0.000 -2.510 -0.000 0.000
8 2 1 6 3 0.628 1.883 0.000 -2.510 -0.000 0.000
12 4 3 10 3 0.628 1.883 0.000 -2.510 -0.000 0.000
11 3 2 8 3 0.628 1.883 0.000 -2.510 -0.000 0.000
11 3 2 9 3 0.628 1.883 0.000 -2.510 -0.000 0.000
10 3 4 5 3 0.979 2.937 0.000 -3.916 -0.000 0.000
11 3 4 5 3 0.979 2.937 0.000 -3.916 -0.000 0.000
9 2 1 5 3 0.979 2.937 0.000 -3.916 -0.000 0.000
8 2 1 5 3 0.979 2.937 0.000 -3.916 -0.000 0.000
7 1 5 4 3 1.590 4.770 0.000 -6.360 -0.000 0.000
6 1 5 4 3 1.590 4.770 0.000 -6.360 -0.000 0.000
13 4 5 1 3 1.590 4.770 0.000 -6.360 -0.000 0.000
12 4 5 1 3 1.590 4.770 0.000 -6.360 -0.000 0.000
5 4 3 2 3 2.874 0.582 2.092 -5.548 -0.000 0.000
5 1 2 3 3 2.874 0.582 2.092 -5.548 -0.000 0.000
[ pairs ]
3 6 1
4 6 1
3 7 1
4 7 1
1 10 1
4 8 1
1 11 1
5 8 1
4 9 1
1 12 1
6 8 1
5 9 1
2 12 1
1 13 1
7 8 1
6 9 1
5 10 1
2 13 1
7 9 1
5 11 1
8 10 1
9 10 1
8 11 1
9 11 1
10 12 1
11 12 1
10 13 1
11 13 1
topol.top
; Topology for THF sII hydrate
#include “oplsaa.ff/forcefield.itp”
; THF topology
#include “oplsaa.ff/thf.itp”
#ifdef POSRES
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
2 1 1000 1000 1000
3 1 1000 1000 1000
4 1 1000 1000 1000
5 1 1000 1000 1000
6 1 1000 1000 1000
7 1 1000 1000 1000
8 1 1000 1000 1000
9 1 1000 1000 1000
10 1 1000 1000 1000
11 1 1000 1000 1000
12 1 1000 1000 1000
13 1 1000 1000 1000
#endif
; Water topology
#include “oplsaa.ff/tip4p.itp”
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
[ system ]
; Name
THF sII hydrate 3x3x3
[ molecules ]
; Compound #mols
THF 216
SOL 3672
em.mdp
; MDP for THF hydrate simulation
define = -DPOSRES
; Run control
integrator = steep ; The algorithm being used
emtol = 1000.0 ; Stops minimization at Fmax < 1000 kj/mol/nm
emstep = 0.01 ; Just step size
nsteps = 50000 ; Maximum number of steps
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor search
ns_type = grid ; Method to determine neighbor list
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
rlist = 1.0 ; Cut-off distance for short-range neighbor list
; Electrostatics
coulombtype = PME ; Model used for electrostaic interactions
rcoulomb = 2.5 ; Distance cut-off from Coulomb (electrostatic) interactions
; van der Waals
vdwtype = cutoff ; As I understand it, this is simply saving computational work. PME would be more accurate because it is more dynamic (?) - I’d have to read more.
rvdw = 2.5
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid (Need to do research)
fourierspacing = 0.12
; EWALD/PME/PPPM parameters (Need to do research)
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
; Temperature and pressure coupling are off during EM
tcoupl = no
pcoupl = no
; No velocities during EM
gen_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4
;lincs-iter = 2