[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