main.cc-new

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

/*____________________________________________________________________________*/
/*
** Initialize torsion arrays and constants.
*/

(void) strcpy( torfmt, "%*s\0" ); /* len(torfmt) is 4 chars */

for (j = 0;  j < MAX_ATOMS;  j++ ) {
    Nnbonds[j] = type[j] = 0;
} 

for (i = 0; i  < MAX_TORS;  i++ ) {
    for (j = 0;  j < MAX_ATOMS;  j++ ) {
	tlist[i][j] = 0;
    }
}

for (i = 0; i  < MAX_TORS;  i++ ) {
    B_isTorConstrained[i] = 0;
    US_torProfile[i][0] = 0;
    N_con[i] = 0;
}

initialiseState( &sInit );
initialiseState( &(mol.S) );

F_W      =  360. / NTORDIVS;
F_hW     =  F_W  / 2.;
F_A_from = -360. + F_hW;
F_A_to   =  360. + F_hW;

for (k = 0; k < MAX_RUNS; k++) {
    for (i = 0; i  < MAX_TORS;  i++ ) {
        sHist[k].tor[i] = 0.;
    }
}

for (i = 0; i < MAX_TORS;  i++ ) {
    if ( (ltorfmt += 3) > LINE_LEN ) {
	prStr( error_message, "ERROR: MAX_TORS = %d torsions declared in \"constants.h\";\n\t LINE_LEN = %d, Therefore you must change \"LINE_LEN\" to exceed %d...\n", MAX_TORS, LINE_LEN, 4+3*MAX_TORS );
	stop( error_message );
	exit( -1 );
    } else {
	(void) strcat( torfmt, " %f\0" );  /* add on 3 chars  for each new torsion... */
    }
} /* len(torfmt) is 4+3*MAX_TORS chars */

for (j = 0; j < MAX_NONBONDS; j++) {
    nonbondlist[j][ATM1] = nonbondlist[j][ATM2] = 0;
}

for (j = 0; j < MAX_RUNS; j++) {
    // isort[j] = j;
    econf[j] = 0.;
}

B_constrain_dist = B_haveCharges = FALSE;
ntor1 = ntor = atomC1 = atomC2 = 0;
sqlower = squpper = 0.;

timeSeedIsSet[0] = 'F';
timeSeedIsSet[1] = 'F';

if (clktck == 0) {        /* fetch clock ticks per second first time */
    if ( (clktck = sysconf(_SC_CLK_TCK)) < (FourByteLong)0L) {
	stop("\"sysconf(_SC_CLK_TCK)\" command failed in \"main.c\"\n");
	exit( -1 );
    } else {
	idct = (float)1. / (float)clktck;
    }
}

(void) strcpy(FN_rms_ref_crds,"unspecified filename\0");

/*____________________________________________________________________________*/
/*
** log(x): compute the natural (base e) logarithm of x,
*/

F_lnH = ((float)log(0.5));

/*____________________________________________________________________________*/
/*
** Output banner...
*/

banner( version );

/*____________________________________________________________________________*/
/*
** Print the time and date when the file was created...
*/

pr( logFile, "This file was created at:\t\t\t" );
printdate( logFile, 1 );

(void) strcpy(hostnm, "unknown host\0");

if (gethostname( hostnm, MAX_CHARS ) == 0) {
    pr( logFile, "                   using:\t\t\t\"%s\"\n", hostnm );
}

pr( logFile, "\nNOTE: \"rus\" stands for:\n\n      r = Real, wall-clock or elapsed time;\n      u = User or cpu-usage time;\n      s = System time\n\nAll timings are in seconds, unless otherwise stated.\n\n\n" );

/*____________________________________________________________________________*/
/*
** (Note: "dock_param_fn" set in "setflags.c"...)
*/
pr( logFile, "Parameter file used for this docking:\t\t%s\n\n", dock_param_fn );


/*____________________________________________________________________________*/
/*
** Start reading in the parameter/run-control file,
*/

