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.
