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

[Justin A. Lemkul]

Cenfeng Fu wrote:
> Hello all,
>  I am simulating liquid benzene. Initially, I generate a cubic box with 
> 600 benzene molecules in it. The length of the box is 4.6 nm. Then I do 
> the energy minimization and a NPT simulation for 20 ns. The force field 
> I used for benzene is OPLSAA (JACS,1990,112,4768). And here is my .mdp 
> file for the NPT simulation.
> 
> title               =  benz
> cpp                 =  /usr/bin/cpp
> define              =
> integrator          =  md
> dt                  =  0.002    ; ps !
> nsteps              =  10000000
> comm_mode           =  Linear
> comm_grps           =  System
> nstxout             =  0
> nstvout             =  2500
> nstlog              =  500
> nstenergy           =  50
> nstxtcout           =  500
> nstlist             =  10
> ns_type             =  grid
> pbc                 =  xyz
> rlist               =  1.0
> coulombtype         =  PME
> rcoulomb            =  1.0
> vdwtype             =  cut-off
> rvdw                =  1.0
> DispCorr            =  EnerPres
> fourierspacing          =  0.12
> fourier_nx              =  0
> fourier_ny              =  0
> fourier_nz              =  0
> pme_order               =  4
> ewald_rtol              =  1e-5
> optimize_fft            =  yes
> tcoupl                  =  nose-hoover
> ref_t                   =  300
> tau_t                   =  0.5
> tc_grps                 =  System
> pcoupl                  =  Parrinello-Rahman
> pcoupltype              =  isotropic
> tau_p                   =  3.0
> compressibility         =  4.5e-5
> ref_p                   =  1.0
> gen_vel                 =  yes
> gen_temp                =  300
> gen_seed                =  175298
> constraints             =  hbonds
> constraint_algorithm    =  LINCS
> ;shake_tol               =  0.0001
> lincs_order             =  4
> lincs_iter              =  1
> 
> To my surprise, the density of the system is too large to the 
> experimental value (0.874 g/cm^3). Here are some results from the 
> g_energy command.
> 
> Statistics over 7500001 steps [ 5000.0000 through 20000.0000 ps ], 8 
> data sets
> All statistics are over 750001 points
> 
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
> -------------------------------------------------------------------------------
> Potential                   13996.9         15    302.701   -83.3695 
>  (kJ/mol)
> Temperature                 299.976    0.00021    3.38239 -0.000582958  (K)
> Pressure                    1.08499     0.0072    1812.66  -0.026224  (bar)
> Box-X                       4.33432    0.00086  0.0102765 -0.00488475  (nm)
> Box-Y                       4.33432    0.00086  0.0102765 -0.00488475  (nm)
> Box-Z                       4.33432    0.00086  0.0102765 -0.00488475  (nm)
> Volume                      81.4276      0.049   0.579125  -0.275291  (nm^3)
> Density                     955.828       0.57    6.80024    3.23201 
>  (kg/m^3)
> 
> Why the density is so large ? Could anyone give me some suggestions?
> 

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

> Thank you in advance!
> 
> Sincerely
> 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

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