[gmx-users] CG Lincs errors
Alex
nedomacho at gmail.com
Thu Dec 15 22:54:23 CET 2016
Mark is right, no two ways about it. For initial equilibration and
assessing preexisting structural strains try vacuum, _much_ smaller
timesteps and possibly low temperatures in vacuum, only then transfer to
solvent, etc. Algorithmically, LINCS requires convergence and you already
are using a pretty high LINCS order... From what I see, dt = 15 fs at 310K
looks like a cowboy mode simulation in this case.
Alex
On Thu, Dec 15, 2016 at 2:32 PM, Mark Abraham <mark.j.abraham at gmail.com>
wrote:
> 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.
> >
> --
> 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