[gmx-users] CG Lincs errors
Nash, Anthony
a.nash at ucl.ac.uk
Fri Dec 16 00:21:44 CET 2016
Alex and Mark, thanks for the information. I’ll drop dt down,
significantly, drop the temperature, and run it for a long time.
Thanks for the ideas.
Anthony
On 15/12/2016 21:54, "gromacs.org_gmx-users-bounces at maillist.sys.kth.se on
behalf of Alex" <gromacs.org_gmx-users-bounces at maillist.sys.kth.se on
behalf of nedomacho at gmail.com> wrote:
>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.
>--
>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