Finished F12Fast notebook

This commit is contained in:
Anthony Scemama 2020-01-07 18:47:54 +01:00
parent 42f8edc594
commit 85a4425a6e
1 changed files with 463 additions and 203 deletions

View File

@ -1308,7 +1308,6 @@
" | _ -> assert false\n",
" in\n",
" \n",
" (* MOs unoccupied in both I and J *)\n",
" let mos_novirt, mos_novirt' =\n",
" let alfa = \n",
" let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
@ -1410,103 +1409,424 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 129,
"metadata": {},
"outputs": [],
"outputs": [
{
"data": {
"text/plain": [
"val m_2112_1H_2Fa :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2112_1H_2Fb :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2112_2Ha_2Fa :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2112_2Hb_2Fa :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2112_2Ha_2Fb :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2121_2Ha_2Fa :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2121_2Hb_2Fa :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2121_2Ha_2Fb :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2122_2Ha_2Fa_ij :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2122_2Hb_2Fb_ij :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2122_2Hb_2Fb_ij2 :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2122_2Ha_2Fa_ij2 :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2122_2Ha_2Fa_nv :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2122_2Hb_2Fb_nv :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val m_2122_2Hb_2Fb_nv2 :\n",
" (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t =\n",
" <abstr>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val integral_value : Determinant.t -> Determinant.t -> float = <fun>\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Checking ... \n",
" - 10 %\n",
" - 20 %\n",
" - 30 %\n",
" - 40 %\n",
" - 50 %\n",
" - 60 %\n",
" - 70 %\n",
" - 80 %\n",
" - 90 %\n",
" - 100 %\n",
"OK: (2, 45) | 2.332967e-02 2.332967e-02 | 2.428613e-17 | cpu : 109.679883 0.039116\n"
]
},
{
"data": {
"text/plain": [
"- : unit = ()\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val ki : Determinant.t =\n",
" phase:+1\n",
" a:+1 +--++-----------------------------------------------------------\n",
" b:+1 +++-------------------------------------------------------------\n",
"\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"val kj : Determinant.t =\n",
" phase:+1\n",
" a:+1 +-+-+-----------------------------------------------------------\n",
" b:+1 ++-+------------------------------------------------------------\n",
"\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"- : float = 0.000896378780435275082\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"- : float = 0.000896378780435274107\n"
]
},
"execution_count": 129,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"let m_2112_1H_2Fa =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_one i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +.\n",
" h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Alfa )\n",
" h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Alfa +. \n",
" h_two i j a l Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.\n",
" h_two i j k a Spin.Alfa Spin.Alfa *. f_one a l Spin.Alfa +.\n",
" sum mos_in (fun m -> -. h_two i j a m Spin.Alfa Spin.Alfa *.\n",
" f_two m a k l Spin.Alfa Spin.Alfa) +.\n",
" sum mos_cabs (fun b -> if b >= a then 0. else\n",
" h_two i j a b Spin.Alfa Spin.Alfa *. f_two a b k l Spin.Alfa Spin.Alfa \n",
" )\n",
" )\n",
" )\n",
"\n",
"let m_2112_1H_2Fb =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_one i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.\n",
" h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Beta)\n",
" h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Beta +.\n",
" h_two i j a l Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.\n",
" h_two i j k a Spin.Alfa Spin.Beta *. f_one a l Spin.Alfa +.\n",
" sum mos_in (fun m ->\n",
" h_two i j a m Spin.Alfa Spin.Beta *. f_two a m k l Spin.Alfa Spin.Beta +.\n",
" h_two i j m a Spin.Alfa Spin.Beta *. f_two m a k l Spin.Alfa Spin.Beta ) +.\n",
" sum mos_cabs (fun b -> \n",
" h_two i j a b Spin.Alfa Spin.Beta *. f_two a b k l Spin.Alfa Spin.Beta \n",
" )\n",
" )\n",
" )\n",
"\n",
"let m_2112_2Ha_2Fa =\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i n a n Spin.Alfa Spin.Alfa *.\n",
" f_two a j k l Spin.Alfa Spin.Alfa )\n",
" h_two i n a n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +. \n",
" h_two j n a n Spin.Alfa Spin.Alfa *. f_two a i l k Spin.Alfa Spin.Alfa \n",
" )\n",
" )\n",
" \n",
"let m_2112_2Hb_2Fa =\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i n a n Spin.Alfa Spin.Beta *.\n",
" f_two a j k l Spin.Alfa Spin.Alfa )\n",
" h_two i n a n Spin.Alfa Spin.Beta *. f_two a j k l Spin.Alfa Spin.Alfa +.\n",
" h_two j n a n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Alfa \n",
" )\n",
" )\n",
" \n",
"let m_2112_2Ha_2Fb =\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i n a n Spin.Alfa Spin.Alfa *.\n",
" f_two a j k l Spin.Alfa Spin.Beta )\n",
" h_two i n a n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.\n",
" h_two j n a n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Beta )\n",
" )\n",
" \n",
"let m_2112_2Hb_2Fb =\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i n a n Spin.Alfa Spin.Beta *.\n",
" f_two a j k l Spin.Alfa Spin.Beta )\n",
" )\n",
" \n",
"\n",
"let m_2121_2Ha_1F =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j a l Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.\n",
" h_two i j k a Spin.Alfa Spin.Alfa *. f_one a l Spin.Alfa)\n",
" )\n",
"\n",
"let m_2121_2Hb_1F =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j a l Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.\n",
" h_two i j k a Spin.Alfa Spin.Beta *. f_one a l Spin.Alfa)\n",
" )\n",
"\n",
"let m_2121_2Ha_2Fa =\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j a l Spin.Alfa Spin.Alfa *.\n",
" f_two a n k n Spin.Alfa Spin.Alfa )\n",
" h_two i j a l Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Alfa +.\n",
" h_two j i a k Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Alfa \n",
" )\n",
" )\n",
" \n",
"\n",
"let m_2121_2Hb_2Fa =\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j a l Spin.Alfa Spin.Beta *.\n",
" f_two a n k n Spin.Alfa Spin.Alfa )\n",
" h_two i j a l Spin.Alfa Spin.Beta *. f_two a n k n Spin.Alfa Spin.Alfa +.\n",
" h_two j i a k Spin.Alfa Spin.Beta *. f_two a n l n Spin.Alfa Spin.Beta \n",
" )\n",
" )\n",
" \n",
"let m_2121_2Ha_2Fb =\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j a l Spin.Alfa Spin.Alfa *.\n",
" f_two a n k n Spin.Alfa Spin.Beta )\n",
" h_two i j a l Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Beta +.\n",
" h_two j i a k Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Beta \n",
" )\n",
" )\n",
" \n",
"let m_2121_2Hb_2Fb =\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
"\n",
"let m_2122_2Ha_2Fa_ij =\n",
" let s = Spin.Alfa in\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j a l Spin.Alfa Spin.Beta *.\n",
" f_two a n k n Spin.Alfa Spin.Beta )\n",
" h_two i n a k s s *. f_two j a n l s s \n",
" +. h_two i n a l s s *. f_two j a k n s s \n",
" -. h_two j n a k s s *. f_two i a n l s s \n",
" -. h_two j n a l s s *. f_two i a k n s s\n",
" )\n",
" )\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
" \n",
"let m_2122_2Hb_2Fb_ij =\n",
" let s, s' = Spin.(Alfa, Beta) in\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two n i a k s s *. f_two a j n l s s'\n",
" +. h_two n j a l s s' *. f_two i a k n s s \n",
" -. h_two n j k a s s' *. f_two i a n l s s'\n",
" )\n",
" )\n",
" \n",
"let m_2122_2Hb_2Fb_ij2 =\n",
" let s, s' = Spin.(Alfa, Beta) in\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" -. h_two i n a l s s' *. f_two a j k n s s' +.\n",
" (if n < j then\n",
" h_two i n k a s s' *. f_two j a l n s' s'\n",
" +. h_two n j a l s' s' *. f_two i a k n s s'\n",
" else \n",
" -. h_two i n k a s s' *. f_two j a n l s' s' \n",
" -. h_two j n a l s' s' *. f_two i a k n s s'\n",
" )\n",
" )\n",
" )\n",
" \n",
"let m_2122_2Ha_2Fa_ij2 =\n",
" let s, s' = Spin.(Alfa, Beta) in\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i n k a s s' *. f_two j a l n s s'\n",
" +. h_two j n l a s s' *. f_two i a k n s s'\n",
" -. h_two i n l a s s' *. f_two j a k n s s'\n",
" -. h_two j n k a s s' *. f_two i a l n s s'\n",
" )\n",
" ) \n",
" \n",
"let m_2122_2Ha_2Fa_nv =\n",
" let s = Spin.Alfa in\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a -> h_two i j a n s s *. f_two n a k l s s ) ) \n",
" \n",
" \n",
" \n",
"let m_2122_2Hb_2Fb_nv =\n",
" let s, s' = Spin.(Alfa, Beta) in\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a -> -. h_two i j a n s s' *. f_two a n k l s s' ) ) \n",
" \n",
"let m_2122_2Hb_2Fb_nv2 =\n",
" let s, s' = Spin.(Alfa, Beta) in\n",
" array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->\n",
" sum mos_cabs (fun a -> -. h_two i j n a s s' *. f_two n a k l s s' ) ) \n",
" \n",
" \n",
" \n",
"let integral_value ki kj = \n",
" let i, j, k, l, s, s', phase =\n",
" match Excitation.of_det ki kj with\n",
@ -1541,158 +1861,59 @@
" | Spin.Alfa -> alfa, beta\n",
" | Spin.Beta -> beta, alfa\n",
" in\n",
"\n",
" \n",
" let mos_virt_a, mos_virt_b = \n",
" Array.init mo_num (fun i -> Some (i+1)) , \n",
" Array.init mo_num (fun i -> Some (i+1)) \n",
" in\n",
" List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);\n",
" List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);\n",
" List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);\n",
" List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);\n",
"\n",
" let mos_virt_a, mos_virt_b = \n",
" Array.to_list mos_virt_a |> Util.list_some,\n",
" Array.to_list mos_virt_b |> Util.list_some\n",
" \n",
" let mos_novirt, mos_novirt' =\n",
" let alfa = \n",
" let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
" let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
" Bitstring.to_list (Bitstring.logor i j)\n",
" in\n",
" let beta = \n",
" let i = Spindeterminant.bitstring @@ Determinant.beta ki in\n",
" let j = Spindeterminant.bitstring @@ Determinant.beta kj in\n",
" Bitstring.to_list (Bitstring.logor i j)\n",
" in\n",
" match s with\n",
" | Spin.Alfa -> alfa, beta\n",
" | Spin.Beta -> beta, alfa\n",
" in\n",
" \n",
" let mos_virt, mos_virt' =\n",
" match s with\n",
" | Spin.Alfa -> mos_virt_a, mos_virt_b\n",
" | Spin.Beta -> mos_virt_b, mos_virt_a\n",
" in\n",
" \n",
" \n",
" \n",
" let result = \n",
" let result = \n",
" if s = s' then\n",
" let s'' = Spin.other s in\n",
" m_2112_1H_2Fa.{i,j,k,l} +. \n",
" sum mos_i (fun n ->\n",
" m_2112_2Ha_2Fa.{i,j,k,l,n} +. m_2112_2Ha_2Fa.{j,i,l,k,n}\n",
" ) +.\n",
" sum mos_i' (fun n ->\n",
" m_2112_2Hb_2Fa.{i,j,k,l,n} +. m_2112_2Hb_2Fa.{j,i,l,k,n}\n",
" ) +.\n",
" m_2121_2Ha_1F.{i,j,k,l} +. \n",
" sum mos_j (fun n ->\n",
" m_2121_2Ha_2Fa.{i,j,k,l,n} +. m_2121_2Ha_2Fa.{j,i,l,k,n}\n",
" ) +.\n",
" sum mos_j' (fun n ->\n",
" m_2121_2Ha_2Fb.{i,j,k,l,n} +. m_2121_2Ha_2Fb.{j,i,l,k,n}\n",
" ) +.\n",
" sum mos_cabs (fun a -> \n",
" sum mos_ij (fun n ->\n",
" h_two i n a k s s *. f_two j a n l s s \n",
" +. h_two i n a l s s *. f_two j a k n s s \n",
" -. h_two j n a k s s *. f_two i a n l s s \n",
" -. h_two j n a l s s *. f_two i a k n s s\n",
" )\n",
" +. sum mos_virt (fun m ->\n",
" -. h_two i j a m s s *. f_two m a k l s s )\n",
" +. sum mos_ij' (fun n ->\n",
" h_two i n k a s s'' *. f_two j a l n s s''\n",
" +. h_two j n l a s s'' *. f_two i a k n s s''\n",
" -. h_two i n l a s s'' *. f_two j a k n s s''\n",
" -. h_two j n k a s s'' *. f_two i a l n s s''\n",
" )\n",
" )\n",
" sum mos_i (fun n -> m_2112_2Ha_2Fa.{n,i,j,k,l} ) +.\n",
" sum mos_i' (fun n -> m_2112_2Hb_2Fa.{n,i,j,k,l} ) +.\n",
" sum mos_j (fun n -> m_2121_2Ha_2Fa.{n,i,j,k,l} ) +.\n",
" sum mos_j' (fun n -> m_2121_2Ha_2Fb.{n,i,j,k,l} ) +.\n",
" sum mos_ij (fun n -> m_2122_2Ha_2Fa_ij.{n,i,j,k,l} ) +. \n",
" sum mos_ij' (fun n -> m_2122_2Ha_2Fa_ij2.{n,i,j,k,l} ) +. \n",
" sum mos_novirt (fun n -> m_2122_2Ha_2Fa_nv.{n,i,j,k,l} ) \n",
" else\n",
" m_2112_1H_2Fb.{i,j,k,l} +.\n",
" sum mos_i (fun n ->\n",
" m_2112_2Ha_2Fb.{i,j,k,l,n} +. m_2112_2Hb_2Fb.{j,i,l,k,n}\n",
" ) +.\n",
" sum mos_i' (fun n ->\n",
" m_2112_2Hb_2Fb.{i,j,k,l,n} +. m_2112_2Ha_2Fb.{j,i,l,k,n}\n",
" ) +.\n",
" m_2121_2Hb_1F.{i,j,k,l} +.\n",
" sum mos_j (fun n ->\n",
" m_2121_2Hb_2Fa.{i,j,k,l,n} +. m_2121_2Hb_2Fb.{j,i,l,k,n}\n",
" ) +.\n",
" sum mos_j' (fun n ->\n",
" m_2121_2Hb_2Fb.{i,j,k,l,n} +. m_2121_2Hb_2Fa.{j,i,l,k,n}\n",
" ) +.\n",
" sum mos_cabs (fun a ->\n",
" sum mos_virt' (fun m ->\n",
" h_two i j a m s s' *. f_two a m k l s s' ) +.\n",
" sum mos_virt (fun m ->\n",
" h_two i j m a s s' *. f_two m a k l s s' ) +.\n",
" sum mos_ij (fun n ->\n",
" h_two n i a k s s *. f_two a j n l s s'\n",
" +. h_two n j a l s s' *. f_two i a k n s s \n",
" -. h_two n j k a s s' *. f_two i a n l s s'\n",
" ) +.\n",
" sum mos_ij' (fun n -> if n >= j then 0. else \n",
" h_two i n k a s s' *. f_two j a l n s' s'\n",
" +. h_two n j a l s' s' *. f_two i a k n s s' ) +. \n",
" sum mos_ij' (fun n -> if n <= j then 0. else \n",
" -. h_two i n k a s s' *. f_two j a n l s' s' \n",
" -. h_two j n a l s' s' *. f_two i a k n s s' ) +.\n",
" sum mos_ij' (fun n -> \n",
" -. h_two i n a l s s' *. f_two a j k n s s'\n",
" )\n",
" )\n",
" m_2112_1H_2Fb.{i,j,k,l} +. \n",
" sum mos_i (fun n -> m_2112_2Ha_2Fb.{n,i,j,k,l} ) +.\n",
" sum mos_i' (fun n -> m_2112_2Ha_2Fb.{n,j,i,l,k} ) +.\n",
" sum mos_j (fun n -> m_2121_2Hb_2Fa.{n,i,j,k,l} ) +.\n",
" sum mos_j' (fun n -> m_2121_2Hb_2Fa.{n,j,i,l,k} ) +.\n",
" sum mos_novirt'(fun n -> m_2122_2Hb_2Fb_nv.{n,i,j,k,l} ) +.\n",
" sum mos_novirt (fun n -> m_2122_2Hb_2Fb_nv2.{n,i,j,k,l} )+.\n",
" sum mos_ij (fun n -> m_2122_2Hb_2Fb_ij.{n,i,j,k,l} ) +. \n",
" sum mos_ij' (fun n -> m_2122_2Hb_2Fb_ij2.{n,i,j,k,l} ) \n",
" in\n",
"\n",
" match phase with\n",
" | Phase.Pos -> result\n",
" | Phase.Neg -> -. result\n",
" \n",
"let _ = check 100 integral_value 2 \n",
"\n",
"let _ = check 100 integral_value 2 1 \n",
"(*\n",
"\n",
"let ki = det_I.(2)\n",
"let kj = det_I.(33)\n",
"let ki = det_I.(7)\n",
"let kj = det_I.(89)\n",
"\n",
"let _ = integral_value ki kj\n",
"let _ = \n",
" let alphas = generate_alphas ki kj 2 1\n",
" let alphas = generate_alphas ki kj 2 \n",
" in compute_HaaF ki alphas kj\n",
"*)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 12. 2 2 "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let integral_value ki kj = \n",
" let h, h', p, p', s, s', phase =\n",
" match Excitation.of_det ki kj with\n",
" | Excitation.(Double (phase,\n",
" { hole=h ; particle=p ; spin=s },\n",
" { hole=h'; particle=p'; spin=s'}) ) -> h, h', p, p', s, s', phase\n",
" | _ -> assert false\n",
" in\n",
"\n",
" let result = \n",
" if s <> s' then (* Alpha-Beta *)\n",
" sum mos_cabs (fun b -> \n",
" sum mos_cabs (fun a -> \n",
" h_two h h' a b s s' *. f_two a b p p' s s'\n",
" )) \n",
" else (* Alpha-Alpha / Beta-Beta *)\n",
" sum mos_cabs (fun b -> \n",
" sum mos_cabs (fun a -> if b >= a then 0. else\n",
" h_two h h' a b s s' *. f_two a b p p' s s'\n",
" )) \n",
" in\n",
" match phase with\n",
" | Phase.Pos -> result\n",
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 100 integral_value 2 2 2 2"
"\n"
]
},
{
@ -1704,9 +1925,48 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 130,
"metadata": {},
"outputs": [],
"outputs": [
{
"data": {
"text/plain": [
"val integral_value : Determinant.t -> Determinant.t -> float = <fun>\n"
]
},
"execution_count": 130,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Checking ... \n",
" - 10 %\n",
" - 20 %\n",
" - 30 %\n",
" - 40 %\n",
" - 50 %\n",
" - 60 %\n",
" - 70 %\n",
" - 80 %\n",
" - 90 %\n",
" - 100 %\n",
"OK: (13, 52) | 6.122906e-04 6.122906e-04 | 3.252607e-19 | cpu : 77.238524 0.119682\n"
]
},
{
"data": {
"text/plain": [
"- : unit = ()\n"
]
},
"execution_count": 130,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"let integral_value ki kj = \n",
" let i, j, m, k, l, n, s1, s2, s3, phase =\n",
@ -1765,7 +2025,14 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 200 integral_value 3 "
"let _ = check 100 integral_value 3 "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---------------------------"
]
},
{
@ -1909,13 +2176,6 @@
" let x = (integral_value ki kj) in\n",
" Printf.printf \"%20.8e %20.8e %20.8e\\n%!\" x v (v-. x)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {