diff --git a/Work.ipynb b/Work.ipynb index 64efcd0..622cbbf 100644 --- a/Work.ipynb +++ b/Work.ipynb @@ -1593,9 +1593,10 @@ " if n == 0 then\n", " prev_iter\n", " else\n", - "(*\n", + " let m_Phi = prev_iter.coef\n", + " in\n", " (* Calcul de tous les alpha -> Matrice *)\n", - " let m_alpha m_Phi =\n", + " let f_alpha m_Phi =\n", "\n", " (* Calcul de toutes les intégrales B_12 -> Matrice *)\n", " let m_b12 = Mat.init_cols n_ao n_ao (fun i j -> \n", @@ -1612,15 +1613,17 @@ " ) i j\n", " )\n", " in\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", + " 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", "\n", " in\n", " (* Matrice de rotation 2 par 2 *)\n", - " let m_R m_alpha = \n", + " let m_R = \n", "\n", + " let m_alpha = f_alpha m_Phi\n", + " in\n", " (* Valeur du plus grand angle alpha*)\n", " let alpha = \n", "\n", @@ -1634,20 +1637,19 @@ " in\n", " let indice_j = ij_max_alpha / n_ao + 1\n", "\n", - " in\n", - " m_alpha.{indice_i,indice_j}\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", + " 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", + " let m_Ksi_tilde = \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", @@ -1656,42 +1658,43 @@ " 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", + " 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", + " \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 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", - " Mat.sub m_Phi m_Psi\n", + " Mat.add m_Psi_tilde m_interm\n", + " in \n", + " let coef = m_Phi_tilde m_Ksi m_Ksi_tilde\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 \n", + " if abs_float (loc -. prev_iter.loc) < 1.e-5 \n", " then\n", - " result\n", + " { coef ; loc }\n", " else\n", - " let result = prev_iter\n", + " prev_iter\n", " in\n", " new_Phi result (n-1)\n", "\n", @@ -1706,159 +1709,7 @@ "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, - "metadata": {}, - "outputs": [], - "source": [ - " new_Phi prev_iter 10;;" + " new_Phi prev_iter 3;;" ] }, {