[gmx-users] Energies in simulation
Payman Pirzadeh
ppirzade at ucalgary.ca
Mon Jun 8 19:54:33 CEST 2009
Well, when I set the box, I used 'editconf' command to rescale the box to
have the right density which was ~997. After simulation, I got the
following:
Energy Average RMSD Fluct. Drift
Tot-Drift
----------------------------------------------------------------------------
---
Density (SI) 956.765 150.266 130.596 0.257477
257.477
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:29 AM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] Energies in simulation
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-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