#pragma warning(disable:4996)£»
#include<vector>
#include <algorithm>
#include"read_write.h"
#include"variant_command.h"
using namespace std;
namespace EVO {
	void Exchange_command(
		int    nant_o[],                  // total number of atoms and types
		int    typenum_o[],               // number of each type 
		double xyz_o[][3],                // atomic coordinates
		double xyz_e[][3],
		int NumExg,
		int Num1[],
		int Num2[]
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		//Exchange operator which exchanges the coordinates of two atoms of different types
		fprintf(fp_olog, ">>>Exchange operator \n");
		double temp[3];
		fprintf(fp_olog, ">>>the num of Exchange is\n ");
		fprintf(fp_olog, "%d\n", NumExg);
		fprintf(fp_olog, ">>>The atoms selected for exchange are\n");
		for (int i = 0; i < NumExg; i++)
		{
			for (int j = 0; j < 3; j++) {
				temp[j] = xyz_o[Num1[i]][j];
				xyz_o[Num1[i]][j] = xyz_o[Num2[i]][j];
				xyz_o[Num2[i]][j] = temp[j];
			}
			fprintf(fp_olog, "%d  %d\n", Num1[i], Num2[i]);
		}
		for (int i = 0; i < nant_o[0]; i++)
		{
			for (int j = 0; j < 3; j++)
			{
				xyz_e[i][j] = xyz_o[i][j];
			}
		}
		fclose(fp_olog);
	}

	void Strain_command(
		char   mode[],
		double vec_o[3][3],               // lattice vectors
		double vec_e[3][3],
		double strain0, double strain1, double strain2,
		double strain3, double strain4, double strain5
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		fprintf(fp_olog, ">>>Strain operator \n");
		int i, j;
		double defvect[3][3] = { {1.0000, 0.0000, 0.0000},{0.0000, 1.0000, 0.0000},{0.0000, 0.0000, 1.0000} };
		//tensile strain along x axis
		if (!strcmp(mode, "xx"))
		{
			fprintf(fp_olog, ">>>the values of tensile strain along xx axis is \n%1.4f \n", strain0);
			defvect[0][0] = 1.0000 + strain0;
		}
		//tensile strain along y axis
		if (!strcmp(mode, "yy"))
		{
			fprintf(fp_olog, ">>> the values of tensile strain along yy axis is \n%1.4f \n", strain1);
			defvect[1][1] = 1.0000 + strain1;
		}
		//tensile strain along z axis
		if (!strcmp(mode, "zz"))
		{
			fprintf(fp_olog, ">>>the values of tensile strain along zz axis is \n%1.4f \n", strain2);
			defvect[2][2] = 1.0000 + strain2;
		}
		// shear strain in x-y plane
		if (!strcmp(mode, "xy"))
		{
			fprintf(fp_olog, ">>>the values of the shear strain in xy plane is \n%1.4f \n", strain3);
			defvect[0][1] = strain3 / 2.0000;
			defvect[1][0] = strain3 / 2.0000;
		}
		//shear strain in x-z plane
		if (!strcmp(mode, "xz"))
		{
			fprintf(fp_olog, ">>> the values of the shear strain in xz plane is \n%1.4f \n", strain4);
			defvect[0][2] = strain4 / 2.0000;
			defvect[2][0] = strain4 / 2.0000;
		}
		//shear strain in y-z plane
		if (!strcmp(mode, "yz"))
		{
			fprintf(fp_olog, ">>>the values of the shear strain in yz plane is \n%1.4f \n", strain5);
			defvect[1][2] = strain5 / 2.0000;
			defvect[2][1] = strain5 / 2.0000;
		}
		//strain
		if (!strcmp(mode, "xyz"))
		{
			fprintf(fp_olog, ">>>the values of the pure strain in xyz is that: \n");
			fprintf(fp_olog, "%1.4f  %1.4f  %1.4f  %1.4f  %1.4f  %1.4f \n",
				strain0, strain1, strain2, strain3, strain4, strain5);
			defvect[0][0] = 1.0000 + strain0; defvect[0][1] = strain3 / 2.0000; defvect[0][2] = strain4 / 2.0000;
			defvect[1][0] = strain3 / 2.0000; defvect[1][1] = 1.0000 + strain1; defvect[1][2] = strain5 / 2.0000;
			defvect[2][0] = strain4 / 2.0000; defvect[2][1] = strain5 / 2.0000; defvect[2][2] = 1.0000 + strain2;
		}

		for (i = 0; i < 3; i++) {
			for (j = 0; j < 3; j++) {
				vec_e[i][j] = vec_o[i][0] * defvect[0][j] + vec_o[i][1] * defvect[1][j] + vec_o[i][2] * defvect[2][j];
			}
		}
		fclose(fp_olog);
	}

