Dear all,
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
[ 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
--
-----------------------------------
Dr. Henri K. Ervasti
Delft University of Technology
Process & Energy Laboratory
Leeghwaterstraat 44 k. 0.30
2628CA Delft
Netherlands
e-mail: h.k.erva...@tudelft.nl
phone: +31 (0)15 27 86662
fax: +31 (0)15 27 82460
--
gmx-users mailing list gmx-users@gromacs.org
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 gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
