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

[Justin A. Lemkul]

Pablo Englebienne wrote:
> 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:
> 

If bonds are supposed to be defined between, i.e. H-Cl and Cl-Cl, then this may 
not be appropriate.  Instead, you could perhaps set the covalent bonds as 
[constraints] in the topology, leaving the other "nonbonded" bonds as [bonds] 
(and then use "constraints = none" in the .mdp file, so only the [constraints] 
are used).

> --[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)
> 

What is the purpose of position restraints here?

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

Probably a consequence of restraining the starting structure.  I'd say you're 
not equilibrating appropriately with position restraints imposed.

> 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).
> 

Insufficient equilibration, and a variety of factors within the .mdp file.  If 
you want feedback on run parameters, you'll have to post the .mdp.

-Justin

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

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

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

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