[gmx-users] Wrong bond and Lennard-Jones potentials in cyclohexane
Carl Caleman
calle at xray.bmc.uu.se
Fri Dec 11 21:44:00 CET 2009
Hi,
I have simulated density and for a couple of cyclohex*** molecules,
using both opls and gaff parameters, and that worked pretty ok for
example (DHvap values):
molecule: gaff: opls: exp:
cyclohexanone 51.33 50.45 45.06 (298.15 K)
cyclohexylamine 59.07 49.84 43.67 (298.15 K)
However I have not worked with cyclohexane. I dont know if this helps
you. I put in the topol.top (opls) file I have used for cyclohexanone,
I hope that helps you. We have generated the opls parameters using
x2top, and then manually made sure that the topologies made sense.
Good luck,
/Carl
#include "ffoplsaa.itp"
[ moleculetype ]
; Name nrexcl
ICE 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
typeB chargeB massB
1 opls_136 1 LIG C 1 -0.12 12.011 ;
2 opls_140 1 LIG H 1 0.06 1.008 ;
3 opls_140 1 LIG H 1 0.06 1.008 ;
4 opls_136 1 LIG C 2 -0.12 12.011 ;
5 opls_140 1 LIG H 2 0.06 1.008 ;
6 opls_140 1 LIG H 2 0.06 1.008 ;
7 opls_136 1 LIG C 3 -0.12 12.011 ;
8 opls_140 1 LIG H 3 0.06 1.008 ;
9 opls_140 1 LIG H 3 0.06 1.008 ;
10 opls_136 1 LIG C 4 -0.12 12.011 ;
11 opls_140 1 LIG H 4 0.06 1.008 ;
12 opls_140 1 LIG H 4 0.06 1.008 ;
13 opls_136 1 LIG C 5 -0.12 12.011 ;
14 opls_140 1 LIG H 5 0.06 1.008 ;
15 opls_140 1 LIG H 5 0.06 1.008 ;
16 opls_280 1 LIG C 6 0.47 12.011 ;
17 opls_281 1 LIG O 6 -0.47 15.9994 ;
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 16 1
4 5 1
4 6 1
4 7 1
7 8 1
7 9 1
7 10 1
10 11 1
10 12 1
10 13 1
13 14 1
13 15 1
13 16 1
16 17 1
[ pairs ]
; ai aj funct c0 c1 c2 c3
1 8 1
1 9 1
1 10 1
1 14 1
1 15 1
2 5 1
2 6 1
2 7 1
2 13 1
2 17 1
3 5 1
3 6 1
3 7 1
3 13 1
3 17 1
4 11 1
4 12 1
4 13 1
4 17 1
5 8 1
5 9 1
5 10 1
5 16 1
6 8 1
6 9 1
6 10 1
6 16 1
7 14 1
7 15 1
7 16 1
8 11 1
8 12 1
8 13 1
9 11 1
9 12 1
9 13 1
10 17 1
11 14 1
11 15 1
11 16 1
12 14 1
12 15 1
12 16 1
14 17 1
15 17 1
[ angles ]
; ai aj ak funct c0 c1 c2
c3
2 1 3 1
2 1 4 1
2 1 16 1
3 1 4 1
3 1 16 1
4 1 16 1
1 4 5 1
1 4 6 1
1 4 7 1
5 4 6 1
5 4 7 1
6 4 7 1
4 7 8 1
4 7 9 1
4 7 10 1
8 7 9 1
8 7 10 1
9 7 10 1
7 10 11 1
7 10 12 1
7 10 13 1
11 10 12 1
11 10 13 1
12 10 13 1
10 13 14 1
10 13 15 1
10 13 16 1
14 13 15 1
14 13 16 1
15 13 16 1
1 16 13 1
1 16 17 1
13 16 17 1
[ dihedrals ]
; ai aj ak al funct c0 c1 c2
c3 c4 c5
2 1 4 5 3
2 1 4 6 3
2 1 4 7 3
3 1 4 5 3
3 1 4 6 3
3 1 4 7 3
16 1 4 5 3
16 1 4 6 3
16 1 4 7 3
2 1 16 13 3
2 1 16 17 3
3 1 16 13 3
3 1 16 17 3
4 1 16 13 3
4 1 16 17 3
1 4 7 8 3
1 4 7 9 3
1 4 7 10 3
5 4 7 8 3
5 4 7 9 3
5 4 7 10 3
6 4 7 8 3
6 4 7 9 3
6 4 7 10 3
4 7 10 11 3
4 7 10 12 3
4 7 10 13 3
8 7 10 11 3
8 7 10 12 3
8 7 10 13 3
9 7 10 11 3
9 7 10 12 3
9 7 10 13 3
7 10 13 14 3
7 10 13 15 3
7 10 13 16 3
11 10 13 14 3
11 10 13 15 3
11 10 13 16 3
12 10 13 14 3
12 10 13 15 3
12 10 13 16 3
10 13 16 1 3
10 13 16 17 3
14 13 16 1 3
14 13 16 17 3
15 13 16 1 3
15 13 16 17 3
[ system ]
; Name
ICE
[ molecules ]
; Compound #mols
ICE 216
Henri Ervasti wrote:
> 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
>
>
>
--
_______________________________________________
Carl Caleman, Ph.D.
Physik Department E17
Technische Universität München
email: carl.caleman at tum.de
More information about the gromacs.org_gmx-users
mailing list