[gmx-users] CG Lincs errors
Nash, Anthony
a.nash at ucl.ac.uk
Fri Dec 16 13:59:33 CET 2016
Hi Peter,
Thanks for the reply, I know we spoke in length on this mater only just
recently. Many thanks for that.
I’ve taken the time step of collagen in vacu down to 0.0001 and I’ve
dropped the temp down to 280. I hope, running over 16 cores for two days
that this should relieve any tension in the backdone before in gradually
increase to 20-40fs. Again, all in vacu with no solvent.
I’ve actually thought about writing a script to modify the equilibrium
bond angles in the CG.itp file for the backbone, using the atomistic
structure as a template. After all, true collagen does not replicate
‘ideal’ collagen along the stretch of the protein e.g., the MMP1 binding
site is not very tightly coiled. Perhaps, there lies my problem.
I haven’t thought too much of the actual full fibril structure, I want to
capture the D-band gap in type I collagen. When I give it some more
thought I will probably look into semi-isotropic pressure coupling first.
Thanks
Anthony
Dr Anthony Nash
Department of Chemistry
University College London
Skeletal Tissue Dynamics Group
Committee member of London Matrix Group @LondonMatrixGrp
On 16/12/2016 09:02, "gromacs.org_gmx-users-bounces at maillist.sys.kth.se on
behalf of Peter Kroon" <gromacs.org_gmx-users-bounces at maillist.sys.kth.se
on behalf of p.c.kroon at rug.nl> wrote:
>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