From 26507ec7a0dfc1d7b392c56c8cf09f63ddefa977 Mon Sep 17 00:00:00 2001 From: yann Date: Mon, 4 May 2020 14:42:00 +0200 Subject: [PATCH] Work: modif boucle test alpha\n", + " let b12 g = Mat.init_cols n_ao n_ao ( fun i j ->\n", " (g.{i,i} -. g.{j,j}) *. g.{i,j})\n", "\n", " in \n", " Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi))\n", " in \n", " let m_a12 =\n", - " let a12 g = Mat.init_cols n_mo n_mo ( fun i j -> \n", + " let a12 g = Mat.init_cols n_ao n_ao ( fun i j -> \n", " g.{i,j} *. g.{i,j} \n", " -. (1. /. 4.) *. (g.{i,i} -. g.{j,j} *. g.{i,i} -. g.{j,j}))\n", " in\n", " Mat.add (a12 phi_x_phi) ( Mat.add (a12 phi_y_phi) (a12 phi_z_phi))\n", " in\n", - " Mat.init_cols n_mo n_mo ( fun i j ->\n", + " Mat.init_cols n_ao n_ao ( fun i j ->\n", " asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )));;\n", - "\n", - "\n", - "(* let m_a12 =\n", - " \n", - " let m_a12_x = Mat.init_cols n_mo n_mo ( fun i j -> \n", - " phi_x_phi.{i,j} *. phi_x_phi.{i,j} \n", - " -. (1. /. 4.)(phi_x_phi.{i,i} -. phi_x_phi.{j,j} *. phi_x_phi.{i,i} -. phi_x_phi.{j,j}))\n", - " in\n", - " let m_a12_y = Mat.init_cols n_mo n_mo ( fun i j -> \n", - " phi_y_phi.{i,j} *. phi_y_phi.{i,j} \n", - " -. (1. /. 4.)(phi_y_phi.{i,i} -. phi_y_phi.{j,j} *. phi_y_phi.{i,i} -. phi_y_phi.{j,j}))\n", - " in\n", - " let m_a12_z = Mat.init_cols n_mo n_mo ( fun i j -> \n", - " phi_z_phi.{i,j} *. phi_z_phi.{i,j} \n", - " -. (1. /. 4.)(phi_z_phi.{i,i} -. phi_z_phi.{j,j} *. phi_z_phi.{i,i} -. phi_z_phi.{j,j}))\n", - " in \n", - " Mat.add m_a12_x Mat.add (m_a12_y m_a12_z)\n", - "\n", - " let m_b12 = \n", - "\n", - " let m_b12_x = Mat.init_cols n_mo n_mo ( fun i j ->\n", - " (phi_x_phi.{i,i} -. phi_x_phi.{j,j}) *. phi_x_phi.{i,j})\n", - " in \n", - " let m_b12_y = Mat.init_cols n_mo n_mo ( fun i j ->\n", - " (phi_y_phi.{i,i} -. phi_y_phi.{j,j}) *. phi_y_phi.{i,j})\n", - " in \n", - " let m_b12_z = Mat.init_cols n_mo n_mo ( fun i j ->\n", - " (phi_z_phi.{i,i} -. phi_z_phi.{j,j}) *. phi_z_phi.{i,j})\n", - " in \n", - " Mat.add m_b12_x Mat.add (m_b12_y m_b12_z)\n", - "\n", - "in\n", - "Mat.init_cols n_mo n_mo ( fun i j ->\n", - " asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )));;*)\n" + " \n" ] }, { @@ -1358,14 +1325,14 @@ "outputs": [], "source": [ "(* Calcul de toutes les intégrales B_12 -> Matrice *)\n", - "let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> \n", + "let m_b12 = Mat.init_cols n_ao n_ao (fun i j -> \n", "integral_general (fun a b e f i j ->\n", " ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}\n", " ) i j\n", " )\n", "\n", "(* Calcul de toutes les intégrales A_12 -> Matrice *)\n", - "let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->\n", + "let m_a12 = Mat.init_cols n_ao n_ao (fun i j ->\n", "integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,j} \n", " -. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} ) \n", " *. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )\n", @@ -1373,7 +1340,7 @@ " )\n", "\n", "(* Calcul de tous les alpha -> Matrice *)\n", - "let m_alpha = Mat.init_cols n_mo n_mo ( fun i j ->\n", + "let m_alpha = Mat.init_cols n_ao n_ao ( fun i j ->\n", " asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n", " )\n", " )" @@ -1390,13 +1357,104 @@ "|> iamax;;\n", "\n", "(* Indices du plus grand angle alpha *)\n", - "let indice_i = ij_max_alpha mod n_mo;;\n", - "let indice_j = ij_max_alpha / n_mo + 1;;\n", + "let indice_i = ij_max_alpha mod n_ao;;\n", + "let indice_j = ij_max_alpha / n_ao + 1;;\n", "\n", "(* Valeur du plus grand angle alpha*)\n", "let alpha = m_alpha.{indice_i,indice_j};;" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* Test algo alpha < pi/. 2. *)\n", + "let prev_alpha = m_alpha\n", + "\n", + "let rec m_new_alpha m_alpha n = \n", + " \n", + " Printf.printf \"%d\\n%!\" n;\n", + " if n == 0 then\n", + " prev_alpha\n", + " else\n", + " \n", + " \n", + " let m_alpha =\n", + "\n", + " Mat.init_cols n_ao n_ao (fun i j -> if m_alpha.{i,j} > (3.14159265359 /. 2.) \n", + " then (m_alpha.{i,j} -. (3.14159265359 /. 2.))\n", + " else m_alpha.{i,j})\n", + "\n", + " in\n", + " (* Extraction du plus grand angle alpha*)\n", + " let ij_max_alpha = Mat.as_vec m_alpha\n", + " |> iamax\n", + " in\n", + " (* Indices du plus grand angle alpha *)\n", + " let indice_i = ij_max_alpha mod n_ao\n", + " in\n", + " let indice_j = ij_max_alpha / n_ao + 1\n", + " in\n", + " (* Valeur du plus grand angle alpha*)\n", + " let alpha = m_alpha.{indice_i,indice_j}\n", + " in\n", + "\n", + " if alpha < (3.14159265359 /. 2.)\n", + " then\n", + " m_alpha\n", + " else\n", + " let m_alpha = prev_alpha\n", + " in\n", + " m_new_alpha m_alpha (n-1);;" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m_new_alpha m_alpha 2;;\n", + "alpha;;" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* Test méthode de calcul de alpha *)\n", + "let calcul_alpha_2 methode = \n", + "\n", + "let alpha_boys = 1.\n", + "in\n", + "let alpha_er = 2.\n", + "in\n", + "let alpha methode =\n", + " match methode with \n", + " | \"Boys\"\n", + " | \"boys\" -> alpha_boys\n", + " | \"ER\"\n", + " | \"er\" -> alpha_er\n", + " | _ -> invalid_arg \"Unknown method, please enter Boys or ER\"\n", + "\n", + "in \n", + "(alpha methode)**2.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "calcul_alpha_2 \"ER\";;\n" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -1535,26 +1593,26 @@ " if n == 0 then\n", " prev_iter\n", " else\n", - "\n", + "(*\n", " (* Calcul de tous les alpha -> Matrice *)\n", " let m_alpha m_Phi =\n", "\n", " (* Calcul de toutes les intégrales B_12 -> Matrice *)\n", - " let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> \n", + " let m_b12 = Mat.init_cols n_ao n_ao (fun i j -> \n", " integral_general (fun a b e f i j ->\n", " ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} ) *. m_Phi.{e,i} *. m_Phi.{f,j}\n", " ) i j\n", " )\n", " in\n", " (* Calcul de toutes les intégrales A_12 -> Matrice *)\n", - " let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->\n", + " let m_a12 = Mat.init_cols n_ao n_ao (fun i j ->\n", " integral_general (fun a b e f i j -> m_Phi.{a,i} *. m_Phi.{b,i} *. m_Phi.{e,i} *. m_Phi.{f,j} \n", " -. 0.25 *. ( m_Phi.{e,i} *. m_Phi.{f,i} -. m_Phi.{e,j} *. m_Phi.{f,j} ) \n", " *. ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} )\n", " ) i j\n", " )\n", " in\n", - " Mat.init_cols n_mo n_mo ( fun i j ->\n", + " Mat.init_cols n_ao n_ao ( fun i j ->\n", " asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n", " )\n", " )\n", @@ -1572,9 +1630,9 @@ " in\n", "\n", " (* Indices du plus grand angle alpha *)\n", - " let indice_i = ij_max_alpha mod n_mo\n", + " let indice_i = ij_max_alpha mod n_ao\n", " in\n", - " let indice_j = ij_max_alpha / n_mo + 1\n", + " let indice_j = ij_max_alpha / n_ao + 1\n", "\n", " in\n", " m_alpha.{indice_i,indice_j}\n", @@ -1620,19 +1678,20 @@ " else 0.)\n", "\n", " in\n", - " Mat.sub m_C m_Psi\n", + " Mat.sub m_Phi m_Psi\n", " in\n", - " \n", + " *)\n", " let coef = Mat.add m_Psi_tilde m_interm in\n", " let loc = alpha \n", " in\n", " let result = \n", " { coef ; loc }\n", " in\n", - " if loc -. prev_iter.loc < 1.e-5 then\n", - " result\n", + " if loc -. prev_iter.loc < 1.e-5 \n", + " then\n", + " result\n", " else\n", - " let result = prev_iter\n", + " let result = prev_iter\n", " in\n", " new_Phi result (n-1)\n", "\n", @@ -1641,6 +1700,158 @@ "\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* Fonction pour itérer , cc critère de convergence*) \n", + "(* Localisation de Edminstion *)\n", + "let ee_ints = AOBasis.ee_ints ao_basis;;\n", + "let m_C = MOBasis.mo_coef mo_basis;;\n", + "let n_ao = Mat.dim1 m_C ;;\n", + "let n_mo = n_ao - 1;;\n", + "let sum a = \n", + " Array.fold_left (fun accu x -> accu +. x) 0. a\n", + ";;\n", + "\n", + "(*Calcul des intégrales*) \n", + "let integral_general g i j =\n", + "Array.map (fun a ->\n", + " let v = \n", + " Array.map (fun b ->\n", + " let u = \n", + " Array.map (fun e ->\n", + " let t = Array.map (fun f ->\n", + " (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n", + " ) (Util.array_range 1 n_ao)\n", + " in sum t\n", + " ) (Util.array_range 1 n_ao)\n", + " in sum u\n", + " ) (Util.array_range 1 n_ao)\n", + " in sum v\n", + ") (Util.array_range 1 n_ao)\n", + "|> sum \n", + "\n", + "\n", + "let m_Phi = m_C;;\n", + "\n", + "type iteration = \n", + "{ coef : Mat.t ;\n", + " loc : float ;\n", + "}\n", + "\n", + "let prev_iter = {coef = m_Phi ;loc = 0. }\n", + "\n", + "(* Calcul de la nouvelle matrice Phi après rotations *)\n", + "let rec new_Phi prev_iter n = \n", + " Printf.printf \"%d\\n%!\" n;\n", + " if n == 0 then\n", + " prev_iter\n", + " else\n", + "\n", + " (* Calcul de tous les alpha -> Matrice *)\n", + " let m_alpha m_Phi =\n", + "\n", + " (* Calcul de toutes les intégrales B_12 -> Matrice *)\n", + " let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> \n", + " integral_general (fun a b e f i j ->\n", + " ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} ) *. m_Phi.{e,i} *. m_Phi.{f,j}\n", + " ) i j\n", + " )\n", + " in\n", + " (* Calcul de toutes les intégrales A_12 -> Matrice *)\n", + " let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->\n", + " integral_general (fun a b e f i j -> m_Phi.{a,i} *. m_Phi.{b,i} *. m_Phi.{e,i} *. m_Phi.{f,j} \n", + " -. 0.25 *. ( m_Phi.{e,i} *. m_Phi.{f,i} -. m_Phi.{e,j} *. m_Phi.{f,j} ) \n", + " *. ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} )\n", + " ) i j\n", + " )\n", + " in\n", + " Mat.init_cols n_mo n_mo ( fun i j ->\n", + " asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n", + " )\n", + " )\n", + "\n", + " in\n", + " (* Matrice de rotation 2 par 2 *)\n", + " let m_R m_alpha = \n", + "\n", + " (* Valeur du plus grand angle alpha*)\n", + " let alpha = \n", + "\n", + " (* Extraction du plus grand angle alpha*)\n", + " let ij_max_alpha = Mat.as_vec m_alpha\n", + " |> iamax\n", + " in\n", + "\n", + " (* Indices du plus grand angle alpha *)\n", + " let indice_i = ij_max_alpha mod n_mo\n", + " in\n", + " let indice_j = ij_max_alpha / n_mo + 1\n", + "\n", + " in\n", + " m_alpha.{indice_i,indice_j}\n", + " in\n", + " Printf.printf \"%f\\n%!\" alpha;\n", + "\n", + "\n", + " Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha\n", + " else if i>j then sin alpha \n", + " else -. sin alpha )\n", + "\n", + " in\n", + " (* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle \n", + " on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n", + " let m_Ksi_tilde m_R = \n", + "\n", + " (* Fonction d'extraction des 2 vecteurs propres i et j de la matrice des OMs pour les mettres dans la matrice Ksi \n", + " (n par 2) pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n", + " let p = indice_i\n", + " in\n", + " let q = indice_j\n", + " in\n", + " let m_Ksi = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_Phi.{i,p} else m_Phi.{i,q} )\n", + " in\n", + " gemm m_Ksi m_R\n", + "\n", + " in\n", + " (* Construction de la nouvelle matrice d'OMs phi~ *)\n", + " let m_Phi_tilde m_Ksi m_Ksi_tilde = \n", + "\n", + " (* Matrice intérmédiaire où les orbitales i et j ont été supprimées et remplacées par des 0, par soustraction de la \n", + " matrice Phi par la matrice *)\n", + " let m_interm = \n", + "\n", + " (* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n", + " let m_Psi_tilde = Mat.init_cols n_ao n_ao (fun i j -> if j=q then m_Ksi_tilde.{i,2}\n", + " else if j=p then m_Ksi_tilde.{i,1}\n", + " else 0.)\n", + " in\n", + " (* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n", + " let m_Psi = Mat.init_cols n_ao n_ao (fun i j -> if j=q then m_Ksi.{i,2}\n", + " else if j=p then m_Ksi.{i,1}\n", + " else 0.)\n", + "\n", + " in\n", + " Mat.sub m_Phi m_Psi\n", + " in\n", + "\n", + " let coef = Mat.add m_Psi_tilde m_interm in\n", + " let loc = alpha \n", + " in\n", + " let result = \n", + " { coef ; loc }\n", + " in\n", + " if loc -. prev_iter.loc < 1.e-5 then\n", + " result\n", + " else\n", + " let result = prev_iter\n", + " in\n", + " new_Phi result (n-1)" + ] + }, { "cell_type": "code", "execution_count": null,