[gmx-users] Re: Chloroform (CHCl3) solvent box for G53a5 force field

[Pablo Englebienne]

Following up on the previous message, I noticed that the topology I 
previously sent (including 10 bonds) works for a minimization, but not 
for an MD simulation. grompp issues the following warning:

WARNING 1 [file topol.top, line 29]:
  Molecule type 'CHCL3' has 10 constraints.
  For stability and efficiency there should not be more constraints than
  internal number of degrees of freedom: 9.

I therefore used the following .itp file, with only the C-H and C-Cl bonds:

--[chcl3.itp]--
[ moleculetype ]
; Name            nrexcl
CHCL3             1

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  
typeB    chargeB      massB
     1       HCHL      1  CHCL3   HChL      1      0.082      1.008
     2       CCHL      1  CHCL3   CChL      1      0.179     12.011
     3      CLCHL      1  CHCL3  CLCh1      1     -0.087     35.453
     4      CLCHL      1  CHCL3  CLCh2      1     -0.087     35.453
     5      CLCHL      1  CHCL3  CLCh3      1     -0.087     35.453

[ bonds ]
;  ai    aj funct            c0            c1            c2            c3
    1     2     2    gb_39
    2     3     2    gb_40
    2     4     2    gb_40
    2     5     2    gb_40

[ angles ]
;  ai    aj    ak funct            c0            c1            
c2            c3
    1     2     3     2  ga_43
    1     2     4     2  ga_43
    1     2     5     2  ga_43
    3     2     4     2  ga_44
    3     2     5     2  ga_44
    4     2     5     2  ga_44

[ exclusions ]
1 2 3 4 5
2 1 3 4 5
3 1 2 4 5
4 1 2 3 5
5 1 2 3 4
--[chcl3.itp]--

These are the steps I took:
- build a box with 216 CHCL3 molecules with genbox
- adjusted the density to 1479 with editconf
- minimized the box to F<100 (steep, 50000 steps)
- equilibrated NVT (position restrained C), 100ps, 300K, tau_t 0.1
- equilibrated NPT (position restrained C), 100ps, 300K, 1bar, tau_t 
0.1, tau_p 2.0, compressibility 1e-4 (from CRC handbook)

Until here, everything looks decent, except for relatively large 
fluctuations in T (RMSD ~9K) and density (RMSD ~10kg m-3).

I then performed an unconstrained MD, 1ns, otherwise identical 
parameters to NPT equilibration. Temp (300K) and density (1450 kg m-3) 
stable, but fluctuating (RMSD 8 and 20 respectively).

In order to compare the results with Tironi and van Gunsteren, Molecular 
Physics 1994, 83, 381, who used the same GROMOS parameters: Epot = -28.6 
+/- 0.3 kJ/mol, density 1520 +/- 12 kg m-3.

This is the output I get from g_energy:

Energy                      Average       RMSD     Fluct.      Drift  
Tot-Drift
-------------------------------------------------------------------------------
Potential                  -4720.16    119.141    118.978 -0.0215895   
-21.5896

Therefore molar Epot = -(-4720)/216 = 21.9 kJ/mol. What factors could be 
accountable for the decrease wrt the reported value of 28.4? 
Experimental deltaHv is 31.4 kJ/mol (same reference).

Would there be any other parameters I should check before using this 
solvent box in production runs?

Regards,
Pablo

-- 
Pablo Englebienne, PhD
Institute of Complex Molecular Systems (ICMS)
Eindhoven University of Technology, TU/e
PO Box 513, HG -1.26
5600 MB Eindhoven, The Netherlands
Tel +31 40 247 5349