[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