[gmx-users] analysis of MD energy terms
Justin A. Lemkul
jalemkul at vt.edu
Fri May 28 03:44:57 CEST 2010
Moeed wrote:
> Hello Justin,
>
>
> [Actually, the reason I posted some parts of my mdp file was because you
> said I am sending redundant information so I just deleted the first
> lines of mdp file :) ]
>
What I was referring to was the meaningless grompp outputs, EM results, etc. An
.mdp file should always be posted in full if it is relevant.
>
> To get the xtc file initially I used nrexcl=3 in the top file. Then I
> used mdrun -rerun to breakdown energies on this xtc file with new top
> file nrexcl=7 to exclude all intramolecular nonbonded interactions (you
> definitely recall all these details :)
>
>
> 1- Does it matter in the first run (to get xtc) what option I use for
> nrexcl in top file. Actually, I thought maybe in the normal run I have
> to set this option to zero. I did that and I got lincs warning (the same
> applies for nrexcl=1). So I tried nrexcl=2 and I noticed the results are
> not the same as nrexcl 3. However, this means I am still having
> nonbonded intramolecular interaction between atoms that are less than 2
> bonds away. Am I right?
>
Haphazardly setting different values for nrexcl will break down the physics!
You must understand the intrinsic assumptions in most MD force fields before
playing with these values. If you set nrexcl = 1 or 2 then nonbonded
interactions are being calculated for atom pairs that are normally governed only
by bonded interactions. Anything less than 3 (or greater than three in most
cases, for that matter) is going to lead to a breakdown and unreasonable
results, if you manage to get the simulation running at all. That's why your
nrexcl = 7 has to be done as a rerun. If you tried this kind of simulation on
its own, the results would be gibberish.
>
> 2- I am wondering why in this breakdown I still see LJ and coloumb 1-4
> interactions while I have nrexcl=2.
>
Because 1-4 interactions are separated by three bonds. Using nrexcl = 2 has no
effect here.
>
> 3- In the second set of results (nrexcl=3) potential energy is 2766 >
> 2156 (the value for nrexcl==2). Logically, when 3 bonds away are
> excluded potential energy should be less I think!..
>
Hard to say. You're completely changing the balance of power between bonded and
nonbonded interactions.
>
> 4- Regarding the units, in chapter 2 of manual the units for energy is
> in KJ/mol. Here is this table I see heat capacity is given in J/mol K.
> Could you please kindly tell me the units what is the unit for energy
> temrs LJ or coulomb?
>
The heat capacity you get is not useful. I believe the method for calculating
it has been updated in the development code, but basically you'll get the same
result (or nearly so) for just about every simulation that you run, regardless
of what's in it :)
>
> 5- Very important question: If I want to compute the interaction
> energies BETWEEN molecuels, do I have to sum up LJ (SR) =(-3026 e.g) and
> coloumb (SR)=74 ? and can I conclude that there is a net attraction
> energy in the system?
>
That depends. Are there long range interactions? PME? Dispersion correction?
>
> 6- If I want to assure that equilibrium is reached, can I make the
> simulation time long enough such that T reaches 300 K? now it is 297 K.
>
You can make run for as much time as you want. Sounds to me like you're
relatively close to where you need to be, but you've got to reach the target
temperature, ensure that it's stable, and then and only then should you be
collecting data.
-Justin
>
> I appreciate your help.
>
> Thanks
>
> Moeed
>
>
>
> _mdrun -rerun on trajectory for nrexcl=2 in top file_
>
>
> Energy Average RMSD Fluct. Drift
> Tot-Drift
>
>
>
> -------------------------------------------------------------------------------
>
>
>
> Angle 4286.19 258.765 231.378 40.1269
> 401.349
>
>
>
> Ryckaert-Bell. 722.966 69.0616 61.1677 11.1046
> 111.068
>
>
>
> *LJ-14 498.432 31.2891 30.9138 1.67344
> 16.7378*
>
> * *
>
> *Coulomb-14 -322.264 11.5056 10.1507
> 1.87612 18.765*
>
> * *
>
> *LJ (SR) -3026.03 45.5198 38.6106 8.35011
> 83.5178*
>
> * *
>
> *Coulomb (SR) 74.4564 2.53605 2.53537 0.02034
> 0.203441*
>
>
>
> Coul. recip. -77.5276 2.3709 1.31777 -0.682623
> -6.8276
>
>
>
> Potential 2156.22 335.429 282.807 62.4689
> 624.814
>
>
>
> Kinetic En. 0 0 0
> 0 0
>
>
>
> Total Energy 2156.22 335.429 282.807 62.4689
> 624.814
>
>
>
> Temperature 0 0 0
> 0 0
>
>
>
> Pressure (bar) 1563.28 161.924 161.475 4.17097
> 41.718
>
>
>
> Cons. rmsd () 0 0 0
> 0 0
>
>
>
> Mu-X -0.00385475 0.893851 0.893682 -0.00602335
> -0.0602456
>
>
>
> Heat Capacity Cv: nan J/mol K (factor = nan)
>
>
>
>
> _
> _
>
> _mdrun -rerun on trajectory for nrexcl=3 in top file_
>
>
>
> Statistics over 5001 steps [ 0.0000 thru 10.0000 ps ], 14 data sets
>
> All averages are over 501 frames
>
>
>
>
>
> Energy Average RMSD Fluct. Drift
> Tot-Drift
>
> -------------------------------------------------------------------------------
>
> Angle 4400.09 284.259 237.268 54.2184
> 542.293
>
> Ryckaert-Bell. 1009.62 124.468 94.0379 28.2417
> 282.474
>
> *LJ-14 648.411 40.6844 39.926 2.70793
> 27.0848*
>
> **
>
> *Coulomb-14 -278.654 34.4227 11.8071
> 11.1987 112.01*
>
> **
>
> *LJ (SR) -3010.16 47.968 42.5388
> 7.67736 76.789*
>
> **
>
> *Coulomb (SR) 74.7846 2.5186 2.47007 0.170417
> 1.70451*
>
> Coul. recip. -78.0436 2.58606 1.54725 -0.717662
> -7.17805
>
> Potential 2766.05 454.268 342.141 103.497
> 1035.18
>
> Kinetic En. 0 0 0
> 0 0
>
> Total Energy 2766.05 454.268 342.141 103.497
> 1035.18
>
> Temperature 0 0 0
> 0 0
>
> Pressure (bar) 1973.01 177.068 172.803 13.3792
> 133.819
>
> Cons. rmsd () 0 0 0
> 0 0
>
> Mu-X -0.0523974 0.826473 0.826318 -0.00553723
> -0.0553834
>
> Heat Capacity Cv: nan J/mol K (factor = nan)
>
>
>
>
> ; Run control
> integrator = md
> dt = 0.002 ; ps !
> nsteps = 5000 ; total 1.0 ps.
> nstcomm = 1 ; frequency for center of mass motion
> removal
>
> ; Output control
> nstenergy = 10 ; frequency to write energies to energy
> file. i.e., energies and other statistical data are stored every 10 steps
> nstxout = 10 ; frequency to write
> coordinates/velocity/force to output trajectory file
> nstvout = 0
> nstfout = 10
> nstlog = 10 ; frequency to write energies to log file
> nstxtcout = 10
>
> ; Neighbor searching
> nstlist = 10 ; neighborlist will be updated at least
> every 10 steps
> ;ns_type = grid
>
> ; Electrostatics/VdW
> coulombtype = PME
> vdw-type = cut-off
> ; Cut-offs
> rlist = 1.0
> rcoulomb = 1.0
> rvdw = 1.0
>
> ; Temperature coupling Berendsen temperature coupling is on in
> two groups
> Tcoupl = berendsen
> tc-grps = HEX ;sol
> tau_t = 0.1 ;0.1
> ref_t = 300 ;300
>
> ; Pressure coupling: Pressure coupling is not on
> Pcoupl = no
> tau_p = 0.5
> compressibility = 4.5e-5
> ref_p = 1.0
>
> ; Velocity generation Generate velocites is on at 300 K.
> Manual p155
> gen_vel = yes
> gen_temp = 300.0
> gen_seed = 173529
>
> ; Bonds
> constraints = all-bonds
> constraint-algorithm = lincs
>
> pbc=xyz
>
>
>
>
>
--
========================================
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