[gmx-users] Energies in simulation

David van der Spoel spoel at xray.bmc.uu.se
Mon Jun 8 19:58:15 CEST 2009


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
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



More information about the gromacs.org_gmx-users mailing list