eval.cc

c++编写的并行拉马克遗传算法的程序。实现分析对接程序

/********************************************************************
     These are the functions associated with the evaluation object.

                                rsh 9/95
********************************************************************/
// #include <iostream.h>
// #include <fstream.h>
#include <math.h>
#include "eval.h"

extern FILE *logFile;

#ifdef sgi
    #include <stdio.h>
    #include <string.h>
    #include "structs.h"
    #include <ieeefp.h>
#else
    extern "C"
    {
        #include <stdio.h>
        #include <string.h>
        #include "structs.h"
    }
#endif

#ifdef sun
    #include <ieeefp.h>
#endif

/*  The chromosome is assumed to have a layout like this -

       | x | y | z | nx | ny | nz | ang | tor1 | ... | tor N |

    where
       x is the x translation
       y is the y translation
       z is the z translation
       nx is the x component of the quaternion
       ny is the y component of the quaternion
       nz is the z component of the quaternion
       ang is the angle portion of the quaternion
       tor 1, ..., tor N are the ntor torsion angles
*/

void make_state_from_rep(Representation **rep, State *stateNow)
/*
    This routine modifies the various components of stateNow to correspond
    to the chromosome.  
*/
{
   register int i;

#ifdef DEBUG
   (void)fprintf(logFile, "eval.cc/make_state_from_rep(Representation **rep, State *stateNow)\n");
#endif /* DEBUG */

   //  Do the translations
   stateNow->T.x = rep[0]->gene(0).real;
   stateNow->T.y = rep[1]->gene(0).real;
   stateNow->T.z = rep[2]->gene(0).real;

   //  Set up the quaternion
   stateNow->Q.nx = rep[3]->gene(0).real;
   stateNow->Q.ny = rep[3]->gene(1).real;
   stateNow->Q.nz = rep[3]->gene(2).real;
   stateNow->Q.ang = rep[4]->gene(0).real;
   
   //  Copy the angles
   for (i=1; i<=stateNow->ntor; i++) {
      stateNow->tor[i-1] = rep[4]->gene(i).real;
   }

   mkUnitQuat(&(stateNow->Q));
}

