Mdif, boyls, algo ER
This commit is contained in:
parent
4d456be51a
commit
2151c8927c
267
Work.ipynb
267
Work.ipynb
@ -161,6 +161,7 @@
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(*\n",
|
||||
"let xyz_string = \n",
|
||||
"\"3\n",
|
||||
"Water\n",
|
||||
@ -168,7 +169,28 @@
|
||||
"H -0.756950272703377558 0. -0.585882234512562827\n",
|
||||
"H 0.756950272703377558 0. -0.585882234512562827\n",
|
||||
"\"\n",
|
||||
" "
|
||||
" *) "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### $H_2$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"let xyz_string = \n",
|
||||
"\"2\n",
|
||||
"H2\n",
|
||||
"H 0. 0. 0.\n",
|
||||
"H 1.8 0. 0.\n",
|
||||
"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -199,6 +221,19 @@
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"let basis_string = \"\n",
|
||||
"HYDROGEN\n",
|
||||
"S 4\n",
|
||||
"1 1.301000E+01 1.968500E-02\n",
|
||||
"2 1.962000E+00 1.379770E-01\n",
|
||||
"3 4.446000E-01 4.781480E-01\n",
|
||||
"4 1.220000E-01 5.012400E-01\n",
|
||||
"S 1\n",
|
||||
"1 1.220000E-01 1.000000E+00\n",
|
||||
"P 1\n",
|
||||
"1 7.270000E-01 1.0000000\n",
|
||||
"\"\n",
|
||||
"(*\n",
|
||||
"let basis_string = \"\n",
|
||||
"HYDROGEN\n",
|
||||
"S 4\n",
|
||||
@ -244,7 +279,8 @@
|
||||
"D 1\n",
|
||||
"1 1.185000E+00 1.0000000\n",
|
||||
"\n",
|
||||
"\""
|
||||
"\"\n",
|
||||
"*)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -950,19 +986,41 @@
|
||||
"\n",
|
||||
"$D(\\theta) = A_{ii} + A_{jj} - (cos^2\\theta B_{ii} + sin^2\\theta B_{jj} - sin2\\theta B_{ij})^2 - (sin^2\\theta B_{ii} + cos^2\\theta B_{jj} + sin2\\theta B_{ij})^2$\n",
|
||||
"\n",
|
||||
"Le but est de minimiser la matrice $D(\\theta)$ il faut donc faire une rotation des orbitales i et j donnants le plus grand élément de la matrice $D(\\theta)$\n",
|
||||
"\n",
|
||||
"Si je comprends bien, on doit prendre le plus grand élément $D(i,j)$ de la matrice $D(\\theta)$ pour $\\theta = 0$ au départ et d'appliquer des rotation en faisant varier $\\theta$ pour minimiser $D(i,j)$ sachant que les extrema de $D(\\theta)$ sont de la forme :\n",
|
||||
"Les extrema de $D(\\theta)$ sont de la forme :\n",
|
||||
"\n",
|
||||
"$Tan(\\theta) = - \\frac{\\gamma}{\\beta}$\n",
|
||||
"\n",
|
||||
"$D(\\theta]$ à 4 extrema :\n",
|
||||
"$D(\\theta)$ à 4 extrema :\n",
|
||||
"\n",
|
||||
"$4\\theta; \\;\\; 4\\theta +\\pi; \\;\\; 4\\theta+ 2\\pi; \\;\\; 4\\theta+ 3\\pi$\n",
|
||||
"\n",
|
||||
"Puis recommencer avec le nouveau $D(i,j)max$\n",
|
||||
"$\\theta_{ij}= Atan(- \\frac{\\gamma}{\\beta})$\n",
|
||||
"\n",
|
||||
"\n"
|
||||
"Ainsi on peut calculer la matrice des $\\theta$ et effectuer la rotation pour le couple d'orbitale $(i,j)$ ayant le $\\theta_{max}$\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Sinon on procède comme Edminson Ruedenberg mais avec les intégrales $A_{12}$ et $B_{12}$ définies comme :\n",
|
||||
"\n",
|
||||
"$A^r_{12} = \\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n",
|
||||
"$*\\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n",
|
||||
"$- \\frac {1}{4}(\\langle \\phi_1 | \\bar r | \\phi_1 \\rangle $\n",
|
||||
"$- \\langle \\phi_2 | \\bar r | \\phi_2 \\rangle * \\langle \\phi_1 | \\bar r | \\phi_1 \\rangle$\n",
|
||||
"$- \\langle \\phi_2 | \\bar r | \\phi_2 \\rangle)$\n",
|
||||
"\n",
|
||||
"Et \n",
|
||||
"\n",
|
||||
"$B^r_{12} = (\\langle \\phi_1 | \\bar r | \\phi_1 \\rangle - \\langle \\phi_2 | \\bar r | \\phi_2 \\rangle)$\n",
|
||||
"$ * \\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n",
|
||||
"\n",
|
||||
"Avec \n",
|
||||
"\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}$"
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -971,7 +1029,14 @@
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(*let multipoles = \n",
|
||||
"(* Test Boys *)\n",
|
||||
"\n",
|
||||
"let n_ao = Mat.dim1 m_C ;;\n",
|
||||
"let n_mo = n_ao - 1;;\n",
|
||||
"\n",
|
||||
"let m_alpha_boys multipoles = \n",
|
||||
"\n",
|
||||
" let multipoles = \n",
|
||||
" AOBasis.multipole ao_basis\n",
|
||||
" in\n",
|
||||
"\n",
|
||||
@ -979,12 +1044,65 @@
|
||||
" Multipole.matrix_x multipoles \n",
|
||||
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef \n",
|
||||
" in \n",
|
||||
" \n",
|
||||
" let phi_y_phi =\n",
|
||||
" Multipole.matrix_y multipoles \n",
|
||||
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef \n",
|
||||
" in\n",
|
||||
"*)"
|
||||
" let phi_z_phi =\n",
|
||||
" Multipole.matrix_y multipoles \n",
|
||||
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef\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} \n",
|
||||
" -. (1. /. 4.) *. (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",
|
||||
" asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )));;\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"(* let m_a12 =\n",
|
||||
" \n",
|
||||
" let m_a12_x = Mat.init_cols n_mo n_mo ( fun i j -> \n",
|
||||
" phi_x_phi.{i,j} *. phi_x_phi.{i,j} \n",
|
||||
" -. (1. /. 4.)(phi_x_phi.{i,i} -. phi_x_phi.{j,j} *. phi_x_phi.{i,i} -. phi_x_phi.{j,j}))\n",
|
||||
" in\n",
|
||||
" let m_a12_y = Mat.init_cols n_mo n_mo ( fun i j -> \n",
|
||||
" phi_y_phi.{i,j} *. phi_y_phi.{i,j} \n",
|
||||
" -. (1. /. 4.)(phi_y_phi.{i,i} -. phi_y_phi.{j,j} *. phi_y_phi.{i,i} -. phi_y_phi.{j,j}))\n",
|
||||
" in\n",
|
||||
" let m_a12_z = Mat.init_cols n_mo n_mo ( fun i j -> \n",
|
||||
" phi_z_phi.{i,j} *. phi_z_phi.{i,j} \n",
|
||||
" -. (1. /. 4.)(phi_z_phi.{i,i} -. phi_z_phi.{j,j} *. phi_z_phi.{i,i} -. phi_z_phi.{j,j}))\n",
|
||||
" in \n",
|
||||
" Mat.add m_a12_x Mat.add (m_a12_y m_a12_z)\n",
|
||||
"\n",
|
||||
" let m_b12 = \n",
|
||||
"\n",
|
||||
" let m_b12_x = Mat.init_cols n_mo n_mo ( fun i j ->\n",
|
||||
" (phi_x_phi.{i,i} -. phi_x_phi.{j,j}) *. phi_x_phi.{i,j})\n",
|
||||
" in \n",
|
||||
" let m_b12_y = Mat.init_cols n_mo n_mo ( fun i j ->\n",
|
||||
" (phi_y_phi.{i,i} -. phi_y_phi.{j,j}) *. phi_y_phi.{i,j})\n",
|
||||
" in \n",
|
||||
" let m_b12_z = Mat.init_cols n_mo n_mo ( fun i j ->\n",
|
||||
" (phi_z_phi.{i,i} -. phi_z_phi.{j,j}) *. phi_z_phi.{i,j})\n",
|
||||
" in \n",
|
||||
" Mat.add m_b12_x Mat.add (m_b12_y m_b12_z)\n",
|
||||
"\n",
|
||||
"in\n",
|
||||
"Mat.init_cols n_mo n_mo ( fun i j ->\n",
|
||||
" asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )));;*)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -1202,11 +1320,7 @@
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": true
|
||||
}
|
||||
},
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(* Localisation de Edminstion *)\n",
|
||||
@ -1240,11 +1354,7 @@
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": true
|
||||
}
|
||||
},
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(* Calcul de toutes les intégrales B_12 -> Matrice *)\n",
|
||||
@ -1272,11 +1382,7 @@
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": true
|
||||
}
|
||||
},
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(* Extraction du plus grand angle alpha*)\n",
|
||||
@ -1341,18 +1447,18 @@
|
||||
"celle contenant i~ et j~ *)\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_mo (fun i j -> if j=q then m_Ksi_tilde.{i,1}\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",
|
||||
" else if j=p then m_Ksi_tilde.{i,2}\n",
|
||||
" else 0.);;\n",
|
||||
" (* p q verif*) \n",
|
||||
"(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
|
||||
"let m_Psi = Mat.init_cols n_ao n_mo (fun i j -> if j=q then m_Ksi.{i,1}\n",
|
||||
" else if j=p then m_Ksi.{i,2}\n",
|
||||
"let m_Psi = Mat.init_cols n_ao n_ao (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",
|
||||
"(* 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 m_interm = Mat.sub m_C m_Psi;;\n",
|
||||
"let m_interm = Mat.sub m_Phi m_Psi;;\n",
|
||||
"\n",
|
||||
"(* Construction de la nouvelle matrice d'OMs phi~ *)\n",
|
||||
"let m_Phi_tilde = Mat.add m_Psi_tilde m_interm;; "
|
||||
@ -1421,29 +1527,30 @@
|
||||
" loc : float ;\n",
|
||||
"}\n",
|
||||
"\n",
|
||||
"let prev_iter = {coef = m_Phi ;loc = 0. }\n",
|
||||
"\n",
|
||||
"(* Calcul de la nouvelle matrice Phi après rotations *)\n",
|
||||
"let rec new_Phi prev_iter n = \n",
|
||||
"\n",
|
||||
" Printf.printf \"%d\\n%!\" n;\n",
|
||||
" if n == 0 then\n",
|
||||
" m_Phi\n",
|
||||
" prev_iter\n",
|
||||
" else\n",
|
||||
"\n",
|
||||
" (* Calcul de tous les alpha -> Matrice *)\n",
|
||||
" let m_alpha m_C =\n",
|
||||
" let m_alpha m_Phi =\n",
|
||||
"\n",
|
||||
" (* Calcul de toutes les intégrales B_12 -> Matrice *)\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_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,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_mo n_mo (fun i j ->\n",
|
||||
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,j} \n",
|
||||
" -. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} ) \n",
|
||||
" *. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,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",
|
||||
@ -1474,7 +1581,7 @@
|
||||
" in\n",
|
||||
" Printf.printf \"%f\\n%!\" alpha;\n",
|
||||
"\n",
|
||||
" in\n",
|
||||
" \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",
|
||||
@ -1503,13 +1610,13 @@
|
||||
" let m_interm = \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_mo (fun i j -> if j=q then m_Ksi_tilde.{i,1}\n",
|
||||
" else if j=p then m_Ksi_tilde.{i,2}\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",
|
||||
" else if j=p then m_Ksi_tilde.{i,1}\n",
|
||||
" else 0.)\n",
|
||||
" in\n",
|
||||
" (* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
|
||||
" let m_Psi = Mat.init_cols n_ao n_mo (fun i j -> if j=q then m_Ksi.{i,1}\n",
|
||||
" else if j=p then m_Ksi.{i,2}\n",
|
||||
" let m_Psi = Mat.init_cols n_ao n_ao (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",
|
||||
" in\n",
|
||||
@ -1517,22 +1624,55 @@
|
||||
" in\n",
|
||||
" \n",
|
||||
" let coef = Mat.add m_Psi_tilde m_interm in\n",
|
||||
" let loc = .... in\n",
|
||||
" let loc = alpha \n",
|
||||
" in\n",
|
||||
" let result = \n",
|
||||
" { coef ; loc }\n",
|
||||
" in\n",
|
||||
" if loc - prev_iter.loc < 1.e-5 then\n",
|
||||
" if loc -. prev_iter.loc < 1.e-5 then\n",
|
||||
" result\n",
|
||||
" else\n",
|
||||
" new_Phi result (n-1)\n",
|
||||
" let result = prev_iter\n",
|
||||
" in\n",
|
||||
"(*-----*)\n",
|
||||
"new_Phi m_C 10;;\n",
|
||||
" new_Phi result (n-1)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
" new_Phi prev_iter 10;;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(*Test critère de convergence -> norme de alpha *)\n",
|
||||
"\n",
|
||||
"let norme m_alpha =\n",
|
||||
" let vec_alpha = \n",
|
||||
" Mat.as_vec m_alpha\n",
|
||||
" in\n",
|
||||
" let vec2_alpha = \n",
|
||||
" Vec.sqr vec_alpha\n",
|
||||
" in\n",
|
||||
" let norme2 = Vec.sum vec2_alpha\n",
|
||||
"in\n",
|
||||
"sqrt(norme2);;\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"norme m_alpha;;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
@ -1540,6 +1680,39 @@
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(* Test choix méthode *)\n",
|
||||
"type localisation = \n",
|
||||
" | Boys\n",
|
||||
" | Edminston \n",
|
||||
"\n",
|
||||
"let toto localisation = \n",
|
||||
" match localisation with \n",
|
||||
" | Boys -> 1. +. 1.\n",
|
||||
" | Edminston -> 1. +. 2.;;\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"toto Boys ;;\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"let titi lo=\n",
|
||||
"if lo=\"Boys\" then 1. else 2.;;\n",
|
||||
"titi \"Boys\";;\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"let alpha localisation toto = \n",
|
||||
"if localisation = 1 \n",
|
||||
"then \n",
|
||||
"\n",
|
||||
"else \n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
|
||||
Loading…
Reference in New Issue
Block a user