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

[pascal.baillod at epfl.ch]

Thank you very much David, thank you very much Mark! 

I actually already have a gromacs REMD implementation, that is an external
script which iteratively lauches jobs, computes protein potential energies (by
summing terms obtained with g_energy protein energy group output), computes swap
probabilities in fct of these energies, dose the possible swaps and relaunches
the new runs. But that is of course rather slow, as a call for grompp is
necessary every new run. 

1) Would there be a way to bypass grompp in the context of the external script?

In order to get prot-energy REMD in the gromacs source code: If I understand you
well, David, I could actually create a new routine "sum_epotPROT" in the file
src/mdlib/tgroup.c?.. It would look as follows:


void sum_epotPROT(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][0]);
  epot[F_COUL_SR]  = sum_v(grps->estat.nn,grps->estat.ee[egCOULSR][1]);
  epot[F_COUL_SR]  = sum_v(grps->estat.nn,grps->estat.ee[egCOULSR][2]);
  epot[F_LJ]       = sum_v(grps->estat.nn,grps->estat.ee[egLJSR][0]);
  epot[F_LJ]       = sum_v(grps->estat.nn,grps->estat.ee[egLJSR][1]);
  epot[F_LJ]       = sum_v(grps->estat.nn,grps->estat.ee[egLJSR][2]);
  epot[F_LJ14]     = sum_v(grps->estat.nn,grps->estat.ee[egLJ14][0]);
  epot[F_LJ14]     = sum_v(grps->estat.nn,grps->estat.ee[egLJ14][1]);
  epot[F_LJ14]     = sum_v(grps->estat.nn,grps->estat.ee[egLJ14][2]);
  epot[F_COUL14]   = sum_v(grps->estat.nn,grps->estat.ee[egCOUL14][0]);
  epot[F_COUL14]   = sum_v(grps->estat.nn,grps->estat.ee[egCOUL14][1]);
  epot[F_COUL14]   = sum_v(grps->estat.nn,grps->estat.ee[egCOUL14][2]);
  epot[F_COUL_LR] += sum_v(grps->estat.nn,grps->estat.ee[egCOULLR][0]);
  epot[F_COUL_LR] += sum_v(grps->estat.nn,grps->estat.ee[egCOULLR][1]);
  epot[F_COUL_LR] += sum_v(grps->estat.nn,grps->estat.ee[egCOULLR][2]);
  epot[F_LJ_LR]   += sum_v(grps->estat.nn,grps->estat.ee[egLJLR][0]);
  epot[F_LJ_LR]   += sum_v(grps->estat.nn,grps->estat.ee[egLJLR][1]);
  epot[F_LJ_LR]   += sum_v(grps->estat.nn,grps->estat.ee[egLJLR][2]);
/* 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][0]);
  epot[F_BHAM]     = sum_v(grps->estat.nn,grps->estat.ee[egBHAMSR][1]);
  epot[F_BHAM]     = sum_v(grps->estat.nn,grps->estat.ee[egBHAMSR][2]);
  epot[F_BHAM_LR]  = sum_v(grps->estat.nn,grps->estat.ee[egBHAMLR][0]);
  epot[F_BHAM_LR]  = sum_v(grps->estat.nn,grps->estat.ee[egBHAMLR][1]);
  epot[F_BHAM_LR]  = sum_v(grps->estat.nn,grps->estat.ee[egBHAMLR][2]);

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

2) Could that work? 

3) You invoke the hassle it will be to consider PME interactions of
protein-protein and protein-solute only.. So that is only done partially in the
"F_COUL_LR" terms above, since

/* lattice part of LR doesnt belong to any group
 * and has been added earlier
 */

Where would "earlier" be and how could I modify it?


4) I would then have to get the replica_exchange routine to use protein
potential energies, instead of total system potential energies. Would this be
possible changing the following line in /src/kernel/md.c

bExchanged = replica_exchange(log,mcr,repl_ex,state,ener[F_EPOT],step,t);

by

bExchanged = replica_exchange(log,mcr,repl_ex,state,ener[F_EPOTprot],step,t);


Thank you very much for your 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	
*******************************************************************************