[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
> unfortunately) to read all your MD novels at ones.
>
> 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.
>
> I have also some other inquiries about this system:
>
> - 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
> am asking
> 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
========================================
More information about the gromacs.org_gmx-users
mailing list