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
========================================
_______________________________________________
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
