From ed290c7f4400f8fa70ec184eabd66c58de260303 Mon Sep 17 00:00:00 2001 From: yann Date: Tue, 2 Jun 2020 10:45:02 +0200 Subject: [PATCH] Work : test BOYS --- Boys.ipynb | 46 ++++---------- Work.ipynb | 153 +++++++++++++++++++++++++++++++++++++--------- Work_test.ipynb | 23 +------ run_cyanine.ipynb | 2 +- 4 files changed, 141 insertions(+), 83 deletions(-) diff --git a/Boys.ipynb b/Boys.ipynb index 4c1b967..d3c134e 100644 --- a/Boys.ipynb +++ b/Boys.ipynb @@ -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", "$$\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", "- 3 \\langle x \\rangle^4_i \\big] + \\big[ ...y...] + \\big[ ...z...] \n", "$$\n", @@ -236,11 +236,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "jupyter": { - "outputs_hidden": true - } - }, + "metadata": {}, "outputs": [], "source": [ "(*\n", @@ -251,7 +247,7 @@ "*)\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", "\n", "(* Antisymétrisation de la matrice *)\n", @@ -270,11 +266,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "jupyter": { - "outputs_hidden": true - } - }, + "metadata": {}, "outputs": [], "source": [ "(* Mise au carré de la matrice : X -> X² *)\n", @@ -294,11 +286,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "jupyter": { - "outputs_hidden": true - } - }, + "metadata": {}, "outputs": [], "source": [ "(* Passage de -tau² à tau² *)\n", @@ -314,11 +302,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "jupyter": { - "outputs_hidden": true - } - }, + "metadata": {}, "outputs": [], "source": [ "(*Test opération sur tau, fonction générale*)\n", @@ -331,11 +315,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "jupyter": { - "outputs_hidden": true - } - }, + "metadata": {}, "outputs": [], "source": [ "(* 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", "\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", "execution_count": null, - "metadata": { - "jupyter": { - "outputs_hidden": true - } - }, + "metadata": {}, "outputs": [], "source": [ "let a= 1.;;\n", diff --git a/Work.ipynb b/Work.ipynb index 5d1ad77..8b2e2e9 100644 --- a/Work.ipynb +++ b/Work.ipynb @@ -1338,11 +1338,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "jupyter": { - "source_hidden": true - } - }, + "metadata": {}, "outputs": [], "source": [ "(*\n", @@ -1490,10 +1486,9 @@ " ),Vec.sum v_d);;\n", "\n", "(*********************)\n", - "(*\n", + "\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", "\n", "(****************************)\n", - "(*\n", - "f_alpha_boys m_C;;\n", - "*)" + "\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", " | \"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\" -> 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", @@ -1624,9 +1684,9 @@ " if n_rec_alpha == 0 \n", " then m_alpha \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", - " 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", " else if m_alpha.{i,j} < 0. \n", " then -. m_alpha.{i,j}\n", @@ -1696,11 +1756,6 @@ "(*\n", "let alpha = alphaij.alpha_max;; (* Fonction -> constante *) \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", - "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", "(*************************)\n", "(*\n", - "let m_Psi = f_Psi m_Ksi indice_i indice_j;; \n", - "let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j;; \n", + "let m_Psi = f_k m_Ksi indice_i indice_j m_C;; \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", " \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", " let moyenne_alpha = norme_alpha /. (float_of_int((Mat.dim1 m_C)* (Mat.dim2 m_C)))\n", " in\n", " let alpha = alpha_max /. epsilon\n", - " in\n", + " in*)\n", " \n", " (* D critère à maximiser *)\n", " let critere_D = alphad.d \n", @@ -1916,7 +1971,7 @@ " in\n", " \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", " (* Stabilisation de la convergence *)\n", "\n", @@ -1925,18 +1980,18 @@ " then max_D\n", " else critere_D\n", " in\n", - " \n", + " (*\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", " else let epsilon = 1.\n", " in\n", - " \n", + " *)\n", " (*Util.debug_matrix \"new_alpha_m\" m_alpha;*)\n", " (*Util.debug_matrix \"m_new_m_C\" m_new_m_C;*)\n", "\n", "\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", " else localisation m_new_m_C methode epsilon (n-1) critere_D m_alpha cc max_D;;" ] @@ -1958,7 +2013,7 @@ "in\n", "localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D;;\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" ] }, + { + "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", "execution_count": null, diff --git a/Work_test.ipynb b/Work_test.ipynb index f36ff51..b30881c 100644 --- a/Work_test.ipynb +++ b/Work_test.ipynb @@ -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", "\n", "(****************************)\n", - "(*\n", - "f_alpha_boys m_C;;\n", - "*)" + "\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", "execution_count": null, diff --git a/run_cyanine.ipynb b/run_cyanine.ipynb index 54ec556..537b3b7 100644 --- a/run_cyanine.ipynb +++ b/run_cyanine.ipynb @@ -69,7 +69,7 @@ "source": [ "(*I : Position et base atomique*)\n", "\n", - "(* 1. Coordonnées *)\n", + "(* 1. Coordonnées *) (* !!! Attention à bien mettre \"nbatome sur la même ligne !!! *)\n", "let charge = 1;;\n", "let xyz_string = \n", "\"8\n",