Nilesh Dhumal wrote:
I have general question. If this parameters are not giving the proper
results then How come I get the proper results with same parameter if I
use 1/2 value of Kbond with the same parameters.
Random chance, I imagine. Your .mdp settings are not correct and the cutoffs
are unreasonably short for any modern force field. If you're seeking to
reproduce others' work, then you have to use the same protocol. You're not
doing that.
I checked the temp when I use half value of Kbond and it goes max. ~320.
This would also indicate instability. You managed to get a "correct" result for
one observable from an incorrect simulation wherein the thermodynamic
observables do not correspond to what you have set.
The way I found all of this out was first and foremost by watching the
trajectory. After about 30 ps, the water molecules freeze in place and vibrate
as if they are about to shear apart. At this exact time, the value of epsilon
(from epsilon.xvg, output by g_dipoles) plummets down, after having steadily
increased for the first 30 ps, as expected. The temperature spikes and the
pressure begins to oscillate by +/- 10^5, far higher than what it should for a
system of a few hundred water molecules.
Watch your trajectories, make sure you're getting the right average temperature
and pressure, and that all your other energetic terms stabilize, then analyze
the result of the simulation. Otherwise, any inferred output is meaningless.
-Justin
Nilesh
On Wed, May 11, 2011 4:03 pm, Justin A. Lemkul wrote:
Nilesh Dhumal wrote:
Hello Justin,
I used 0.5 fs time step and still I got dielectric constant ~2.
This is the md.mdp file I used. I checked the temperature it didn't go
more than ~ 320 K.
Still, that's unacceptable if your ref_t is 295. The fact that the
temperature is still rocketing up suggests the same instability I found.
title = cpeptide MD cpp = /usr/bin/cpp
integrator = md dt = 0.001 ; ps ! nsteps
= 6500000 ; total 5 ps.
nstcomm = 1000 nstxout = 1000 nstvout
= 1000
nstfout = 1 nstlist = 1 ns_type =
grid rlist = 0.6 rcoulomb = 0.6 rvdw
= 0.7
These cutoffs make no sense. The paper you've cited used 0.9 nm. If
you're trying to reproduce previous work, use the same settings.
coulombtype = PME vdwtype = cut-off
Similarly here, you haven't applied dispersion correction, which you
should, per the methods used in the JCP article.
-Justin
pbc = xyz fourierspacing = 0.12 fourier_nx = 0
fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol
= 1e-5
optimize_fft = yes ; Berendsen temperature coupling is on
Tcoupl = v-rescale
tau_t = 0.1 tc-grps =system ref_t = 295 ; Pressure coupling is on
Pcoupl = parrinello-rahman
pcoupltype = isotropic tau_p = 0.5 compressibility
= 4.5e-5
ref_p = 1.0 ; Generate velocites is on at 300 K.
gen_vel = yes gen_temp = 295.0 gen_seed
= 173529
Nilesh
On Tue, May 10, 2011 6:14 pm, Justin A. Lemkul wrote:
The gmx-developers list is not the forum for these types of
questions. I am replying via gmx-users, which is where the discussion
should stay.
I took a few minutes to dig into this. My conclusion is that your
system is not stable. I would encourage you to analyze the
temperature and pressure of your systems that are giving odd results.
When I used the
topology you provided, the temperature spiked to over 600 K and the
box began to oscillate extensively as if the system were about to
explode. The resulting epsilon value was about 1.6.
If I reduce the timestep to 0.5 fs, the results are much more
reasonable, and using the Ka and Kb values in the paper (not halved!)
I get a much
more sensible result, even after just 100 ps. It looks to me as if
the combination of these parameters and a 1-fs timestep is not
entirely stable. I know the original authors used dt = 1 fs, but the
point remains.
-Justin
Nilesh Dhumal wrote:
Hello,
I am using flexible water model for my system. I am referring a
paper J.
Chem. Phys. 124, 024503 (2006). Author have used Amber type force
field.
i.e. 1/2 Kbond (r-req)**2 + 1/2 Kangle(theta-thetaeq)**2.
Kbond = 443153.3808 kJ/mol nm**2
Kangle= 317.5656 kJ/mol rad**2.
I am using olss-aa force field parameters in Gromacs VERSION
4.0.7.
I checked some papers in which author have used oplsaa force field
in Gromacs. 1/2 factor is not in opls force field if I compare opls
and amber.
I didn’t get the proper dielectric constant for water when I used
the parameters reported in paper
(Kbond = 443153.3808 kJ/mol nm**2,Kangle= 317.5656 kJ/mol rad**2)
I half the value of Kb and get the proper dielectric constant (~80)
for water reported in paper. If I half Kangle then I don’t get
proper value. Below are the results for the dielectric constant of
water. Bond length is nm.
Here I have done some analysis. The original value reported in J.
Chem.
Phys. 124, 024503 2006, paper are
Kbond = 443153.3808 kJ/mol nm**2
Kangle = 317.5656 kJ/mol rad**2.
bond length Kbond angle Kangle dielectric constant
0.1012
443153.3808 113.24 317.5656 ~1.9 : orginal value
0.1012 221576.6904 113.24 317.5656 ~80 : 1/2
(Kbond)
0.1012 443153.3808 113.24 158.7828 ~1.58 : 1/2
(kangle)
0.1012 221576.6904 113.24 317.5656 ~1.9 : 1/2
(Kbond)&(Kangle)
Here I pasted spc_fw.itp file used for the simulation.
; Derived from parsing of runfiles/alat.top.orig
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
;1 3 yes 0.5 0.5
; comb-rule 3 is square-root sigma, the OPLSAA version
[ atomtypes ]
; full atom descriptions are available in ffoplsaa.atp
; name bond_type mass charge ptype sigma
epsilon opls_1001 OW 8 15.99940 -0.820 A
3.1655e-01
6.503e-01
opls_1002 HW 1 1.00800 0.410 A 0.00e+00
0.00e+00
[ bondtypes ]
; i j func b0 kb
OW HW 1 0.1012 443153.3808 ; J. Chem. Phys.
(2006),124,024503
[ angletypes ]
; i j k func th0 cth
HW OW HW 1 113.24 317.5656 ; J. Chem. Phys.
(2006),124,024503
[ moleculetype ]
; Name nrexcl
WAT 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
typeB chargeB massB 1 opls_1001 1 WAT OW 1
-0.82
15.99940 ;
2 opls_1002 1 WAT HW1 1 0.41 1.008 ;
3 opls_1002 1 WAT HW2 1 0.41 1.008 ;
[ bonds ]
; i j funct
1 2 1
1 3 1
[ angles ]
; i j k funct
2 1 3 1
In water.top file, I included spc_fw.itp file.
; Include water topology
#include "spc_fw.itp"
I run the simulation 6.5 ns for collecting data and I have total
256
water molecules.
Is there anything wrong with my .itp file?
I am using Gromacs VERSION 4.0.7.
Thanks
Nilesh
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
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/Support/Mailing_Lists/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/Support/Mailing_Lists
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
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/Support/Mailing_Lists/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/Support/Mailing_Lists
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
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/Support/Mailing_Lists/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/Support/Mailing_Lists
