[gmx-users] CG Lincs errors
Peter Kroon
p.c.kroon at rug.nl
Fri Dec 16 10:02:55 CET 2016
As a note to Alex (and the rest of the list), the coarse-grained Martini
forcefield is usually run with timesteps between 20-40 fs. 15fs is
already rather low. I do agree that longer equilibration at low timestep
(5 or 10) might help.
Alternatively, Do you think a semiisotropic pressure coupling might be
applicable in this case, since it's an infinite collagen polymer?
Peter (PhD in the Martini group)
On 16-12-16 00:21, Nash, Anthony wrote:
> 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