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

Suppression notebook anciens/ mise à jour/ ajout nouveaux

This commit is contained in:
Yann Damour 2020-06-09 11:37:40 +02:00
parent ed290c7f44
commit 6f1e3cea5d
7 changed files with 1369 additions and 4555 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

754
Calcul_test.ipynb
View File

@ -1,754 +0,0 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* Algorithme de localisation des orbitales : Méthode de Boys ou Edminston-Ruedenberg *)\n",
"\n",
"(*I : Position et base atomique*)\n",
" (* 1. Fonction création chaîne linéaire de n H *)\n",
" (* 2. Base atomique *)\n",
"\n",
"(* II : Hartree-Fock *)\n",
" (* 1. Def pour HF *)\n",
" (* 2. Calcul de Hartree-Fock*)\n",
" \n",
"(* III : Définition des fonctions pour la localisation *)\n",
" (* 1. Définitions de base nécessaire pour la suite *) \n",
" (* 2. Fonction de calcul de tous les alpha ER *)\n",
" (* 3. Fonction de calcul de tous les alpha Boys *)\n",
" (* 4. Pattern matching : calcul de alpha et de D *)\n",
" (* 5. Norme de alpha *)\n",
" (* 6. Détermination alpha_max et ses indices i et j *) \n",
" (* 7. Matrice de rotation 2 par 2 *)\n",
" (* 8. Fonction : i,j -> matrice Ksi (n par 2) *)\n",
" (* 9. Fonction : i~,j~ -> matrice Ksi~ (n par 2) *)\n",
" (* 10. Matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)\n",
" (* 11. Matrice intérmédiaire où les orbitales i et j ont été supprimées et remplacées par des 0*)\n",
" \n",
"(* IV : Boucle de localisation des OMs *) \n",
"\n",
"(* V : Fonctions pour la séparation et le rassemblement des matrices *)\n",
" (* 1. Fonction de création d'une list d'entier à partir d'un vecteur de float*)\n",
" (* 2. Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)\n",
" (* 3. Fonction de séparation d'une matrice en 2 sous matrice *)\n",
" (* 4. Liste des OMs occupées *)\n",
" (* 5. Fonction de rassemblement de 2 matrices *) \n",
" (* 6. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)\n",
" \n",
"(* VI : Calcul et assemblage des OMs occupées/virtuels *) \n",
" \n",
"(* VII : Calcul Coef CI *)\n",
" (* 1. Def de la mo_basis pour le HF pre CI *)\n",
" (* 2. HF pour CI *)\n",
" (* 3. Calcul de coef CI *)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"open Lacaml.D\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*I : Position et base atomique*)\n",
"\n",
"(* 1. Fonction création chaîne linéaire de n H *)\n",
"\n",
"let xyz d n = \n",
" let accu = \"\"\n",
" in\n",
" let rec toto accu d n =\n",
" let accu = \n",
" if n=0 \n",
" then accu ^ \"\"\n",
" else\n",
" match n with \n",
" | 1 -> \"H 0. 0. 0.\\n\" ^ accu\n",
" | x -> toto (\"H\" ^ \" \" ^ string_of_float( d *. float_of_int(n-1)) ^ \" 0. 0.\\n\" ^ accu ) d (n-1)\n",
" in\n",
" accu\n",
" in string_of_int(n) ^ \"\\nH\" ^ string_of_int(n) ^ \"\\n\" ^ toto accu d n;;\n",
" \n",
"let xyz_string = xyz 1.8 20;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* 2. Base atomique *)\n",
"\n",
"let basis_string = \n",
"\"\n",
"HYDROGEN\n",
"S 6\n",
"1 0.3552322122E+02 0.9163596281E-02\n",
"2 0.6513143725E+01 0.4936149294E-01\n",
"3 0.1822142904E+01 0.1685383049E+00\n",
"4 0.6259552659E+00 0.3705627997E+00\n",
"5 0.2430767471E+00 0.4164915298E+00\n",
"6 0.1001124280E+00 0.1303340841E+00\n",
"\"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* II : Hartree-Fock *)\n",
"\n",
"(* 1. Def pour HF *)\n",
"\n",
"let nuclei =\n",
" Nuclei.of_xyz_string xyz_string\n",
" \n",
"let basis = \n",
" Basis.of_nuclei_and_basis_string nuclei basis_string\n",
" \n",
"let simulation = \n",
" Simulation.make ~charge:0 ~multiplicity:1 ~nuclei basis\n",
" \n",
"let ao_basis = \n",
" Simulation.ao_basis simulation\n",
" \n",
"let nocc =\n",
" let elec = Simulation.electrons simulation in\n",
" Electrons.n_alfa elec\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* 2. Calcul de Hartree-Fock*)\n",
"\n",
"let hf = HartreeFock.make simulation\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let mo_basis = MOBasis.of_hartree_fock hf\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let mo_coef = MOBasis.mo_coef mo_basis\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* III : Définition des fonctions pour la localisation *)\n",
"\n",
"(* 1. 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 multipoles = \n",
" AOBasis.multipole ao_basis;;\n",
" \n",
"let sum a = \n",
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
" \n",
"let pi = 3.14159265358979323846264338;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 2. Fonction de calcul de tous les alpha ER *)\n",
"let f_alpha m_C =\n",
"\n",
" let n_mo = Mat.dim2 m_C in\n",
" (*let t0 = Sys.time () in*)\n",
" \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",
" let v_d = Vec.init n_mo (fun i -> 0.) in\n",
" \n",
" (* Tableaux temporaires *)\n",
" let m_pqr =\n",
" Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)\n",
" in\n",
" let m_qr_i = Mat.create (n_ao*n_ao) n_mo in\n",
" let m_ri_j = Mat.create (n_ao*n_mo) n_mo in\n",
" let m_ij_k = Mat.create (n_mo*n_mo) n_mo in\n",
" \n",
" Array.iter (fun s ->\n",
" (* Grosse boucle externe sur s *)\n",
" Array.iter (fun r ->\n",
" Array.iter (fun q ->\n",
" Array.iter (fun p ->\n",
" m_pqr.{p,q,r} <- ERI.get_phys ee_ints p q r s\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" (* Conversion d'un tableau a 3 indices en une matrice nao x nao^2 *)\n",
" let m_p_qr =\n",
" Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]\n",
" |> Bigarray.array2_of_genarray\n",
" in\n",
" \n",
" let m_qr_i =\n",
" (* (qr,i) = <i r|q s> = \\sum_p <p r | q s> C_{pi} *)\n",
" gemm ~transa:`T ~c:m_qr_i m_p_qr m_C\n",
" in\n",
" \n",
" let m_q_ri =\n",
" (* Transformation de la matrice (qr,i) en (q,ri) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo)\n",
" in\n",
" \n",
" let m_ri_j =\n",
" (* (ri,j) = <i r | j s> = \\sum_q <i r | q s> C_{bj} *)\n",
" gemm ~transa:`T ~c:m_ri_j m_q_ri m_C\n",
" in\n",
" \n",
" let m_r_ij =\n",
" (* Transformation de la matrice (ri,j) en (r,ij) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo)\n",
" in\n",
" \n",
" let m_ij_k =\n",
" (* (ij,k) = <i k | j s> = \\sum_r <i r | j s> C_{rk} *)\n",
" gemm ~transa:`T ~c:m_ij_k m_r_ij m_C\n",
" in\n",
" \n",
" let m_ijk =\n",
" (* Transformation de la matrice (ei,j) en (e,ij) *)\n",
" Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |]\n",
" |> Bigarray.array3_of_genarray\n",
" in\n",
" \n",
" Array.iter (fun j ->\n",
" Array.iter (fun i ->\n",
" m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});\n",
" m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j} -.\n",
" 0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_C.{s,i} +.\n",
" (m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_C.{s,j})\n",
" ) (Util.array_range 1 n_mo);\n",
" v_d.{j} <- v_d.{j} +. m_ijk.{j,j,j} *. m_C.{s,j}\n",
" ) (Util.array_range 1 n_mo)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" (*let t1 = Sys.time () in\n",
" Printf.printf \"t = %f s\\n%!\" (t1 -. t0);*)\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",
" ),Vec.sum v_d);;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 3. Fonction de calcul de tous les alpha Boys *)\n",
"let f_alpha_boys 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",
"\n",
" let m_b12= \n",
" let b12 g = Mat.init_cols n_mo n_mo ( 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",
" g.{i,j} *. g.{i,j} -. 0.25 *. ((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",
" 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",
" 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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 4. Pattern matching : calcul de alpha et de D *)\n",
"type alphad = {\n",
"m_alpha : Mat.t;\n",
"d : float;\n",
"}\n",
"\n",
"let m_alpha_d methode m_C = \n",
" let alpha methode =\n",
" \n",
" 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",
" | \"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",
"\n",
"in \n",
"alpha methode;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 5. Norme de alpha *)\n",
"\n",
"let norme m = \n",
" let vec_m = Mat.as_vec m\n",
" in\n",
" let vec2 = Vec.sqr vec_m\n",
"in sqrt(Vec.sum vec2);;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 6. 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",
"type alphaij = {\n",
" alpha_max : float;\n",
" indice_ii : int;\n",
" indice_jj : int;};;\n",
"\n",
"let rec new_m_alpha m_alpha m_C n_rec_alpha=\n",
"\n",
" let n_mo = Mat.dim2 m_C\n",
" in\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_mo n_mo (fun i j -> \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",
" then (m_alpha.{i,j} +. ( pi /. 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 et j du alpha max *)\n",
" let indice_ii, indice_jj = \n",
" let max = max_element3 alpha_m \n",
" in\n",
" (max - 1) mod n_mo +1, (max - 1) / n_mo +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",
" alpha_m.{i,j}\n",
" \n",
" (*Printf.printf \"%i %i\\n%!\" i j;*)\n",
" \n",
" in\n",
" let alpha_max = alpha alpha_m \n",
" in\n",
" \n",
" (*Printf.printf \"%f\\n%!\" alpha_max;*)\n",
" \n",
" if alpha_max < pi /. 2.\n",
" then {alpha_max; indice_ii; indice_jj}\n",
" else new_m_alpha alpha_m m_C (n_rec_alpha-1);;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 7. Matrice de rotation 2 par 2 *)\n",
"let f_R alpha =\n",
" Mat.init_cols 2 2 (fun i j -> \n",
" if i=j \n",
" then cos alpha\n",
" else if i>j \n",
" then sin alpha \n",
" else -. sin alpha )\n",
";;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 8. 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 =\n",
" let n_ao = Mat.dim1 m_C\n",
"in\n",
"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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 9. 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 m_C = gemm m_Ksi m_R;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 10. Fonction pour la création de matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)\n",
"let f_k mat indice_i indice_j m_C = \n",
"let n_mo = Mat.dim2 m_C\n",
" in\n",
" let n_ao = Mat.dim1 m_C\n",
"in\n",
"Mat.init_cols n_ao n_mo (fun i j -> \n",
" if j=indice_i \n",
" then mat.{i,1}\n",
" else if j=indice_j \n",
" then mat.{i,2}\n",
" else 0.)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 11. 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\n",
"des coefs par la matrice *)\n",
"let f_interm m_C m_Psi = Mat.sub m_C m_Psi;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)\n",
"let new_m_C m_C methode m_alpha epsilon =\n",
" \n",
" (* alphaij contient le alpha max ainsi que ses indices i et j *)\n",
" let n_rec_alpha = 10 (* Nombre ditération max pour réduire les valeurs de alpha *)\n",
" in\n",
" let alphaij = new_m_alpha m_alpha m_C n_rec_alpha \n",
" in\n",
"\n",
" (* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)\n",
" let alpha = (alphaij.alpha_max) *. epsilon \n",
" in\n",
"\n",
" (*Printf.printf \"%f\\n%!\" alpha;*)\n",
" \n",
" (* Indice i et j du alpha max après calcul *)\n",
" let indice_i = alphaij.indice_ii \n",
" in\n",
" let indice_j = alphaij.indice_jj \n",
" in\n",
"\n",
" (*Printf.printf \"%i %i\\n%!\" indice_i indice_j;*)\n",
" \n",
" (* Matrice de rotation *)\n",
" let m_R = f_R alpha \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 \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 m_C \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_k m_Ksi indice_i indice_j m_C\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_k m_Ksi_tilde indice_i indice_j m_C \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 \n",
" in\n",
" \n",
" (*Util.debug_matrix \"m_interm\" m_interm;*)\n",
" \n",
"(* Matrice après rotation *)\n",
"Mat.add m_Psi_tilde m_interm"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"(* Localisation de Edminstion ou de Boys *)\n",
"\n",
"(* Calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)\n",
"let rec localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D=\n",
"\n",
" (*Printf.printf \"%i\\n%!\" n;\n",
" Printf.printf \"%f\\n%!\" epsilon;*)\n",
"\n",
" (*Util.debug_matrix \"m_C\" m_C;*)\n",
"\n",
" (*Util.debug_matrix \"new_alpha_m\" prev_m_alpha;*)\n",
"\n",
"if n == 0 \n",
" then m_C\n",
" else\n",
"\n",
" let m_new_m_C = new_m_C m_C methode prev_m_alpha epsilon\n",
" in\n",
" (* Fonction de pattern matching en fonction de la méthode *)\n",
" let alphad = m_alpha_d methode m_new_m_C\n",
" in\n",
" \n",
" (* D critère à maximiser *)\n",
" let critere_D = alphad.d \n",
" in\n",
" let diff = max_D -. critere_D \n",
" in\n",
" let max_D =\n",
" if diff > 0.\n",
" then max_D\n",
" else critere_D\n",
" in\n",
" (*Printf.printf \"%f\\n%!\" max_D;*)\n",
" Printf.printf \"%e\\n%!\" critere_D;\n",
" (*Printf.printf \"%f\\n%!\" diff;*)\n",
" \n",
" (* Matrice des alphas *)\n",
" let m_alpha = alphad.m_alpha\n",
" in\n",
" (*let norme_alpha = norme m_alpha\n",
" in \n",
" Printf.printf \"%f\\n%!\" norme_alpha;*)\n",
"if diff > 0. && epsilon >= 0.1\n",
" then localisation m_C methode (epsilon /. 2.) (n-1) critere_D prev_m_alpha cc max_D\n",
" else\n",
" let epsilon = 1.\n",
" in\n",
" \n",
"(*Util.debug_matrix \"new_alpha_m\" m_alpha;*)\n",
"(*Util.debug_matrix \"m_new_m_C\" m_new_m_C;*)\n",
"\n",
"if diff**2. < cc**2.\n",
" then m_new_m_C\n",
" else\n",
"localisation m_new_m_C methode epsilon (n-1) critere_D m_alpha cc max_D;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let localise localisation m_C methode epsilon n cc=\n",
" let alphad = m_alpha_d methode m_C \n",
" in\n",
" let prev_m_alpha = alphad.m_alpha\n",
" in \n",
" let prev_critere_D = 0.\n",
" in\n",
" let max_D = 0.\n",
"in\n",
"localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* V : Fonctions pour la séparation et le rassemblement des matrices *)\n",
"\n",
"(* 1. Fonction de création d'une list d'entier à partir d'un vecteur de float*)\n",
"let int_list vec = \n",
" let float_list = Vec.to_list vec\n",
" in\n",
" let g a = int_of_float(a)\n",
"in List.map g float_list;;\n",
"\n",
"(* 2. Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)\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 list_int = int_list vec \n",
"in\n",
"List.filter (fun x -> x > 0) list_int;;\n",
"\n",
"(* 3. Fonction de séparation d'une matrice en 2 sous matrice, la première matrice correspondant aux colonnes de la matrice dont le numéro est présent\n",
"dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)\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",
"(* 4. Liste des OMs occupées *)\n",
"let list_occ = \n",
" let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))\n",
"in int_list vec_occ;;\n",
"\n",
"(* 5. Fonction de rassemblement de 2 matrices *)\n",
"let assemble m_occ m_vir = Mat.init_cols n_ao n_mo ( fun i j ->\n",
" if j > nocc \n",
" then m_vir.{i,j-(n_mo - nocc)}\n",
" else m_occ.{i,j});;\n",
" \n",
"(* 6. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)\n",
"let m_occ , m_vir = split_mat m_C list_occ;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"localise localisation m_occ \"Boys\" 1. 1000 10e-5;;"
]
}
],
"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
}

