[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