getinitialstate.cc

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

/* getInitialState.cc */

#include <math.h>

#ifdef sgi
    #include <stdio.h>
    #include <sys/types.h>
    #include <sys/times.h>
    #include <sys/param.h>
    #include <time.h>
    #include "getInitialState.h"
#else
    extern "C"
    {
        #include <stdio.h>
        #include <sys/types.h>
        #include <sys/times.h>
        #include <sys/param.h>
        #include <time.h>
        #include "getInitialState.h"
    }
#endif


extern FILE *logFile;
extern char *programname;


void getInitialState(   

            float *Addr_e0total,
            float e0max,

            State *sInit, /* was qtn0[QUAT] and tor0[MAX_TORS] */
            State *sMinm, /* was qtnMin[QUAT] and torMin[MAX_TORS] */
            State *sLast, /* was qtnLast[QUAT] and torLast[MAX_TORS] */

            Boole B_RandomTran0,
            Boole B_RandomQuat0,
            Boole B_RandomDihe0,

            float charge[MAX_ATOMS],
            float q1q2[MAX_NONBONDS],
            float crd[MAX_ATOMS][SPACE],
            float crdpdb[MAX_ATOMS][SPACE],
            char  atomstuff[MAX_ATOMS][MAX_CHARS],
            float elec[MAX_ATOMS],
            float emap[MAX_ATOMS],
            float e_internal[NEINT][ATOM_MAPS][ATOM_MAPS],
            Boole B_calcIntElec,
            float xhi,
            float yhi,
            float zhi,
            float xlo,
            float ylo,
            float zlo,
            float inv_spacing,
            float map[MAX_GRID_PTS][MAX_GRID_PTS][MAX_GRID_PTS][MAX_MAPS],
            int   natom,
            int   Nnb,
            int   nonbondlist[MAX_NONBONDS][2],
            int   ntor,
            int   tlist[MAX_TORS][MAX_ATOMS],
            int   type[MAX_ATOMS],
            float vt[MAX_TORS][SPACE],
            int   irun1,
            int   outlev,
            int   MaxRetries,
            float torsFreeEnergy)

{
    float e0total = 0.;
    float e0inter = 0.;
    float e0intra = 0.;
    float e0min = BIG;
    int   retries = 0;
    register int i = 0;
    Clock  initStart;
    Clock  initEnd;
    struct tms tms_initStart;
    struct tms tms_initEnd;


    /* Store time at start of initialization process... */

    initStart = times( &tms_initStart );

    retries = 0;
    do {
        /*
        ** while e0total, initial energy, is too high,
        */

        /*
        ** Initialize all state variables...
        */
        if (B_RandomTran0) {
            sInit->T.x = random_range( xlo, xhi );
            sInit->T.y = random_range( ylo, yhi );
            sInit->T.z = random_range( zlo, zhi );
            if (outlev > 1) {
                pr( logFile, "Random initial translation,  tran0 %.3f %.3f %.3f\n", sInit->T.x, sInit->T.y, sInit->T.z);
            }
        }/*if*/
        if (B_RandomQuat0) {
            sInit->Q.nx  = random_range( -1., 1. );
            sInit->Q.ny  = random_range( -1., 1. );
            sInit->Q.nz  = random_range( -1., 1. );
            sInit->Q.ang = Rad( random_range( 0., 360.) );/*convert to radians*/

            mkUnitQuat( &(sInit->Q) );

            if (outlev > 1) {
                pr( logFile, "Random initial quaternion,  quat0 %.3f %.3f %.3f %.1f\n", sInit->Q.nx, sInit->Q.ny, sInit->Q.nz, Deg( sInit->Q.ang ) );
            }
        }/*if*/
        if ( B_RandomDihe0 && (ntor > 0) ) {
            if (outlev > 1) {
                pr( logFile, "Random initial torsions, ndihe = %d\ndihe0 = ", ntor);
            }
            sInit->ntor = ntor;
            for (i=0; i<ntor; i++) {
                sInit->tor[i] = random_range(-180.,180.);
                if (outlev > 1) {
                    pr( logFile, "%7.2f ", sInit->tor[i] ); /*in degrees*/
                }
                sInit->tor[i] = Rad( sInit->tor[i] ); /*now in radians*/
            }
            if (outlev > 1) {
                pr( logFile, "\n");
            }
        }/*if*/

        copyState( sMinm, *sInit );
        copyState( sLast, *sInit );

/* _________________________________________________________________________
**
** Initialize the automated docking simulation,
** _________________________________________________________________________
*/
        initautodock( atomstuff, crd, crdpdb, xhi, yhi, zhi, xlo, ylo, zlo, 
            natom, ntor, sInit, tlist, vt, outlev);
        
        e0inter = trilinterp( crd, charge, type, natom, map, inv_spacing, elec, emap, xlo, ylo, zlo );
        e0intra = eintcal( nonbondlist, e_internal, crd, type, Nnb, B_calcIntElec, q1q2);
        e0total = e0inter + e0intra;

        if (e0total < e0min) {
            /*
            ** This energy is lower so update
            ** the initialization minimum-energy state variables,
            */
            e0min = e0total;
            copyState( sMinm, *sInit );
        }

        if ( (e0total > e0max) && ((!B_RandomTran0)||(!B_RandomQuat0)) ) {
            B_RandomTran0 = TRUE;
            B_RandomQuat0 = TRUE;
        }

        ++retries;
        if ((retries > 0) && (retries < MaxRetries) &&
            (e0total > e0max) && (outlev > 1)) {

            pr(logFile, "Initial total energy, e0total = %.3f, too high!\n", e0total);
            pr(logFile, "Number of attempts = %d (run %d)\n\n", retries, irun1);
            pr(logFile, "Will try again...\n");

        } else if (retries >= MaxRetries) {

            pr( logFile, "Sorry, too many retries (%d).  Continuing...\n\nWill use the state with the lowest energy found, %.2f\n\n", MaxRetries, e0min);
            e0total = e0min;
            copyState( sInit, *sMinm );

            break;
        }
        fflush( logFile );
    } while ( e0total > e0max );

    cnv_state_to_coords( *sInit, vt, tlist, ntor, crdpdb, crd, natom );

    e0inter = trilinterp( crd, charge, type, natom, map, inv_spacing, elec, emap, xlo, ylo, zlo );
    e0intra = eintcal( nonbondlist, e_internal, crd, type, Nnb, B_calcIntElec, q1q2);
    e0total = e0inter + e0intra;

    copyState( sMinm, *sInit );
    copyState( sLast, *sInit );

    prInitialState(e0inter, e0intra, torsFreeEnergy, natom, crd, atomstuff, type, emap, elec, charge);

    initEnd = times( &tms_initEnd );
    pr(logFile, "Number of initialization attempts = %d (run %d)\n", retries, irun1);
    pr(logFile, "Time spent initializing: (Real, CPU, System): ");
    timesys( initEnd - initStart, &tms_initStart, &tms_initEnd );

    *Addr_e0total = e0total;
    fflush( logFile );

}
/* EOF */