simanneal.cc

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

/* simanneal.cc */

#include <math.h>

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



extern FILE *logFile;
extern char *programname;


void simanneal( int   *Addr_nconf,
		int   Nnb,
		float WallEnergy,
		char  atomstuff[MAX_ATOMS][MAX_CHARS],
		float charge[MAX_ATOMS],
		Boole B_calcIntElec,
		float q1q2[MAX_NONBONDS],
		float crd[MAX_ATOMS][SPACE],
		float crdpdb[MAX_ATOMS][SPACE],
		char  FN_dpf[MAX_CHARS],
		float e_internal[NEINT][ATOM_MAPS][ATOM_MAPS],
		float econf[MAX_RUNS],
		Boole B_either,
		float elec[MAX_ATOMS],
		float emap[MAX_ATOMS],
		float xhi,
		float yhi,
		float zhi,
		int   NcycMax,
		float inv_spacing,
		int   irunmax,
		Clock jobStart,
		float xlo,
		float ylo,
		float zlo,
		float map[MAX_GRID_PTS][MAX_GRID_PTS][MAX_GRID_PTS][MAX_MAPS],
		int   naccmax,
		int   natom,
		int   nonbondlist[MAX_NONBONDS][2],
		int   nrejmax,
		int   ntor1,
		int   ntor,
		int   outlev,

		State sInit, /* tor0, qtn0 */
		State sHist[MAX_RUNS], /* was qtnHist, torHist */

		float qtwFac,
		Boole B_qtwReduc,
		float qtwStep0,
		Boole B_selectmin,
		char  FN_ligand[MAX_CHARS],
		float sml_center[SPACE],
		float RT0,
		Boole B_RTChange,
		float RTFac,
		struct tms tms_jobStart,
		int   tlist[MAX_TORS][MAX_ATOMS],
		float torFac,
		Boole B_torReduc,
		float torStep0,
		char  FN_trj[MAX_CHARS],
		int   trj_cyc_max,
		int   trj_cyc_min,
		int   trj_freq,
		float trnFac,
		Boole B_trnReduc,
		float trnStep0,
		int   type[MAX_ATOMS],
		float vt[MAX_TORS][SPACE],
		Boole B_writeTrj,
		Boole B_constrain,
		int   atomC1,
		int   atomC2,
		float sqlower,
		float squpper,
		Boole B_linear_schedule,
		float RTreduc,
		/*float maxrad,*/
		Boole B_watch,
		char  FN_watch[MAX_CHARS],
		Boole B_isGaussTorCon,
		unsigned short US_torProfile[MAX_TORS][NTORDIVS],
		Boole B_isTorConstrained[MAX_TORS],
		Boole B_ShowTorE,
		unsigned short US_TorE[MAX_TORS],
		float F_TorConRange[MAX_TORS][MAX_TOR_CON][2],
		int N_con[MAX_TORS],
		Boole B_RandomTran0,
		Boole B_RandomQuat0,
		Boole B_RandomDihe0,
		float e0max,
		float torsFreeEnergy,
		int   MaxRetries)

