Work: problème index boucle
This commit is contained in:
parent
759ea51b6e
commit
0b36d1ea8c
511
Work.ipynb
511
Work.ipynb
@ -78,7 +78,7 @@
|
||||
"(* --------- *)\n",
|
||||
"\n",
|
||||
"(*Mettre le bon chemin ici *)\n",
|
||||
"#cd \"/home/scemama/QCaml\";;\n",
|
||||
"#cd \"/home/ydamour/QCaml\";;\n",
|
||||
"\n",
|
||||
"#use \"topfind\";;\n",
|
||||
"#require \"jupyter.notebook\";;\n",
|
||||
@ -615,7 +615,7 @@
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"let nocc = 2 (* On a 10 electrons, donc 5 orbitales occupees. *)\n",
|
||||
"let nocc = 1 (* On a 10 electrons, donc 5 orbitales occupees. *)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"let c_of_h m_H = \n",
|
||||
@ -1029,6 +1029,7 @@
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(*\n",
|
||||
"(* Test Boys *)\n",
|
||||
"\n",
|
||||
"let n_ao = Mat.dim1 m_C ;;\n",
|
||||
@ -1070,7 +1071,8 @@
|
||||
" asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )));;\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"m_alpha_boys multipoles;;\n"
|
||||
"m_alpha_boys multipoles;;\n",
|
||||
"*)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -1244,6 +1246,24 @@
|
||||
"\n",
|
||||
"$= \\left(\\sum_e c_{ei} \\sum_f c_{fi} - \\sum_e c_{ej} \\sum_f c_{fj} \\right) \\left(\\sum_a c_{ai} \\sum_b c_{bi} - \\sum_a c_{aj} \\sum_b c_{bj} \\right) [ \\chi_a \\chi_b| \\chi_e \\chi_f ] $\n",
|
||||
"\n",
|
||||
"(\n",
|
||||
"Ainsi comme :\n",
|
||||
"\n",
|
||||
"$\\phi_i ^2 = \\left( \\sum_a c_{ai} \\chi_a \\right)^2$\n",
|
||||
"$= \\sum_a c_{ai} \\sum_b c_{bi} \\chi_a \\chi_b$\n",
|
||||
"\n",
|
||||
"On aura :\n",
|
||||
"\n",
|
||||
"$D=\\sum_n [\\phi_n^2|\\phi_n^2]$\n",
|
||||
"\n",
|
||||
"$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} [\\chi_a \\chi_b| \\phi_n^2])$\n",
|
||||
"\n",
|
||||
"$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} \\sum_e c_{en} \\sum_f c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$\n",
|
||||
"\n",
|
||||
"$= \\sum_n (\\sum_a \\sum_b \\sum_e \\sum_f c_{an} c_{bn} c_{en} c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$\n",
|
||||
"\n",
|
||||
")\n",
|
||||
"\n",
|
||||
"Mais aussi :\n",
|
||||
"\n",
|
||||
"$B_{ij} = [\\phi_i ^2 - \\phi_j ^2 | \\phi_i \\phi_j ] = [\\left( \\sum_a c_{ai} \\chi_i \\right)^2 - \\left( \\sum_a c_{aj} \\chi_j \\right)^2| \\phi_i \\phi_j ] $\n",
|
||||
@ -1285,86 +1305,6 @@
|
||||
"## Calculs : Localisation de Edminston"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(* Localisation de Edminstion *)\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 = n_ao - 1;;\n",
|
||||
"let sum a = \n",
|
||||
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
|
||||
";;\n",
|
||||
"\n",
|
||||
"(*Calcul des intégrales*) \n",
|
||||
"let integral_general g i j =\n",
|
||||
"Array.map (fun a ->\n",
|
||||
" let v = \n",
|
||||
" Array.map (fun b ->\n",
|
||||
" let u = \n",
|
||||
" Array.map (fun e ->\n",
|
||||
" let t = Array.map (fun f ->\n",
|
||||
" (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
|
||||
" ) (Util.array_range 1 n_ao)\n",
|
||||
" in sum t\n",
|
||||
" ) (Util.array_range 1 n_ao)\n",
|
||||
" in sum u\n",
|
||||
" ) (Util.array_range 1 n_ao)\n",
|
||||
" in sum v\n",
|
||||
") (Util.array_range 1 n_ao)\n",
|
||||
"|> sum"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(* Calcul de toutes les intégrales B_12 -> Matrice *)\n",
|
||||
"let m_b12 = Mat.init_cols n_ao n_ao (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",
|
||||
" ) i j\n",
|
||||
" )\n",
|
||||
"\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",
|
||||
"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",
|
||||
" ) i j\n",
|
||||
" )\n",
|
||||
"\n",
|
||||
"(* Calcul de tous les alpha -> Matrice *)\n",
|
||||
"let m_alpha = Mat.init_cols n_ao n_ao ( fun i j ->\n",
|
||||
" asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n",
|
||||
" )\n",
|
||||
" )"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"(* Extraction du plus grand angle alpha*)\n",
|
||||
"let ij_max_alpha = Mat.as_vec m_alpha\n",
|
||||
"|> iamax;;\n",
|
||||
"\n",
|
||||
"(* Indices du plus grand angle alpha *)\n",
|
||||
"let indice_i = ij_max_alpha mod n_ao;;\n",
|
||||
"let indice_j = ij_max_alpha / n_ao + 1;;\n",
|
||||
"\n",
|
||||
"(* Valeur du plus grand angle alpha*)\n",
|
||||
"let alpha = m_alpha.{indice_i,indice_j};;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
@ -1468,7 +1408,7 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Rotation des orbitales"
|
||||
"### ER-V1"
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -1477,6 +1417,100 @@
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"(* Localisation de Edminstion *)\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 sum a = \n",
|
||||
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
|
||||
";;\n",
|
||||
"\n",
|
||||
"(*Calcul des intégrales*) \n",
|
||||
"let integral_general g i j =\n",
|
||||
"Array.map (fun a ->\n",
|
||||
" let v = \n",
|
||||
" Array.map (fun b ->\n",
|
||||
" let u = \n",
|
||||
" Array.map (fun e ->\n",
|
||||
" let t = Array.map (fun f ->\n",
|
||||
" (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
|
||||
" ) (Util.array_range 1 n_ao)\n",
|
||||
" in sum t\n",
|
||||
" ) (Util.array_range 1 n_ao)\n",
|
||||
" in sum u\n",
|
||||
" ) (Util.array_range 1 n_ao)\n",
|
||||
" in sum v\n",
|
||||
") (Util.array_range 1 n_ao)\n",
|
||||
"|> sum\n",
|
||||
";;\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",
|
||||
"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",
|
||||
" ) i j\n",
|
||||
" )\n",
|
||||
"\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,j} *. 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",
|
||||
" ) i j\n",
|
||||
" )\n",
|
||||
"\n",
|
||||
"(* Calcul de tous les alpha -> Matrice *)\n",
|
||||
"let m_alpha = 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",
|
||||
"\n",
|
||||
"let s_D = \n",
|
||||
" let v_D = \n",
|
||||
" let m_D = 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,i} \n",
|
||||
" ) i j\n",
|
||||
" )\n",
|
||||
" in Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
|
||||
"in Vec.sum v_D \n",
|
||||
"\n",
|
||||
"(*\n",
|
||||
"(* Calcul de D *)\n",
|
||||
"let m_D = 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,i} \n",
|
||||
" ) i j\n",
|
||||
" )\n",
|
||||
" \n",
|
||||
"let v_D = Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
|
||||
"\n",
|
||||
"let s_D = Vec.sum v_D\n",
|
||||
"*)\n",
|
||||
"\n",
|
||||
" (* \n",
|
||||
"(* Calcul de D *)\n",
|
||||
"let v_D = Vec.init n_ao (fun i ->\n",
|
||||
"integral_general (fun a b e f i -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
|
||||
" ) i \n",
|
||||
" )\n",
|
||||
"\n",
|
||||
"*)\n",
|
||||
"(* Rotation des orbitales *)\n",
|
||||
"(*\n",
|
||||
"(* Extraction du plus grand angle alpha*)\n",
|
||||
"let ij_max_alpha = Mat.as_vec m_alpha\n",
|
||||
"|> iamax;;\n",
|
||||
"\n",
|
||||
"(* Indices du plus grand angle alpha *)\n",
|
||||
"let indice_i = ij_max_alpha mod n_ao;;\n",
|
||||
"let indice_j = ij_max_alpha / n_ao + 1;;\n",
|
||||
"\n",
|
||||
"(* Valeur du plus grand angle alpha*)\n",
|
||||
"let alpha = m_alpha.{indice_i,indice_j};;\n",
|
||||
"*)\n",
|
||||
"\n",
|
||||
"(*\n",
|
||||
"let m_Phi = m_C\n",
|
||||
"\n",
|
||||
"(* Matrice de rotation 2 par 2 *)\n",
|
||||
@ -1492,22 +1526,12 @@
|
||||
"\n",
|
||||
"(* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle on\n",
|
||||
"obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n",
|
||||
"let m_Ksi_tilde = gemm m_Ksi m_R;;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Réinjection"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"let m_Ksi_tilde = gemm m_Ksi m_R;;\n",
|
||||
"*)\n",
|
||||
"\n",
|
||||
"(* Réinjection*)\n",
|
||||
"\n",
|
||||
"(*\n",
|
||||
"(* 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",
|
||||
@ -1528,7 +1552,8 @@
|
||||
"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;; "
|
||||
"let m_Phi_tilde = Mat.add m_Psi_tilde m_interm;; \n",
|
||||
"*)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -1538,13 +1563,6 @@
|
||||
"# TESTS"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
@ -1558,6 +1576,7 @@
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"(* Fonction pour itérer , cc critère de convergence*) \n",
|
||||
"(* Localisation de Edminstion *)\n",
|
||||
"let ee_ints = AOBasis.ee_ints ao_basis;;\n",
|
||||
@ -1592,16 +1611,23 @@
|
||||
" loc : float ;\n",
|
||||
"}\n",
|
||||
"\n",
|
||||
"let prev_iter = {coef = m_C ;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",
|
||||
" Printf.printf \"%d\\n%!\" n;\n",
|
||||
" (*\n",
|
||||
" Util.debug_matrix \"prev_iter\" prev_iter.coef;\n",
|
||||
" *)\n",
|
||||
" if n == 0 then\n",
|
||||
" prev_iter\n",
|
||||
" else\n",
|
||||
" let m_Phi = prev_iter.coef\n",
|
||||
" in\n",
|
||||
" (*\n",
|
||||
" Util.debug_matrix \"m_coef\" prev_iter.coef;\n",
|
||||
" Util.debug_matrix \"m_Phi\" m_Phi;\n",
|
||||
" *)\n",
|
||||
" (* Calcul de tous les alpha -> Matrice *)\n",
|
||||
" let f_alpha m_Phi =\n",
|
||||
"\n",
|
||||
@ -1614,49 +1640,70 @@
|
||||
" 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_Phi.{a,i} *. m_Phi.{b,i} *. m_Phi.{e,i} *. m_Phi.{f,j} \n",
|
||||
" integral_general (fun a b e f i j -> m_Phi.{a,i} *. m_Phi.{b,j} *. 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_mo ( fun i j ->\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",
|
||||
"\n",
|
||||
" in\n",
|
||||
" in \n",
|
||||
" (* Calcul de D *)\n",
|
||||
" let s_D = \n",
|
||||
" let v_D = \n",
|
||||
" let m_D = 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,i} \n",
|
||||
" ) i j\n",
|
||||
" )\n",
|
||||
" in Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
|
||||
" in Vec.sum v_D \n",
|
||||
" \n",
|
||||
" in \n",
|
||||
" (* Matrice de rotation 2 par 2 *)\n",
|
||||
" let indice_i, indice_j, m_R = \n",
|
||||
" let indice_i, indice_j, m_R = \n",
|
||||
"\n",
|
||||
" \n",
|
||||
" let m_alpha = f_alpha m_Phi\n",
|
||||
" in\n",
|
||||
" (*\n",
|
||||
" Util.debug_matrix \"m_alpha\" m_alpha;\n",
|
||||
" *)\n",
|
||||
" (* Valeur du plus grand angle alpha*)\n",
|
||||
" let alpha = \n",
|
||||
"\n",
|
||||
" (* Extraction du plus grand angle alpha*)\n",
|
||||
" (* Extraction du plus grand angle alpha*)\n",
|
||||
" let ij_max_alpha = Mat.as_vec m_alpha\n",
|
||||
" |> iamax\n",
|
||||
" in\n",
|
||||
"\n",
|
||||
" (* Indices du plus grand angle alpha *)\n",
|
||||
" let indice_i = ij_max_alpha mod n_ao\n",
|
||||
" let indice_i = ij_max_alpha mod n_mo\n",
|
||||
" in\n",
|
||||
" let indice_j = ij_max_alpha / n_ao + 1\n",
|
||||
" let indice_j = ij_max_alpha / n_mo + 1\n",
|
||||
"\n",
|
||||
" in\n",
|
||||
" m_alpha.{indice_i,indice_j}\n",
|
||||
" (*\n",
|
||||
" Printf.printf \"%i\\n%!\" indice_i;\n",
|
||||
" Printf.printf \"%i\\n%!\" indice_j;\n",
|
||||
" *)\n",
|
||||
" let alpha = \n",
|
||||
" m_alpha.{indice_i,indice_j}\n",
|
||||
" in\n",
|
||||
" \n",
|
||||
" (*\n",
|
||||
" Printf.printf \"%f\\n%!\" alpha;\n",
|
||||
" \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",
|
||||
" (*\n",
|
||||
" Util.debug_matrix \"m_R\" m_R;\n",
|
||||
" Printf.printf \"%d %d\\n%!\" indice_i indice_j;\n",
|
||||
" *)\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, m_Ksi_tilde = \n",
|
||||
@ -1672,9 +1719,10 @@
|
||||
" m_Ksi, gemm m_Ksi m_R\n",
|
||||
"\n",
|
||||
" in\n",
|
||||
" \n",
|
||||
" (*\n",
|
||||
" Util.debug_matrix \"m_Ksi\" m_Ksi;\n",
|
||||
" Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n",
|
||||
" *)\n",
|
||||
" (* Construction de la nouvelle matrice d'OMs phi~ *)\n",
|
||||
" let m_Phi_tilde m_Ksi m_Ksi_tilde = \n",
|
||||
"\n",
|
||||
@ -1700,20 +1748,31 @@
|
||||
" in \n",
|
||||
" let coef = m_Phi_tilde m_Ksi m_Ksi_tilde\n",
|
||||
" in\n",
|
||||
" \n",
|
||||
" Util.debug_matrix \"coef\" coef;\n",
|
||||
" let loc = alpha \n",
|
||||
" \n",
|
||||
" let loc = s_D \n",
|
||||
" in\n",
|
||||
" let result = \n",
|
||||
" \n",
|
||||
" Printf.printf \"%f\\n%!\" loc;\n",
|
||||
" Printf.printf \"%f\\n%!\" prev_iter.loc;\n",
|
||||
" \n",
|
||||
" if abs_float (loc -. prev_iter.loc) < 1.e-5 \n",
|
||||
" then\n",
|
||||
" { coef ; loc }\n",
|
||||
" else\n",
|
||||
" prev_iter\n",
|
||||
" in\n",
|
||||
" new_Phi result (n-1)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" new_Phi {coef; loc;} (n-1)\n",
|
||||
"\n",
|
||||
" \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"new_Phi {coef = m_C ;loc = 0. } 10;;\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
@ -1723,7 +1782,199 @@
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
" new_Phi prev_iter 3;;"
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"(* Fonction pour itérer , cc critère de convergence*) \n",
|
||||
"(* Localisation de Edminstion *)\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 sum a = \n",
|
||||
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
|
||||
";;\n",
|
||||
"\n",
|
||||
"(*Calcul des intégrales*) \n",
|
||||
"let integral_general g i j =\n",
|
||||
"Array.map (fun a ->\n",
|
||||
" let v = \n",
|
||||
" Array.map (fun b ->\n",
|
||||
" let u = \n",
|
||||
" Array.map (fun e ->\n",
|
||||
" let t = Array.map (fun f ->\n",
|
||||
" (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
|
||||
" ) (Util.array_range 1 n_ao)\n",
|
||||
" in sum t\n",
|
||||
" ) (Util.array_range 1 n_ao)\n",
|
||||
" in sum u\n",
|
||||
" ) (Util.array_range 1 n_ao)\n",
|
||||
" in sum v\n",
|
||||
") (Util.array_range 1 n_ao)\n",
|
||||
"|> sum \n",
|
||||
"\n",
|
||||
"let m_Phi = None;;\n",
|
||||
"type iteration = \n",
|
||||
"{ coef : Mat.t ;\n",
|
||||
" loc : float ;\n",
|
||||
"}\n",
|
||||
"\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",
|
||||
" Printf.printf \"%d\\n%!\" n;\n",
|
||||
" (*\n",
|
||||
" Util.debug_matrix \"prev_iter\" prev_iter.coef;\n",
|
||||
" *)\n",
|
||||
" if n == 0 then\n",
|
||||
" prev_iter\n",
|
||||
" else\n",
|
||||
" let m_Phi = prev_iter.coef\n",
|
||||
" in\n",
|
||||
" (*\n",
|
||||
" Util.debug_matrix \"m_coef\" prev_iter.coef;\n",
|
||||
" Util.debug_matrix \"m_Phi\" m_Phi;\n",
|
||||
" *)\n",
|
||||
" (* Calcul de tous les alpha -> Matrice *)\n",
|
||||
" let f_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_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_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_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",
|
||||
" in \n",
|
||||
" (* Matrice de rotation 2 par 2 *)\n",
|
||||
" let indice_i, indice_j, m_R = \n",
|
||||
" let n = 3\n",
|
||||
" in\n",
|
||||
" let rec new_alpha m_alpha n =\n",
|
||||
" \n",
|
||||
" let m_alpha = f_alpha m_Phi\n",
|
||||
" in\n",
|
||||
" (*\n",
|
||||
" Util.debug_matrix \"m_alpha\" m_alpha;\n",
|
||||
" *)\n",
|
||||
" (* Valeur du plus grand angle alpha*)\n",
|
||||
" (* Extraction du plus grand angle alpha*)\n",
|
||||
" let ij_max_alpha = Mat.as_vec m_alpha\n",
|
||||
" |> iamax\n",
|
||||
" in\n",
|
||||
"\n",
|
||||
" (* Indices du plus grand angle alpha *)\n",
|
||||
" let indice_i = ij_max_alpha mod n_mo\n",
|
||||
" in\n",
|
||||
" let indice_j = ij_max_alpha / n_mo + 1\n",
|
||||
"\n",
|
||||
" \n",
|
||||
" (*\n",
|
||||
" Printf.printf \"%i\\n%!\" indice_i;\n",
|
||||
" Printf.printf \"%i\\n%!\" indice_j;\n",
|
||||
" *)\n",
|
||||
" in\n",
|
||||
" if m_alpha.{indice_i,indice_j} < (3.14 /. 2.)\n",
|
||||
" then m_alpha\n",
|
||||
" else new_alpha (Mat.sub m_alpha (Mat.make n_ao n_ao (3.14 /. 2.)))\n",
|
||||
" \n",
|
||||
" in\n",
|
||||
" let alpha = m_alpha.{indice_i, indice_j}\n",
|
||||
" \n",
|
||||
" \n",
|
||||
" (*\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",
|
||||
" (*\n",
|
||||
" Util.debug_matrix \"m_R\" m_R;\n",
|
||||
" Printf.printf \"%d %d\\n%!\" indice_i indice_j;\n",
|
||||
" *)\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, 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",
|
||||
" let p = indice_i\n",
|
||||
" in\n",
|
||||
" let q = indice_j\n",
|
||||
" 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",
|
||||
" 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",
|
||||
" *)\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_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",
|
||||
" (* Matrice intérmédiaire où les orbitales i et j ont été supprimées et remplacées par des 0, par soustraction de la \n",
|
||||
" matrice Phi par la matrice *)\n",
|
||||
" 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_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",
|
||||
" in\n",
|
||||
" Mat.sub m_Phi m_Psi\n",
|
||||
" in\n",
|
||||
" Mat.add m_Psi_tilde m_interm\n",
|
||||
" in \n",
|
||||
" let coef = m_Phi_tilde m_Ksi m_Ksi_tilde\n",
|
||||
" in\n",
|
||||
" \n",
|
||||
" Util.debug_matrix \"coef\" coef;\n",
|
||||
" \n",
|
||||
" let loc = alpha \n",
|
||||
" in\n",
|
||||
" (*\n",
|
||||
" Printf.printf \"%f\\n%!\" loc;\n",
|
||||
" Printf.printf \"%f\\n%!\" prev_iter.loc;\n",
|
||||
" *)\n",
|
||||
" if abs_float (loc -. prev_iter.loc) < 1.e-10 \n",
|
||||
" then\n",
|
||||
" { coef ; loc }\n",
|
||||
" else\n",
|
||||
" new_Phi {coef; loc;} (n-1)\n",
|
||||
"\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
|
||||
Loading…
Reference in New Issue
Block a user