main.cc-new

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

/*____________________________________________________________________________*/

    case DPF_TRJFRQ:
        /*
        **  trjfrq
        **  Trajectory frequency,
        */
        (void) sscanf( line, "%*s %d", &trj_freq);
        B_write_trj = (trj_freq > 0);
        pr( logFile, UnderLine );
        pr( logFile, "\t\tTRAJECTORY INFORMATION\n" );
        pr( logFile, "\t\t______________________\n\n\n" );
        if (B_write_trj) {
	    pr( logFile, "Output frequency for trajectory frames =\tevery %d step%s\n", trj_freq, (trj_freq > 1)?"s.":"." );
        } else {
	    pr( logFile, "No trajectory of states will be written.\n\n" );
	    pr( logFile, "Subsequent \"trjbeg\", \"trjend\", \"trjout\" and \"trjsel\" parameters will be ignored.\n\n" );
        }
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_TRJBEG:
        /*
        **  trjbeg
        **  Trajectory begin cycle,
        */
        (void) sscanf( line, "%*s %d", &trj_begin_cyc ); 
        pr( logFile, "Begin outputting trajectory of states at cycle:\t%d\n", trj_begin_cyc );
        if (trj_begin_cyc < 0) {
	    trj_begin_cyc = 0;
        } else if (trj_begin_cyc > ncycles) {
	    trj_begin_cyc = trj_end_cyc = ncycles;
        }
        --trj_begin_cyc;
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_TRJEND:
        /*
        **  trjend
        **  Trajectory end cycle,
        */
        (void) sscanf( line, "%*s %d", &trj_end_cyc ); 
        pr( logFile, "Cease outputting trajectory of states at cycle:\t%d\n", trj_end_cyc );
        if (trj_end_cyc > ncycles) {
	    trj_end_cyc = ncycles;
        } else if (trj_end_cyc < 0) {
	    trj_end_cyc = 1;
        }
        --trj_end_cyc;
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_TRJOUT:
        /*
        **  trjout
        **  Trajectory file,
        */
        (void) sscanf( line, "%*s %s", FN_trj );
        pr( logFile, "\nWrite trajectory of state variables to file: \"%s\"\n", FN_trj);
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_TRJSEL:
        /*
        **  trjsel
        **  Trajectory select,
        */
        (void) sscanf( line, "%*s %1s", &out_acc_rej );
        B_acconly = (out_acc_rej == 'A');
        B_either  = (out_acc_rej == 'E');
        if (B_acconly) {
	    pr( logFile, "Output *accepted* states only.\n" );
        } else if (B_either) {
	    pr( logFile, "Output *either* accepted or rejected states.\n" );
        } else {
	    pr( logFile, "WARNING: Missing or unknown accepted/rejected output flag.\n" );
        }
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_EXTNRG:
        /*
        **  extnrg
        **  Wall Energy,
        */
        (void) sscanf( line, "%*s %f", &WallEnergy );
        pr( logFile, "External grid energy (beyond grid map walls) = %.2f\n\n", WallEnergy );
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_CLUSTER:
        /*
        **  cluster
        **  Cluster mode,
        */
        (void) sscanf( line, "%*s %s", FN_clus );
        B_cluster_mode = TRUE;
        pr( logFile, "Cluster mode is now set.\n\n" );
	clmode( atm_typ_str, num_atm_maps, clus_rms_tol,
	  hostnm, jobStart, tms_jobStart, 
	  B_write_all_clusmem, FN_clus, crdpdb, lig_center, 
	  B_symmetry_flag, FN_rms_ref_crds );
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_CLUSALL:
        /*
        ** write_all_clusmem
        ** Write all cluster members...
        */
        B_write_all_clusmem = TRUE;
        pr( logFile, "All members of each cluster will be written out after the clustering histogram.\n(This is instead of outputting just the lowest energy member in each.)\n\n" );
        break; 

