Here's an excellent learning experience.
1. Simplify the problem: minimize and attempt to equilibrate one single hexane
molecule before trying to build your stacked structure. If the system blows up
here, then the topology is likely to blame. More on this in a few moments.
2. Watch the trajectory of the explosion, setting nstxout = 1 if necessary to
catch all of the details. If you do this (with your input files below) you will
see hydrogen atoms crashing in on themselves, severe angle distortions, and
ultimately, collapse of the system. I would strongly suggest you do this, for
your own education - knowing what to look for is half the battle.
3. Consider why this might be happening. Could x2top have given you a flawed
output file? Quite possibly. Have a look at your [dihedrals] directive. Are
all possible dihedrals covered? Clearly not. Every dihedral (H-C-C-H, C-C-C-C,
C-C-C-H) needs to be accounted for in OPLS. You have only 5, when you should
have 45 dihedral terms!
Here's how I figured this out. With a molecule as simple as hexane, writing an
.rtp entry is trivial, like so:
[ HEX ]
[ atoms ]
CAA opls_157 -0.180 1
HA1 opls_140 0.060 1
HA2 opls_140 0.060 1
HA3 opls_140 0.060 1
CAC opls_158 -0.120 2
HC1 opls_140 0.060 2
HC2 opls_140 0.060 2
CAE opls_158 -0.120 3
HE1 opls_140 0.060 3
HE2 opls_140 0.060 3
CAF opls_158 -0.120 4
HF1 opls_140 0.060 4
HF2 opls_140 0.060 4
CAD opls_158 -0.120 5
HD1 opls_140 0.060 5
HD2 opls_140 0.060 5
CAB opls_157 -0.180 6
HB1 opls_140 0.060 6
HB2 opls_140 0.060 6
HB3 opls_140 0.060 6
[ bonds ]
CAA HA1
CAA HA2
CAA HA3
CAA CAC
CAC HC1
CAC HC2
CAC CAE
CAE HE1
CAE HE2
CAE CAF
CAF HF1
CAF HF2
CAF CAD
CAD HD1
CAD HD2
CAD CAB
CAB HB1
CAB HB2
CAB HB3
(in conjunction with the following .pdb file):
HETATM 1 CAA HEX 1 8.330 1.510 -0.010 1.00 20.00 C
HETATM 2 HA1 HEX 1 9.281 1.200 -0.024 1.00 20.00 H
HETATM 3 HA2 HEX 1 8.154 2.080 -0.813 1.00 20.00 H
HETATM 4 HA3 HEX 1 8.166 2.044 0.820 1.00 20.00 H
HETATM 5 CAC HEX 1 7.400 0.300 -0.030 1.00 20.00 C
HETATM 6 HC1 HEX 1 7.584 -0.267 0.773 1.00 20.00 H
HETATM 7 HC2 HEX 1 7.573 -0.231 -0.860 1.00 20.00 H
HETATM 8 CAE HEX 1 5.940 0.730 -0.010 1.00 20.00 C
HETATM 9 HE1 HEX 1 5.754 1.291 -0.816 1.00 20.00 H
HETATM 10 HE2 HEX 1 5.769 1.266 0.817 1.00 20.00 H
HETATM 11 CAF HEX 1 5.010 -0.480 -0.020 1.00 20.00 C
HETATM 12 HF1 HEX 1 5.192 -1.038 0.790 1.00 20.00 H
HETATM 13 HF2 HEX 1 5.186 -1.020 -0.843 1.00 20.00 H
HETATM 14 CAD HEX 1 3.540 -0.050 -0.010 1.00 20.00 C
HETATM 15 HD1 HEX 1 3.357 0.507 -0.820 1.00 20.00 H
HETATM 16 HD2 HEX 1 3.363 0.489 0.813 1.00 20.00 H
HETATM 17 CAB HEX 1 2.610 -1.270 -0.020 1.00 20.00 C
HETATM 18 HB1 HEX 1 1.658 -0.964 -0.013 1.00 20.00 H
HETATM 19 HB2 HEX 1 2.780 -1.812 -0.843 1.00 20.00 H
HETATM 20 HB3 HEX 1 2.785 -1.830 0.790 1.00 20.00 H
Then let pdb2gmx do the hard work for you. It will write all the necessary
bonded terms, simply by specifying the bonds.
Energy minimization still yields a positive potential, but it is quite low. In
this case, that's alright, since there are numerous unsatisfied hydrophobic
interactions that will likely be satisfied once there are a few other hexane
molecules around. Equilibration works just fine after that, although I would
seriously recommend *against* using cutoff electrostatics. The artefacts are
well-established.
Hopefully this has given you (and others who may come upon this post) some
insight into how to effectively diagnose problems like this one. With more
complex molecules, writing .rtp entries is not so trivial, but parameterization
in general is a very advanced concept. Knowing what the force field requires is
the biggest battle of all.
-Justin
Moeed wrote:
Dear experts,
Could you please take a look at energy values. The system contains only
stack of hexane molecules. MD sun gives lincs warning.
1- How can I get rid of high repulsive potential?
*genconf_d -f Hexane.gro -nbox 4.0 8.0 8.0 -dist 1.1 0.55 0.55 -o
Hexane-stack.gro*
********------------------------***************************************************************************
energy minimization
my output.mdrun_em:
Steepest Descents converged to Fmax < 1000 in 30 steps
Potential Energy = 4.76092783832156e+05
Maximum force = 9.86600079729483e+02 on atom 3115
Norm of force = 5.33725886931530e+02
*********************************************************************************************************
g_energy file:
Statistics over 61 steps [ 0.0000 thru 0.1200 ps ], 12 data sets
All averages are exact over 61 steps
Energy Average RMSD Fluct. Drift
Tot-Drift
------------------------------
-------------------------------------------------
Angle 78416.5 38731.2 25083.3 838187
102231
LJ-14 1.16009e+08 8.87401e+08 8.87401e+08
3.06025e+06 373250
Coulomb-14 734.065 480.683 140.481 13056.3
1592.44
LJ (SR) 5482.29 5991.07 4445.16 114080
13914.1
Coulomb (SR) 1711.47 732.288 228.708 -19758
-2409.83
Potential 1.16326e+08 8.8739e+08 8.8739e+08
921124 112347
Kinetic En. 1.77327e+18 5.4851e+18 4.28286e+18 9.73294e+19
1.1871e+19
Total Energy 1.77327e+18 5.4851e+18 4.28286e+18 9.73294e+19
1.1871e+19
Temperature 4.06507e+16 1.25741e+17 9.81811e+16 2.2312e+18
2.72133e+17
Pressure (bar) 1.48406e+15 4.59052e+15 3.58436e+15 8.14557e+16
9.93493e+15
T-HEX 4.06507e+16 1.25741e+17 9.81811e+16 2.2312e+18
2.72133e+17
Lamb-HEX 0.99144 0.00206848 0 -0.0617392
-0.00753016
Heat Capacity Cv: -0.934076 J/mol K (factor = 9.56799)
***************
*************************************************************************
title = Hexane
cpp = /lib/cpp
; Run control
integrator = md
dt = 0.002 ; ps !
nsteps = 5000 ; total 1.0 ps.
nstcomm = 1 ; frequency for center of mass motion
removal
; Output control
nstenergy = 10 ; frequency to write energies to energy
file. i.e., energies and other statistical data are stored every 10 steps
nstxout = 10 ; frequency to write
coordinates/velocity/force to output trajectory file
nstvout = 0
nstfout = 10
nstlog = 10 ; frequency to write energies to log file
; Neighbor searching
nstlist = 10 ; neighborlist will be updated at least
every 10 steps
;ns_type = grid
; Electrostatics/VdW
coulombtype = cut-off
vdw-type = cut-off
; Cut-offs
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Temperature coupling Berendsen temperature coupling is on in
two groups
Tcoupl = berendsen
tc-grps = HEX ;sol
tau_t = 0.1 ;0.1
ref_t = 300 ;300
; Pressure coupling: Pressure coupling is not on
Pcoupl = no
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Velocity generation Generate velocites is on at 300 K.
Manual p155
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
; Bonds
constraints = all-bonds
constraint-algorithm = lincs
pbc=xyz
****************************************************************************************
;
; File 'Hexane.top' was generated
; By user: moeed (500)
; On host: moeed-desktop
; At date: Thu Apr 8 13:51:19 2010
;
; This is your include topology file
; Generated by x2top
;
; Include forcefield parameters
#include "ffoplsaa.itp"
[ moleculetype ]
; Name nrexcl
HEX 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
typeB chargeB massB
1 opls_157 1 HEX C1 1 -0.18 12.011 ;
qtot -0.18
2 opls_158 1 HEX C2 2 -0.12 12.011 ;
qtot -0.3
3 opls_158 1 HEX C3 3 -0.12 12.011 ;
qtot -0.42
4 opls_158 1 HEX C4 4 -0.12 12.011 ;
qtot -0.54
5 opls_158 1 HEX C5 5 -0.12 12.011 ;
qtot -0.66
6 opls_157 1 HEX C6 6 -0.18 12.011 ;
qtot -0.84
7 opls_140 1 HEX H1 6 0.06 1.008 ;
qtot -0.78
8 opls_140 1 HEX H2 6 0.06 1.008 ;
qtot -0.72
9 opls_140 1 HEX H3 6 0.06 1.008 ;
qtot -0.66
10 opls_140 1 HEX H4 5 0.06 1.008 ;
qtot -0.6
11 opls_140 1 HEX H5 5 0.06 1.008 ;
qtot -0.54
12 opls_140 1 HEX H6 4 0.06 1.008 ;
qtot -0.48
13 opls_140 1 HEX H7 4 0.06 1.008 ;
qtot -0.42
14 opls_140 1 HEX H8 3 0.06 1.008 ;
qtot -0.36
15 opls_140 1 HEX H9 3 0.06 1.008 ;
qtot -0.3
16 opls_140 1 HEX H10 2 0.06 1.008 ;
qtot -0.24
17 opls_140 1 HEX H11 2 0.06 1.008 ;
qtot -0.18
18 opls_140 1 HEX H12 1 0.06 1.008 ;
qtot -0.12
19 opls_140 1 HEX H13 1 0.06 1.008 ;
qtot -0.06
20 opls_140 1 HEX H14 1 0.06 1.008 ;
qtot 0
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1 1.530000e-01 4.000000e+05 1.530000e-01 4.000000e+05
1 18 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
1 19 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
1 20 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
2 3 1 1.530000e-01 4.000000e+05 1.530000e-01 4.000000e+05
2 16 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
2 17 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
3 4 1 1.540000e-01 4.000000e+05 1.540000e-01 4.000000e+05
3 14 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
3 15 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
4 5 1 1.530000e-01 4.000000e+05 1.530000e-01 4.000000e+05
4 12 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
4 13 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
5 6 1 1.530000e-01 4.000000e+05 1.530000e-01 4.000000e+05
5 10 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
5 11 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
6 7 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
6 8 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
6 9 1 1.090000e-01 4.000000e+05 1.090000e-01 4.000000e+05
[ pairs ]
; ai aj funct c0 c1 c2 c3
1 4 1
1 14 1
1 15 1
2 5 1
2 12 1
2 13 1
3 6 1
3 10 1
3 11 1
3 18 1
3 19 1
3 20 1
4 7 1
4 8 1
4 9 1
4 16 1
4 17 1
5 14 1
5 15 1
6 12 1
6 13 1
7 10 1
7 11 1
8 10 1
8 11 1
9 10 1
9 11 1
10 12 1
10 13 1
11 12 1
11 13 1
12 14 1
12 15 1
13 14 1
13 15 1
14 16 1
14 17 1
15 16 1
15 17 1
16 18 1
16 19 1
16 20 1
17 18 1
17 19 1
17 20 1
[ angles ]
; ai aj ak funct c0 c1
c2 c3
2 1 18 1 1.120000e+02 4.000000e+02 1.120000e+02
4.000000e+02
2 1 19 1 1.110000e+02 4.000000e+02 1.110000e+02
4.000000e+02
2 1 20 1 1.110000e+02 4.000000e+02 1.110000e+02
4.000000e+02
18 1 19 1 1.070000e+02 4.000000e+02 1.070000e+02
4.000000e+02
18 1 20 1 1.070000e+02 4.000000e+02 1.070000e+02
4.000000e+02
19 1 20 1 1.080000e+02 4.000000e+02 1.080000e+02
4.000000e+02
1 2 3 1 1.130000e+02 4.000000e+02 1.130000e+02
4.000000e+02
1 2 16 1 1.090000e+02 4.000000e+02 1.090000e+02
4.000000e+02
1 2 17 1 1.090000e+02 4.000000e+02 1.090000e+02
4.000000e+02
3 2 16 1 1.100000e+02 4.000000e+02 1.100000e+02
4.000000e+02
3 2 17 1 1.100000e+02 4.000000e+02 1.100000e+02
4.000000e+02
16 2 17 1 1.060000e+02 4.000000e+02 1.060000e+02
4.000000e+02
2 3 4 1 1.130000e+02 4.000000e+02 1.130000e+02
4.000000e+02
2 3 14 1 1.100000e+02 4.000000e+02 1.100000e+02
4.000000e+02
2 3 15 1 1.100000e+02 4.000000e+02 1.100000e+02
4.000000e+02
4 3 14 1 1.090000e+02 4.000000e+02 1.090000e+02
4.000000e+02
4 3 15 1 1.090000e+02 4.000000e+02 1.090000e+02
4.000000e+02
14 3 15 1 1.060000e+02 4.000000e+02 1.060000e+02
4.000000e+02
3 4 5 1 1.130000e+02 4.000000e+02 1.130000e+02
4.000000e+02
3 4 12 1 1.090000e+02 4.000000e+02 1.090000e+02
4.000000e+02
3 4 13 1 1.090000e+02 4.000000e+02 1.090000e+02
4.000000e+02
5 4 12 1 1.100000e+02 4.000000e+02 1.100000e+02
4.000000e+02
5 4 13 1 1.100000e+02 4.000000e+02 1.100000e+02
4.000000e+02
12 4 13 1 1.060000e+02 4.000000e+02 1.060000e+02
4.000000e+02
4 5 6 1 1.130000e+02 4.000000e+02 1.130000e+02
4.000000e+02
4 5 10 1 1.100000e+02 4.000000e+02 1.100000e+02
4.000000e+02
4 5 11 1 1.100000e+02 4.000000e+02 1.100000e+02
4.000000e+02
6 5 10 1 1.090000e+02 4.000000e+02 1.090000e+02
4.000000e+02
6 5 11 1 1.090000e+02 4.000000e+02 1.090000e+02
4.000000e+02
10 5 11 1 1.060000e+02 4.000000e+02 1.060000e+02
4.000000e+02
5 6 7 1 1.110000e+02 4.000000e+02 1.110000e+02
4.000000e+02
5 6 8 1 1.110000e+02 4.000000e+02 1.110000e+02
4.000000e+02
5 6 9 1 1.120000e+02 4.000000e+02 1.120000e+02
4.000000e+02
7 6 8 1 1.080000e+02 4.000000e+02 1.080000e+02
4.000000e+02
7 6 9 1 1.070000e+02 4.000000e+02 1.070000e+02
4.000000e+02
8 6 9 1 1.070000e+02 4.000000e+02 1.070000e+02
4.000000e+02
[ dihedrals ]
; ai aj ak al funct c0 c1
c2 c3 c4 c5
18 1 2 3 3 3.600000e+02 5.000000e+00 3.000000e+00
3.600000e+02 5.000000e+00 3.000000e+00
1 2 3 4 3 3.600000e+02 5.000000e+00 3.000000e+00
3.600000e+02 5.000000e+00 3.000000e+00
2 3 4 5 3 3.600000e+02 5.000000e+00 3.000000e+00
3.600000e+02 5.000000e+00 3.000000e+00
3 4 5 6 3 3.600000e+02 5.000000e+00 3.000000e+00
3.600000e+02 5.000000e+00 3.000000e+00
4 5 6 7 3 3.600000e+02 5.000000e+00 3.000000e+00
3.600000e+02 5.000000e+00 3.000000e+00
[ system ]
; Name
HEX
[ molecules ]
; Compound #mols
HEX 256
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
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 [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
