[gmx-users] CG Lincs errors

Mark Abraham mark.j.abraham at gmail.com
Thu Dec 15 22:32:51 CET 2016


Hi,

If a simulation isn't stable with a small time step (as I think you are
saying) then moving to a larger time step is guaranteed to make that worse.
Try an even smaller time step, for a long time, and see what happens. Or
take a subset of your protein and see what happens. Or simulate in vacuo
for a while. Your topology could be unsuited to your starting structure,
e.g. some part is under a lot of tension that gets released at some point
and no finite time step can in practice deal with the velocity of the
recoil...

Mark

On Thu, 15 Dec 2016 23:06 Nash, Anthony <a.nash at ucl.ac.uk> wrote:

> 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
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>


More information about the gromacs.org_gmx-users mailing list