[gmx-users] Re: the density of liquid benzene is too large

[Justin A. Lemkul]

Cenfeng Fu wrote:
> 
>     Several possible reasons:
> 
>     1. The force field parameters aren't perfect, so there is some inherent
>     disagreement between simulation and reality.  What is the expected
>     value for
>     this force field? 
> 
> 
>     2. You're using the isothermal compressibility for water.  If your
>     system is
>     pure liquid benzene, I'd think you would want to use the
>     compressibility for
>     benzene.  I don't know how big the difference would be off-hand, but
>     at least
>     you'd be convinced that your simulation was set up properly.
> 
>     3. Although it won't matter a huge amount, to what temperature does the
>     experimental density correspond?  Usually these parameters are given
>     at 25C (298
>     K).  The temperature you've used is 300 K.  Again, a minor point,
>     but one worth
>     doing correctly in a simulation even though the such a change would
>     not account
>     for the magnitude of difference you're currently seeing.
> 
>     -Justin
> 
> Hi Justin,
>   Thanks for your suggestions!
>   The OPLSAA model for liquid benzene could get a density of 
> 0.873+/-0.001 g/cm^3 (JACS,1990,112,4768) at 298K and 1 atm. I have done 
> some test and maybe I have found what is the problem. In the previous 
> simulation, I applied long range dispersion corrections for energy and 
> pressure with "DispCorr = EnerPres". When I apply the long range 
> dispersion corrections only for energy with "DispCorr = Ener" (the 
> temperature is set to be 298K and the pressure is set to be 1 bar), the 
> density of the system is 0.883, which is closely to the experimental 
> value and the expect value of this model.
>   Now, I have anther question. After the simulation, I want to calculate 
> the hear capacity of liquid benzene. So I using this command:
> g_energy -f *.edr -s *.tpr -o energy.xvg -b 10000 -nmol 600 -nconstr 12
>  And these are the results:
> Statistics over 5000001 steps [ 10000.0000 through 20000.0000 ps ], 11 
> data sets
> All statistics are over 500001 points
> 
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
> -------------------------------------------------------------------------------
> Potential                   19.7926      0.036   0.435187 -0.0870648  
> (kJ/mol)
> Kinetic En.                 29.7264    1.7e-05   0.360236 -7.68964e-05  
> (kJ/mol)
> Total Energy                 49.519      0.036   0.575006 -0.0871418  
> (kJ/mol)
> Temperature                 297.999    0.00017    3.61127 -0.000770665  (K)
> Pressure                    1.15978     0.0036    171.048 -0.00963121  (bar)
> Box-X                       4.45033     0.0015  0.0105071 -0.0038387  (nm)
> Box-Y                       4.45033     0.0015  0.0105071 -0.0038387  (nm)
> Box-Z                       4.45033     0.0015  0.0105071 -0.0038387  (nm)
> Volume                      88.1423      0.089   0.624413  -0.228234  (nm^3)
> Density                     883.013       0.89    6.25236    2.28197  
> (kg/m^3)
> Enthalpy                    29712.2         21    345.005   -52.2857  
> (kJ/mol)
> 
> Temperature dependent fluctuation properties at T = 297.999. #constr/mol 
> = 12
> Isothermal Compressibility: 0.000107512 /bar
> Adiabatic bulk modulus:        9301.25  bar
> Heat capacity at constant pressure Cp:    218.791 J/mol K
> Thermal expansion coefficient alphaP: 0.000135136 1/K
> 
> I got a hear capacity at constant pressure with 218.791 J/(mol K). 
> However, the experimental value is 135.98 J/(mol K), and the expect 
> value of the model is 130.54 J/(mol K) (JACS,1990,112,4768). In the new 
> simulation, I used LINCS for all-bonds. So I think the #nconstr should 
> be 12 in the g_energy command. Is this value for #nconstr right? If it 
> is wrong, what value should I use. Or I should not use long range 
> dispersion correction for energy? Are there other mistakes with my 
> parameters?
> 

The original derivation of OPLS benzene parameters (in the paper you cite) 
applied dispersion correction to the energy term, so I would suspect that you 
should do the same.  No constraints are mentioned in that paper, so it is hard 
to say whether or not you should be using them.

You can easily see the effects of constraints by running a few more simulations 
(no constraints, then hbonds only, etc).

Another challenge pertains to the use of cutoffs.  The paper you cite used a 
1.3-nm ring-ring cutoff.  I don't know how you would implement that in Gromacs, 
but in your original post 
(http://lists.gromacs.org/pipermail/gmx-users/2010-October/054499.html) you were 
using 1.0-nm cutoffs.  Differences in the way nonbonded interactions are 
calculated can (and quite often, will) make a difference in your results.

-Justin

> Best regards!
> Cenfeng Fu
> 

-- 
========================================

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

========================================