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. Davide -- gmx-users mailing list gmx-users@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-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists