#include "Meta1D.cuh"
float BSpline_Interpolation_4(int i, float x)// in fact need 1.-x
{
	if (i == -1)
	{
		return 1. / 6.*x*x*x;
	}
	else if (i == 0)
	{
		return -0.5*x*x*x + 0.5*x*x + 0.5*x + 1. / 6.;
	}
	else if (i == 1)
	{
		return 0.5*x*x*x - x*x + 2. / 3.;
	}
	else if (i == 2)
	{
		return -1. / 6.*x*x*x + 0.5*x*x - 0.5*x + 1. / 6.;
	}
	else
	{
		printf("fuck\n");
		getchar();
		return 0.;
	}
}
float DBSpline_Interpolation_4(int i, float x)// in fact need 1.-x, and is -D
{
	if (i == -1)
	{
		return 3. / 6.*x*x;
	}
	else if (i == 0)
	{
		return -0.5*3.*x*x + 0.5*2.*x + 0.5;
	}
	else if (i == 1)
	{
		return 0.5*3.*x*x - 2.*x;
	}
	else if (i == 2)
	{
		return -3. / 6.*x*x + 0.5*2.*x - 0.5;
	}
	else
	{
		printf("fuck\n");
		getchar();
		return 0.;
	}
}
void META::Initial(CONTROLLER *controller, char * module_name)
{
	if (module_name == NULL)
	{
		strcpy(this->module_name, "meta");
	}
	else
	{
		strcpy(this->module_name, module_name);
	}
	sprintf(read_potential_file_name, "Meta_Potential.txt");
	sprintf(write_potential_file_name, "Meta_Potential.txt");

	printf("START INITIALIZING META:\n");

	if (controller[0].Command_Exist(this->module_name, "Read_Potential_File"))
	{
		sscanf(controller[0].Command(this->module_name, "Read_Potential_File"), "%s", read_potential_file_name);
		printf("\tSTART READING Potential_File: %s\n", read_potential_file_name);
		Read_Potential();
	}
	else
	{
		int check_necessary_inpu_exist = 0;
		if (controller[0].Command_Exist(this->module_name, "Cv_Minimal"))
		{
			cv_min = atof(controller[0].Command(this->module_name, "Cv_Minimal"));
			check_necessary_inpu_exist = check_necessary_inpu_exist + 1;
		}
		if (controller[0].Command_Exist(this->module_name, "Cv_Maximum"))
		{
			cv_max = atof(controller[0].Command(this->module_name, "Cv_Maximum"));
			check_necessary_inpu_exist = check_necessary_inpu_exist + 1;
		}
		if (controller[0].Command_Exist(this->module_name, "dCv"))
		{
			dcv = atof(controller[0].Command(this->module_name, "dCv"));
		}
		
		if ((float)(cv_max - cv_min) / dcv < 1.0)
		{
			check_necessary_inpu_exist = -1;
		}
		if (check_necessary_inpu_exist != 2)
		{
			printf("	Please check if Meta necessary input is complete and correct!\n");
			getchar();
		}

		grid_numbers = (float)(cv_max - cv_min) / dcv + 1.;
		dcv = (float)(cv_max - cv_min) / (grid_numbers - 1);

		Malloc_Safely((void**)&potential_list, sizeof(float)*grid_numbers);
		for (int i = 0; i < grid_numbers; i = i + 1)
		{
			//potential_list[i] = 10.*dcv*(grid_numbers-i);
			potential_list[i] = 0;
		}

	}
	if (controller[0].Command_Exist(this->module_name, "Height"))
	{
		height = atof(controller[0].Command(this->module_name, "Height"));
	}
	if (controller[0].Command_Exist(this->module_name, "Sigma"))
	{
		sigma = atof(controller[0].Command(this->module_name, "Sigma"));
	}
	if (controller[0].Command_Exist(this->module_name, "Wall_height"))
	{
		border_potential_height = atof(controller[0].Command(this->module_name, "Wall_height"));
	}
	if (controller[0].Command_Exist(this->module_name, "WellTemp_Factor"))
	{
		welltemp_factor = atof(controller[0].Command(this->module_name, "WellTemp_Factor"));
	}
	if (controller[0].Command_Exist(this->module_name, "Write_Potential_File"))
	{
		sscanf(controller[0].Command(this->module_name, "Write_Potential_File"), "%s", write_potential_file_name);
	}
	if (controller[0].Command_Exist(this->module_name, "Potential_Update_Interval"))
	{
		potential_update_interval = atof(controller[0].Command(this->module_name, "Potential_Update_Interval"));
	}
	else
	{
		printf("\tNo Potential_Update_Intervel input, use write_information_interval instead!\n");
		potential_update_interval = atof(controller[0].Command("write_information_interval")); //Potential更新频率，未定义时默认设置为体系轨迹更新频率
	}
	controller[0].Step_Print_Initial("metad", "%.2f");

	printf("\tPotential_Update_Interval: %d\n", potential_update_interval);
	printf("\tCv_Minimal: %f\n", cv_min);
	printf("\tCv_Maximum: %f\n", cv_max);
	printf("\tdCv: %f\n", dcv);
	printf("\tBias_Height: %f\n", height);
	printf("\tBias_Sigma: %f\n", sigma);
	printf("\tBondary_Height: %f\n", border_potential_height);
	printf("\tWellTemp_Factor: %f\n", welltemp_factor);
	printf("\tWrite_Potential_File: %s\n", write_potential_file_name);
	is_initialized = 1;
	printf("END (INITIALIZING META)\n");
}
void META::Write_Potential()
{
	FILE *temp_file = NULL;
	Open_File_Safely(&temp_file, write_potential_file_name, "w");
	fprintf(temp_file, "system name and description\n");
	fprintf(temp_file, "%f %f %f\n",cv_min,cv_max,dcv);
	fprintf(temp_file, "%d\n", grid_numbers);
	for (int i = 0; i < grid_numbers; i = i + 1)
	{
		fprintf(temp_file, "%f %f\n", (float)dcv*i + cv_min, potential_list[i]);
	}
	fclose(temp_file);
}
void META::Read_Potential()
{
	FILE *temp_file = NULL;
	Open_File_Safely(&temp_file, read_potential_file_name, "r");
	char temp_char[256];
	fgets(temp_char, 256, temp_file);
	fscanf(temp_file, "%f %f %f\n", &cv_min, &cv_max, &dcv);
	fscanf(temp_file, "%d\n", &grid_numbers);
	Malloc_Safely((void**)&potential_list, sizeof(float)*grid_numbers);
	float temp_float;
	for (int i = 0; i < grid_numbers; i = i + 1)
	{
		fscanf(temp_file, "%f %f\n", &temp_float, &potential_list[i]);
	}
	fclose(temp_file);
}
//这一部分修改过了，20210524
__global__ void Add_Frc(const int atom_numbers, VECTOR *frc, VECTOR *cv_grad, float dheight_dcv)
{
	for (int i = threadIdx.x; i < atom_numbers; i = i + blockDim.x)
	{
		frc[i] = frc[i] + dheight_dcv*cv_grad[i];
	}
}
//__global__ void Meta_Force_CUDA
//(int atom_numbers, VECTOR *frc, float *cv, VECTOR *cv_grad,
//int grid_numbers,float cv_min,float cv_max,float dcv_inverse,float *potential_list
//)
//{
//	float cv_temp = cv[0];
//	cv_temp = fminf(fmaxf(cv_temp, cv_min),cv_max);
//	int grid_i=(cv[0] - cv_min)*dcv_inverse;
//	float dheight_dcv = (potential_list[grid_i + 1] - potential_list[grid_i])*dcv_inverse;
//	Add_Frc << <1, 256 >> > (atom_numbers, frc, cv_grad, fb, dheight_dcv);
//}
void META::Meta_Force(int atom_numbers, VECTOR *frc, VECTOR *cv_grad, float dheight_dcv)
{
	Add_Frc << <1, 256 >> >	(atom_numbers, frc, cv_grad, dheight_dcv);
}
float META::Potential(float x)
{
	float temp_pos = (float)(x - cv_min) / dcv;
	int pos = (float)temp_pos;
		if (pos >= 1 && pos <= grid_numbers - 3)
	{
		float scale = 1. - temp_pos + (float)pos;
		return BSpline_Interpolation_4(-1, scale)*potential_list[pos - 1]
			+ BSpline_Interpolation_4(0, scale)*potential_list[pos]
			+ BSpline_Interpolation_4(1, scale)*potential_list[pos + 1]
			+ BSpline_Interpolation_4(2, scale)*potential_list[pos + 2];
	}
	else if (pos == 0)
	{
		float scale = 1. - temp_pos + (float)pos;
		return BSpline_Interpolation_4(-1, scale)*border_potential_height
			+ BSpline_Interpolation_4(0, scale)*potential_list[0]
			+ BSpline_Interpolation_4(1, scale)*potential_list[1]
			+ BSpline_Interpolation_4(2, scale)*potential_list[2];
	}
	else if (pos == grid_numbers - 2)
	{
		float scale = 1. - temp_pos + (float)pos;
		return BSpline_Interpolation_4(-1, scale)*potential_list[grid_numbers - 3]
			+ BSpline_Interpolation_4(0, scale)*potential_list[grid_numbers - 2]
			+ BSpline_Interpolation_4(1, scale)*potential_list[grid_numbers - 1]
			+ BSpline_Interpolation_4(2, scale)*border_potential_height;
	}
	else
	{
		printf("POTENTIAL::Potential Warning: Input of x=%f out of range!\n", x);
		return 0.;
	}
}
float META::DPotential(float x)
{
	float temp_pos = (float)(x - cv_min) / dcv;
	int pos = (float)temp_pos;
	if (pos >= 1 && pos <= grid_numbers - 3)
	{
		float scale = 1. - temp_pos + (float)pos;
		return DBSpline_Interpolation_4(-1, scale)*potential_list[pos - 1]
			+ DBSpline_Interpolation_4(0, scale)*potential_list[pos]
			+ DBSpline_Interpolation_4(1, scale)*potential_list[pos + 1]
			+ DBSpline_Interpolation_4(2, scale)*potential_list[pos + 2];
	}
	else if (pos == 0)
	{
		float scale = 1. - temp_pos + (float)pos;
		return DBSpline_Interpolation_4(-1, scale)*border_potential_height
			+ DBSpline_Interpolation_4(0, scale)*potential_list[0]
			+ DBSpline_Interpolation_4(1, scale)*potential_list[1]
			+ DBSpline_Interpolation_4(2, scale)*potential_list[2];
	}
	else if (pos == grid_numbers - 2)
	{
		float scale = 1. - temp_pos + (float)pos;
		return DBSpline_Interpolation_4(-1, scale)*potential_list[grid_numbers - 3]
			+ DBSpline_Interpolation_4(0, scale)*potential_list[grid_numbers - 2]
			+ DBSpline_Interpolation_4(1, scale)*potential_list[grid_numbers - 1]
			+ DBSpline_Interpolation_4(2, scale)*border_potential_height;
	}
	else
	{
		printf("POTENTIAL::Potential Warning: Input of x=%f out of range!\n", x);
		return 0.;
	}
}
void META::AddPotential(float h_cv, int steps)
{
	if (steps % potential_update_interval == 0)
	{
		for (int i = 0; i < grid_numbers; i = i + 1)
		{
			float pos = (float)i*dcv + cv_min;
			float add = height *1. / sqrtf(2.*3.141592654) / sigma * expf(-(pos - h_cv)*(pos - h_cv) / 2. / sigma / sigma);
			potential_list[i] = potential_list[i] + add*expf(-potential_list[i] / welltemp_factor);
		}
		float baseline = potential_list[grid_numbers - 1];
		for (int i = 0; i < grid_numbers; i = i + 1)
		{
			potential_list[i] = potential_list[i] - baseline;
		}
	}
}