[gmx-users] How to reduce high repulsion from system?
Justin A. Lemkul
jalemkul at vt.edu
Fri Apr 30 01:14:00 CEST 2010
Justin A. Lemkul wrote:
>
> 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!
>
I should also add in here that x2top does produce the correct output if one
employs the -alldih flag (which is off by default).
> 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