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

[pascal.baillod at epfl.ch]

Hello,

I thank you very much, David and Yang! Combining your 2 pieces of advice, I
finally managed to get something running.. It is still incorrect, as I have, for
every step, null protein potential energies (and a swapping probability of 100%,
obviously). Could that have something to do with some PME terms I should still
include, David? Let me post you my modifications, in 4 different files, as well
as an output sample, here below.. Thank you very much for your help!!

Regards,

Pascal

---------------------------------------------------------------------
1) src/mdlib/tgroup.c : summing prot epot 
--------------------------------------------------------------------

void sum_epotPROT(t_grpopts *opts,t_groups *grps,real epot[])
{
  int i;
                                                                               
                                                                             
  /* Accumulate energies */
  epot[F_COUL_SR]  = grps->estat.ee[egCOULSR][0] +
0.5*(grps->estat.ee[egCOULSR][1] + grps->estat.ee[egCOULSR][2]);
  epot[F_LJ]       = grps->estat.ee[egLJSR][0] + 0.5*(grps->estat.ee[egLJSR][1]
+ grps->estat.ee[egLJSR][2]);
  epot[F_LJ14]     = grps->estat.ee[egLJ14][0] + 0.5*(grps->estat.ee[egLJ14][1]
+ grps->estat.ee[egLJ14][2]);
  epot[F_COUL14]   = grps->estat.ee[egCOUL14][0] +
0.5*(grps->estat.ee[egCOUL14][1] + grps->estat.ee[egCOUL14][2]);
  epot[F_COUL_LR]  = grps->estat.ee[egCOULLR][0] +
0.5*(grps->estat.ee[egCOULLR][1] + grps->estat.ee[egCOULLR][2]);
  epot[F_LJ_LR]    = grps->estat.ee[egLJLR][0] + 0.5*(grps->estat.ee[egLJLR][1]
+ grps->estat.ee[egLJLR][2]);
/* lattice part of LR doesnt belong to any group
 * and has been added earlier
 */
  epot[F_BHAM]     = grps->estat.ee[egBHAMSR][0] +
0.5*(grps->estat.ee[egBHAMSR][1] + grps->estat.ee[egBHAMSR][2]);
  epot[F_BHAM_LR]     = grps->estat.ee[egBHAMLR][0] +
0.5*(grps->estat.ee[egBHAMLR][1] + grps->estat.ee[egBHAMLR][2]);
                                                                               
                                                                             
  epot[F_EPOTprot] = 0;
  for(i=0; (i<F_EPOTprot); i++)
    if (i != F_DISRESVIOL && i != F_ORIRESDEV && i !=
F_DIHRESVIOL)
      epot[F_EPOTprot] += epot[i];
}


...added at end of the file.


---------------------------------------------------------------------
2) include/types/idef.h 
--------------------------------------------------------------------

typedef atom_id t_iatom;

/* this MUST correspond to the 
   t_interaction_function[F_NRE] in gmxlib/ifunc.c */
enum {
  F_BONDS,
  F_G96BONDS,

/* other contributions.... */

  F_DVDL,
  F_DVDLKIN,
  F_NRE,       /* This number is for the total number of energies	*/
  F_EPOTprot   /* my addition, after F_NRE as suggested by Yang	*/
};


---------------------------------------------------------------------
3) src/gmxlib/ifunc.c (update interaction data with new prot epot, as indicated
by David)
-----------------------------------------------------------------------


/* this MUST correspond to the enum in include/types/idef.h */
const t_interaction_function interaction_function[F_NRE]=
{
  def_bond   ("BONDS",    "Bond",            2, 2, 2,  eNR_BONDS,  bonds         ),

/* other contributions.... */

  def_bond   ("G96BONDS", "G96Bond",         2, 2, 2,  eNR_BONDS,  g96bonds      
  def_nofc   ("DV/DL",    "dVpot/dlambda"    ),
  def_nofc   ("DK/DL",    "dEkin/dlambda"    ),
  def_nofc   ("EPOTprot",     "Potential Protein"        )
};




---------------------------------------------------------------------
4) src/kernel/md.c : Printing protein energies
---------------------------------------------------------------------
    
    if ( MASTER(cr)){
      fprintf(log,"Epot protein : %12.8f\n", ener[F_EPOTprot] );
    }

..added at line 935. This prints out the protein potential energy 9 times, I
guess I have to put it elsewhere..

---------------------------------------------------------------------
5) Output obtained:
---------------------------------------------------------------------


           Step           Time         Lambda
             50        0.07500        0.00000
                                                                               
                                                                             
   Energies (kJ/mol)
        G96Bond       G96Angle    Proper Dih.  Improper Dih.          LJ-14
    1.52789e+03    1.63171e+03    6.80293e+02    5.20036e+02    2.75589e+02
     Coulomb-14        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential
    1.04654e+04    8.45580e+04   -5.65189e+05   -6.30367e+04   -5.28567e+05
    Kinetic En.   Total Energy    Temperature Pressure (bar)
    9.60138e+04   -4.32553e+05    2.96506e+02    5.24819e+01
                                                                               
                                                                             
Epot protein :   0.00000000
Replica exchange at step 50 time 0.075
Repl 0 <-> 1  dE =  1.697e-46  dpV =  4.220e-07  d =  4.220e-07
Repl ex  0 x  1
Repl pr   1.0
                                                                               
                                                                             
Epot protein :   0.00000000
Epot protein :   0.00000000
Epot protein :   0.00000000
Epot protein :   0.00000000
Epot protein :   0.00000000
Epot protein :   0.00000000
Epot protein :   0.00000000
Epot protein :   0.00000000
Epot protein :   0.00000000




*******************************************************************************
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	
*******************************************************************************