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

Work : correction boucle et nettoyage

Browse Source
This commit is contained in:
Yann Damour 2020-05-06 15:39:23 +02:00
parent 0b36d1ea8c
commit 14fc81a642
1 changed files with 93 additions and 310 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

403
Work.ipynb
View File

@ -1199,7 +1199,7 @@
"\n",
"$j~ (x) = -sin(\\gamma) i(x) + cos(\\gamma) j(x)$\n",
"\n",
"On part de la matrices de orbitales moléculaires Phi et on cherche la paire d'orbitale (i,j) ayant le plus grand angle de rotation alpha, avec alpha défini comme :\n",
"On part de la matrices de orbitales moléculaires Phi et on cherche la paire d'orbitale (i,j) ayant le plus grand angle de rotation alpha, avec alpha défini comme (si ce dernier est supétieur à $\\frac{\\pi}{2}$ on devra soustraire $\\frac{\\pi}{2}$ aux éléments de la matrices :\n",
"\n",
"\n",
"$Cos(4 \\alpha)= -A_{12} / (A_{12}^2 + B_{12}^2)^{1/2}$\n",
@ -1246,24 +1246,6 @@
"\n",
"$= \\left(\\sum_e c_{ei} \\sum_f c_{fi} - \\sum_e c_{ej} \\sum_f c_{fj} \\right) \\left(\\sum_a c_{ai} \\sum_b c_{bi} - \\sum_a c_{aj} \\sum_b c_{bj} \\right) [ \\chi_a \\chi_b| \\chi_e \\chi_f ] $\n",
"\n",
"(\n",
"Ainsi comme :\n",
"\n",
"$\\phi_i ^2 = \\left( \\sum_a c_{ai} \\chi_a \\right)^2$\n",
"$= \\sum_a c_{ai} \\sum_b c_{bi} \\chi_a \\chi_b$\n",
"\n",
"On aura :\n",
"\n",
"$D=\\sum_n [\\phi_n^2|\\phi_n^2]$\n",
"\n",
"$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} [\\chi_a \\chi_b| \\phi_n^2])$\n",
"\n",
"$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} \\sum_e c_{en} \\sum_f c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$\n",
"\n",
"$= \\sum_n (\\sum_a \\sum_b \\sum_e \\sum_f c_{an} c_{bn} c_{en} c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$\n",
"\n",
")\n",
"\n",
"Mais aussi :\n",
"\n",
"$B_{ij} = [\\phi_i ^2 - \\phi_j ^2 | \\phi_i \\phi_j ] = [\\left( \\sum_a c_{ai} \\chi_i \\right)^2 - \\left( \\sum_a c_{aj} \\chi_j \\right)^2| \\phi_i \\phi_j ] $\n",
@ -1295,7 +1277,22 @@
"\n",
"Ainsi, on soustrait à la matrice des OMs Phi la matrice m_ksi pour supprimer les OMs i et j de celle ci puis on additionne cette nouvelle matrice à la matrice m_ksi_tilde pour créer la nouvelle matrice des OMs m_Phi_tilde.\n",
"\n",
"On repart de cette nouvelle matrice m_Phi_tilde et on cherche la paire d'orbitale (k,l) ayant le plus grand angle de rotation alpha. Et on procède comme nous l'avons précedemment de manière intérative. Le but étant de maximiser D, c'est à dire maximiser le terme de répulsion.\n"
"On repart de cette nouvelle matrice m_Phi_tilde et on cherche la paire d'orbitale (k,l) ayant le plus grand angle de rotation alpha. Et on procède comme nous l'avons précedemment de manière intérative. Le but étant de maximiser D, c'est à dire maximiser le terme de répulsion.\n",
"\n",
"Pour le calcul de D (le critère de localisation qu'on maximise) , sachant que :\n",
"\n",
"$\\phi_i ^2 = \\left( \\sum_a c_{ai} \\chi_a \\right)^2$\n",
"$= \\sum_a c_{ai} \\sum_b c_{bi} \\chi_a \\chi_b$\n",
"\n",
"On aura :\n",
"\n",
"$D=\\sum_n [\\phi_n^2|\\phi_n^2]$\n",
"\n",
"$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} [\\chi_a \\chi_b| \\phi_n^2])$\n",
"\n",
"$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} \\sum_e c_{en} \\sum_f c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$\n",
"\n",
"$= \\sum_n (\\sum_a \\sum_b \\sum_e \\sum_f c_{an} c_{bn} c_{en} c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$"
]
},
{
@ -1311,10 +1308,12 @@
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"(* Test algo alpha < pi/. 2. *)\n",
"(*)\n",
"let m_alpha = f_alpha m_Phi\n",
"in\n",
"let prev_alpha = m_alpha\n",
"\n",
"in\n",
"let rec m_new_alpha m_alpha n = \n",
" \n",
" Printf.printf \"%d\\n%!\" n;\n",
@ -1353,18 +1352,6 @@
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"m_new_alpha m_alpha 2;;\n",
"alpha;; \n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
@ -1388,17 +1375,9 @@
" | _ -> invalid_arg \"Unknown method, please enter Boys or ER\"\n",
"\n",
"in \n",
"(alpha methode)**2.\n",
"(alpha methode)**2.;;\n",
"\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"*)\n",
"(*\n",
"calcul_alpha_2 \"ER\";;\n",
"*)"
@ -1408,7 +1387,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### ER-V1"
"### ER-V1 \n",
"-> Fonctionne, il n'itère pas => Point de départ"
]
},
{
@ -1417,6 +1397,7 @@
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"\n",
"(* Localisation de Edminstion *)\n",
"let ee_ints = AOBasis.ee_ints ao_basis;;\n",
@ -1467,6 +1448,7 @@
" )\n",
" )\n",
"\n",
"(*\n",
"let s_D = \n",
" let v_D = \n",
" let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
@ -1475,42 +1457,44 @@
" )\n",
" in Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
"in Vec.sum v_D \n",
"*)\n",
"\n",
"(*\n",
"(* Calcul de D *)\n",
"let m_D = Mat.init_cols n_mo n_mo (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,i} \n",
" ) i j\n",
" )\n",
" \n",
" \n",
"\n",
"let v_D = Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
"\n",
"let s_D = Vec.sum v_D\n",
"*)\n",
"\n",
" (* \n",
"\n",
"(* \n",
"(* Calcul de D *)\n",
"let v_D = Vec.init n_ao (fun i ->\n",
"integral_general (fun a b e f i -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
" ) i \n",
"integral_general (fun a b e f i _ -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
" )\n",
" )\n",
"\n",
"let s_D = Vec.sum v_D\n",
"*)\n",
"\n",
"(* Rotation des orbitales *)\n",
"(*\n",
"\n",
"(* Extraction du plus grand angle alpha*)\n",
"let ij_max_alpha = Mat.as_vec m_alpha\n",
"|> iamax;;\n",
"\n",
"(* Indices du plus grand angle alpha *)\n",
"let indice_i = ij_max_alpha mod n_ao;;\n",
"let indice_j = ij_max_alpha / n_ao + 1;;\n",
"let indice_i = (ij_max_alpha - 1) mod n_ao + 1;;\n",
"let indice_j = (ij_max_alpha - 1) / n_mo + 1;;\n",
"\n",
"(* Valeur du plus grand angle alpha*)\n",
"let alpha = m_alpha.{indice_i,indice_j};;\n",
"*)\n",
"\n",
"(*\n",
"\n",
"\n",
"let m_Phi = m_C\n",
"\n",
"(* Matrice de rotation 2 par 2 *)\n",
@ -1527,11 +1511,11 @@
"(* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle on\n",
"obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n",
"let m_Ksi_tilde = gemm m_Ksi m_R;;\n",
"*)\n",
"\n",
"\n",
"(* Réinjection*)\n",
"\n",
"(*\n",
"\n",
"(* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur \n",
"les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice\n",
"des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec \n",
@ -1539,7 +1523,7 @@
"\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,2}\n",
" else if j=p then m_Ksi_tilde.{i,1}\n",
" else 0.);;\n",
" (* p q verif*) \n",
"(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
@ -1553,6 +1537,7 @@
"\n",
"(* Construction de la nouvelle matrice d'OMs phi~ *)\n",
"let m_Phi_tilde = Mat.add m_Psi_tilde m_interm;; \n",
"\n",
"*)"
]
},
@ -1567,7 +1552,8 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Réécriture pour itérerer ER"
"## Réécriture pour itérerer ER\n",
"-> Il fonctionne apparemment mais je suis sceptique pour certains résultats, il faut réécrire toutes les fonctions en fonction de leur arguments pour le refaire mieux et éviter les problèmes rencontrés sur les variables loc"
]
},
{
@ -1576,7 +1562,6 @@
"metadata": {},
"outputs": [],
"source": [
"\n",
"(* Fonction pour itérer , cc critère de convergence*) \n",
"(* Localisation de Edminstion *)\n",
"let ee_ints = AOBasis.ee_ints ao_basis;;\n",
@ -1609,24 +1594,24 @@
"type iteration = \n",
"{ coef : Mat.t ;\n",
" loc : float ;\n",
"}\n",
" }\n",
"\n",
"(*let prev_iter = {coef = m_C ;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",
" (*\n",
" \n",
" Util.debug_matrix \"prev_iter\" prev_iter.coef;\n",
" *)\n",
" \n",
" if n == 0 then\n",
" prev_iter\n",
" else\n",
" let m_Phi = prev_iter.coef\n",
" in\n",
" (*\n",
" \n",
" Util.debug_matrix \"m_coef\" prev_iter.coef;\n",
" Util.debug_matrix \"m_Phi\" m_Phi;\n",
" (*Util.debug_matrix \"m_Phi\" m_Phi;\n",
" *)\n",
" (* Calcul de tous les alpha -> Matrice *)\n",
" let f_alpha m_Phi =\n",
@ -1668,32 +1653,34 @@
" \n",
" let m_alpha = f_alpha m_Phi\n",
" in\n",
" (*\n",
" \n",
" Util.debug_matrix \"m_alpha\" m_alpha;\n",
" *)\n",
" \n",
" (* Valeur du plus grand angle alpha*)\n",
" (* Extraction du plus grand angle alpha*)\n",
" let ij_max_alpha = Mat.as_vec m_alpha\n",
" |> iamax\n",
" in\n",
"\n",
" \n",
" Printf.printf \"%i\\n%!\" ij_max_alpha;\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 - 1) mod n_mo + 1\n",
" in\n",
" let indice_j = ij_max_alpha / n_mo + 1\n",
" let indice_j = (ij_max_alpha - 1) / n_mo + 1\n",
"\n",
" in\n",
" (*\n",
" \n",
" Printf.printf \"%i\\n%!\" indice_i;\n",
" Printf.printf \"%i\\n%!\" indice_j;\n",
" *)\n",
" \n",
" let alpha = \n",
" m_alpha.{indice_i,indice_j}\n",
" in\n",
" \n",
" (*\n",
" \n",
" Printf.printf \"%f\\n%!\" alpha;\n",
" *)\n",
" \n",
" indice_i, indice_j, \n",
" Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha\n",
" else if i>j then sin alpha \n",
@ -1701,55 +1688,59 @@
"\n",
" in\n",
" (*\n",
" Util.debug_matrix \"m_R\" m_R;\n",
" Util.debug_matrix \"m_R\" m_R;*)\n",
" Printf.printf \"%d %d\\n%!\" indice_i indice_j;\n",
" *)\n",
" \n",
" let p = indice_i\n",
" in\n",
" let q = indice_j\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, 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",
" let p = indice_i\n",
" in\n",
" let q = indice_j\n",
" in\n",
" Printf.printf \"%d %d\\n%!\" indice_i indice_j;\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",
" m_Ksi, gemm m_Ksi m_R\n",
"\n",
" in\n",
" (*\n",
" Util.debug_matrix \"m_Ksi\" m_Ksi;\n",
" Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n",
" *)\n",
" in\n",
" Printf.printf \"%d %d\\n%!\" p q;\n",
" Util.debug_matrix \"m_Ksi\" m_Ksi;\n",
" Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n",
" \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_mo (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 intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
" let m_Psi = Mat.init_cols n_ao n_mo (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",
" (* 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_mo (fun i j -> if j=p then m_Ksi_tilde.{i,1}\n",
" else if j=q then m_Ksi_tilde.{i,2}\n",
" else 0.)\n",
" in\n",
" Printf.printf \"%d %d\\n%!\" p q;\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 intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
" let m_Psi = Mat.init_cols n_ao n_mo (fun i j -> if j=p then m_Ksi.{i,1}\n",
" else if j=q then m_Ksi.{i,2}\n",
" else 0.)\n",
"\n",
" in\n",
" (*Util.debug_matrix \"m_Psi\" m_Psi;*)\n",
" Mat.sub m_Phi m_Psi\n",
" in\n",
" (*Util.debug_matrix \"m_interm\" m_interm;*)\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",
" Util.debug_matrix \"coef\" coef;\n",
" Util.debug_matrix \"Nouvelle_matrice_des_coef\" coef;\n",
" \n",
" let loc = s_D \n",
" in\n",
@ -1772,211 +1763,10 @@
"metadata": {},
"outputs": [],
"source": [
"new_Phi {coef = m_C ;loc = 0. } 10;;\n",
"new_Phi {coef = m_C ;loc = 0. } 1;;\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"(* 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 = Mat.dim2 m_C ;;\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",
"let m_Phi = None;;\n",
"type iteration = \n",
"{ coef : Mat.t ;\n",
" loc : float ;\n",
"}\n",
"\n",
"(*let prev_iter = {coef = m_C ;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",
" (*\n",
" Util.debug_matrix \"prev_iter\" prev_iter.coef;\n",
" *)\n",
" if n == 0 then\n",
" prev_iter\n",
" else\n",
" let m_Phi = prev_iter.coef\n",
" in\n",
" (*\n",
" Util.debug_matrix \"m_coef\" prev_iter.coef;\n",
" Util.debug_matrix \"m_Phi\" m_Phi;\n",
" *)\n",
" (* Calcul de tous les alpha -> Matrice *)\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_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",
" in \n",
" (* Matrice de rotation 2 par 2 *)\n",
" let indice_i, indice_j, m_R = \n",
" let n = 3\n",
" in\n",
" let rec new_alpha m_alpha n =\n",
" \n",
" let m_alpha = f_alpha m_Phi\n",
" in\n",
" (*\n",
" Util.debug_matrix \"m_alpha\" m_alpha;\n",
" *)\n",
" (* Valeur du plus grand angle alpha*)\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",
" \n",
" (*\n",
" Printf.printf \"%i\\n%!\" indice_i;\n",
" Printf.printf \"%i\\n%!\" indice_j;\n",
" *)\n",
" in\n",
" if m_alpha.{indice_i,indice_j} < (3.14 /. 2.)\n",
" then m_alpha\n",
" else new_alpha (Mat.sub m_alpha (Mat.make n_ao n_ao (3.14 /. 2.)))\n",
" \n",
" in\n",
" let alpha = m_alpha.{indice_i, indice_j}\n",
" \n",
" \n",
" (*\n",
" Printf.printf \"%f\\n%!\" alpha;\n",
" *)\n",
" indice_i, indice_j, \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",
" (*\n",
" Util.debug_matrix \"m_R\" m_R;\n",
" Printf.printf \"%d %d\\n%!\" indice_i indice_j;\n",
" *)\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, 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",
" 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",
" m_Ksi, gemm m_Ksi m_R\n",
"\n",
" in\n",
" (*\n",
" Util.debug_matrix \"m_Ksi\" m_Ksi;\n",
" Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n",
" *)\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_mo (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 intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
" let m_Psi = Mat.init_cols n_ao n_mo (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.add m_Psi_tilde m_interm\n",
" in \n",
" let coef = m_Phi_tilde m_Ksi m_Ksi_tilde\n",
" in\n",
" \n",
" Util.debug_matrix \"coef\" coef;\n",
" \n",
" let loc = alpha \n",
" in\n",
" (*\n",
" Printf.printf \"%f\\n%!\" loc;\n",
" Printf.printf \"%f\\n%!\" prev_iter.loc;\n",
" *)\n",
" if abs_float (loc -. prev_iter.loc) < 1.e-10 \n",
" then\n",
" { coef ; loc }\n",
" else\n",
" new_Phi {coef; loc;} (n-1)\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
@ -2000,13 +1790,6 @@
"norme m_alpha;;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,