scemama/StageYann
1
0
Fork 0
Code Pull Requests Releases Activity

Work : test BOYS

Browse Source
This commit is contained in:
Yann Damour 2020-06-02 10:45:02 +02:00
parent 1fcf7d2e98
commit ed290c7f44
4 changed files with 141 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
Boys.ipynb
View File

@ -75,7 +75,7 @@
"{\\cal B}_2= \\sum_{i=1}^N \\big[ \\langle x^2 \\rangle_i - \\langle x \\rangle^2_i + \\langle y^2 \\rangle_i - …\\langle y \\rangle^2_i + \\langle z^2 \\rangle_i - \\langle z \\rangle^2_i \\big]\n", "{\\cal B}_2= \\sum_{i=1}^N \\big[ \\langle x^2 \\rangle_i - \\langle x \\rangle^2_i + \\langle y^2 \\rangle_i - …\\langle y \\rangle^2_i + \\langle z^2 \\rangle_i - \\langle z \\rangle^2_i \\big]\n",
"$$\n", "$$\n",
"$$\n", "$$\n",
"{\\cal B}_2 = \\sum_{i=1}^N \\big[ \\langle x^4 \\rangle_i - 4 \\langle x^3 \\rangle_i \\langle x \\rangle_i\n", "{\\cal B}_4= \\sum_{i=1}^N \\big[ \\langle x^4 \\rangle_i - 4 \\langle x^3 \\rangle_i \\langle x \\rangle_i\n",
" + 6 \\langle x^2 \\rangle_i \\langle x \\rangle^2_i\n", " + 6 \\langle x^2 \\rangle_i \\langle x \\rangle^2_i\n",
"- 3 \\langle x \\rangle^4_i \\big] + \\big[ ...y...] + \\big[ ...z...] \n", "- 3 \\langle x \\rangle^4_i \\big] + \\big[ ...y...] + \\big[ ...z...] \n",
"$$\n", "$$\n",
@ -236,11 +236,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"metadata": { "metadata": {},
"jupyter": {
"outputs_hidden": true
}
},
"outputs": [], "outputs": [],
"source": [ "source": [
"(*\n", "(*\n",
@ -251,7 +247,7 @@
"*)\n", "*)\n",
"\n", "\n",
"(* Construction de la matrice X n par n *)\n", "(* Construction de la matrice X n par n *)\n",
"let n =5 \n", "let n =2\n",
"let ran=Mat.random ~range:1. ~from:0. n n;; (* Construction d'une matrice random n*n *)\n", "let ran=Mat.random ~range:1. ~from:0. n n;; (* Construction d'une matrice random n*n *)\n",
"\n", "\n",
"(* Antisymétrisation de la matrice *)\n", "(* Antisymétrisation de la matrice *)\n",
@ -270,11 +266,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"metadata": { "metadata": {},
"jupyter": {
"outputs_hidden": true
}
},
"outputs": [], "outputs": [],
"source": [ "source": [
"(* Mise au carré de la matrice : X -> X² *)\n", "(* Mise au carré de la matrice : X -> X² *)\n",
@ -294,11 +286,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"metadata": { "metadata": {},
"jupyter": {
"outputs_hidden": true
}
},
"outputs": [], "outputs": [],
"source": [ "source": [
"(* Passage de -tau² à tau² *)\n", "(* Passage de -tau² à tau² *)\n",
@ -314,11 +302,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"metadata": { "metadata": {},
"jupyter": {
"outputs_hidden": true
}
},
"outputs": [], "outputs": [],
"source": [ "source": [
"(*Test opération sur tau, fonction générale*)\n", "(*Test opération sur tau, fonction générale*)\n",
@ -331,11 +315,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"metadata": { "metadata": {},
"jupyter": {
"outputs_hidden": true
}
},
"outputs": [], "outputs": [],
"source": [ "source": [
"(* Calcul de cos tau à partir du vecteur tau *)\n", "(* Calcul de cos tau à partir du vecteur tau *)\n",
@ -387,7 +367,11 @@
"let b = gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m)));;\n", "let b = gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m)));;\n",
"\n", "\n",
"(* Somme de 2 matrices, ici -> a + b -> R *)\n", "(* Somme de 2 matrices, ici -> a + b -> R *)\n",
"let m_r = Mat.add a b;; (* Matrice de rotation R *)" "let m_r = Mat.add a b;; (* Matrice de rotation R *)\n",
"\n",
"gesvd m_r;;\n",
"\n",
"\n"
] ]
}, },
{ {
@ -494,11 +478,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"metadata": { "metadata": {},
"jupyter": {
"outputs_hidden": true
}
},
"outputs": [], "outputs": [],
"source": [ "source": [
"let a= 1.;;\n", "let a= 1.;;\n",

153
Work.ipynb
View File

@ -1338,11 +1338,7 @@
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,
"metadata": { "metadata": {},
"jupyter": {
"source_hidden": true
}
},
"outputs": [], "outputs": [],
"source": [ "source": [
"(*\n", "(*\n",
@ -1490,10 +1486,9 @@
" ),Vec.sum v_d);;\n", " ),Vec.sum v_d);;\n",
"\n", "\n",
"(*********************)\n", "(*********************)\n",
"(*\n", "\n",
"f_alpha m_C;;\n", "f_alpha m_C;;\n",
"let m_alpha , s_D = f_alpha m_C;;\n", "let m_alpha , s_D = f_alpha m_C;;\n"
"*)"
] ]
}, },
{ {
@ -1539,9 +1534,73 @@
" Vec.sum(Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)));;\n", " Vec.sum(Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)));;\n",
"\n", "\n",
"(****************************)\n", "(****************************)\n",
"(*\n", "\n",
"f_alpha_boys m_C;;\n", "f_alpha_boys m_C;;\n"
"*)" ]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Test méthode de boys *)\n",
"\n",
"let f_theta m_C = \n",
" let n_mo = Mat.dim2 m_C\n",
" in\n",
" let phi_x_phi =\n",
" Multipole.matrix_x multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
" in \n",
" let phi_y_phi =\n",
" Multipole.matrix_y multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
" in\n",
" let phi_z_phi =\n",
" Multipole.matrix_z multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
" in\n",
" let phi_x2_phi =\n",
" Multipole.matrix_x2 multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
" in \n",
" let phi_y2_phi =\n",
" Multipole.matrix_y2 multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
" in\n",
" let phi_z2_phi =\n",
" Multipole.matrix_z2 multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
" in\n",
" let m_b_0 =\n",
" let b_0 g h = Mat.init_cols n_mo n_mo (fun i j ->\n",
" h.{i,i} +. h.{j,j} -. (g.{i,i}**2. +. g.{j,j}**2.))\n",
" in Mat.add (b_0 phi_x_phi phi_x2_phi) (Mat.add (b_0 phi_y_phi phi_y2_phi) (b_0 phi_z_phi phi_z2_phi)) \n",
" in\n",
" let m_beta = \n",
" let beta g = Mat.init_cols n_mo n_mo (fun i j ->\n",
" (g.{i,i} -. g.{j,j})**2. -. 4. *. g.{i,j}**2.)\n",
" in Mat.add (beta phi_x_phi) (Mat.add (beta phi_y_phi) (beta phi_z_phi))\n",
" (*in Util.debug_matrix \"m_beta\" m_beta;*)\n",
" in\n",
" let m_gamma = \n",
" let gamma g = Mat.init_cols n_mo n_mo (fun i j ->\n",
" 4. *. g.{i,j} *. (g.{i,i} -. g.{j,j}) )\n",
" in Mat.add (gamma phi_x_phi) (Mat.add (gamma phi_y_phi) (gamma phi_z_phi))\n",
" (*in Util.debug_matrix \"m_gamma\" m_gamma;*)\n",
" \n",
" in\n",
" let m_theta = Mat.init_cols n_mo n_mo (fun i j ->\n",
" -. 0.25 *. atan(m_gamma.{i,j} /. m_beta.{i,j}))\n",
" in\n",
" let m_critere_B = Mat.init_cols n_mo n_mo (fun i j ->\n",
" m_b_0.{i,j} +. 0.25 *. ((1. -. cos(4. *. m_theta.{i,j})) *. m_beta.{i,j} +. sin(4. *. m_theta.{i,j}) *. m_gamma.{i,j}))\n",
" in m_theta, Vec.sum(Mat.as_vec m_critere_B)\n",
";;\n",
"\n",
"f_theta m_C;;\n",
"f_theta m_B;;"
] ]
}, },
{ {
@ -1562,6 +1621,7 @@
" match methode with \n", " match methode with \n",
" | \"Boys\"\n", " | \"Boys\"\n",
" | \"boys\" -> let alpha_boys , d_boys = f_alpha_boys m_C in {m_alpha = alpha_boys; d = d_boys}\n", " | \"boys\" -> let alpha_boys , d_boys = f_alpha_boys m_C in {m_alpha = alpha_boys; d = d_boys}\n",
" | \"BOYS\" -> let theta_boys, b_boys = f_theta m_C in {m_alpha = theta_boys; d = b_boys}\n",
" | \"ER\"\n", " | \"ER\"\n",
" | \"er\" -> let alpha_er , d_er = f_alpha m_C in {m_alpha = alpha_er; d = d_er}\n", " | \"er\" -> let alpha_er , d_er = f_alpha m_C in {m_alpha = alpha_er; d = d_er}\n",
" | _ -> invalid_arg \"Unknown method, please enter Boys or ER\"\n", " | _ -> invalid_arg \"Unknown method, please enter Boys or ER\"\n",
@ -1624,9 +1684,9 @@
" if n_rec_alpha == 0 \n", " if n_rec_alpha == 0 \n",
" then m_alpha \n", " then m_alpha \n",
" else Mat.init_cols n_mo n_mo (fun i j -> \n", " else Mat.init_cols n_mo n_mo (fun i j -> \n",
" if (m_alpha.{i,j}) > (pi /. 2.) \n", " if (m_alpha.{i,j}) >= (pi /. 2.) \n",
" then (m_alpha.{i,j} -. ( pi /. 2.))\n", " then (m_alpha.{i,j} -. ( pi /. 2.))\n",
" else if m_alpha.{i,j} < -. pi /. 2.\n", " else if m_alpha.{i,j} <= -. pi /. 2.\n",
" then (m_alpha.{i,j} +. ( pi /. 2.))\n", " then (m_alpha.{i,j} +. ( pi /. 2.))\n",
" else if m_alpha.{i,j} < 0. \n", " else if m_alpha.{i,j} < 0. \n",
" then -. m_alpha.{i,j}\n", " then -. m_alpha.{i,j}\n",
@ -1696,11 +1756,6 @@
"(*\n", "(*\n",
"let alpha = alphaij.alpha_max;; (* Fonction -> constante *) \n", "let alpha = alphaij.alpha_max;; (* Fonction -> constante *) \n",
"f_R alpha;;\n", "f_R alpha;;\n",
"*)\n",
"(*\n",
"alpha_v \"deloc\" alphaij;;\n",
"let alpha = (alpha_v \"loc\" alphaij) ;;\n",
"f_R alpha ;;\n",
"*)" "*)"
] ]
}, },
@ -1750,7 +1805,7 @@
"\n", "\n",
"(*************************)\n", "(*************************)\n",
"(*\n", "(*\n",
"let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi;; \n", "let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi m_C;; \n",
"*)" "*)"
] ]
}, },
@ -1774,8 +1829,8 @@
" else 0.)\n", " else 0.)\n",
"(*************************)\n", "(*************************)\n",
"(*\n", "(*\n",
"let m_Psi = f_Psi m_Ksi indice_i indice_j;; \n", "let m_Psi = f_k m_Ksi indice_i indice_j m_C;; \n",
"let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j;; \n", "let m_Psi_tilde = f_k m_Ksi_tilde indice_i indice_j m_C;; \n",
"*)" "*)"
] ]
}, },
@ -1902,12 +1957,12 @@
" (* Critères possibles de localisation -> norme_alpha, moyenne_alpha, alpha, diff *)\n", " (* Critères possibles de localisation -> norme_alpha, moyenne_alpha, alpha, diff *)\n",
" \n", " \n",
" (* Norme de m_alpha, moyenne et max *)\n", " (* Norme de m_alpha, moyenne et max *)\n",
" let norme_alpha = norme m_alpha\n", " (*let norme_alpha = norme m_alpha\n",
" in\n", " in\n",
" let moyenne_alpha = norme_alpha /. (float_of_int((Mat.dim1 m_C)* (Mat.dim2 m_C)))\n", " let moyenne_alpha = norme_alpha /. (float_of_int((Mat.dim1 m_C)* (Mat.dim2 m_C)))\n",
" in\n", " in\n",
" let alpha = alpha_max /. epsilon\n", " let alpha = alpha_max /. epsilon\n",
" in\n", " in*)\n",
" \n", " \n",
" (* D critère à maximiser *)\n", " (* D critère à maximiser *)\n",
" let critere_D = alphad.d \n", " let critere_D = alphad.d \n",
@ -1916,7 +1971,7 @@
" in\n", " in\n",
" \n", " \n",
" (*Printf.printf \"%f\\n%!\" diff;*)\n", " (*Printf.printf \"%f\\n%!\" diff;*)\n",
" Printf.printf \"%i %f %f %f %f %f\\n%!\" n max_D critere_D alpha norme_alpha moyenne_alpha;\n", " (*Printf.printf \"%i %f %f %f %f %f\\n%!\" n max_D critere_D alpha norme_alpha moyenne_alpha;*)\n",
" \n", " \n",
" (* Stabilisation de la convergence *)\n", " (* Stabilisation de la convergence *)\n",
"\n", "\n",
@ -1925,18 +1980,18 @@
" then max_D\n", " then max_D\n",
" else critere_D\n", " else critere_D\n",
" in\n", " in\n",
" \n", " (*\n",
" if diff > 0. && epsilon > 0.001\n", " if diff > 0. && epsilon > 0.001\n",
" then localisation m_C methode (epsilon /. 1000.) (n-1) critere_D prev_m_alpha cc max_D\n", " then localisation m_C methode (epsilon /. 1000.) (n-1) critere_D prev_m_alpha cc max_D\n",
" else let epsilon = 1.\n", " else let epsilon = 1.\n",
" in\n", " in\n",
" \n", " *)\n",
" (*Util.debug_matrix \"new_alpha_m\" m_alpha;*)\n", " (*Util.debug_matrix \"new_alpha_m\" m_alpha;*)\n",
" (*Util.debug_matrix \"m_new_m_C\" m_new_m_C;*)\n", " (*Util.debug_matrix \"m_new_m_C\" m_new_m_C;*)\n",
"\n", "\n",
"\n", "\n",
" (* Critère de conv pour sortir *)\n", " (* Critère de conv pour sortir *)\n",
" if diff**2. < cc**2.\n", " if alpha_max**2. +. 1. < cc**2.\n",
" then m_new_m_C\n", " then m_new_m_C\n",
" else localisation m_new_m_C methode epsilon (n-1) critere_D m_alpha cc max_D;;" " else localisation m_new_m_C methode epsilon (n-1) critere_D m_alpha cc max_D;;"
] ]
@ -1958,7 +2013,7 @@
"in\n", "in\n",
"localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D;;\n", "localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D;;\n",
" \n", " \n",
"localise localisation m_C \"Boys\" 1. 1000 10e-5;;" "let m_B = localise localisation m_C \"BOYS\" 1. 100 10e-10;;"
] ]
}, },
{ {
@ -2201,6 +2256,48 @@
"*)\n" "*)\n"
] ]
}, },
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let a = [1;2;3;4;5];;\n",
"List.iter print_int a;;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let print_list f lst =\n",
" let rec print_elements = function\n",
" | [] -> ()\n",
" | h::t -> f h; print_string \";\"; print_elements t\n",
" in\n",
" print_string \"[\";\n",
" print_elements lst;\n",
" print_string \"]\";;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print_list print_int [3;6;78;5;2;34;7];;"
]
},
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,

