Work_test : nouvelle boucle avec variation du pas

This commit is contained in:
Yann Damour 2020-05-21 17:10:28 +02:00
parent a44126a743
commit e29ba94200
4 changed files with 857 additions and 42 deletions

754
Calcul_test.ipynb Normal file
View File

@ -0,0 +1,754 @@
{
"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
}

View File

@ -161,14 +161,40 @@
},
"outputs": [],
"source": [
"let xyz_string = \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": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 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 10;;"
]
},
{
@ -199,7 +225,24 @@
},
"outputs": [],
"source": [
"let basis_string = \"\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",
"\n",
"\n",
"\n",
"\n",
"\n",
"(*let basis_string = \"\n",
"HYDROGEN\n",
"S 4\n",
"1 1.301000E+01 1.968500E-02\n",
@ -244,7 +287,8 @@
"D 1\n",
"1 1.185000E+00 1.0000000\n",
"\n",
"\""
"\"\n",
"*)"
]
},
{

View File

@ -320,6 +320,29 @@
"hidden": true
},
"outputs": [],
"source": [
"(*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*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let nuclei =\n",
" Nuclei.of_xyz_string xyz_string\n",
@ -1818,7 +1841,7 @@
" in\n",
" (*let norme_alpha = norme m_alpha\n",
" in*)\n",
" Util.debug_matrix \"new_alpha_m\" m_alpha;\n",
" (*Util.debug_matrix \"new_alpha_m\" m_alpha;*)\n",
" (*Printf.printf \"%f\\n%!\" norme_alpha;*)\n",
" \n",
" (*Util.debug_matrix \"m_alpha\" m_alpha;*)\n",
@ -1833,7 +1856,7 @@
" let alpha = (alphaij.alpha_max) *. epsilon \n",
" in\n",
"\n",
" Printf.printf \"%f\\n%!\" alpha;\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",
@ -1890,7 +1913,7 @@
"in\n",
"\n",
"(*Util.debug_matrix \"new_alpha_m\" (f_alpha m_C);*)\n",
"Util.debug_matrix \"m_new_m_C\" m_new_m_C;\n",
"(*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",
@ -1908,7 +1931,7 @@
"(* Fonction / Matrice des coef / Méthode(\"Boys\" ou \"ER\") / Pas(<=1.) \n",
"/ Nombre d'itérations max / 0. (valeur de D pour initier la boucle) / critère de convergence sur D*)\n",
"\n",
"let new_m_boys = localisation m_C \"boys\" 1. 5 0. 10e-7;;\n",
"let new_m_boys = localisation m_C \"boys\" 1. 20 0. 10e-7;;\n",
"(*let new_m_er = localisation m_C \"ER\" 1. 100 0. 10e-7;;\n",
"*)\n"
]

View File

@ -204,7 +204,7 @@
" 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 20;;\n"
]
},
{
@ -305,15 +305,6 @@
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let nuclei =\n",
" Nuclei.of_xyz_string xyz_string\n",
@ -1467,7 +1458,7 @@
" let alpha = (alphaij.alpha_max) *. epsilon \n",
" in\n",
"\n",
" Printf.printf \"%f\\n%!\" alpha;\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",
@ -1528,13 +1519,14 @@
"(* 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=\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 \"%i\\n%!\" n;\n",
" Printf.printf \"%f\\n%!\" epsilon;*)\n",
"\n",
" Util.debug_matrix \"m_C\" m_C;\n",
" (*Util.debug_matrix \"m_C\" m_C;*)\n",
"\n",
" Util.debug_matrix \"new_alpha_m\" prev_m_alpha;\n",
" (*Util.debug_matrix \"new_alpha_m\" prev_m_alpha;*)\n",
"\n",
"if n == 0 \n",
" then m_C\n",
@ -1549,8 +1541,16 @@
" (* D critère à maximiser *)\n",
" let critere_D = alphad.d \n",
" in\n",
" \n",
" Printf.printf \"%f\\n%!\" critere_D;\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",
@ -1558,18 +1558,19 @@
" (*let norme_alpha = norme m_alpha\n",
" in \n",
" Printf.printf \"%f\\n%!\" norme_alpha;*)\n",
"if diff > 0. && epsilon > 0.09\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",
" let diff = prev_critere_D -. critere_D\n",
"\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;;\n"
"localisation m_new_m_C methode epsilon (n-1) critere_D m_alpha cc max_D;;\n"
]
},
{
@ -1578,26 +1579,19 @@
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"\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;;\n",
"localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D;;\n",
"\n",
"localise localisation m_C \"Boys\" 1. 100 10e-7;;\n",
"*)\n",
"m_C;;\n",
"\n",
"let alphad = m_alpha_d \"boys\" m_C \n",
"\n",
"let prev_m_alpha = alphad.m_alpha\n",
"let prev_critere_D = 0.;;\n",
"\n",
"localisation m_C \"Boys\" 1. 5 prev_critere_D prev_m_alpha 10e-7;;"
"localise localisation m_C \"Boys\" 1. 1000 10e-5;;\n"
]
},
{