[gmx-users] Computing prot pot ener in the gmx code

[pascal.baillod at epfl.ch]

Hello David,

Thanks for your indications concerning the summing of potential energy terms in
gromacs. (For reminder, I am trying to get them in order to carry out REMD
simulations in function of the potential energy of a solute and not of the
entire system.) I checked the tgroup.c file and found the following:


void sum_epot(t_grpopts *opts,t_groups *grps,real epot[])
{
  int i;

  /* Accumulate energies */
  epot[F_COUL_SR]  = sum_v(grps->estat.nn,grps->estat.ee[egCOULSR]);
  epot[F_LJ]       = sum_v(grps->estat.nn,grps->estat.ee[egLJSR]);
  epot[F_LJ14]     = sum_v(grps->estat.nn,grps->estat.ee[egLJ14]);
  epot[F_COUL14]   = sum_v(grps->estat.nn,grps->estat.ee[egCOUL14]);
  epot[F_COUL_LR] += sum_v(grps->estat.nn,grps->estat.ee[egCOULLR]);
  epot[F_LJ_LR]   += sum_v(grps->estat.nn,grps->estat.ee[egLJLR]);
/* lattice part of LR doesnt belong to any group
 * and has been added earlier
 */
  epot[F_BHAM]     = sum_v(grps->estat.nn,grps->estat.ee[egBHAMSR]);
  epot[F_BHAM_LR]  = sum_v(grps->estat.nn,grps->estat.ee[egBHAMLR]);

  epot[F_EPOT] = 0;
  for(i=0; (i<F_EPOT); i++)
    if (i != F_DISRESVIOL && i != F_ORIRESDEV && i != F_DIHRESVIOL)
      epot[F_EPOT] += epot[i];
}


However, these components apparently refer to the entire system, and not to a
subpart, as can be defined in the .mdp files with the keyword

energygrps               = Protein Other

The .edr file obtained when running a simulation with this keywords contains the
following potential energy components:



 1=       G96Bond   2=      G96Angle   3=   Proper Dih.   4= Improper Dih.
   5=         LJ-14   6=    Coulomb-14   7=       LJ (SR)   8=  Coulomb (SR)
   9=  Coul. recip.  10=     Potential  11=   Kinetic En.  12=  Total Energy
  13=   Temperature  14=Pressure (bar)  15=         Box-X  16=         Box-Y
  17=         Box-Z  18=        Volume  19=  Density (SI)  20=            pV
  21=        Vir-XX  22=        Vir-XY  23=        Vir-XZ  24=        Vir-YX
  25=        Vir-YY  26=        Vir-YZ  27=        Vir-ZX  28=        Vir-ZY
  29=        Vir-ZZ  30= Pres-XX (bar)  31= Pres-XY (bar)  32= Pres-XZ (bar)
  33= Pres-YX (bar)  34= Pres-YY (bar)  35= Pres-YZ (bar)  36= Pres-ZX (bar)
  37= Pres-ZY (bar)  38= Pres-ZZ (bar)  39= #Surf*SurfTen  40=  Pcoupl-Mu-XX
  41=  Pcoupl-Mu-YY  42=  Pcoupl-Mu-ZZ  43=          Mu-X  44=          Mu-Y
  45=          Mu-Z  46=Coul-SR:Protein-Protein  47=LJ-SR:Protein-Protein 
48=Coul-14:Protein-Protein
  49=LJ-14:Protein-Protein  50=Coul-SR:Protein-Other  51=LJ-SR:Protein-Other 
52=Coul-14:Protein-Other
  53=LJ-14:Protein-Other  54=Coul-SR:Other-Other  55=LJ-SR:Other-Other 
56=Coul-14:Other-Other
  57=LJ-14:Other-Other  58=     T-Protein  59=       T-Other  60=  Lamb-Protein
  61=    Lamb-Other


To obtain the potential energy due to protein-protein and protein-water
interactions (excluding all water to water contributions), one could sum up the
energy components obtained with the input numbers 48-53 given to g_energy. I
would like to do this summation directly during the gromacs run (changing
slightly the source code) and feed this protein potential energy to the
replica-exchange subroutine. Where could I find this in the source code?

Thank you very much for any help!

Pascal

*******************************************************************************
Pascal Baillod (PhD student) 
*******************************************************************************
Swiss Federal Institute of Technology EPFL	        Tel: +41-(0)21-693-0322
Institute of Chemical Sciences and Engineering ,	Fax: +41-(0)21-693-0320
Laboratory of Computational Chemistry and Biochemistry	pascal.baillod at epfl.ch
Room BCH 4121, Avenue Forel,	                        http://lcbcpc21.epfl.ch
CH-1015 Lausanne	
*******************************************************************************