# [gmx-users] Re: NPT simulation

Justin A. Lemkul jalemkul at vt.edu
Sun Aug 29 20:39:21 CEST 2010

```
Moeed wrote:
> Hello Dr. Chaban,
>
> System size is 5 *5*3 nm when density reaches around 600 SI. Can you
> please explain what do you mean by "What about dependence heat of
> evaporation vs system size ?". I dont know how to compute vap. heat in
> gromacs!
>
> Unfortunately, I am still unclear about energy units. :(
> Yes, only when I use
> g_energy -f ener.edr -o energy.xvg -nmol 8       I get pretty similar
> values as you obtained (around 2200 KJ/molecule !? ) but then units
> should be energy per molecule. I am confused because Before you had said
> units are per mole. Without using -nmol 8 option total energy is around
> 20,000 (apparently g_energy by default gives  KJ/mole)
>
> Another reason I am confused is because
>
> g_energy -f ener.edr -o energy.xvg according to manual chapter 2 should
> give energies in KJ/mol and not system. I did some simulations for
> different number of polymer chains in a big box. (1, 2, 4, 8 chains).
> The magnitude of energy values increase as number of chains increase..My
> question is if values are in the units of KJ/mol why let's say nonbonded
> interaction energy values are increasing in value as number of molecuels
> goes up? Should not the values be the same if we are referring to MOLE
> number of molecules?
>
> I know manual can not be wrong but from what I see it makes more sense
> if unit is energy/system!
>

When dividing the energy of a system by the number of molecules (in a
homogeneous system), you are extracting what I believe is commonly referred to
as "configurational energy" which, for relatively simple systems, should
converge fairly quickly.  The reason your total energy values increase with
system size is a simple matter of potential energy.  More interactions mean that
the magnitude of the potential will increase, and likewise with the kinetic
energy, more particles that have velocity imply a greater sum.

You can convince yourself of this fact by running relatively short simulations
of water boxes, using, i.e. spc216.gro and then a larger construct from it.
Within a few hundred ps, you should get reasonably converged configurational
energies, though of course the total energy will be larger simply by virtue of
the system size.

Energies are in kJ/mol, no question.  The extrapolation would be that a mole of
a given species in the given configuration would have this energy in kJ.  For a
single molecule, the energy is whatever value you obtain divided by Avogadro's
number, which is a relatively unimportant quantity if you are interested in bulk
dynamics.

And for the record, I do *not* enjoy novel-length posts, rather complete ones
with sufficient detail to diagnose the issue :)

-Justin

> Thanks,
>
>
> 2010/8/27 Vitaly Chaban <vvchaban at gmail.com <mailto:vvchaban at gmail.com>>
>
>     Dear Moeed:
>
>     1. Please note, this is Justin here who likes long descriptions but
>     I like those which are as short as possible. :-) Please formulate
>     your questions more specifically because I am too lazy (and busy
>
>     2. If I said that the energy was given per mole, than it is given
>     per mole of molecules, not systems, be sure. Use g_energy -nmol  XXX
>     where it is reasonable. I think your energy vlue before was correct,
>     just divide it by NMOL.
>
>     3. Use NVT (no barostat) if you really need some speficif density.
>     Calculate the box size by hand.
>
>     4. Polymers are rather specific MD systems. I think your 1 PE + 343
>     hexanes should be a minimal good composition. By the way, what is
>     the size of MD box in nmxnmxnm? What about dependence heat of
>     evaporation vs system size ?
>
>     5. If I am not mistaken "LJ-SR" is a total Lennard-Jones energy of
>     system except 1-4 interactions. If I am not true, please somebody
>     correct me.
>
>     6. Just equilibrate the system to avoid volume dances in the
>     production run.
>
>     Good luck!
>
>     --
>     Dr. Vitaly V. Chaban
>     Department of Chemistry
>     University of Rochester,
>     Rochester, NY, U.S.A.
>
>
>
>     Dear Dr.Chaban,
>
>     Thanks for your message.  Now I  have to groups: PE is a single
>     chain and
>     HEX (343 hexane molecules). I am interested in nonbonded interaction
>     energies between PE and solvent as well as between PE-PE and HEX
>     -HEX. I ran
>     the simulation for 300ps and the system is equilibrating (300K and
>     30bar and
>     also total energy relaxes to a 46000KJ/mol). Actually I can not
>     trust these
>     results since if you remember you got much less energy terms.. I was
>     shameful of bothering you and discussing about that at that time since I
>     thought it does not make sense getting different results with a same
>     parameter file and gerometry.  You said it might be because I am getting
>     results in a different unit but in fact to get the energy break down
>     I just
>     issue g_energy command and dont touch anything. So energy values listed
>     below should be in KJ/mol (MOLE by default and not molecule). Since
>     H and C
>     are the only atoms present in the system I have switched off
>     electrostatics
>     and total energy is about 46000KJ/mol. I would be thankful if you
>     could take
>     a look at my system in a proper time.
>
>
>     - the density I am getting is about 650 SI and I need sth about 400. In
>     general the key parameter to change density in NPT is pressure
>     right? so I
>     could use P=20bar to get less density and it is kind of try and
>     error . am I
>     right?
>
>     - I have still a foggy image of the system size in simulations. How
>     can I
>     decide on number of polymer and solvent particles in system? (of
>     course the
>     smaller the better in terms of simulation time..) but the reason I
>     this is because I ran simulations for different number of polymer chains
>     (1,2,4,8) and I saw energy terms as total energy increase in value as
>     systmes becomes larger. If units are in energy/mol so why is this
>     happening?
>     It does not make sense to calculate energy for a certain (MOLE)
>     amount of
>     molecules and get different  values for different system size. (in
>     brief: I
>     can no t decide on number of polymer and hexane molecules)
>
>     -now in this system there is only one polymer chain and energy breakdown
>     shows energy for PE-PE LJ SR and 1-4. Does it make sanse to have SR
>     within a
>     single chain? (does this come from 1-5,1-6....interactions on chain?
>     in case
>     of hexane this can be only 1-5,1-6,1-7 becase there are no more than
>     7 bonds
>     away each atom) LJ-SR:HEX-HEX    = -7651.06,  & LJ-SR:PE-PE =-96.7831
>
>     -at the beginning of the simulation I see the system is decreasing
>     in size
>     and after about 10ps box expands again and after 20th ps starts
>     contracting
>     again. (volume starts from 8000 nm3 and reaches around 500 then
>     increases to
>     1600 nm3 and after that converges nicely to a plateau..). I tried both
>     semiisotropic and isotropic options but the same behaviour has been
>     observed.
>
>
>     I appreciate any help and comment...
>     Moeed :)
>
>     **********************************************************************
>
>     Energy                      Average       RMSD     Fluct.      Drift
>     Tot-Drift
>     -------------------------------------------------------------------------------
>     Bond                         8536.6    562.056    441.004    4.02362
>     1207.09
>     Angle                       13076.1    682.442    576.028     4.2256
>     1267.69
>     Ryckaert-Bell.              3276.27    329.245    302.614   -1.49785
>     -449.355
>     LJ-14                       1762.95    58.3619    57.9592  0.0790304
>     23.7092
>     Coulomb-14                        0          0          0
>     0          0
>     LJ (SR)                    -8457.99    2676.83    1915.36   -21.5925
>     -6477.77
>     Coulomb (SR)                      0          0          0
>     0          0
>     Potential                   18193.9    1946.39    1467.66   -14.7621
>     -4428.64
>     Kinetic En.                 27009.8    325.523     325.51  0.0343522
>     10.3057
>     Total Energy                45203.7       1970    1501.36   -14.7277
>     -4418.33
>     Temperature                 299.914    3.61457    3.61442 0.000381443
>     0.114433
>     Pressure (bar)              30.8942    533.069    533.059  0.0378739
>     11.3622
>     Box-X                       5.41674    2.51691    2.19253 -0.0142718
>     -4.28156
>     Box-Y                       5.41674    2.51691    2.19253 -0.0142718
>     -4.28156
>     Box-Z                       5.02579    2.60455    2.00013 -0.0192635
>     -5.77908
>     Volume                      306.955    929.783    856.612   -4.17472
>     -1252.42
>     Density (SI)                530.808    191.837      139.2    1.52423
>     457.271
>     #Surf*SurfTen              -3.00385    3396.99    3396.96   0.165929
>     49.779
>     Coul-SR:PE-PE                     0          0          0
>     0          0
>     LJ-SR:PE-PE                -96.7831    12.0843     11.971  0.0190581
>     5.71744
>     Coul-14:PE-PE                     0          0          0
>     0          0
>     LJ-14:PE-PE                 145.957    10.2536    10.2529 0.00134961
>     0.404885
>     Coul-SR:PE-HEX                    0          0          0
>     0          0
>     LJ-SR:PE-HEX               -710.145    226.458    156.934   -1.88519
>     -565.559
>     Coul-14:PE-HEX                    0          0          0
>     0          0
>     LJ-14:PE-HEX                      0          0          0
>     0          0
>     Coul-SR:HEX-HEX                   0          0          0
>     0          0
>     LJ-SR:HEX-HEX              -7651.06    2451.43    1758.13   -19.7263
>     -5917.92
>     Coul-14:HEX-HEX                   0          0          0
>     0          0
>     LJ-14:HEX-HEX               1616.99    56.0996    55.6947  0.0776808
>     23.3043
>
>     *************************************************************************
>     pbc              =  xyz                   ; use priodic BCs in all
>     directions
>     energygrps          =  PE HEX
>
>     ;        Run control
>     integrator          =  md
>     dt                  =  0.001
>     nsteps              =  300000
>     nstcomm             =  1
>
>     ;        Output control
>     nstenergy           =  100
>     nstxout             =  100
>     nstvout             =  0
>     nstfout             =  0
>     nstlog              =  1000
>     nstxtcout          =  0
>
>     ;        Neighbor searching
>     nstlist             =  10                 ;
>     ns_type             =  grid
>
>     ;        Electrostatics/VdW
>     ;coulombtype         =  PME
>     vdw-type            =  Shift
>     rcoulomb-switch     =  0
>     rvdw-switch         =  0.9
>
>     ;        Cut-offs
>     rlist               =  1.1
>     rcoulomb            =  1.1 ;1.0
>     rvdw                =  1.0
>
>     ;        Temperature coupling
>     Tcoupl              =  v-rescale
>     tc-grps             =  System
>     tau_t               =  0.1
>     ref_t               =  300
>
>     ;        Pressure coupling
>     Pcoupl              =  Parrinello-Rahman;berendsen
>     Pcoupltype          =  semiisotropic    ;isotropic ;
>     tau_p               =  1      1
>     compressibility     =  4.5e-5 4.5e-5
>     ref_p               =  30.0    30.0
>
>     ;        Velocity generation
>     gen_vel             =  yes
>     gen_temp            =  300.0
>     gen_seed            =  173529
>
>     ;        Bonds
>     constraints         =  none
>     constraint-algorithm = lincs
>
>
>
>
>
> --
> Moeed Shahamat
> Graduate Student (Materials Modeling Research Group)
> McGill University- Department of Chemical Engineering
> Montreal, Quebec H3A 2B2, Canada
> Web:http://mmrg.chemeng.mcgill.ca/pages/current-group-members/moeed-shahamat.php
> Web:http://mmrg.chemeng.mcgill.ca/
>

--
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================

```