simanneal.cc

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

		if (B_within_constraint) {
		    /* 
		    ** Normally, this is true.
		    ** If the distance-constraint was set,
		    ** and was just violated, this is false.
		    */

		    for (i = 0;  i < natom;  i++) {
			B_outside= is_out_grid(crd[i][X], crd[i][Y], crd[i][Z]);
			if ( B_outside ) {
			    /*
			    ** Outside grid!
			    */
			    ++nedge;
			    ++nrej;
			    /*pr(logFile, "e"); / *###*/
			    /* etot += WallEnergy; ++ntot;*/
			    /*
			    ** Undo this move,
			    */
			    copyState( &sNow, sLast );
			    /*
			    ** Subtract just the translation;
			    ** hence "doubling back" from the edge.
			    */
			    sNow.T.x -= sChange.T.x;
			    sNow.T.y -= sChange.T.y;
			    sNow.T.z -= sChange.T.z;

			    if ( B_writeTrj && B_inRange && B_either && (++count == trj_freq)) {
				count = 0;
				e = eintra = WallEnergy;
				output_state( FP_trj, sNow, ntor, nacc+nrej, e,
				    eintra, (char)'e', B_watch, FN_watch, 
				    atomstuff, natom, crd);
			    }/*writeTrj*/

			    break;/*...out of i*/

			}/*outside*/
		    }/*for atoms i*/

		    if ( !B_outside ){ /*inside grid maps*/

 			/* Calculate Energy of System, =======================*/

			/*
			** MORE ACCURATE METHOD, (SLOWER):
			*/
			e = quicktrilinterp( crd, charge, type, natom, map, inv_spacing, xlo, ylo, zlo) + (eintra = eintcal( nonbondlist, e_internal, crd, type, Nnb, B_calcIntElec, q1q2));

			/*
			** LESS ACCURATE  METHOD (FASTER):
			** (alternative to above line).
			** 
			** Only calculate internal energy if 
			** trilinterp energy is low enough.
			** 
			** if ( (e = quicktrilinterp( crd, charge, type, natom, 
			** map, inv_spacing, xlo, ylo, zlo)) < 
			** ENERGY_CUTOFF) { 
			**   e += (eintra = eintcal( nonbondlist, e_internal,
			**   crd, type, Nnb, B_calcIntElec, q1q2 ));
			** }
			*/

			if (B_isGaussTorCon) {
			    /*** This looks wrong... for (Itor = 0; Itor <= ntor; Itor++) { ***/
			    for (Itor = 0; Itor < ntor; Itor++) {
				if (B_isTorConstrained[Itor] == 1) {
				    indx = Rad2Div( sNow.tor[Itor] );
				    if (B_ShowTorE) {
					e += (float)( US_TorE[Itor] 
						  = US_torProfile[Itor][indx] );
				    } else {
					e += (float)US_torProfile[Itor][indx];
				    }
				}
			    }
			}

			/* Apply the Metropolis energy test... ===============*/

			if (e <= eLast) {
			    /*
			    **  Accept this move immediately.
			    */
			    ++nacc;
			    ++nAcc;
			    etot += (eLast = e);
			    ++ntot;
			    copyState( &sLast, sNow );

			    /* pr(logFile, "A"); / *###*/
			    if (e < eMin) {
				/*
				** Update minimum-energy state variables,
				*/
				eMin = e;
				copyState( &sMin,  sNow );
			    }
			    if ( B_writeTrj && B_inRange && (++count == trj_freq)){
				count = 0;
				output_state( FP_trj, sNow, ntor, nacc+nrej, e,
				    eintra, (char)'A', B_watch, FN_watch, 
				    atomstuff, natom, crd);
			    }/*write trajectory*/

			} else {

                            /* Probabilistic move. ===========================*/

			    if (exp((double)((eLast-e)*inv_RT))<local_random()){
				/*
				** Failed the probability test. 
				** Reject this move.
				*/
				++nrej;
				/* pr(logFile, "R"); / *###*/
				copyState( &sNow, sLast );

				if ( B_writeTrj && B_inRange && B_either && 
				    (++count == trj_freq)) {
				    count = 0;
				    output_state( FP_trj, sNow, ntor, nacc+nrej,
				        e, eintra, (char)'R', B_watch, FN_watch,
					atomstuff, natom, crd);
				}/*write trajectory*/

			    } else {
				/*
				** Passed the probability test.  
				** Accept this move.
				** A chance to escape a local minimum...
				*/
				++nacc;
				++nAccProb;
				etot += e;
				++ntot;
				/*pr(logFile, "a"); / *###*/

				eLast = e;
				copyState( &sLast, sNow );
				if ( B_writeTrj && B_inRange && (++count == trj_freq))  {
				    count = 0;
				    output_state(FP_trj, sNow, ntor, nacc+nrej,
					e, eintra, (char)'a', B_watch, FN_watch,
					atomstuff, natom, crd);
				}/*write trajectory*/
			    }/*passed Monte Carlo probablility test*/
			}/*e > eLast, Probabilistic move...*/
		    }/*inside grid maps*/
		}/*within_constraint*/

            } while ( nacc < naccmax  &&  nrej < nrejmax );