double Eval::operator()(Representation **rep)
{
   register int i;
   int   B_outside = 0;
   int   I_tor = 0;
   int   indx = 0;
   double energy = 0.0;

#ifdef DEBUG
    (void)fprintf(logFile,"eval.cc/double Eval::operator()(Representation **rep)\n");
#endif /* DEBUG */

   make_state_from_rep(rep, &stateNow);

#ifdef DEBUG
    if (is_out_grid(stateNow.T.x, stateNow.T.y, stateNow.T.z)) {
       (void)fprintf(logFile,"eval.cc/stateNow.T is outside grid!\n");
    }
#endif /* DEBUG */

#ifdef DEBUG
    (void)fprintf(logFile,"eval.cc/Converting state to coordinates...\n");
#endif /* DEBUG */
 
   // Ligand could be inside or could still be outside, check all the atoms...
   cnv_state_to_coords(stateNow, vt, tlist, stateNow.ntor, crdpdb, crd, natom);

#ifdef DEBUG
(void)fprintf(logFile,"eval.cc/Checking to see if all coordinates are inside grid...\n");
#endif /* DEBUG */

   //  Check to see if crd is valid
   for (i=0; (i<natom)&&(!B_outside); i++) {
      B_outside = is_out_grid(crd[i][0], crd[i][1], crd[i][2]);
   } // i

   if (!B_outside) {

#ifdef DEBUG
(void)fprintf(logFile,"eval.cc/All coordinates are inside grid...\n");
#endif /* DEBUG */

/*---------------------------------------------------------------------------
* I have removed this call to save a stack-push, and inlined the code 
* instead...
*                 -- Garrett
*
*    energy = evaluate_energy(crd, charge, type, natom, map, inv_spacing, 
*                              xlo, ylo, zlo, nonbondlist, 
*                              e_internal, Nnb, B_calcIntElec, q1q2, 
*                              B_isGaussTorCon, B_isTorConstrained, stateNow, 
*                              B_ShowTorE, US_TorE, US_torProfile);
* --------------------------------------------------------------------------*/

        energy = quicktrilinterp( crd, charge, type, natom, map, 
                                  inv_spacing, xlo, ylo, zlo)
                 + eintcal( nonbondlist, e_internal, crd, type, Nnb, 
                            B_calcIntElec, q1q2);
        /*
        energy = trilinterp( crd, charge, type, natom, map, 
                                  inv_spacing, eval_elec, eval_emap, xlo, ylo, zlo)
                 + eintcal( nonbondlist, e_internal, crd, type, Nnb, 
                            B_calcIntElec, q1q2);
                            */
     
        if (B_isGaussTorCon) {
            for (I_tor = 0; I_tor <= stateNow.ntor; I_tor++) {
                if (B_isTorConstrained[I_tor] == 1) {
                    indx = Rad2Div( WrpModRad(stateNow.tor[I_tor]) );
                    if (B_ShowTorE) {
                        energy += (double)(US_TorE[I_tor] = US_torProfile[I_tor][indx]);
                    } else {
                        energy += (double)US_torProfile[I_tor][indx];
                    }
                }
            } // I_tor
        }/*if*/
   } else {
        /*
         * This confuses the GA and GA-LS, because there is no gradient
         * information when all outside conformations are given the same
         * energy.
         *
         * energy = BIG_ENERGY;  / / A really big number defined in autocomm.h
         */
        /*
         * Instead...
         *
         * Penalise atoms outside grid based on the square of the 
         * distance from centre of grid map, otherwise use the normal 
         * trilinear interpolation.
         */
        energy = outsidetrilinterp( crd, charge, type, natom, map,
                                    inv_spacing, // eval_elec, eval_emap, 
                                    xlo, ylo, zlo,
                                    xhi, yhi, zhi,  xcen, ycen, zcen )
                 + eintcal( nonbondlist, e_internal, crd, type, Nnb,
                            B_calcIntElec, q1q2);
        if (B_isGaussTorCon) {
            for (I_tor = 0; I_tor <= stateNow.ntor; I_tor++) {
                if (B_isTorConstrained[I_tor] == 1) {
                    indx = Rad2Div( WrpModRad(stateNow.tor[I_tor]) );
                    if (B_ShowTorE) {
                        energy += (double)(US_TorE[I_tor] = US_torProfile[I_tor][indx
]);
                    } else {
                        energy += (double)US_torProfile[I_tor][indx];
                    }
                }
            } // I_tor
        } // if
   }

   num_evals++;

   if (!finite(energy)) {
      (void)fprintf( logFile, "eval.cc:  ERROR!  energy is infinite!\n\n");
      for (i=0; i<natom; i++) {
          // (void)fprintf( logFile, "ATOM  %5d  C   INF     1    %8.3f%8.3f%8.3f %+8.2f %+6.2f  %+6.3f\n", i+1, crd[i][X], crd[i][Y], crd[i][Z], eval_emap[i], eval_elec[i], charge[i]); 
          (void)fprintf( logFile, "ATOM  %5d  C   INF     1    %8.3f%8.3f%8.3f %+8.2f %+6.2f  %+6.3f\n", i+1, crd[i][X], crd[i][Y], crd[i][Z], 0.0, 0.0, charge[i]); 
      } // i
   }
   if (ISNAN(energy)) {
      (void)fprintf( logFile, "eval.cc:  ERROR!  energy is not a number!\n\n");
      for (i=0; i<natom; i++) {
          // (void)fprintf( logFile, "ATOM  %5d  C   NaN     1    %8.3f%8.3f%8.3f %+8.2f %+6.2f  %+6.3f\n", i+1, crd[i][X], crd[i][Y], crd[i][Z], eval_emap[i], eval_elec[i], charge[i]); 
          (void)fprintf( logFile, "ATOM  %5d  C   NaN     1    %8.3f%8.3f%8.3f %+8.2f %+6.2f  %+6.3f\n", i+1, crd[i][X], crd[i][Y], crd[i][Z], 0.0, 0.0, charge[i]); 
      } // i
   }
   return(energy);
}

int Eval::write(FILE *out_file, Representation **rep)
{
    int i, retval;
    //char rec14[14];

#ifdef DEBUG
    (void)fprintf(logFile,"eval.cc/int Eval::write(FILE *out_file, Representation **rep)\n");
#endif /*DEBUG*/

    make_state_from_rep(rep, &stateNow);
    cnv_state_to_coords(stateNow, vt, tlist, stateNow.ntor, crdpdb, crd, natom);
    for (i=0; i<natom; i++) {
        // strncpy( rec14, &atomstuff[i][13], (size_t)13);
        // rec14[13]='\0';
        //strncpy(rec14, "C   RES     1\0", (size_t)14);
        //retval = fprintf( out_file, "ATOM  %5d  %13s    %8.3f%8.3f%8.3f %+8.2f %+6.2f  %+6.3f\n", i+1, rec14, crd[i][X], crd[i][Y], crd[i][Z], 0., 0., charge[i]); 
        retval = fprintf( out_file, "ATOM  %5d  C   RES     1    %8.3f%8.3f%8.3f %+8.2f %+6.2f  %+6.3f\n", i+1, crd[i][X], crd[i][Y], crd[i][Z], 0., 0., charge[i]); 
    } // i
    return retval;
}