{
    char message[LINE_LEN];


    FILE *FP_trj;

    State sNow; /* qtnNow, torNow */
    State sChange; /* qtnChange, torChange */
    State sLast; /* qtnLast, torLast */
    State sMin; /* qtnMin, torMin */
    State sSave; /* qtnSave, torSave */

    float d[SPACE];
    float e = 0.;
    float e0 = 0.;
    float einter = 0.;
    float eintra = 0.;
    float eLast = 0.;
    float eMin = BIG_ENERGY;
    float etot = 0.0;
    float inv_RT = 0.;
    float qtwStep;
    float rsqC1C2;
    float RT = 616.;
    float torTmp;
    float torStep;
    float trnStep;
    /* ** float xloTrn; ** float xhiTrn; ** float yloTrn; ** float yhiTrn; ** float zloTrn; ** float zhiTrn; ** float lo[SPACE]; ** float trnStepHi; ** float qtwStepHi; ** float torStepHi; */

    Boole B_inRange = FALSE;
    Boole B_outside = FALSE;
    Boole B_within_constraint = TRUE;

    int count = 0;
    int Itor = 0;
    int icycle = -1;
    int icycle1 = -1;
    int irun = -1;
    int irun1 = -1;
    int indx;

    register int i = 0;
    register int nacc = 0;
    register int nAcc = 0;
    register int nAccProb = 0;
    register int nedge = 0;
    register int nrej = 0;
    register int ntot = 1;
    register int XYZ = 0;

    struct tms tms_cycStart;
    struct tms tms_cycEnd;
    struct tms tms_jobEnd;

    Clock cycStart;
    Clock cycEnd;
    Clock jobEnd;

    time_t time_seed;

    /* trnStepHi = HI_NRG_JUMP_FACTOR * trnStep; ** qtwStepHi = HI_NRG_JUMP_FACTOR * qtwStep; ** torStepHi = HI_NRG_JUMP_FACTOR * torStep; */
    /* lo[X] = xlo;  lo[Y] = ylo;  lo[Z] = zlo;*/
    /* xloTrn = xlo + maxrad; ** xhiTrn = xhi - maxrad; ** yloTrn = ylo + maxrad; ** yhiTrn = yhi - maxrad; ** zloTrn = zlo + maxrad; ** zhiTrn = zhi - maxrad; */

/* Open the trajectory file for writing, =====================================*/

    if ( B_writeTrj ) { 
        if ( (FP_trj = fopen(FN_trj, "w")) == NULL ) {
            prStr( message, "\n%s: can't create trajectory file %s\n", programname, FN_trj);
            pr_2x( stderr, logFile, message );
            prStr( message, "\n%s: Unsuccessful Completion.\n\n", programname);
            pr_2x( stderr, logFile, message );
            jobEnd = times( &tms_jobEnd );
            timesys( jobEnd - jobStart, &tms_jobStart, &tms_jobEnd );
            pr_2x( logFile, stderr, UnderLine );
            exit( -1 );
        }/*END PROGRAM*/
    }/*endif*/

/* Begin the automated docking simulation, ===================================*/

    pr( logFile, "\n\n\t\tBEGINNING MONTE CARLO SIMULATED ANNEALING\n");
    pr( logFile, "     \t\t_________________________________________\n\n\n\n" );


    for ( irun = 0;  irun < irunmax;  irun++ ) { /*===========================*/

        irun1 = 1 + irun;

	/*
	** Initialize random number generator with a time-dependent seed...
	*/
	time_seed = time( &time_seed );
	seed_random( time_seed );

	if (outlev > 0) {
	    pr(logFile, "\n\tINITIALIZING AUTOMATED DOCKING SIMULATION\n" );
	    pr(logFile, "\t_________________________________________\n\n" );
	    pr( logFile, "RUN %d...\n\n", irun1);
	    pr( logFile, "Time-dependent Seed:  %ld\n\n", time_seed);
	}

        if ( B_writeTrj ) {
            pr( FP_trj, "ntorsions %d\nrun %d\n", ntor, irun1 );
            fflush( FP_trj );
        }

	getInitialState( &e0, e0max,
			 &sInit, &sMin, &sLast, 
			 B_RandomTran0, B_RandomQuat0, B_RandomDihe0, 
			 charge, q1q2, crd, crdpdb, atomstuff,
			 elec, emap, e_internal, B_calcIntElec,
			 xhi, yhi, zhi, xlo, ylo, zlo,
			 inv_spacing, map, natom, Nnb, nonbondlist,
			 ntor, tlist, type, vt, irun1, outlev, MaxRetries,
			 torsFreeEnergy);

        RT = RT0;		/* Initialize the "annealing" temperature */
	if (RT <= APPROX_ZERO) { RT = 616.; }
	inv_RT = 1. / RT;

        eMin    = min(BIG_ENERGY, e0);
        eLast   = e0;
        trnStep = trnStep0;	/* translation*/
        qtwStep = qtwStep0;	/* quaternion angle, w*/
        torStep = torStep0;	/* torsion angles*/

        if (outlev > 0) {
	    pr( logFile, "\n\n\t\tBEGINNING SIMULATED ANNEALING");
	    pr( logFile, "\n\t\t_____________________________\n\n");
	    pr( logFile, "\n      \t      \tMinimum     Average     | Acc/    Accepted:    Rejected:     |          |  XYZ-Translation  |        Time:        \n");
	      pr( logFile, "Run:  \tCycle:\tEnergy:     Energy:     |   /Rej: Down:  Up:   Total: Edge:  |   RT:    |   of Min.Energy   |  Real, CPU, System  \n" );
	      pr( logFile, "______\t______\t___________ ___________ | ______ ______ ______ ______ ______ | ________ | _________________ | ____________________\n" );
/*	                                  12345678901 12345678901   123456 123456 123456 123456 123456   12345678   12345 12345 12345
**	                 "%D /%D\T%D /%D\T%+11.2F     %+11.2F       %6.2F  %6D    %6D    %6D    %6D      %8.1F     %5.2F %5.2F %5.2F   ",
**	                 IRUN1, IRUNMAX, ICYCLE1, ICYCLEMAX, EmIN, ETOT/NTOT, QTNmIN[x], QTNmIN[y], QTNmIN[z], (NREJ!=0) ? (FLOAT)NACC/NREJ : -999., NACC, NREJ, NEDGE, rt
*/
        }
        fflush(logFile);
/*____________________________________________________________________________*/

        for ( icycle = 0;  icycle < NcycMax;  icycle++ ) {

            cycStart = times( &tms_cycStart );

            icycle1 = icycle + (ntot = 1);
            B_inRange = (icycle >= trj_cyc_min) && (icycle <= trj_cyc_max);
            if ( B_writeTrj && B_inRange ) {
		pr( FP_trj, "cycle %d\ntemp %f\n", icycle1, RT );
		fflush(  FP_trj  );
            }
            nAcc = nAccProb = nacc = nrej = nedge = 0;
	    etot = eLast;

/*____________________________________________________________________________*/

            do {
		/*
		** Do one Monte Carlo step,
                ** while the number of accepted steps < naccmax 
		** and number of rejected steps < nrejmax...
		*/

		mkNewState( &sNow, &sLast, &sChange,
                    vt, tlist, ntor, crd, crdpdb, natom,
		    trnStep,
		    /*qtwStep,*/
		    torStep,
		    F_TorConRange, N_con);

		if (B_constrain) {
		    for (XYZ = 0;  XYZ < SPACE;  XYZ++) {
			d[XYZ] = crd[atomC1][XYZ] - crd[atomC2][XYZ];
		    }
		    rsqC1C2 = sqhypotenuse(d[X],  d[Y], d[Z]);
		    if (! (B_within_constraint = (rsqC1C2 > sqlower) && 
						 (rsqC1C2 < squpper))) {
			copyState( &sNow, sLast );
		    }
		}/*if B_constrain*/