/*____________________________________________________________________________*/

	    if ((nacc == 0) && (nedge == nrejmax)) {
		/*
		**  Clear indication that ligand got stuck on an edge...
		*/
		pr(logFile, "\n\n>>> Ligand appears to be stuck on an edge: forced to re-initialize. <<<\n\n");
		--icycle;
		getInitialState( &e0, e0max,
			 &sInit, &sMin, &sLast, 
			 TRUE, TRUE, TRUE, 
			 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);

	    } else {

		if ( B_trnReduc )  trnStep *= trnFac;       
		if ( B_qtwReduc )  qtwStep *= qtwFac;
		if ( B_torReduc )  torStep *= torFac;
		/*
		**  Output-level dependent diagnostics...
		*/
		if (outlev > 0) {
		    /*pr(logFile, "\n"); / *###*/
		    pr( logFile, "%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, NcycMax, eMin, etot/ntot, (nrej!=0) ? (float)nacc/nrej : 999.99, nAcc, nAccProb, nrej, nedge, RT, sMin.T.x, sMin.T.y, sMin.T.z );
		    cycEnd = times( &tms_cycEnd );
		    timesys( cycEnd - cycStart, &tms_cycStart, &tms_cycEnd );
		    if (outlev > 1) {
			pr( logFile, "\tEnergy:   \tState:\n\t__________\t____________________________________________________________\nMinimum\t%+6.2f\t(%+.2f,%+.2f,%+.2f), q = [w,(x,y,z)] = [%5.1f deg, (%+.2f,%+.2f,%+.2f)],\n", eMin, sMin.T.x, sMin.T.y, sMin.T.z, Deg(sMin.Q.ang) , sMin.Q.nx, sMin.Q.ny, sMin.Q.nz );
			pr( logFile, "\nLast\t%+6.2f\t(%+.2f,%+.2f,%+.2f), q = [w,(x,y,z)] = [%5.1f deg, (%+.2f,%+.2f,%+.2f)],\n", eLast, sLast.T.x, sLast.T.y, sLast.T.z, Deg(sLast.Q.ang) , sLast.Q.nx, sLast.Q.ny, sLast.Q.nz );
			if (ntor > 0) {
			    pr( logFile, "Minimum:\t(" );
			    for (i=0; i<ntor; i++) {
				pr( logFile, "%.1f%s ", Deg(sMin.tor[i]), (i < ntor1)?",":" deg)" );
			    }
			    pr( logFile, "\nLast:\t(" );
			    for (i=0; i<ntor; i++) {
				pr( logFile, "%.1f%s ", Deg(sLast.tor[i]), (i < ntor1)?",":" deg)" );
			    }
			    pr( logFile, "\n" );
			}
			if ( B_trnReduc )
			    pr( logFile, "\nTranslation step size reduced; now =\t\t +/- %.2f A\n", trnFac);
			if ( B_qtwReduc )
			    pr( logFile, "\nQuaternion Rotation step size reduced; now =\t +/- %.2f deg\n", qtwFac);
			if ( B_torReduc )
			    pr( logFile, "\nTorsion step size reduced; now =\t\t +/- %.2f deg\n", torFac);
		    }/*outlev > 1*/
		    flushLog;
		}/*outlev > 0*/
		/*
		** Reduce temperature,
		*/
		if ( B_linear_schedule ) {
		    RT -= RTreduc;

		    if (RT <= APPROX_ZERO) {
			inv_RT = 1. / APPROX_ZERO;
		    } else {
			inv_RT = 1. / RT;
		    }
		} else if ( B_RTChange ) {
		    inv_RT = 1./( RT *= RTFac );
		}
		/*
		** Start next cycle at minimum state?
		*/
		if ( B_selectmin ) {
		    eLast = eMin;
		    copyState( &sLast, sMin );
		}
	    } /* Prepares for next cycle */

        } /* icycle */
