main.cc-new

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

        */
        retval = sscanf( line, "%*s %f %d", &e0max, &MaxRetries );
	if (retval == 0) {
	    pr( logFile, "Could not read any arguments!\n" );
	} else if (retval == EOF) {
	    pr( logFile, "End of file encountered!\n");
	} else if (retval == 1) {
	    pr(logFile,"Using the default maximum number of retries for initialization, %d retries.\n\n", MaxRetries);
	} else if (retval == 2) {
	    pr(logFile,"Using user-specified maximum number of retries for initialization, %d retries.\n\n", MaxRetries);
	}
        if (e0max < 0.) {
	    e0max = 1000.;
        }
        pr(logFile,"\nIf the initial energy is greater than e0max, %.3f,\nthen a new, random initial state will be created.\n\n",e0max);
        (void) fflush(logFile);
        break;

/*____________________________________________________________________________*/

    case DPF_SIMANNEAL:
        /*
        ** simanneal
        */
        /*
        ** Calculate reduction factor based on initial and final step values,
        ** and number of cycles...
        */
	if (!command_mode) {
	    ncycm1 = ncycles - 1;
	    if (ncycm1 < 0) {
		ncycm1 = 1;
	    }
	    if (B_CalcTrnRF) {
		trnFac = RedFac(trnStep0, trnStepFinal, ncycm1);
		pr( logFile, "Calculated reduction factor for translations     = %-.3f /cycle\n", trnFac);
		B_trnReduc = (trnFac != 1.);       
	    }
	    if (B_CalcQtwRF) {
		qtwFac = RedFac(qtwStep0, qtwStepFinal, ncycm1);
		pr( logFile, "Calculated reduction factor for quaternion angle = %-.3f /cycle\n", qtwFac );
		B_qtwReduc = (qtwFac != 1.);
	    }
	    if (B_CalcTorRF) {
		torFac    = RedFac(torStep0, torStepFinal, ncycm1);
		pr( logFile, "Calculated reduction factor for torsion angles   = %-.3f /cycle\n", torFac );
		B_torReduc = (torFac != 1.);
	    }
	    pr(logFile, "\n");
	    /*
	    ** Number of ligands read in...
	    */
	    if (nlig > 0) {
		pr( logFile, "\nTotal number of ligands read in by the DPF \"move\" command = %d\n\n", nlig );
	    }
	    if (nres > 0) {
		pr( logFile, "\nTotal number of residues read in by the DPF \"flex\" command = %d\n\n", nres );
	    }
	    pr(logFile, "\n");
	    if (B_havenbp) {
		nbe( atm_typ_str, e_internal, num_atm_maps );
	    }
	    if (B_cluster_mode) {
		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 );
	    }
	    for (j = 0; j < MAX_RUNS; j++) {
		econf[j] = torsFreeEnergy;
	    }
	    /* ___________________________________________________________________
	    ** 
	    ** Begin the automated docking simulation,
	    ** ___________________________________________________________________
	    */
	    simanneal(&nconf, Nnb, WallEnergy, atomstuff, charge,B_calcIntElec,
		    q1q2,crd,crdpdb,dock_param_fn,e_internal,econf,B_either,
		    elec,emap,xhi,yhi,zhi,ncycles,inv_spacing,nruns,jobStart,
		    xlo,ylo,zlo,map, 
		    naccmax, natom, nonbondlist, nrejmax, ntor1, ntor, outlev, 
		    sInit, sHist,   qtwFac, B_qtwReduc, qtwStep0,
		    B_selectmin,FN_ligand,lig_center,RT0,B_tempChange,RTFac, 
		    tms_jobStart, tlist, torFac, B_torReduc, torStep0, 
		    FN_trj, trj_end_cyc, trj_begin_cyc, trj_freq, trnFac,  
		    B_trnReduc, trnStep0, type, vt, B_write_trj,
		    B_constrain_dist, atomC1, atomC2, sqlower, squpper, 
		    B_linear_schedule, RTreduc, 
		    /*maxrad,*/
		    B_watch, FN_watch, 
		    B_isGaussTorCon, US_torProfile, B_isTorConstrained,
		    B_ShowTorE, US_TorE, F_TorConRange, N_con, 
		    B_RandomTran0, B_RandomQuat0, B_RandomDihe0,
		    e0max, torsFreeEnergy, MaxRetries);
	    (void) fflush(logFile);
	} else {
	    (void)fprintf(logFile, "NOTE: Command mode has been set, so simulated annealing cannot be performed.\n\n");
	}
        break;

