All F12 diagrams OK

This commit is contained in:
Anthony Scemama 2019-12-18 12:56:35 +01:00
parent 3b22e014b4
commit 3edbf9ca78
1 changed files with 335 additions and 169 deletions

View File

@ -1523,7 +1523,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"# 8. 1 2 22a"
"# 8. 1 2 22"
]
},
{
@ -1599,7 +1599,7 @@
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 50 integral_value 1 2 2 2"
"let _ = check 0 integral_value 1 2 2 2"
]
},
{
@ -1688,13 +1688,345 @@
"let _ = check 0 integral_value 2 1 1 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 10. 2 1 21"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"png_image \"2_1_21.png\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
"\n",
"$m$ : orbital indices of MOs unoccupied both in $|I\\rangle$ and $|J\\rangle$.\n",
"\n",
"$n$ : orbital indices of MOs occupied in $|I\\rangle$\n",
"\n",
"$\\sum_{a} \n",
" \\hat{T}_{ij}^{al} \\hat{T}_{a}^{k} +\n",
" \\hat{T}_{ij}^{ka} \\hat{T}_{a}^{l} +\n",
" \\hat{T}_{i\\bar{j}}^{a\\bar{l}} \\hat{T}_{a}^{k} +\n",
" \\hat{T}_{i\\bar{j}}^{k\\bar{a}} \\hat{T}_{\\bar{a}}^{\\bar{l}}\n",
"$\n",
"\n",
"\n",
"$$\n",
"\\sum_{a}\n",
" \\langle i j || a l \\rangle \\left( f_{ak} + \\sum_{n} \\left[ a n || k n \\right] + \\sum_{\\bar{n}} \\left[ a \\bar{n} | k \\bar{n} \\right] \\right) + \\\\\n",
"\\sum_{a}\n",
" \\langle i j || k a \\rangle \\left( f_{al} + \\sum_{n} \\left[ n a || n l \\right] + \\sum_{\\bar{n}} \\left[ \\bar{n} a | \\bar{n} l \\right] \\right) \n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let integral_value ki kj = \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",
"\n",
" let result = \n",
" \n",
" let mos, mos', s'' =\n",
" match s with\n",
" | Spin.Alfa -> mos_a kj, mos_b kj, Spin.Beta\n",
" | Spin.Beta -> mos_b kj, mos_a kj, Spin.Alfa\n",
" in\n",
" sum mos_cabs (fun a -> \n",
" h_two i j a l s s' *.\n",
" ( f_one a k s +.\n",
" sum mos (fun n -> f_two a n k n s s) +.\n",
" sum mos' (fun n -> f_two a n k n s s'') \n",
" ) ) +.\n",
"\n",
" let mos, mos', s'' =\n",
" match s' with\n",
" | Spin.Alfa -> mos_a kj, mos_b kj, Spin.Beta\n",
" | Spin.Beta -> mos_b kj, mos_a kj, Spin.Alfa\n",
" in\n",
" sum mos_cabs (fun a -> \n",
" h_two i j k a s s' *.\n",
" ( f_one a l s' +.\n",
" sum mos (fun n -> f_two n a n l s' s') +.\n",
" sum mos' (fun n -> f_two n a n l s'' s') \n",
" ) \n",
" ) \n",
" in\n",
"\n",
" match phase with\n",
" | Phase.Pos -> result\n",
" | Phase.Neg -> -. result\n",
"\n",
"\n",
"let _ = check 0 integral_value 2 1 2 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 11. *2 1 22*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"png_image \"2_1_22.png\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
"\n",
"$m$ : orbital indices of MOs unoccupied both in $|I\\rangle$ and $|J\\rangle$.\n",
"\n",
"$n$ : orbital indices of MOs occupied in $|I\\rangle$\n",
"\n",
"$\\sum_{a} \n",
" \\hat{T}_{ij}^{am} \\hat{T}_{am}^{kl} +\n",
" \\hat{T}_{ij}^{ma} \\hat{T}_{ma}^{kl} -\n",
" \\hat{T}_{in}^{al} \\hat{T}_{aj}^{kn} -\n",
" \\hat{T}_{nj}^{ka} \\hat{T}_{ia}^{nl} \n",
"$\n",
"\n",
"\n",
"$$\n",
"\\sum_{a}\n",
"\\sum_{m} \\langle i j || a m \\rangle \\left[ a m || k l \\right] + \\langle i j || m a \\rangle \\left[ m a || k l \\right] - \\\\\n",
"\\sum_{n} \\langle i n || a l \\rangle \\left[ a j || k n \\right] + \\langle n j || k a \\rangle \\left[ i a || n l \\right] \n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let integral_value ki kj = \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",
"\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 result = \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",
" 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",
" 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",
" 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 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",
" in\n",
"\n",
" match phase with\n",
" | Phase.Pos -> result\n",
" | Phase.Neg -> -. result\n",
"\n",
"\n",
"let _ = check 0 integral_value 2 1 2 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 12. 2 2 22"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"png_image \"2_2_22.png\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
"\n",
"$$\n",
"\\sum_{b}\\sum_{a<b}\n",
" \\langle i j || a b \\rangle \n",
" \\left[ ab || kj \\right]\n",
"$$"
]
},
{
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 13. *3 1 2 2*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"png_image \"3_1_22.png\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$|J\\rangle = \\hat{T}_{ijm}^{kln}|I\\rangle$\n",
"\n",
"$\\sum_{a}\n",
" \\hat{T}_{ij}^{al} \\hat{T}_{am}^{kn} \n",
" + \\hat{T}_{ij}^{ka} \\hat{T}_{am}^{ln} \n",
" + \\hat{T}_{im}^{ka} \\hat{T}_{ja}^{ln} \n",
" + \\hat{T}_{im}^{an} \\hat{T}_{aj}^{kl} \n",
" + \\hat{T}_{jm}^{an} \\hat{T}_{ia}^{kl} \n",
" + \\hat{T}_{jm}^{la} \\hat{T}_{ia}^{kn} \n",
"$\n",
"\n",
"$$\n",
"\\sum_{a} \\langle i j || a l \\rangle \\left[ a m || k n \\right] +\n",
" \\langle i j || k a \\rangle \\left[ a m || l n \\right] +\n",
" \\langle i m || k a \\rangle \\left[ j a || l n \\right] +\n",
" \\langle i m || a n \\rangle \\left[ a j || k l \\right] +\n",
" \\langle j m || a n \\rangle \\left[ i a || k l \\right] +\n",
" \\langle j m || l a \\rangle \\left[ i a || k n \\right] \n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let integral_value ki kj = \n",
" let i, j, m, k, l, n, s1, s2, s3, phase =\n",
" match Excitation.of_det ki kj with\n",
" | Excitation.(Triple (phase,\n",
" { hole=h1 ; particle=p1 ; spin=s1 },\n",
" { hole=h2 ; particle=p2 ; spin=s2 },\n",
" { hole=h3 ; particle=p3 ; spin=s3 }) ) -> h1, h2, h3, p1, p2, p3, s1, s2, s3, phase\n",
" | _ -> assert false\n",
" in\n",
"\n",
" let result = \n",
" sum mos_cabs (fun a -> \n",
" h_two i j a l s1 s2 *. f_two a m k n s1 s3 +.\n",
" h_two i j k a s1 s2 *. f_two a m l n s2 s3 +.\n",
" h_two i m k a s1 s3 *. f_two j a l n s2 s3 +.\n",
" h_two i m a n s1 s3 *. f_two a j k l s1 s2 +.\n",
" h_two j m a n s2 s3 *. f_two i a k l s1 s2 +.\n",
" h_two j m l a s2 s3 *. f_two i a k n s1 s3 \n",
" ) \n",
" in\n",
" match phase with\n",
" | Phase.Pos -> result\n",
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 0 integral_value 3 1 2 2"
]
},
{
"cell_type": "raw",
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"let ki = det_I.(0)\n",
"let kj = det_I.(11)\n",
@ -1787,172 +2119,6 @@
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 10. *2 1 21*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"png_image \"2_1_21.png\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 11. *2 1 22*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"png_image \"2_1_22.png\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 12. 2 2 22"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"png_image \"2_2_22.png\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
"\n",
"$$\n",
"\\sum_{b}\\sum_{a<b}\n",
" \\langle i j || a b \\rangle \n",
" \\left[ ab || kj \\right]\n",
"$$"
]
},
{
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 13. *3 1 2 2*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"png_image \"3_1_22.png\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$|J\\rangle = \\hat{T}_{ijm}^{kln}|I\\rangle$\n",
"\n",
"$$\n",
"\\sum_{a} \\langle i j || k a \\rangle \\left[ m a || n l \\right] \n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let integral_value ki kj = \n",
" let h1, h2, h3, p1, p2, p3, s1, s2, s3, phase =\n",
" match Excitation.of_det ki kj with\n",
" | Excitation.(Triple (phase,\n",
" { hole=h1 ; particle=p1 ; spin=s1 },\n",
" { hole=h2 ; particle=p2 ; spin=s2 },\n",
" { hole=h3 ; particle=p3 ; spin=s3 }) ) -> h1, h2, h3, p1, p2, p3, s1, s2, s3, phase\n",
" | _ -> assert false\n",
" in\n",
"\n",
"(*\n",
"123\n",
"213\n",
"231\n",
"321\n",
"312\n",
"132\n",
"*)\n",
"\n",
" let result = \n",
" List.fold_left (fun accu a -> accu +.\n",
" h_two h1 h2 p1 a s1 s2 *. f_two h3 a p3 p2 s3 s2 -.\n",
" h_two h2 h1 p2 a s2 s1 *. f_two h3 a p3 p1 s3 s1 +.\n",
" h_two h2 h3 p2 a s2 s3 *. f_two h1 a p1 p3 s1 s3 -.\n",
" h_two h3 h2 p3 a s3 s2 *. f_two h1 a p1 p2 s1 s2 +.\n",
" h_two h3 h1 p3 a s3 s1 *. f_two h2 a p2 p1 s2 s1 -.\n",
" h_two h1 h3 p1 a s1 s3 *. f_two h2 a p2 p3 s2 s3 \n",
" ) 0. mos_cabs\n",
" in\n",
" match phase with\n",
" | Phase.Pos -> result\n",
" | Phase.Neg -> -. result\n",
" \n",
"\n",
"let _ = check 100 integral_value 3 1 2 2"
]
},
{
"cell_type": "code",
"execution_count": null,