From e66d790cb7f7dcb6cb74f2bd2da57056f1d5ec42 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 12 Dec 2019 16:28:30 +0100 Subject: [PATCH] Working on notebook --- CI/Determinant.ml | 3 + CI/Excitation.ml | 2 +- Notebooks/F12_matrix.ipynb | 279 +++++++++++++++++++++++++++++++++---- 3 files changed, 253 insertions(+), 31 deletions(-) diff --git a/CI/Determinant.ml b/CI/Determinant.ml index a75e01d..a35aba5 100644 --- a/CI/Determinant.ml +++ b/CI/Determinant.ml @@ -84,12 +84,15 @@ let annihilation spin h t = let single_excitation spin h p t = + assert (h <> p); match spin with | Spin.Alfa -> { t with alfa = Spindeterminant.single_excitation h p t.alfa } | Spin.Beta -> { t with beta = Spindeterminant.single_excitation h p t.beta } let double_excitation spin h p spin' h' p' t = + assert (h <> p); + assert (h' <> p'); match spin, spin' with | Spin.(Alfa, Beta) -> { alfa = Spindeterminant.single_excitation h p t.alfa ; beta = Spindeterminant.single_excitation h' p' t.beta } diff --git a/CI/Excitation.ml b/CI/Excitation.ml index 09b40da..013e259 100644 --- a/CI/Excitation.ml +++ b/CI/Excitation.ml @@ -133,7 +133,7 @@ let pp_s_exc ppf t = | Spin.Beta -> "\\beta " ) t.hole t.particle -let pp_exc ppf t = +let pp ppf t = let phase, l = match t with | Identity p -> p, [] diff --git a/Notebooks/F12_matrix.ipynb b/Notebooks/F12_matrix.ipynb index 7ba560b..4e2aeca 100644 --- a/Notebooks/F12_matrix.ipynb +++ b/Notebooks/F12_matrix.ipynb @@ -655,18 +655,21 @@ " let mos_a, mos_b = Determinant.to_lists ki in\n", " [ List.map (fun hole -> \n", " List.map (fun particle ->\n", - " Determinant.single_excitation Spin.Alfa hole particle ki\n", + " if hole = particle then None else\n", + " Some (Determinant.single_excitation Spin.Alfa hole particle ki)\n", " ) mos_abs\n", " ) mos_a\n", " ;\n", " List.map (fun hole -> \n", " List.map (fun particle ->\n", - " Determinant.single_excitation Spin.Beta hole particle ki\n", + " if hole = particle then None else\n", + " Some (Determinant.single_excitation Spin.Beta hole particle ki)\n", " ) mos_abs\n", " ) mos_b \n", " ]\n", " |> List.concat\n", " |> List.concat\n", + " |> Util.list_some\n", " |> List.filter (fun x -> Determinant.is_none x = false)\n", " |> List.map (fun x -> Determinant.set_phase Phase.Pos x)\n", " in\n", @@ -678,6 +681,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", " if hole' > hole && particle' < particle then\n", " Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Alfa hole' particle' ki)\n", " else None\n", @@ -1293,7 +1297,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# 6. *1 1 21*" + "# 6. 1 1 21" ] }, { @@ -1371,16 +1375,16 @@ " List.fold_left (fun accu j -> accu +.\n", " (h_two i j k a s s ) *.\n", " (f_one a j s +. \n", - " (List.fold_left (fun accu m -> accu +. h_two m a m j s s) 0. mos_i ) +.\n", - " (List.fold_left (fun accu m -> accu +. h_two m a m j s' s) 0. mos_i') ) \n", - " ) 0. mos_i\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", " (h_two i j k a s s') *.\n", " (f_one a j s' +. \n", - " (List.fold_left (fun accu m -> accu +. h_two m a m j s s') 0. mos_i ) +.\n", - " (List.fold_left (fun accu m -> accu +. h_two m a m j s' s') 0. mos_i') ) \n", - " ) 0. mos_i'\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", " in\n", " match phase with\n", @@ -1388,26 +1392,7 @@ " | Phase.Neg -> -. result\n", " \n", "\n", - "let _ = check 11 integral_value 1 1 2 1" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "let ki = det_I.(0)\n", - "let kj = det_I.(10)\n", - "\n", - "\n", - "let _ =\n", - " let alphas =\n", - " generate_alphas ki kj 1 1 2 1\n", - " in\n", - " compute_HaaF ki alphas kj\n", - "\n", - "let _ = integral_value ki kj" + "let _ = check 0 integral_value 1 1 2 1" ] }, { @@ -1430,7 +1415,239 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# 8. 1 2 22 " + "$|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", + "\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", + "\\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", + "\\end{align}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "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 =\n", + " match s with\n", + " | Spin.Alfa -> mos_virt_a\n", + " | Spin.Beta -> mos_virt_b\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", + "\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", + " 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", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "let ki = det_I.(0)\n", + "let kj = det_I.(11)\n", + "\n", + "let alpha = \n", + " generate_alphas ki kj 1 1 2 2\n", + " |> List.filter (fun alpha -> abs_float (compute_HaaF ki [alpha] kj) < 2.0e-05)\n", + " |> List.fold_left (fun (a,mx) alpha ->\n", + " let x = abs_float (compute_HaaF ki [alpha] kj) in\n", + " if x > mx then (alpha, x) else (a,mx) ) (ki,-1.)\n", + " |> fst\n", + " \n", + "\n", + "let _ = Determinant.to_lists ki\n", + "let _ = Determinant.to_lists alpha\n", + "let _ = Determinant.to_lists kj\n", + "\n", + "let _ = \n", + " Format.printf \"Excitation: %a\\n@.\" Excitation.pp_exc (Excitation.of_det ki alpha);\n", + " Format.printf \"Excitation: %a\\n@.\" Excitation.pp_exc (Excitation.of_det alpha kj);\n", + " Format.printf \"Excitation: %a\\n@.\" Excitation.pp_exc (Excitation.of_det ki kj)\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 =\n", + " match s with\n", + " | Spin.Alfa -> mos_virt_a\n", + " | Spin.Beta -> mos_virt_b\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", + "\n", + " let result = \n", + " let s' = Spin.other s in\n", + " let a = 13 in\n", + " (*\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", + " *)\n", + " let a = 16\n", + " and j = 1\n", + " and m = 7\n", + " in\n", + " (*\n", + " h_two i j a m s s' *. f_two a m k j s s'\n", + " *)\n", + " h_two i j m a s s' *. f_two m a k j s s'\n", + " (*\n", + " let a = 18\n", + " and j = 1\n", + " and l = 2\n", + " in\n", + " -. h_two l j k a s s' *. f_two i a l j s s' \n", + " *)\n", + "\n", + " in\n", + " \n", + " match phase with\n", + " | Phase.Pos -> result\n", + " | Phase.Neg -> -. result\n", + " \n", + "\n", + "let _ = compute_HaaF ki [alpha] kj\n", + "let _ = integral_value ki kj\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 8. 1 2 22a" ] }, { @@ -1477,6 +1694,8 @@ " | Spin.Alfa -> mos_a ki, mos_b ki\n", " | Spin.Beta -> mos_b ki, mos_a ki\n", " in\n", + " \n", + " \n", "\n", " let result = \n", " (* Alpha-Beta *)\n",