/*____________________________________________________________________________*/

    case DPF_SET_GA:

      if (GlobalSearchMethod != NULL) {
          pr(logFile, "Deleting the previous settings for the Genetic Algorithm.\n");
          delete GlobalSearchMethod;
      }

      pr(logFile, "Passing the current parameters to the Genetic Algorithm.\n");
      GlobalSearchMethod = new Genetic_Algorithm(e_mode, s_mode, c_mode, w_mode, elitism, c_rate, m_rate, window_size, num_generations);
      ((Genetic_Algorithm *)GlobalSearchMethod)->mutation_values(low, high, alpha, beta);
      ((Genetic_Algorithm *)GlobalSearchMethod)->initialize(pop_size, 7+sInit.ntor);

      (void) fflush(logFile);
      break;
/*____________________________________________________________________________*/

    case DPF_SET_SW1:

      if (LocalSearchMethod != NULL) {
          pr(logFile, "Deleting the previous settings for the local search Solis-Wets algorithm (SW1 object).\n");
          delete LocalSearchMethod;
      }

      pr(logFile, "Passing the current settings to the local search Solis-Wets algorithm (SW1 object).\n");
      LocalSearchMethod = new Solis_Wets1(7+sInit.ntor, max_its, max_succ, max_fail, rho, lb_rho, 2.0, 0.5, search_freq);

      (void) fflush(logFile);
      break;
/*____________________________________________________________________________*/

    case DPF_SET_PSW1:

      if (LocalSearchMethod != NULL) {
          pr(logFile, "Deleting the previous settings for the local search pseudo-Solis-Wets algorithm (pSW1 object).\n");
          delete LocalSearchMethod;
      }

      pr(logFile, "Passing the current settings to the local search pseudo-Solis-Wets algorithm (pSW1 object).\n");

      //  Allocate space for the variable rho's
      rho_ptr = new float[7+sInit.ntor];
      lb_rho_ptr = new float[7+sInit.ntor];

      //  Initialize the rho's corresponding to the translation
      for (j=0; j<3; j++) {
         rho_ptr[j] = trnStep0;
         lb_rho_ptr[j] = trnStepFinal;
      }

      //  Initialize the rho's corresponding to the quaterion
      for (; j<7; j++) {
         rho_ptr[j] = qtwStep0;
         lb_rho_ptr[j] = qtwStepFinal;
      }

      //  Initialize the rho's corresponding to the torsions
      for (; j<7+sInit.ntor; j++) {
         rho_ptr[j] = torStep0;
         lb_rho_ptr[j] = torStepFinal;
      }

      LocalSearchMethod = new Pseudo_Solis_Wets1(7+sInit.ntor, max_its, max_succ, max_fail, 2.0, 0.5, search_freq, rho_ptr, lb_rho_ptr);

      (void) fflush(logFile);
      break;

