[gmx-users] CG Lincs errors

[Nash, Anthony]

Hi all,

I¹m hoping for some help. I¹m very sorry, this is a bit of a long one.

I¹ve been struggling for almost a month trying to run a CG representation
of our all-atom model of a collagen protein (3 polypeptide chains in a
protein). Our original AMBER all-atom model had been successful modelling
using MD. I went on to use the latest version of Martinize.py with the
latest version of the MARTINI forcefield fields.

After a little tweaking (the way AMBER names histidine residues), I
successful converted the molecule (approx 3100 amino acids) into a CG
representation. I successfully energy min the protein in vacuum to a
threshold of 500, and in solvent to a threshold of 750 using steepest
descent. Looking for a system at an energy min of a threshold around 300 I
begin to see LINCS warnings. Observing the initial structure, there is
nothing obviously wrong with the bond network (both protein and polarized
CG water).

I take the system that energy mins at 750 (protein-water mix, with no
fault reported), and went straight to NPT, 20fs step. Blew up. After a bit
of chatting with the MARTINI community, I¹ve started with an NVT ensemble,
beginning at 5s then through 10fs, 15fs, and 20fs. I only run for 1000
steps before switching. Keeping any of the simulations running for longer
throws lincs warnings followed by a segmentation fault from the warning:

"3 particles communicated to PME rank 7 are more than 2/3 times the
cut-off out of the domain decomposition cell of their charge group in
dimension x."

Observing the trajectories of any of the extended simulations shows the
protein snapping like a rope, and always at the same place. I have watched
every trajectory at this point, using numerous energy min start points, to
try and understand why it is blowing up. I can¹t see any obvious reason. I
was told to consider how the temperature is changing. Below is an example
of the temperature and pressure from an NPT of 20fs step continued from
the very short 20fs step NVT simulation (hoping that perhaps CG without
pressure just doesn¹t behave happily; I was wrong).


TEMP:
Š
6.630000  311.000336
    6.645000  311.371643
    6.660000  311.724213
    6.675000  313.878693
    6.690000  681558.937500


PRESSURE:
Š
6.630000    3.559879
    6.645000    3.901433
    6.660000    3.589078
    6.675000    4.158611
    6.690000  81762.437500

The final LINCS warning from this same run:

Step 300, time 4.5 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000035, max 0.003386 (between atoms 2125 and 2126)
bonds that rotated more than 45 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
   2125   2126   68.3    0.2781   0.2691      0.2700
   2125   2127   45.9    0.2789   0.2701      0.2700


At this stage the structure ruptures as described above.


My NVT settings (with NPT included to save space) are:

-----------------
title                    = Martini

integrator               = md
dt                       = 0.015
nsteps                   = 1000
nstcomm                  = 100
;comm-grps                       =

nstxout                  = 1000
nstvout                  = 1000
nstfout                  = 0
nstlog                   = 1
nstenergy                = 1
nstxout-compressed       = 0
compressed-x-precision   = 0
;compressed-x-grps        =
energygrps               = collagen solvent

cutoff-scheme            = Verlet
nstlist                  = 20
ns_type                  = grid

pbc                      = xyz
verlet-buffer-tolerance  = 0.005

coulombtype              = PME ;reaction-field
rcoulomb                 = 1.1
fourierspacing       = 0.16 ;0.2  ;0.12

epsilon_r                = 2.5 ;15      ; 2.5 (with polarizable water)
epsilon_rf               = 0
vdw_type                 = cutoff
vdw-modifier             = Potential-shift-verlet
rvdw                     = 1.1

tcoupl                   = v-rescale ;berendsen ;v-rescale
tc-grps                  = collagen solvent
tau_t                    = 0.5 0.5 ;1.0 1.0
ref_t                    = 310 310

Pcoupl                   = berendsen   ;parrinello-rahman
Pcoupltype               = isotropic
tau_p                    = 12.0  ; parrinello-rahman is more stable with
larger tau-p, DdJ, 20130422
compressibility          = 10e-4
ref_p                    = 1.0
refcoord_scaling         = com

gen_vel                  = no
gen_temp                 = 310
gen_seed                 = 473529

continuation = yes
constraints              = none
constraint_algorithm     = lincs
lincs-warnangle = 45
lincs-order=8
lincs-iter=4


‹‹‹‹‹‹‹‹‹‹

Every setting bar the lincs iter, order, warnangle were supplied with the
latest version of MARTINI. During many NVT runs I have adjusted the tau-t
to try and keep the thermostat from oscillating its way into infinity.

I¹m curious, will an out of control thermostat break a structure, or will
a structure breaking (for what ever reason this structure is breaking)
cause the thermostat to go out of control?

My only thought thought is the initial .itp file that Martinize created. I
informed the script that this was collagen, therefor it sets ³F² and the
corresponding bonded parameters. A human collagen is not perfect in its
helical structure. Could there be underlying forces contributed from badly
bonded backbone-backbone arrangements?

[ atoms ]
    1    N0     1   GLY    BB     1  0.0000 ; F
    2    N0     2   LEU    BB     2  0.0000 ; F
    3    C1     2   LEU   SC1     3  0.0000 ; F
    4    N0     3   SER    BB     4  0.0000 ; F


Many thanks
Anthony

Thanks
Anthony