[gmx-developers] REMD simulations based on protein potential energy

pascal.baillod at epfl.ch pascal.baillod at epfl.ch
Wed Mar 8 11:20:45 CET 2006


Hi!

I'm back to fixing up a protein-potential-energy based gromacs REMD
implementation.. I have more or less everything fixed so that 

1) I sum up all prot bonded and prot-prot and prot-solvent interaction terms, in
the routine sum_epotPROT I added in the file tgroup.c (see code below) and store
the result in the new var ener(F_EPOTprot).

2) I use the resulting new energy term (ener(F_EPOTprot), 59th el of the "ener"
data structure) as potential energy to be used in the Monte-Carlo probability
defining whether 2 REMD replicas will be swapped or not.

It works fine and I obtain, as expected, small numbers for the potential
energies and high swapping probabilities (about 90%). 

I now only have to make sure that I have these protein potential energies
correct (I get something between -200 and 0 kJ/mol), which I suspect is not the
case. For reminder, I work with a 103 residue protein solvated in spc water,
resulting in a 38373 atom system. 

1) In order to cross check, I am having mdrun write protein and solvent energy
terms to the .edr file. With g_energy, I can obtain all prot bonded and
prot-prot and prot-solvent interaction terms. Summing them up (see the detail of
the terms I sum up below), I obtain the potential energy of the protein and of
its interaction with solvent molecules, which is around +200 to +600 kJ/mol,
higher than the one computed in my ener(F_EPOTprot) term (something between -200
and 0 kJ/mol)..

2) If I remember well, David once told me about missing PME terms I would have
to get elsewhere in the code (elsewhere than in tgroup.c).. Where would that be?

Many thanks for any help!!!

Pascal



----------------------------------------------------------------
MY NEW TGROUP.C ROUTINE AND G_ENERGY OPTIONS:
----------------------------------------------------------------



1) sum_epotPROT routine, computing new prot epot term ener(F_EPOTprot) in tgroup.c
--------------------------------------------------------------------


void sum_epotPROT(t_grpopts *opts,t_groups *grps,real epot[])
{
  int i;
  epot[F_EPOTprot]=0;
  
  /* 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_EPOTprot] += epot[F_COUL_SR];
  epot[F_LJ]       = grps->estat.ee[egLJSR][0] + 0.5*(grps->estat.ee[egLJSR][1]
+ grps->estat.ee[egLJSR][2]);
  epot[F_EPOTprot] += epot[F_LJ];
  epot[F_LJ14]     = grps->estat.ee[egLJ14][0] + 0.5*(grps->estat.ee[egLJ14][1]
+ grps->estat.ee[egLJ14][2]);
  epot[F_EPOTprot] += epot[F_LJ14];
  epot[F_COUL14]   = grps->estat.ee[egCOUL14][0] +
0.5*(grps->estat.ee[egCOUL14][1] + grps->estat.ee[egCOUL14][2]);
  epot[F_EPOTprot] += epot[F_COUL14];
  epot[F_COUL_LR]  = grps->estat.ee[egCOULLR][0] +
0.5*(grps->estat.ee[egCOULLR][1] + grps->estat.ee[egCOULLR][2]);
  epot[F_EPOTprot] += epot[F_COUL_LR];
  epot[F_LJ_LR]    = grps->estat.ee[egLJLR][0] + 0.5*(grps->estat.ee[egLJLR][1]
+ grps->estat.ee[egLJLR][2]);
  epot[F_EPOTprot] += epot[F_LJ_LR];
  
  
/*  fprintf(stderr,"grps->estat.nn : %d\n",grps->estat.nn);        */
/*  fprintf(stderr,"before loop : F_EPOTprot : %d\n",F_EPOTprot);  */
  
/* 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) {
/*      fprintf(stderr,"in loop : F_EPOTprot : %d, i= %d\n",F_EPOTprot,i);       */
/*      fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[i],
grps->estat.ee[i][0], grps->estat.ee[i][1],grps->estat.ee[i][2]);    */
/*       epot[F_EPOTprot] += epot[i];  */
/*      fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[i],
grps->estat.ee[i][0], grps->estat.ee[i][1],grps->estat.ee[i][2]);    */
/*      fprintf(stderr,"epot prot: %d \t %f\n", i,epot[i]);        */          
            
    }
  }
  
  fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[F_COUL_SR],
grps->estat.ee[egCOULSR][0],
grps->estat.ee[egCOULSR][1],grps->estat.ee[egCOULSR][2]); 
  fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[F_LJ],
grps->estat.ee[egLJSR][0], grps->estat.ee[egLJSR][1],grps->estat.ee[egLJSR][2]);
  fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[F_LJ14],
grps->estat.ee[egLJ14][0], grps->estat.ee[egLJ14][1],grps->estat.ee[egLJ14][2]);
  fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[F_COUL14],
grps->estat.ee[egCOUL14][0],
grps->estat.ee[egCOUL14][1],grps->estat.ee[egCOUL14][2]);
  fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[F_COUL_LR],
grps->estat.ee[egCOULLR][0],
grps->estat.ee[egCOULLR][1],grps->estat.ee[egCOULLR][2]);
  fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[F_LJ_LR],
grps->estat.ee[egLJLR][0], grps->estat.ee[egLJLR][1],grps->estat.ee[egLJLR][2]);
  fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[F_BHAM],
grps->estat.ee[egBHAMSR][0],
grps->estat.ee[egBHAMSR][1],grps->estat.ee[egBHAMSR][2]);
  fprintf(stderr,"epot prot: %d \t %f \t %f \t %f \t %f \n", i,epot[F_BHAM_LR],
grps->estat.ee[egBHAMLR][0],
grps->estat.ee[egBHAMLR][1],grps->estat.ee[egBHAMLR][2]);    
  fprintf(stderr,"epot prot from tgroup: \t %f \n", epot[F_EPOTprot]);         
                                                                            
}






2) g_energy terms I sum up in order to get a control value (that should be the
same) for the protein potential energy
--------------------------------------------------------------------


Select the terms you want from the following list
End your selection with 0
   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

TERMS I SUM UP:
1 2 3 4 46 47 48 49 50 51 52 53 0




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



More information about the gromacs.org_gmx-developers mailing list