/*____________________________________________________________________________*/

    case DPF_GALS:
	(void) fflush( logFile );
        /*
        ** Genetic Algorithm-Local search
        */
	if (!command_mode) {
	    (void) sscanf( line, "%*s %d", &nruns );

	    if ( nruns > MAX_RUNS ) {

		prStr( error_message, "ERROR: %d runs requested, but only dimensioned for %d.\nChange \"MAX_RUNS\" in \"constants.h\".", nruns, MAX_RUNS);
		stop( error_message );
		exit( -1 );

	    } else if ((GlobalSearchMethod==NULL)||(LocalSearchMethod==NULL)) {

		prStr(error_message, "ERROR:  Need to use \"set_ga\" to allocate both Global Optimization object AND Local Optimization object.\n");
		stop(error_message);
		exit(-1);

	    }

	    pr( logFile, "Number of requested Hybrid Global-Local Search dockings = %d run%c\n", nruns, (nruns > 1)?'s':' ');

	    evaluate.setup(crd, charge, type, natom, map, inv_spacing, 
	      elec, emap,
	      xlo, xhi, ylo, yhi, zlo, zhi, nonbondlist, e_internal, Nnb, 
	      B_calcIntElec, q1q2, B_isGaussTorCon, B_isTorConstrained,
	      B_ShowTorE, US_TorE, US_torProfile, vt, tlist, crdpdb, sInit, mol);

	    for (j=0; j<nruns; j++) {
		j1 = j + 1;

		(void) fprintf( logFile, "\n\n\tBEGINNING GENETIC ALGORITHM-LOCAL SEARCH DOCKING\n");
		(void) fflush( logFile );
		pr( logFile, "\nRun:\t%d / %d\n", j1, nruns );

		// Update time-based RNG seeds...
		if (timeSeedIsSet[0] == 'T') {
                    pr(logFile, "Updating First Time-dependent Seed Now.\n");
		    seed[0] = (FourByteLong)time( &time_seed );
		}
		if (timeSeedIsSet[1] == 'T') {
                    pr(logFile, "Updating Second Time-dependent Seed Now.\n");
		    seed[1] = (FourByteLong)time( &time_seed );
		}
		setall(seed[0], seed[1]);
		initgn(-1); // Reinitializes the state of the current random number generator
		pr(logFile, "Seeds:\t%ld %ld\nDate:\t", seed[0], seed[1]);
		printdate( logFile, 2 );
		(void) fflush( logFile );

		gaStart = times( &tms_gaStart );

		//  Can get rid of the following line
		((Genetic_Algorithm *)GlobalSearchMethod)->initialize(pop_size, 7+sInit.ntor);

		sHist[nconf] = call_glss( GlobalSearchMethod, LocalSearchMethod, sInit, 
					  num_evals, pop_size, xlo, xhi, 
					  ylo, yhi, zlo, zhi, outlev, &mol);

		gaEnd = times( &tms_gaEnd );
		pr( logFile, "\nRun completed;  time taken for this run:\n");
		timesyshms( gaEnd - gaStart, &tms_gaStart, &tms_gaEnd );
		pr( logFile, "\n");
		printdate( logFile, 1 );
		(void) fflush( logFile );

		pr(logFile, "Total number of Energy Evaluations: %lu\n", evaluate.evals() );
		pr(logFile, "Total number of Generations:        %u\n", ((Genetic_Algorithm *)GlobalSearchMethod)->num_generations());
 
		pr( logFile, UnderLine );
		pr( logFile, "\n\n\tFINAL GENETIC ALGORITHM-LOCAL SEARCH DOCKED STATE\n" );
		pr( logFile,     "\t_________________________________________________\n\n\n" );

		writePDBQState( j, FN_ligand, dock_param_fn, lig_center, 
		    &(sHist[nconf]), ntor, &eintra, &einter, natom, atomstuff, 
		    crd, emap, elec, charge, torsFreeEnergy,
		    vt, tlist, crdpdb, nonbondlist, e_internal,
		    type, Nnb, B_calcIntElec, q1q2,
		    map, inv_spacing, xlo, ylo, zlo, xhi, yhi, zhi);
		econf[nconf] = eintra + einter; // new2

		++nconf;
	    } // Next run
	    (void) fflush(logFile);
	} else {
	    (void)fprintf(logFile, "NOTE: Command mode has been set.  Sorry, genetic algorithm-local search cannot be performed.\n\n");
	}
	break;

