[gmx-users] Protein in SDS, GROMOS96 ff, MD crash - SDS molecule comes apart, LINCS error
Venky Krishna
venky.agas at gmail.com
Wed Apr 25 04:15:54 CEST 2007
Hi gmx-users,
System Setup:
My System consists of a one helix of a transmembrane protein + 58 SDS
molecules + water + ions. (added in that order)
I energy minimized the system (using l-bfgs algorithm) after addition
of each set of molecules to the box, starting from just the protein
in the box. The potential energy was minimized well (to a negative
value) and the Fmax on the atoms was < 50 as stated in the em.mdp file.
Next I equilibrated the system for 1 ns by doing a position
restrained MD (protein restrained), to relax the SDS molecules around
the protein (NVT equilibration) using berendsen thermostat.
Later NPT md simulation was started by linearly heating the system
thus obtained from NVT equilibration from 0 K to 300 K in 50 ps. the
rest of the 10 ns the temperature was kept at 300 K using nose-hoover
thermostat. Pressure coupling using parrinello-rahman.
Problem:
So here is the problem... after 500 ps, a SDS molecule rips apart and
the simulation crashes with LINCS error. Quite surprisingly, there
was only one SDS molecule showing such instability... all other
molecules are intact. any inputs on how to correct the problem... for
now I just deleted the SDS molecule, equilibrated the system for 200
ps and rerun the MD simulation... I have a bad feeling it will crash
again.
Parameters:
I generated the SDS topology using PRODRG2 BETA server.. since I
wanted GROMOS96 parameters. Here is how the SDS top file looks like.
~~~~~~~~~~~~~~~~~~~~
[ moleculetype ]
; Name nrexcl
SDS 3
[ atoms ]
; nr type resnr resid atom cgnr charge mass
1 OA 1 SDS OAD 1 -0.726 15.9994
2 SDMSO 1 SDS SAQ 1 2.226 32.0600
3 OM 1 SDS OAB 1 -1.179 15.9994
4 OM 1 SDS OAC 1 -1.179 15.9994
5 OA 1 SDS OAP 1 -0.207 15.9994
6 CH2 1 SDS CAO 1 0.033 14.0270
7 CH2 1 SDS CAN 1 0.032 14.0270
8 CH2 1 SDS CAM 2 0.000 14.0270
9 CH2 1 SDS CAL 2 0.000 14.0270
10 CH2 1 SDS CAK 2 0.000 14.0270
11 CH2 1 SDS CAJ 2 0.000 14.0270
12 CH2 1 SDS CAI 3 0.007 14.0270
13 CH2 1 SDS CAH 3 0.006 14.0270
14 CH2 1 SDS CAG 3 0.008 14.0270
15 CH2 1 SDS CAF 3 0.007 14.0270
16 CH2 1 SDS CAE 3 0.007 14.0270
17 CH3 1 SDS CAA 3 -0.035 15.0350
[ bonds ]
; ai aj fu c0, c1, ...
1 2 2 0.144 6100000.0 0.144 6100000.0 ; OAD SAQ
2 3 2 0.153 8040000.0 0.153 8040000.0 ; SAQ OAB
2 4 2 0.153 8040000.0 0.153 8040000.0 ; SAQ OAC
2 5 2 0.144 6100000.0 0.144 6100000.0 ; SAQ OAP
5 6 2 0.143 8180000.0 0.143 8180000.0 ; OAP CAO
6 7 2 0.153 7150000.0 0.153 7150000.0 ; CAO CAN
7 8 2 0.153 7150000.0 0.153 7150000.0 ; CAN CAM
8 9 2 0.153 7150000.0 0.153 7150000.0 ; CAM CAL
9 10 2 0.153 7150000.0 0.153 7150000.0 ; CAL CAK
10 11 2 0.153 7150000.0 0.153 7150000.0 ; CAK CAJ
11 12 2 0.153 7150000.0 0.153 7150000.0 ; CAJ CAI
12 13 2 0.153 7150000.0 0.153 7150000.0 ; CAI CAH
13 14 2 0.153 7150000.0 0.153 7150000.0 ; CAH CAG
14 15 2 0.153 7150000.0 0.153 7150000.0 ; CAG CAF
15 16 2 0.153 7150000.0 0.153 7150000.0 ; CAF CAE
16 17 2 0.153 7150000.0 0.153 7150000.0 ; CAE CAA
[ pairs ]
; ai aj fu c0, c1, ...
1 6 1 ; OAD CAO
2 7 1 ; SAQ CAN
3 6 1 ; OAB CAO
4 6 1 ; OAC CAO
5 8 1 ; OAP CAM
6 9 1 ; CAO CAL
7 10 1 ; CAN CAK
8 11 1 ; CAM CAJ
9 12 1 ; CAL CAI
10 13 1 ; CAK CAH
11 14 1 ; CAJ CAG
12 15 1 ; CAI CAF
13 16 1 ; CAH CAE
14 17 1 ; CAG CAA
[ angles ]
; ai aj ak fu c0, c1, ...
1 2 3 2 109.5 518.0 109.5 518.0 ; OAD
SAQ OAB
1 2 4 2 109.5 518.0 109.5 518.0 ; OAD
SAQ OAC
1 2 5 2 109.5 518.0 109.5 518.0 ; OAD
SAQ OAP
3 2 4 2 109.5 518.0 109.5 518.0 ; OAB
SAQ OAC
3 2 5 2 109.5 518.0 109.5 518.0 ; OAB
SAQ OAP
4 2 5 2 109.5 518.0 109.5 518.0 ; OAC
SAQ OAP
2 5 6 2 120.0 530.0 120.0 530.0 ; SAQ
OAP CAO
5 6 7 2 109.5 520.0 109.5 520.0 ; OAP
CAO CAN
6 7 8 2 109.5 520.0 109.5 520.0 ; CAO
CAN CAM
7 8 9 2 109.5 520.0 109.5 520.0 ; CAN
CAM CAL
8 9 10 2 109.5 520.0 109.5 520.0 ; CAM
CAL CAK
9 10 11 2 109.5 520.0 109.5 520.0 ; CAL
CAK CAJ
10 11 12 2 109.5 520.0 109.5 520.0 ; CAK
CAJ CAI
11 12 13 2 109.5 520.0 109.5 520.0 ; CAJ
CAI CAH
12 13 14 2 109.5 520.0 109.5 520.0 ; CAI
CAH CAG
13 14 15 2 109.5 520.0 109.5 520.0 ; CAH
CAG CAF
14 15 16 2 109.5 520.0 109.5 520.0 ; CAG
CAF CAE
15 16 17 2 109.5 520.0 109.5 520.0 ; CAF
CAE CAA
[ dihedrals ]
; ai aj ak al fu c0, c1, m, ...
2 1 3 4 2 35.3 334.8 35.3 334.8 ; imp
SAQ OAD OAB OAC
1 2 5 6 1 0.0 1.3 3 0.0 1.3 3 ; dih
OAD SAQ OAP CAO
7 6 5 2 1 0.0 1.3 3 0.0 1.3 3 ; dih
CAN CAO OAP SAQ
8 7 6 5 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAM CAN CAO OAP
9 8 7 6 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAL CAM CAN CAO
10 9 8 7 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAK CAL CAM CAN
11 10 9 8 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAJ CAK CAL CAM
12 11 10 9 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAI CAJ CAK CAL
13 12 11 10 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAH CAI CAJ CAK
14 13 12 11 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAG CAH CAI CAJ
15 14 13 12 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAF CAG CAH CAI
16 15 14 13 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAE CAF CAG CAH
17 16 15 14 1 0.0 5.9 3 0.0 5.9 3 ; dih
CAA CAE CAF CAG
~~~~~~~~~~~~~~~~~~~~
and here are my relevant MD mdp parameters.
####################
integrator = md
tinit = 0
dt = 0.002
nsteps = 5000000
init_step = 0
nstlist = 5
ns-type = grid
pbc = xyz
rlist = 0.9
domain-decomposition = no
coulombtype = pme
rcoulomb = 0.9
epsilon-r = 1
vdw-type = Cut-off
rvdw-switch = 0
rvdw = 1.4
DispCorr = No
table-extension = 1
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = yes
; OPTIONS FOR WEAK COUPLING ALGORITHMS
tcoupl = nose-hoover
tc-grps = Protein SDS SOL Na+ Cl-
tau-t = 0.5 0.5 0.5 0.5 0.5
ref-t = 300 300 300 300 300
Pcoupl = Parrinello-Rahman
Pcoupltype = Isotropic
tau-p = 5
compressibility = 4.5e-5
ref-p = 1
andersen_seed = 815131
; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing = single single single single single
; Number of time points to use for specifying annealing in each group
annealing_npoints = 2 2 2 2 2
; List of times at the annealing points for each group
annealing_time = 0 20 0 20 0 20 0 20 0 20
; Temp. at each annealing point, for each group.
annealing_temp = 0 300 0 300 0 300 0 300 0 300
; GENERATE VELOCITIES FOR STARTUP RUN
gen-vel = yes
gen-temp = 0
gen-seed = 3529
; OPTIONS FOR BONDS
constraints = all-bonds
constraint-algorithm = Lincs
unconstrained-start = no
lincs-order = 4
lincs-iter = 1
lincs-warnangle = 30
morse = no
####################
Here is the ERROR....
$$$$$$$$$$$$$$$$$$$$$$$$$$$
Step 257327, time 514.654 (ps) LINCS WARNING
relative constraint deviation after LINCS:
max 25.765928 (between atoms 920 and 921) rms 1.109525
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
920 921 92.7 0.1507 3.8543 0.1440
921 922 119.1 0.1556 0.3664 0.1530
921 923 78.8 0.1442 0.1825 0.1530
921 924 100.1 0.1613 1.1378 0.1440
924 925 92.9 0.1390 3.1287 0.1430
925 926 92.6 0.1870 3.1463 0.1530
926 927 157.8 0.1521 0.1948 0.1530
927 928 128.2 0.1545 0.2347 0.1530
928 929 63.1 0.1529 0.1001 0.1530
Constraint error in algorithm Lincs at step 257327
Step 257328, time 514.656 (ps) LINCS WARNING
relative constraint deviation after LINCS:
max 390842875904.000000 (between atoms 924 and 925) rms
20729815040.000000
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
699 700 118.2 0.1440 0.3004 0.1440
700 701 99.2 0.1530 0.9467 0.1530
700 702 62.8 0.1530 0.1665 0.1530
700 703 77.7 0.1440 0.1298 0.1440
703 704 31.6 0.1430 0.1368 0.1430
716 717 90.5 0.1440 34.8807 0.1440
717 718 90.1 0.1530 147.2678 0.1530
717 719 90.1 0.1530 110.9749 0.1530
717 720 90.9 0.1440 36.1724 0.1440
720 721 91.5 0.1430 10.8505 0.1430
721 722 97.2 0.1530 2.9375 0.1530
722 723 145.7 0.1530 0.2025 0.1530
723 724 42.0 0.1530 0.1993 0.1530
886 887 57.6 0.1440 0.1784 0.1440
887 888 67.3 0.1530 0.2050 0.1530
887 889 63.1 0.1530 0.1974 0.1530
887 890 100.5 0.1440 0.8066 0.1440
890 891 100.0 0.1430 0.7412 0.1430
891 892 155.6 0.1530 0.1709 0.1530
892 893 42.1 0.1530 0.2004 0.1530
903 904 90.9 0.1440 15469.8779 0.1440
904 905 90.1 0.1530 51128.5117 0.1530
904 906 90.2 0.1530 46422.8555 0.1530
904 907 90.1 0.1440 134111.4688 0.1440
907 908 90.1 0.1430 168917.5156 0.1430
908 909 90.3 0.1530 16148.1357 0.1530
909 910 93.8 0.1530 687.8440 0.1530
910 911 93.1 0.1530 101.2424 0.1530
911 912 108.5 0.1530 14.7890 0.1530
913 914 130.7 0.1530 0.6334 0.1530
914 915 59.3 0.1530 0.2331 0.1530
920 921 124.5 3.8543 590567680.0000 0.1440
921 922 115.9 0.3664 507195392.0000 0.1530
921 923 127.3 0.1825 597398592.0000 0.1530
921 924 89.6 1.1378 54504730624.0000 0.1440
924 925 94.0 3.1287 55890534400.0000 0.1430
925 926 94.6 3.1463 52505100288.0000 0.1530
926 927 91.7 0.1948 53183660032.0000 0.1530
927 928 50.6 0.2347 1845543040.0000 0.1530
928 929 58.1 0.1001 2104577408.0000 0.1530
929 930 162.9 0.1506 350578240.0000 0.1530
930 931 31.0 0.1543 339091040.0000 0.1530
931 932 86.4 0.1530 0.3569 0.1530
932 933 44.6 0.1531 0.2167 0.1530
1117 1118 65.7 0.1530 0.2308 0.1530
1118 1119 106.3 0.1530 0.5409 0.1530
1119 1120 104.2 0.1530 0.6235 0.1530
1120 1121 98.6 0.1530 1.0182 0.1530
1121 1122 95.9 0.1530 1.4904 0.1530
1122 1123 95.0 0.1530 1.7516 0.1530
1124 1125 107.3 0.1440 0.4820 0.1440
1125 1126 107.7 0.1530 0.5021 0.1530
1125 1127 107.1 0.1530 0.5176 0.1530
1125 1128 104.5 0.1440 0.5717 0.1440
1128 1129 78.4 0.1430 0.2178 0.1430
1129 1130 63.6 0.1530 0.1758 0.1530
1130 1131 53.0 0.1530 0.1784 0.1530
1131 1132 40.3 0.1530 0.1644 0.1530
Constraint error in algorithm Lincs at step 257328
t = 514.656 ps: Water molecule starting at atom 18916 can not be
settled.
Check for bad contacts and/or reduce the timestep.Wrote pdb files
with previous and current coordinates
Large VCM(group rest): 24544858112.00000, -56862011392.00000,
610747392.00000, ekin-cm: 2.45315e+26
Group rest with mass 1.27898e+05, Ekrot 1.46688e+15 Det(I)
= nan
COM: 49089732.00000 -113724048.00000 1221496.62500
P: 3139232935706624.00000 -7272524777783296.00000
78113232912384.00000
V: 24544858112.00000 -56862011392.00000 610747392.00000
J: -162155133050814464.00000 83557368560651075584.00000
-184521431156065828864.00000
w: -0.00001 0.00002 -0.00001
Inertia tensor (3x3):
Inertia tensor[ 0]={ 1.59138e+25, -2.53183e+25, 6.93924e+24}
Inertia tensor[ 1]={-2.53183e+25, 4.38301e+25, -1.26821e+25}
Inertia tensor[ 2]={ 6.93924e+24, -1.26821e+25, 8.36827e+24}
$$$$$$$$$$$$$$$$$$$$$$$$$$
Regards,
Venky.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20070424/d4197129/attachment.html>
More information about the gromacs.org_gmx-users
mailing list