[gmx-users] MARTINI force-field in pulling simulations
Justin Lemkul
jalemkul at vt.edu
Thu Feb 7 20:09:47 CET 2013
On 2/7/13 1:44 PM, Davide Mercadante wrote:
> Dear Justin,
>
> Thank you for your reply. I decreased the time step from 0.02 to 0.005 and
> run the simulation again. The simulation still crashes giving LINCS
> warning on the same atoms but does it later.
> Do you advice to keep reducing the time step in order to reach a simulated
> time where the pull of the whole molecule occurs and I see the force
> peaking? I am still not sure why this happens...
>
No, I think the issue is probably more fundamental. MARTINI uses fixed
secondary structure (bonds in the topology to preserve geometry). You're
probably just pulling against those and the algorithms that work on bonded
interactions are failing.
-Justin
> Thanks.
>
> Davide
>
> On 7/02/13 1:05 PM, "Justin Lemkul" <jalemkul at vt.edu> wrote:
>
>>
>>
>> On 2/7/13 5:20 AM, Davide Mercadante wrote:
>>> Dear All,
>>>
>>> I am trying to run a pulling simulation on a small protein (18 aa) using
>>> the GC forcefield MARTINI (v2.2). I have energy minimized and
>>> equilibrated (NPT) my system and everything seems fine. My system
>>> consists of the protein + water + ions NA+ and CL-.
>>>
>>> After the equilibration I start a constant velocity pulling along the
>>> z-direction of the last atom of the chain while the first atom is
>>> positionally restrained in xyz (basically I am stretching the protein).
>>> At some point from the start of the pulling simulation and before
>>> reaching the full extension of the chain (force profile is still
>>> steadily increasing without peaking) the simulation crashes giving me
>>> these LINCS warnings:
>>>
>>> Step 293252, time 5865.04 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.157274, max 0.291135 (between atoms 10 and 11)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>>
>>> Step 293252, time 5865.04 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.128562, max 0.230219 (between atoms 7 and 8)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>>
>>> Step 293253, time 5865.06 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.413291, max 0.828515 (between atoms 7 and 8)
>>>
>>> Step 293253, time 5865.06 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.587032, max 0.986890 (between atoms 11 and 13)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 7 8 46.6 0.2686 0.0579 0.3500
>>> 10 11 121.3 0.2481 0.0550 0.3500
>>> 13 15 168.5 0.2851 0.0278 0.3500
>>>
>>> Step 293253, time 5865.06 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.108373, max 0.322423 (between atoms 19 and 21)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 7 8 45.0 0.2686 0.0600 0.3500
>>> 10 11 88.1 0.2481 0.0497 0.3500
>>> Wrote pdb files with previous and current coordinates
>>> Wrote pdb files with previous and current coordinates
>>>
>>> Step 293254, time 5865.08 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.330732, max 0.654588 (between atoms 23 and 25)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 21 23 42.0 0.2905 0.1426 0.3500
>>> 17 18 92.9 0.3146 0.1560 0.3000
>>> 15 17 35.3 0.5839 0.4003 0.3500
>>> 13 15 141.6 0.0278 0.2742 0.3500
>>>
>>> Step 293254, time 5865.08 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.312540, max 0.555587 (between atoms 19 and 21)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 8 9 48.0 0.3406 0.1435 0.3100
>>> 7 8 82.3 0.0579 0.4596 0.3500
>>> 10 11 102.9 0.0550 0.5130 0.3500
>>> 11 12 51.9 0.5425 0.3643 0.4000
>>> 11 13 38.5 0.6954 0.1898 0.3500
>>> 13 15 143.4 0.0278 0.2879 0.3500
>>> 15 17 37.3 0.5839 0.3730 0.3500
>>>
>>> Step 293254, time 5865.08 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.343084, max 0.751867 (between atoms 3 and 5)
>>> 17 18 93.2 0.3146 0.1563 0.3000
>>> 21 23 41.6 0.2905 0.1443 0.3500
>>> 23 24 33.0 0.3784 0.4007 0.4000
>>> 11 13 58.2 0.6954 0.3283 0.3500
>>> 23 24 33.1 0.3784 0.3999 0.4000
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 7 8 82.2 0.0579 0.4605 0.3500
>>> 8 9 48.4 0.3406 0.1439 0.3100
>>> 10 11 104.0 0.0550 0.5496 0.3500
>>> 11 12 57.6 0.5425 0.3679 0.4000
>>> 11 13 61.9 0.6954 0.1402 0.3500
>>> 13 15 141.2 0.0278 0.4588 0.3500
>>> Wrote pdb files with previous and current coordinates
>>> Wrote pdb files with previous and current coordinates
>>> Wrote pdb files with previous and current coordinates
>>>
>>> Step 293255, time 5865.1 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 32.745333, max 105.228383 (between atoms 17 and 18)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 5 6 30.3 0.3950 0.3971 0.4000
>>> 13 15 172.9 0.2879 4.6287 0.3500
>>> 15 16 33.0 0.2435 5.0674 0.3300
>>> 15 17 51.2 0.3730 11.1793 0.3500
>>> 17 18 107.1 0.1563 31.8685 0.3000
>>> 17 19 143.4 0.3380 10.1576 0.3500
>>> 19 20 86.5 0.3397 6.2047 0.3300
>>> 21 22 103.7 0.3750 6.6220 0.4000
>>> 21 23 150.3 0.1426 6.8152 0.3500
>>>
>>> Step 293255, time 5865.1 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 10.273318, max 30.025748 (between atoms 19 and 21)
>>> 23 24 149.4 0.3999 1.3487 0.4000
>>> 25 26 79.9 0.3039 2.4910 0.3000
>>> 25 27 165.9 0.2225 2.4307 0.3500
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 21 22 103.7 0.3750 6.6180 0.4000
>>> 21 23 150.2 0.1426 6.8083 0.3500
>>> 19 20 86.5 0.3397 6.2048 0.3300
>>> 17 19 143.4 0.3380 10.1586 0.3500
>>> 17 18 107.1 0.1563 31.8685 0.3000
>>> 15 17 51.2 0.3730 11.1771 0.3500
>>> 15 16 33.0 0.2435 5.0657 0.3300
>>> 13 15 173.1 0.2879 4.6383 0.3500
>>> 23 24 149.5 0.3999 1.3474 0.4000
>>> 25 26 79.9 0.3039 2.4804 0.3000
>>> 25 27 165.7 0.2225 2.3691 0.3500
>>> 28 29 44.1 0.3191 0.0921 0.3100
>>> 28 30 39.4 0.3101 0.2122 0.3500
>>> Wrote pdb files with previous and current coordinates
>>> Wrote pdb files with previous and current coordinates
>>>
>>> Step 293256 Warning: pressure scaling more than 1%, mu: 2.21304 2.21304
>>> 2.21304
>>>
>>> Step 293256 Warning: pressure scaling more than 1%, mu: 2.21304 2.21304
>>> 2.21304
>>>
>>> Step 293256 Warning: pressure scaling more than 1%, mu: 2.21304 2.21304
>>> 2.21304
>>>
>>> Step 293256 Warning: pressure scaling more than 1%, mu: 2.21304 2.21304
>>> 2.21304
>>>
>>> Step 293256 Warning: pressure scaling more than 1%, mu: 2.21304 2.21304
>>> 2.21304
>>>
>>> Step 293256 Warning: pressure scaling more than 1%, mu: 2.21304 2.21304
>>> 2.21304
>>>
>>> Step 293256 Warning: pressure scaling more than 1%, mu: 2.21304 2.21304
>>> 2.21304
>>>
>>> Step 293256 Warning: pressure scaling more than 1%, mu: 2.21304 2.21304
>>> 2.21304
>>>
>>> Step 293256, time 5865.12 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 3192.715629, max 6000.654094 (between atoms 25 and 26)
>>>
>>> Step 293256, time 5865.12 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.212149, max 0.517573 (between atoms 5 and 7)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 5 7 37.9 0.3543 0.5312 0.3500
>>> 7 8 33.5 0.3821 0.3834 0.3500
>>> 8 9 56.8 0.3361 0.3502 0.3100
>>> 8 10 165.7 0.3548 1.0526 0.3500
>>>
>>> Step 293256, time 5865.12 (ps) LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 95.799929, max 246.403254 (between atoms 19 and 21)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 8 9 56.7 0.3361 0.3566 0.3100
>>> 8 10 165.7 0.3548 1.0708 0.3500
>>> 7 8 33.4 0.3821 0.3897 0.3500
>>> 5 7 37.9 0.3543 0.5312 0.3500
>>> 11 12 121.6 0.2312 2.6433 0.4000
>>> 11 13 149.8 0.3060 15.6610 0.3500
>>> 13 14 123.6 0.4129 14.6116 0.4000
>>> 13 14 124.2 0.4129 14.9085 0.4000
>>> 13 15 169.6 4.6287 44.1184 0.3500
>>> 11 12 121.1 0.2312 2.7328 0.4000
>>> 11 13 149.6 0.3060 15.9678 0.3500
>>> 15 16 168.1 5.0674 39.6694 0.3300
>>> 15 17 159.6 11.1793 63.3062 0.3500
>>> 17 18 114.1 31.8685 33.0371 0.3000
>>> 17 19 112.0 10.1576 36.5026 0.3500
>>> 19 20 71.0 6.2047 9.2619 0.3300
>>> 23 24 62.3 1.3474 352.2061 0.4000
>>> 25 26 63.9 2.4804 1800.5086 0.3000
>>> 25 27 67.1 2.3691 2007.2339 0.3500
>>> 13 15 169.4 4.6287 43.0156 0.3500
>>> 27 28 66.0 0.6878 1992.3756 0.3500
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2 angle previous, current, constraint length
>>> 19 20 71.0 6.2047 9.2610 0.3300
>>> 17 19 112.0 10.1576 36.4983 0.3500
>>> 17 18 114.1 31.8685 33.0303 0.3000
>>> 15 17 159.6 11.1793 63.3215 0.3500
>>> 15 16 168.1 5.0674 39.6801 0.3300
>>> 13 15 169.6 4.6287 44.0779 0.3500
>>> 13 14 124.6 0.4129 14.7983 0.4000
>>> 11 13 148.8 0.3060 16.1472 0.3500
>>> 23 24 62.3 1.3474 352.1905 0.4000
>>> 25 26 63.9 2.4804 1800.4962 0.3000
>>> 25 27 67.1 2.3691 2007.1944 0.3500
>>> 27 28 66.0 0.6878 1992.3052 0.3500
>>> 28 29 44.6 0.0921 600.7400 0.3100
>>> 28 30 137.0 0.2122 696.2760 0.3500
>>> 30 31 38.7 0.2803 155.6825 0.2650
>>> Wrote pdb files with previous and current coordinates
>>> Wrote pdb files with previous and current coordinates
>>>
>>> I have a feeling that this happens when I am reaching the full extension
>>> of the chain because if in the input file (pull.mdp - copied below) I
>>> increase the pull rate the simulation simply crashes before.
>>> This crash does not happen in a full-atom simulation of the same protein
>>> (using for example OPLS) and the force profile reaches the peak wanted
>>> when the protein approaches the full extension.
>>>
>>> I was wondering if you please could help me to understand why this
>>> happens, if this is a force-field related issue and how I could possibly
>>> solve the problem.
>>> The force costant I apply in the posre.itp file for atom number 1 is
>>> 5000 kj mol^1 nm^2 while the pulling force constant (harmonic potential)
>>> is 500 kj mol^1 nm^2.
>>> The pull.mdp input file is appended below.
>>>
>>> title = Umbrella pulling simulation
>>> define = -DPOSRES
>>> ; Run parameters
>>> integrator = md
>>> dt = 0.02
>>> tinit = 0
>>> nsteps = 500000 ; 10000.00 ps
>>> nstcomm = 10
>>> ; Output parameters
>>> nstxout = 100 ; every 2 ps
>>> nstvout = 100
>>> nstfout = 100
>>> nstxtcout = 100 ; every 2 ps
>>> nstenergy = 100
>>> ; Bond parameters
>>> constraint_algorithm = lincs
>>> constraints = all-bonds
>>> continuation = yes ; continuing from NPT
>>> ; Single-range cutoff scheme
>>> nstlist = 5
>>> ns_type = grid
>>> rlist = 1.4
>>> rcoulomb = 1.4
>>> rvdw = 1.4
>>> ; PME electrostatics parameters
>>> coulombtype = PME
>>> fourierspacing = 0.12
>>> fourier_nx = 0
>>> fourier_ny = 0
>>> fourier_nz = 0
>>> pme_order = 4
>>> ewald_rtol = 1e-5
>>> optimize_fft = yes
>>> ; Berendsen temperature coupling is on in two groups
>>> Tcoupl = V-rescale
>>> tc_grps = Protein Non-Protein
>>> tau_t = 2 2
>>> ref_t = 310 310
>>> ; Pressure coupling is on
>>> Pcoupl = Berendsen
>>> pcoupltype = isotropic
>>> tau_p = 2.0
>>> compressibility = 4.5e-5
>>> ref_p = 1.0
>>> refcoord_scaling = com
>>> ; Generate velocities is off
>>> gen_vel = no
>>> ; Periodic boundary conditions are on in all directions
>>> pbc = xyz
>>> ; Long-range dispersion correction
>>> DispCorr = EnerPres
>>>
>>> ; PULL-CODE
>>>
>>> pull = umbrella
>>> pull_geometry = direction
>>> pull_start = yes ; (we don't want to calculate the initial distance
>>> between 'pull' and 'freeze' by hand)
>>> pull_ngroups = 1 ; (because the reference group doesn't count to this
>>> value)
>>> pull_dim = N N Y ; (because we want to pull only along the z-axis)
>>> pull_group0 = freeze
>>> pull_group1 = pull
>>> pull_vec1 = 0 0 1
>>> pull_init1 = 0.0 0.0 0.0
>>> pull_rate1 = 0.001
>>> pull_k1 = 500
>>> pull_constr_tol = 1e-6
>>> pull_pbc_atom0 = 0
>>> pull_kB1 = 500
>>> pull_nstxout = 100 ; every 2 ps
>>> pull_nstfout = 100 ; every 2 ps
>>>
>>> Thank you, any suggestions will be greatly appreciated.
>>
>> Try reducing the time step.
>>
>> -Justin
>>
>> --
>> ========================================
>>
>> Justin A. Lemkul, Ph.D.
>> Research Scientist
>> Department of Biochemistry
>> Virginia Tech
>> Blacksburg, VA
>> jalemkul[at]vt.edu | (540) 231-9080
>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>>
>> ========================================
>> --
>> gmx-users mailing list gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-request at gromacs.org.
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>
--
========================================
Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list