	void Ripple_command(
		int    nant_o[],                  // total number of atoms and types
		double xyz_o[][3],                // atomic coordinates
		double xyz_e[][3],
		int temp,
		double angle_1,
		double angle_2                   // randomly chosen from the interval 0~2pi to vary the regions of the structure that are strongly displaced
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		//Ripple opetator which shifts the coordinates of each atom long a random axis
		fprintf(fp_olog, ">>>Ripple operator \n");
		int i;
		char axis;
		int x, y, z;
		double MaxDisp = 0.5000;                     // the maximum possible displacement of any atom in the displaced(z) direction
		//temp = floor(6 * RandDigit());
		switch (temp)
		{
		case 0: x = 0, y = 1, z = 2; axis = 'z'; break;
		case 1: x = 0, y = 2, z = 1; axis = 'y'; break;
		case 2: x = 1, y = 0, z = 2; axis = 'z'; break;
		case 3: x = 1, y = 2, z = 0; axis = 'x'; break;
		case 4: x = 2, y = 0, z = 1; axis = 'y'; break;
		case 5: x = 2, y = 1, z = 0; axis = 'x'; break;
		}                                          // choosing a random axis
		fprintf(fp_olog, ">>>The displaced(z) direction is \n%c \n", axis);
		fprintf(fp_olog, ">>>angle_1, angle_2=\n%1.4f  %1.4f \n", angle_1, angle_2);
		double theta_1 = angle_1 * Pi / 180;
		double theta_2 = angle_2 * Pi / 180;
		for (i = 0; i < nant_o[0]; i++) {
			xyz_e[i][x] = xyz_o[i][x];
			xyz_e[i][y] = xyz_o[i][y];
			xyz_e[i][z] = xyz_o[i][z] + MaxDisp * cos(2 * Pi * xyz_o[i][x] + theta_1) * cos(2 * Pi * xyz_o[i][y] + theta_2);
			if (xyz_e[i][z] >= 1) { --xyz_e[i][z]; }
			if (xyz_e[i][z] < 0) { ++xyz_e[i][z]; }
		}
		fclose(fp_olog);
	}

	//The random move operator moves all atoms¡¯ positions randomly, keeping the lattice vectors unchanged
	void Randmove_command(
		int    nant_o[],
		double xyz_o[][3],
		double xyz_e[][3],
		int num,
		double dxyz
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		fprintf(fp_olog, ">>>Random move perator \n");
		fprintf(fp_olog, ">>>The distance of atom then moves is \n%lf\n", dxyz);
		if (num > nant_o[0])
		{
			fprintf(fp_olog, ">>>wrong,the randmove atom is higher than total atom!");
			return;
		}
		int count = 0;
		vector<int> dis(num);
		fprintf(fp_olog, ">>>The serial number of the moved atom is\n");
		for (int k = 0; k < num; k++)
		{
			dis[k] = rand() % nant_o[0];
			for (int l = 0; l < k; l++)
			{
				while (dis[k] == dis[l])
				{
					dis[k] = rand() % num;
				}
			}
			fprintf(fp_olog, "%d ", dis[k]);
			count++;
			if (count == num)
				break;
		}
		fprintf(fp_olog, "\n");
		int i, j;
		for (i = 0; i < nant_o[0]; i++)
		{
			for (j = 0; j < 3; j++)
			{
				if (find(dis.begin(), dis.end(), i) != dis.end())
				{
					xyz_e[i][j] = xyz_o[i][j] + dxyz;
				}
				else
					xyz_e[i][j] = xyz_o[i][j];
				if (xyz_e[i][j] >= 1) { --xyz_e[i][j]; }
				if (xyz_e[i][j] < 0) { ++xyz_e[i][j]; }
			}
		}

		fclose(fp_olog);
	}

