[gmx-users] How to reduce high repulsion from system?
Justin A. Lemkul
jalemkul at vt.edu
Fri Apr 30 00:49:19 CEST 2010
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
========================================
More information about the gromacs.org_gmx-users
mailing list