800
Localisation.ml Normal file
View File

@ -0,0 +1,800 @@
(* Algorithme de localisation des orbitales *)
(* SOMMAIRE *)
(* I. Définition des types *)
(* II : Définition des fonctions pour la localisation *)
(* 1. Définitions de base nécessaire pour la suite *)
(* 2. Fonction de calcul de la matrice des angles selon la méthode de Edminston *)
(* 3. Fonction calcul de la matrice des angles selon la méthode de Boys façon Edminston *)
(* 4. Fonction de calcul de la matrice des angles selon la méthode de Boys *)
(* 5. Pattern matching : calcul de alpha et du critere D *)
(* 6. Norme d'une matrice *)
(* 7. Détermination de alpha_max et ses indices i et j.
Si alpha max > pi/2 on soustrait pi/2 à la matrice des alphas de manière récursive
Si alpha max < pi/2 on ajoute pi/2 à la matrice des alphas de manière récursive *)
(* 8. Matrice de rotation 2 par 2 *)
(* 9. Matrice de rotation n par n avec un angle fixe *)
(* 10. 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)
pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)
(* 11. Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle
on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)
(* 12. Fonction pour la création de matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)
(* 13. 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
par la matrice *)
(* III. Localisation des OMs *)
(* 1. Boucle de calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)
(* 2. Fonction de localisation *)
(* IV : Fonctions pour la séparation et le rassemblement des matrices *)
(* 1. Fonction de création d'une liste d'entier à partir d'un vecteur de float*)
(* 2. Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)
(* 3. Fonction de séparation d'une matrice en 2 sous matrice, la première matrice correspondant aux colonnes de la matrice dont le numéro est présent
dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)
(* 4. Liste des OMs occupées *)
(* 5. Liste des OMs de coeur *)
(* 6. Liste des OMs pi au sein des OMs virtuelles *)
(* 7. Fonction de rassemblement de 2 matrices *)
(* V. Calcul *)
(* 1. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)
(* 2. Application d'une légère rotation d'un angle fixé aux OMs pour améliorer la convergence *) (* !!! Réglage de l'angle de rotation !!! *)
(* 3. Légende *)
(* 4. Rotation des orbitales : localise localisation new_X / "methode" / epsilon = 1. si <1 diminue l'angle des rotations / nb itération / critère de convergence *) (* !!! Paramètres du calcul !!! *)
(* 5. Rassemblement des occupées et des virtuelles *)
let methode = "BOYS";;
let iteration = 100000;;
let cc = 10e-3;;
let epsilon = 1.;;
(*let list_split_vir2 = [];;*)
open Lacaml.D
(* I. Définition des types *)
type alphaij = {
alpha_max : float;
indice_ii : int;
indice_jj : int;
}
type alphad = {
m_alpha : Mat.t;
d : float;
}
let localize mo_basis =
let simulation = MOBasis.simulation mo_basis in
let nocc =
let elec = Simulation.electrons simulation
in
Electrons.n_alfa elec
in
let ao_basis = MOBasis.ao_basis mo_basis in
let n_core =
(Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
in
(* II : Définition des fonctions pour la localisation *)
(* 1. Définitions de base nécessaire pour la suite *)
let m_C = MOBasis.mo_coef mo_basis in
let n_ao = Mat.dim1 m_C in
let pi = acos(-1.) in
(* 2. Fonction de calcul de la matrice des angles selon la méthode de Edminston *)
let f_alpha m_C =
let ee_ints = AOBasis.ee_ints ao_basis in
let n_mo = Mat.dim2 m_C in
(*let t0 = Sys.time () in*)
let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in
let m_a12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in
let v_d = Vec.init n_mo (fun i -> 0.) in
(* Tableaux temporaires *)
let m_pqr =
Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)
in
let m_qr_i = Mat.create (n_ao*n_ao) n_mo in
let m_ri_j = Mat.create (n_ao*n_mo) n_mo in
let m_ij_k = Mat.create (n_mo*n_mo) n_mo in
Array.iter (fun s ->
(* Grosse boucle externe sur s *)
Array.iter (fun r ->
Array.iter (fun q ->
Array.iter (fun p ->
m_pqr.{p,q,r} <- ERI.get_phys ee_ints p q r s
) (Util.array_range 1 n_ao)
) (Util.array_range 1 n_ao)
) (Util.array_range 1 n_ao);
(* Conversion d'un tableau a 3 indices en une matrice nao x nao^2 *)
let m_p_qr =
Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]
|> Bigarray.array2_of_genarray
in
let m_qr_i =
(* (qr,i) = <i r|q s> = \sum_p <p r | q s> C_{pi} *)
gemm ~transa:`T ~c:m_qr_i m_p_qr m_C
in
let m_q_ri =
(* Transformation de la matrice (qr,i) en (q,ri) *)
Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo)
in
let m_ri_j =
(* (ri,j) = <i r | j s> = \sum_q <i r | q s> C_{bj} *)
gemm ~transa:`T ~c:m_ri_j m_q_ri m_C
in
let m_r_ij =
(* Transformation de la matrice (ri,j) en (r,ij) *)
Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo)
in
let m_ij_k =
(* (ij,k) = <i k | j s> = \sum_r <i r | j s> C_{rk} *)
gemm ~transa:`T ~c:m_ij_k m_r_ij m_C
in
let m_ijk =
(* Transformation de la matrice (ei,j) en (e,ij) *)
Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |]
|> Bigarray.array3_of_genarray
in
Array.iter (fun j ->
Array.iter (fun i ->
m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});
m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j} -.
0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_C.{s,i} +.
(m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_C.{s,j})
) (Util.array_range 1 n_mo);
v_d.{j} <- v_d.{j} +. m_ijk.{j,j,j} *. m_C.{s,j}
) (Util.array_range 1 n_mo)
) (Util.array_range 1 n_ao);
(*let t1 = Sys.time () in
Printf.printf "t = %f s\n%!" (t1 -. t0);*)
(Mat.init_cols n_mo n_mo ( fun i j ->
if i= j
then 0.
else if +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} <= 0.
then pi /. 4. +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j}
else -. pi /. 4. +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j})
,Vec.sum v_d)
in
(* 3. Fonction calcul de la matrice des angles selon la méthode de Boys façon Edminston *)
let f_alpha_boys m_C =
let multipoles = AOBasis.multipole ao_basis in
let n_mo = Mat.dim2 m_C
in
let phi_x_phi =
Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let phi_y_phi =
Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let phi_z_phi =
Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let m_b12=
let b12 g = Mat.init_cols n_mo n_mo ( fun i j ->
(g.{i,i} -. g.{j,j}) *. g.{i,j})
in
Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi))
in
let m_a12 =
let a12 g = Mat.init_cols n_mo n_mo ( fun i j ->
g.{i,j} *. g.{i,j} -. 0.25 *. ((g.{i,i} -. g.{j,j}) *. (g.{i,i} -. g.{j,j})))
in
Mat.add (a12 phi_x_phi) ( Mat.add (a12 phi_y_phi) (a12 phi_z_phi))
in
(Mat.init_cols n_mo n_mo ( fun i j ->
if i=j
then 0.
else if +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} >= 0.
then pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j}
else -. pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j})
,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.)))
in
(* 4. Fonction de calcul de la matrice des angles selon la méthode de Boys *)
let f_theta m_C =
let multipoles = AOBasis.multipole ao_basis in
let n_mo = Mat.dim2 m_C
in
let phi_x_phi =
Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let phi_y_phi =
Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let phi_z_phi =
Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let phi_x2_phi =
Multipole.matrix_x2 multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let phi_y2_phi =
Multipole.matrix_y2 multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let phi_z2_phi =
Multipole.matrix_z2 multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let m_b_0 =
let b_0 g h = Mat.init_cols n_mo n_mo (fun i j ->
h.{i,i} +. h.{j,j} -. (g.{i,i}**2. +. g.{j,j}**2.))
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))
in
let m_beta =
let beta g = Mat.init_cols n_mo n_mo (fun i j ->
(g.{i,i} -. g.{j,j})**2. -. 4. *. g.{i,j}**2.)
in
Mat.add (beta phi_x_phi) (Mat.add (beta phi_y_phi) (beta phi_z_phi))
in
(*Util.debug_matrix "m_beta" m_beta;*)
let m_gamma =
let gamma g = Mat.init_cols n_mo n_mo (fun i j ->
4. *. g.{i,j} *. (g.{i,i} -. g.{j,j}) )
in
Mat.add (gamma phi_x_phi) (Mat.add (gamma phi_y_phi) (gamma phi_z_phi))
in
(*Util.debug_matrix "m_gamma" m_gamma;*)
let m_theta = Mat.init_cols n_mo n_mo (fun i j ->
if i = j
then 0.
else +. 0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j})
in
(*Util.debug_matrix "m_theta" m_theta;*)
let m_critere_B = Mat.init_cols n_mo n_mo (fun i j ->
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}))
in
m_theta, Vec.sum(Mat.as_vec m_critere_B)
in
(* 5. Pattern matching : calcul de alpha et du critere D *)
let m_alpha_d methode m_C =
let alpha methode =
match methode with
| "Boys_er"
| "BOYS_ER"
| "boys_er" -> let alpha_boys , d_boys = f_alpha_boys m_C in {m_alpha = alpha_boys; d = d_boys}
| "Boys"
| "boys"
| "BOYS" -> let theta_boys, b_boys = f_theta m_C in {m_alpha = theta_boys; d = 1. /. b_boys}
| "ER"
| "er" -> let alpha_er , d_er = f_alpha m_C in {m_alpha = alpha_er; d = d_er}
| _ -> invalid_arg "Unknown method, please enter Boys or ER"
in
alpha methode
in
(* 6. Norme d'une matrice *)
let norme m =
let vec_m = Mat.as_vec m in
let vec2 = Vec.sqr vec_m
in
sqrt(Vec.sum vec2)
in
(* 7. Détermination de alpha_max et ses indices i et j.
Si alpha max > pi/2 on soustrait pi/2 à la matrice des alphas de manière récursive
Si alpha max < pi/2 on ajoute pi/2 à la matrice des alphas de manière récursive *)
let rec new_m_alpha m_alpha m_C n_rec_alpha=
let n_mo = Mat.dim2 m_C in
let alpha_m =
(*Printf.printf "%i\n%!" n_rec_alpha;*)
if n_rec_alpha == 0
then m_alpha
else
Mat.init_cols n_mo n_mo (fun i j ->
if (m_alpha.{i,j}) >= (pi /. 2.)
then (m_alpha.{i,j} -. ( pi /. 2.))
else if m_alpha.{i,j} <= -. pi /. 2.
then (m_alpha.{i,j} +. ( pi /. 2.))
else m_alpha.{i,j} )
in
(*Util.debug_matrix "alpha_m" alpha_m;*)
(* Détermination de l'emplacement du alpha max *)
let max_element alpha_m =
Mat.as_vec alpha_m
|> iamax
in
(* indice i et j du alpha max *)
let indice_ii, indice_jj =
let max = max_element alpha_m
in
(max - 1) mod n_mo +1, (max - 1) / n_mo +1
in
(* Valeur du alpha max*)
let alpha alpha_m =
let i = indice_ii in
let j = indice_jj
in
(*Printf.printf "%i %i\n%!" i j;*)
alpha_m.{i,j}
in
let alpha_max = alpha alpha_m in
(*Printf.printf "%f\n%!" alpha_max;*)
if alpha_max < pi /. 2.
then {alpha_max; indice_ii; indice_jj}
else new_m_alpha alpha_m m_C (n_rec_alpha-1)
in
(* 8. Matrice de rotation 2 par 2 *)
let f_R alpha =
Mat.init_cols 2 2 (fun i j ->
if i=j
then cos alpha
else if i>j
then sin alpha
else -. sin alpha )
in
(* 9. Matrice de rotation n par n avec un angle fixe *)
let rotate m_C angle=
let n_mo = Mat.dim2 m_C in
(* Antisymétrisation et 0 sur la diagonale*)
let m_angle angle=
Mat.init_cols n_mo n_mo (fun i j ->
if i = j
then 0.
else if i>j
then -. angle
else angle)
in
(*Util.debug_matrix "m" m;*)
let m = m_angle angle in
(* Mise au carré -> X² *)
let mm = gemm m m in
(*Util.debug_matrix "mm" mm;*)
(* Diagonalisation de X² -> diag = vecteur contenant les eigenvalues (-tau²),
m_angle2 -> Matrice avec les eigenvectors (W) *)
let diag = syev ~vectors:true mm in
(* Passage de -tau² à tau² *)
let square_tau= Vec.abs diag in
(* Passage de tau² à tau *)
let tau = Vec.sqrt square_tau in
(* Calcul de cos tau à partir du vecteur tau *)
let cos_tau = Vec.cos tau in
(* Matrice cos tau *)
let m_cos_tau =
Mat.init_cols n_mo n_mo (fun i j ->
if i=j
then cos_tau.{i}
else 0.)
in
(*Util.debug_matrix "m_cos_tau" m_cos_tau;*)
(* Calcul de sin tau à partir du vecteur tau *)
let sin_tau = Vec.sin tau in
(* Matrice de sin tau *)
let m_sin_tau =
Mat.init_cols n_mo n_mo (fun i j ->
if i=j
then sin_tau.{i}
else 0.)
in
(*Util.debug_matrix "m_sin_tau" m_sin_tau;*)
(* Calcul de la transposée de W (m_angles2) *)
let transpose_mm = Mat.transpose_copy mm in
(*Util.debug_matrix "transpose_mm" transpose_mm;*)
(* Calcul du vecteur tau^-1 *)
let tau_1 =
Vec.init n_mo (fun i ->
if tau.{i}<= 0.001
then 1.
else 1. /. tau.{i})
in
(* Calcul de la matrice tau^⁻1 *)
let m_tau_1 =
Mat.init_cols n_mo n_mo (fun i j ->
if i=j
then tau_1.{i}
else 0.)
in
(*Util.debug_matrix "m_tau_1" m_tau_1;*)
(* gemm mm (gemm m_cos_tau transpose_mm) -> W cos(tau) Wt *)
let a = gemm mm (gemm m_cos_tau transpose_mm) in
(*Util.debug_matrix "a" a;*)
(* gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) -> W tau^-1 sin(tau) Wt X *)
let b = gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) in
(*Util.debug_matrix "b" b;*)
(* Somme de 2 matrices, ici -> a + b -> R *)
let m_r = Mat.add a b (* Matrice de rotation R *)
in
gemm m_C m_r
in
(* 10. 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)
pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)
let f_Ksi indice_i indice_j m_C =
let n_ao = Mat.dim1 m_C
in
Mat.init_cols n_ao 2 (fun i j ->
if j=1
then m_C.{i,indice_i}
else m_C.{i,indice_j} )
in
(* 11. Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle
on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)
let f_Ksi_tilde m_R m_Ksi m_C =
gemm m_Ksi m_R
in
(* 12. Fonction pour la création de matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)
let f_k mat indice_i indice_j m_C =
let n_mo = Mat.dim2 m_C in
let n_ao = Mat.dim1 m_C
in
Mat.init_cols n_ao n_mo (fun i j ->
if j=indice_i
then mat.{i,1}
else if j=indice_j
then mat.{i,2}
else 0.)
in
(* 13. 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
par la matrice *)
let f_interm m_C m_Psi =
Mat.sub m_C m_Psi
in
(* III. Localisation des OMs *)
(* 1. Boucle de calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)
let rec localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D=
(*Printf.printf "%i\n%!" n;*)
(*Printf.printf "%f\n%!" epsilon;*)
(*Util.debug_matrix "m_C" m_C;
Util.debug_matrix "new_alpha_m" prev_m_alpha;*)
if n == 0
then m_C
else
(* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)
let new_m_C m_C methode m_alpha epsilon =
(* alphaij contient le alpha max ainsi que ses indices i et j *)
let n_rec_alpha = 10 in (* Nombre ditération max pour réduire les valeurs de alpha *)
let alphaij = new_m_alpha m_alpha m_C n_rec_alpha in
(* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)
let alpha = (alphaij.alpha_max) *. epsilon in
(*Printf.printf "%f\n%!" alpha;*)
(* Indice i et j du alpha max après calcul *)
let indice_i = alphaij.indice_ii in
let indice_j = alphaij.indice_jj in
(*Printf.printf "%i %i\n%!" indice_i indice_j;*)
(* Matrice de rotation *)
let m_R = f_R alpha in
(*Util.debug_matrix "m_R" m_R;*)
(* Matrice qui va subir la rotation *)
let m_Ksi = f_Ksi indice_i indice_j m_C in
(*Util.debug_matrix "m_Ksi" m_Ksi;*)
(* Matrice ayant subit la rotation *)
let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi m_C in
(*Util.debug_matrix "m_Ksi_tilde" m_Ksi_tilde;*)
(* Matrice pour supprimerles coef des orbitales i et j dans la matrice des coef *)
let m_Psi = f_k m_Ksi indice_i indice_j m_C in
(*Util.debug_matrix "m_Psi" m_Psi;*)
(* Matrice pour ajouter les coef des orbitales i~ et j~ dans la matrice des coef *)
let m_Psi_tilde = f_k m_Ksi_tilde indice_i indice_j m_C in
(*Util.debug_matrix "m_Psi_tilde" m_Psi_tilde;*)
(* Matrice avec les coef des orbitales i et j remplacés par 0 *)
let m_interm = f_interm m_C m_Psi
(*Util.debug_matrix "m_interm" m_interm;*)
in
(* Matrice après rotation *)
(Mat.add m_Psi_tilde m_interm, alpha)
in
(* Extraction de la nouvelle matrice des coef et du alpha_max *)
let m_new_m_C, alpha_max = new_m_C m_C methode prev_m_alpha epsilon in
(* Fonction de pattern matching en fonction de la méthode *)
let alphad = m_alpha_d methode m_new_m_C in
(* Matrice des alphas *)
let m_alpha = alphad.m_alpha in
(* Critères possibles de localisation -> norme_alpha, moyenne_alpha, alpha, diff *)
(* Norme de m_alpha, moyenne et max *)
let norme_alpha = norme m_alpha in
let moyenne_alpha = norme_alpha /. (float_of_int((Mat.dim1 m_C)* (Mat.dim2 m_C))) in
let alpha = alpha_max /. epsilon in
(* D critère à maximiser *)
let critere_D = alphad.d in
let diff = max_D -. critere_D in
(*Printf.printf "%f\n%!" diff;*)
Printf.printf "%i %f %f %f %f %f\n%!" n max_D critere_D alpha norme_alpha moyenne_alpha;
let max_D =
if diff > 0.
then max_D
else critere_D
in
(*Util.debug_matrix "new_alpha_m" m_alpha;*)
(*Util.debug_matrix "m_new_m_C" m_new_m_C;*)
(* Critère de conv pour sortir *)
if alpha_max**2. < cc**2.
then m_new_m_C
else localisation m_new_m_C methode epsilon (n-1) critere_D m_alpha cc max_D
in
(* 2. Fonction de localisation *)
let localise localisation m_C methode epsilon n cc =
let alphad = m_alpha_d methode m_C in
let prev_m_alpha = alphad.m_alpha in
let prev_critere_D = 0. in
let max_D = 0.
in
localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D
in
(* IV : Fonctions pour la séparation et le rassemblement des matrices *)
(* 1. Fonction de création d'une liste d'entier à partir d'un vecteur de float*)
let int_list vec =
let float_list = Vec.to_list vec in
let g a = int_of_float(a)
in
List.map g float_list
in
(* 2. Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)
let miss_elem mat list =
let n_mo = Mat.dim2 mat in
let vec = Vec.init (n_mo) (fun i ->
if List.mem i list
then 0.
else float_of_int(i))
in
let list_int = int_list vec
in
List.filter (fun x -> x > 0) list_int
in
(* 3. Fonction de séparation d'une matrice en 2 sous matrice, la première matrice correspondant aux colonnes de la matrice dont le numéro est présent
dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)
let split_mat mat list =
let vec_of_mat = Mat.to_col_vecs mat in
let f a = vec_of_mat.(a-1) in
let vec_list_1 = List.map f list in
let list_2 = miss_elem mat list in
let vec_list_2 = List.map f list_2
in
(Mat.of_col_vecs_list vec_list_1,Mat.of_col_vecs_list vec_list_2)
in
(* 4. Liste des OMs occupées *)
let list_occ =
let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))
in
int_list vec_occ
in
(* 5. Liste des OMs de coeur *)
let list_core =
let vec_core = Vec.init (n_core) (fun i -> float_of_int(i))
in int_list vec_core
in
(* 6. Liste des OMs pi au sein des OMs virtuelles *)
let list_pi mat =
let n_ao = Mat.dim1 m_C in
let n_mo = Mat.dim2 m_C in
let m_pi = Mat.init_cols n_ao n_mo ( fun i j ->
if m_C.{i,j}**2. < 10e-8**2.
then 0.
else m_C.{i,j})
in
(*Printf.printf "%i\n" nocc;
Util.debug_matrix "m_pi" m_pi;*)
let vec = Mat.to_col_vecs m_pi in
let vec_dot = Vec.init (n_mo-nocc) (fun i -> dot vec.(0) vec.(i-1+nocc))
in
let vec_pi = Vec.init (n_mo-nocc) (fun i ->
if vec_dot.{i} = 0.
then float_of_int(i)
else 0.)
in
let list_pi = int_list vec_pi in
let rec remove x list =
match list with
| [] -> []
| var :: tail -> if var = x
then remove x tail
else var :: (remove x tail)
in
remove 0 list_pi
in
(* 7. Fonction de rassemblement de 2 matrices *)
let assemble m_occ m_vir =
let v_occ = Mat.to_col_vecs m_occ in
let v_vir = Mat.to_col_vecs m_vir in
let m = Array.append v_occ v_vir
in
Mat.of_col_vecs m
in
(* V. Calcul *)
(* 1. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes, puis des orbitales de coeur : m_core et occupées : m_occu *)
let m_occ , m_vir = split_mat m_C list_occ in
(* Util.debug_matrix "m_vir" m_vir;*)
let m_core, m_occu = split_mat m_occ list_core in
let dim_m_core = Mat.dim1 m_core in
let list_split_vir2 = list_pi m_vir in
let m_vir2, m_vir1 = split_mat m_vir list_split_vir2 in
(* Util.debug_matrix "m_vir1" m_vir1;
Util.debug_matrix "m_vir2" m_vir2;*)
let dim_m_vir2 = Mat.dim1 m_vir2 in
(*Util.debug_matrix "m_vir2" m_vir2;*)
(* 2. Application d'une légère rotation d'un angle fixé aux OMs pour améliorer la convergence *)
let new_core =
if dim_m_core = 0
then m_core
else rotate m_core 0.001
in
let new_occu = rotate m_occu 0.001 in
let new_vir2 =
if dim_m_vir2 = 0
then m_vir2
else rotate m_vir2 0.001
in
let new_vir1 = rotate m_vir1 0.001 in
(* 3. Légende *)
let texte = "Methode : " ^ methode ^ "; Nombre iterarion : " ^ string_of_int(iteration) ^ "; cc : " ^ string_of_float(cc)^string_of_float(epsilon) in
let legende = "n max_D Critere_D alpha_max norme_alpha moyenne_alpha"
in
Printf.printf "%s\n" "";
Printf.printf "%s\n" "Localisation des orbitales";
Printf.printf "%s\n" "";
Printf.printf "%s\n" texte;
Printf.printf "%s\n" "";
Printf.printf "%s\n" legende;
Printf.printf "%s\n" "";
(* 4. Rotation des orbitales : localise localisation new_X / "methode" / epsilon = 1. si <1 diminue l'angle des rotations / nb itération / critère de convergence *)
Printf.printf "%s\n" "Localisation orbitales de coeurs, 1/4";
let loc_m_core =
if dim_m_core= 0
then new_core
else localise localisation new_core methode epsilon iteration cc
in
Printf.printf "%s\n" "";
Printf.printf "%s\n" "Localisation orbitales occupées, 2/4";
let loc_m_occu = localise localisation new_occu methode epsilon iteration cc
in
Printf.printf "%s\n" "";
Printf.printf "%s\n" "Localisation orbitales virtuelles pi, 3/4";
let loc_m_vir2 =
if dim_m_vir2 = 0
then new_vir2
else localise localisation new_vir2 methode epsilon iteration cc
in
Printf.printf "%s\n" "";
Printf.printf "%s\n" "Localisation orbitales virtuelles sigma, 4/4";
let loc_m_vir1 = localise localisation new_vir1 methode epsilon iteration cc
in
(* 5. Rassemblement des occupées et des virtuelles *)
let m_assemble_occ = assemble loc_m_core loc_m_occu
in
let m_assemble_vir = assemble loc_m_vir2 loc_m_vir1
in
let m_assemble_loc = assemble m_assemble_occ m_assemble_vir
in
m_assemble_loc

682
Work.ipynb
View File

@ -152,33 +152,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## H$_2$O"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(*\n",
"let xyz_string = \n",
"\"3\n",
"Water\n",
"O 0. 0. 0.\n",
"H -0.756950272703377558 0. -0.585882234512562827\n",
"H 0.756950272703377558 0. -0.585882234512562827\n",
"\"\n",
"*)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### H$_4$"
"### H$_n$"
]
},
{
@ -189,7 +163,7 @@
"source": [
"(* Fonction création chaîne linéaire de n H *)\n",
"\n",
"let xyz d n = \n",
"(*let xyz d n = \n",
" let accu = \"\"\n",
" in\n",
" let rec toto accu d n =\n",
@ -204,7 +178,20 @@
" accu\n",
" in string_of_int(n) ^ \"\\nH\" ^ string_of_int(n) ^ \"\\n\" ^ toto accu d n;;\n",
" \n",
"let xyz_string = xyz 1.8 10;;\n"
"let xyz_string = xyz 1.8 6;;\n",
"*)\n",
"let xyz_string = \n",
"\"8 \n",
"C1\n",
"C 0.000000 0.000000 0.424165\n",
"H 0.000000 0.000000 1.508704\n",
"N 0.000000 1.158087 -0.174215\n",
"N 0.000000 -1.158087 -0.174215\n",
"H 0.000000 1.258360 -1.178487\n",
"H 0.000000 2.005112 0.371143\n",
"H 0.000000 -2.005112 0.371143\n",
"H 0.000000 -1.258360 -1.178487\n",
"\""
]
},
{
@ -235,67 +222,65 @@
},
"outputs": [],
"source": [
"\n",
"let basis_string = \n",
"\"HYDROGEN\n",
"S 3\n",
"1 0.1873113696E+02 0.3349460434E-01\n",
"2 0.2825394365E+01 0.2347269535E+00\n",
"3 0.6401216923E+00 0.8137573261E+00\n",
"S 1\n",
"1 0.1612777588E+00 1.0000000\n",
"\n",
"CARBON\n",
"S 6\n",
"1 0.3552322122E+02 0.9163596281E-02\n",
"2 0.6513143725E+01 0.4936149294E-01\n",
"3 0.1822142904E+01 0.1685383049E+00\n",
"4 0.6259552659E+00 0.3705627997E+00\n",
"5 0.2430767471E+00 0.4164915298E+00\n",
"6 0.1001124280E+00 0.1303340841E+00\n",
"\"\n",
"\n",
"\n",
"(*\n",
"let basis_string = \n",
"\"\n",
"HYDROGEN\n",
"S 4\n",
"1 1.301000E+01 1.968500E-02\n",
"2 1.962000E+00 1.379770E-01\n",
"3 4.446000E-01 4.781480E-01\n",
"4 1.220000E-01 5.012400E-01\n",
"1 0.3047524880E+04 0.1834737132E-02\n",
"2 0.4573695180E+03 0.1403732281E-01\n",
"3 0.1039486850E+03 0.6884262226E-01\n",
"4 0.2921015530E+02 0.2321844432E+00\n",
"5 0.9286662960E+01 0.4679413484E+00\n",
"6 0.3163926960E+01 0.3623119853E+00\n",
"S 3\n",
"1 0.7868272350E+01 -0.1193324198E+00\n",
"2 0.1881288540E+01 -0.1608541517E+00\n",
"3 0.5442492580E+00 0.1143456438E+01\n",
"S 1\n",
"1 1.220000E-01 1.000000E+00\n",
"1 0.1687144782E+00 0.1000000000E+01\n",
"P 3\n",
"1 0.7868272350E+01 0.6899906659E-01\n",
"2 0.1881288540E+01 0.3164239610E+00\n",
"3 0.5442492580E+00 0.7443082909E+00\n",
"P 1\n",
"1 7.270000E-01 1.0000000\n",
"1 0.1687144782E+00 0.1000000000E+01\n",
"\n",
"OXYGEN\n",
"S 9\n",
"1 1.172000E+04 7.100000E-04\n",
"2 1.759000E+03 5.470000E-03\n",
"3 4.008000E+02 2.783700E-02\n",
"4 1.137000E+02 1.048000E-01\n",
"5 3.703000E+01 2.830620E-01\n",
"6 1.327000E+01 4.487190E-01\n",
"7 5.025000E+00 2.709520E-01\n",
"8 1.013000E+00 1.545800E-02\n",
"9 3.023000E-01 -2.585000E-03\n",
"S 9\n",
"1 1.172000E+04 -1.600000E-04\n",
"2 1.759000E+03 -1.263000E-03\n",
"3 4.008000E+02 -6.267000E-03\n",
"4 1.137000E+02 -2.571600E-02\n",
"5 3.703000E+01 -7.092400E-02\n",
"6 1.327000E+01 -1.654110E-01\n",
"7 5.025000E+00 -1.169550E-01\n",
"8 1.013000E+00 5.573680E-01\n",
"9 3.023000E-01 5.727590E-01\n",
"NITROGEN\n",
"S 6\n",
"1 0.4173511460E+04 0.1834772160E-02\n",
"2 0.6274579110E+03 0.1399462700E-01\n",
"3 0.1429020930E+03 0.6858655181E-01\n",
"4 0.4023432930E+02 0.2322408730E+00\n",
"5 0.1282021290E+02 0.4690699481E+00\n",
"6 0.4390437010E+01 0.3604551991E+00\n",
"S 3\n",
"1 0.1162636186E+02 -0.1149611817E+00\n",
"2 0.2716279807E+01 -0.1691174786E+00\n",
"3 0.7722183966E+00 0.1145851947E+01\n",
"S 1\n",
"1 3.023000E-01 1.000000E+00\n",
"P 4\n",
"1 1.770000E+01 4.301800E-02\n",
"2 3.854000E+00 2.289130E-01\n",
"3 1.046000E+00 5.087280E-01\n",
"4 2.753000E-01 4.605310E-01\n",
"1 0.2120314975E+00 0.1000000000E+01\n",
"P 3\n",
"1 0.1162636186E+02 0.6757974388E-01\n",
"2 0.2716279807E+01 0.3239072959E+00\n",
"3 0.7722183966E+00 0.7408951398E+00\n",
"P 1\n",
"1 2.753000E-01 1.000000E+00\n",
"D 1\n",
"1 1.185000E+00 1.0000000\n",
"1 0.2120314975E+00 0.1000000000E+01\n",
"\"\n",
"*)"
"(*\"HYDROGEN\n",
"S 3\n",
"1 0.1873113696E+02 0.3349460434E-01\n",
"2 0.2825394365E+01 0.2347269535E+00\n",
"3 0.6401216923E+00 0.8137573261E+00\n",
"S 1\n",
"1 0.1612777588E+00 1.0000000\n",
"\"*)\n"
]
},
{
@ -326,14 +311,17 @@
" Basis.of_nuclei_and_basis_string nuclei basis_string\n",
" \n",
"let simulation = \n",
" Simulation.make ~charge:0 ~multiplicity:1 ~nuclei basis\n",
" Simulation.make ~charge:1 ~multiplicity:1 ~nuclei basis\n",
" \n",
"let ao_basis = \n",
" Simulation.ao_basis simulation\n",
" \n",
"let nocc =\n",
" let elec = Simulation.electrons simulation in\n",
" Electrons.n_alfa elec"
" Electrons.n_alfa elec\n",
" \n",
"let n_core =\n",
" (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2"
]
},
{
@ -677,11 +665,13 @@
" gemm m_X m_C' (* C = X.C' *)\n",
" \n",
" \n",
"let m_C = \n",
"let m_C = m_assemble_boys\n",
"(*\n",
" match Guess.make ~nocc ~guess:`Hcore ao_basis with\n",
" | Hcore m_H -> c_of_h m_H\n",
" | _ -> assert false\n",
"*)"
" *) \n",
" *)"
]
},
{
@ -762,12 +752,13 @@
"outputs": [],
"source": [
"(*\n",
"let m_C =\n",
" let m_C', _ = \n",
"let m_C = m_assemble_boys\n",
"(* let m_C', _ = \n",
" Util.xt_o_x m_F m_X\n",
" |> Util.diagonalize_symm\n",
" in\n",
" gemm m_X m_C'\n",
"*)\n",
"*)"
]
},
@ -1326,13 +1317,11 @@
"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 multipoles = \n",
" AOBasis.multipole ao_basis;;\n",
" \n",
"let sum a = \n",
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
" \n",
"let pi = 3.14159265358979323846264338;;\n"
"let pi = acos(-1.);;\n"
]
},
{
@ -1481,14 +1470,17 @@
" Printf.printf \"t = %f s\\n%!\" (t1 -. t0);*)\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",
" ),Vec.sum v_d);;\n",
" if i= j \n",
" then 0.\n",
" else if +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} <= 0. \n",
" then pi /. 4. +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j}\n",
" else -. pi /. 4. +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j})\n",
" ,Vec.sum v_d);;\n",
"\n",
"(*********************)\n",
"\n",
"f_alpha m_C;;\n",
"let m_alpha , s_D = f_alpha m_C;;\n"
"(*\n",
"f_alpha m_C;;*)\n",
"(*let m_alpha , s_D = f_alpha m_C;;*)\n"
]
},
{
@ -1499,6 +1491,7 @@
"source": [
"(* Fonction calcul alpha Boys *)\n",
"let f_alpha_boys m_C = \n",
" let multipoles = AOBasis.multipole ao_basis in\n",
" let n_mo = Mat.dim2 m_C\n",
" in\n",
" let phi_x_phi =\n",
@ -1529,13 +1522,16 @@
" in\n",
" (Mat.init_cols n_mo n_mo ( 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",
" 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",
" then 0.\n",
" else if +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} >= 0.\n",
" then pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j}\n",
" else -. pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j})\n",
" ,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",
"*)"
]
},
{
@ -1546,61 +1542,64 @@
"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",
" let f_theta m_C = \n",
" let multipoles = AOBasis.multipole ao_basis in\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",
"\n",
"f_theta m_C;;\n",
"f_theta m_B;;"
"\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",
" let m_theta = Mat.init_cols n_mo n_mo (fun i j ->\n",
" if i = j \n",
" then 0.\n",
" else +. 0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j})\n",
" in (*Util.debug_matrix \"m_theta\" m_theta;*)\n",
"\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;;*)"
]
},
{
@ -1619,9 +1618,9 @@
" let alpha methode =\n",
" \n",
" 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",
" | \"Boys_er\"\n",
" | \"boys_er\" -> 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 = 1. /. 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",
@ -1633,7 +1632,7 @@
"(*\n",
"m_alpha_d \"ER\" ;;\n",
"m_alpha_d \"Boys\" ;;\n",
"let methode = \"ER\";;\n",
"let methode = \"BOYS\"\n",
"let alphad = m_alpha_d methode m_C;;\n",
"let m_alpha = alphad.m_alpha;;\n",
"let d = alphad.d;;\n",
@ -1688,8 +1687,6 @@
" then (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",
" else m_alpha.{i,j} )\n",
" in \n",
" \n",
@ -1731,7 +1728,7 @@
"\n",
"(*************************)\n",
"(*\n",
"let m_alpha,d = f_alpha m_C\n",
"let m_alpha,d = f_theta m_C\n",
"let alphaij = new_m_alpha m_alpha m_C 3;;\n",
"alphaij.alpha_max;;\n",
"*)"
@ -1759,6 +1756,182 @@
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Matrice de rotation n par n *)\n",
" let new_c m_C angle= \n",
" let n_mo = Mat.dim2 m_C in\n",
"\n",
" (* Antisymétrisation et 0 sur la diagonale*)\n",
" let m_angle angle= Mat.init_cols n_mo n_mo (fun i j ->\n",
" if i = j \n",
" then 0.\n",
" else if i>j \n",
" then -. angle\n",
" else angle) in\n",
"\n",
" (*Util.debug_matrix \"m\" m;*)\n",
" let m = m_angle angle in\n",
" (* Mise au carré -> X² *)\n",
" let mm = gemm m m in\n",
"\n",
" (*Util.debug_matrix \"mm\" mm;*)\n",
"\n",
" (* Diagonalisation de X² -> diag = vecteur contenant les eigenvalues (-tau²),\n",
" m_angle2 -> Matrice avec les eigenvectors (W) *)\n",
" let diag = syev ~vectors:true mm in\n",
" (* Passage de -tau² à tau² *)\n",
" let square_tau= Vec.abs diag in\n",
" (* Passage de tau² à tau *)\n",
" let tau = Vec.sqrt square_tau in \n",
" (* Calcul de cos tau à partir du vecteur tau *)\n",
" let cos_tau = Vec.cos tau in\n",
" (* Matrice cos tau *)\n",
" let m_cos_tau = Mat.init_cols n_mo n_mo (fun i j -> if i=j then cos_tau.{i} else 0.) in\n",
"\n",
" (*Util.debug_matrix \"m_cos_tau\" m_cos_tau;*)\n",
"\n",
" (* Calcul de sin tau à partir du vecteur tau *)\n",
" let sin_tau = Vec.sin tau in \n",
" (* Matrice de sin tau *)\n",
" let m_sin_tau = Mat.init_cols n_mo n_mo (fun i j -> if i=j then sin_tau.{i} else 0.) in\n",
"\n",
" (*Util.debug_matrix \"m_sin_tau\" m_sin_tau;*)\n",
"\n",
" (* Calcul de la transposée de W (m_angles2) *)\n",
" let transpose_mm = Mat.transpose_copy mm in\n",
"\n",
" (*Util.debug_matrix \"transpose_mm\" transpose_mm;*)\n",
"\n",
" (* Calcul du vecteur tau^-1 *)\n",
" let tau_1 = Vec.init n_mo (fun i -> if tau.{i}<= 0.001 then 1.\n",
" else 1. /. tau.{i}) in\n",
" (* Calcul de la matrice tau^⁻1 *)\n",
" let m_tau_1 = Mat.init_cols n_mo n_mo (fun i j -> if i=j then tau_1.{i} else 0.) in\n",
"\n",
" (*Util.debug_matrix \"m_tau_1\" m_tau_1;*)\n",
"\n",
" (* gemm mm (gemm m_cos_tau transpose_mm) -> W cos(tau) Wt *)\n",
" let a = gemm mm (gemm m_cos_tau transpose_mm) in\n",
"\n",
" (*Util.debug_matrix \"a\" a;*)\n",
"\n",
" (* gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) -> W tau^-1 sin(tau) Wt X *)\n",
" let b = gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) in\n",
"\n",
" (*Util.debug_matrix \"b\" b;*)\n",
"\n",
" (* Somme de 2 matrices, ici -> a + b -> R *)\n",
" let m_r = Mat.add a b (* Matrice de rotation R *)\n",
" in gemm m_C m_r\n",
";;\n",
"\n",
"(************************************************)\n",
"\n",
"(*\n",
"f_alpha_boys m_C;;\n",
"f_alpha_boys new_c;;\n",
"*)\n",
"\n",
"(********************************************)\n",
"\n",
"let f_r alphaij m_C = \n",
" let n_mo = Mat.dim2 m_C in\n",
" \n",
" let indice_i = alphaij.indice_ii in\n",
" let indice_j = alphaij.indice_jj in\n",
" let alpha = alphaij.alpha_max in\n",
" (* Matrice des angles -> X *)\n",
" let m_angle = Mat.init_cols n_mo n_mo (fun i j -> \n",
" if (i = indice_i && j = indice_j) || (j = indice_i && i = indice_j) \n",
" then alpha\n",
" else 0.) in\n",
" \n",
" (*Util.debug_matrix \"m_angle\" m_angle;*)\n",
" \n",
" (*\n",
" let m_alpha = alphad.m_alpha in\n",
" let m_angle = Mat.init_cols n_mo n_mo (fun i j -> \n",
" 0.1 *. m_alpha.{i,j})\n",
" in\n",
" *)\n",
" (* Antisymétrisation et 0 sur la diagonale*)\n",
" let m =Mat.init_cols n_mo n_mo (fun i j ->\n",
" if i = j \n",
" then 0.\n",
" else if i>j \n",
" then m_angle.{i,j} \n",
" else -. m_angle.{i,j})\n",
" \n",
" in\n",
" \n",
" (*Util.debug_matrix \"m\" m;*)\n",
" \n",
" (* Mise au carré -> X² *)\n",
" let mm = gemm m m in\n",
" \n",
" (*Util.debug_matrix \"mm\" mm;*)\n",
" \n",
" (* Diagonalisation de X² -> diag = vecteur contenant les eigenvalues (-tau²),\n",
" m_angle2 -> Matrice avec les eigenvectors (W) *)\n",
" let diag = syev ~vectors:true mm in\n",
" (* Passage de -tau² à tau² *)\n",
" let square_tau= Vec.abs diag in\n",
" (* Passage de tau² à tau *)\n",
" let tau = Vec.sqrt square_tau in \n",
" (* Calcul de cos tau à partir du vecteur tau *)\n",
" let cos_tau = Vec.cos tau in\n",
" (* Matrice cos tau *)\n",
" let m_cos_tau = Mat.init_cols n_mo n_mo (fun i j -> if i=j then cos_tau.{i} else 0.) in\n",
" \n",
" (*Util.debug_matrix \"m_cos_tau\" m_cos_tau;*)\n",
" \n",
" (* Calcul de sin tau à partir du vecteur tau *)\n",
" let sin_tau = Vec.sin tau in \n",
" (* Matrice de sin tau *)\n",
" let m_sin_tau = Mat.init_cols n_mo n_mo (fun i j -> if i=j then sin_tau.{i} else 0.) in\n",
" \n",
" (*Util.debug_matrix \"m_sin_tau\" m_sin_tau;*)\n",
" \n",
" (* Calcul de la transposée de W (m_angles2) *)\n",
" let transpose_mm = Mat.transpose_copy mm in\n",
" \n",
" (*Util.debug_matrix \"transpose_mm\" transpose_mm;*)\n",
" \n",
" (* Calcul du vecteur tau^-1 *)\n",
" let tau_1 = Vec.init n_mo (fun i -> if tau.{i}<= 0.001 then 1.\n",
" else 1. /. tau.{i}) in\n",
" (* Calcul de la matrice tau^⁻1 *)\n",
" let m_tau_1 = Mat.init_cols n_mo n_mo (fun i j -> if i=j then tau_1.{i} else 0.) in\n",
" \n",
" (*Util.debug_matrix \"m_tau_1\" m_tau_1;*)\n",
" \n",
" (* gemm mm (gemm m_cos_tau transpose_mm) -> W cos(tau) Wt *)\n",
" let a = gemm mm (gemm m_cos_tau transpose_mm) in\n",
" \n",
" (*Util.debug_matrix \"a\" a;*)\n",
" \n",
" (* gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) -> W tau^-1 sin(tau) Wt X *)\n",
" let b = gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) in\n",
" \n",
" (*Util.debug_matrix \"b\" b;*)\n",
" \n",
" (* Somme de 2 matrices, ici -> a + b -> R *)\n",
"Mat.add a b (* Matrice de rotation R *)\n",
";;\n",
"\n",
"(*******************************************************************) \n",
"\n",
"(*\n",
"let m_r = f_r alphaij m_C;;\n",
"gemm m_C m_r;;\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
@ -1872,8 +2045,8 @@
"\n",
" (*Printf.printf \"%i\\n%!\" n;*)\n",
" (*Printf.printf \"%f\\n%!\" epsilon;*)\n",
" (*Util.debug_matrix \"m_C\" m_C;*)\n",
" (*Util.debug_matrix \"new_alpha_m\" prev_m_alpha;*)\n",
" (*Util.debug_matrix \"m_C\" m_C;\n",
" Util.debug_matrix \"new_alpha_m\" prev_m_alpha;*)\n",
"\n",
"if n == 0 \n",
" then m_C\n",
@ -1881,97 +2054,89 @@
"\n",
" (* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)\n",
" let new_m_C m_C methode m_alpha epsilon =\n",
"\n",
" (*\n",
" let alphaij = new_m_alpha m_alpha m_C 3 in\n",
" let m_r = f_r alphaij m_C in\n",
" (*Util.debug_matrix \"m_r\" m_r;*)\n",
" let new_m = gemm m_C m_r in\n",
" (*Util.debug_matrix \"new_m\" new_m;*)\n",
" \n",
" new_m, alphaij.alpha_max\n",
" *)\n",
" \n",
" (* alphaij contient le alpha max ainsi que ses indices i et j *)\n",
" let n_rec_alpha = 10 (* Nombre ditération max pour réduire les valeurs de alpha *)\n",
" in\n",
" let alphaij = new_m_alpha m_alpha m_C n_rec_alpha \n",
" in\n",
" let n_rec_alpha = 10 in (* Nombre ditération max pour réduire les valeurs de alpha *)\n",
" let alphaij = new_m_alpha m_alpha m_C n_rec_alpha in\n",
"\n",
" (* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)\n",
" let alpha = (alphaij.alpha_max) *. epsilon \n",
" in\n",
"\n",
" let alpha = (alphaij.alpha_max) *. epsilon in\n",
" \n",
" (*Printf.printf \"%f\\n%!\" alpha;*)\n",
"\n",
" (* Indice i et j du alpha max après calcul *)\n",
" let indice_i = alphaij.indice_ii \n",
" in\n",
" let indice_j = alphaij.indice_jj \n",
" in\n",
" let indice_i = alphaij.indice_ii in\n",
" let indice_j = alphaij.indice_jj in\n",
"\n",
" (*Printf.printf \"%i %i\\n%!\" indice_i indice_j;*)\n",
"\n",
" \n",
" (* Matrice de rotation *)\n",
" let m_R = f_R alpha \n",
" in\n",
" let m_R = f_R alpha in\n",
"\n",
" (*Util.debug_matrix \"m_R\" m_R;*)\n",
"\n",
" \n",
" (* Matrice qui va subir la rotation *)\n",
" let m_Ksi = f_Ksi indice_i indice_j m_C \n",
" in\n",
" let m_Ksi = f_Ksi indice_i indice_j m_C 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 m_C \n",
" in\n",
" let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi m_C 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_k m_Ksi indice_i indice_j m_C\n",
" in\n",
" let m_Psi = f_k m_Ksi indice_i indice_j m_C 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_k m_Ksi_tilde indice_i indice_j m_C \n",
" in\n",
" let m_Psi_tilde = f_k m_Ksi_tilde indice_i indice_j m_C 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 \n",
" in\n",
" let m_interm = f_interm m_C m_Psi in\n",
"\n",
" (*Util.debug_matrix \"m_interm\" m_interm;*)\n",
"\n",
" (* Matrice après rotation *)\n",
" (Mat.add m_Psi_tilde m_interm, alpha)\n",
"\n",
" \n",
" in\n",
" (* Extraction de la nouvelle matrice des coef et du alpha_max *)\n",
" let m_new_m_C, alpha_max = new_m_C m_C methode prev_m_alpha epsilon\n",
" in\n",
" let m_new_m_C, alpha_max = new_m_C m_C methode prev_m_alpha epsilon in\n",
"\n",
" (* Fonction de pattern matching en fonction de la méthode *)\n",
" let alphad = m_alpha_d methode m_new_m_C\n",
" in\n",
" let alphad = m_alpha_d methode m_new_m_C in\n",
" \n",
" (* Matrice des alphas *)\n",
" let m_alpha = alphad.m_alpha\n",
" in\n",
" let m_alpha = alphad.m_alpha in\n",
" \n",
" (* 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",
" 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",
" \n",
" let norme_alpha = norme m_alpha in\n",
" let moyenne_alpha = norme_alpha /. (float_of_int((Mat.dim1 m_C)* (Mat.dim2 m_C))) in\n",
" let alpha = alpha_max /. epsilon in\n",
" \n",
" \n",
" (* D critère à maximiser *)\n",
" let critere_D = alphad.d \n",
" in\n",
" let diff = max_D -. critere_D \n",
" in\n",
" let critere_D = alphad.d in\n",
" let diff = max_D -. critere_D 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",
@ -1991,7 +2156,7 @@
"\n",
"\n",
" (* Critère de conv pour sortir *)\n",
" if alpha_max**2. +. 1. < cc**2.\n",
" if alpha_max**2. < 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;;"
]
@ -2003,17 +2168,14 @@
"outputs": [],
"source": [
"let localise localisation m_C methode epsilon n cc=\n",
" let alphad = m_alpha_d methode m_C \n",
" in\n",
" let prev_m_alpha = alphad.m_alpha\n",
" in \n",
" let prev_critere_D = 0.\n",
" in\n",
" let alphad = m_alpha_d methode m_C in\n",
" let prev_m_alpha = alphad.m_alpha in \n",
" let prev_critere_D = 0. in\n",
" let max_D = 0.\n",
"in\n",
"localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D;;\n",
" \n",
"let m_B = localise localisation m_C \"BOYS\" 1. 100 10e-10;;"
"(*let new_mc = new_c m_C 0.001;;\n",
"let m_B = localise localisation new_mc \"Boys\" 1. 200 10e-10;;*)"
]
},
{
@ -2044,7 +2206,7 @@
"\n",
"(* Fonction de séparation d'une matrice en 2 sous matrice, la première matrice correspondant aux colonnes de la matrice dont le numéro est présent\n",
"dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)\n",
"let split_mat mat list =\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",
@ -2061,6 +2223,11 @@
" let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))\n",
"in int_list vec_occ;;\n",
"\n",
"(* Liste des OMs de coeur *)\n",
"let list_core = \n",
" let vec_core = Vec.init (n_core) (fun i -> float_of_int(i))\n",
"in int_list vec_core;;\n",
"\n",
"(* Fonction de rassemblement de 2 matrices *)\n",
"let assemble m_occ m_vir =\n",
" let v_occ = Mat.to_col_vecs m_occ\n",
@ -2089,7 +2256,15 @@
"metadata": {},
"outputs": [],
"source": [
"let m_occ , m_vir = split_mat m_C list_occ;;"
"let m_occ , m_vir = split_mat m_C list_occ;;\n",
"let dim_m_core = Mat.dim1 m_core;;\n",
"let m_core, m_occu = split_mat m_occ list_core;;\n",
"let new_core = \n",
" if dim_m_core = 0 \n",
" then m_core \n",
" else new_c m_core 0.001;;\n",
"let new_occu = new_c m_occu 0.001;;\n",
"let new_vir = new_c m_vir 0.001;;\n"
]
},
{
@ -2104,13 +2279,13 @@
"(*let new_m_boys = localise localisation m_C \"boys\" 1. 100 10e-7;;\n",
"let new_m_er = localise localisation m_C \"ER\" 1. 100 10e-7;;\n",
"*)\n",
"\n",
"let loc_m_occ_boys = localise localisation m_occ \"Boys\" 1. 1000 10e-7;;\n",
"\n",
"(*\n",
"let loc_m_vir_boys = localise localisation m_vir \"boys\" 1. 100 10e-7;;\n",
"let m_assemble_boys = assemble loc_m_occ_boys loc_m_vir_boys;;\n",
"*)\n",
"let loc_m_core_boys = \n",
" if dim_m_core=0 \n",
" then m_core \n",
" else localise localisation m_core \"boys_er\" 1. 10 10e-3;;\n",
"let loc_m_occ_boys = localise localisation new_occ \"boys_er\" 1. 10 10e-3;;\n",
"let loc_m_vir_boys = localise localisation new_vir \"BOYS\" 1. 10 10e-3;;\n",
"(*let m_assemble_boys = assemble loc_m_occ_boys loc_m_vir_boys;;*)\n",
"\n",
"(*\n",
"let loc_m_occ_er = localise localisation m_occ \"ER\" 1. 100 10e-7;;\n",
@ -2213,11 +2388,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"jupyter": {
"outputs_hidden": true
}
},
"metadata": {},
"outputs": [],
"source": [
"(*\n",
@ -2298,6 +2469,55 @@
"print_list print_int [3;6;78;5;2;34;7];;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let mat = mo_coef;;\n",
" let n_ao = Mat.dim1 mat ;;\n",
" let n_mo = Mat.dim2 mat ;;\n",
" let m_pi = Mat.init_cols n_ao n_mo ( fun i j ->\n",
" if mat.{i,j}**2. < 10e-8**2.\n",
" then 0.\n",
" else mat.{i,j})\n",
" ;;\n",
" Printf.printf \"%i\\n\" nocc;;\n",
" Util.debug_matrix \"m_pi\" m_pi;;\n",
" let vec = Mat.to_col_vecs m_pi;;\n",
" let vec_dot = Vec.init (n_mo-nocc) (fun i -> dot vec.(0) vec.(i-1+nocc))\n",
" ;;\n",
" let vec_pi = Vec.init (n_mo-nocc) (fun i ->\n",
" if vec_dot.{i} = 0.\n",
" then float_of_int(i)\n",
" else 0.);;\n",
" \n",
" let list_pi = int_list vec_pi;;\n",
" let rec remove x list =\n",
" match list with\n",
" | [] -> []\n",
" | var :: tail -> if var = x\n",
" then remove x tail\n",
" else var :: (remove x tail);;\n",
" remove 0 list_pi\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,

1754
Work_test.ipynb
View File

File diff suppressed because it is too large Load Diff

118
run_Localisation.ml Normal file
View File

@ -0,0 +1,118 @@
open Lacaml.D
let read_qp_mo dirname =
let ic = Unix.open_process_in ("zcat "^dirname^"/mo_basis/mo_coef.gz") in
let check = String.trim (input_line ic) in
assert (check = "2");
let int_list =
input_line ic
|> String.split_on_char ' '
|> List.filter (fun x -> x <> "")
|> List.map int_of_string
in
let n_ao, n_mo =
match int_list with
| [ x ; y ] -> x, y
| _ -> assert false
in
let result =
Mat.init_cols n_ao n_mo (fun i j ->
let s = input_line ic in
Scanf.sscanf s " %f " (fun x -> x)
)
in
let exit_code = Unix.close_process_in ic in
assert (exit_code = Unix.WEXITED 0);
result
let () =
let open Command_line in
begin
set_header_doc (Sys.argv.(0) ^ " - QCaml command");
set_description_doc "Localizes MOs";
set_specs
[ { short='b' ; long="basis" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the basis set"; } ;
{ short='x' ; long="xyz" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
{ short='m' ; long="multiplicity" ; opt=Optional;
arg=With_arg "<int>";
doc="Spin multiplicity (2S+1). Default is singlet"; } ;
{ short='c' ; long="charge" ; opt=Optional;
arg=With_arg "<int>";
doc="Total charge of the molecule. Default is 0"; } ;
{ short='o' ; long="output" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the localized MOs"; } ;
{ short='i' ; long="import" ; opt=Optional;
arg=With_arg "<string>";
doc="Name of the EZFIO directory containing MOs"; } ;
]
end;
let basis_file = Util.of_some @@ Command_line.get "basis" in
let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
let ezfio_file = Command_line.get "import" in
let charge =
match Command_line.get "charge" with
| Some x -> int_of_string x
| None -> 0
in
let multiplicity =
match Command_line.get "multiplicity" with
| Some x -> int_of_string x
| None -> 1
in
let output_filename = Util.of_some @@ Command_line.get "output" in
(* II : Hartree-Fock *)
(* 1. Def pour HF *)
let simulation =
Simulation.of_filenames
~charge ~multiplicity ~nuclei:nuclei_file basis_file
in
(* 2. Calcul de Hartree-Fock*)
let mo_basis =
match ezfio_file with
| Some ezfio_file ->
let mo_coef = read_qp_mo ezfio_file in
let mo_type = MOBasis.Localized "Boys" in
let elec = Simulation.electrons simulation in
let n_mo = Mat.dim2 mo_coef in
let mo_occupation = Vec.init n_mo (fun i ->
if i <= Electrons.n_beta elec then 2.
else if i <= Electrons.n_alfa elec then 1.
else 0.) in
MOBasis.make ~simulation ~mo_type ~mo_occupation ~mo_coef ()
| None -> HartreeFock.make simulation
|> MOBasis.of_hartree_fock
in
let local_mos = Localisation.localize mo_basis in
let oc = open_out output_filename in
Printf.fprintf oc "[\n";
Mat.as_vec local_mos
|> Vec.iter (fun x -> Printf.fprintf oc "%20.15e,\n" x);
Printf.fprintf oc "]\n";
close_out oc

913
run_cyanine.ipynb
View File

@ -1,913 +0,0 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* Algorithme de localisation des orbitales : Méthode de Boys ou Edminston-Ruedenberg *)\n",
"\n",
"(*I : Position et base atomique*)\n",
" (* 1. Coordonnées *)\n",
" (* 2. Base atomique *)\n",
"\n",
"(* II : Hartree-Fock *)\n",
" (* 1. Def pour HF *) (* !!! -> Charge !!! *)\n",
" (* 2. Calcul de Hartree-Fock*)\n",
" \n",
"(* III : Définition des fonctions pour la localisation *)\n",
" (* 1. Définitions de base nécessaire pour la suite *) \n",
" (* 2. Fonction de calcul de tous les alpha ER *)\n",
" (* 3. Fonction de calcul de tous les alpha Boys *)\n",
" (* 4. Pattern matching : calcul de alpha et de D *)\n",
" (* 5. Norme de alpha *)\n",
" (* 6. Détermination alpha_max et ses indices i et j *) \n",
" (* 7. Matrice de rotation 2 par 2 *)\n",
" (* 8. Fonction : i,j -> matrice Ksi (n par 2) *)\n",
" (* 9. Fonction : i~,j~ -> matrice Ksi~ (n par 2) *)\n",
" (* 10. Matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)\n",
" (* 11. Matrice intérmédiaire où les orbitales i et j ont été supprimées et remplacées par des 0*)\n",
" \n",
"(* IV : Localisation des OMs *)\n",
" (* 1. Boucle de localisation des OMs *) (* !!! -> Choix paramètres et critères localisation !!! *)\n",
" (* 2. Fonction de localisation *)\n",
"\n",
"\n",
"(* V : Fonctions pour la séparation et le rassemblement des matrices *)\n",
" (* 1. Fonction de création d'une list d'entier à partir d'un vecteur de float*)\n",
" (* 2. Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)\n",
" (* 3. Fonction de séparation d'une matrice en 2 sous matrice *)\n",
" (* 4. Liste des OMs occupées *)\n",
" (* 5. Fonction de rassemblement de 2 matrices *) \n",
" (* 6. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)\n",
" \n",
"(* VI : Calcul et assemblage des OMs occupées/virtuels *) (* !!! -> Choix paramètres de calcul !!! *)\n",
" \n",
"(* VII : Calcul Coef CI *)\n",
" (* 1. Def de la mo_basis pour le HF pre CI *)\n",
" (* 2. HF pour CI *)\n",
" (* 3. Calcul de coef CI *)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"open Lacaml.D\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*I : Position et base atomique*)\n",
"\n",
"(* 1. Coordonnées *) (* !!! Attention à bien mettre \"nbatome sur la même ligne !!! *)\n",
"let charge = 1;;\n",
"let xyz_string = \n",
"\"8\n",
"C1\n",
"C 0.000000 0.000000 0.424165\n",
"H 0.000000 0.000000 1.508704\n",
"N 0.000000 1.158087 -0.174215\n",
"N 0.000000 -1.158087 -0.174215\n",
"H 0.000000 1.258360 -1.178487\n",
"H 0.000000 2.005112 0.371143\n",
"H 0.000000 -2.005112 0.371143\n",
"H 0.000000 -1.258360 -1.178487\n",
"\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* 2. Base atomique -> 6-31G *)\n",
"\n",
"let basis_string = \n",
"\"\n",
"HYDROGEN\n",
"S 3\n",
"1 0.1873113696E+02 0.3349460434E-01\n",
"2 0.2825394365E+01 0.2347269535E+00\n",
"3 0.6401216923E+00 0.8137573261E+00\n",
"S 1\n",
"1 0.1612777588E+00 1.0000000\n",
"\n",
"CARBON\n",
"S 6\n",
"1 0.3047524880E+04 0.1834737132E-02\n",
"2 0.4573695180E+03 0.1403732281E-01\n",
"3 0.1039486850E+03 0.6884262226E-01\n",
"4 0.2921015530E+02 0.2321844432E+00\n",
"5 0.9286662960E+01 0.4679413484E+00\n",
"6 0.3163926960E+01 0.3623119853E+00\n",
"S 3\n",
"1 0.7868272350E+01 -0.1193324198E+00\n",
"2 0.1881288540E+01 -0.1608541517E+00\n",
"3 0.5442492580E+00 0.1143456438E+01\n",
"S 1\n",
"1 0.1687144782E+00 0.1000000000E+01\n",
"P 3\n",
"1 0.7868272350E+01 0.6899906659E-01\n",
"2 0.1881288540E+01 0.3164239610E+00\n",
"3 0.5442492580E+00 0.7443082909E+00\n",
"P 1\n",
"1 0.1687144782E+00 0.1000000000E+01\n",
"\n",
"NITROGEN\n",
"S 6\n",
"1 0.4173511460E+04 0.1834772160E-02\n",
"2 0.6274579110E+03 0.1399462700E-01\n",
"3 0.1429020930E+03 0.6858655181E-01\n",
"4 0.4023432930E+02 0.2322408730E+00\n",
"5 0.1282021290E+02 0.4690699481E+00\n",
"6 0.4390437010E+01 0.3604551991E+00\n",
"S 3\n",
"1 0.1162636186E+02 -0.1149611817E+00\n",
"2 0.2716279807E+01 -0.1691174786E+00\n",
"3 0.7722183966E+00 0.1145851947E+01\n",
"S 1\n",
"1 0.2120314975E+00 0.1000000000E+01\n",
"P 3\n",
"1 0.1162636186E+02 0.6757974388E-01\n",
"2 0.2716279807E+01 0.3239072959E+00\n",
"3 0.7722183966E+00 0.7408951398E+00\n",
"P 1\n",
"1 0.2120314975E+00 0.1000000000E+01\n",
"\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* II : Hartree-Fock *)\n",
"\n",
"(* 1. Def pour HF *)\n",
"\n",
"let nuclei =\n",
" Nuclei.of_xyz_string xyz_string\n",
" \n",
"let basis = \n",
" Basis.of_nuclei_and_basis_string nuclei basis_string\n",
" \n",
"let simulation = \n",
" Simulation.make ~charge ~multiplicity:1 ~nuclei basis \n",
" \n",
"let ao_basis = \n",
" Simulation.ao_basis simulation\n",
" \n",
"let nocc =\n",
" let elec = Simulation.electrons simulation in\n",
" Electrons.n_alfa elec\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* 2. Calcul de Hartree-Fock*)\n",
"\n",
"let hf = HartreeFock.make simulation\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let mo_basis = MOBasis.of_hartree_fock hf\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let mo_coef = MOBasis.mo_coef mo_basis\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* III : Définition des fonctions pour la localisation *)\n",
"\n",
"(* 1. 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 multipoles = \n",
" AOBasis.multipole ao_basis;;\n",
" \n",
"let sum a = \n",
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
" \n",
"let pi = 3.14159265358979323846264338;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 2. Fonction de calcul de tous les alpha ER *)\n",
"let f_alpha m_C =\n",
"\n",
" let n_mo = Mat.dim2 m_C in\n",
" (*let t0 = Sys.time () in*)\n",
" \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",
" let v_d = Vec.init n_mo (fun i -> 0.) in\n",
" \n",
" (* Tableaux temporaires *)\n",
" let m_pqr =\n",
" Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)\n",
" in\n",
" let m_qr_i = Mat.create (n_ao*n_ao) n_mo in\n",
" let m_ri_j = Mat.create (n_ao*n_mo) n_mo in\n",
" let m_ij_k = Mat.create (n_mo*n_mo) n_mo in\n",
" \n",
" Array.iter (fun s ->\n",
" (* Grosse boucle externe sur s *)\n",
" Array.iter (fun r ->\n",
" Array.iter (fun q ->\n",
" Array.iter (fun p ->\n",
" m_pqr.{p,q,r} <- ERI.get_phys ee_ints p q r s\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" (* Conversion d'un tableau a 3 indices en une matrice nao x nao^2 *)\n",
" let m_p_qr =\n",
" Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]\n",
" |> Bigarray.array2_of_genarray\n",
" in\n",
" \n",
" let m_qr_i =\n",
" (* (qr,i) = <i r|q s> = \\sum_p <p r | q s> C_{pi} *)\n",
" gemm ~transa:`T ~c:m_qr_i m_p_qr m_C\n",
" in\n",
" \n",
" let m_q_ri =\n",
" (* Transformation de la matrice (qr,i) en (q,ri) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo)\n",
" in\n",
" \n",
" let m_ri_j =\n",
" (* (ri,j) = <i r | j s> = \\sum_q <i r | q s> C_{bj} *)\n",
" gemm ~transa:`T ~c:m_ri_j m_q_ri m_C\n",
" in\n",
" \n",
" let m_r_ij =\n",
" (* Transformation de la matrice (ri,j) en (r,ij) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo)\n",
" in\n",
" \n",
" let m_ij_k =\n",
" (* (ij,k) = <i k | j s> = \\sum_r <i r | j s> C_{rk} *)\n",
" gemm ~transa:`T ~c:m_ij_k m_r_ij m_C\n",
" in\n",
" \n",
" let m_ijk =\n",
" (* Transformation de la matrice (ei,j) en (e,ij) *)\n",
" Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |]\n",
" |> Bigarray.array3_of_genarray\n",
" in\n",
" \n",
" Array.iter (fun j ->\n",
" Array.iter (fun i ->\n",
" m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});\n",
" m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j} -.\n",
" 0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_C.{s,i} +.\n",
" (m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_C.{s,j})\n",
" ) (Util.array_range 1 n_mo);\n",
" v_d.{j} <- v_d.{j} +. m_ijk.{j,j,j} *. m_C.{s,j}\n",
" ) (Util.array_range 1 n_mo)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" (*let t1 = Sys.time () in\n",
" Printf.printf \"t = %f s\\n%!\" (t1 -. t0);*)\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",
" ),Vec.sum v_d);;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 3. Fonction de calcul de tous les alpha Boys *)\n",
"\n",
"let f_alpha_boys 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",
"\n",
" let m_b12= \n",
" let b12 g = Mat.init_cols n_mo n_mo ( 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",
" g.{i,j} *. g.{i,j} -. 0.25 *. ((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",
" 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",
" 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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 4. Pattern matching : calcul de alpha et de D *)\n",
"\n",
"type alphad = {\n",
"m_alpha : Mat.t;\n",
"d : float;\n",
"}\n",
"\n",
"let m_alpha_d methode m_C = \n",
" let alpha methode =\n",
" \n",
" 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",
" | \"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",
"\n",
"in \n",
"alpha methode;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 5. Norme de alpha *)\n",
"\n",
"let norme m = \n",
" let vec_m = Mat.as_vec m\n",
" in\n",
" let vec2 = Vec.sqr vec_m\n",
"in sqrt(Vec.sum vec2);;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 6. 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",
"type alphaij = {\n",
" alpha_max : float;\n",
" indice_ii : int;\n",
" indice_jj : int;};;\n",
"\n",
"let rec new_m_alpha m_alpha m_C n_rec_alpha=\n",
"\n",
" let n_mo = Mat.dim2 m_C\n",
" in\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_mo n_mo (fun i j -> \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",
" then (m_alpha.{i,j} +. ( pi /. 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 et j du alpha max *)\n",
" let indice_ii, indice_jj = \n",
" let max = max_element3 alpha_m \n",
" in\n",
" (max - 1) mod n_mo +1, (max - 1) / n_mo +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",
" alpha_m.{i,j}\n",
" \n",
" (*Printf.printf \"%i %i\\n%!\" i j;*)\n",
" \n",
" in\n",
" let alpha_max = alpha alpha_m \n",
" in\n",
" \n",
" (*Printf.printf \"%f\\n%!\" alpha_max;*)\n",
" \n",
" if alpha_max < pi /. 2.\n",
" then {alpha_max; indice_ii; indice_jj}\n",
" else new_m_alpha alpha_m m_C (n_rec_alpha-1);;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 7. Matrice de rotation 2 par 2 *)\n",
"let f_R alpha =\n",
" Mat.init_cols 2 2 (fun i j -> \n",
" if i=j \n",
" then cos alpha\n",
" else if i>j \n",
" then sin alpha \n",
" else -. sin alpha )\n",
";;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 8. 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 =\n",
" let n_ao = Mat.dim1 m_C\n",
"in\n",
"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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 9. 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 m_C = gemm m_Ksi m_R;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 10. Fonction pour la création de matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)\n",
"let f_k mat indice_i indice_j m_C = \n",
"let n_mo = Mat.dim2 m_C\n",
" in\n",
" let n_ao = Mat.dim1 m_C\n",
"in\n",
"Mat.init_cols n_ao n_mo (fun i j -> \n",
" if j=indice_i \n",
" then mat.{i,1}\n",
" else if j=indice_j \n",
" then mat.{i,2}\n",
" else 0.)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 11. 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\n",
"des coefs par la matrice *)\n",
"let f_interm m_C m_Psi = Mat.sub m_C m_Psi;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* IV : Localisation des OMs *)\n",
"\n",
"\n",
"(* 1. Boucle de calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)\n",
"let rec localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D=\n",
"\n",
" (*Printf.printf \"%i\\n%!\" n;*)\n",
" (*Printf.printf \"%f\\n%!\" epsilon;*)\n",
" (*Util.debug_matrix \"m_C\" m_C;*)\n",
" (*Util.debug_matrix \"new_alpha_m\" prev_m_alpha;*)\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 methode m_alpha epsilon =\n",
"\n",
" (* alphaij contient le alpha max ainsi que ses indices i et j *)\n",
" let n_rec_alpha = 10 (* Nombre ditération max pour réduire les valeurs de alpha *)\n",
" in\n",
" let alphaij = new_m_alpha m_alpha m_C n_rec_alpha \n",
" in\n",
"\n",
" (* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)\n",
" let alpha = (alphaij.alpha_max) *. epsilon \n",
" in\n",
"\n",
" (*Printf.printf \"%f\\n%!\" alpha;*)\n",
"\n",
" (* Indice i et j du alpha max après calcul *)\n",
" let indice_i = alphaij.indice_ii \n",
" in\n",
" let indice_j = alphaij.indice_jj \n",
" in\n",
"\n",
" (*Printf.printf \"%i %i\\n%!\" indice_i indice_j;*)\n",
"\n",
" (* Matrice de rotation *)\n",
" let m_R = f_R alpha \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 \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 m_C \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_k m_Ksi indice_i indice_j m_C\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_k m_Ksi_tilde indice_i indice_j m_C \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 \n",
" in\n",
"\n",
" (*Util.debug_matrix \"m_interm\" m_interm;*)\n",
"\n",
" (* Matrice après rotation *)\n",
" (Mat.add m_Psi_tilde m_interm, alpha)\n",
"\n",
" in\n",
" (* Extraction de la nouvelle matrice des coef et du alpha_max *)\n",
" let m_new_m_C, alpha_max = new_m_C m_C methode prev_m_alpha epsilon\n",
" in\n",
"\n",
" (* Fonction de pattern matching en fonction de la méthode *)\n",
" let alphad = m_alpha_d methode m_new_m_C\n",
" in\n",
" \n",
" (* Matrice des alphas *)\n",
" let m_alpha = alphad.m_alpha\n",
" in\n",
" \n",
" (* 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",
" 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",
" \n",
" (* D critère à maximiser *)\n",
" let critere_D = alphad.d \n",
" in\n",
" let diff = max_D -. critere_D \n",
" 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",
" \n",
" (* Stabilisation de la convergence *)\n",
"\n",
" let max_D =\n",
" if diff > 0.\n",
" then max_D\n",
" else critere_D\n",
" in\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",
" (*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 alpha**2. < 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;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 2. Fonction de localisation *)\n",
"let localise localisation m_C methode epsilon n cc=\n",
" let alphad = m_alpha_d methode m_C \n",
" in\n",
" let prev_m_alpha = alphad.m_alpha\n",
" in \n",
" let prev_critere_D = 0.\n",
" in\n",
" let max_D = 0.\n",
"in\n",
"localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D;;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* V : Fonctions pour la séparation et le rassemblement des matrices *)\n",
"\n",
"(* 1. Fonction de création d'une list d'entier à partir d'un vecteur de float*)\n",
"let int_list vec = \n",
" let float_list = Vec.to_list vec\n",
" in\n",
" let g a = int_of_float(a)\n",
"in List.map g float_list;;\n",
"\n",
"(* 2. Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)\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 list_int = int_list vec \n",
"in\n",
"List.filter (fun x -> x > 0) list_int;;\n",
"\n",
"(* 3. Fonction de séparation d'une matrice en 2 sous matrice, la première matrice correspondant aux colonnes de la matrice \n",
"dont le numéro est présent dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)\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",
"(* 4. Liste des OMs occupées *)\n",
"let list_occ = \n",
" let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))\n",
"in int_list vec_occ;;\n",
"\n",
"(* 5. Fonction de rassemblement de 2 matrices *)\n",
"let assemble m_occ m_vir =\n",
" let occ = Mat.to_col_vecs m_occ in\n",
" let vir = Mat.to_col_vecs m_vir in\n",
" Array.concat [ occ ; vir ]\n",
" |> Mat.of_col_vecs\n",
" \n",
"(* 6. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)\n",
"let m_occ , m_vir = split_mat m_C list_occ;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* VI : Calcul et assemblage des OMs occupées/virtuels *)\n",
"\n",
"(* Calcul : \n",
"localise localisation / Matrice des coef / \"Methode\" (boys/er) / nb itération max / critère de convergence \n",
"*)\n",
"\n",
"(*let new_m_boys = localise localisation m_C \"boys\" 1. 100 10e-7;;\n",
"let new_m_er = localise localisation m_C \"ER\" 1. 100 10e-7;;\n",
"*)\n",
"(* Legende au dessus des valeurs print *)\n",
"let texte = \"n max_D Critere_D alpha_max norme_alpha moyenne_alpha\"\n",
"in Printf.printf \"%s\\n\" texte;;\n",
"\n",
"let loc_m_occ_boys = localise localisation m_occ \"Boys\" 1. 10000 10e-3;;\n",
"let loc_m_vir_boys = localise localisation m_vir \"boys\" 1. 10000 10e-3;;\n",
"let m_assemble_boys = assemble loc_m_occ_boys loc_m_vir_boys;;\n",
"\n",
"Printf.printf \"[\\n\";;\n",
"Mat.as_vec m_assemble_boys\n",
"|> Vec.iter (fun x -> Printf.printf \"%20.15e,\\n\" x);;\n",
"Printf.printf \"]\\n\";;\n",
"\n",
"(*\n",
"let loc_m_occ_er = localise localisation m_occ \"ER\" 1. 100 10e-7;;\n",
"let loc_m_vir_er = localise localisation m_vir \"ER\" 1. 100 10e-7;;\n",
"let m_assemble_er = assemble loc_m_occ_er loc_m_vir_er;;*)\n",
"(*\n",
"let loc_assemble_boys = localise localisation m_assemble_boys \"boys\" 1. 100 10e-7;;\n",
"let loc_assemble_er = localise localisation m_assemble_er \"er\" 1. 100 10e-7;;*)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"(* VII : Calcul Coef CI *)\n",
"\n",
"(* 1. Def de la mo_basis pour le HF pre CI *)\n",
"\n",
"(*\n",
"let mo_base1 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"Boys\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:new_m_boys ();;\n",
"\n",
"let mo_base2 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"ER\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:new_m_er ();;\n",
"*)\n",
"\n",
"let mo_base3 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"Boys\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:m_assemble_boys ();;\n",
"\n",
"(*\n",
"let mo_base4 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"ER\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:m_assemble_er ();;\n",
"\n",
"let mo_base5 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"Boys\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:loc_assemble_boys ();;\n",
"\n",
"let mo_base6 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"ER\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:loc_assemble_er ();;\n",
"*)\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"(* 2. HF pour CI *)\n",
"\n",
"let nuclei =\n",
" Nuclei.of_xyz_string xyz_string\n",
" \n",
"let basis = \n",
" Basis.of_nuclei_and_basis_string nuclei basis_string\n",
" \n",
"let simulation = \n",
" Simulation.make ~charge:0 ~multiplicity:1 ~nuclei basis\n",
" \n",
"let ao_basis = \n",
" Simulation.ao_basis simulation\n",
" \n",
"let hf = HartreeFock.make simulation\n",
"\n",
"(*let mo_basis = MOBasis.of_hartree_fock hf*)\n",
"let mo_basis = mo_base3\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"(* 3. Calcul de coef CI *)\n",
"\n",
"let ci =\n",
" DeterminantSpace.fci_of_mo_basis mo_basis ~frozen_core:false\n",
" |> CI.make\n",
" \n",
"let ci_coef, ci_energy = Lazy.force ci.eigensystem ;;\n",
"let m_ci_coef = Mat.as_vec ci_coef;;\n",
"Vec.iteri (fun i x -> Printf.printf \"%d %e\\n\" i x) m_ci_coef;;\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
}

903
run_loc.ipynb
View File

@ -1,903 +0,0 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* Algorithme de localisation des orbitales : Méthode de Boys ou Edminston-Ruedenberg *)\n",
"\n",
"(*I : Position et base atomique*)\n",
" (* 1. Coordonnées *)\n",
" (* 2. Base atomique *)\n",
"\n",
"(* II : Hartree-Fock *)\n",
" (* 1. Def pour HF *) (* !!! -> Charge !!! *)\n",
" (* 2. Calcul de Hartree-Fock*)\n",
" \n",
"(* III : Définition des fonctions pour la localisation *)\n",
" (* 1. Définitions de base nécessaire pour la suite *) \n",
" (* 2. Fonction de calcul de tous les alpha ER *)\n",
" (* 3. Fonction de calcul de tous les alpha Boys *)\n",
" (* 4. Pattern matching : calcul de alpha et de D *)\n",
" (* 5. Norme de alpha *)\n",
" (* 6. Détermination alpha_max et ses indices i et j *) \n",
" (* 7. Matrice de rotation 2 par 2 *)\n",
" (* 8. Fonction : i,j -> matrice Ksi (n par 2) *)\n",
" (* 9. Fonction : i~,j~ -> matrice Ksi~ (n par 2) *)\n",
" (* 10. Matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)\n",
" (* 11. Matrice intérmédiaire où les orbitales i et j ont été supprimées et remplacées par des 0*)\n",
" \n",
"(* IV : Localisation des OMs *)\n",
" (* 1. Boucle de localisation des OMs *) (* !!! -> Choix paramètres et critères localisation !!! *)\n",
" (* 2. Fonction de localisation *)\n",
"\n",
"\n",
"(* V : Fonctions pour la séparation et le rassemblement des matrices *)\n",
" (* 1. Fonction de création d'une list d'entier à partir d'un vecteur de float*)\n",
" (* 2. Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)\n",
" (* 3. Fonction de séparation d'une matrice en 2 sous matrice *)\n",
" (* 4. Liste des OMs occupées *)\n",
" (* 5. Fonction de rassemblement de 2 matrices *) \n",
" (* 6. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)\n",
" \n",
"(* VI : Calcul et assemblage des OMs occupées/virtuels *) (* !!! -> Choix paramètres de calcul !!! *)\n",
" \n",
"(* VII : Calcul Coef CI *)\n",
" (* 1. Def de la mo_basis pour le HF pre CI *)\n",
" (* 2. HF pour CI *)\n",
" (* 3. Calcul de coef CI *)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"open Lacaml.D\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*I : Position et base atomique*)\n",
"\n",
"(* 1. Coordonnées *)\n",
"let xyz_string =\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* 2. Base atomique -> 6-31G *)\n",
"\n",
"let basis_string = \n",
"\"\n",
"HYDROGEN\n",
"S 3\n",
"1 0.1873113696E+02 0.3349460434E-01\n",
"2 0.2825394365E+01 0.2347269535E+00\n",
"3 0.6401216923E+00 0.8137573261E+00\n",
"S 1\n",
"1 0.1612777588E+00 1.0000000\n",
"\n",
"CARBON\n",
"S 6\n",
"1 0.3047524880E+04 0.1834737132E-02\n",
"2 0.4573695180E+03 0.1403732281E-01\n",
"3 0.1039486850E+03 0.6884262226E-01\n",
"4 0.2921015530E+02 0.2321844432E+00\n",
"5 0.9286662960E+01 0.4679413484E+00\n",
"6 0.3163926960E+01 0.3623119853E+00\n",
"S 3\n",
"1 0.7868272350E+01 -0.1193324198E+00\n",
"2 0.1881288540E+01 -0.1608541517E+00\n",
"3 0.5442492580E+00 0.1143456438E+01\n",
"S 1\n",
"1 0.1687144782E+00 0.1000000000E+01\n",
"P 3\n",
"1 0.7868272350E+01 0.6899906659E-01\n",
"2 0.1881288540E+01 0.3164239610E+00\n",
"3 0.5442492580E+00 0.7443082909E+00\n",
"P 1\n",
"1 0.1687144782E+00 0.1000000000E+01\n",
"\n",
"NITROGEN\n",
"S 6\n",
"1 0.4173511460E+04 0.1834772160E-02\n",
"2 0.6274579110E+03 0.1399462700E-01\n",
"3 0.1429020930E+03 0.6858655181E-01\n",
"4 0.4023432930E+02 0.2322408730E+00\n",
"5 0.1282021290E+02 0.4690699481E+00\n",
"6 0.4390437010E+01 0.3604551991E+00\n",
"S 3\n",
"1 0.1162636186E+02 -0.1149611817E+00\n",
"2 0.2716279807E+01 -0.1691174786E+00\n",
"3 0.7722183966E+00 0.1145851947E+01\n",
"S 1\n",
"1 0.2120314975E+00 0.1000000000E+01\n",
"P 3\n",
"1 0.1162636186E+02 0.6757974388E-01\n",
"2 0.2716279807E+01 0.3239072959E+00\n",
"3 0.7722183966E+00 0.7408951398E+00\n",
"P 1\n",
"1 0.2120314975E+00 0.1000000000E+01\n",
"\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* II : Hartree-Fock *)\n",
"\n",
"(* 1. Def pour HF *)\n",
"\n",
"let nuclei =\n",
" Nuclei.of_xyz_string xyz_string\n",
" \n",
"let basis = \n",
" Basis.of_nuclei_and_basis_string nuclei basis_string\n",
" \n",
"let simulation = \n",
" Simulation.make ~charge:0 ~multiplicity:1 ~nuclei basis \n",
" \n",
"let ao_basis = \n",
" Simulation.ao_basis simulation\n",
" \n",
"let nocc =\n",
" let elec = Simulation.electrons simulation in\n",
" Electrons.n_alfa elec\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"(* 2. Calcul de Hartree-Fock*)\n",
"\n",
"let hf = HartreeFock.make simulation\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let mo_basis = MOBasis.of_hartree_fock hf\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let mo_coef = MOBasis.mo_coef mo_basis\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* III : Définition des fonctions pour la localisation *)\n",
"\n",
"(* 1. 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 multipoles = \n",
" AOBasis.multipole ao_basis;;\n",
" \n",
"let sum a = \n",
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
" \n",
"let pi = 3.14159265358979323846264338;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 2. Fonction de calcul de tous les alpha ER *)\n",
"let f_alpha m_C =\n",
"\n",
" let n_mo = Mat.dim2 m_C in\n",
" (*let t0 = Sys.time () in*)\n",
" \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",
" let v_d = Vec.init n_mo (fun i -> 0.) in\n",
" \n",
" (* Tableaux temporaires *)\n",
" let m_pqr =\n",
" Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)\n",
" in\n",
" let m_qr_i = Mat.create (n_ao*n_ao) n_mo in\n",
" let m_ri_j = Mat.create (n_ao*n_mo) n_mo in\n",
" let m_ij_k = Mat.create (n_mo*n_mo) n_mo in\n",
" \n",
" Array.iter (fun s ->\n",
" (* Grosse boucle externe sur s *)\n",
" Array.iter (fun r ->\n",
" Array.iter (fun q ->\n",
" Array.iter (fun p ->\n",
" m_pqr.{p,q,r} <- ERI.get_phys ee_ints p q r s\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" (* Conversion d'un tableau a 3 indices en une matrice nao x nao^2 *)\n",
" let m_p_qr =\n",
" Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]\n",
" |> Bigarray.array2_of_genarray\n",
" in\n",
" \n",
" let m_qr_i =\n",
" (* (qr,i) = <i r|q s> = \\sum_p <p r | q s> C_{pi} *)\n",
" gemm ~transa:`T ~c:m_qr_i m_p_qr m_C\n",
" in\n",
" \n",
" let m_q_ri =\n",
" (* Transformation de la matrice (qr,i) en (q,ri) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo)\n",
" in\n",
" \n",
" let m_ri_j =\n",
" (* (ri,j) = <i r | j s> = \\sum_q <i r | q s> C_{bj} *)\n",
" gemm ~transa:`T ~c:m_ri_j m_q_ri m_C\n",
" in\n",
" \n",
" let m_r_ij =\n",
" (* Transformation de la matrice (ri,j) en (r,ij) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo)\n",
" in\n",
" \n",
" let m_ij_k =\n",
" (* (ij,k) = <i k | j s> = \\sum_r <i r | j s> C_{rk} *)\n",
" gemm ~transa:`T ~c:m_ij_k m_r_ij m_C\n",
" in\n",
" \n",
" let m_ijk =\n",
" (* Transformation de la matrice (ei,j) en (e,ij) *)\n",
" Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |]\n",
" |> Bigarray.array3_of_genarray\n",
" in\n",
" \n",
" Array.iter (fun j ->\n",
" Array.iter (fun i ->\n",
" m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});\n",
" m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j} -.\n",
" 0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_C.{s,i} +.\n",
" (m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_C.{s,j})\n",
" ) (Util.array_range 1 n_mo);\n",
" v_d.{j} <- v_d.{j} +. m_ijk.{j,j,j} *. m_C.{s,j}\n",
" ) (Util.array_range 1 n_mo)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" (*let t1 = Sys.time () in\n",
" Printf.printf \"t = %f s\\n%!\" (t1 -. t0);*)\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",
" ),Vec.sum v_d);;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 3. Fonction de calcul de tous les alpha Boys *)\n",
"\n",
"let f_alpha_boys 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",
"\n",
" let m_b12= \n",
" let b12 g = Mat.init_cols n_mo n_mo ( 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",
" g.{i,j} *. g.{i,j} -. 0.25 *. ((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",
" 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",
" 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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 4. Pattern matching : calcul de alpha et de D *)\n",
"\n",
"type alphad = {\n",
"m_alpha : Mat.t;\n",
"d : float;\n",
"}\n",
"\n",
"let m_alpha_d methode m_C = \n",
" let alpha methode =\n",
" \n",
" 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",
" | \"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",
"\n",
"in \n",
"alpha methode;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 5. Norme de alpha *)\n",
"\n",
"let norme m = \n",
" let vec_m = Mat.as_vec m\n",
" in\n",
" let vec2 = Vec.sqr vec_m\n",
"in sqrt(Vec.sum vec2);;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 6. 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",
"type alphaij = {\n",
" alpha_max : float;\n",
" indice_ii : int;\n",
" indice_jj : int;};;\n",
"\n",
"let rec new_m_alpha m_alpha m_C n_rec_alpha=\n",
"\n",
" let n_mo = Mat.dim2 m_C\n",
" in\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_mo n_mo (fun i j -> \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",
" then (m_alpha.{i,j} +. ( pi /. 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 et j du alpha max *)\n",
" let indice_ii, indice_jj = \n",
" let max = max_element3 alpha_m \n",
" in\n",
" (max - 1) mod n_mo +1, (max - 1) / n_mo +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",
" alpha_m.{i,j}\n",
" \n",
" (*Printf.printf \"%i %i\\n%!\" i j;*)\n",
" \n",
" in\n",
" let alpha_max = alpha alpha_m \n",
" in\n",
" \n",
" (*Printf.printf \"%f\\n%!\" alpha_max;*)\n",
" \n",
" if alpha_max < pi /. 2.\n",
" then {alpha_max; indice_ii; indice_jj}\n",
" else new_m_alpha alpha_m m_C (n_rec_alpha-1);;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 7. Matrice de rotation 2 par 2 *)\n",
"let f_R alpha =\n",
" Mat.init_cols 2 2 (fun i j -> \n",
" if i=j \n",
" then cos alpha\n",
" else if i>j \n",
" then sin alpha \n",
" else -. sin alpha )\n",
";;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 8. 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 =\n",
" let n_ao = Mat.dim1 m_C\n",
"in\n",
"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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 9. 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 m_C = gemm m_Ksi m_R;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 10. Fonction pour la création de matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)\n",
"let f_k mat indice_i indice_j m_C = \n",
"let n_mo = Mat.dim2 m_C\n",
" in\n",
" let n_ao = Mat.dim1 m_C\n",
"in\n",
"Mat.init_cols n_ao n_mo (fun i j -> \n",
" if j=indice_i \n",
" then mat.{i,1}\n",
" else if j=indice_j \n",
" then mat.{i,2}\n",
" else 0.)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 11. 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\n",
"des coefs par la matrice *)\n",
"let f_interm m_C m_Psi = Mat.sub m_C m_Psi;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* IV : Localisation des OMs *)\n",
"\n",
"\n",
"(* 1. Boucle de calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)\n",
"let rec localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D=\n",
"\n",
" (*Printf.printf \"%i\\n%!\" n;*)\n",
" (*Printf.printf \"%f\\n%!\" epsilon;*)\n",
" (*Util.debug_matrix \"m_C\" m_C;*)\n",
" (*Util.debug_matrix \"new_alpha_m\" prev_m_alpha;*)\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 methode m_alpha epsilon =\n",
"\n",
" (* alphaij contient le alpha max ainsi que ses indices i et j *)\n",
" let n_rec_alpha = 10 (* Nombre ditération max pour réduire les valeurs de alpha *)\n",
" in\n",
" let alphaij = new_m_alpha m_alpha m_C n_rec_alpha \n",
" in\n",
"\n",
" (* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)\n",
" let alpha = (alphaij.alpha_max) *. epsilon \n",
" in\n",
"\n",
" (*Printf.printf \"%f\\n%!\" alpha;*)\n",
"\n",
" (* Indice i et j du alpha max après calcul *)\n",
" let indice_i = alphaij.indice_ii \n",
" in\n",
" let indice_j = alphaij.indice_jj \n",
" in\n",
"\n",
" (*Printf.printf \"%i %i\\n%!\" indice_i indice_j;*)\n",
"\n",
" (* Matrice de rotation *)\n",
" let m_R = f_R alpha \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 \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 m_C \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_k m_Ksi indice_i indice_j m_C\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_k m_Ksi_tilde indice_i indice_j m_C \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 \n",
" in\n",
"\n",
" (*Util.debug_matrix \"m_interm\" m_interm;*)\n",
"\n",
" (* Matrice après rotation *)\n",
" (Mat.add m_Psi_tilde m_interm, alpha)\n",
"\n",
" in\n",
" (* Extraction de la nouvelle matrice des coef et du alpha_max *)\n",
" let m_new_m_C, alpha_max = new_m_C m_C methode prev_m_alpha epsilon\n",
" in\n",
"\n",
" (* Fonction de pattern matching en fonction de la méthode *)\n",
" let alphad = m_alpha_d methode m_new_m_C\n",
" in\n",
" \n",
" (* Matrice des alphas *)\n",
" let m_alpha = alphad.m_alpha\n",
" in\n",
" \n",
" (* 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",
" 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",
" \n",
" (* D critère à maximiser *)\n",
" let critere_D = alphad.d \n",
" in\n",
" let diff = max_D -. critere_D \n",
" 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",
" \n",
" (* Stabilisation de la convergence *)\n",
"\n",
" let max_D =\n",
" if diff > 0.\n",
" then max_D\n",
" else critere_D\n",
" in\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",
" (*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 alpha**2. < 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;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 2. Fonction de localisation *)\n",
"let localise localisation m_C methode epsilon n cc=\n",
" let alphad = m_alpha_d methode m_C \n",
" in\n",
" let prev_m_alpha = alphad.m_alpha\n",
" in \n",
" let prev_critere_D = 0.\n",
" in\n",
" let max_D = 0.\n",
"in\n",
"localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D;;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* V : Fonctions pour la séparation et le rassemblement des matrices *)\n",
"\n",
"(* 1. Fonction de création d'une list d'entier à partir d'un vecteur de float*)\n",
"let int_list vec = \n",
" let float_list = Vec.to_list vec\n",
" in\n",
" let g a = int_of_float(a)\n",
"in List.map g float_list;;\n",
"\n",
"(* 2. Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)\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 list_int = int_list vec \n",
"in\n",
"List.filter (fun x -> x > 0) list_int;;\n",
"\n",
"(* 3. Fonction de séparation d'une matrice en 2 sous matrice, la première matrice correspondant aux colonnes de la matrice \n",
"dont le numéro est présent dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)\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",
"(* 4. Liste des OMs occupées *)\n",
"let list_occ = \n",
" let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))\n",
"in int_list vec_occ;;\n",
"\n",
"(* 5. Fonction de rassemblement de 2 matrices *)\n",
"let assemble m_occ m_vir =\n",
" let occ = Mat.to_col_vecs m_occ in\n",
" let vir = Mat.to_col_vecs m_vir in\n",
" Array.concat [ occ ; vir ]\n",
" |> Mat.of_col_vecs\n",
" \n",
"(* 6. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)\n",
"let m_occ , m_vir = split_mat m_C list_occ;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* VI : Calcul et assemblage des OMs occupées/virtuels *)\n",
"\n",
"(* Calcul : \n",
"localise localisation / Matrice des coef / \"Methode\" (boys/er) / nb itération max / critère de convergence \n",
"*)\n",
"\n",
"(*let new_m_boys = localise localisation m_C \"boys\" 1. 100 10e-7;;\n",
"let new_m_er = localise localisation m_C \"ER\" 1. 100 10e-7;;\n",
"*)\n",
"(* Legende au dessus des valeurs print *)\n",
"let texte = \"n max_D Critere_D alpha_max norme_alpha moyenne_alpha\"\n",
"in Printf.printf \"%s\\n\" texte;;\n",
"\n",
"let loc_m_occ_boys = localise localisation m_occ \"Boys\" 1. 10000 10e-3;;\n",
"let loc_m_vir_boys = localise localisation m_vir \"boys\" 1. 10000 10e-3;;\n",
"let m_assemble_boys = assemble loc_m_occ_boys loc_m_vir_boys;;\n",
"\n",
"\n",
"Printf.printf \"[\\n\";;\n",
"Mat.as_vec m_assemble_boys\n",
"|> Vec.iter (fun x -> Printf.printf \"%20.15e,\\n\" x);;\n",
"Printf.printf \"]\\n\";;\n",
"\n",
"(*\n",
"let loc_m_occ_er = localise localisation m_occ \"ER\" 1. 100 10e-7;;\n",
"let loc_m_vir_er = localise localisation m_vir \"ER\" 1. 100 10e-7;;\n",
"let m_assemble_er = assemble loc_m_occ_er loc_m_vir_er;;*)\n",
"(*\n",
"let loc_assemble_boys = localise localisation m_assemble_boys \"boys\" 1. 100 10e-7;;\n",
"let loc_assemble_er = localise localisation m_assemble_er \"er\" 1. 100 10e-7;;*)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"(* VII : Calcul Coef CI *)\n",
"\n",
"(* 1. Def de la mo_basis pour le HF pre CI *)\n",
"\n",
"(*\n",
"let mo_base1 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"Boys\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:new_m_boys ();;\n",
"\n",
"let mo_base2 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"ER\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:new_m_er ();;\n",
"*)\n",
"\n",
"let mo_base3 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"Boys\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:m_assemble_boys ();;\n",
"\n",
"(*\n",
"let mo_base4 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"ER\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:m_assemble_er ();;\n",
"\n",
"let mo_base5 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"Boys\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:loc_assemble_boys ();;\n",
"\n",
"let mo_base6 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized \"ER\")\n",
"~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:loc_assemble_er ();;\n",
"*)\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"(* 2. HF pour CI *)\n",
"\n",
"let nuclei =\n",
" Nuclei.of_xyz_string xyz_string\n",
" \n",
"let basis = \n",
" Basis.of_nuclei_and_basis_string nuclei basis_string\n",
" \n",
"let simulation = \n",
" Simulation.make ~charge:0 ~multiplicity:1 ~nuclei basis\n",
" \n",
"let ao_basis = \n",
" Simulation.ao_basis simulation\n",
" \n",
"let hf = HartreeFock.make simulation\n",
"\n",
"(*let mo_basis = MOBasis.of_hartree_fock hf*)\n",
"let mo_basis = mo_base3\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"(* 3. Calcul de coef CI *)\n",
"\n",
"let ci =\n",
" DeterminantSpace.fci_of_mo_basis mo_basis ~frozen_core:false\n",
" |> CI.make\n",
" \n",
"let ci_coef, ci_energy = Lazy.force ci.eigensystem ;;\n",
"let m_ci_coef = Mat.as_vec ci_coef;;\n",
"Vec.iteri (fun i x -> Printf.printf \"%d %e\\n\" i x) m_ci_coef;;\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
}