//mod�e 4 : g�m�rie sph�es molles

#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "time.h"


#undef N	

#include "parseciffile.h"
#include "sglite.h"
#include "SpaceGroup.h"


#define N 1000


#ifndef M_PI
	#define M_PI 3.1415926535897932384626433832795
#endif




double drand(){

	return ((double)rand())/(double)RAND_MAX;	// g��ateur de nombres al�toires
}


int main(int argc, char** argv, char** envp){


	double passage[3][3]; // matrice de passage de coordonn�s obliques en coordonn�s cart�iennes
	double inv_passage[3][3]; // // matrice de passage de coordonn�s obliques en coordonn�s cart�iennes inverse

	double alpha;  
	double beta;  //d�inition des angles
	double gamma;
	double a;
	double b;    // d�inition des param�res de maille
	double c;

	double coord_obliq_O[N][3];  // coordonn�s O dans le rep�e quelconque
	double coord_obliq_tmp_O[N][3];  // coordonn�s O temporaires dans le rep�e quelconque
	double coord_orthon_O[N][3];  // coordonn�s O dans le rep�e orthonorm�
	double coord_orthon_new_O[N][3]; // coordonn�s orthonorm�s des O avec Te �l'origine

	double coord_obliq_Te[N][3];  // coordonn�s Te dans le rep�e quelconque
	double coord_orthon_Te[N][3];	//coordonn�s Te dans le rep�e orthonorm�

	double coord_spher_lp[3]; //coordonn�s sph�iques lp
	double coord_spher_lp_min[3]; //coordonn�s sph�iques lp min
	double coord_orthon_new_lp[3]; //coordonn�s orthonorm�s lp min avec Te �l'origine
	double coord_orthon_lp[3]; //coordonn�s orthonorm�s lp min
	double coord_obliq_lp[3]; // coordonn�s obliques lp min
	double coord_spher_lp_min_opt[3]; // coordonn�s sph�iques lp min optims�s
				

	double dist_Olp[N];  //distance O-lp
	double dist2_Telp;	//distance Te-lp
	double K =0.004; //0.003085;	//constante de raideur de la force de rappel entre Te et lp        K = 0.0025;
	//double K = 1.00;	//constante de raideur de la force de rappel entre Te et lp

	double radius_O = 1.4;
	double radius_lp = 1.12;	//= 1.011;	// = 1.15; 
	double r_cut = radius_O + radius_lp ;
	double energy; //�ergie de r�ulsion
	double energy_min; // �ergie de r�ulsion minimale
	
	
	double p=0.33;	 //= 0.325; //facteur de Born-Mayer //0.345;
	
	int A = 1 ;	//constante pr�xponentielle
	int m;
	int i;
	int k; 
	int j;
	int V=10;

	double lp_pos_x; //position x de la paire libre
	double lp_pos_y; //position y de la paire libre
	double lp_pos_z; //position z de la paire libre
				
	double O_pos_x; //position x d'un oxyg�e
	double O_pos_y; //position y d'un oxyg�e
	double O_pos_z; //position z d'un oxyg�e
					
	double d_Olp_x; //diff�ence position O-lp en x
	double d_Olp_y; //diff�ence position O-lp en y
	double d_Olp_z; //diff�ence position O-lp en z




	if (argc != 2 && argc !=3 && argc !=4 ) {	//&& argc !=5) {
		printf("Usage : LPLocalisation 'ciffilename' <O-LP distance> <O-LP refined (0/1)>\n");
		return 0;
	}

	SCrystalInfo crystal;

	double LPOdist = 0.9;
	int isLPOrefined = 0;
	//int model = 1;

	if (argc > 3){
		LPOdist = atof(argv[2]);
		if (argc == 4){
			isLPOrefined = atoi(argv[3]);
		}
	//	if (argc == 5){
	//		model = atoi(argv[4]);
	//	}
	}
	
coord_spher_lp[0]= LPOdist;
coord_spher_lp_min[0] = LPOdist;

	// first load cif file...
	loadFileCIF(argv[1],&crystal);
	// second init spacegroup class		
	SpaceGroup spacegroup(crystal.space_group_HM);
	
	CrystMatrix_REAL coords;

	cout << "equivalent positions for Te :" << endl;
	for (i=0;i<crystal.nbte;i++){


//  #########################################################       	probleme coord obliq Te
		
	
		int mult = spacegroup.GetAllSymmetrics(crystal.atomTe[i].x,crystal.atomTe[i].y,crystal.atomTe[i].z,coords,false,false,true);
		cout << mult << "  " << crystal.atomTe[i].mult << endl;
		
			coord_obliq_Te[i][0]=coords(0,0);
			coord_obliq_Te[i][1]=coords(0,1);
			coord_obliq_Te[i][2]=coords(0,2);
	}		
	cout << coords(0,0) << endl;


	int nb_oxygen = 0; // multiples O des positions &#65533;uivalentes


	cout << "equivalent positions for O :" << endl;
	for (i=0;i<crystal.nbo;i++){
		int mult = spacegroup.GetAllSymmetrics(crystal.atomO[i].x,crystal.atomO[i].y,crystal.atomO[i].z,coords,false,false,true);
		cout << mult << "  " << crystal.atomO[i].mult << endl;
		
		for (j=0;j<mult;j++){
			coord_obliq_tmp_O[nb_oxygen][0]=coords(j,0);
			coord_obliq_tmp_O[nb_oxygen][1]=coords(j,1);
			coord_obliq_tmp_O[nb_oxygen][2]=coords(j,2);
			nb_oxygen++;
		}
	}

	cout << coords(0,0) << endl;


double aa , bb , cc ;
double u[3]; //diff&#65533;ence coordonn&#65533;s obliques Te - O en x,y,z
double metric_tensor[3][3]; //tenseur m&#65533;rique
double metric_tensor_u[3]; //tenseur m&#65533;rique * u
int nombre_O = 0; //nombre de O ayant une forte influence sur lp
int h;	//nombre_Te
double pos_O[1000][3];	// coordonn&#65533;s O 

double dist_TeO ;
double dist2_TeO ;

a = crystal.cell_a;
b = crystal.cell_b;
c = crystal.cell_c;
alpha = crystal.cell_alpha*M_PI/180;
beta = crystal.cell_beta*M_PI/180;
gamma = crystal.cell_gamma*M_PI/180;

// d&#65533;inition des coefficients du tenseur m&#65533;rique
metric_tensor[0][0] = a * a ;
metric_tensor[1][1] = b * b ;
metric_tensor[2][2] = c * c ;
metric_tensor[0][1] = metric_tensor[1][0] = a * b * cos(gamma) ;
metric_tensor[0][2] = metric_tensor[2][0] = a * c * cos(beta) ;
metric_tensor[1][2] = metric_tensor[2][1] = b * c * cos(alpha) ;


for (h=0 ; h < crystal.nbte ; h ++){

	for (m=0; m < nb_oxygen ; m ++){
	
		for (aa = -1.0 ; aa <= 1.0 ; aa += 1.0 ){
			for (bb = -1.0 ; bb <= 1.0 ; bb += 1.0 ){
				for (cc = -1.0 ; cc <= 1.0 ; cc += 1.0 ){
				
					pos_O[m][0] =  coord_obliq_tmp_O[m][0] + aa ;	
					pos_O[m][1] =  coord_obliq_tmp_O[m][1] + bb ;	
					pos_O[m][2] =  coord_obliq_tmp_O[m][2] + cc ;	
		
					u[0] = (coord_obliq_Te[h][0] - pos_O[m][0]);
					u[1] = (coord_obliq_Te[h][1] - pos_O[m][1]);
					u[2] = (coord_obliq_Te[h][2] - pos_O[m][2]);

					for ( i=0 ; i<3 ; i ++){
					
						metric_tensor_u[i] = 0;
									
						for (k=0; k<3; k++){
							metric_tensor_u[i] += metric_tensor[i][k] * u[k];
						}
					}

					dist2_TeO=0;

					for ( i=0; i<3 ; i++ ){

					dist2_TeO += u[i] * metric_tensor_u[i];
					}

					dist_TeO = 0;
					dist_TeO = sqrt(dist2_TeO);
					

					if (dist_TeO < 4.00){
						coord_obliq_O[nombre_O][0] = pos_O[m][0] ;
						coord_obliq_O[nombre_O][1] = pos_O[m][1] ;
						coord_obliq_O[nombre_O][2] = pos_O[m][2] ;
							/*	printf("dist Te-O = %6.4f %1.1i \t",dist_TeO,nombreO);
								printf("x = %6.4f \t",coordobliq[nombreO][0]);
								printf("y = %6.4f \t",coordobliq[nombreO][1]);
								printf("z = %6.4f \n",coordobliq[nombreO][2]);		*/					
						nombre_O++;
					}
				}
			}
		}
	}
			


			


// d&#65533;inition des coefficients de la matrice de passage
passage[1][0] = passage[2][0] = passage[2][1] = 0.0;
passage[0][0]= a;
passage[0][1]= b*cos(gamma);
passage[0][2]= c*cos(beta);
passage[1][1]= b*sin(gamma);
passage[1][2]= c*((cos(alpha)-cos(beta)*cos(gamma))/sin(gamma));
passage[2][2]= c*sqrt(sin(gamma)*sin(gamma)-(cos(alpha)*cos(alpha)+cos(beta)*cos(beta)-2*cos(alpha)*cos(beta)*cos(gamma)))/sin(gamma);

for (i=0;i<3;i++){

	coord_orthon_Te[h][i]=0.0;
	for (k=0;k<3;k++){

		coord_orthon_Te[h][i]+=passage[i][k]*coord_obliq_Te[h][k]; //coord orthonorm&#65533;s du Te
	}
}


for (m = 0 ; m < nombre_O ; m++){

	for (i=0;i<3;i++){

		coord_orthon_O[m][i]=0.0;
		for (k=0;k<3;k++){

		coord_orthon_O[m][i] += passage[i][k] * coord_obliq_O[m][k]; // coord orthonorm&#65533;s des O
		}
	}
	for (i=0;i<3;i++){

		coord_orthon_new_O[m][i]=0.0;
		coord_orthon_new_O[m][i]=coord_orthon_O[m][i] - coord_orthon_Te[h][i]; //coordonn&#65533;s orthonorm&#65533;s des O qd Te &#65533;l'origine
	}
}

energy_min = 1000.0;

//#####################

//for ( coord_spher_lp[0] = 0.8 ; coord_spher_lp[1] <= 1.3 ; coord_spher_lp[1] += 0.05){

for (coord_spher_lp[1] = 0. ; coord_spher_lp[1] <= M_PI ; coord_spher_lp[1] += 0.005){

		for (coord_spher_lp[2] = 0. ; coord_spher_lp[2] <= 2*M_PI; coord_spher_lp[2] += 0.005){

			energy = 0.0;

			lp_pos_x = coord_spher_lp[0]*sin(coord_spher_lp[1])*cos(coord_spher_lp[2]);
			lp_pos_y = coord_spher_lp[0]*sin(coord_spher_lp[1])*sin(coord_spher_lp[2]);
			lp_pos_z = coord_spher_lp[0]*cos(coord_spher_lp[1]);

			dist2_Telp = lp_pos_x * lp_pos_x + lp_pos_y * lp_pos_y + lp_pos_z * lp_pos_z ;
			energy = 0.5*K * dist2_Telp;

			for (m = 0 ; m < nombre_O ; m++){

				 O_pos_x = coord_orthon_new_O[m][0];
				 O_pos_y = coord_orthon_new_O[m][1];
				 O_pos_z = coord_orthon_new_O[m][2];
				
				 d_Olp_x = O_pos_x - lp_pos_x;
				 d_Olp_y = O_pos_y - lp_pos_y;
				 d_Olp_z = O_pos_z - lp_pos_z;
				
				dist_Olp[m] = sqrt(d_Olp_x*d_Olp_x+d_Olp_y*d_Olp_y+d_Olp_z*d_Olp_z);
	
				if (dist_Olp[m]< r_cut)
					energy += A *(exp(-dist_Olp[m]/p)-exp(-r_cut/p));
				
			}
		
			if (energy <= energy_min){

				energy_min = energy;
//				coord_spher_lp_min[0] = coord_spher_lp[0];
				coord_spher_lp_min[1] = coord_spher_lp[1];
				coord_spher_lp_min[2] = coord_spher_lp[2];
				}
			
		}
	} 
//}

//#####################

double delta_r = 0.005;
double delta_theta = 0.001;
double delta_phi = 0.001;

srand((unsigned)time(NULL));
srand(1000);

coord_spher_lp_min_opt[0] = coord_spher_lp_min[0];
coord_spher_lp_min_opt[1] = coord_spher_lp_min[1];
coord_spher_lp_min_opt[2] = coord_spher_lp_min[2];

	
double	nb_cycles = 1000000.0;

printf("\n \nenergie sans optimisation = %9.10f \n", energy_min);

//goto lab1;

energy_min = 10000.0;
if (isLPOrefined == 1){

	for ( i = 0; i < nb_cycles ; i++){

				coord_spher_lp_min[0] = coord_spher_lp_min_opt[0] + delta_r*(drand()-0.5)*2.0;
				coord_spher_lp_min[1] = coord_spher_lp_min_opt[1] + delta_theta*(drand()-0.5)*2.0;
				coord_spher_lp_min[2] = coord_spher_lp_min_opt[2] + delta_phi*(drand()-0.5)*2.0;
				energy=0.0;
			
			lp_pos_x = coord_spher_lp_min[0]*sin(coord_spher_lp_min[1])*cos(coord_spher_lp_min[2]);
			lp_pos_y = coord_spher_lp_min[0]*sin(coord_spher_lp_min[1])*sin(coord_spher_lp_min[2]);
			lp_pos_z = coord_spher_lp_min[0]*cos(coord_spher_lp_min[1]);

			dist2_Telp = lp_pos_x * lp_pos_x + lp_pos_y * lp_pos_y + lp_pos_z * lp_pos_z ;
			energy =0.5* K * dist2_Telp;

			for (m = 0 ; m < nombre_O ; m++){

				 O_pos_x = coord_orthon_new_O[m][0];
				 O_pos_y = coord_orthon_new_O[m][1];
				 O_pos_z = coord_orthon_new_O[m][2];
				
				 d_Olp_x = O_pos_x - lp_pos_x;
				 d_Olp_y = O_pos_y - lp_pos_y;
				 d_Olp_z = O_pos_z - lp_pos_z;
				
				dist_Olp[m] = sqrt(d_Olp_x*d_Olp_x+d_Olp_y*d_Olp_y+d_Olp_z*d_Olp_z);

				if (dist_Olp[m]< r_cut)
					energy+= A *(exp(-dist_Olp[m]/p)-exp(-r_cut/p));
				
			}
		if (energy < energy_min && coord_spher_lp_min[0] < V){
			energy_min = energy;
			coord_spher_lp_min_opt[0] = coord_spher_lp_min[0];
			coord_spher_lp_min_opt[1] = coord_spher_lp_min[1];
			coord_spher_lp_min_opt[2] = coord_spher_lp_min[2];
		}
	}
}
else{
	
	for ( i = 0; i < nb_cycles ; i++){

			
				coord_spher_lp_min[1] = coord_spher_lp_min_opt[1] + delta_theta*(drand()-0.5)*2.0;
				coord_spher_lp_min[2] = coord_spher_lp_min_opt[2] + delta_phi*(drand()-0.5)*2.0;
				energy=0.0;
			
			lp_pos_x = coord_spher_lp[0]*sin(coord_spher_lp[1])*cos(coord_spher_lp[2]);
			lp_pos_y = coord_spher_lp[0]*sin(coord_spher_lp[1])*sin(coord_spher_lp[2]);
			lp_pos_z = coord_spher_lp[0]*cos(coord_spher_lp[1]);

		dist2_Telp = lp_pos_x * lp_pos_x + lp_pos_y * lp_pos_y + lp_pos_z * lp_pos_z ;
			energy =0.5* K * dist2_Telp;

			for (m = 0 ; m < nombre_O ; m++){

				 O_pos_x = coord_orthon_new_O[m][0];
				 O_pos_y = coord_orthon_new_O[m][1];
				 O_pos_z = coord_orthon_new_O[m][2];
				
				 d_Olp_x = O_pos_x - lp_pos_x;
				 d_Olp_y = O_pos_y - lp_pos_y;
				 d_Olp_z = O_pos_z - lp_pos_z;
				
				dist_Olp[m] = sqrt(d_Olp_x*d_Olp_x+d_Olp_y*d_Olp_y+d_Olp_z*d_Olp_z);

				if (dist_Olp[m]< r_cut)
				energy += A*(exp(-dist_Olp[m]/p)-exp(-r_cut/p));
				
			}	

		if (energy < energy_min){
			energy_min = energy;
			coord_spher_lp_min_opt[1] = coord_spher_lp_min[1];
			coord_spher_lp_min_opt[2] = coord_spher_lp_min[2];
		}
	}
}

//###############################################################################################################
// affichage de la distance et de l'&#65533;ergie entre la paire libre et chaque oxygene 

			lp_pos_x = coord_spher_lp_min_opt[0]*sin(coord_spher_lp_min_opt[1])*cos(coord_spher_lp_min_opt[2]);
			lp_pos_y = coord_spher_lp_min_opt[0]*sin(coord_spher_lp_min_opt[1])*sin(coord_spher_lp_min_opt[2]);
			lp_pos_z = coord_spher_lp_min_opt[0]*cos(coord_spher_lp_min_opt[1]);
			dist2_Telp = lp_pos_x * lp_pos_x + lp_pos_y * lp_pos_y + lp_pos_z * lp_pos_z ;

			for (m = 0 ; m < nombre_O ; m++){

				energy =0.5* K * dist2_Telp;
				 O_pos_x = coord_orthon_new_O[m][0];
				 O_pos_y = coord_orthon_new_O[m][1];
				 O_pos_z = coord_orthon_new_O[m][2];
				
				 d_Olp_x = O_pos_x - lp_pos_x;
				 d_Olp_y = O_pos_y - lp_pos_y;
				 d_Olp_z = O_pos_z - lp_pos_z;
				
				dist_Olp[m] = sqrt(d_Olp_x*d_Olp_x+d_Olp_y*d_Olp_y+d_Olp_z*d_Olp_z);


				if (dist_Olp[m] < 4.5){
					if (dist_Olp[m]<= r_cut){
						energy = A*(exp(-dist_Olp[m]/p)-exp(-r_cut/p)) + 0.5* K * dist2_Telp;
					}
				printf(" dist_Olp = %5.3f \t", dist_Olp[m]);
				printf(" energie = %9.10f \n", energy);
				}
			}

//###############################################################################################################

//lab1:

printf("energie avec optimisation = %9.10f \n \n", energy_min);



coord_orthon_new_lp[0] = coord_spher_lp_min_opt[0] * sin(coord_spher_lp_min_opt[1]) * cos(coord_spher_lp_min_opt[2]); //
coord_orthon_new_lp[1] = coord_spher_lp_min_opt[0] * sin(coord_spher_lp_min_opt[1]) * sin(coord_spher_lp_min_opt[2]); // coord orthonorm&#65533;s de la paire libre qd Te &#65533;l'origine
coord_orthon_new_lp[2] = coord_spher_lp_min_opt[0] * cos(coord_spher_lp_min_opt[1]);								  //


for (i=0; i<3 ; i++){

	coord_orthon_lp[i]=0;
	coord_orthon_lp[i]= coord_orthon_new_lp[i] + coord_orthon_Te[h][i]; // coord orthonorm&#65533;s de la paire libre
}


// d&#65533;inition des coefficients de la matrice de passage inverse
inv_passage[1][0] = inv_passage[2][0] = inv_passage[2][1] = 0;
inv_passage[0][0]=1/passage[0][0];
inv_passage[1][1]=1/passage[1][1];
inv_passage[2][2]=1/passage[2][2];
inv_passage[0][1]=-passage[0][1]*inv_passage[0][0]*inv_passage[1][1];
inv_passage[1][2]=-passage[1][2]*inv_passage[1][1]*inv_passage[2][2];
inv_passage[0][2]=-(passage[0][1]*inv_passage[1][2]+passage[0][2]*inv_passage[2][2])*inv_passage[0][0];


for (i=0;i<3;i++){

	coord_obliq_lp[i]=0.0;
	for (k=0;k<3;k++){

		coord_obliq_lp[i]+=inv_passage[i][k]*coord_orthon_lp[k]; //coord obliques de la paire libre
	}
}


printf("pour l'atome de Te place en \n");
printf(" %6.4f \n", coord_obliq_Te[h][0]);
printf(" %6.4f \n", coord_obliq_Te[h][1]);
printf(" %6.4f \n", coord_obliq_Te[h][2]);

printf("coordonnees de la paire libre : \n");
printf(" %6.4f \n", coord_obliq_lp[0]);
printf(" %6.4f \n", coord_obliq_lp[1]);
printf(" %6.4f \n", coord_obliq_lp[2]);

printf("longueur de la liaison Te-lp = \n %5.3f \n", coord_spher_lp_min_opt[0]);

}

	return 0;

}