/*____________________________________________________________________________*/

    case DPF_RMSNOSYM:
        /*
        **  rmsnosym
        **  Calculate RMS values in the normal way,
        **  ignoring any atom-type equivalences...
        */
        B_symmetry_flag = FALSE;
        pr( logFile, "Symmetry will be ignored in RMS calculations.\n\n" );
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_SCHEDLIN:
        /*
        **  linear_schedule
        **  Use a linear (arithmetic) temperature
        **  reduction schedule.  This is necessary for
        **  more accurate entropy estimations...
        */
        B_linear_schedule = TRUE;
        pr( logFile, "A linear temperature reduction schedule will be used...\n\n" );
        if (ncycles == -1) {
	    pr( logFile, "\nPlease specify the number of cycles first!\n\n" );
        } else {
	    RTreduc = RT0 / ncycles;
	    pr( logFile, "Annealing temperature will be reduced by %.3f cal mol per cycle.\n\n", RTreduc );
        }
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_INTELEC:
        /*
        **  intelec
        **  Calculate internal electrostatic energies...
        */
        B_calcIntElec = TRUE;
        pr( logFile, "Internal electrostatic energies will be calculated.\n\n");
        if (!B_haveCharges) {
	    pr( logFile, "%s: WARNING! No partial atomic charges have been supplied yet.\n\n",programname);
	} else {
	    pr(logFile,"Calculating the product of the partial atomic charges q1*q2 for all %d non-bonded pairs...\n\n\n",Nnb);
	    pr(logFile,"Non-bonded                           Scaled\n");
	    pr(logFile,"   Pair     Atom1-Atom2    q1*q2      q1*q2\n");
	    pr(logFile,"__________  ___________  _________  _________\n");
	    for (i = 0;  i < Nnb;  i++) {
		atm1 = nonbondlist[i][ATM1];
		atm2 = nonbondlist[i][ATM2];
		q1q2[i] = charge[atm1] * charge[atm2];
		pr(logFile,"   %4d     %5d-%-5d    %5.2f",i+1,atm1+1,atm2+1,q1q2[i]);
		q1q2[i] *= ELECSCALE;
		pr(logFile,"    %5.2f\n",q1q2[i]);
	    }
	    pr(logFile,"\n");
	}
        (void) fflush(logFile);
        break; 