/*____________________________________________________________________________*/

    case DPF_LS:
       (void) sscanf(line, "%*s %d", &nruns);

	if (!command_mode) {
	    if ( nruns > MAX_RUNS ) {

	       prStr( error_message, "ERROR: %d runs requested, but only dimensioned for %d.\nChange \"MAX_RUNS\" in \"constants.h\".", nruns, MAX_RUNS);
	       stop( error_message );
	       exit( -1 );

	   } else if (LocalSearchMethod==NULL) {

	       prStr(error_message, "ERROR:  Need to use \"set_sw1\" or \"set_psw1\" to allocate a Local Optimization object.\n");
	       stop(error_message);
	       exit(-1);
	   }

	   pr( logFile, "Number of Local Search (LS) only dockings = %d run%c\n", nruns, (nruns > 1)?'s':' ');


	   evaluate.setup(crd, charge, type, natom, map, inv_spacing, 
	      elec, emap,
	      xlo, xhi, ylo, yhi, zlo, zhi, nonbondlist,
	      e_internal, Nnb, B_calcIntElec, q1q2,B_isGaussTorCon,B_isTorConstrained,
	      B_ShowTorE, US_TorE, US_torProfile, vt, tlist, crdpdb, sInit, mol);

	   for (j=0; j<nruns; j++) {

	       pr( logFile, "\nDoing Local Search run: %d / %d.\n", j+1, nruns );
	       if (timeSeedIsSet[0] == 'T') {
	           seed[0] = (FourByteLong)time( &time_seed );
	       }
	       if (timeSeedIsSet[1] == 'T') {
	           seed[1] = (FourByteLong)time( &time_seed );
	       }
	       setall(seed[0], seed[1]);
	       initgn(-1); // Reinitializes the state of the current generator

	       pr(logFile, "Seeds: %ld %ld\n", seed[0], seed[1]);
	       (void) fflush( logFile );
	       gaStart = times( &tms_gaStart );

	       sHist[nconf] = call_ls(LocalSearchMethod, sInit, pop_size, &mol);

	       pr(logFile, "There were %lu Energy Evaluations.\n", evaluate.evals());

	       gaEnd = times( &tms_gaEnd );
	       pr( logFile, "Time taken for this Local Search (LS) run:\n");
	       timesyshms( gaEnd - gaStart, &tms_gaStart, &tms_gaEnd );
	       pr( logFile, "\n");
	       (void) fflush( logFile );
		
	       pr( logFile, UnderLine );
	       pr( logFile, "\n\n\tFINAL LOCAL SEARCH DOCKED STATE\n" );
	       pr( logFile,     "\t_______________________________\n\n\n" );
	       
	       writePDBQState( j, FN_ligand, dock_param_fn, lig_center, 
		    &(sHist[nconf]), ntor, &eintra, &einter, natom, atomstuff, 
		    crd, emap, elec, charge, torsFreeEnergy,
		    vt, tlist, crdpdb, nonbondlist, e_internal,
		    type, Nnb, B_calcIntElec, q1q2,
		    map, inv_spacing, xlo, ylo, zlo, xhi, yhi, zhi);

	       econf[nconf] = eintra + einter; // new2
	       
	       ++nconf;

	   } // Next run
	   (void) fflush(logFile);
       } else {
	    (void)fprintf(logFile, "NOTE: Command mode has been set, so local search cannot be performed.\n\n");
       }
       break;

/*____________________________________________________________________________*/

    case DPF_GS:
      (void) sscanf(line, "%*s %d", &nruns);

      if (!command_mode) {
	  if (nruns>MAX_RUNS) {

	      prStr(error_message, "ERROR:  %d runs requested, but only dimensioned for %d.\nChange \"MAX_RUNS\" in \"constants.h\".", nruns, MAX_RUNS);
	      stop(error_message);
	      exit(-1);

	  } else if (GlobalSearchMethod==NULL) {

	      prStr(error_message, "ERROR:  Need to use \"set_ga\" to allocate a Global Optimization object.\n");
	      stop(error_message);
	      exit(-1);
	  }

	  pr(logFile, "Number of Genetic Algorithm (GA) only dockings = %d run%c\n", nruns, (nruns>1)?'s':' ');


	  evaluate.setup(crd, charge, type, natom, map, inv_spacing, 
	     elec, emap,
	     xlo, xhi, ylo, yhi, zlo, zhi, nonbondlist,
	     e_internal, Nnb, B_calcIntElec, q1q2, B_isGaussTorCon,B_isTorConstrained,
	     B_ShowTorE, US_TorE, US_torProfile, vt, tlist, crdpdb, sInit, mol);

	  for (j=0; j<nruns; j++) {

	      fprintf( logFile, "\n\n\tBEGINNING GENETIC ALGORITHM DOCKING\n");
	      pr(logFile, "\nDoing Genetic Algorithm run:  %d/%d.\n", j+1, nruns);

	      if (timeSeedIsSet[0] == 'T') {
		  seed[0] = (FourByteLong)time( &time_seed );
	      }
	      if (timeSeedIsSet[1] == 'T') {
		  seed[1] = (FourByteLong)time( &time_seed );
	      }
	      setall(seed[0], seed[1]);
	      initgn(-1); // Reinitializes the state of the current generator

	      pr(logFile, "Seeds:  %ld %ld\n", seed[0], seed[1]);
	      (void) fflush(logFile);

	      gaStart = times(&tms_gaStart);

	      sHist[nconf] = call_gs(GlobalSearchMethod, sInit, num_evals, pop_s