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

[Yang Ye]

pascal.baillod at epfl.ch wrote:

>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];
>}
>  
>
you may do so in include/type/idef.h. Add one more item after F_NRE
This is one trick so that you just need to change the memory allocation 
in md.c to snew(ener,F_NRE+x); x is 1 for this case. If you put new 
energy part before F_NRE, there will be compile error (or running error? 
I couldn't remember).

This is what I could help in coding.

-- 
/Regards,/
Yang Ye
/Computational Biology Lab
School of Biological Sciences
Nanyang Technological University
Singapore
Tel: 6316-2884
/