diff --git a/Sauvegarde.ipynb b/Sauvegarde.ipynb new file mode 100644 index 0000000..5b77777 --- /dev/null +++ b/Sauvegarde.ipynb @@ -0,0 +1,737 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sauvegarde des anciennes fonctions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(*\n", + "(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n", + "let f_alpha m_C eps=\n", + " let n_mo = Mat.dim2 m_C\n", + " in\n", + " (* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)\n", + " let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n", + " let m_a12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n", + " Array.iter (fun a ->\n", + " Array.iter (fun b ->\n", + " let mca = Vec.init n_mo (fun i -> m_C.{a,i} *. m_C.{b,i}) in\n", + " Array.iter (fun e ->\n", + " Array.iter (fun f ->\n", + " let integral = ERI.get_phys ee_ints f b e a in\n", + " if abs_float integral > eps then\n", + " Array.iter ( fun i -> \n", + " let mcei = m_C.{e,i} *. m_C.{f,i} in\n", + " Array.iter ( fun j -> \n", + " let x = m_C.{e,i} *. m_C.{f,j} *. integral in\n", + " let mcaij = ( mca.{i} -. mca.{j} ) in\n", + " m_b12.{i,j} <- m_b12.{i,j} +. mcaij *. x;\n", + " m_a12.{i,j} <- m_a12.{i,j} +. m_C.{a,i} *. m_C.{b,j} *. x\n", + " -. 0.25 *. ( mcei -. m_C.{e,j} *. m_C.{f,j} ) *. mcaij *. integral \n", + " ) (Util.array_range 1 n_mo)\n", + " ) (Util.array_range 1 n_mo)\n", + " ) (Util.array_range 1 n_ao)\n", + " ) (Util.array_range 1 n_ao)\n", + " ) (Util.array_range 1 n_ao)\n", + " ) (Util.array_range 1 n_ao);\n", + "\n", + " Mat.init_cols n_mo n_mo ( fun i j -> \n", + " if i= j then 0. \n", + " else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))));;\n", + "\n", + "(*********************)\n", + "\n", + "f_alpha m_C 1.e-5;; \n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(*\n", + "(* Calcul de D -> critère à maximiser dans ER*)\n", + "let s_D m_C = \n", + " let v_D = \n", + " let n_mo = Mat.dim2 m_C\n", + " in\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", + " in Vec.init n_mo ( fun i -> m_D.{i,i} )\n", + " in Vec.sum v_D ;;\n", + "\n", + "(******************)\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", + "let toto = s_D m_C;;\n", + "\n", + "*)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Boucle localisation avec fonctions définies à l'extérieur" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* Ancienne boucle\n", + "\n", + "(* Localisation de Edminstion *)\n", + "let n_rec_alpha = 10;;\n", + "(* Calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)\n", + "let rec final_m_C m_C n =\n", + "\n", + " Printf.printf \"%i\\n%!\" n;\n", + "\n", + " (*Util.debug_matrix \"m_C\" m_C;*)\n", + "\n", + "if n == 0 \n", + " then m_C\n", + " else\n", + " \n", + " (* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)\n", + " let new_m_C m_C =\n", + " \n", + " (* D critère à maximiser *)\n", + " let critere_D = s_D m_C (* Fonction -> constante *)\n", + " in\n", + " Printf.printf \"%f\\n%!\" critere_D;\n", + " \n", + " (* Matrice des alphas *)\n", + " let m_alpha = f_alpha m_C (* Fonction -> constante *)\n", + " in\n", + "\n", + " (*Util.debug_matrix \"m_alpha\" m_alpha;*)\n", + "\n", + " (* alphaij contient le alpha max ainsi que ses indices i et j *)\n", + " let alphaij = new_m_alpha m_alpha n_rec_alpha (* Fonction -> constante *)\n", + " in\n", + "\n", + " (* Valeur de alpha max après calcul *)\n", + " let alpha = alphaij.alpha_max (* Fonction -> constante *)\n", + " in\n", + "\n", + " (* Indice i et j du alpha max après calcul *)\n", + " let indice_i = alphaij.indice_ii (* Fonction -> constante *)\n", + " in\n", + " let indice_j = alphaij.indice_jj (* Fonction -> constante *)\n", + " in\n", + "\n", + " (*Printf.printf \"%i %i\\n%!\" indice_i indice_j;*)\n", + "\n", + " (* Matrice de rotaion *)\n", + " let m_R = f_R alpha (* Fonction -> constante *)\n", + " in\n", + "\n", + " (*Util.debug_matrix \"m_R\" m_R;*)\n", + "\n", + " (* Matrice qui va subir la rotation *)\n", + " let m_Ksi = f_Ksi indice_i indice_j m_C (* Fonction -> constante *)\n", + " in\n", + "\n", + " (*Util.debug_matrix \"m_Ksi\" m_Ksi;*)\n", + "\n", + " (* Matrice ayant subit la rotation *)\n", + " let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi (* Fonction -> constante *)\n", + " in\n", + "\n", + " (*Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;*)\n", + "\n", + " (* Matrice pour supprimerles coef des orbitales i et j dans la matrice des coef *)\n", + " let m_Psi = f_Psi m_Ksi indice_i indice_j (* Fonction -> constante *)\n", + " in\n", + "\n", + " (*Util.debug_matrix \"m_Psi\" m_Psi;*)\n", + "\n", + " (* Matrice pour ajouter les coef des orbitales i~ et j~ dans la matrice des coef *)\n", + " let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j (* Fonction -> constante *)\n", + " in\n", + "\n", + " (*Util.debug_matrix \"m_Psi_tilde\" m_Psi_tilde;*)\n", + "\n", + " (* Matrice avec les coef des orbitales i et j remplacés par 0 *)\n", + " let m_interm = f_interm m_C m_Psi (* Fonction -> constante *)\n", + " in\n", + " (* Matrice après rotation *)\n", + " Mat.add m_Psi_tilde m_interm\n", + " in\n", + " let m_new_m_C = new_m_C m_C (* Fonction -> constante *)\n", + "\n", + "in\n", + "\n", + "(*Util.debug_matrix \"new_alpha_m\" (f_alpha m_C);*)\n", + "(*Util.debug_matrix \"m_new_m_C\" m_new_m_C;*)\n", + "\n", + "final_m_C m_new_m_C (n-1);;\n", + "\n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(*\n", + "(*Cellule de calcul*)\n", + "final_m_C m_C 10;;\n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(*Définition des fonctions teste*)\n", + "(*\n", + "(* Localisation de Edminstion *)\n", + "\n", + "(* Définitions de base nécessaire pour la suite *)\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", + "\n", + "(*Fonction général de 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", + "(******************************************************************************************************)\n", + "\n", + "(* Fonction de calcul de tous les alpha -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n", + "let f_alpha m_C =\n", + "\n", + "(* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)\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", + "in\n", + "(* Fonction de calcul de toutes les intégrales A_12 -> Matrice, dépend de m_C *)\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,j} *. 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", + " ) i j\n", + " )\n", + "in\n", + " Mat.init_cols n_ao n_ao ( fun i j ->\n", + " if i= j \n", + " then 0. \n", + " else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))))\n", + "\n", + ";;\n", + "\n", + "\n", + "f_alpha m_C;;\n", + " \n", + "(******************************************************************************************************)\n", + "\n", + "(* Calcul de D *)\n", + " let s_D m_C = \n", + " let v_D = \n", + " let m_D = 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,i} \n", + " ) i j\n", + " )\n", + " in Vec.init n_ao ( fun i -> m_D.{i,1} )\n", + " in Vec.sum v_D ;;\n", + " \n", + " s_D m_C;;\n", + " \n", + "(******************************************************************************************************)\n", + "\n", + "(* Si alpha max > pi/2 on doit soustraire pi/2 à la matrice des alphas *)\n", + "let m_alpha = f_alpha m_C;;\n", + "\n", + "type alphaij = {\n", + "alpha_max : float;\n", + "indice_ii : int;\n", + "indice_jj : int;}\n", + "\n", + "let rec new_m_alpha m_alpha n_rec_alpha=\n", + " let alpha_m =\n", + " Printf.printf \"%i\\n%!\" n_rec_alpha;\n", + " if n_rec_alpha == 0 \n", + " then m_alpha \n", + " else Mat.init_cols n_ao n_ao (fun i j -> \n", + " if (m_alpha.{i,j}) > 3.14 /. 2. \n", + " then (m_alpha.{i,j} -. ( 3.14 /. 2.))\n", + " else if m_alpha.{i,j} < -. 3.14 /. 2.\n", + " then (m_alpha.{i,j} +. ( 3.14 /. 2.))\n", + " else if m_alpha.{i,j} < 0. \n", + " then -. m_alpha.{i,j}\n", + " else m_alpha.{i,j} )\n", + " \n", + " in \n", + " Util.debug_matrix \"alpha_m\" alpha_m;\n", + " let max_element3 alpha_m = \n", + " Mat.as_vec alpha_m\n", + " |> iamax\n", + " in\n", + " let indice_ii = \n", + " let max = max_element3 alpha_m\n", + " in\n", + " Printf.printf \"%i\\n%!\" max;\n", + " (max - 1) mod n_ao +1 \n", + " in\n", + " let indice_jj = \n", + " let max = max_element3 alpha_m\n", + " in\n", + " (max - 1) / n_ao +1\n", + " in\n", + " let alpha alpha_m = \n", + " let i = indice_ii \n", + " in\n", + " let j = indice_jj \n", + " in\n", + " alpha_m.{i,j}\n", + " in\n", + " let alpha_max = alpha alpha_m \n", + " in\n", + " Printf.printf \"%f\\n%!\" alpha_max;\n", + " if alpha_max < 3.14 /. 2.\n", + " then {alpha_max; indice_ii; indice_jj}\n", + " else new_m_alpha alpha_m (n_rec_alpha-1);;\n", + " \n", + "\n", + "new_m_alpha m_alpha 10 ;;\n", + "\n", + "let alphaij = new_m_alpha m_alpha 10 ;;\n", + "\n", + "(******************************************************************************************)\n", + "let alpha = alphaij.alpha_max\n", + "\n", + "(* Matrice de rotation 2 par 2 *)\n", + "let m_R 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", + "\n", + "(******************************************************************************************)\n", + "\n", + "let indice_i = alphaij.indice_ii;;\n", + "let indice_j = alphaij.indice_jj;;\n", + "let m_R = m_R alpha;;\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 par 2)\n", + "pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n", + "let m_Ksi indice_i indice_j m_C = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_C.{i,indice_i} else m_C.{i,indice_j} );;\n", + "\n", + "m_Ksi indice_i indice_j m_C;;\n", + "\n", + "let m_Ksi = m_Ksi indice_i indice_j m_C;;\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_tilde m_R m_Ksi = gemm m_Ksi m_R;;\n", + "\n", + "m_Ksi_tilde m_R m_Ksi ;;\n", + "\n", + "(******************************************************************************************)\n", + "let m_Ksi_tilde = m_Ksi_tilde m_R m_Ksi\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", + "celle contenant i~ et j~ *)\n", + "\n", + "(* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n", + "let m_Psi_tilde m_Ksi_tilde indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi_tilde.{i,1}\n", + " else if j=indice_j then m_Ksi_tilde.{i,2}\n", + " else 0.);;\n", + "\n", + "m_Psi_tilde m_Ksi_tilde indice_i indice_j;;\n", + "\n", + "(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n", + "let m_Psi m_Ksi indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi.{i,1}\n", + " else if j=indice_j then m_Ksi.{i,2}\n", + " else 0.);;\n", + " \n", + "m_Psi m_Ksi indice_i indice_j;;\n", + "\n", + "(******************************************************************************************)\n", + "let m_Psi = m_Psi m_Ksi indice_i indice_j;;\n", + "let m_Psi_tilde = m_Psi_tilde m_Ksi_tilde indice_i indice_j;; \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 matrice Phi\n", + "par la matrice *)\n", + "let m_interm m_C m_Psi = Mat.sub m_C m_Psi;;\n", + "\n", + "let m_interm = m_interm m_C m_Psi;;\n", + "\n", + "(* Construction de la nouvelle matrice d'OMs phi~ *)\n", + "let m_Phi_tilde m_Psi_tilde m_interm = Mat.add m_Psi_tilde m_interm;;\n", + "\n", + "m_Phi_tilde m_Psi_tilde m_interm;;\n", + "\n", + "*)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Programmes avec fonctions définies à l'intérieurs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* Localisation de Edminstion *)\n", + "\n", + "(*Fonctionne -> programme avec les fonctions définies à l'intérieur*)\n", + "\n", + "(*\n", + "\n", + "(* Définitions de base nécessaire pour la suite *)\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", + "\n", + "(*Fonction général de 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", + "\n", + "type alphaij = {\n", + " alpha_max : float;\n", + " indice_ii : int;\n", + " indice_jj : int;};;\n", + " \n", + "let n_rec_alpha = 10;;\n", + "\n", + "(*\n", + "let alpha = 1.;;\n", + "let y_R 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", + "let y_R = y_R alpha;; \n", + "let m_C = gemm m_C y_R;;\n", + "*)\n", + "(******************************************************************************************************)\n", + "\n", + "let rec final_m_C m_C n =\n", + "\n", + "Printf.printf \"%i\\n%!\" n;\n", + "(*\n", + "Util.debug_matrix \"m_C\" m_C;\n", + "*)\n", + "\n", + "if n == 0 then m_C\n", + " else\n", + "\n", + "let new_m_C m_C =\n", + "\n", + " (* Fonction de calcul de tous les alpha -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n", + " let f_alpha m_C =\n", + "\n", + " (* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)\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", + " in\n", + " (* Fonction de calcul de toutes les intégrales A_12 -> Matrice, dépend de m_C *)\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,j} *. 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", + " ) i j\n", + " )\n", + " in\n", + " Mat.init_cols n_ao n_ao ( fun i j -> if i= j \n", + " then 0. \n", + " else\n", + " 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))))\n", + " in\n", + " \n", + " (* Calcul de D *)\n", + " let s_D m_C = \n", + " let v_D = \n", + " let m_D = 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,i} \n", + " ) i j\n", + " )\n", + " in Vec.init n_ao ( fun i -> m_D.{i,1} )\n", + " in Vec.sum v_D \n", + " in\n", + " let critere_D = s_D m_C (* Fonction -> constante *)\n", + " in\n", + " Printf.printf \"%f\\n%!\" critere_D;\n", + " \n", + " let m_alpha = f_alpha m_C (* Fonction -> constante *)\n", + " in\n", + " (*\n", + " Util.debug_matrix \"m_alpha\" m_alpha;\n", + " *)\n", + " \n", + " (* Détermination alpha_max et ses indices i et j.\n", + " Si alpha max > pi/2 on soustrait pi/2 à la matrice des alphas de manière récursive *)\n", + " let rec new_m_alpha m_alpha n_rec_alpha=\n", + " let alpha_m =\n", + " \n", + " Printf.printf \"%i\\n%!\" n_rec_alpha;\n", + " \n", + " if n_rec_alpha == 0 \n", + " then m_alpha \n", + " else Mat.init_cols n_ao n_ao (fun i j -> \n", + " if (m_alpha.{i,j}) > (3.14 /. 2.) \n", + " then (m_alpha.{i,j} -. ( 3.14 /. 2.))\n", + " else if m_alpha.{i,j} < -. 3.14 /. 2.\n", + " then (m_alpha.{i,j} +. ( 3.14 /. 2.))\n", + " else if m_alpha.{i,j} < 0. \n", + " then -. m_alpha.{i,j}\n", + " else m_alpha.{i,j} )\n", + " in \n", + " \n", + " Util.debug_matrix \"alpha_m\" alpha_m;\n", + " \n", + " (* Détermination de l'emplacement du alpha max *)\n", + " let max_element3 alpha_m = \n", + " Mat.as_vec alpha_m\n", + " |> iamax\n", + " in\n", + " \n", + " (* indice i du alpha max *)\n", + " let indice_ii = \n", + " let max = max_element3 alpha_m (* Fonction -> constante *)\n", + " in\n", + " (*\n", + " Printf.printf \"%i\\n%!\" max;\n", + " *)\n", + " (max - 1) mod n_ao +1 \n", + " in\n", + " \n", + " (* indice j du alpha max *)\n", + " let indice_jj = \n", + " let max = max_element3 alpha_m (* Fonction -> constante *)\n", + " in\n", + " (max - 1) / n_ao +1\n", + " in\n", + " \n", + " (* Valeur du alpha max*)\n", + " let alpha alpha_m = \n", + " let i = indice_ii \n", + " in\n", + " let j = indice_jj \n", + " in\n", + " (*\n", + " Printf.printf \"%i %i\\n%!\" i j;\n", + " *)\n", + " alpha_m.{i,j}\n", + " \n", + " in\n", + " let alpha_max = alpha alpha_m (* Fonction -> constante *)\n", + " in\n", + " Printf.printf \"%f\\n%!\" alpha_max;\n", + " if alpha_max < 3.14 /. 2.\n", + " then {alpha_max; indice_ii; indice_jj}\n", + " else new_m_alpha alpha_m (n_rec_alpha-1)\n", + " in\n", + " let alphaij = new_m_alpha m_alpha n_rec_alpha (* Fonction -> constante *)\n", + " in\n", + " \n", + " (* Valeur de alpha max après calcul *)\n", + " let alpha = alphaij.alpha_max (* Fonction -> constante *)\n", + " in\n", + " \n", + " (* Matrice de rotation 2 par 2 *)\n", + " let f_R 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", + " in\n", + " (* Indice i et j du alpha max après calcul *)\n", + " let indice_i = alphaij.indice_ii (* Fonction -> constante *)\n", + " in\n", + " let indice_j = alphaij.indice_jj (* Fonction -> constante *)\n", + " in\n", + " \n", + " Printf.printf \"%i %i\\n%!\" indice_i indice_j;\n", + " \n", + " let m_R = f_R alpha (* Fonction -> constante *)\n", + " in\n", + " Util.debug_matrix \"m_R\" m_R;\n", + " (* Fonction d'extraction des 2 vecteurs propres i et j de la matrice des OMs pour les mettres dans la matrice Ksi (n par 2)\n", + " pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n", + " let f_Ksi indice_i indice_j m_C = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_C.{i,indice_i} else m_C.{i,indice_j} )\n", + " in\n", + " let m_Ksi = f_Ksi indice_i indice_j m_C (* Fonction -> constante *)\n", + " in\n", + " \n", + " Util.debug_matrix \"m_Ksi\" m_Ksi;\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 f_Ksi_tilde m_R m_Ksi = gemm m_Ksi m_R\n", + " in\n", + " let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi (* Fonction -> constante *)\n", + " in\n", + " \n", + " Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\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", + " celle contenant i~ et j~ *)\n", + " \n", + " (* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n", + " let f_Psi_tilde m_Ksi_tilde indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi_tilde.{i,1}\n", + " else if j=indice_j then m_Ksi_tilde.{i,2}\n", + " else 0.)\n", + "\n", + " in\n", + "\n", + " (* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n", + " let f_Psi m_Ksi indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi.{i,1}\n", + " else if j=indice_j then m_Ksi.{i,2}\n", + " else 0.)\n", + " in\n", + " let m_Psi = f_Psi m_Ksi indice_i indice_j (* Fonction -> constante *)\n", + " in\n", + " \n", + " Util.debug_matrix \"m_Psi\" m_Psi;\n", + " \n", + " let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j (* Fonction -> constante *)\n", + " in\n", + " \n", + " Util.debug_matrix \"m_Psi_tilde\" m_Psi_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 matrice Phi\n", + " par la matrice *)\n", + " let f_interm m_C m_Psi = Mat.sub m_C m_Psi\n", + " in\n", + " let m_interm = f_interm m_C m_Psi (* Fonction -> constante *)\n", + "in\n", + "Mat.add m_Psi_tilde m_interm\n", + "in\n", + "let m_new_m_C = new_m_C m_C (* Fonction -> constante *)\n", + "in\n", + "Util.debug_matrix \"new_alpha_m\" (f_alpha m_C);\n", + "\n", + "Util.debug_matrix \"m_new_m_C\" m_new_m_C;\n", + "\n", + "final_m_C m_new_m_C (n-1);;\n", + "\n", + "(*****************************)\n", + "\n", + "final_m_C m_C 10;;\n", + "\n", + "*)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "OCaml 4.10.0", + "language": "OCaml", + "name": "ocaml-jupyter" + }, + "language_info": { + "codemirror_mode": "text/x-ocaml", + "file_extension": ".ml", + "mimetype": "text/x-ocaml", + "name": "OCaml", + "nbconverter_exporter": null, + "pygments_lexer": "OCaml", + "version": "4.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/Test.ipynb b/Test.ipynb new file mode 100644 index 0000000..c0b7c2e --- /dev/null +++ b/Test.ipynb @@ -0,0 +1,193 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#use \"topfind\";;\n", + "#require \"jupyter.notebook\";;\n", + "#require \"lacaml.top\";;\n", + "#require \"alcotest\";;\n", + "#require \"str\";;\n", + "#require \"bigarray\";;\n", + "#require \"zarith\";;\n", + "#require \"getopt\";;\n", + "open Lacaml.D" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "let n_ao = 8;;\n", + "let n_mo = 6;; \n", + "let ran=\n", + "Mat.random ~range:1. ~from:0. n_ao n_mo;; " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* Extraction des colonnes i,j,k,... d'une matrice à partir d'une liste [i,j,k,...] *)\n", + "\n", + "let toto = [2;3;6];;\n", + "(*\n", + "ran;;\n", + "let new_ran = Mat.transpose_copy ran ;;\n", + "\n", + "let vec = Mat.to_col_vecs ran;;\n", + "\n", + "vec.(0);;\n", + "\n", + "let f a = vec.(a);;\n", + "\n", + "let n = List.map f toto;;\n", + "\n", + "let vec_list = [vec.(0);vec.(1)];;\n", + "\n", + "\n", + "let mat_n = Mat.of_col_vecs_list vec_list;;\n", + "\n", + "\n", + "let occ ran toto = \n", + " let new_ran = Mat.transpose_copy ran\n", + " in\n", + " let vec = Mat.to_col_vecs new_ran\n", + " in \n", + " let f a = vec.(a+1)\n", + " in \n", + " let vec_list = List.map f toto \n", + "in Mat.of_col_vecs_list vec_list;;\n", + "\n", + "occ ran toto;;\n", + "\n", + "\n", + "\n", + "let list_a = [0;0;2;0;5];;\n", + "\n", + "let new_list list= List.filter (fun x -> x > 0.) list;;\n", + "\n", + "let vecto = Vec.init 8 (fun i -> if List.mem i toto\n", + "then 0.\n", + "else float_of_int(i) );;\n", + "\n", + "let vecto_list = Vec.to_list vecto;;\n", + "\n", + "\n", + "new_list vecto_list;;\n", + "*)\n", + "(******************)\n", + "(*\n", + "ran;;\n", + "let vec_of_mat = Mat.to_col_vecs ran;;\n", + "let f a = vec_of_mat.(a-1);;\n", + "vec_of_mat.(5);;\n", + "let vec_list = List.map f toto;;\n", + "let new_m = Mat.of_col_vecs_list vec_list;;\n", + "*)\n", + "\n", + "(* Fp *)\n", + "let miss_elem mat list = \n", + " let n_mo = Mat.dim2 mat\n", + " in\n", + " let vec = Vec.init (n_mo) (fun i ->\n", + " if List.mem i list\n", + " then 0.\n", + " else float_of_int(i))\n", + " in\n", + " let vec_list = Vec.to_list vec\n", + " in\n", + " let g a = int_of_float(a)\n", + " in\n", + " let vec_list_int = List.map g vec_list\n", + " \n", + " in\n", + " List.filter (fun x -> x > 0) vec_list_int;;\n", + "\n", + "let split_mat mat list =\n", + " let vec_of_mat = Mat.to_col_vecs mat\n", + " in\n", + " let f a = vec_of_mat.(a-1)\n", + " in\n", + " let vec_list_1 = List.map f list\n", + " in\n", + " let list_2 = miss_elem mat list\n", + " in\n", + " let vec_list_2 = List.map f list_2\n", + "in (Mat.of_col_vecs_list vec_list_1,Mat.of_col_vecs_list vec_list_2);; \n", + "\n", + "let m_occ , m_vir = split_mat ran toto;;\n", + "\n", + "(*\n", + "let miss_elem mat list = \n", + " let n_mo = Mat.dim2 mat\n", + " in\n", + " let vec = Vec.init (n_mo) (fun i ->\n", + " if List.mem i list\n", + " then 0.\n", + " else float_of_int(i))\n", + " in\n", + " let vec_list = Vec.to_list vec\n", + " in\n", + " let g a = int_of_float(a)\n", + " in\n", + " let vec_list_int = List.map g vec_list\n", + " \n", + " in\n", + " List.filter (fun x -> x > 0) vec_list_int;;\n", + "*) \n", + "\n", + " \n", + "\n", + "(*let tata = miss_elem toto;;\n", + "\n", + "let titi = [1.;2.]\n", + "\n", + "let f a = int_of_float(a);;\n", + "let test list = List.map f list;;\n", + "\n", + "test titi;;\n", + "*)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "OCaml 4.10.0", + "language": "OCaml", + "name": "ocaml-jupyter" + }, + "language_info": { + "codemirror_mode": "text/x-ocaml", + "file_extension": ".ml", + "mimetype": "text/x-ocaml", + "name": "OCaml", + "nbconverter_exporter": null, + "pygments_lexer": "OCaml", + "version": "4.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/Work.ipynb b/Work.ipynb index 1e28875..cfe15ce 100644 --- a/Work.ipynb +++ b/Work.ipynb @@ -1809,7 +1809,7 @@ " ) i j\n", " );;\n", "let toto = s_D m_C;;\n", - "toto *. 6.;;\n", + "\n", "*)" ] },