[gmx-users] Energies in simulation
Payman Pirzadeh
ppirzade at ucalgary.ca
Mon Jun 8 19:24:51 CEST 2009
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.
Payman
-----Original Message-----
From: gmx-users-bounces at gromacs.org [mailto:gmx-users-bounces at 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-bounces at gromacs.org [mailto:gmx-users-bounces at 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 at 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-request at 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
spoel at xray.bmc.uu.se spoel at gromacs.org http://folding.bmc.uu.se
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
_______________________________________________
gmx-users mailing list gmx-users at 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-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
More information about the gromacs.org_gmx-users
mailing list