	//This operator slips a group of atoms by a random distance along a random direction that is parallel to a specific plane
	void Slip_command(
		int    nant_o[],                  // total number of atoms and types
		double xyz_o[][3],                // atomic coordinates
		double xyz_e[][3],
		double zlower, double ydisp, double xlower,
		double zdisp, double ylower, double xdisp
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		fprintf(fp_olog, ">>>Slip operator \n");
		int i;
		fprintf(fp_olog, ">>>slip on y-z plane: zlower,slip on z-x plane: xlower,slip on x-y plane: ylower=\n%1.2f %1.2f %1.2f\n", zlower, xlower, ylower);
		fprintf(fp_olog, ">>>displace along y axis,displace along z axis,displace along x axis=\n%1.4f %1.4f %1.4f\n", ydisp, zdisp, xdisp);
		fclose(fp_olog);
		for (i = 0; i < nant_o[0]; i++) {
			//slip on y-z plane
			if (xyz_o[i][2] > zlower)
			{
				xyz_e[i][1] = xyz_o[i][1] + ydisp;
				if (xyz_e[i][1] > 1) { --xyz_e[i][1]; }//
			}
			else
			{
				xyz_e[i][1] = xyz_o[i][1];
			}

			//slip on z-x plane
			if (xyz_o[i][0] > xlower)
			{
				xyz_e[i][2] = xyz_o[i][2] + zdisp;
				if (xyz_e[i][2] > 1) { --xyz_e[i][2]; }
			}
			else
			{
				xyz_e[i][2] = xyz_o[i][2];
			}

			//slip on x-y plane
			if (xyz_o[i][1] > ylower)
			{
				xyz_e[i][0] = xyz_o[i][0] + xdisp;
				if (xyz_e[i][0] > 1) { --xyz_e[i][0]; }
			}
			else
			{
				xyz_e[i][0] = xyz_o[i][0];
			}
		}
	}

	void Twist_command(
		int    nant_o[],                  // total number of atoms and types
		double xyz_o[][3],                // atomic coordinates
		double xyz_e[][3],
		char axis,
		double angletemp
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		fprintf(fp_olog, ">>>Twist operator \n");
		int i, j;
		fprintf(fp_olog, ">>>Selecting a randomly chosen axis: \n%c\n>>>The rotation angle is \n%3.2f \n", axis, angletemp);
		fclose(fp_olog);
		double angle = angletemp * Pi / 180;
		//twist along x axis
		if (axis == 'x')
		{
			for (i = 0; i < nant_o[0]; i++) {
				xyz_e[i][0] = xyz_o[i][0];
				xyz_e[i][1] = xyz_o[i][1] * cos(angle * xyz_o[i][0]) + xyz_o[i][2] * sin(angle * xyz_o[i][0]);
				xyz_e[i][2] = -xyz_o[i][1] * sin(angle * xyz_o[i][0]) + xyz_o[i][2] * cos(angle * xyz_o[i][0]);
			}
		}
		//twsit along y axis
		if (axis == 'y')
		{
			for (i = 0; i < nant_o[0]; i++) {
				xyz_e[i][0] = xyz_o[i][0] * cos(angle * xyz_o[i][1]) - xyz_o[i][2] * sin(angle * xyz_o[i][1]);
				xyz_e[i][1] = xyz_o[i][1];
				xyz_e[i][2] = xyz_o[i][0] * sin(angle * xyz_o[i][1]) + xyz_o[i][2] * cos(angle * xyz_o[i][1]);
			}
		}
		//twist along z axis
		if (axis == 'z')
		{
			for (i = 0; i < nant_o[0]; i++) {
				xyz_e[i][0] = xyz_o[i][0] * cos(angle * xyz_o[i][2]) + xyz_o[i][1] * sin(angle * xyz_o[i][2]);
				xyz_e[i][1] = -xyz_o[i][0] * sin(angle * xyz_o[i][2]) + xyz_o[i][1] * cos(angle * xyz_o[i][2]);
				xyz_e[i][2] = xyz_o[i][2];
			}
		}
		for (i = 0; i < nant_o[0]; i++) {
			for (j = 0; j < 3; j++) {
				if (xyz_e[i][j] >= 1) { --xyz_e[i][j]; }
				if (xyz_e[i][j] < 0) { ++xyz_e[i][j]; }
			}
		}
	}

