void wham_E_Histogram()

in Documents/NamdSample/namd210b/lib/replica/FEP_wca/wham/wham-fep.cpp [229:359]
106 lines of code
20 McCabe index (conditional complexity)

void wham_E_Histogram(void)
{
	int i, j, k, t, Converged=0, icycle, Idx;
	double arg1, arg2, argmax=0.0, Bot_i, Top_i, diff, *pDataTmp, kbt_Inv, E_Local;
	double *pEmin=NULL, *pEmax=NULL, **pE_Hist=NULL;
	int *pN_Bin=NULL, **nHist_E=NULL;

	//start	to build the histogram of E for each window and free the memory for all data points, epprtd[][]
	pEmin = new double[nb_wind+1];
	CHECK_POINTER(pEmin)
	pEmax = new double[nb_wind+1];
	CHECK_POINTER(pEmax)
	pN_Bin = new int[nb_wind+1];
	CHECK_POINTER(pN_Bin)
	nHist_E = new int*[nb_wind+1];
	CHECK_POINTER(nHist_E)
	pE_Hist = new double*[nb_wind+1];
	CHECK_POINTER(pE_Hist)
	for(i=0; i<=nb_wind; i++)	{
		pEmin[i] = 1.0E100;
		pEmax[i] = -1.0E100;
		for(j=0; j<NPoints[i]; j++)	{
			if(epprtd[i][j] > pEmax[i])	{
				pEmax[i] = epprtd[i][j];
			}
			if(epprtd[i][j] < pEmin[i])	{
				pEmin[i] = epprtd[i][j];
			}
		}
		pN_Bin[i] = (int)((pEmax[i]-pEmin[i])/E_BIN_SIZE)+1;
		nHist_E[i] = new int[pN_Bin[i]];
		CHECK_POINTER(nHist_E[i])
		pE_Hist[i] = new double[pN_Bin[i]];
		CHECK_POINTER(pE_Hist[i])
		memset(nHist_E[i], 0, sizeof(int)*pN_Bin[i]);

		for(j=0; j<pN_Bin[i]; j++)	{	// assign the list of energies histogram
			pE_Hist[i][j] = pEmin[i] + j*E_BIN_SIZE;
		}

		for(j=0; j<NPoints[i]; j++)	{
			Idx = (int)((epprtd[i][j] - pEmin[i])/E_BIN_SIZE);
			nHist_E[i][Idx]++;
		}
	}
	Deallocate_Mem_2D_Array(epprtd, max_wind, max_data);
	epprtd = NULL;
	//end	to build the histogram of E for each window and free the memory for all data points, epprtd[][]


	kbt_Inv = 1.0/kbt;
	pDataTmp = new double[nb_wind+1];

	for(icycle=1; icycle<=NCYCLE; icycle++)	{
		for(k=0; k<=nb_wind; k++)	{
			F[k]=0.0;
			
			for(i=0; i<=nb_wind; i++)	{
				for(t=0; t<pN_Bin[i]; t++)	{
					if(nHist_E[i][t] == 0)	{
						continue;
					}
					E_Local = pE_Hist[i][t];
					arg1   = (-lambda[k]*E_Local)*kbt_Inv;
					argmax = arg1;
					
					for(j=0; j<=nb_wind; j++)	{
						arg2  = (F2[j]-lambda[j]*E_Local)*kbt_Inv;
						
						if(arg2 > argmax)	argmax = arg2;
						pDataTmp[j] = arg2;
					}
					Top_i = exp(arg1-argmax) * nHist_E[i][t];
					
					Bot_i  = 0.0;
					for(j=0; j<=nb_wind; j++)	{
						arg2  = pDataTmp[j]-argmax;
						Bot_i  = Bot_i  + NPoints[j] * exp(arg2);
					}
					F[k] = F[k] + Top_i / Bot_i;
				}
			}
		}
		
		Converged = 1;
		for(k=0; k<=nb_wind; k++)	{
			F[k] = -kbt*log( F[k] );
//			printf("window %d: F(k) = %.5lf\n", k+1, F[k]);
			diff = fabs( F2[k] - F[k] );
			
			F2[k] = F[k] - F[0];
			
			if(diff > TOL) Converged = 0;
		}
		
		if(Converged) {
			printf("\nConverged at step: %d\n", icycle);
			delete []pDataTmp;

			for(i=0; i<=nb_wind; i++)	{
				delete [](nHist_E[i]);
				delete [](pE_Hist[i]);
			}
			delete []nHist_E;
			delete []pE_Hist;

			delete []pEmin;
			delete []pEmax;
			delete []pN_Bin;

			return;
		}
		
	}
	
	printf("\nFail to converge in %d steps.\n", NCYCLE);
	delete []pDataTmp;

	for(i=0; i<=nb_wind; i++)	{
		delete [](nHist_E[i]);
		delete [](pE_Hist[i]);
	}
	delete []nHist_E;
	delete []pE_Hist;
	
	delete []pEmin;
	delete []pEmax;
	delete []pN_Bin;

	return;
}