work : calcul indice j avec i

This commit is contained in:
Yann Damour 2020-05-18 16:04:18 +02:00
parent 19ffb39db5
commit 961954547a
2 changed files with 67 additions and 93 deletions

View File

@ -125,7 +125,7 @@
" \n",
"let hf = HartreeFock.make simulation\n",
"\n",
"let mo_basis = MOBasis.of_hartree_fock hf"
"let mo_basis = MOBasis.of_hartree_fock hf;;"
]
},
{

View File

@ -1299,11 +1299,6 @@
" \n",
"let sum a = \n",
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
" \n",
"type alphaij = {\n",
" alpha_max : float;\n",
" indice_ii : int;\n",
" indice_jj : int;};;\n",
" \n",
"let pi = 3.14159265358979323846264338;;\n"
]
@ -1528,14 +1523,8 @@
"\n",
" let alpha_boys , d_boys = f_alpha_boys m_C\n",
" in\n",
" (*let d_boys = d_boys m_C\n",
" in*)\n",
" let alpha_er , d_er = f_alpha m_C\n",
" in\n",
" (*let alpha_er = mat_alpha m_C\n",
" in\n",
" let d_er = s_D m_C\n",
" in*)\n",
" let alpha methode =\n",
" match methode with \n",
" | \"Boys\"\n",
@ -1584,6 +1573,11 @@
"metadata": {},
"outputs": [],
"source": [
"type alphaij = {\n",
" alpha_max : float;\n",
" indice_ii : int;\n",
" indice_jj : int;};;\n",
"\n",
"(* Détermination alpha_max et ses indices i et j.\n",
"Si alpha max > pi/2 on soustrait pi/2 à la matrice des alphas de manière récursive *)\n",
"let rec new_m_alpha m_alpha m_C n_rec_alpha=\n",
@ -1614,23 +1608,13 @@
" |> iamax\n",
" in\n",
" \n",
" (* indice i du alpha max *)\n",
" let indice_ii = \n",
" (* indice i et j du alpha max *)\n",
" let indice_ii, indice_jj = \n",
" let max = max_element3 alpha_m \n",
" in\n",
" \n",
" (*Printf.printf \"%i\\n%!\" max;*)\n",
" \n",
" (max - 1) mod n_mo +1 \n",
" (max - 1) mod n_mo +1, (max - 1) / n_mo +1\n",
" in\n",
" \n",
" (* indice j du alpha max *)\n",
" let indice_jj = \n",
" let max = max_element3 alpha_m \n",
" in\n",
" (max - 1) / n_mo +1\n",
" in\n",
" \n",
" \n",
" (* Valeur du alpha max*)\n",
" let alpha alpha_m = \n",
" let i = indice_ii \n",
@ -1756,43 +1740,7 @@
"metadata": {},
"outputs": [],
"source": [
"(* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur \n",
"les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice\n",
"des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec \n",
"celle contenant i~ et j~ *)\n",
" \n",
"(* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n",
"(*\n",
"let f_Psi_tilde m_Ksi_tilde 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 m_Ksi_tilde.{i,1}\n",
" else if j=indice_j \n",
" then m_Ksi_tilde.{i,2}\n",
" else 0.)\n",
" \n",
";;\n",
"*)\n",
"(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
"(*\n",
"let f_Psi m_Ksi 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 m_Ksi.{i,1}\n",
" else if j=indice_j \n",
" then m_Ksi.{i,2}\n",
" else 0.)\n",
";;\n",
"*)\n",
"\n",
"(* 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",
@ -1805,10 +1753,8 @@
" then mat.{i,2}\n",
" else 0.)\n",
"(*************************)\n",
"\n",
"(*\n",
"let m_Psi = f_Psi m_Ksi indice_i indice_j;; \n",
"\n",
"let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j;; \n",
"*)"
]
@ -1824,12 +1770,9 @@
"let f_interm m_C m_Psi = Mat.sub m_C m_Psi;;\n",
"\n",
"(*************************)\n",
"\n",
"(*\n",
"let m_interm = f_interm m_C m_Psi;; \n",
"\n",
"let new_m_C m_C= Mat.add m_Psi_tilde m_interm;;\n",
"\n",
"let m_new_m_C = new_m_C m_C;;\n",
"*)"
]
@ -1893,8 +1836,8 @@
" (* Matrice des alphas *)\n",
" let m_alpha = alphad.m_alpha\n",
" in\n",
" (*let norme_alpha = norme m_alpha\n",
" in*)\n",
" let norme_alpha = norme m_alpha\n",
" in\n",
" \n",
" (*Printf.printf \"%f\\n%!\" norme_alpha;*)\n",
" \n",
@ -1957,12 +1900,12 @@
" (*Util.debug_matrix \"m_interm\" m_interm;*)\n",
" \n",
" (* Matrice après rotation *)\n",
" ( Mat.add m_Psi_tilde m_interm, critere_D)\n",
" ( Mat.add m_Psi_tilde m_interm, critere_D, norme_alpha)\n",
" in\n",
" let m_new_m_C , critere_D = new_m_C m_C methode loc_deloc \n",
" let m_new_m_C , critere_D, norme_alpha = new_m_C m_C methode loc_deloc \n",
" in\n",
" let diff = prev_critere_D -. critere_D +. 1.\n",
" \n",
" let diff = prev_critere_D -. critere_D \n",
"\n",
"in\n",
"\n",
"(*Util.debug_matrix \"new_alpha_m\" (f_alpha m_C);*)\n",
@ -1971,7 +1914,6 @@
"if diff**2. < cc**2.\n",
" then m_new_m_C\n",
" else\n",
"\n",
"localisation m_new_m_C methode loc_deloc epsilon (n-1) critere_D cc;;"
]
},
@ -2066,32 +2008,64 @@
"metadata": {},
"outputs": [],
"source": [
"let new_m_boys = localisation m_C \"boys\" \"loc\" 1. 100 0. 10e-7;;\n",
"let new_m_er = localisation m_C \"ER\" \"loc\" 1. 100 0. 10e-7;;\n",
"\n",
"let loc_m_occ_boys = localisation m_occ \"boys\" \"loc\" 1. 100 0. 10e-7;;\n",
"let loc_m_occ_er = localisation m_occ \"ER\" \"loc\" 1. 100 0. 10e-7;;\n",
"let loc_m_vir_boys = localisation m_vir \"boys\" \"loc\" 1. 100 0. 10e-7;;\n",
"let loc_m_vir_er = localisation m_vir \"ER\" \"loc\" 1. 100 0. 10e-7;;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let loc_m_vir_er = localisation m_vir \"ER\" \"loc\" 1. 100 0. 10e-7;;\n",
"\n",
"let m_assemble_boys = assemble loc_m_occ_boys loc_m_vir_boys;;\n",
"let m_assemble_er = assemble loc_m_occ_er loc_m_vir_er;;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let m_assemble_er = assemble loc_m_occ_er loc_m_vir_er;;\n",
"\n",
"let loc_assemble_boys = localisation m_assemble_boys \"boys\" \"loc\" 1. 100 0. 10e-7;;\n",
"let loc_assemble_er = localisation m_assemble_er \"er\" \"loc\" 1. 100 0. 10e-7;;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"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 hf = HartreeFock.make simulation\n",
"\n",
"let mo_basis = MOBasis.of_hartree_fock hf"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"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 "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,