[gmx-users] Wrong bond and Lennard-Jones potentials in cyclohexane
Mark Abraham
Mark.Abraham at anu.edu.au
Fri Dec 11 21:16:44 CET 2009
Henri Ervasti wrote:
> Dear all,
See comment below.
> I have tried to simulate cyclohexane to produce vaporization enthalpy
> for it. I have tried it with two approaches:
>
> 1) using geometry, torsion potentials (harmonic) + charges obtained from
> B3LYP/6-311++G** (DFT) calculations and OPLS-AA Lennard-Jones, bond and
> angle potentials.
>
> 2) Using only OPLS-AA force field values (slightly modified for
> torsions).
>
> In both cases I do achieve a reasonable accuracy for density, but the
> energetics go totally wrong and even the sign of the result is the
> opposite than it should be. I have noticed that this results from very
> large bond potentials and L-J potentials in the averaged results. From a
> 512 molecules NPT MD simulation (minimization using steep, 400ps
> relaxation, 2ns run, 2fs time-step) using OPLS-AA parameters only, I get
> an average energy (the total divided by 512) for a bond to be 115.9
> kJ/mol, while for a single cyclohexane in vacuum this is 24.4 kJ/mol.
> The L-J results into an average (the total divided by 512) 14.03 kJ/mol,
> which suggests that the L-J part would be repulsive! If I approximate
> that the intramolecular energies would be the same (thus ignore bond,
> angle, torsion terms) I get vaporization enthalpy of about -20.6 kJ/mol,
> while the experimental result is around +33 kJ/mol. Thus, a wrong sign
> indicating that the system should be in the gas phase, enlarging (which
> it does not do). The results are even slightly worse using the DFT
> geometries, torsions, charges. Also, the creators of OPLS-AA have
> simulated cyclohexane and obtained DvapH very close to the experimental
> results [Jorgensen et al, JACS 1996, 118, 11225-11236 + the supporting
> info].
>
> I have tried using different time-step, different L-J for hydrogen and
> carbon, different torsion potentials, but the results are always the
> same: very large potentials for bonds and L-J.
>
> This discrepancy baffles me, especially that this problem does not
> appear when I have modelled f.ex. methanol, ethanol, hexane, acetone...
> Even N-formyl morpholine did not show this problem, although it also has
> a six-membered ring. Could this be some kind of problem in the code for
> (more) symmetric (non-aromatic) ring/cyclic molecules?
>
> See the topology file and a pdb file for the OPLS-AA cyclohexane
> simulation at the end of the message.
>
> Thank you for your help!
>
> Kind regards,
> Henri Ervasti
>
> PDB:
>
> HEADER
> HETATM 1 CAA CHE 1 0.140 2.640 0.160 1.00 20.00
> C
> HETATM 2 CAB CHE 1 -0.030 1.160 -0.190 1.00 20.00
> C
> HETATM 3 CAC CHE 1 1.450 3.180 -0.410 1.00 20.00
> C
> HETATM 4 HAA CHE 1 0.148 2.753 1.244 1.00 20.00
> H
> HETATM 5 CAD CHE 1 1.160 0.360 0.350 1.00 20.00
> C
> HETATM 6 HAC CHE 1 -0.080 1.048 -1.273 1.00 20.00
> H
> HETATM 7 HAD CHE 1 -0.916 0.808 0.338 1.00 20.00
> H
> HETATM 8 CAF CHE 1 2.470 0.890 -0.220 1.00 20.00
> C
> HETATM 9 HAG CHE 1 1.045 -0.687 0.068 1.00 20.00
> H
> HETATM 10 HAH CHE 1 1.191 0.528 1.426 1.00 20.00
> H
> HETATM 11 CAE CHE 1 2.630 2.380 0.130 1.00 20.00
> C
> HETATM 12 HAK CHE 1 2.466 0.772 -1.304 1.00 20.00
> H
> HETATM 13 HAL CHE 1 3.277 0.358 0.285 1.00 20.00
> H
> HETATM 14 HAI CHE 1 2.595 2.455 1.217 1.00 20.00
> H
> HETATM 15 HAJ CHE 1 3.552 2.756 -0.313 1.00 20.00
> H
> HETATM 16 HAB CHE 1 -0.668 3.178 -0.337 1.00 20.00
> H
> HETATM 17 HAF CHE 1 1.420 3.015 -1.487 1.00 20.00
> H
> HETATM 18 HAE CHE 1 1.563 4.226 -0.123 1.00 20.00
> H
>
> TOPOLOGY:
>
> [ defaults ]
> ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> 1 2 yes 0.5 0.5
See manual 5.3.3. ffoplsaa.itp uses combination rule 3.
Mark
> [ atomtypes ]
> ; name at.num mass charge ptype sigma epsilon
> CH1 6 12.01100 0.000 A 3.50000e-01 2.76144e-01
> HC 1 1.00790 0.000 A 2.50000e-01 1.25000e-01
>
> [ nonbond_params ]
> ; i j func c6 c12
>
> [ pairtypes ]
> ; i j func cs6 cs12
>
> [ moleculetype ]
> ; Name nrexcl
> CHE 5
>
> [ atoms ]
> ; nr type resnr resid atom cgnr charge mass
> 1 CH1 1 CHE CAA 1 -0.120 12.0110
> 2 CH1 1 CHE CAB 2 -0.120 12.0110
> 3 CH1 1 CHE CAC 3 -0.120 12.0110
> 4 HC 1 CHE HAA 1 0.060 1.0079
> 5 CH1 1 CHE CAD 4 -0.120 12.0110
> 6 HC 1 CHE HAC 2 0.060 1.0079
> 7 HC 1 CHE HAD 2 0.060 1.0079
> 8 CH1 1 CHE CAF 5 -0.120 12.0110
> 9 HC 1 CHE HAG 4 0.060 1.0079
> 10 HC 1 CHE HAH 4 0.060 1.0079
> 11 CH1 1 CHE CAE 6 -0.120 12.0110
> 12 HC 1 CHE HAK 5 0.060 1.0079
> 13 HC 1 CHE HAL 5 0.060 1.0079
> 14 HC 1 CHE HAI 6 0.060 1.0079
> 15 HC 1 CHE HAJ 6 0.060 1.0079
> 16 HC 1 CHE HAB 1 0.060 1.0079
> 17 HC 1 CHE HAF 3 0.060 1.0079
> 18 HC 1 CHE HAE 3 0.060 1.0079
>
> [ bonds ]
> ; ai aj fu c0, c1, ... CHARMM/22 ; B3LYP/CBSB7
> 1 2 1 0.1529 224262.4 ; 0.154
> 1 3 1 0.1529 224262.4 ; 0.154
> 1 4 1 0.1090 284512.0 ; 0.110
> 1 16 1 0.1090 284512.0 ; 0.110
> 2 5 1 0.1529 224262.4 ; 0.154
> 2 6 1 0.1090 284512.0 ; 0.110
> 2 7 1 0.1090 284512.0 ; 0.110
> 3 11 1 0.1529 224262.4 ; 0.154
> 3 17 1 0.1090 284512.0 ; 0.110
> 3 18 1 0.1090 284512.0 ; 0.110
> 5 8 1 0.1529 224262.4 ; 0.154
> 5 9 1 0.1090 284512.0 ; 0.110
> 5 10 1 0.1090 284512.0 ; 0.110
> 8 11 1 0.1529 224262.4 ; 0.154
> 8 12 1 0.1090 284512.0 ; 0.110
> 8 13 1 0.1090 284512.0 ; 0.110
> 11 14 1 0.1090 284512.0 ; 0.110
> 11 15 1 0.1090 284512.0 ; 0.110
>
> [ pairs ]
> ; ai aj fu c0, c1, ...
>
> [ exclusions ]
> ; ai aj ak ...
>
> [ angles ]
> ; ai aj ak fu c0, c1, ... CHARMM/22 ; B3LYP/CBSB7
> 1 2 6 1 110.7 313.800 ; 110.2
> 1 2 7 1 110.7 313.800 ; 109.1
> 1 3 11 1 112.7 488.273 ; 111.6
> 1 3 17 1 110.7 313.800 ; 110.2
> 1 3 18 1 110.7 313.800 ; 109.1
> 2 1 3 1 112.7 488.273 ; 111.5
> 2 1 4 1 110.7 313.800 ; 109.1
> 2 1 16 1 110.7 313.800 ; 110.2
> 2 5 8 1 112.7 488.273 ; 111.6
> 2 5 9 1 110.7 313.800 ; 109.1
> 2 5 10 1 110.7 313.800 ; 110.2
> 3 1 4 1 110.7 313.800 ; 109.1
> 3 1 16 1 110.7 313.800 ; 110.2
> 3 11 8 1 112.7 488.273 ; 111.6
> 3 11 14 1 110.7 313.800 ; 109.1
> 3 11 15 1 110.7 313.800 ; 110.2
> 4 1 16 1 107.8 276.144 ; 106.4
> 5 2 6 1 110.7 313.800 ; 110.2
> 5 2 7 1 110.7 313.800 ; 109.1
> 5 8 11 1 112.7 488.273 ; 111.5
> 5 8 12 1 110.7 313.800 ; 110.2
> 5 8 13 1 110.7 313.800 ; 109.2
> 6 2 7 1 107.8 276.144 ; 106.5
> 8 5 9 1 110.7 313.800 ; 109.1
> 8 5 10 1 110.7 313.800 ; 110.2
> 8 11 14 1 110.7 313.800 ; 109.1
> 8 11 15 1 110.7 313.800 ; 110.2
> 9 5 10 1 107.8 276.144 ; 106.5
> 11 3 17 1 110.7 313.800 ; 110.2
> 11 3 18 1 110.7 313.800 ; 109.1
> 11 8 12 1 110.7 313.800 ; 110.2
> 11 8 13 1 110.7 313.800 ; 109.2
> 12 8 13 1 107.8 276.144 ; 106.4
> 14 11 15 1 107.8 276.144 ; 106.5
>
> [ dihedrals ]
> ; ai aj ak al fu c0, c1, m, ...
> 1 2 5 8 5 0.0000 0.0000 40.0000 0.0000
> 1 2 5 9 5 0.0000 0.0000 1.2552 0.0000
> 1 2 5 10 5 0.0000 0.0000 1.2552 0.0000
> 1 3 11 8 5 0.0000 0.0000 40.0000 0.0000
> 1 3 11 14 5 0.0000 0.0000 1.2552 0.0000
> 1 3 11 15 5 0.0000 0.0000 1.2552 0.0000
> 2 1 3 11 5 0.0000 0.0000 40.0000 0.0000
> 2 1 3 17 5 0.0000 0.0000 1.2552 0.0000
> 2 1 3 18 5 0.0000 0.0000 1.2552 0.0000
> 2 5 8 11 5 0.0000 0.0000 40.0000 0.0000
> 2 5 8 12 5 0.0000 0.0000 1.2552 0.0000
> 2 5 8 13 5 0.0000 0.0000 1.2552 0.0000
> 3 1 2 5 5 0.0000 0.0000 40.0000 0.0000
> 3 1 2 6 5 0.0000 0.0000 1.2552 0.0000
> 3 1 2 7 5 0.0000 0.0000 1.2552 0.0000
> 3 11 8 5 5 0.0000 0.0000 40.0000 0.0000
> 3 11 8 12 5 0.0000 0.0000 1.2552 0.0000
> 3 11 8 13 5 0.0000 0.0000 1.2552 0.0000
> 4 1 2 5 5 0.0000 0.0000 1.2552 0.0000
> 4 1 2 6 5 0.0000 0.0000 1.2552 0.0000
> 4 1 2 7 5 0.0000 0.0000 1.2552 0.0000
> 4 1 3 11 5 0.0000 0.0000 1.2552 0.0000
> 4 1 3 17 5 0.0000 0.0000 1.2552 0.0000
> 4 1 3 18 5 0.0000 0.0000 1.2552 0.0000
> 5 2 1 16 5 0.0000 0.0000 1.2552 0.0000
> 5 8 11 14 5 0.0000 0.0000 1.2552 0.0000
> 5 8 11 15 5 0.0000 0.0000 1.2552 0.0000
> 6 2 1 16 5 0.0000 0.0000 1.2552 0.0000
> 6 2 5 8 5 0.0000 0.0000 1.2552 0.0000
> 6 2 5 9 5 0.0000 0.0000 1.2552 0.0000
> 6 2 5 10 5 0.0000 0.0000 1.2552 0.0000
> 7 2 1 16 5 0.0000 0.0000 1.2552 0.0000
> 7 2 5 8 5 0.0000 0.0000 1.2552 0.0000
> 7 2 5 9 5 0.0000 0.0000 1.2552 0.0000
> 7 2 5 10 5 0.0000 0.0000 1.2552 0.0000
> 8 11 3 17 5 0.0000 0.0000 1.2552 0.0000
> 8 11 3 18 5 0.0000 0.0000 1.2552 0.0000
> 9 5 8 11 5 0.0000 0.0000 1.2552 0.0000
> 9 5 8 12 5 0.0000 0.0000 1.2552 0.0000
> 9 5 8 13 5 0.0000 0.0000 1.2552 0.0000
> 10 5 8 11 5 0.0000 0.0000 1.2552 0.0000
> 10 5 8 12 5 0.0000 0.0000 1.2552 0.0000
> 10 5 8 13 5 0.0000 0.0000 1.2552 0.0000
> 11 3 1 16 5 0.0000 0.0000 1.2552 0.0000
> 12 8 11 14 5 0.0000 0.0000 1.2552 0.0000
> 12 8 11 15 5 0.0000 0.0000 1.2552 0.0000
> 13 8 11 14 5 0.0000 0.0000 1.2552 0.0000
> 13 8 11 15 5 0.0000 0.0000 1.2552 0.0000
> 14 11 3 17 5 0.0000 0.0000 1.2552 0.0000
> 14 11 3 18 5 0.0000 0.0000 1.2552 0.0000
> 15 11 3 17 5 0.0000 0.0000 1.2552 0.0000
> 15 11 3 18 5 0.0000 0.0000 1.2552 0.0000
>
> [ system ]
> cyclohexane
>
> [ molecules ]
> CHE 512
>
>
>
More information about the gromacs.org_gmx-users
mailing list