Work : calcul D_boys et match D

This commit is contained in:
Yann Damour 2020-05-08 10:49:21 +02:00
parent 9e39302fcf
commit 7b24683ea3

View File

@ -203,24 +203,24 @@
"metadata": {},
"outputs": [],
"source": [
"(* Fonction création chaîne de n H *)\n",
"(* Fonction création chaîne linéaire de n H *)\n",
"\n",
"let xyz n = \n",
"let xyz d n = \n",
" let accu = \"\"\n",
" in\n",
" let rec toto accu n =\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\" ^ \" \" ^ string_of_float(1.8 *. float_of_int(n)) ^ \" 0. 0.\\n\" ^ accu\n",
" | x -> toto (\"H\" ^ \" \" ^ string_of_float(1.8 *. float_of_int(n)) ^ \" 0. 0.\\n\" ^ accu ) (n-1)\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 n;;\n",
" in string_of_int(n) ^ \"\\nH\" ^ string_of_int(n) ^ \"\\n\" ^ toto accu d n;;\n",
" \n",
"let xyz_string = xyz 4;;"
"let xyz_string = xyz 1.8 4;;"
]
},
{
@ -1054,7 +1054,15 @@
"\n",
"$A^r_{12}=A^x_{12} + A^y_{12} + A^z_{12}$\n",
"\n",
"$B^r_{12}=B^x_{12} + B^y_{12} + B^z_{12}$"
"$B^r_{12}=B^x_{12} + B^y_{12} + B^z_{12}$\n",
"\n",
"Et le critère à minimiser est :\n",
"\n",
"$D= \\sum_i < \\phi_i | r | \\phi_i >$\n",
"\n",
"Avec \n",
"\n",
"$< \\phi_i | r | \\phi_i > = < \\phi_i | x | \\phi_i > + < \\phi_i | y | \\phi_i > + < \\phi_i | z | \\phi_i >$\n"
]
},
{
@ -1247,7 +1255,8 @@
"\n",
"\n",
"Avec : \n",
" $A_{ij} = [ \\phi_i \\phi_j | \\phi_i \\phi_j] - 1/4 [\\phi_i^2 - \\phi_j^2 | \\phi_i^2 - \\phi_j^2] $\n",
"\n",
"$A_{ij} = [ \\phi_i \\phi_j | \\phi_i \\phi_j] - 1/4 [\\phi_i^2 - \\phi_j^2 | \\phi_i^2 - \\phi_j^2] $\n",
" \n",
"Où :\n",
"\n",
@ -1627,8 +1636,9 @@
" then 0. \n",
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))));;\n",
"\n",
"f_alpha m_C;; \n",
" "
"(*********************)\n",
"\n",
"f_alpha m_C;; "
]
},
{
@ -1644,11 +1654,12 @@
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
" ) i j\n",
" )\n",
" in Vec.init n_ao ( fun i -> m_D.{i,1} )\n",
" in Vec.init n_ao ( fun i -> m_D.{i,i} )\n",
" in Vec.sum v_D ;;\n",
"\n",
"(******************)\n",
"\n",
"s_D m_C;;\n"
"s_D m_C;;"
]
},
{
@ -1669,7 +1680,7 @@
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef \n",
" in\n",
" let phi_z_phi =\n",
" Multipole.matrix_y multipoles \n",
" Multipole.matrix_z multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef\n",
" in \n",
"\n",
@ -1689,10 +1700,54 @@
" Mat.init_cols n_ao n_ao ( fun i j ->\n",
" 0.25 *. asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )));;\n",
"\n",
"(****************************)\n",
"\n",
"f_alpha_boys multipoles;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Fonction de calcul de D Boys *)\n",
"\n",
"let d_boys multipoles = \n",
"\n",
" let phi_x2_phi =\n",
" Multipole.matrix_x2 multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef \n",
" in\n",
" \n",
" (*Util.debug_matrix \"phi_x2_phi\" phi_x2_phi;*)\n",
" \n",
" let phi_y2_phi =\n",
" Multipole.matrix_y2 multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef \n",
" in\n",
" \n",
" (*Util.debug_matrix \"phi_y2_phi\" phi_y2_phi;*)\n",
" \n",
" let phi_z2_phi =\n",
" Multipole.matrix_z2 multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef\n",
"\n",
" in\n",
" \n",
" (*Util.debug_matrix \"phi_z2_phi\" phi_z2_phi;*)\n",
" \n",
" let v_D_boys = \n",
" Vec.init n_ao ( fun i -> phi_x2_phi.{i,i} +. phi_y2_phi.{i,i} +. phi_z2_phi.{i,i})\n",
"\n",
"in\n",
"Vec.sum v_D_boys;;\n",
"\n",
"(*************************)\n",
"\n",
"d_boys multipoles;;"
]
},
{
"cell_type": "code",
"execution_count": null,
@ -1717,6 +1772,8 @@
"in \n",
"alpha methode;;\n",
"\n",
"(*************************)\n",
"\n",
"calcul_alpha_2 \"ER\" ;;\n",
"calcul_alpha_2 \"Boys\" ;;\n",
"let methode = \"ER\";;\n",
@ -1732,7 +1789,7 @@
"(* Test méthode de calcul du critère de convergence *)\n",
"let calcul_cc methode = \n",
"\n",
"let b_boys = 2.\n",
"let b_boys = d_boys multipoles\n",
"in\n",
"let d_er = s_D m_C\n",
"in\n",
@ -1747,7 +1804,9 @@
"in \n",
"cc methode;;\n",
"\n",
"calcul_cc \"ER\";;\n"
"(*************************)\n",
"\n",
"calcul_cc \"boys\";;"
]
},
{
@ -1821,6 +1880,8 @@
" then {alpha_max; indice_ii; indice_jj}\n",
" else new_m_alpha alpha_m (n_rec_alpha-1);;\n",
"\n",
"(*************************)\n",
"\n",
"(*\n",
"let m_alpha = f_alpha m_C;;\n",
"new_m_alpha m_alpha 3;;\n",
@ -1842,6 +1903,8 @@
" then sin alpha \n",
" else -. sin alpha );;\n",
"\n",
"(*************************)\n",
"\n",
"(*\n",
"let alpha = alphaij.alpha_max;; (* Fonction -> constante *) \n",
"f_R alpha;;\n",
@ -1876,6 +1939,8 @@
"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 = 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",
"(*************************)\n",
"\n",
"(*\n",
"let m_Ksi = f_Ksi indice_i indice_j m_C;; (* Fonction -> constante *)\n",
"*)"
@ -1890,6 +1955,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 = gemm m_Ksi m_R;;\n",
"\n",
"(*************************)\n",
"\n",
"(*\n",
"let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi;; (* Fonction -> constante *)\n",
"*)"
@ -1922,6 +1990,9 @@
" then m_Ksi.{i,2}\n",
" else 0.)\n",
";;\n",
"\n",
"(*************************)\n",
"\n",
"(*\n",
"let m_Psi = f_Psi m_Ksi indice_i indice_j;; (* Fonction -> constante *)\n",
"\n",
@ -1938,6 +2009,9 @@
"(* 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\n",
"par la matrice *)\n",
"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;; (* Fonction -> constante *)\n",
"\n",