diff --git a/Calcul_test.ipynb b/Calcul_test.ipynb new file mode 100644 index 0000000..398ed53 --- /dev/null +++ b/Calcul_test.ipynb @@ -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) = = \\sum_p

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) = = \\sum_q 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) = = \\sum_r 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 +} diff --git a/HartreeFock.ipynb b/HartreeFock.ipynb index aee73a4..766b23e 100644 --- a/HartreeFock.ipynb +++ b/HartreeFock.ipynb @@ -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", + "*)" ] }, { diff --git a/Work.ipynb b/Work.ipynb index 0c14d29..1476d3e 100644 --- a/Work.ipynb +++ b/Work.ipynb @@ -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" ] diff --git a/Copy_test.ipynb b/Work_test.ipynb similarity index 98% rename from Copy_test.ipynb rename to Work_test.ipynb index 0fe4b17..5a93336 100644 --- a/Copy_test.ipynb +++ b/Work_test.ipynb @@ -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" ] }, {