void wham()

in Documents/NamdSample/namd210/lib/replica/FEP_wca/wham/wham-fep.cpp [167:224]

• 44 lines of code
• 11 McCabe index (conditional complexity)

void wham(void)
{
	int i, j, k, t, Converged=0, icycle;
	double arg1, arg2, argmax=0.0, Bot_i, Top_i, diff, *pDataTmp, kbt_Inv;
	
	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<NPoints[i]; t++)	{
					arg1   = (-lambda[k]*epprtd[i][t])*kbt_Inv;
					argmax = arg1;
					
					for(j=0; j<=nb_wind; j++)	{
						arg2  = (F2[j]-lambda[j]*epprtd[i][t])*kbt_Inv;
						
						if(arg2 > argmax)	argmax = arg2;
						pDataTmp[j] = arg2;
					}
					Top_i = exp(arg1-argmax);
					
					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;
			return;
		}
		
	}
	
	printf("\nFail to converge in %d steps.\n", NCYCLE);
	delete []pDataTmp;
}