	void Rotation_command(
		double vec_o[3][3],     // lattice vectors
		double vec_e[3][3],
		double angle1,
		double angle2,
		double angle3
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		fprintf(fp_olog, ">>>Rotation operator \n");
		int i;
		double theta1;
		double theta2;
		double theta3;
		double vec_t[3][3];
		double vec_m[3][3];
		theta1 = angle1 * Pi / 180.00;
		theta2 = angle2 * Pi / 180.00;
		theta3 = angle3 * Pi / 180.00;
		fprintf(fp_olog, ">>>the rotation angles in x, y and z axis are \n%3.2f %3.2f %3.2f\n", angle1, angle2, angle3);
		//fprintf(fp_olog,">>> the rotation angles in x, y and z axis are %3.2f, %3.2f and %3.2f degree, respectively \n", theta1, theta2, theta3);
		fclose(fp_olog);
		//rotate along x axis
		for (i = 0; i < 3; i++) {
			vec_t[i][0] = vec_o[i][0];
			vec_t[i][1] = vec_o[i][1] * cos(theta1) + vec_o[i][2] * sin(theta1);
			vec_t[i][2] = -vec_o[i][1] * sin(theta1) + vec_o[i][2] * cos(theta1);
		}
		//rotate along y axis
		for (i = 0; i < 3; i++) {
			vec_m[i][0] = vec_t[i][0] * cos(theta2) - vec_t[i][2] * sin(theta2);
			vec_m[i][1] = vec_t[i][1];
			vec_m[i][2] = vec_t[i][0] * sin(theta2) + vec_t[i][2] * cos(theta2);
		}
		//rotate along z axis
		for (i = 0; i < 3; i++) {
			vec_e[i][0] = vec_m[i][0] * cos(theta3) + vec_m[i][1] * sin(theta3);
			vec_e[i][1] = -vec_m[i][0] * sin(theta3) + vec_m[i][1] * cos(theta3);
			vec_e[i][2] = vec_m[i][2];
		}
	}

	void Crossover_command(double& latt1_o,
		double& latt2_o,                  // Scaling Lattice parameter 
		double& latt_e,
		int    nant_o[],                    // total number of atoms and types
		int    typenum_o[],               // number of each type
		double vec1_o[3][3],
		double vec2_o[3][3],              // lattice vectors
		double vec_e[3][3],
		double xyz1_o[][3],
		double xyz2_o[][3],               // atomic coordinates
		double xyz_e[][3],
		char axis,
		double segment
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		fprintf(fp_olog, ">>> Crossover operator \n");
		fprintf(fp_olog, ">>>>>>>>>>>>>>>>>>>>>>> NOTE <<<<<<<<<<<<<<<<<<<<<<<< \n");
		fprintf(fp_olog, ">>>The POSCAR file for Crossover must have the same total number and type of atoms and the same sequence in POSCAR! \n");
		int i, j, k;
		latt_e = 1.0;
		for (i = 0; i < 3; i++) {
			for (j = 0; j < 3; j++) {
				vec_e[i][j] = (latt1_o * vec1_o[i][j] + latt2_o * vec2_o[i][j]) / 2.0000;
			}
		}
		fprintf(fp_olog, "Selecting a random axis: %c \n", axis); // select a random axis
		if (segment < 0.25 || segment>0.75)
		{
			fprintf(fp_olog, "segement is error!\n");
			return;
		}
		fprintf(fp_olog, "POSCAR_1: %1.2f; POSCAR_2: %1.2f \n", segment, 1.0000 - segment);
		int countnum;
		double temp[1000][3];
		int Nlower = 0;
		int Nuper = typenum_o[0];
		int zz;
		switch (axis)
		{
		case 'x': zz = 0; break;
		case 'y': zz = 1; break;
		case 'z': zz = 2; break;
		}
		for (i = 0; i < nant_o[1]; i++) {
			countnum = 0;
			for (j = Nlower; j < Nuper; j++) {
				if (xyz1_o[j][zz] < segment)
				{
					for (k = 0; k < 3; k++) { temp[countnum][k] = xyz1_o[j][k]; }
					countnum++;
				}
			}
			for (j = Nlower; j < Nuper; j++) {
				if (xyz2_o[j][zz] >= segment)
				{
					for (k = 0; k < 3; k++) { temp[countnum][k] = xyz2_o[j][k]; }
					countnum++;
				}
			}
			for (countnum; countnum < typenum_o[i]; countnum++) {
				temp[countnum][j] = RandDigit();
				fprintf(fp_olog, "%lf ", temp[countnum][j]);
			}
			for (j = Nlower; j < Nuper; j++) {
				for (k = 0; k < 3; k++) {
					xyz_e[j][k] = temp[j - Nlower][k];
				}
			}
			Nlower = Nlower + typenum_o[i];
			Nuper = Nuper + typenum_o[i + 1];
		}
		fclose(fp_olog);
	}

