[gmx-users] Energies in simulation
Payman Pirzadeh
ppirzade at ucalgary.ca
Mon Jun 8 20:59:16 CEST 2009
Well, I tried to equilibrate the system, and I think the program has
included all the distortion that the system had been through. When I checked
the trajectory, the box was highly distorted as the simulation started, then
returned to its cubic shape. Also, I am a bit confused by the word
"strange". Well, based on manual the rvdw>rlist and it is. Also, authors in
the paper have tested U with various cutoff radii from 8-11A. The results
were very close! To make sure we are on the same page, the paper address is
J. Chem. Phys., Vol118, No.16 pp7401 (2003).
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:58 AM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] Energies in simulation
Payman Pirzadeh wrote:
> 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
The RMSD is very large. Are you sure this is in equilibrium? It could be
that your box is exploding. Please check as well that the strange twin
range cut-off that you are using is what the original authors used.
>
>
>
> 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