This commit is contained in:
Anthony Scemama 2020-05-04 16:14:06 +02:00
parent 1d473819c9
commit 759ea51b6e
2 changed files with 75 additions and 16 deletions

View File

@ -923,6 +923,60 @@
"ERI.get_phys ints 1 2 1 2;;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"type localisation = \n",
"| Boys \n",
"| Edminston\n",
"\n",
"let alpha_boys ao_basis m_C =\n",
" ao_basis + m_C\n",
" \n",
"let alpha_er ao_basis m_C =\n",
" ao_basis + m_C\n",
"\n",
"let make_alpha loc =\n",
" match loc with\n",
" | Boys -> alpha_boys \n",
" | Edminston -> alpha_er \n",
" \n",
" \n",
"let l_method_of_string s =\n",
" match s with\n",
" | \"Boys\"\n",
" | \"boys\" -> Boys\n",
" | \"ER\"\n",
" | \"Er\"\n",
" | \"er\" -> Edminston\n",
" | _ -> invalid_arg \"message d'erreur\"\n",
"\n",
"\n",
"(*\n",
"let localize m =\n",
" let alpha = make_alpha m\n",
";; \n",
"\n",
"localize (l_method_of_string \"Boys\")\n",
";;\n",
"*)\n",
"\n",
"\n",
"let f c a b = \n",
" let add a b = a+b in\n",
" let sub a b = a-b in\n",
" let g =\n",
" match c with \n",
" | Boys -> add \n",
" | _ -> sub \n",
" in\n",
" g a b\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,

View File

@ -78,7 +78,7 @@
"(* --------- *)\n",
"\n",
"(*Mettre le bon chemin ici *)\n",
"#cd \"/home/ydamour/QCaml\";;\n",
"#cd \"/home/scemama/QCaml\";;\n",
"\n",
"#use \"topfind\";;\n",
"#require \"jupyter.notebook\";;\n",
@ -615,7 +615,7 @@
"metadata": {},
"outputs": [],
"source": [
"let nocc = 5 (* On a 10 electrons, donc 5 orbitales occupees. *)\n",
"let nocc = 2 (* On a 10 electrons, donc 5 orbitales occupees. *)\n",
"\n",
"\n",
"let c_of_h m_H = \n",
@ -884,7 +884,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Localisation de Boyls"
"## Localisation de Boys"
]
},
{
@ -1563,7 +1563,7 @@
"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 = n_ao - 1;;\n",
"let n_mo = Mat.dim2 m_C ;;\n",
"let sum a = \n",
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
";;\n",
@ -1586,15 +1586,13 @@
") (Util.array_range 1 n_ao)\n",
"|> sum \n",
"\n",
"\n",
"let m_Phi = m_C;;\n",
"\n",
"let m_Phi = None;;\n",
"type iteration = \n",
"{ coef : Mat.t ;\n",
" loc : float ;\n",
"}\n",
"\n",
"let prev_iter = {coef = m_Phi ;loc = 0. }\n",
"let prev_iter = {coef = m_C ;loc = 0. }\n",
"\n",
"(* Calcul de la nouvelle matrice Phi après rotations *)\n",
"let rec new_Phi prev_iter n = \n",
@ -1608,28 +1606,28 @@
" let f_alpha m_Phi =\n",
"\n",
" (* Calcul de toutes les intégrales B_12 -> Matrice *)\n",
" let m_b12 = Mat.init_cols n_ao n_ao (fun i j -> \n",
" let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> \n",
" integral_general (fun a b e f i j ->\n",
" ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} ) *. m_Phi.{e,i} *. m_Phi.{f,j}\n",
" ) i j\n",
" )\n",
" in\n",
" (* Calcul de toutes les intégrales A_12 -> Matrice *)\n",
" let m_a12 = Mat.init_cols n_ao n_ao (fun i j ->\n",
" let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->\n",
" integral_general (fun a b e f i j -> m_Phi.{a,i} *. m_Phi.{b,i} *. m_Phi.{e,i} *. m_Phi.{f,j} \n",
" -. 0.25 *. ( m_Phi.{e,i} *. m_Phi.{f,i} -. m_Phi.{e,j} *. m_Phi.{f,j} ) \n",
" *. ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} )\n",
" ) i j\n",
" )\n",
" in\n",
" Mat.init_cols n_ao n_ao ( fun i j ->\n",
" Mat.init_cols n_ao n_mo ( fun i j ->\n",
" asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n",
" )\n",
" )\n",
"\n",
" in\n",
" (* Matrice de rotation 2 par 2 *)\n",
" let m_R = \n",
" let indice_i, indice_j, m_R = \n",
"\n",
" let m_alpha = f_alpha m_Phi\n",
" in\n",
@ -1651,14 +1649,17 @@
" in\n",
" Printf.printf \"%f\\n%!\" alpha;\n",
" \n",
" indice_i, indice_j, \n",
" Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha\n",
" else if i>j then sin alpha \n",
" else -. sin alpha )\n",
"\n",
" in\n",
" Util.debug_matrix \"m_R\" m_R;\n",
" Printf.printf \"%d %d\\n%!\" indice_i indice_j;\n",
" (* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle \n",
" on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n",
" let m_Ksi_tilde = \n",
" let m_Ksi, m_Ksi_tilde = \n",
"\n",
" (* Fonction d'extraction des 2 vecteurs propres i et j de la matrice des OMs pour les mettres dans la matrice Ksi \n",
" (n par 2) pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n",
@ -1668,15 +1669,18 @@
" in\n",
" let m_Ksi = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_Phi.{i,p} else m_Phi.{i,q} )\n",
" in\n",
" gemm m_Ksi m_R\n",
" m_Ksi, gemm m_Ksi m_R\n",
"\n",
" in\n",
" \n",
" Util.debug_matrix \"m_Ksi\" m_Ksi;\n",
" Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n",
" (* Construction de la nouvelle matrice d'OMs phi~ *)\n",
" let m_Phi_tilde m_Ksi m_Ksi_tilde = \n",
"\n",
" \n",
" (* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n",
" let m_Psi_tilde = Mat.init_cols n_ao n_ao (fun i j -> if j=q then m_Ksi_tilde.{i,2}\n",
" let m_Psi_tilde = Mat.init_cols n_ao n_mo (fun i j -> if j=q then m_Ksi_tilde.{i,2}\n",
" else if j=p then m_Ksi_tilde.{i,1}\n",
" else 0.)\n",
" in\n",
@ -1685,7 +1689,7 @@
" let m_interm = \n",
"\n",
" (* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
" let m_Psi = Mat.init_cols n_ao n_ao (fun i j -> if j=q then m_Ksi.{i,2}\n",
" let m_Psi = Mat.init_cols n_ao n_mo (fun i j -> if j=q then m_Ksi.{i,2}\n",
" else if j=p then m_Ksi.{i,1}\n",
" else 0.)\n",
"\n",
@ -1696,6 +1700,7 @@
" in \n",
" let coef = m_Phi_tilde m_Ksi m_Ksi_tilde\n",
" in\n",
" Util.debug_matrix \"coef\" coef;\n",
" let loc = alpha \n",
" in\n",
" let result = \n",