[gmx-users] Energies in simulation

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


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