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
_______________________________________________
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
