Working on F12 nb

This commit is contained in:
Anthony Scemama 2019-12-20 16:09:02 +01:00
parent e7ff462c9b
commit d02620655c
1 changed files with 502 additions and 190 deletions

View File

@ -98,7 +98,9 @@
"#install_printer MOClass.pp ;;\n",
"let pp_mo ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;\n",
"#install_printer pp_mo;;\n",
"(*\n",
"#install_printer DeterminantSpace.pp;;\n",
"*)\n",
"#install_printer SpindeterminantSpace.pp;;\n",
"#install_printer Bitstring.pp;;\n"
]
@ -123,8 +125,8 @@
"metadata": {},
"outputs": [],
"source": [
"let basis_filename = \"/home/scemama/qp2/data/basis/sto-3g\" \n",
"let aux_basis_filename = \"/home/scemama/qp2/data/basis/6-31g\" \n",
"let basis_filename = \"/home/scemama/qp2/data/basis/6-31g\" \n",
"let aux_basis_filename = \"/home/scemama/qp2/data/basis/cc-pvdz\" \n",
"let nuclei = Nuclei.of_zmt_string \"c\" \n",
"let frozen_core = false\n",
"let multiplicity = 1\n",
@ -513,7 +515,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"let ci = CI.make ~n_states:state in_space\n",
@ -803,6 +807,7 @@
"outputs": [],
"source": [
"let check n integral_value exc cabs lexc rexc =\n",
" let cpudet, cpuint = ref 0., ref 0. in\n",
" let det_list =\n",
" match n with\n",
" | 0 -> det_I\n",
@ -822,10 +827,16 @@
" begin\n",
" let alphas = generate_alphas ki kj exc cabs lexc rexc in\n",
" let det_value =\n",
" compute_HaaF ki alphas kj\n",
" let t0 = Unix.gettimeofday () in\n",
" let result = compute_HaaF ki alphas kj in\n",
" cpudet := !cpudet +. Unix.gettimeofday () -. t0;\n",
" result\n",
" in\n",
" let int_value = \n",
" integral_value ki kj\n",
" let t0 = Unix.gettimeofday () in\n",
" let result = integral_value ki kj in\n",
" cpuint := !cpuint +. Unix.gettimeofday () -. t0;\n",
" result\n",
" in\n",
"(* Printf.printf \"%d %d %e %e\\n%!\" i j det_value int_value; *)\n",
" (i,j,det_value,int_value)\n",
@ -854,9 +865,9 @@
" (*\n",
" Printf.printf \"OK: %e\\n%!\" error\n",
" *)\n",
" Printf.printf \"OK: (%d, %d) | %e %e | %e\\n%!\" i j d v error\n",
" Printf.printf \"OK: (%d, %d) | %e %e | %e | cpu : %f %f\\n%!\" i j d v error !cpudet !cpuint\n",
" else\n",
" Printf.printf \"Failed: (%d, %d) | %e %e | %e\\n%!\" i j d v error\n",
" Printf.printf \"Failed: (%d, %d) | %e %e | %e | cpu : %f %f\\n%!\" i j d v error !cpudet !cpuint\n",
"\n",
" \n",
"let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l"
@ -1006,7 +1017,9 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"let m_0111_1H_1F = \n",
@ -1197,7 +1210,7 @@
" ) 0. [ (mos_virt_a,mos_a ki,Spin.Alfa) ; (mos_virt_b,mos_b ki,Spin.Beta)]\n",
" \n",
"let _ =\n",
" check 0 integral_value 0 1 2 2\n",
" check 100 integral_value 0 1 2 2\n",
" "
]
},
@ -1256,7 +1269,7 @@
" result\n",
" \n",
"let _ =\n",
" check 0 integral_value 0 1 2 2\n",
" check 100 integral_value 0 1 2 2\n",
" "
]
},
@ -1330,7 +1343,7 @@
" ) 0. [ (mos_a ki,Spin.Alfa) ; (mos_b ki,Spin.Beta)]\n",
" \n",
"let _ =\n",
" check 0 integral_value 0 2 2 2\n",
" check 100 integral_value 0 2 2 2\n",
" "
]
},
@ -1373,7 +1386,7 @@
" ) 0. [ (mos_a ki) ; (mos_b ki)]\n",
" \n",
"let _ =\n",
" check 0 integral_value 0 2 2 2\n",
" check 100 integral_value 0 2 2 2\n",
" "
]
},
@ -1454,14 +1467,14 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 0 integral_value 1 1 1 1"
"let _ = check 100 integral_value 1 1 1 1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
"scrolled": false
},
"outputs": [],
"source": [
@ -1562,7 +1575,7 @@
" | Phase.Neg -> -. result\n",
" \n",
"let _ =\n",
" check 0 integral_value 1 1 1 1\n",
" check 100 integral_value 1 1 1 1\n",
" "
]
},
@ -1663,7 +1676,7 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 0 integral_value 1 1 1 2"
"let _ = check 100 integral_value 1 1 1 2"
]
},
{
@ -1766,7 +1779,7 @@
" \n",
" \n",
"\n",
"let _ = check 0 integral_value 1 1 1 2"
"let _ = check 100 integral_value 1 1 1 2"
]
},
{
@ -1867,7 +1880,7 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 0 integral_value 1 1 2 1"
"let _ = check 100 integral_value 1 1 2 1"
]
},
{
@ -1970,7 +1983,7 @@
" \n",
" \n",
"\n",
"let _ = check 0 integral_value 1 1 2 1"
"let _ = check 100 integral_value 1 1 2 1"
]
},
{
@ -2025,30 +2038,292 @@
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"let ki = det_I.(6)\n",
"let kj = det_I.(58)\n",
"\n",
"let _ = \n",
" Format.printf \"|I> -> |J> : %a |\\n@.\" Excitation.pp (Excitation.of_det ki kj) ;\n",
"generate_alphas ki kj 1 1 2 2\n",
"|> Array.of_list\n",
"|> Array.mapi (fun kk alpha -> \n",
"\n",
"(*\n",
"let _ = Determinant.to_lists ki\n",
"let _ = Determinant.to_lists alpha\n",
"let _ = Determinant.to_lists kj\n",
"\n",
" \n",
"let _ = \n",
" Format.printf \"|I> -> |J> : %a\\n@.\" Excitation.pp (Excitation.of_det ki kj);\n",
" Format.printf \"|I> -> |a> : %a\\n@.\" Excitation.pp (Excitation.of_det ki alpha);\n",
" Format.printf \"|a> -> |J> : %a\\n@.\" Excitation.pp (Excitation.of_det alpha kj)\n",
"*)\n",
"\n",
" let integral_value ki kj = \n",
" \n",
" let exc0 = Array.init (aux_num+1) (fun _ -> [|\"-\";\"-\"|]) in\n",
" let exc1 = Array.init (aux_num+1) (fun _ -> [|\"-\";\"-\"|]) in \n",
" let exc2 = Array.init (aux_num+1) (fun _ -> [|\"-\";\"-\"|]) in \n",
" \n",
" let i, k, s, phase =\n",
" match Excitation.of_det ki kj with\n",
" | Excitation.(Single (phase, { hole ; particle ; spin } )) ->\n",
" hole, particle, spin, phase\n",
" | _ -> assert false\n",
" in\n",
" let spin = function\n",
" | Spin.Alfa -> 0\n",
" | _ -> 1\n",
" in\n",
" exc0.(i).(spin s ) <- \"i\" \n",
" ; exc0.(k).(spin s ) <- \"k\" \n",
" ;\n",
" let s0 = s in\n",
" \n",
" let i, j, k, l, s, s', p1 =\n",
" match Excitation.of_det ki alpha with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
" \n",
" if exc0.(i).(spin s ) = \"-\" then exc0.(i).(spin s ) <- \"n\";\n",
" if exc0.(j).(spin s') = \"-\" then exc0.(j).(spin s') <- \"n\";\n",
" if exc0.(k).(spin s ) = \"-\" then exc0.(k).(spin s ) <- if k > mo_num then \"a\" else \"m\";\n",
" if exc0.(l).(spin s') = \"-\" then exc0.(l).(spin s') <- if l > mo_num then \"a\" else \"m\";\n",
" \n",
" let string_h = \n",
" Printf.sprintf \"h_two %s %s %s %s %s %s *. \" \n",
" exc0.(i).(spin s )\n",
" exc0.(j).(spin s')\n",
" exc0.(k).(spin s )\n",
" exc0.(l).(spin s')\n",
" (if s = s0 then \"s \" else \"s'\")\n",
" (if s' = s0 then \"s \" else \"s'\")\n",
" (*\n",
" (if exc0.(i).(spin s ) = \"i\" || exc0.(k).(spin s ) = \"k\" then \"s \" else \n",
" if exc0.(i).(spin s ) = \"j\" || exc0.(k).(spin s ) = \"l\" then \"s' \" else \"s''\")\n",
" (if exc0.(j).(spin s') = \"i\" || exc0.(l).(spin s') = \"k\" then \"s \" else \n",
" if exc0.(j).(spin s') = \"j\" || exc0.(l).(spin s') = \"l\" then \"s' \" else \"s''\")\n",
" *)\n",
" in \n",
" \n",
" let i, j, k, l, s, s', p2 =\n",
" match Excitation.of_det alpha kj with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
" \n",
" let string_f = \n",
" Printf.sprintf \"f_two %s %s %s %s %s %s\" \n",
" exc0.(i).(spin s )\n",
" exc0.(j).(spin s')\n",
" exc0.(k).(spin s )\n",
" exc0.(l).(spin s')\n",
" (*\n",
" (if exc0.(i).(spin s ) = \"i\" || exc0.(k).(spin s ) = \"k\" then \"s \" else \n",
" if exc0.(i).(spin s ) = \"j\" || exc0.(k).(spin s ) = \"l\" then \"s' \" else \n",
" if exc0.(i).(spin s ) = \"n\" || exc0.(k).(spin s ) = \"n\" then \"s''\" else\n",
" if s = s0 then \"s\" else \"s'\")\n",
" (if exc0.(j).(spin s') = \"i\" || exc0.(l).(spin s') = \"k\" then \"s \" else \n",
" if exc0.(j).(spin s') = \"j\" || exc0.(l).(spin s') = \"l\" then \"s' \" else \n",
" if exc0.(j).(spin s') = \"n\" || exc0.(l).(spin s') = \"n\" then \"s''\" else\n",
" if s' = s0 then \"s\" else \"s'\")\n",
" *)\n",
" (if s = s0 then \"s \" else \"s'\")\n",
" (if s' = s0 then \"s \" else \"s'\")\n",
" in\n",
" \n",
" Format.printf \"|I> -> |a> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det ki alpha) string_h ;\n",
" Format.printf \"|a> -> |J> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det alpha kj) string_f ;\n",
" \n",
" let sign = \n",
" if Phase.add p1 p2 = phase then \"+. \" else \"-. \"\n",
" in\n",
"sign ^ string_h ^ string_f\n",
"(*\n",
" \n",
" let mos, mos' =\n",
" match s with\n",
" | Spin.Alfa -> mos_a ki, mos_b ki\n",
" | Spin.Beta -> mos_b ki, mos_a ki\n",
" in\n",
"\n",
" let result = \n",
" let k=4 and l=5 in\n",
" let i=1 and j=2 in\n",
" let n=2 and a=7 in\n",
" let s = Spin.Alfa\n",
" and s' = Spin.Alfa\n",
" and s'' = Spin.Beta\n",
" in\n",
" (* \n",
" h_two j n l a s' s'' *. f_two a i n k s s''\n",
" *)\n",
" h_two i n k a s s'' *. f_two j a l n s s''\n",
" in\n",
"\n",
" match phase with\n",
" | Phase.Pos -> result\n",
" | Phase.Neg -> -. result\n",
" *)\n",
"\n",
"in\n",
"integral_value ki kj\n",
"(*\n",
"let a = (compute_HaaF ki [alpha] kj)\n",
"and b = (integral_value ki kj)\n",
"in\n",
"if kk = 31 then\n",
" Format.printf \"%6d %e %e@.@.\" (kk) a b\n",
" *)\n",
")\n",
"|> Array.to_list\n",
"|> List.sort_uniq compare\n",
"\n",
"let _ = compute_HaaF ki (generate_alphas ki kj 1 1 2 2) kj\n",
"\n",
"let _ = integral_value ki kj\n",
"\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
"scrolled": false
},
"outputs": [],
"source": [
"let m_1122_oa =\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 l j k a Spin.Alfa Spin.Alfa *.\n",
" f_two i a l j Spin.Alfa Spin.Alfa \n",
" )\n",
" )\n",
"let integral_value ki kj = \n",
" let i, k, s, phase =\n",
" match Excitation.of_det ki kj with\n",
" | Excitation.(Single (phase, { hole ; particle ; spin })) ->\n",
" hole, particle, spin, phase\n",
" | _ -> assert false\n",
" in\n",
" \n",
" (* MOs unoccupied in both I and J *)\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",
" 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",
" let mos_j, mos_j' =\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.logand 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.logand i j)\n",
" in\n",
" match s with\n",
" | Spin.Alfa -> alfa, beta\n",
" | Spin.Beta -> beta, alfa\n",
" in\n",
" let result = \n",
" let s' = Spin.other s in\n",
" sum mos_cabs (fun a -> \n",
" sum mos_j (fun n -> \n",
" sum mos_virt (fun m -> \n",
" h_two n i a m s s *. f_two m a k n s s \n",
" ))) +.\n",
" sum mos_cabs (fun a' -> \n",
" sum mos_j' (fun n' -> \n",
" sum mos_virt (fun m -> \n",
" h_two i n' m a' s s' *. f_two m a' k n' s s'\n",
" ))) +.\n",
" sum mos_cabs (fun a -> \n",
" sum mos_j' (fun n' -> \n",
" sum mos_virt' (fun m' -> \n",
" h_two i n' a m' s s' *. f_two a m' k n' s s'\n",
" ))) -. \n",
" sum mos_cabs (fun a -> \n",
" sum mos_j (fun n -> \n",
" sum mos_j (fun l -> if l <=n then 0. else\n",
" h_two n l a k s s *. f_two i a l n s s \n",
" ))) -. \n",
" sum mos_cabs (fun a' -> \n",
" sum mos_j' (fun n' -> \n",
" sum mos_j (fun l -> \n",
" h_two l n' k a' s s' *. f_two i a' l n' s s'\n",
" ))) \n",
" in\n",
" match phase with\n",
" | Phase.Pos -> result\n",
" | Phase.Neg -> -. result\n",
" \n",
"let _ = check 10 integral_value 1 1 2 2\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"let m_1122_va =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j m a Spin.Alfa Spin.Alfa *.\n",
" h_two j i a m Spin.Alfa Spin.Alfa *.\n",
" f_two m a k j Spin.Alfa Spin.Alfa \n",
" )\n",
" )\n",
"\n",
"let m_1122_v2 =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j m a Spin.Alfa Spin.Beta *.\n",
" f_two m a k j Spin.Alfa Spin.Beta \n",
" )\n",
" )\n",
"\n",
"let m_1122_v3 =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j a m Spin.Alfa Spin.Beta *.\n",
" f_two a m k j Spin.Alfa Spin.Beta \n",
" )\n",
" )\n",
"\n",
"let m_1122_oa =\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 j l a k Spin.Alfa Spin.Alfa *.\n",
" f_two i a l j Spin.Alfa Spin.Alfa \n",
" )\n",
" )\n",
"\n",
"let m_1122_o =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
" sum mos_cabs (fun a ->\n",
@ -2057,22 +2332,6 @@
" )\n",
" )\n",
"\n",
"let m_1122_v =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j m a Spin.Alfa Spin.Beta *.\n",
" f_two m a k j Spin.Alfa Spin.Beta \n",
" )\n",
" )\n",
"\n",
"let m_1122_v2 =\n",
" array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
" sum mos_cabs (fun a ->\n",
" h_two i j a m Spin.Alfa Spin.Beta *.\n",
" f_two a m k j Spin.Alfa Spin.Beta \n",
" )\n",
" )\n",
"\n",
"\n",
"let integral_value ki kj = \n",
" let i, k, s, phase =\n",
@ -2097,13 +2356,13 @@
" Array.to_list mos_virt_b |> Util.list_some\n",
" in\n",
"\n",
" let mos_virt, mos_virt' =\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",
" in\n",
" \n",
" let mos_j, mos_j' =\n",
" let mos_j, mos_j' =\n",
" let alfa = \n",
" let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
" let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
@ -2117,8 +2376,8 @@
" match s with\n",
" | Spin.Alfa -> alfa, beta\n",
" | Spin.Beta -> beta, alfa\n",
" in\n",
"\n",
" in\n",
" \n",
" let result = \n",
" sum mos_j (fun j -> \n",
" sum mos_virt (fun m -> \n",
@ -2126,18 +2385,18 @@
" )) +.\n",
" sum mos_j' (fun j -> \n",
" sum mos_virt (fun m -> \n",
" m_1122_v.{i,j,k,m}\n",
" m_1122_v2.{i,j,k,m}\n",
" )) +.\n",
" sum mos_j' (fun j -> \n",
" sum mos_virt' (fun m -> \n",
" m_1122_v2.{i,j,k,m}\n",
" m_1122_v3.{i,j,k,m}\n",
" )) -.\n",
" sum mos_j (fun j -> \n",
" sum mos_j (fun l -> if l = i then 0. else\n",
" sum mos_j (fun l -> if l <= j then 0. else\n",
" m_1122_oa.{i,j,k,l}\n",
" )) -. \n",
" sum mos_j' (fun j -> \n",
" sum mos_j (fun l -> if l = i then 0. else\n",
" sum mos_j (fun l -> \n",
" m_1122_o.{i,j,k,l}\n",
" )) \n",
" in\n",
@ -2146,7 +2405,7 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 0 integral_value 1 1 2 2\n",
"let _ = check 100 integral_value 1 1 2 2\n",
"\n",
"\n"
]
@ -2231,14 +2490,14 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 0 integral_value 1 2 2 2"
"let _ = check 100 integral_value 1 2 2 2"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
"scrolled": false
},
"outputs": [],
"source": [
@ -2297,7 +2556,7 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 0 integral_value 1 2 2 2\n",
"let _ = check 100 integral_value 1 2 2 2\n",
"\n",
"\n"
]
@ -2392,7 +2651,7 @@
" | Phase.Neg -> -. result\n",
"\n",
"\n",
"let _ = check 0 integral_value 2 1 1 2"
"let _ = check 100 integral_value 2 1 1 2"
]
},
{
@ -2423,7 +2682,7 @@
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
"scrolled": false
},
"outputs": [],
"source": [
@ -2506,7 +2765,7 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 0 integral_value 2 1 1 2\n",
"let _ = check 100 integral_value 2 1 1 2\n",
"\n",
"\n"
]
@ -2600,14 +2859,14 @@
" | Phase.Neg -> -. result\n",
"\n",
"\n",
"let _ = check 0 integral_value 2 1 2 1"
"let _ = check 100 integral_value 2 1 2 1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
"scrolled": false
},
"outputs": [],
"source": [
@ -2690,7 +2949,7 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 0 integral_value 2 1 2 1\n",
"let _ = check 100 integral_value 2 1 2 1\n",
"\n",
"\n"
]
@ -2699,7 +2958,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"# *11. 2 1 22*"
"# 11. 2 1 22"
]
},
{
@ -2768,41 +3027,67 @@
" \n",
" let result = \n",
" \n",
" let mos, mos_virt =\n",
" let mos_virt, mos_virt' =\n",
" match s with\n",
" | Spin.Alfa -> mos_a ki, mos_virt_a\n",
" | Spin.Beta -> mos_b ki, mos_virt_b\n",
" | Spin.Alfa -> mos_virt_a, mos_virt_b\n",
" | Spin.Beta -> mos_virt_b, mos_virt_a\n",
" in\n",
" \n",
" let mos', mos_virt' =\n",
" match s' with\n",
" | Spin.Alfa -> mos_a ki, mos_virt_a\n",
" | Spin.Beta -> mos_b ki, mos_virt_b\n",
" let mos, mos' =\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.logand 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.logand i j)\n",
" in\n",
" match s with\n",
" | Spin.Alfa -> alfa, beta\n",
" | Spin.Beta -> beta, alfa\n",
" in\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",
"\n",
" if s = s' then \n",
" let s'' = Spin.other s' in\n",
" sum mos_cabs (fun a -> \n",
" sum mos (fun n ->\n",
" h_two i n k a s s *. f_two a j n l s s' +. \n",
" h_two j n l a s s *. f_two a i n k s s' ) +.\n",
" sum mos' (fun n ->\n",
" h_two i n k a s s' *. f_two a j n l s' s') \n",
" (*\n",
" sum mos_virt (fun m -> \n",
" h_two i j m a s s' *. f_two m a k l s s' )\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' (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",
" else\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 (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' (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' (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' (fun n -> \n",
" h_two i n a l s s' *. f_two a j k n s s' ) -.\n",
" sum mos (fun n -> \n",
" h_two n j k a s s' *. f_two i a n l s s' )\n",
" +.\n",
" sum mos' (fun n -> \n",
" h_two n j a l s s' *. f_two a i k n s s ) +.\n",
" sum mos (fun n -> \n",
" h_two n i a k s s' *. f_two a j l n s' s' )\n",
" *)\n",
" )\n",
" -. h_two i n a l s s' *. f_two a j k n s s'\n",
" )\n",
" )\n",
" in\n",
"\n",
" match phase with\n",
@ -2810,117 +3095,144 @@
" | Phase.Neg -> -. result\n",
"\n",
"\n",
"let _ = check 10 integral_value 2 1 2 2"
"let _ = check 300 integral_value 2 1 2 2"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"let ki = det_I.(0)\n",
"let kj = det_I.(9)\n",
"\n",
"let _ = \n",
"generate_alphas ki kj 2 1 2 2\n",
"|> Array.of_list\n",
"|> Array.mapi (fun kk alpha -> \n",
"\n",
"(*\n",
"let _ = Determinant.to_lists ki\n",
"let _ = Determinant.to_lists alpha\n",
"let _ = Determinant.to_lists kj\n",
"let ki = det_I.(99)\n",
"let kj = det_I.(193)\n",
"\n",
"let alpha_to_string alpha = \n",
" \n",
"let _ = \n",
" Format.printf \"|I> -> |J> : %a\\n@.\" Excitation.pp (Excitation.of_det ki kj);\n",
" Format.printf \"|I> -> |a> : %a\\n@.\" Excitation.pp (Excitation.of_det ki alpha);\n",
" Format.printf \"|a> -> |J> : %a\\n@.\" Excitation.pp (Excitation.of_det alpha kj)\n",
"*)\n",
" let i, j, k, l, s, s', phase =\n",
" let exc0 = Array.init (aux_num+1) (fun _ -> [|\"-\";\"-\"|]) in\n",
" \n",
" let i, j, k, l, s, s', phase =\n",
" match Excitation.of_det ki kj with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
"\n",
"\n",
" let integral_value ki kj = \n",
" \n",
" (* \n",
" let i, j, k, l, s, s', phase =\n",
" match Excitation.of_det ki kj with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
" let i, a, s, _ = \n",
" match Excitation.of_det ki alpha with\n",
" | Excitation.(Single (phase, { hole ; particle ; spin })) ->\n",
" hole, particle, spin, phase\n",
" | _ -> assert false\n",
" in\n",
"\n",
" let i', j', k', l', t, u, _ =\n",
" match Excitation.of_det ki alpha with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
" \n",
" let i'', j'', k'', l'', t', u', _ =\n",
" match Excitation.of_det alpha kj with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
" *)\n",
" (*\n",
" let j, k, l, s' = \n",
" if i' = a then\n",
" (assert (k' = k) ; j', k', l', t')\n",
" else\n",
" (assert (j'=a) ; assert (l'=k) ; i', l', k', t)\n",
" in\n",
" *)\n",
" \n",
" Format.printf \"|I> -> |J> : %a |\\n@.\" Excitation.pp (Excitation.of_det ki kj) ;\n",
" Format.printf \"|I> -> |a> : %a |\\n@.\" Excitation.pp (Excitation.of_det ki alpha) ;\n",
" Format.printf \"|a> -> |J> : %a |\\n@.\" Excitation.pp (Excitation.of_det alpha kj);\n",
" in\n",
" let spin = function\n",
" | Spin.Alfa -> 0\n",
" | _ -> 1\n",
" in\n",
" exc0.(i).(spin s ) <- \"i\" \n",
" ; exc0.(j).(spin s') <- \"j\" \n",
" ; exc0.(k).(spin s ) <- \"k\" \n",
" ; exc0.(l).(spin s') <- \"l\" \n",
" ;\n",
" let s0, s0' = s, s' in\n",
" \n",
" \n",
" let mos, mos' =\n",
" match s with\n",
" | Spin.Alfa -> mos_a ki, mos_b ki\n",
" | Spin.Beta -> mos_b ki, mos_a ki\n",
" in\n",
"\n",
" let result = \n",
" let k=5 and l=4 in\n",
" let i=1 and j=2 in\n",
" let n=2 and a=6 in\n",
" let s = Spin.Alfa\n",
" and s' = Spin.Alfa\n",
" and s'' = Spin.Beta\n",
" in\n",
" h_two j n l a s s *. f_two a i n k s s'\n",
" in\n",
"\n",
" match phase with\n",
" | Phase.Pos -> result\n",
" | Phase.Neg -> -. result\n",
"\n",
"in\n",
"let a = (compute_HaaF ki [alpha] kj)\n",
"and b = (integral_value ki kj)\n",
"in\n",
" let i, j, k, l, s, s', p1 =\n",
" match Excitation.of_det ki alpha with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
" \n",
" if exc0.(i).(spin s ) = \"-\" then exc0.(i).(spin s ) <- \"n\";\n",
" if exc0.(j).(spin s') = \"-\" then exc0.(j).(spin s') <- \"n\";\n",
" if exc0.(k).(spin s ) = \"-\" then exc0.(k).(spin s ) <- if k > mo_num then \"a\" else \"m\";\n",
" if exc0.(l).(spin s') = \"-\" then exc0.(l).(spin s') <- if l > mo_num then \"a\" else \"m\";\n",
" \n",
" let string_h = \n",
" Printf.sprintf \"h_two %s %s %s %s %s %s *. \" \n",
" exc0.(i).(spin s )\n",
" exc0.(j).(spin s')\n",
" exc0.(k).(spin s )\n",
" exc0.(l).(spin s')\n",
" (if s = s0 then \"s \" else \"s'\")\n",
" (if s' = s0 then \"s \" else \"s'\")\n",
" in \n",
" \n",
" let i, j, k, l, s, s', p2 =\n",
" match Excitation.of_det alpha kj with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
" \n",
" let string_f = \n",
" Printf.sprintf \"f_two %s %s %s %s %s %s\" \n",
" exc0.(i).(spin s )\n",
" exc0.(j).(spin s')\n",
" exc0.(k).(spin s )\n",
" exc0.(l).(spin s')\n",
" (if s = s0 then \"s \" else \"s'\")\n",
" (if s' = s0 then \"s \" else \"s'\")\n",
" in\n",
"(*\n",
"if a <> b then\n",
"*)\n",
"Format.printf \"%6d %e %e@.\" (kk) a b\n",
")\n",
"\n"
" Format.printf \"|I> -> |a> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det ki alpha) string_h ;\n",
" Format.printf \"|a> -> |J> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det alpha kj) string_f ;\n",
" *)\n",
"\n",
" \n",
" let sign = \n",
" if Phase.add p1 p2 = phase then \"+.\" else \"-.\"\n",
" in\n",
" sign ^ string_h ^ string_f\n",
" \n",
"\n",
"let alpha_debug alpha = \n",
" \n",
" let i, j, k, l, s, s', p1 =\n",
" match Excitation.of_det ki alpha with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
" Printf.printf \"%d %d %d %d \" i j k l;\n",
" \n",
" let i, j, k, l, s, s', p2 =\n",
" match Excitation.of_det alpha kj with\n",
" | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
" hole, hole', particle, particle', spin, spin', phase\n",
" | _ -> assert false\n",
" in\n",
" (*\n",
" Format.printf \"|I> -> |a> : %a | \\n@.\" Excitation.pp (Excitation.of_det ki alpha) ;\n",
" Format.printf \"|a> -> |J> : %a | \\n@.\" Excitation.pp (Excitation.of_det alpha kj) ;\n",
" *)\n",
" Printf.printf \"%d %d %d %d \\n%!\" i j k l\n",
"\n",
" \n",
"let strings = \n",
" Format.printf \"|I> -> |J> : %a |\\n@.\" Excitation.pp (Excitation.of_det ki kj) ;\n",
" generate_alphas ki kj 2 1 2 2\n",
" |> Array.of_list\n",
" |> Array.mapi (fun kk alpha -> alpha_to_string alpha)\n",
" |> Array.to_list\n",
" |> List.sort_uniq compare\n",
" |> Array.of_list\n",
"\n",
"let _ = Array.iteri (fun i x -> Printf.printf \"%d %s \\n%!\" i x) strings\n",
"\n",
"let _ =\n",
" let v =\n",
" let alphas =\n",
" generate_alphas ki kj 2 1 2 2\n",
" (*\n",
" |> List.filter (fun alpha ->\n",
" let x = alpha_to_string alpha in\n",
" x = strings.(6)\n",
" )\n",
" *)\n",
" in\n",
" (*\n",
" List.iter alpha_debug alphas ;\n",
" Printf.printf \"\\n%!\";\n",
" *)\n",
" compute_HaaF ki alphas kj \n",
" in\n",
" let x = (integral_value ki kj) in\n",
" Printf.printf \"%20.8e %20.8e %20.8e\\n%!\" x v (v-. x)\n"
]
},
{