while( fgets(line, LINE_LEN, parFile) != NULL ) { /* PARSING-DPF parFile */

    dpf_keyword = parse_dpf_line( line );

    switch( dpf_keyword ) {
    case -1:
        pr( logFile, "DPF> %s", line );
        prStr( error_message,"%s: WARNING: Unrecognized keyword in docking parameter file.\n", programname );
        pr_2x( logFile, stderr, error_message );
        continue;
        /* break; */

    case DPF_NULL:
    case DPF_COMMENT:
        pr( logFile, "DPF> %s", line );
        (void) fflush(logFile);
        break;

    default:
        pr( logFile, "\n\nDPF> %s\n", line );
        indcom = strindex( line, "#" );
        if (indcom != -1) {
	    line[ indcom ] = '\0'; /* Truncate "line" at the comment */
        }
        (void) fflush(logFile);
        break;
    } /* switch */

    switch( dpf_keyword ) {

/*____________________________________________________________________________*/

    case DPF_NULL:
    case DPF_COMMENT:
        break;

/*____________________________________________________________________________*/
    case DPF_TYPES:
        /*
        ** types
        ** ATOM_TYPE_NAMES
        */
        (void) strncpy( atm_typ_str, "????????", (size_t)ATOM_MAPS);
        dpftypes( &Htype,&num_all_maps,&num_atm_maps,atm_typ_str,line );
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_FLD:
        /*
        ** fld
        ** GRID_DATA_FILE
        ** Read the (AVS-format) grid data file, .fld
        */
        readfield( &inv_spacing, &spacing, FN_gdfld, 
            FN_gpf, gridpts1, gridpts, &xhi,&yhi,&zhi, 
            jobStart, line, &xlo,&ylo,&zlo, FN_receptor, 
            map_center, tms_jobStart );
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/
 
    case DPF_MAP:
        /*
        ** map
        ** ATOMIC AFFINITY or ELECTROSTATIC GRID MAP
        ** Read in active site grid map...
        */
        B_charMap = FALSE;
        readmap( &B_havemap, &imap, &num_atm_maps, &spacing, 
            atm_typ_str, FN_gdfld, gridpts1, gridpts, 
            jobStart, line, FN_receptor, map, map_center, 
            mapmax, mapmin, tms_jobStart, B_charMap );
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/
 
    case DPF_CHARMAP:
        /*
        ** charmap
        ** ATOMIC AFFINITY or ELECTROSTATIC GRID MAP
        ** Read in active site grid map...
        */
        B_charMap = TRUE;
        readmap( &B_havemap, &imap, &num_atm_maps, &spacing, 
            atm_typ_str, FN_gdfld, gridpts1, gridpts, 
            jobStart, line, FN_receptor, map, map_center, 
            mapmax, mapmin, tms_jobStart, B_charMap );
        (void) fflush(logFile);
        break;
/*____________________________________________________________________________*/

    case DPF_MOVE:
        /*
        ** move
        ** Movable ligand,
        */
        mol = readPDBQ( line,
             atm_typ_str, num_atm_maps,
             &natom, 
             crdpdb, charge, &B_haveCharges, type,
             pdbaname, FN_ligand, atomstuff, Htype,
             &B_constrain_dist, &atomC1, &atomC2, 
             &sqlower, &squpper,
             &ntor1, &ntor, tlist, vt, 
             &Nnb, Nnbonds, nonbondlist,
             jobStart, tms_jobStart, hostnm);

        sInit.ntor = mol.S.ntor;
        ++nmol;
        ++nlig;

        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_ABOUT:
        /*
        **  about
        **  Rotation center for current ligand,
        */
        (void) sscanf( line, "%*s %f %f %f", &lig_center[X], &lig_center[Y], &lig_center[Z]);
        pr( logFile, "Small molecule center of rotation =\t" );
        pr( logFile, "(%+.3f, %+.3f, %+.3f)\n\n", lig_center[X], lig_center[Y], lig_center[Z]);
        /*
        **  Center the ligand,
        */
        if ( nmol == 0 ) {
	    pr( logFile, "Must specify a ligand PDBQ file, using the \"move\" command.\n");
	} else {
	    pr( logFile, "Translating small molecule by:\t" );
	    pr( logFile, "(%+.3f, %+.3f, %+.3f)\n\n", -lig_center[X], -lig_center[Y], -lig_center[Z]);
	    /*
	    **  Zero-out on central point...
	    */
	    maxrad = -1.;
	    
	    for ( i=0; i<natom; i++ ) {
		r2sum=0.;
		for (XYZ = 0;  XYZ < SPACE;  XYZ++) {
		    c = crd[i][XYZ] = (crdpdb[i][XYZ] -= lig_center[XYZ]);
		    r2sum += c*c;
		} /* XYZ */
		maxrad = max(maxrad,sqrt(r2sum));
	    } /* i */
	    pr( logFile, "Furthest ligand atom from \"about\" center is %.3f Angstroms (maxrad).\n\n",maxrad);

	}
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_TRAN0:
        /* 
        **  tran0
        **  Initial_translation,
        */
        (void) sscanf( line, "%*s %s", param[0]);
        for (i=0; i<6; i++) {
	    param[0][i] = (char)tolower( (int)param[0][i] );
        }
        if (equal(param[0],"random",6)) {
	    B_RandomTran0 = TRUE;
	    mol.S.T.x = sInit.T.x = random_range( xlo, xhi );
	    mol.S.T.y = sInit.T.y = random_range( ylo, yhi );
	    mol.S.T.z = sInit.T.z = random_range( zlo, zhi );
        } else {
	    B_RandomTran0 = FALSE;
	    (void) sscanf( line,"%*s %lf %lf %lf", &(sInit.T.x), &(sInit.T.y), &(sInit.T.z));
	    mol.S.T.x = sInit.T.x;
	    mol.S.T.y = sInit.T.y;
	    mol.S.T.z = sInit.T.z;
	    pr( logFile, "Initial translation =\t\t\t(%.3f, %.3f, %.3f) Angstroms\n", sInit.T.x, sInit.T.y, sInit.T.z );
	    (void) fflush(logFile);
        }
        break;

/*____________________________________________________________________________*/

    case DPF_QUAT0:
        /*
        **  quat0
        **  Initial_quaternion,
        */
        (void) sscanf( line, "%*s %s", param[0]);
        for (i=0; i<6; i++) {
	    param[0][i] = (char)tolower( (int)param[0][i] );
        }
        if (equal(param[0],"random",6)) {
	    B_RandomQuat0 = TRUE;
	    sInit.Q.nx  = random_range( -1.,   1. );
	    sInit.Q.ny  = random_range( -1.,   1. );
	    sInit.Q.nz  = random_range( -1.,   1. );
	    sInit.Q.ang = random_range(  0., 360. );

	    pr( logFile, "Each run will begin with a new, random initial quaternion.\n");
        } else {
	    B_RandomQuat0 = FALSE;

        }
	(void) sscanf( line, "%*s %lf %lf %lf %lf", &sInit.Q.nx, &sInit.Q.ny, &sInit.Q.nz, &sInit.Q.ang);
	pr( logFile, "Initial quaternion,  q = [(x,y,z),w] =\t[ (%.3f, %.3f, %.3f), %.1f deg ],\n", sInit.Q.nx, sInit.Q.ny, sInit.Q.nz, sInit.Q.ang);
	mol.S.Q.nx  = sInit.Q.nx;
	mol.S.Q.ny  = sInit.Q.ny;
	mol.S.Q.nz  = sInit.Q.nz;
	mol.S.Q.ang = sInit.Q.ang;

	mol.S.Q.ang = sInit.Q.ang = Rad( sInit.Q.ang ); /*convert to radians*/
	mkUnitQuat( &sInit.Q );
	mkUnitQuat( &(mol.S.Q) );

	pr( logFile, "Quaternion Vector Normalized to:  %.3f %.3f %.3f\n", sInit.Q.nx, sInit.Q.ny, sInit.Q.nz );
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_NDIHE:
        /*
        **  ndihe
        **  Number of dihedral angles to be specified by "dihe0"
        */
        (void) sscanf( line, "%*s %d", &ndihed );
        if ( nmol == 0 ) {
	    pr( logFile, "Must specify a ligand PDBQ file, using the \"move\" command.\n");
	} else {
	    pr( logFile, "Number of torsions = %d\n", ndihed);
	    if ( ndihed != ntor ) {
		pr( logFile, "WARNING! You requested %d torsions, but I found %d in PDBQ-file specifications.\n", ndihed, ntor );
	    } /* if */
	}
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_DIHE0:
        /*
        **  dihe0
        **  Initial dihedral angles, input in degrees,
        */
        (void) sscanf( line, "%*s %s", param[0]);
        for (i=0; i<6; i++) {
	    param[0][i] = (char)tolower( (int)param[0][i] );
        }
	if (equal(param[0],"random",6)) {
	    B_RandomDihe0 = TRUE;
	    sInit.ntor = nval = ndihed;
	    for ( i=0; i<nval; i++ ) {
		sInit.tor[i] = random_range( -180., 180. );
	    }
	} else {
	    B_RandomDihe0 = FALSE;
	    retval = (int)sscanf( line, torfmt, TOR_ARG_LIST );
	    if (retval == 0) {
		pr( logFile, "Could not read any torsions!\n" );
	    } else if (retval == EOF) {
		pr( logFile, "End of file encountered while reading dihe0 line\n");
	    } else if (retval < ndihed) {
		pr( logFile, "WARNING!  Only %d initial torsion angles were detected on input line.\n",retval);
		pr( logFile, "I am sorry, you set 'ndihe', the number of dihedrals, to %d torsions.\n", ndihed);
	    } else {
		pr( logFile, "%d initial torsion angles were detected on input line.\n", retval );
	    }