From efe0be4bb2050baf8f524060ff02cb585e853869 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 12 Dec 2019 18:41:14 +0100 Subject: [PATCH] Working on f12 notebook --- Notebooks/F12_matrix.ipynb | 310 +++++++++++++++++-------------------- 1 file changed, 145 insertions(+), 165 deletions(-) diff --git a/Notebooks/F12_matrix.ipynb b/Notebooks/F12_matrix.ipynb index 4e2aeca..c52f8a8 100644 --- a/Notebooks/F12_matrix.ipynb +++ b/Notebooks/F12_matrix.ipynb @@ -694,7 +694,7 @@ " List.map (fun particle ->\n", " List.map (fun hole' -> \n", " List.map (fun particle' ->\n", - " if hole' > hole && particle' < particle then\n", + " if hole' > hole && particle' < particle && hole <> particle && hole' <> particle' then\n", " Some (Determinant.double_excitation Spin.Beta hole particle Spin.Beta hole' particle' ki)\n", " else None\n", " ) mos_abs\n", @@ -706,6 +706,7 @@ " List.map (fun particle ->\n", " List.map (fun hole' -> \n", " List.map (fun particle' ->\n", + " if hole = particle || hole' = particle' then None else\n", " Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Beta hole' particle' ki)\n", " ) mos_abs\n", " ) mos_b \n", @@ -836,7 +837,8 @@ " else\n", " Printf.printf \"Failed: (%d, %d) | %e %e | %e\\n%!\" i j d v error\n", "\n", - " " + " \n", + "let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l" ] }, { @@ -898,27 +900,26 @@ " (* Alpha-Beta *)\n", " let s, s' = Spin.(Alfa, Beta) in \n", " let mos_a, mos_b = mos_a ki, mos_b ki in\n", - " List.fold_left (fun accu a -> accu +.\n", - " List.fold_left (fun accu i -> accu +. \n", + " \n", + " sum mos_cabs (fun a ->\n", + " sum mos_a (fun i -> \n", " (h_one i a s +.\n", - " List.fold_left (fun accu j -> accu +. h_two i j a j s s ) 0. mos_a +. \n", - " List.fold_left (fun accu j -> accu +. h_two i j a j s s') 0. mos_b ) *. \n", + " sum mos_a (fun j -> h_two i j a j s s ) +.\n", + " sum mos_b (fun j -> h_two i j a j s s') ) *. \n", " (f_one a i s +.\n", - " List.fold_left (fun accu j -> accu +. f_two a j i j s s ) 0. mos_a +.\n", - " List.fold_left (fun accu j -> accu +. f_two a j i j s s') 0. mos_b )\n", - " ) 0. mos_a\n", - " ) 0. mos_cabs\n", - " +.\n", - " List.fold_left (fun accu a -> accu +.\n", - " List.fold_left (fun accu i -> accu +. \n", + " sum mos_a (fun j -> f_two a j i j s s ) +.\n", + " sum mos_b (fun j -> f_two a j i j s s') )\n", + " )\n", + " +.\n", + " sum mos_b (fun i -> \n", " (h_one i a s +.\n", - " List.fold_left (fun accu j -> accu +. h_two i j a j s' s') 0. mos_b +. \n", - " List.fold_left (fun accu j -> accu +. h_two i j a j s' s ) 0. mos_a ) *. \n", + " sum mos_b (fun j -> h_two i j a j s' s') +.\n", + " sum mos_a (fun j -> h_two i j a j s' s ) ) *. \n", " (f_one a i s +.\n", - " List.fold_left (fun accu j -> accu +. f_two a j i j s' s ) 0. mos_a +.\n", - " List.fold_left (fun accu j -> accu +. f_two a j i j s' s') 0. mos_b )\n", - " ) 0. mos_b\n", - " ) 0. mos_cabs\n", + " sum mos_a (fun j -> f_two a j i j s' s ) +.\n", + " sum mos_b (fun j -> f_two a j i j s' s') )\n", + " ) \n", + " ) \n", " \n", "let _ =\n", " check 0 integral_value 0 1 1 1\n", @@ -992,40 +993,29 @@ " \n", " (* Alpha-Beta *)\n", " let s, s' = Spin.(Alfa, Beta) in\n", - " List.fold_left (fun accu j -> accu +. \n", - " List.fold_left (fun accu i -> accu +. \n", - " List.fold_left (fun accu k -> accu +. \n", - " List.fold_left (fun accu a ->\n", - " accu +. h_two i j a k s s' *. f_two a k i j s s'\n", - " ) 0. mos_cabs\n", - " ) 0. mos_virt_b\n", - " ) 0. (mos_a ki)\n", - " ) 0. (mos_b ki)\n", + " sum mos_cabs (fun a ->\n", + " sum (mos_b ki) (fun j -> \n", + " sum (mos_a ki) (fun i -> \n", + " sum mos_virt_b (fun k -> \n", + " h_two i j a k s s' *. f_two a k i j s s'\n", + " )))\n", " +.\n", - " List.fold_left (fun accu j -> accu +. \n", - " List.fold_left (fun accu i -> accu +. \n", - " List.fold_left (fun accu k -> accu +. \n", - " List.fold_left (fun accu a ->\n", - " accu +. h_two i j a k s s' *. f_two a k i j s s'\n", - " ) 0. mos_cabs\n", - " ) 0. mos_virt_a\n", - " ) 0. (mos_b ki)\n", - " ) 0. (mos_a ki)\n", + " sum (mos_a ki) (fun j -> \n", + " sum (mos_b ki) (fun i -> \n", + " sum mos_virt_a (fun k -> \n", + " h_two i j a k s s' *. f_two a k i j s s'\n", + " )))) \n", " +.\n", " (* Alpha-Alpha / Beta-Beta *)\n", " List.fold_left (fun accu (mos_virt,mos,s) -> \n", " let s' = s in accu +.\n", - " List.fold_left (fun accu j -> accu +. \n", - " List.fold_left (fun accu i -> accu +. \n", - " if i < j then accu else\n", - " List.fold_left (fun accu k -> accu +. \n", - " List.fold_left (fun accu a -> \n", - " accu +. h_two i j a k s s' *. f_two a k i j s s'\n", - " ) 0. mos_cabs\n", - " ) 0. mos_virt\n", - " ) 0. mos\n", - " ) 0. mos\n", - " ) 0. [ (mos_virt_a,mos_a ki,Spin.Alfa) ; (mos_virt_b,mos_b ki,Spin.Beta)]\n", + " sum mos (fun j -> \n", + " sum mos (fun i -> if i < j then 0. else\n", + " sum mos_virt (fun k -> \n", + " sum mos_cabs (fun a -> \n", + " h_two i j a k s s' *. f_two a k i j s s'\n", + " )))) \n", + " ) 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", @@ -1083,31 +1073,23 @@ "let integral_value ki kj = \n", " (* Alpha-Beta *)\n", " let s, s' = Spin.(Alfa, Beta) in\n", - " List.fold_left (fun accu j -> accu +. \n", - " List.fold_left (fun accu i -> accu +. \n", - " List.fold_left (fun accu b -> accu +. \n", - " List.fold_left (fun accu a ->\n", - " accu +. h_two i j a b s s' *. f_two a b i j s s'\n", - " ) 0. mos_cabs\n", - " ) 0. mos_cabs\n", - " ) 0. (mos_a ki)\n", - " ) 0. (mos_b ki)\n", + " sum mos_cabs (fun a -> \n", + " sum mos_cabs (fun b -> \n", + " sum (mos_a ki) (fun i -> \n", + " sum (mos_b ki) (fun j -> \n", + " h_two i j a b s s' *. f_two a b i j s s'\n", + " ))))\n", " +.\n", " (* Alpha-Alpha / Beta-Beta *)\n", " List.fold_left (fun accu (mos,s) -> \n", " let s' = s in accu +.\n", - " List.fold_left (fun accu j -> accu +. \n", - " List.fold_left (fun accu i -> accu +. \n", - " if i < j then accu else\n", - " List.fold_left (fun accu b -> accu +. \n", - " List.fold_left (fun accu a -> \n", - " if b > a then accu else\n", - " accu +. h_two i j a b s s' *. f_two a b i j s s'\n", - " ) 0. mos_cabs\n", - " ) 0. mos_cabs\n", - " ) 0. mos\n", - " ) 0. mos\n", - " ) 0. [ (mos_a ki,Spin.Alfa) ; (mos_b ki,Spin.Beta)]\n", + " sum mos_cabs (fun b -> \n", + " sum mos_cabs (fun a -> if b > a then 0. else\n", + " sum mos (fun j -> \n", + " sum mos (fun i -> if i < j then 0. else\n", + " h_two i j a b s s' *. f_two a b i j s s'\n", + " ))))\n", + " ) 0. [ (mos_a ki,Spin.Alfa) ; (mos_b ki,Spin.Beta)]\n", " \n", "let _ =\n", " check 0 integral_value 0 2 2 2\n", @@ -1176,14 +1158,15 @@ "\n", " let result = \n", " let s' = Spin.other s in\n", - " List.fold_left (fun accu a -> accu +.\n", + " \n", + " sum mos_cabs (fun a -> \n", " (h_one i a s +. \n", - " (List.fold_left (fun accu j -> accu +. h_two i j a j s s ) 0. mos_i ) +.\n", - " (List.fold_left (fun accu j -> accu +. h_two i j a j s s') 0. mos_i') ) *.\n", + " sum mos_i (fun j -> h_two i j a j s s ) +.\n", + " sum mos_i' (fun j -> h_two i j a j s s') ) *.\n", " (f_one a k s +. \n", - " (List.fold_left (fun accu j -> accu +. f_two a j k j s s ) 0. mos_j ) +.\n", - " (List.fold_left (fun accu j -> accu +. f_two a j k j s s') 0. mos_j') ) \n", - " ) 0. mos_cabs\n", + " sum mos_j (fun j -> f_two a j k j s s ) +.\n", + " sum mos_j' (fun j -> f_two a j k j s s') ) \n", + " )\n", " in\n", " match phase with\n", " | Phase.Pos -> result\n", @@ -1269,21 +1252,21 @@ "\n", " let result = \n", " let s' = Spin.other s in\n", - " List.fold_left (fun accu a -> accu +.\n", - " List.fold_left (fun accu j -> accu +.\n", + " sum mos_cabs (fun a -> \n", + " sum mos_j (fun j -> \n", " (h_one j a s +. \n", - " (List.fold_left (fun accu m -> accu +. h_two j m a m s s ) 0. mos_i ) +.\n", - " (List.fold_left (fun accu m -> accu +. h_two j m a m s s') 0. mos_i') ) *.\n", + " sum mos_i (fun m -> h_two j m a m s s ) +.\n", + " sum mos_i' (fun m -> h_two j m a m s s') ) *.\n", " (f_two a i j k s s )\n", - " ) 0. mos_j\n", + " ) \n", " +.\n", - " List.fold_left (fun accu j -> accu +.\n", + " sum mos_j' (fun j -> \n", " (h_one j a s +. \n", - " (List.fold_left (fun accu m -> accu +. h_two j m a m s' s ) 0. mos_i ) +.\n", - " (List.fold_left (fun accu m -> accu +. h_two j m a m s' s') 0. mos_i') ) *.\n", + " sum mos_i (fun m -> h_two j m a m s' s ) +.\n", + " sum mos_i' (fun m -> h_two j m a m s' s') ) *.\n", " (f_two a i j k s' s )\n", - " ) 0. mos_j'\n", - " ) 0. mos_cabs\n", + " ) \n", + " ) \n", " in\n", " match phase with\n", " | Phase.Pos -> result\n", @@ -1371,21 +1354,20 @@ "\n", " let result = \n", " let s' = Spin.other s in\n", - " List.fold_left (fun accu a -> accu +.\n", - " List.fold_left (fun accu j -> accu +.\n", + " sum mos_cabs (fun a -> \n", + " sum mos_j (fun j -> \n", " (h_two i j k a s s ) *.\n", " (f_one a j s +. \n", - " (List.fold_left (fun accu m -> accu +. f_two m a m j s s) 0. mos_i ) +.\n", - " (List.fold_left (fun accu m -> accu +. f_two m a m j s' s) 0. mos_i') ) \n", - " ) 0. mos_j\n", - " +.\n", - " List.fold_left (fun accu j -> accu +.\n", + " sum mos_i (fun m -> f_two m a m j s s) +.\n", + " sum mos_i' (fun m -> f_two m a m j s' s) ) \n", + " ) \n", + " +.\n", + " sum mos_j' (fun j -> \n", " (h_two i j k a s s') *.\n", " (f_one a j s' +. \n", - " (List.fold_left (fun accu m -> accu +. f_two m a m j s s') 0. mos_i ) +.\n", - " (List.fold_left (fun accu m -> accu +. f_two m a m j s' s') 0. mos_i') ) \n", - " ) 0. mos_j'\n", - " ) 0. mos_cabs\n", + " sum mos_i (fun m -> f_two m a m j s s') +.\n", + " sum mos_i' (fun m -> f_two m a m j s' s') ) \n", + " )) \n", " in\n", " match phase with\n", " | Phase.Pos -> result\n", @@ -1417,26 +1399,32 @@ "source": [ "$|J\\rangle = \\hat{T}_{i}^{k}|I\\rangle$\n", "\n", - "$\\sum_{j} \\sum_{a} \\sum_{m} \\hat{T}_{ij}^{ma} \\hat{T}_{ma}^{kj} + \n", - " \\hat{T}_{lj}^{ka} \\hat{T}_{ia}^{lj} +\n", - " \\hat{T}_{ij}^{am} \\hat{T}_{am}^{kj} \n", - " $\n", + "$\\sum_{j} \\sum_{a} \\sum_{m}\n", + " \\hat{T}_{ij}^{ma} \\hat{T}_{ma}^{kj} + \n", + " \\hat{T}_{lj}^{ka} \\hat{T}_{ia}^{lj} + \\\\\n", + "\\sum_{\\bar{j}} \\sum_{\\bar{a}} \\sum_{m}\n", + " \\hat{T}_{i\\bar{j}}^{m\\bar{a}} \\hat{T}_{m\\bar{a}}^{k\\bar{j}} + \n", + " \\hat{T}_{i\\bar{j}}^{a\\bar{m}} \\hat{T}_{a\\bar{m}}^{k\\bar{j}} +\n", + " \\hat{T}_{l\\bar{j}}^{k\\bar{a}} \\hat{T}_{i\\bar{a}}^{l\\bar{j}}\n", + "$\n", "\n", "$j,l$ : orbital indices of MOs occupied both in $|I\\rangle$ and $|J\\rangle$.\n", "\n", "$m$ : orbital indices of virtual MOs, unoccupied both in $|I\\rangle$ and $|J\\rangle$.\n", "\n", + "$a$ : orbital indices of CABS MOs\n", + "\n", "\\begin{align}\n", - " \\sum_{a} \\sum_{j} \\sum_{m}\n", - " \\langle i j || m a \\rangle \\left[ m a || k j \\right] + \n", - " \\sum_{\\bar{a}} \\sum_{\\bar{j}} \\sum_{m} \n", - " \\langle i \\bar{j} | m \\bar{a} \\rangle \\left[ m \\bar{a} | k \\bar{j} \\right] + \\\\\n", - " \\sum_{a} \\sum_{j} \\sum_{l}\n", - " \\langle l j || k a \\rangle \\left[ i a || l j \\right] + \n", - " \\sum_{\\bar{a}} \\sum_{\\bar{j}} \\sum_{l} \n", - " \\langle l \\bar{j} | k \\bar{a} \\rangle \\left[ i \\bar{a} | l \\bar{j} \\right] + \\\\\n", - " \\sum_{a} \\sum_{\\bar{j}} \\sum_{\\bar{m}} \n", - " \\langle i \\bar{j} | a \\bar{m} \\rangle \\left[ a \\bar{m} | k \\bar{j} \\right] \\\\\n", + "% \\sum_{a} \\sum_{j} \\sum_{m}\n", + "% \\langle i j || m a \\rangle \\left[ m a || k j \\right] + \n", + "% \\sum_{\\bar{a}} \\sum_{\\bar{j}} \\sum_{m} \n", + "% \\langle i \\bar{j} | m \\bar{a} \\rangle \\left[ m \\bar{a} | k \\bar{j} \\right] + \\\\\n", + "% \\sum_{a} \\sum_{j} \\sum_{l}\n", + "% \\langle l j || k a \\rangle \\left[ i a || l j \\right] + \n", + "% \\sum_{\\bar{a}} \\sum_{\\bar{j}} \\sum_{l} \n", + "% \\langle l \\bar{j} | k \\bar{a} \\rangle \\left[ i \\bar{a} | l \\bar{j} \\right] + \\\\\n", + "% \\sum_{a} \\sum_{\\bar{j}} \\sum_{\\bar{m}} \n", + "% \\langle i \\bar{j} | a \\bar{m} \\rangle \\left[ a \\bar{m} | k \\bar{j} \\right] \\\\\n", "\\end{align}" ] }, @@ -1456,7 +1444,7 @@ " | _ -> assert false\n", " in\n", " \n", - " (* mos unoccupied in both I and J *)\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", @@ -1471,10 +1459,10 @@ " Array.to_list mos_virt_b |> Util.list_some\n", " in\n", "\n", - " let mos_virt =\n", + " let mos_virt, mos_virt' =\n", " match s with\n", - " | Spin.Alfa -> mos_virt_a\n", - " | Spin.Beta -> 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_j, mos_j' =\n", @@ -1495,35 +1483,38 @@ "\n", " let result = \n", " let s' = Spin.other s in\n", - " List.fold_left (fun accu a -> accu +.\n", - " List.fold_left (fun accu m -> accu +.\n", - " List.fold_left (fun accu j -> accu +.\n", - " h_two i j m a s s *. f_two m a k j s s ) 0. mos_j\n", - " +.\n", - " List.fold_left (fun accu j -> accu +.\n", - " h_two i j m a s s' *. f_two m a k j s s' ) 0. mos_j'\n", - " ) 0. mos_virt\n", - " -.\n", - " List.fold_left (fun accu l -> accu +.\n", - " List.fold_left (fun accu j -> accu +.\n", - " h_two l j k a s s *. f_two i a l j s s ) 0. mos_j\n", - " +.\n", - " List.fold_left (fun accu j -> accu +.\n", - " h_two l j k a s s' *. f_two i a l j s s' ) 0. mos_j'\n", - " ) 0. mos_j\n", - " +.\n", - " List.fold_left (fun accu m -> accu +.\n", - " List.fold_left (fun accu j -> accu +.\n", - " h_two i j a m s s' *. f_two a m k j s s' ) 0. mos_j'\n", - " ) 0. mos_virt\n", - " ) 0. mos_cabs\n", + " sum mos_cabs (fun a -> \n", + " sum mos_j (fun j -> \n", + " sum mos_virt (fun m -> \n", + " h_two i j m a s s *. f_two m a k j s s \n", + " ))) +.\n", + " sum mos_cabs (fun a -> \n", + " sum mos_j' (fun j -> \n", + " sum mos_virt (fun m -> \n", + " h_two i j m a s s' *. f_two m a k j s s' \n", + " ))) +.\n", + " sum mos_cabs (fun a -> \n", + " sum mos_j' (fun j -> \n", + " sum mos_virt' (fun m -> \n", + " h_two i j a m s s' *. f_two a m k j s s'\n", + " ))) +.\n", + " sum mos_cabs (fun a -> \n", + " sum mos_j (fun j -> \n", + " sum mos_j (fun l -> if l = i then 0. else\n", + " h_two l j k a s s *. f_two i a l j s s \n", + " ))) +. \n", + " sum mos_cabs (fun a -> \n", + " sum mos_j' (fun j -> \n", + " sum mos_j (fun l -> if l = i then 0. else\n", + " h_two l j k a s s' *. f_two i a l j s s'\n", + " ))) \n", " in\n", " match phase with\n", " | Phase.Pos -> result\n", " | Phase.Neg -> -. result\n", " \n", "\n", - "let _ = check 20 integral_value 1 1 2 2\n", + "let _ = check 11 integral_value 1 1 2 2\n", "\n", "\n" ] @@ -1699,25 +1690,17 @@ "\n", " let result = \n", " (* Alpha-Beta *)\n", - " List.fold_left (fun accu j ->\n", - " let s' = Spin.other s in\n", - " accu +. List.fold_left (fun accu b ->\n", - " accu +. List.fold_left (fun accu a ->\n", - " accu +. h_two h j a b s s' *. f_two a b p j s s'\n", - " ) 0. mos_cabs\n", - " ) 0. mos_cabs\n", - " ) 0. mos' \n", + " let s' = Spin.other s in\n", + " sum mos_cabs (fun b -> \n", + " sum mos_cabs (fun a -> \n", + " sum mos' (fun j -> h_two h j a b s s' *. f_two a b p j s s')\n", + " ))\n", " +.\n", " (* Alpha-Alpha / Beta-Beta *)\n", - " List.fold_left (fun accu j ->\n", - " let s' = s in \n", - " accu +. List.fold_left (fun accu b ->\n", - " accu +. List.fold_left (fun accu a -> \n", - " if b > a then accu else\n", - " accu +. h_two h j a b s s' *. f_two a b p j s s'\n", - " ) 0. mos_cabs\n", - " ) 0. mos_cabs\n", - " ) 0. mos \n", + " sum mos_cabs (fun b -> \n", + " sum mos_cabs (fun a -> if b >= a then 0. else\n", + " sum mos (fun j -> h_two h j a b s s *. f_two a b p j s s)\n", + " )) \n", " in\n", " match phase with\n", " | Phase.Pos -> result\n", @@ -1821,18 +1804,15 @@ "\n", " let result = \n", " if s <> s' then (* Alpha-Beta *)\n", - " List.fold_left (fun accu b ->\n", - " accu +. List.fold_left (fun accu a ->\n", - " accu +. h_two h h' a b s s' *. f_two a b p p' s s'\n", - " ) 0. mos_cabs\n", - " ) 0. mos_cabs\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", - " List.fold_left (fun accu b ->\n", - " accu +. List.fold_left (fun accu a -> \n", - " if b > a then accu else\n", - " accu +. h_two h h' a b s s' *. f_two a b p p' s s'\n", - " ) 0. mos_cabs\n", - " ) 0. mos_cabs\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",