[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