[gmx-users] A simulation of pure chloroform with AMBER compatible parameters never gives the right density

[anaome at fundp.ac.be]

Dear Users,

I am trying to reproduce the results obtained by Fox & Kollman on  
chloroform presented in the following papers but can never achieve the  
proper density of 1.459 they obtain with their own rigid CHCL3 model.  
The best I obtain with the parameters presented below is 1.390. I  
tried with both version 4.5.3 and 4.0.7.
Can anyone help me on that ?
Thank you

Aymeric Naômé

The Journal of Physical Chemistry, 1996, 100, 25, 10779-10783
The Journal of Physical Chemistry B, 1998, 102, 41, 8070-8079

The parameters are:
qC=-0.345 eps=0.1094 kcal/mol R*=1.908 Ang.
qCl=0.012 eps=0.2550 R*=2.000
qH=0.309  eps=0.0157 R*=1.187
C-Cl=1.758 Ang.
C-H=1.100
H-C-Cl=107.7°
Cl-C-Cl=112.2°

And my corresponding topology file is the following:
sigma=2*(R*/10)/(2^(1/6))
epsilon=eps*4.184

++++++++++++++++++++
[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               no              0.5     0.8333

[ atomtypes ]
; name      at.num  mass     charge ptype  sigma      epsilon
CT           6      12.01    0.0000  A   3.39967e-01  4.57730e-01
H3           1       1.008   0.0000  A   2.11499e-01  6.56888e-02
Cl          17      35.453   0.0000  A   3.56359e-01  1.10669e+00


[ moleculetype ]
; Name            nrexcl
CHCL3             3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass   
typeB    chargeB      massB
      1         CT      1  CHCL3      C      1     -0.345     12.011    
; qtot -0.386
      2         H3      1  CHCL3      H      1      0.309      1.008    
; qtot -0.120
      3         Cl      1  CHCL3    Cl1      1      0.012     35.453    
; qtot -0.080
      4         Cl      1  CHCL3    Cl2      1      0.012     35.453    
; qtot -0.040
      5         Cl      1  CHCL3    Cl3      1      0.012     35.453    
; qtot  0
[ constraints ]
     1     3     1    0.1758
     1     4     1    0.1758
     1     5     1    0.1758
     1     2     1    0.1100
;    2     3     2    0.2340
     2     4     2    0.2340
     2     5     2    0.2340
     3     4     2    0.2901
     4     5     2    0.2901
     5     3     2    0.2901

[ system ]
; Name
CHCl3_256

[ molecules ]
; Compound        #mols
CHCL3             256
++++++++++++++++++++

The simulation parameters in the first paper:
MD program: AMBER 4.1
200-300 molecules
dt=2fs (rigid model)
cut-off=12 Ang.
long range correction according to Allen&Tildesley
Berendsen thermostat with coupling constant 0.4 ps-1
temperature=300 K
inverse compressibility 1.086 10-4 bar-1
Berendsen barostat with coupling constant 0.2 ps-1
pressure=1 atm
SHAKE constraint solver

Below is my mdp file with the corresponding parameters:
++++++++++++++++++++
title		    =  CHCl3_box
cpp                 =  /lib/cpp -traditional
define              =
constraints         =
integrator          =  md
dt                  =  0.002	; fs !
nsteps              =  400000	; total 10000 ps.
nstcomm             =  1
nstxout             =  500
nstvout             =  1000
nstfout             =  0
nstlog              =  500
nstenergy           =  500
nstlist             =  10
ns_type             =  grid
rlist               =  1.2
coulombtype		=  cut-off
rcoulomb            =  1.2
rvdw                =  1.2
DispCorr            =  EnerPres
fourierspacing		=  0.12
fourier_nx		=  0
fourier_ny		=  0
fourier_nz		=  0
pme_order		=  6
ewald_rtol		=  1e-5
optimize_fft		=  yes
energygrps		=  CHCL3
; Berendsen temperature coupling is on in three groups
Tcoupl              =  berendsen
tau_t               = 5
tc_grps		    = CHCL3
ref_t               = 300
; Pressure coupling is on
Pcoupl              =  Berendsen
pcoupltype          =  isotropic
tau_p               =  2.5
compressibility     =  10.86e-5
ref_p               =  1.0
; Generate velocites is on at 300 K.
gen_vel             =  yes
gen_temp            =  300.0
gen_seed            =  173529
; Energy group exclusion
energygrp_excl      =
freezegrps          =
freezedim           =
; Non-equilibrium Thermodynamics
acc_grps	    =
accelarate	    =
; center of mass
comm_mode           = linear
comm_grps           =
; Type of constraint algorithm
constraint-algorithm     = Shake
++++++++++++++++++++

-- 
Aymeric Naômé
Ph. D. Student
UCPTS Division
University of Namur
61 Rue de Bruxelles
5000 Namur
BELGIUM