	void Stripple_command(int    nant_o[],                    // total number of atoms and types
		double vec_o[3][3],               // lattice vectors
		double vec_e[3][3],
		double xyz_o[][3],                // atomic coordinates
		double xyz_e[][3],
		double strain0, double strain1, double strain2,
		double strain3, double strain4, double strain5,
		int temp,
		double angle_1, double angle_2
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		fprintf(fp_olog, ">>>Stripple operator which merges strain with ripple \n");
		fclose(fp_olog);
		char mode[] = { "xyz" };
		//Strain operator
		Strain_command(mode, vec_o, vec_e, strain0, strain1, strain2, strain3, strain4, strain5);
		// Ripple opetator which shifts the coordinates of each atom long a random axis
		Ripple_command(nant_o, xyz_o, xyz_e, temp, angle_1, angle_2);
	}

	void Permustrain_command(int nant_o[],                    // total number of atoms and types
		int    typenum_o[],
		double vec_o[3][3],               // lattice vectors
		double vec_e[3][3],
		double xyz_o[][3],                // atomic coordinates
		double xyz_e[][3],
		double strain0, double strain1, double strain2,
		double strain3, double strain4, double strain5,
		int ExgNum,
		int num1[], int num2[]
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		fprintf(fp_olog, ">>>Permustrain operator which combines strain with exchange \n");
		fclose(fp_olog);
		char mode[] = { "xyz" };
		//Strain operator
		double strain[6];
		Strain_command(mode, vec_o, vec_e, strain0, strain1, strain2, strain3, strain4, strain5);
		//Exchange operator which exchanges the coordinates of two atoms of different types
		Exchange_command(nant_o, typenum_o, xyz_o, xyz_e, ExgNum, num1, num2);
	}

	void Rotstrain_command(double vec_o[3][3],               // lattice vectors
		double vec_e[3][3],
		double angle1,
		double angle2,
		double angle3,
		char mode[],
		double strain
	) {
		FILE* fp_olog = fopen("command_olog", "at+");
		fprintf(fp_olog, ">>>Rotstrain operator which combines rotation with strain \n");
		fclose(fp_olog);
		double vec_t[3][3];
		//Rotation operator
		Rotation_command(vec_o, vec_t, angle1, angle2, angle3);
		//Strain operator
	   /*if (!strcmp(mode, "xx"))
		   Strain_command(mode, vec_t, vec_e, strain, 0, 0, 0, 0, 0);
	   if (!strcmp(mode, "yy"))
		   Strain_command(mode, vec_t, vec_e, 0, strain, 0, 0, 0, 0);
	   if (!strcmp(mode, "zz"))
		   Strain_command(mode, vec_t, vec_e, 0, 0, strain, 0, 0, 0);*/
		if (!strcmp(mode, "xy"))
			Strain_command(mode, vec_t, vec_e, 0, 0, 0, strain, 0, 0);
		if (!strcmp(mode, "xz"))
			Strain_command(mode, vec_t, vec_e, 0, 0, 0, 0, strain, 0);
		if (!strcmp(mode, "yz"))
			Strain_command(mode, vec_t, vec_e, 0, 0, 0, 0, 0, strain);
		else
		{
			FILE* fp_olog = fopen("command_olog", "at+");
			fprintf(fp_olog, ">>>the mode of strain is wrong! \n");
			fclose(fp_olog);
		}
	}
}



