rendered paste body//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 �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�ence coordonn�s obliques Te - O en x,y,z
double metric_tensor[3][3]; //tenseur m�rique
double metric_tensor_u[3]; //tenseur m�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�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�inition des coefficients du tenseur m�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�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�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�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�s orthonorm�s des O qd Te �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'�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�s de la paire libre qd Te �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�s de la paire libre
}
// d�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;
}