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

Yang Ye leafyoung81-group at yahoo.com
Thu Jan 19 20:33:30 CET 2006


I think, with this modification

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

F_EPOTprot can be added before F_NRE. David, you are right!

Go back to your problem, I think that you are just one step away from "Wow".

First, setup a debug environment. A note to you is that always write 
test and debug code along with your code.
Set nstep=2 and nstenergy=1. Write code in sum_epotPROT to print the 
values inside epot[].

Next, I just realized, where do you call sum_epotPROT in md.c?

=)

Yang Ye

pascal.baillod at epfl.ch wrote:
> 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	
> *******************************************************************************
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-developers-request at gromacs.org.
>
>   

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



More information about the gromacs.org_gmx-developers mailing list