23
Work_test.ipynb
View File

@ -1053,9 +1053,8 @@
" Vec.sum(Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)));;\n", " Vec.sum(Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)));;\n",
"\n", "\n",
"(****************************)\n", "(****************************)\n",
"(*\n", "\n",
"f_alpha_boys m_C;;\n", "f_alpha_boys m_C;;\n"
"*)"
] ]
}, },
{ {
@ -1094,24 +1093,6 @@
"*)" "*)"
] ]
}, },
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Test choix cc *)\n",
"\n",
"let choix_cc cc = \n",
" let toto cc = \n",
" match cc with \n",
" | \"alpha_max\" -> let toto = alpha_max /. epsilon in {tutu = toto}\n",
" | \"norme_alpha\" -> let toto = norme_alpha in {tutu = toto}\n",
" | \"critere_D\" -> let toto = critere_D in {tutu = toto}\n",
" | _ -> invalid_arg \"Unknown method, please enter alpha_max\"\n",
"in toto_cc;;\n"
]
},
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,

2
run_cyanine.ipynb
View File

@ -69,7 +69,7 @@
"source": [ "source": [
"(*I : Position et base atomique*)\n", "(*I : Position et base atomique*)\n",
"\n", "\n",
"(* 1. Coordonnées *)\n", "(* 1. Coordonnées *) (* !!! Attention à bien mettre \"nbatome sur la même ligne !!! *)\n",
"let charge = 1;;\n", "let charge = 1;;\n",
"let xyz_string = \n", "let xyz_string = \n",
"\"8\n", "\"8\n",