/*____________________________________________________________________________*/

    case DPF_SEED:
        /*
        **  seed
        **  Set the random-number gerator's seed value,
        */
        retval = (int)sscanf( line, "%*s %s %s", param[0], param[1]);
	timeSeedIsSet[0] = 'F';
	timeSeedIsSet[1] = 'F';
        pr(logFile, "%d seed%c found.\n", retval, ((retval==1)? ' ' : 's'));
        for (j=0; j<retval; j++) {
            for (i=0; i<(int)strlen(param[j]); i++) {
                param[j][i] = (char)tolower( (int)param[j][i] );
            }
            pr(logFile, "argument \"%s\" found\n", param[j]);
        }
        if ((retval==2) || (retval==1)) {
            for (i=0; i<retval ; i++ ) {
                if (equal(param[i], "tim", 3)) {
		    timeSeedIsSet[i] = 'T';
                    seed[i] = (FourByteLong)time( &time_seed );
                    seed_random(seed[i]);
	            pr(logFile,"Random number generator was seeded with the current time, value = %ld\n",seed[i]);
                } else if (equal(param[i], "pid", 3)) {
		    timeSeedIsSet[i] = 'F';
                    seed[i] = getpid();
                    seed_random(seed[i]);
	            pr(logFile,"Random number generator was seeded with the process ID, value   = %ld\n",seed[i]);
                } else {
		    timeSeedIsSet[i] = 'F';
                    seed[i] = atol(param[i]);
                    seed_random(seed[i]);
	            pr(logFile,"Random number generator was seeded with the user-specified value  %ld\n",seed[i]);
                }
            }/*i*/
            pr(logFile, "\n");
        } else {
            pr(logFile, "Error encountered reading seeds!\n");
        }

        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_WATCH:
        /*
        **  watch
        **  for watching a job's progress PDBQ file in AVS,
        */
        (void) sscanf( line, "%*s %s", FN_watch);
        if (B_write_trj) {
            pr(logFile,"\nAutoDock will create the watch-file \"%s\", for real-time monitoring of runs.\n\n", FN_watch);
            pr(logFile,"\nThe watch-file will be updated every %d moves, in accordance with the trajectory parameters..\n\n", trj_freq);
            B_watch = TRUE;
        } else {
	    pr(logFile,"\nYou must set \"trjfrq\" to be greater than zero. No watch-file will be created.\n\n");
            B_watch = FALSE;
        }
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_GAUSSTORCON: 
    case DPF_HARDTORCON: 
        /*
        ** "gausstorcon" Add Gaussian torsion contraints,
        ** "hardtorcon"  Add Hard torsion contraints,
        */
        (void) sscanf( line, "%*s %d %f %f", &I_tor, &F_torPref, &F_torHWdth);
        if (I_tor <= 0) {
	    pr( logFile, "\nTorsion IDs less than 1 (%d) are not allowed!\n\n", I_tor);
        } else if (I_tor > ntor) {
	    pr( logFile, "\nRequested torsion ID (%d) is larger than the number of torsions found (%d)!\n\n", I_tor, ntor);
	} else { /* torsion-ID accepted */
	    --I_tor;    /* Because humans start at 1, and C at 0... */
	    if ( B_isTorConstrained[I_tor] == 0 ) {

		if (dpf_keyword ==  DPF_GAUSSTORCON) {
		    B_isGaussTorCon = TRUE;
		    B_isTorConstrained[I_tor] = 1;
		    /*
		    ** Initialize... Torsion Energy Profile...
		    ** Set energies at every torsion division 
		    ** to the user-defined (maximum) barrier energy,
		    */
		    for (US_tD = 0;  US_tD < NTORDIVS;  US_tD++) {
			US_torProfile[I_tor][US_tD] = US_torBarrier;
		    }
		} else {
		    /*
		    ** DPF_HARDTORCON
		    */
		    B_isTorConstrained[I_tor] = 2;
		}
	    }
	    if (dpf_keyword ==  DPF_GAUSSTORCON) {
		(void) strcpy( S_contype, " half-" );
	    } else {
		(void) strcpy( S_contype, " " );
	    }
		/*
	    ** F_torPref ranges from -180. to +180. degrees...
	    */
	    F_torPref = Wrp(ModDeg(F_torPref));
	    if (F_torHWdth < 0.) {
		pr(logFile,"\nI'm sorry, negative%swidths (%.1f) are not allowed. I will use the default (%.1f) instead.\n\n", S_contype, F_torHWdth, DEFHWDTH);
		F_torHWdth = DEFHWDTH;
	    } else if (F_torHWdth > 90.) {
		pr(logFile,"\nI'm sorry, your requested%swidth (%.1f) is too large. I will use the default (%.1f) instead.\n\n", S_contype, F_torHWdth, DEFHWDTH);
		F_torHWdth = DEFHWDTH;
	    }
	    pr( logFile, "For torsion %d, Adding a constrained-torsion zone centered on %.1f degrees;\n%swidth = %.1f degrees.\n\n", 1+I_tor, F_torPref, S_contype, F_torHWdth);

	    if (dpf_keyword == DPF_GAUSSTORCON) {
		/*
		** Calculate the torsion energy profile;
		** combine this with previous profile without
		** losing any information.
		*/
		for (F_A = F_A_from;  F_A <= F_A_to;  F_A += F_W) {
		    F_Aova = (F_A - F_torPref) / F_torHWdth;
		    US_energy = (unsigned short) (((float)US_torBarrier) * (1. - exp(F_lnH * F_Aova*F_Aova)));
		    /*
		    ** if F_A(<-180.or>180), wrap to -180to180,
		    */
		    F_tor = Wrp(ModDeg(F_A));
		    /*
		    ** Convert from F_tor to US_tD
		    */
		    US_tD = (unsigned short) ((F_tor - F_hW + 180.)/F_W);
		    US_torProfile[I_tor][US_tD] = min(US_energy,US_torProfile[I_tor][US_tD]);
		}/* for F_A */
		/*
		** Ensure that the lowest point(s) in the profile are
		** zero...
		*/
		US_min = TORBARMAX;
		for (US_tD = 0;  US_tD < NTORDIVS;  US_tD++) {
		    US_min = min(US_min,US_torProfile[I_tor][US_tD]);
		}
		for (US_tD = 0;  US_tD < NTORDIVS;  US_tD++) {
		    US_torProfile[I_tor][US_tD] -= US_min;
		}
	    } else { /*DPF_HARDTORCON*/

		iCon = N_con[I_tor] + 1;
		if (iCon < MAX_TOR_CON) {
		    F_TorConRange[I_tor][N_con[I_tor]][LOWER] = F_torPref - 0.5* F_torHWdth;
		    F_TorConRange[I_tor][N_con[I_tor]][UPPER] = F_torPref + 0.5* F_torHWdth;
		    N_con[I_tor] = iCon;
		} else {
		    pr(logFile,"\n\n I'm sorry, you can only have %d (=MAX_TOR_CON) torsion constraints.\nIf you need more, change the \"#define MAX_TOR_CON\" line in \"constants.h\".\n\n",MAX_TOR_CON);
		}/* Still room to add another constraint. */
	    } /*DPF_HARDTORCON*/
	}/* torsion-ID accepted */
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_BARRIER:
        /*
        **  barrier
        **  Define torsion-barrier energy...
        */
        (void) sscanf( line, "%*s %d", &I_torBarrier);
        US_torBarrier = (unsigned short)I_torBarrier;
        US_torBarrier = min(US_torBarrier, TORBARMAX);
        pr(logFile,"\nTorsion barrier energy is set to %uhd\n\n", US_torBarrier);
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_SHOWTORPEN:
        /*
        **  showtorpen
        **  Show torsion's penalty energy.
        */
        B_ShowTorE = TRUE;
        pr(logFile,"\nConstrained torsion penalty energies will be stored during docking, and output after each run\n\n");
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_E0MAX:
        /*
        **  e0max
        **  Set maximum initial energy,