diff --git a/Notebooks/F12_matrix.ipynb b/Notebooks/F12_matrix.ipynb index 903c6ea..99a9fd7 100644 --- a/Notebooks/F12_matrix.ipynb +++ b/Notebooks/F12_matrix.ipynb @@ -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" ] }, {