/*____________________________________________________________________________*/

        if ( B_selectmin ) {
	    copyState( &sSave, sMin );
        } else {
	    copyState( &sSave, sLast );
        }

	sSave.Q.ang = WrpRad( ModRad( sSave.Q.ang ) );

        for (i=0; i<ntor; i++) {
            sSave.tor[i] = WrpRad( ModRad( sSave.tor[i] ) );
        }
   
        pr( logFile, UnderLine );
        pr( logFile, "\n\n\t\tFINAL DOCKED STATE\n" );
        pr( logFile,     "\t\t__________________\n\n\n" );
        pr( logFile, "Run Number %d, \n\nFinal Energy = %+.2f\n", irun1, eLast);

        pr( logFile, "Final Translation = %.2f, %.2f, %.2f\n", sSave.T.x, sSave.T.y, sSave.T.z );
        pr( logFile, "Final Quaternion Rotation Angle = %5.1f deg\n", Deg(sSave.Q.ang) );
        pr( logFile, "Final Quaternion Unit Vector = ( %+.2f, %+.2f, %+.2f )\n", sSave.Q.nx, sSave.Q.ny, sSave.Q.nz );

	copyState( &sHist[ *Addr_nconf ], sSave );

        if (ntor > 0) {
            pr( logFile, "Final Torsions:\n" );
            for (i=0; i<ntor; i++) {
		torTmp = Deg( sSave.tor[i] );
		torTmp = ModDeg( torTmp );
		torTmp = WrpDeg( torTmp );
                sHist[ *Addr_nconf ].tor[i] = sSave.tor[i];
                pr( logFile, "          %2d = %7.2f deg", i+1, torTmp);
		if ((B_isTorConstrained[i] == 1) && B_ShowTorE) {
		    pr(logFile, ", Energetic penalty = %uhd\n", US_TorE[i]);
		} else {
		    pr(logFile, "\n");
		}
            }
        }

	cnv_state_to_coords( sSave, vt, tlist, ntor, crdpdb, crd, natom );

	if (ntor > 0) {
	    eintra = eintcal( nonbondlist, e_internal, crd, type, Nnb, 
		B_calcIntElec, q1q2);
	} else {
	    eintra = 0.0;
	}
        einter = trilinterp( crd, charge, type, natom, map, inv_spacing, 
	    elec, emap, xlo, ylo, zlo);

	writePDBQ( irun, FN_ligand, FN_dpf, sml_center, sSave, ntor,
	  eintra, einter, natom, atomstuff, crd, emap, elec, charge,
	  torsFreeEnergy);

        econf[(*Addr_nconf)] = eLast;
        ++(*Addr_nconf);

    } /* Loop over runs ======================================================*/

    if ( B_writeTrj ) {
    	fclose( FP_trj );
    }
}
/* EOF */