From c0080617ab1f68391efaf633a45ba92d9182e868 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 3 Dec 2019 12:25:31 +0100 Subject: [PATCH] Fixed bytecode C binding --- CI/Determinant.ml | 8 +- CI/Determinant.mli | 2 +- CI/DeterminantSpace.ml | 4 +- CI/DeterminantSpace.mli | 2 +- CI/Phase.ml | 2 +- CI/Phase.mli | 2 +- CI/Spindeterminant.ml | 10 +- CI/Spindeterminant.mli | 2 +- CI/SpindeterminantSpace.ml | 10 +- CI/SpindeterminantSpace.mli | 2 +- MOBasis/MOClass.ml | 11 +- MOBasis/MOClass.mli | 4 +- Notebooks/F12_matrix.ipynb | 1151 +++++++++++++++++++++++++---------- Utils/Bitstring.ml | 17 +- Utils/Electrons.ml | 1 + Utils/Electrons.mli | 3 + Utils/math_functions.c | 6 +- 17 files changed, 901 insertions(+), 336 deletions(-) diff --git a/CI/Determinant.ml b/CI/Determinant.ml index d932db1..2c51d80 100644 --- a/CI/Determinant.ml +++ b/CI/Determinant.ml @@ -90,11 +90,11 @@ let double_excitation spin h p spin' h' p' t = -let pp_det n ppf t = +let pp n ppf t = Format.fprintf ppf "@[@[phase:%a@]@;@[a:%a@]@;@[b:%a@]@]@." - Phase.pp_phase (phase t) - (Spindeterminant.pp_spindet n) t.alfa - (Spindeterminant.pp_spindet n) t.beta + Phase.pp (phase t) + (Spindeterminant.pp n) t.alfa + (Spindeterminant.pp n) t.beta diff --git a/CI/Determinant.mli b/CI/Determinant.mli index aa0836f..eed4e94 100644 --- a/CI/Determinant.mli +++ b/CI/Determinant.mli @@ -76,7 +76,7 @@ val negate_phase : t -> t (** {1 Printers} *) -val pp_det : int -> Format.formatter -> t -> unit +val pp : int -> Format.formatter -> t -> unit (** First [int] is the number of MOs to print. *) (** {1 Unit tests} *) diff --git a/CI/DeterminantSpace.ml b/CI/DeterminantSpace.ml index d9d1e24..e024921 100644 --- a/CI/DeterminantSpace.ml +++ b/CI/DeterminantSpace.ml @@ -348,11 +348,11 @@ let cas_f12_of_mo_basis mo_basis ~frozen_core n m mo_num = -let pp_det_space ppf t = +let pp ppf t = Format.fprintf ppf "@[[ "; let i = ref 0 in determinant_stream t |> Stream.iter (fun d -> Format.fprintf ppf "@[@[%8d@]@;@[%a@]@]@;" !i - (Determinant.pp_det (MOBasis.size (mo_basis t))) d; incr i) ; + (Determinant.pp (MOBasis.size (mo_basis t))) d; incr i) ; Format.fprintf ppf "]@]" diff --git a/CI/DeterminantSpace.mli b/CI/DeterminantSpace.mli index 1c189a3..14afee7 100644 --- a/CI/DeterminantSpace.mli +++ b/CI/DeterminantSpace.mli @@ -71,4 +71,4 @@ val cas_f12_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> int -> (** {2 Printing} *) -val pp_det_space : Format.formatter -> t -> unit +val pp : Format.formatter -> t -> unit diff --git a/CI/Phase.ml b/CI/Phase.ml index cfa9e07..991c7fb 100644 --- a/CI/Phase.ml +++ b/CI/Phase.ml @@ -25,7 +25,7 @@ let add_nperm phase = function | 0 -> phase | nperm -> add phase (of_nperm nperm) -let pp_phase ppf = function +let pp ppf = function | Pos -> Format.fprintf ppf "@[+1@]" | Neg -> Format.fprintf ppf "@[-1@]" diff --git a/CI/Phase.mli b/CI/Phase.mli index f186370..ffa36c1 100644 --- a/CI/Phase.mli +++ b/CI/Phase.mli @@ -19,4 +19,4 @@ val neg : t -> t (** {1 Printers} *) -val pp_phase : Format.formatter -> t -> unit +val pp : Format.formatter -> t -> unit diff --git a/CI/Spindeterminant.ml b/CI/Spindeterminant.ml index 275d869..066f283 100644 --- a/CI/Spindeterminant.ml +++ b/CI/Spindeterminant.ml @@ -128,10 +128,10 @@ let n_electrons = function | None -> 0 -let pp_spindet n ppf = function +let pp n ppf = function | None -> Format.fprintf ppf "@[None@]" | Some s -> - Format.fprintf ppf "@[%a %a@]" Phase.pp_phase s.phase Bitstring.pp + Format.fprintf ppf "@[%a %a@]" Phase.pp s.phase Bitstring.pp s.bitstring @@ -204,9 +204,9 @@ let test_case () = let det = of_list 10 l_a in let l_b = [ 1 ; 7 ; 3 ; 5 ] in let det2 = of_list 10 l_b in - Format.printf "%a@." (pp_spindet 7) det; - Format.printf "%a@." (pp_spindet 7) det2; - Format.printf "%a@." (pp_spindet 7) (single_excitation_reference 2 7 det); + Format.printf "%a@." (pp 7) det; + Format.printf "%a@." (pp 7) det2; + Format.printf "%a@." (pp 7) (single_excitation_reference 2 7 det); Alcotest.(check bool) "single 1" true (single_excitation_reference 2 7 det = det2); Alcotest.(check bool) "single 2" true (single_excitation 2 7 det = single_excitation_reference 2 7 det); Alcotest.(check bool) "single 3" true (single_excitation_reference 4 7 det |> is_none); diff --git a/CI/Spindeterminant.mli b/CI/Spindeterminant.mli index 48651ef..50bb00b 100644 --- a/CI/Spindeterminant.mli +++ b/CI/Spindeterminant.mli @@ -81,7 +81,7 @@ val to_array : t -> int array (** {1 Printers}. *) -val pp_spindet : int -> Format.formatter -> t -> unit +val pp : int -> Format.formatter -> t -> unit (** First [int] is the number of MOs to print *) diff --git a/CI/SpindeterminantSpace.ml b/CI/SpindeterminantSpace.ml index 589687b..06d25fb 100644 --- a/CI/SpindeterminantSpace.ml +++ b/CI/SpindeterminantSpace.ml @@ -61,10 +61,12 @@ let cas_of_mo_basis mo_basis ~frozen_core elec_num n m = -let pp_spindet_space ppf t = - Format.fprintf ppf "@[[ "; - Array.iteri (fun i d -> Format.fprintf ppf "@[@[%8d@] @[%a@]@]@;" i - (Spindeterminant.pp_spindet (MOBasis.size (mo_basis t))) d) (spin_determinants t) ; +let pp ppf t = + Format.fprintf ppf "@[ ["; + let pp = Spindeterminant.pp @@ MOBasis.size (mo_basis t) in + Array.iteri (fun i d -> + Format.fprintf ppf "@[@[%8d@] @[%a@]@]@;" i pp d) + (spin_determinants t) ; Format.fprintf ppf "]@]" diff --git a/CI/SpindeterminantSpace.mli b/CI/SpindeterminantSpace.mli index 555abcd..be955ce 100644 --- a/CI/SpindeterminantSpace.mli +++ b/CI/SpindeterminantSpace.mli @@ -38,7 +38,7 @@ val cas_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> int -> t (** {2 Printing} *) -val pp_spindet_space : Format.formatter -> t -> unit +val pp : Format.formatter -> t -> unit diff --git a/MOBasis/MOClass.ml b/MOBasis/MOClass.ml index 242295d..2276d3f 100644 --- a/MOBasis/MOClass.ml +++ b/MOBasis/MOClass.ml @@ -9,7 +9,7 @@ type mo_class = type t = mo_class list -let pp_mo_occ ppf = function +let pp_mo_class ppf = function | Core i -> Format.fprintf ppf "@[Core %d@]" i | Inactive i -> Format.fprintf ppf "@[Inactive %d@]" i | Active i -> Format.fprintf ppf "@[Active %d@]" i @@ -17,6 +17,15 @@ let pp_mo_occ ppf = function | Deleted i -> Format.fprintf ppf "@[Deleted %d@]" i | Auxiliary i -> Format.fprintf ppf "@[Auxiliary %d@]" i +let pp ppf t = + Format.fprintf ppf "@[[@,"; + let rec aux = function + | [] -> Format.fprintf ppf "]@]" + | x :: [] -> Format.fprintf ppf "%a@,]@]" pp_mo_class x + | x :: rest -> ( Format.fprintf ppf "%a@,;@," pp_mo_class x; aux rest ) + in + aux t + let of_list t = t diff --git a/MOBasis/MOClass.mli b/MOBasis/MOClass.mli index 4e63c64..9f4cd1b 100644 --- a/MOBasis/MOClass.mli +++ b/MOBasis/MOClass.mli @@ -52,5 +52,7 @@ val mo_class_array : t -> mo_class array (** {2 Printers} *) -val pp_mo_occ : Format.formatter -> mo_class -> unit +val pp_mo_class : Format.formatter -> mo_class -> unit + +val pp : Format.formatter -> t -> unit diff --git a/Notebooks/F12_matrix.ipynb b/Notebooks/F12_matrix.ipynb index fef4455..de7ec8d 100644 --- a/Notebooks/F12_matrix.ipynb +++ b/Notebooks/F12_matrix.ipynb @@ -19,72 +19,55 @@ "execution_count": 1, "metadata": {}, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "- : unit = ()\n", - "Findlib has been successfully loaded. Additional directives:\n", - " #require \"package\";; to load a package\n", - " #list;; to list the available packages\n", - " #camlp4o;; to load camlp4 (standard syntax)\n", - " #camlp4r;; to load camlp4 (revised syntax)\n", - " #predicates \"p,q,...\";; to set these predicates\n", - " Topfind.reset();; to force that packages will be reloaded\n", - " #thread;; to enable threads\n", - "\n", - "- : unit = ()\n" - ] - }, { "name": "stderr", "output_type": "stream", "text": [ - "/home/scemama/qp2/external/opam/default/lib/bytes: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/base64: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/base64/base64.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/ocaml/compiler-libs: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/ocaml/compiler-libs/ocamlcommon.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/result: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/result/result.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/ppx_deriving/runtime: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/ppx_deriving/runtime/ppx_deriving_runtime.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/ppx_deriving_yojson/runtime: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/ppx_deriving_yojson/runtime/ppx_deriving_yojson_runtime.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/ocaml/unix.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/uuidm: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/uuidm/uuidm.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/easy-format: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/easy-format/easy_format.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/biniou: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/biniou/biniou.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/yojson: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/yojson/yojson.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/jupyter: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/jupyter/jupyter.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/jupyter/notebook: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/jupyter/notebook/jupyter_notebook.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/ocaml/bigarray.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/lacaml: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/lacaml/lacaml.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/astring: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/astring/astring.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/cmdliner: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/cmdliner/cmdliner.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/seq: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/stdlib-shims: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/stdlib-shims/stdlib_shims.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/fmt: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/fmt/fmt.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/fmt/fmt_cli.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/fmt/fmt_tty.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/alcotest: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/alcotest/alcotest.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/ocaml/str.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/zarith: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/zarith/zarith.cma: loaded\n", - "/home/scemama/qp2/external/opam/default/lib/getopt: added to search path\n", - "/home/scemama/qp2/external/opam/default/lib/getopt/getopt.cma: loaded\n" + "/home/scemama/qp2/external/opam/4.07.1/lib/bytes: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/base64: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/base64/base64.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/compiler-libs: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/compiler-libs/ocamlcommon.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/result: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/result/result.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/ppx_deriving/runtime: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/ppx_deriving/runtime/ppx_deriving_runtime.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/ppx_deriving_yojson/runtime: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/ppx_deriving_yojson/runtime/ppx_deriving_yojson_runtime.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/unix.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/uuidm: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/uuidm/uuidm.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/easy-format: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/easy-format/easy_format.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/biniou: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/biniou/biniou.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/yojson: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/yojson/yojson.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/jupyter: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/jupyter/jupyter.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/jupyter/notebook: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/jupyter/notebook/jupyter_notebook.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/bigarray.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/lacaml: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/lacaml/lacaml.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/astring: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/astring/astring.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/cmdliner: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/cmdliner/cmdliner.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/seq: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/stdlib-shims: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/stdlib-shims/stdlib_shims.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/fmt: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/fmt/fmt.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/fmt/fmt_cli.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/fmt/fmt_tty.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/alcotest: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/alcotest/alcotest.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/str.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/zarith: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/zarith/zarith.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/getopt: added to search path\n", + "/home/scemama/qp2/external/opam/4.07.1/lib/getopt/getopt.cma: loaded\n" ] } ], @@ -127,7 +110,8 @@ "#load_rec \"Util.cma\";;\n", "#load_rec \"Matrix.cmo\";;\n", "#load_rec \"Simulation.cmo\";;\n", - "#load_rec \"Simulation.cmo\";;\n", + "#load_rec \"Spindeterminant.cmo\";;\n", + "#load_rec \"Determinant.cmo\";;\n", "#load_rec \"HartreeFock.cmo\";;\n", "#load_rec \"MOBasis.cmo\";;\n", "#load_rec \"F12CI.cmo\";;\n", @@ -150,7 +134,7 @@ { "data": { "text/plain": [ - "val f : Format.formatter -> MOBasis.t -> unit = \n" + "val pp_mo : Format.formatter -> MOBasis.t -> unit = \n" ] }, "execution_count": 3, @@ -170,8 +154,12 @@ "#install_printer Util.pp_matrix;;\n", "#install_printer HartreeFock.pp ;;\n", "#install_printer Fock.pp ;;\n", - "let f ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;\n", - "#install_printer f;;\n" + "#install_printer MOClass.pp ;;\n", + "let pp_mo ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;\n", + "#install_printer pp_mo;;\n", + "#install_printer DeterminantSpace.pp;;\n", + "#install_printer SpindeterminantSpace.pp;;\n", + "#install_printer Bitstring.pp;;\n" ] }, { @@ -196,7 +184,7 @@ { "data": { "text/plain": [ - "val basis_filename : string = \"/home/scemama/qp2/data/basis/cc-pvdz\"\n" + "val basis_filename : string = \"/home/scemama/qp2/data/basis/6-31g\"\n" ] }, "execution_count": 4, @@ -206,7 +194,7 @@ { "data": { "text/plain": [ - "val aux_basis_filename : string = \"/home/scemama/qp2/data/basis/cc-pvtz\"\n" + "val aux_basis_filename : string = \"/home/scemama/qp2/data/basis/cc-pvdz\"\n" ] }, "execution_count": 4, @@ -223,6 +211,16 @@ "metadata": {}, "output_type": "execute_result" }, + { + "data": { + "text/plain": [ + "val frozen_core : bool = false\n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, { "data": { "text/plain": [ @@ -255,37 +253,26 @@ " # Angular Coordinates (Bohr) Exponents Coefficients\n", " Momentum X Y Z\n", "-----------------------------------------------------------------------\n", - " 1-3 S 0.0000 0.0000 0.0000 2.94000000e+03 6.80000000e-04\n", - " 4.41200000e+02 5.23600000e-03\n", - " 1.00500000e+02 2.66060000e-02\n", - " 2.84300000e+01 9.99930000e-02\n", - " 9.16900000e+00 2.69702000e-01\n", - " 3.19600000e+00 4.51469000e-01\n", - " 1.15900000e+00 2.95074000e-01\n", - " 1.81100000e-01 1.25870000e-02\n", + " 1-3 S 0.0000 0.0000 0.0000 1.26458570e+03 1.94480000e-03\n", + " 1.89936810e+02 1.48351000e-02\n", + " 4.31590890e+01 7.20906000e-02\n", + " 1.20986630e+01 2.37154200e-01\n", + " 3.80632320e+00 4.69198700e-01\n", + " 1.27289030e+00 3.56520200e-01\n", " \n", - " 2.94000000e+03 -1.23000000e-04\n", - " 4.41200000e+02 -9.66000000e-04\n", - " 1.00500000e+02 -4.83100000e-03\n", - " 2.84300000e+01 -1.93140000e-02\n", - " 9.16900000e+00 -5.32800000e-02\n", - " 3.19600000e+00 -1.20723000e-01\n", - " 1.15900000e+00 -1.33435000e-01\n", - " 1.81100000e-01 5.30767000e-01\n", + " 3.19646310e+00 -1.12648700e-01\n", + " 7.47813300e-01 -2.29506400e-01\n", + " 2.19966300e-01 1.18691670e+00\n", " \n", - " 5.89000000e-02 1.00000000e+00\n", + " 8.23099000e-02 1.00000000e+00\n", " \n", " \n", "-----------------------------------------------------------------------\n", - " 4-9 P 0.0000 0.0000 0.0000 3.61900000e+00 2.91110000e-02\n", - " 7.11000000e-01 1.69365000e-01\n", - " 1.95100000e-01 5.13458000e-01\n", + " 3-8 P 0.0000 0.0000 0.0000 3.19646310e+00 5.59802000e-02\n", + " 7.47813300e-01 2.61550600e-01\n", + " 2.19966300e-01 7.93972300e-01\n", " \n", - " 6.01800000e-02 1.00000000e+00\n", - " \n", - " \n", - "-----------------------------------------------------------------------\n", - " 10-15 D 0.0000 0.0000 0.0000 2.38000000e-01 1.00000000e+00\n", + " 8.23099000e-02 1.00000000e+00\n", " \n", " \n", "-----------------------------------------------------------------------\n", @@ -308,7 +295,20 @@ " # Angular Coordinates (Bohr) Exponents Coefficients\n", " Momentum X Y Z\n", "-----------------------------------------------------------------------\n", - " 1-7 S 0.0000 0.0000 0.0000 2.94000000e+03 6.80000000e-04\n", + " 1-6 S 0.0000 0.0000 0.0000 1.26458570e+03 1.94480000e-03\n", + " 1.89936810e+02 1.48351000e-02\n", + " 4.31590890e+01 7.20906000e-02\n", + " 1.20986630e+01 2.37154200e-01\n", + " 3.80632320e+00 4.69198700e-01\n", + " 1.27289030e+00 3.56520200e-01\n", + " \n", + " 3.19646310e+00 -1.12648700e-01\n", + " 7.47813300e-01 -2.29506400e-01\n", + " 2.19966300e-01 1.18691670e+00\n", + " \n", + " 8.23099000e-02 1.00000000e+00\n", + " \n", + " 2.94000000e+03 6.80000000e-04\n", " 4.41200000e+02 5.23600000e-03\n", " 1.00500000e+02 2.66060000e-02\n", " 2.84300000e+01 9.99930000e-02\n", @@ -328,57 +328,23 @@ " \n", " 5.89000000e-02 1.00000000e+00\n", " \n", - " 6.86300000e+03 2.36000000e-04\n", - " 1.03000000e+03 1.82600000e-03\n", - " 2.34700000e+02 9.45200000e-03\n", - " 6.65600000e+01 3.79570000e-02\n", - " 2.16900000e+01 1.19965000e-01\n", - " 7.73400000e+00 2.82162000e-01\n", - " 2.91600000e+00 4.27404000e-01\n", - " 1.13000000e+00 2.66278000e-01\n", - " 1.10100000e-01 -7.27500000e-03\n", - " \n", - " 6.86300000e+03 -4.30000000e-05\n", - " 1.03000000e+03 -3.33000000e-04\n", - " 2.34700000e+02 -1.73600000e-03\n", - " 6.65600000e+01 -7.01200000e-03\n", - " 2.16900000e+01 -2.31260000e-02\n", - " 7.73400000e+00 -5.81380000e-02\n", - " 2.91600000e+00 -1.14556000e-01\n", - " 1.13000000e+00 -1.35908000e-01\n", - " 1.10100000e-01 5.77441000e-01\n", - " \n", - " 2.57700000e-01 1.00000000e+00\n", - " \n", - " 4.40900000e-02 1.00000000e+00\n", - " \n", " \n", "-----------------------------------------------------------------------\n", - " 4-18 P 0.0000 0.0000 0.0000 3.61900000e+00 2.91110000e-02\n", + " 3-14 P 0.0000 0.0000 0.0000 3.19646310e+00 5.59802000e-02\n", + " 7.47813300e-01 2.61550600e-01\n", + " 2.19966300e-01 7.93972300e-01\n", + " \n", + " 8.23099000e-02 1.00000000e+00\n", + " \n", + " 3.61900000e+00 2.91110000e-02\n", " 7.11000000e-01 1.69365000e-01\n", " 1.95100000e-01 5.13458000e-01\n", " \n", " 6.01800000e-02 1.00000000e+00\n", " \n", - " 7.43600000e+00 1.07360000e-02\n", - " 1.57700000e+00 6.28540000e-02\n", - " 4.35200000e-01 2.48180000e-01\n", - " \n", - " 1.43800000e-01 1.00000000e+00\n", - " \n", - " 4.99400000e-02 1.00000000e+00\n", - " \n", " \n", "-----------------------------------------------------------------------\n", - " 10-27 D 0.0000 0.0000 0.0000 2.38000000e-01 1.00000000e+00\n", - " \n", - " 3.48000000e-01 1.00000000e+00\n", - " \n", - " 1.80300000e-01 1.00000000e+00\n", - " \n", - " \n", - "-----------------------------------------------------------------------\n", - " 41-50 F 0.0000 0.0000 0.0000 3.25000000e-01 1.00000000e+00\n", + " 19-24 D 0.0000 0.0000 0.0000 2.38000000e-01 1.00000000e+00\n", " \n", " \n", "-----------------------------------------------------------------------\n", @@ -429,9 +395,10 @@ } ], "source": [ - "let basis_filename = \"/home/scemama/qp2/data/basis/cc-pvdz\" \n", - "let aux_basis_filename = \"/home/scemama/qp2/data/basis/cc-pvtz\" \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 \"be\" \n", + "let frozen_core = false\n", "let multiplicity = 1\n", "let state = 1\n", "\n", @@ -464,44 +431,33 @@ "name": "stdout", "output_type": "stream", "text": [ - "21 significant shell pairs computed in 0.055000 seconds\n", + "15 significant shell pairs computed in 0.004974 seconds\n", "1\n", "2\n", - "3\n", + "5\n", "6\n", - "9\n", - "Computed ERIs in 0.238291 seconds\n", + "Computed ERIs in 0.026928 seconds\n", "MOs =\n", "\n", "\n", - " -- 1 -- -- 2 -- -- 3 -- -- 4 -- -- 5 --\n", - " 1 1.0006 0.00028918 0 0 0\n", - " 2 0.00576322 1.00367 0 1.02041E-15 1.08506E-15\n", - " 3 -0.00181801 0.582215 0 -8.39813E-16 -8.40021E-16\n", - " ... ... ... ... ...\n", - " 13 -0.00209072 -0.00172308 9.60612E-16 0 0\n", - " 14 0 -0 0 0 0\n", - " 15 -0.00209072 -0.00172308 -1.26705E-15 0 4.17104E-16\n", + " -- 1 -- -- 2 -- -- 3 -- -- 4 -- -- 5 --\n", + " 1 0.997996 -0.22266 0 0 -0\n", + " 2 0.0161427 0.288767 2.2404E-15 0 -7.73804E-16\n", + " 3 0 3.58344E-16 0.258415 0.037609 0.0381979\n", + " ... ... ... ... ...\n", + " 7 0 0 0.78727 0.114442 0.116234\n", + " 8 0 0 -0.128237 0.0798944 0.789475\n", + " 9 0 0 -0.100874 0.791587 -0.0964853\n", " \n", "\n", - " -- 6 -- -- 7 -- -- 8 -- -- 9 -- -- 10 --\n", - " 1 0.549926 0 -2.88529E-15 -0 8.58482E-16\n", - " 2 2.77214 0 -1.582E-14 -0 3.81097E-15\n", - " 3 -1.826 0 9.32169E-15 -0 -3.39771E-15\n", - " ... ... ... ... ...\n", - " 13 0.224421 0 -1.20926E-15 -0 -0.0220572\n", - " 14 0 0 -0 -1.20416E-15 0.0718784\n", - " 15 0.224421 0 -0 -0 0.847943\n", - " \n", - "\n", - " -- 11 -- -- 12 -- -- 13 -- -- 14 -- -- 15 --\n", - " 1 -1.28321E-15 0 0 0 0.516133\n", - " 2 -1.06861E-14 0 1.77725E-15 0 4.52711\n", - " 3 -1.84016E-15 0 1.31997E-15 0 1.08073\n", - " ... ... ... ... ...\n", - " 13 -0.964125 -0.0163169 0.262399 0.0292786 -1.45115\n", - " 14 0.0475417 0.0403496 0.0724528 0.992822 -0\n", - " 15 0.406875 0.061849 -0.328775 -0.0593936 -1.45115\n", + " -- 6 -- -- 7 -- -- 8 -- -- 9 --\n", + " 1 0.00338551 0 -0 0\n", + " 2 2.01724 -3.80748E-15 4.70544E-15 -6.99768E-16\n", + " 3 1.56069E-15 1.292 0.249105 0.0981672\n", + " ... ... ... ...\n", + " 7 -2.3622E-15 -1.05652 -0.20374 -0.0802898\n", + " 8 1.64398E-15 0.172096 -0.528064 -0.92522\n", + " 9 2.61906E-15 0.135375 -0.918775 0.549573\n", " \n" ] }, @@ -516,35 +472,33 @@ " ------------------------------------------------------------\n", " # HF energy Convergence HOMO-LUMO\n", " ------------------------------------------------------------\n", - " 1 -14.05142130 8.0958e-01 0.2991\n", - " 2 -14.51982122 2.7933e-01 0.3600\n", - " 3 -14.56519147 8.3990e-02 0.3674\n", - " 4 -14.57117146 2.6217e-02 0.3679\n", - " 5 -14.57214606 1.0744e-02 0.3676\n", - " 6 -14.57232209 5.0392e-03 0.3675\n", - " 7 -14.57235521 2.2956e-03 0.3674\n", - " 8 -14.57236154 1.0293e-03 0.3673\n", - " 9 -14.57236304 1.3825e-08 0.3673\n", - " 10 -14.57236304 1.6054e-09 0.3673\n", + " 1 -14.35428398 3.4264e-01 0.3789\n", + " 2 -14.52767889 1.8460e-01 0.3848\n", + " 3 -14.56058094 7.7052e-02 0.3856\n", + " 4 -14.56582672 3.0375e-02 0.3848\n", + " 5 -14.56662378 1.1794e-02 0.3842\n", + " 6 -14.56676403 1.7857e-06 0.3837\n", + " 7 -14.56676404 1.7598e-08 0.3837\n", + " 8 -14.56676403 1.3644e-09 0.3837\n", " ------------------------------------------------------------\n", "\n", "\n", " ============================================================\n", - " One-electron energy -19.0610711219\n", - " Kinetic 14.5598880943\n", - " Potential -33.6209592162\n", + " One-electron energy -19.1171562840\n", + " Kinetic 14.6787428206\n", + " Potential -33.7958991046\n", " -------------------------------------------------------- \n", - " Two-electron energy 4.4887080836\n", - " Coulomb 7.1543328124\n", - " Exchange -2.6656247288\n", + " Two-electron energy 4.5503922505\n", + " Coulomb 7.2329463384\n", + " Exchange -2.6825540879\n", " -------------------------------------------------------- \n", - " HF HOMO -8.4100296256\n", - " HF LUMO 1.5842692516\n", - " HF LUMO-LUMO 9.9942988772\n", + " HF HOMO -8.1986650549\n", + " HF LUMO 2.2431802579\n", + " HF LUMO-LUMO 10.4418453128\n", " -------------------------------------------------------- \n", - " Electronic energy -14.5723630383\n", + " Electronic energy -14.5667640335\n", " Nuclear repulsion 0.0000000000\n", - " Hartree-Fock energy -14.5723630383\n", + " Hartree-Fock energy -14.5667640335\n", " ============================================================\n", " \n", "\n" @@ -558,35 +512,25 @@ "data": { "text/plain": [ "val mo_basis : MOBasis.t =\n", - " Eigenvalues: -4.732765 -0.309063 0.058221 0.058221 0.058221 \n", - " -- 1 -- -- 2 -- -- 3 -- -- 4 -- -- 5 --\n", - " 1 1.0006 0.00028918 0 0 0\n", - " 2 0.00576322 1.00367 0 1.02041E-15 1.08506E-15\n", - " 3 -0.00181801 0.582215 0 -8.39813E-16 -8.40021E-16\n", - " ... ... ... ... ...\n", - " 13 -0.00209072 -0.00172308 9.60612E-16 0 0\n", - " 14 0 -0 0 0 0\n", - " 15 -0.00209072 -0.00172308 -1.26705E-15 0 4.17104E-16\n", + " Eigenvalues: -4.706891 -0.301295 0.082435 0.082435 0.082435 \n", + " -- 1 -- -- 2 -- -- 3 -- -- 4 -- -- 5 --\n", + " 1 0.997996 -0.22266 0 0 -0\n", + " 2 0.0161427 0.288767 2.2404E-15 0 -7.73804E-16\n", + " 3 0 3.58344E-16 0.258415 0.037609 0.0381979\n", + " ... ... ... ... ...\n", + " 7 0 0 0.78727 0.114442 0.116234\n", + " 8 0 0 -0.128237 0.0798944 0.789475\n", + " 9 0 0 -0.100874 0.791587 -0.0964853\n", " \n", - " Eigenvalues: 0.277480 0.350137 0.350137 0.350137 0.650740 \n", - " -- 6 -- -- 7 -- -- 8 -- -- 9 -- -- 10 --\n", - " 1 0.549926 0 -2.88529E-15 -0 8.58482E-16\n", - " 2 2.77214 0 -1.582E-14 -0 3.81097E-15\n", - " 3 -1.826 0 9.32169E-15 -0 -3.39771E-15\n", - " ... ... ... ... ...\n", - " 13 0.224421 0 -1.20926E-15 -0 -0.0220572\n", - " 14 0 0 -0 -1.20416E-15 0.0718784\n", - " 15 0.224421 0 -0 -0 0.847943\n", - " \n", - " Eigenvalues: 0.650740 0.650740 0.650740 0.650740 1.188790 \n", - " -- 11 -- -- 12 -- -- 13 -- -- 14 -- -- 15 --\n", - " 1 -1.28321E-15 0 0 0 0.516133\n", - " 2 -1.06861E-14 0 1.77725E-15 0 4.52711\n", - " 3 -1.84016E-15 0 1.31997E-15 0 1.08073\n", - " ... ... ... ... ...\n", - " 13 -0.964125 -0.0163169 0.262399 0.0292786 -1.45115\n", - " 14 0.0475417 0.0403496 0.0724528 0.992822 -0\n", - " 15 0.406875 0.061849 -0.328775 -0.0593936 -1.45115\n", + " Eigenvalues: 0.439754 0.464931 0.464931 0.464931 \n", + " -- 6 -- -- 7 -- -- 8 -- -- 9 --\n", + " 1 0.00338551 0 -0 0\n", + " 2 2.01724 -3.80748E-15 4.70544E-15 -6.99768E-16\n", + " 3 1.56069E-15 1.292 0.249105 0.0981672\n", + " ... ... ... ...\n", + " 7 -2.3622E-15 -1.05652 -0.20374 -0.0802898\n", + " 8 1.64398E-15 0.172096 -0.528064 -0.92522\n", + " 9 2.61906E-15 0.135375 -0.918775 0.549573\n", " \n", " \n" ] @@ -606,104 +550,153 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### FCI-F12" + "# FCI-F12" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notations:\n", + "\n", + "* $\\langle ij || kl \\rangle = \\int \\phi_i(r_1) \\phi_j(r_2) \\frac{1}{r_{12}} \\phi_k(r1) \\phi_l(r2) $ \n", + "* $\\left[ ij || kl \\right] = \\int \\phi_i(r_1) \\phi_j(r_2) f_{12} \\phi_k(r1) \\phi_l(r2) $ \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Common functions" ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 6, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Number of determinants : 105 105 11025\n" - ] + "data": { + "text/plain": [ + "val f12 : F12factor.t =\n", + " {F12factor.expo_s = 1.;\n", + " gaussian =\n", + " {GaussianOperator.coef_g =\n", + " [ -0.314400 -0.303700 -0.168100 -0.098110 -0.060240 -0.037260 ];\n", + " expo_sg =\n", + " [ 0.220900 1.004000 3.622000 12.160000 45.870000 254.400000 ];\n", + " expo_sg_inv =\n", + " [ 4.526935 0.996016 0.276091 0.082237 0.021801 0.003931 ]}}\n" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" }, { - "name": "stderr", - "output_type": "stream", - "text": [ - "4-idx transformation \n", - "15 / 15\n" - ] + "data": { + "text/plain": [ + "val mo_num : int = 9\n" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "136 significant shell pairs computed in 2.999540 seconds\n", - "1\n", - "2\n", - "3\n", - "6\n", - "9\n", - "15\n", - "16\n", - "17\n", - "18\n", - "19\n", - "22\n", - "25\n", - "28\n", - "34\n", - "40\n", - "Computed ERIs in 110.013171 seconds\n" - ] + "data": { + "text/plain": [ + "val pp_spindet : Format.formatter -> Spindeterminant.t -> unit = \n" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" }, { - "name": "stderr", - "output_type": "stream", - "text": [ - "4-idx transformation \n", - "44 / 44\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "136 significant shell pairs computed in 0.279177 seconds\n", - "1\n", - "2\n", - "3\n", - "6\n", - "9\n", - "15\n", - "16\n", - "17\n", - "18\n", - "19\n", - "22\n", - "25\n", - "28\n", - "34\n", - "40\n", - "Computed ERIs in 8.198963 seconds\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "4-idx transformation \n", - "44 / 44\n", - "Computing HF12 integrals\n" - ] - }, - { - "ename": "error", - "evalue": "runtime_error", - "output_type": "error", - "traceback": [ - "\u001b[31mException: Invalid_argument \"index out of bounds\".\nRaised by primitive operation at unknown location\nCalled from file \"CI/CIMatrixElement.ml\", line 37, characters 47-66\nCalled from file \"list.ml\", line 117, characters 24-34\nCalled from file \"CI/CIMatrixElement.ml\", line 37, characters 8-76\nCalled from file \"list.ml\", line 88, characters 20-23\nCalled from file \"CI/CI.ml\", line 64, characters 2-54\nCalled from file \"CI/CI.ml\", line 641, characters 12-44\nCalled from file \"src/vec4_D.ml\", line 57, characters 29-32\nCalled from file \"camlinternalLazy.ml\", line 27, characters 17-27\nRe-raised at file \"camlinternalLazy.ml\", line 34, characters 10-11\nCalled from file \"CI/CI.ml\", line 687, characters 8-112\nCalled from file \"camlinternalLazy.ml\", line 27, characters 17-27\nRe-raised at file \"camlinternalLazy.ml\", line 34, characters 10-11\nCalled from unknown location\nCalled from file \"[21]\", line 2, characters 4-85\nCalled from file \"toplevel/toploop.ml\", line 180, characters 17-56\n\u001b[0m" - ] + "data": { + "text/plain": [ + "val pp_det : Format.formatter -> Determinant.t -> unit = \n" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "let fcif12 =\n", - " F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ~state ()\n" + "let f12 = Util.of_some @@ Simulation.f12 simulation \n", + "\n", + "let mo_num = MOBasis.size mo_basis \n", + "\n", + "let pp_spindet = Spindeterminant.pp mo_num\n", + "\n", + "let pp_det = Determinant.pp mo_num\n", + "\n", + ";;\n", + "\n", + "#install_printer pp_spindet ;;\n", + "#install_printer pp_det ;;\n", + "\n" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "let simulation_aux = \n", + " let charge = Charge.to_int @@ Simulation.charge simulation \n", + " and multiplicity = Electrons.multiplicity @@ Simulation.electrons simulation\n", + " and nuclei = Simulation.nuclei simulation\n", + " in\n", + " let general_basis = \n", + " Basis.general_basis @@ Simulation.basis simulation\n", + " in\n", + " GeneralBasis.combine [\n", + " general_basis ; GeneralBasis.read aux_basis_filename\n", + " ]\n", + " |> Basis.of_nuclei_and_general_basis nuclei\n", + " |> Simulation.make ~f12 ~charge ~multiplicity ~nuclei \n", + "\n", + "\n", + "let aux_basis = \n", + " MOBasis.of_mo_basis simulation_aux mo_basis\n", + " \n", + "let () = ignore @@ MOBasis.f12_ints aux_basis\n", + "let () = ignore @@ MOBasis.two_e_ints aux_basis\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Determinant-based functions" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Integrals" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$\\langle I | \\hat{H} | J \\rangle = \\begin{cases}\n", + "\\sum_i h_{ii} + \\frac{1}{2} \\sum_{ij} \\langle ij || ij \\rangle \\text{ when } |J\\rangle = |I\\rangle \\\\\n", + "h_{ik} + \\sum_{j} \\langle ij || kj \\rangle \\text{ when } |J\\rangle = \\hat{T}_i^k |I\\rangle \\\\\n", + "\\langle ij || kl \\rangle \\text{ when } |J\\rangle = \\hat{T}_{ij}^{kl} |I\\rangle \\\\\n", + "\\end{cases}$\n", + "\n", + "\n", + "$\\langle I | \\hat{F} | J \\rangle = \\begin{cases}\n", + "\\sum_i f_{ii} + \\frac{1}{2} \\sum_{ij} \\langle ij || ij \\rangle \\text{ when } |J\\rangle = |I\\rangle \\\\\n", + "f_{ik} + \\sum_{j} \\langle ij || kj \\rangle \\text{ when } |J\\rangle = \\hat{T}_i^k |I\\rangle \\\\\n", + "\\langle ij || kl \\rangle \\text{ when } |J\\rangle = \\hat{T}_{ij}^{kl} |I\\rangle \\\\\n", + "\\end{cases}$" ] }, { @@ -712,16 +705,560 @@ "metadata": {}, "outputs": [ { - "ename": "error", - "evalue": "compile_error", - "output_type": "error", - "traceback": [ - "\u001b[32mFile \"[7]\", line 1, characters 0-1:\n\u001b[31mError: Syntax error\n\u001b[36m 1: \u001b[30m\u001b[4m%\u001b[0m\u001b[30mload_ext itikz\u001b[0m\n" - ] + "data": { + "text/plain": [ + "val cancel_singles : bool = false\n" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val f12_integrals :\n", + " MOBasis.t ->\n", + " ('a -> 'b -> 'c -> float) *\n", + " (int -> int -> int -> int -> 'd -> 'd -> float) * 'e option = \n" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val h_ij : MOBasis.t -> CIMatrixElement.De.t -> CIMatrixElement.De.t -> float =\n", + " \n" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val f_ij : MOBasis.t -> CIMatrixElement.De.t -> CIMatrixElement.De.t -> float =\n", + " \n" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val hf_ij :\n", + " MOBasis.t -> CIMatrixElement.De.t -> CIMatrixElement.De.t -> float list =\n", + " \n" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "%load_ext itikz" + "let cancel_singles = false \n", + "\n", + "let f12_integrals mo_basis =\n", + "\n", + " \n", + " let two_e_ints = MOBasis.f12_ints mo_basis in\n", + " let f2 i j k l s s' =\n", + " let ijkl = F12.get_phys two_e_ints i j k l in\n", + " let ijlk = F12.get_phys two_e_ints i j l k in\n", + " if s' <> s then\n", + " ijkl \n", + " else\n", + " (ijkl -. ijlk) \n", + " in\n", + " ( (fun _ _ _ -> 0.),\n", + " (if cancel_singles then\n", + " (fun i j k l s s' ->\n", + " (* Required to cancel out single excitations *)\n", + " if (i=k && j<>l) || (j=l && i<>k) then\n", + " 0.\n", + " else\n", + " f2 i j k l s s'\n", + " ) \n", + " else f2),\n", + " None\n", + " )\n", + "\n", + "let h_ij mo_basis ki kj =\n", + " let integrals =\n", + " List.map (fun f -> f mo_basis)\n", + " [ CI.h_integrals ]\n", + " in\n", + " CIMatrixElement.make integrals ki kj \n", + " |> List.hd\n", + "\n", + "\n", + "let f_ij mo_basis ki kj =\n", + " let integrals =\n", + " List.map (fun f -> f mo_basis)\n", + " [ f12_integrals ]\n", + " in\n", + " CIMatrixElement.make integrals ki kj \n", + " |> List.hd \n", + "\n", + "\n", + "let hf_ij mo_basis ki kj =\n", + " let integrals =\n", + " List.map (fun f -> f mo_basis)\n", + " [ CI.h_integrals ; f12_integrals ]\n", + " in\n", + " CIMatrixElement.make integrals ki kj\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Determinant space" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "- : Spindeterminant.t =\n", + "+1 +-+-+-+---------------------------------------------------------\n" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Spindeterminant.of_list 10 [1; 3 ; 5; 7]" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "val det_space : SpindeterminantSpace.t =\n", + " [ 0 +1 ++--------------------------------------------------------------\n", + " 1 +1 --+-------------------------------------------------------------\n", + " 2 +1 ---+------------------------------------------------------------\n", + " 3 +1 ----+-----------------------------------------------------------\n", + " 4 +1 +----+----------------------------------------------------------\n", + " 5 +1 -+---+----------------------------------------------------------\n", + " 6 +1 --+--+----------------------------------------------------------\n", + " 7 +1 ---+-+----------------------------------------------------------\n", + " 8 +1 ----++----------------------------------------------------------\n", + " 9 +1 ------+---------------------------------------------------------\n", + " 10 +1 +++----+--------------------------------------------------------\n", + " 11 +1 ---+---+--------------------------------------------------------\n", + " 12 +1 ----+--+--------------------------------------------------------\n", + " 13 +1 -----+-+--------------------------------------------------------\n", + " 14 +1 ++----++--------------------------------------------------------\n", + " 15 +1 --+---++--------------------------------------------------------\n", + " 16 +1 ---+--++--------------------------------------------------------\n", + " 17 +1 ----+-++--------------------------------------------------------\n", + " 18 +1 +----+++--------------------------------------------------------\n", + " 19 +1 -+---+++--------------------------------------------------------\n", + " 20 +1 --+--+++--------------------------------------------------------\n", + " 21 +1 ---+-+++--------------------------------------------------------\n", + " 22 +1 ----++++--------------------------------------------------------\n", + " 23 +1 +++-----+-------------------------------------------------------\n", + " 24 +1 ---+----+-------------------------------------------------------\n", + " 25 +1 ----+---+-------------------------------------------------------\n", + " 26 +1 +----+--+-------------------------------------------------------\n", + " 27 +1 -+---+--+-------------------------------------------------------\n", + " 28 +1 --+--+--+-------------------------------------------------------\n", + " 29 +1 ---+-+--+-------------------------------------------------------\n", + " 30 +1 ----++--+-------------------------------------------------------\n", + " 31 +1 +-----+-+-------------------------------------------------------\n", + " 32 +1 -+----+-+-------------------------------------------------------\n", + " 33 +1 --+---+-+-------------------------------------------------------\n", + " 34 +1 ---+--+-+-------------------------------------------------------\n", + " 35 +1 ----+-+-+-------------------------------------------------------\n", + " ]\n" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "let det_space = \n", + "(*\n", + " DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num\n", + " *)\n", + " SpindeterminantSpace.fci_of_mo_basis mo_basis ~frozen_core (Electrons.n_alfa @@ Simulation.electrons simulation)\n" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "unsigned int v; // current permutation of bits \n", + "unsigned int w; // next permutation of bits\n", + "\n", + "unsigned int t = v | (v - 1); // t gets v's least significant 0 bits set to 1\n", + "// Next set to 1 the most significant bit to change, \n", + "// set to 0 the least significant ones, and add the necessary 1 bits.\n", + "w = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1)); " + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "- : int = 31\n" + ] + }, + "execution_count": 51, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "let _ = \n", + " let open Bitstring in\n", + "(*\n", + "\n", + " let u = Bitstring.of_int 63 in\n", + "\n", + " let t = logor u @@ minus_one u in\n", + " let not_t = lognot t in\n", + " let neg_not_t = neg not_t in\n", + " \n", + " shift_right\n", + " (minus_one @@ logand not_t neg_not_t) \n", + " ( (trailing_zeros u) + 1),\n", + " *)\n", + " popcount (Bitstring.of_int 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "val elec_num : int = 4\n" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val alfa_num : int = 2\n" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val mo_class : MOClass.t =\n", + " [Active 1;Active 2;Active 3;Active 4;Active 5;Active 6;Active 7;Active 8;\n", + " Active 9]\n" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val m : int list -> Bitstring.t = \n" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val occ_mask : Bitstring.t =\n", + " ----------------------------------------------------------------\n" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val active_mask : Bitstring.t =\n", + " +++++++++-------------------------------------------------------\n" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val neg_active_mask : Bitstring.t =\n", + " ---------+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val permtutations : int -> int -> Bitstring.t list = \n" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val spin_determinants : Bitstring.t list =\n", + " [++--------------------------------------------------------------;\n", + " --+-------------------------------------------------------------;\n", + " ---+------------------------------------------------------------;\n", + " ----+-----------------------------------------------------------;\n", + " -----+----------------------------------------------------------;\n", + " ++----+---------------------------------------------------------;\n", + " --+---+---------------------------------------------------------;\n", + " ---+--+---------------------------------------------------------;\n", + " ----+-+---------------------------------------------------------;\n", + " +----++---------------------------------------------------------;\n", + " -+---++---------------------------------------------------------;\n", + " --+--++---------------------------------------------------------;\n", + " ---+-++---------------------------------------------------------;\n", + " ----+++---------------------------------------------------------;\n", + " +------+--------------------------------------------------------;\n", + " -+-----+--------------------------------------------------------;\n", + " --+----+--------------------------------------------------------;\n", + " ---+---+--------------------------------------------------------;\n", + " ----+--+--------------------------------------------------------;\n", + " +----+-+--------------------------------------------------------;\n", + " -+---+-+--------------------------------------------------------;\n", + " --+--+-+--------------------------------------------------------;\n", + " ---+-+-+--------------------------------------------------------;\n", + " ----++-+--------------------------------------------------------;\n", + " +-----++--------------------------------------------------------;\n", + " -+----++--------------------------------------------------------;\n", + " --+---++--------------------------------------------------------;\n", + " ---+--++--------------------------------------------------------;\n", + " ----+-++--------------------------------------------------------;\n", + " +----+++--------------------------------------------------------;\n", + " -+---+++--------------------------------------------------------;\n", + " --+--+++--------------------------------------------------------;\n", + " ---+-+++--------------------------------------------------------;\n", + " ----++++--------------------------------------------------------;\n", + " +-------+-------------------------------------------------------;\n", + " -+------+-------------------------------------------------------]\n" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "let elec_num = Electrons.n_elec @@ Simulation.electrons simulation\n", + "\n", + "let alfa_num = Electrons.n_alfa @@ Simulation.electrons simulation\n", + "\n", + "let mo_class = MOClass.fci ~frozen_core mo_basis \n", + "\n", + "\n", + "let m l = \n", + " List.fold_left (fun accu i -> let j = i-1 in\n", + " Bitstring.logor accu (Bitstring.shift_left_one mo_num j)\n", + " ) (Bitstring.zero mo_num) l\n", + " \n", + "let occ_mask = m (MOClass.core_mos mo_class)\n", + "\n", + "let active_mask = m (MOClass.active_mos mo_class)\n", + "\n", + "\n", + "let neg_active_mask = Bitstring.lognot active_mask \n", + "\n", + " (* Here we generate the FCI space and filter out unwanted determinants\n", + " with excitations involving the core electrons. This should be improved. *)\n", + " \n", + "let permtutations m n = \n", + " let open Bitstring in\n", + " let rec aux k u rest =\n", + " if k=1 then\n", + " List.rev (u :: rest)\n", + " else\n", + " let t = (logor u (minus_one u)) in\n", + " let t' = plus_one t in\n", + " let t'' = shift_right (minus_one (logand (lognot t) (neg @@ lognot t))) (trailing_zeros u + 1) in \n", + " (*\n", + " let t'' = shift_right (minus_one (logand (lognot t) t')) (trailing_zeros u + 1) in\n", + " *)\n", + " (aux [@tailcall]) (k-1) (logor t' t'') (u :: rest)\n", + " in\n", + " aux (Util.binom n m) (minus_one (shift_left_one n m)) []\n", + "\n", + " \n", + "let spin_determinants =\n", + " Bitstring.permtutations alfa_num mo_num\n", + " (*\n", + " |> List.filter (fun b -> Bitstring.logand neg_active_mask b = occ_mask)\n", + " |> List.map (fun b -> Spindeterminant.of_bitstring b)\n", + " |> Array.of_list\n", + " *)\n" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "let ci = CI.make ~n_states:state det_space\n", + "\n", + "let ci_coef, ci_energy = Lazy.force ci.eigensystem " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Creates a function `f : Determinant.t -> bool` that returns `true` when the determinant has one or two electrons in the CABS." + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "let is_a_double det_space =\n", + " let mo_class = DeterminantSpace.mo_class det_space in\n", + " let mo_num = Array.length @@ MOClass.mo_class_array mo_class in\n", + " let m l =\n", + " List.fold_left (fun accu i ->\n", + " let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)\n", + " ) (Bitstring.zero mo_num) l\n", + " in\n", + " let aux_mask = m (MOClass.auxiliary_mos mo_class) in\n", + " fun k -> \n", + " let alfa =\n", + " Determinant.alfa k\n", + " |> Spindeterminant.bitstring\n", + " in\n", + " let beta =\n", + " Determinant.beta k\n", + " |> Spindeterminant.bitstring\n", + " in\n", + " let a = Bitstring.logand aux_mask alfa\n", + " and b = Bitstring.logand aux_mask beta\n", + " in\n", + " match Bitstring.popcount a + Bitstring.popcount b with\n", + " | 2 | 1 -> true\n", + " | 0 | _ -> false\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Permutation operator $p_{12}$ that generates a new determinant with electrons 1 and 2 swapped." + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "let p12 det_space =\n", + " let mo_class = DeterminantSpace.mo_class det_space in\n", + " let mo_num = Array.length @@ MOClass.mo_class_array mo_class in\n", + " let m l =\n", + " List.fold_left (fun accu i ->\n", + " let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)\n", + " ) (Bitstring.zero mo_num) l\n", + " in\n", + " let aux_mask = m (MOClass.auxiliary_mos mo_class) in\n", + " let not_aux_mask =\n", + " Bitstring.(shift_left_one mo_num (mo_num-1) |> minus_one |> logxor aux_mask)\n", + " in\n", + " fun k ->\n", + " let alfa =\n", + " Determinant.alfa k\n", + " |> Spindeterminant.bitstring\n", + " in\n", + " let beta =\n", + " Determinant.beta k\n", + " |> Spindeterminant.bitstring\n", + " in\n", + " let a = Bitstring.logand aux_mask alfa\n", + " and b = Bitstring.logand aux_mask beta\n", + " in\n", + " match Bitstring.popcount a, Bitstring.popcount b with\n", + " | 2, 0 \n", + " | 0, 2 -> Some (Determinant.negate_phase k)\n", + " | 1, 1 -> Some (Determinant.of_spindeterminants\n", + " (Spindeterminant.of_bitstring @@\n", + " Bitstring.(logor b (logand not_aux_mask alfa)) )\n", + " (Spindeterminant.of_bitstring @@\n", + " Bitstring.(logor a (logand not_aux_mask beta))\n", + " ) )\n", + " | 1, 0 \n", + " | 0, 1 -> Some k\n", + " | _ -> None\n", + " \n", + " \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Determinants of the FCI space" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "let in_dets = \n", + " DeterminantSpace.determinant_stream ci.CI.det_space \n", + " |> Util.stream_to_list\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Integral-based functions" ] }, { @@ -734,7 +1271,7 @@ ], "metadata": { "kernelspec": { - "display_name": "OCaml default", + "display_name": "OCaml 4.07.1", "language": "OCaml", "name": "ocaml-jupyter" }, diff --git a/Utils/Bitstring.ml b/Utils/Bitstring.ml index feec615..066e5ed 100644 --- a/Utils/Bitstring.ml +++ b/Utils/Bitstring.ml @@ -13,6 +13,7 @@ module One = struct let shift_left_one i = 1 lsl i let testbit x i = ( (x lsr i) land 1 ) = 1 let logor a b = a lor b + let neg a = - a let logxor a b = a lxor b let logand a b = a land b let lognot a = lnot a @@ -34,7 +35,7 @@ module One = struct let pp ppf s = - Format.fprintf ppf "@[@[%a@]i@]" (Util.pp_bitstring 64) + Format.fprintf ppf "@[@[%a@]@]" (Util.pp_bitstring 64) (Z.of_int s) end @@ -58,6 +59,7 @@ module Many = struct let logxor = Z.logxor let logand = Z.logand let lognot = Z.lognot + let neg = Z.neg let minus_one = Z.pred let plus_one = Z.succ let trailing_zeros = Z.trailing_zeros @@ -68,7 +70,7 @@ module Many = struct if z = Z.zero then 0 else Z.popcount z let pp ppf s = - Format.fprintf ppf "@[@[%a@]m@]" (Util.pp_bitstring (Z.numbits s)) s + Format.fprintf ppf "@[@[%a@]@]" (Util.pp_bitstring (Z.numbits s)) s end @@ -95,6 +97,10 @@ let is_zero = function | One x -> One.is_zero x | Many x -> Many.is_zero x +let neg = function +| One x -> One (One.neg x) +| Many x -> Many (Many.neg x) + let shift_left x i = match x with | One x -> One (One.shift_left x i) | Many x -> Many (Many.shift_left x i) @@ -184,9 +190,14 @@ let permtutations m n = if k=1 then List.rev (u :: rest) else - let t = (logor u (minus_one u)) in + let t = logor u @@ minus_one u in let t' = plus_one t in + let not_t = lognot t in + let neg_not_t = neg not_t in + let t'' = shift_right (minus_one @@ logand not_t neg_not_t) (trailing_zeros u + 1) in + (* let t'' = shift_right (minus_one (logand (lognot t) t')) (trailing_zeros u + 1) in + *) (aux [@tailcall]) (k-1) (logor t' t'') (u :: rest) in aux (Util.binom n m) (minus_one (shift_left_one n m)) [] diff --git a/Utils/Electrons.ml b/Utils/Electrons.ml index 642ae93..e99c2d5 100644 --- a/Utils/Electrons.ml +++ b/Utils/Electrons.ml @@ -27,4 +27,5 @@ let charge e = let n_alfa t = t.n_alfa let n_beta t = t.n_beta +let n_elec t = t.n_alfa + t.n_beta let multiplicity t = t.multiplicity diff --git a/Utils/Electrons.mli b/Utils/Electrons.mli index caf0122..d84244b 100644 --- a/Utils/Electrons.mli +++ b/Utils/Electrons.mli @@ -14,6 +14,9 @@ val make : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t val charge : t -> Charge.t (** Sum of the charges of the electrons. *) +val n_elec : t -> int +(** Number of electrons *) + val n_alfa : t -> int (** Number of alpha electrons *) diff --git a/Utils/math_functions.c b/Utils/math_functions.c index ba4a05a..353141f 100644 --- a/Utils/math_functions.c +++ b/Utils/math_functions.c @@ -46,7 +46,7 @@ CAMLprim int32_t popcnt(int64_t i) CAMLprim value popcnt_bytecode(value i) { - return copy_int32(__builtin_popcountll (i)); + return caml_copy_int32(__builtin_popcountll (Int64_val(i))); } @@ -58,7 +58,7 @@ CAMLprim int32_t trailz(int64_t i) CAMLprim value trailz_bytecode(value i) { - return copy_int32(__builtin_ctzll (i)); + return caml_copy_int32(__builtin_ctzll (Int64_val(i))); } CAMLprim int32_t leadz(int64_t i) @@ -69,7 +69,7 @@ CAMLprim int32_t leadz(int64_t i) CAMLprim value leadz_bytecode(value i) { - return copy_int32(__builtin_clzll (i)); + return caml_copy_int32(__builtin_clzll (Int64_val(i))); }