Payman Pirzadeh wrote:
Here is the content of .itp file which I developed:
; This is an itp file to describe water's six-site model by H. Nada and J.P.
J. M. Van der Eerden, J. Chem. Phys. Vol.118, no.16, pp7401-7413 (2003)
; This model is a combination of TIP4P and TIP5P. It has three LJ sites and
3 Coulomb sites
; O-H bond length is 0.980A, HOH angle is 108.00degrees, LOL angle is 111.00
degrees, O-M and O-L are about 0.230A and 0.8892A respectively
[ defaults ]
; non-bondedtype combrule genpairs FudgeLJ
FudgeQQ N
1 2 NO
[ atomtypes ]
;name mass charge ptype c6 c12
OW 15.9994 0.0 A 0.3115 0.714845562
HW 1.00800 0.477 A 0.0673 0.11541
MW 0.000 -0.866 D 0.00 0.00
LW 0.00 -0.044 D 0.00 0.00
[ moleculetype ]
;molname nrexcl
SOL 1
[ atoms ]
; nr atomtype resnr residuename atom cgnr charge
1 OW 1 SOL OW 1 0.0
2 HW 1 SOL HW1 1 0.477
3 HW 1 SOL HW2 1 0.477
4 MW 1 SOL MW 1 -0.866
5 LW 1 SOL LP1 1 -0.044
6 LW 1 SOL LP2 1 -0.044
[ settles ]
; OW function doh dhh
1 1 0.0980 0.15856
[ dummies3 ]
; These set of parameters are for M site which can be easily calculated
using TIP4P calculations from tip4p.itp
; So, it will be described as dummy site 3: r(v)= r(i) + a*(r(i)-r(j)) +
b*(r(i)-r(k))
; const = |OH|/{|OH|*cos(HOH/2)} => Due to vector algebra a=b=const/2.
Remember that OM is in the same direction of OH bonds.
; Remember this site is in the same plane of OH bonds; so, its function 1
;
; site from function a b
4 1 2 3 1 0.199642536 0.199642536
; Now we define the position of L sites which can be obtained from tip5p.itp
; So, it will be described as dummy site 3out: r(v) = r(i) + a*(r(i)-r(j)) +
b*(r(i)-r(k)) + c*(r(ij)Xr(ik))
; const1 = {|OL|*cos(LOL/2)}/{|OH|*cos(HOH/2)} => Due to vector algebra
|a|=|b|=const/2. since the lone pairs are in opposite direction of OH bonds,
a minus sign is added. This part is similar to M site.
; const2 = {|OL|*sin(LOL/2)}/{|OH|*|OH|*sin(HOH)} => The denominator is the
magnitude of vector product of OH bonds.
; This sites are tetrahedral sites; so, its function 4
;
; site from function a b c
5 1 2 3 4 -0.437172388 -0.437172388
8.022961206
6 1 2 3 4 -0.437172388 -0.437172388
-8.022961206
[ exclusions ]
1 2 3 4 5 6
2 1 3 4 5 6
3 1 2 4 5 6
4 1 2 3 5 6
5 1 2 3 4 6
6 1 2 3 4 5
And here is the mdp file which I used for the simulation run:
cpp = cpp
include = -I../top
define =
; Run control
integrator = md
dt = 0.001 ;1 fs
nsteps = 1000000 ;10 ns
comm_mode = linear
nstcomm = 1
;Output control
nstxout = 5000
nstlog = 5000
nstenergy = 5000
nstxtcout = 1000
xtc_grps =
energygrps =
; Neighbour Searching
nstlist = 10
ns_type = grid
rlist = 0.9
; Electrostatistics
coulombtype = PME
rcoulomb = 0.9
epsilon_r = 1
; Vdw
vdwtype = cut-off
rvdw = 1.2
DispCorr = EnerPres
;Ewald
fourierspacing = 0.12
pme_order = 4
ewald_rtol = 1e-6
optimize_fft = yes
; Temperature coupling
tcoupl = Berendsen
tc-grps = System
tau_t = 0.1
ref_t = 300
; Pressure Coupling
Pcoupl = Berendsen
Pcoupltype = isotropic
tau_p = 1.0
compressibility = 5.5e-5
ref_p = 1.0
gen_vel = yes
The expected Potential energy is supposed to be around -41.5kJ/mol while my
potential is around -22.2kJ/mol. I calculated the energies by g_energy
command.
And do yo have the right density?
Payman
-----Original Message-----
From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org]
On Behalf Of David van der Spoel
Sent: June 8, 2009 11:13 AM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] Energies in simulation
Payman Pirzadeh wrote:
To the best of my knowledge no, but how can I check that?
A. read the original paper: is your topology correct? Are the simulation
parameters the same?
B. post the itp file here and mdp file and specify energy and expected
energy. How about energy units?
-----Original Message-----
From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org]
On Behalf Of David van der Spoel
Sent: June 8, 2009 11:06 AM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] Energies in simulation
Payman Pirzadeh wrote:
Hi,
I am using my own water model which I developed its .itp file. When
simulation is done after 1ns and energy is kinetic and potential
energies are analyzed, the kinetic value is almost OK, but the potential
energy is almost half of the value reported in literature and another MD
code that I am currently using. I double-checked the parameters I gave
in the .itp with TIP4P and TIP5P to make sure everything is correct in
format and unit. But I can not figure out the problem. Any ideas?
Is there any self-energy involved (i.e. a monomer energy that yo have to
subtract)?
Payman
------------------------------------------------------------------------
_______________________________________________
gmx-users mailing list gmx-users@gromacs.org
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 gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
--
David.
________________________________________________________________________
David van der Spoel, PhD, Professor of Biology
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596, 75124 Uppsala, Sweden
phone: 46 18 471 4205 fax: 46 18 511 755
sp...@xray.bmc.uu.se sp...@gromacs.org http://folding.bmc.uu.se
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
_______________________________________________
gmx-users mailing list gmx-users@gromacs.org
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 gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
