From 4dccd121e0da23254321db43a883f089e8b2686f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 30 Mar 2020 18:06:21 +0200 Subject: [PATCH 1/3] Envoi Nico --- Notebooks/F12_fast.ipynb | 900 ++++++++++++++++++++++++++++++++++++--- Simulation.ml | 16 +- 2 files changed, 842 insertions(+), 74 deletions(-) diff --git a/Notebooks/F12_fast.ipynb b/Notebooks/F12_fast.ipynb index ef18f66..81c558e 100644 --- a/Notebooks/F12_fast.ipynb +++ b/Notebooks/F12_fast.ipynb @@ -9,16 +9,126 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "heading_collapsed": true + }, "source": [ "## Initialization" ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 1, + "metadata": { + "ExecuteTime": { + "end_time": "2020-03-30T10:15:59.908500Z", + "start_time": "2020-03-30T10:15:58.731Z" + }, + "hidden": true, + "init_cell": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "val png_image : string -> unit = \n" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + }, + { + "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/4.10.0/lib/bytes: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/base64: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/base64/base64.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ocaml/compiler-libs: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ocaml/compiler-libs/ocamlcommon.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/result: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/result/result.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ppx_deriving/runtime: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ppx_deriving/runtime/ppx_deriving_runtime.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ppx_deriving_yojson/runtime: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ppx_deriving_yojson/runtime/ppx_deriving_yojson_runtime.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ocaml/unix.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/uuidm: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/uuidm/uuidm.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/easy-format: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/easy-format/easy_format.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/biniou: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/biniou/biniou.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/yojson: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/yojson/yojson.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/jupyter: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/jupyter/jupyter.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/jupyter/notebook: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/jupyter/notebook/jupyter_notebook.cma: loaded\n" + ] + }, + { + "data": { + "text/plain": [ + "val png_image : string -> Jupyter_notebook.display_id = \n" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/scemama/qp2/external/opam/4.10.0/lib/ocaml/compiler-libs/ocamlbytecomp.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ocaml/compiler-libs/ocamltoplevel.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ocaml/bigarray.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/lacaml: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/lacaml/lacaml.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/lacaml/top: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/lacaml/top/lacaml_top.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/astring: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/astring/astring.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/cmdliner: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/cmdliner/cmdliner.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/seq: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/stdlib-shims: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/stdlib-shims/stdlib_shims.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/fmt: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/fmt/fmt.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/fmt/fmt_cli.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/fmt/fmt_tty.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/re: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/re/re.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/alcotest: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/alcotest/alcotest.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/ocaml/str.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/zarith: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/zarith/zarith.cma: loaded\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/getopt: added to search path\n", + "/home/scemama/qp2/external/opam/4.10.0/lib/getopt/getopt.cma: loaded\n" + ] + } + ], "source": [ "let png_image = print_endline ;;\n", "\n", @@ -51,15 +161,24 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "### Modules" ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, + "execution_count": 2, + "metadata": { + "ExecuteTime": { + "end_time": "2020-03-30T10:16:00.033900Z", + "start_time": "2020-03-30T10:15:58.750Z" + }, + "hidden": true, + "init_cell": true + }, "outputs": [], "source": [ "#load \"Constants.cmo\";;\n", @@ -75,16 +194,36 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "### Printers" ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 3, + "metadata": { + "ExecuteTime": { + "end_time": "2020-03-30T10:16:00.095400Z", + "start_time": "2020-03-30T10:15:58.760Z" + }, + "hidden": true, + "init_cell": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "val pp_mo : Format.formatter -> MOBasis.t -> unit = \n" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "#install_printer AngularMomentum.pp_string ;;\n", "#install_printer Basis.pp ;;\n", @@ -113,25 +252,251 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "heading_collapsed": true + }, "source": [ "## Run" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "### Simulation\n" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": { + "ExecuteTime": { + "end_time": "2020-03-30T10:16:00.287500Z", + "start_time": "2020-03-30T10:15:58.769Z" + }, + "hidden": true, + "init_cell": true, "scrolled": false }, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "val basis_filename : string = \"/home/scemama/qp2/data/basis/6-31g\"\n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val aux_basis_filename : string = \"/home/scemama/qp2/data/basis/cc-pvdz\"\n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val nuclei : Nuclei.t = [|(Element.C, 0.0000 0.0000 0.0000)|]\n" + ] + }, + "execution_count": 4, + "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": [ + "val multiplicity : int = 1\n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val state : int = 1\n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val basis : Basis.t =\n", + " \n", + " Atomic Basis set\n", + " ----------------\n", + "\n", + "-----------------------------------------------------------------------\n", + " # Angular Coordinates (Bohr) Exponents Coefficients\n", + " Momentum X Y Z\n", + "-----------------------------------------------------------------------\n", + " 1-3 S 0.0000 0.0000 0.0000 3.04752488e+03 1.83473713e-03\n", + " 4.57369518e+02 1.40373228e-02\n", + " 1.03948685e+02 6.88426223e-02\n", + " 2.92101553e+01 2.32184443e-01\n", + " 9.28666296e+00 4.67941348e-01\n", + " 3.16392696e+00 3.62311985e-01\n", + " \n", + " 7.86827235e+00 -1.19332420e-01\n", + " 1.88128854e+00 -1.60854152e-01\n", + " 5.44249258e-01 1.14345644e+00\n", + " \n", + " 1.68714478e-01 1.00000000e+00\n", + " \n", + " \n", + "-----------------------------------------------------------------------\n", + " 4-9 P 0.0000 0.0000 0.0000 7.86827235e+00 6.89990666e-02\n", + " 1.88128854e+00 3.16423961e-01\n", + " 5.44249258e-01 7.44308291e-01\n", + " \n", + " 1.68714478e-01 1.00000000e+00\n", + " \n", + " \n", + "-----------------------------------------------------------------------\n", + "\n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val aux_basis : Basis.t =\n", + " \n", + " Atomic Basis set\n", + " ----------------\n", + "\n", + "-----------------------------------------------------------------------\n", + " # Angular Coordinates (Bohr) Exponents Coefficients\n", + " Momentum X Y Z\n", + "-----------------------------------------------------------------------\n", + " 1-6 S 0.0000 0.0000 0.0000 3.04752488e+03 1.83473713e-03\n", + " 4.57369518e+02 1.40373228e-02\n", + " 1.03948685e+02 6.88426223e-02\n", + " 2.92101553e+01 2.32184443e-01\n", + " 9.28666296e+00 4.67941348e-01\n", + " 3.16392696e+00 3.62311985e-01\n", + " \n", + " 7.86827235e+00 -1.19332420e-01\n", + " 1.88128854e+00 -1.60854152e-01\n", + " 5.44249258e-01 1.14345644e+00\n", + " \n", + " 1.68714478e-01 1.00000000e+00\n", + " \n", + " 6.66500000e+03 6.92000000e-04\n", + " 1.00000000e+03 5.32900000e-03\n", + " 2.28000000e+02 2.70770000e-02\n", + " 6.47100000e+01 1.01718000e-01\n", + " 2.10600000e+01 2.74740000e-01\n", + " 7.49500000e+00 4.48564000e-01\n", + " 2.79700000e+00 2.85074000e-01\n", + " 5.21500000e-01 1.52040000e-02\n", + " 1.59600000e-01 -3.19100000e-03\n", + " \n", + " 6.66500000e+03 -1.46000000e-04\n", + " 1.00000000e+03 -1.15400000e-03\n", + " 2.28000000e+02 -5.72500000e-03\n", + " 6.47100000e+01 -2.33120000e-02\n", + " 2.10600000e+01 -6.39550000e-02\n", + " 7.49500000e+00 -1.49981000e-01\n", + " 2.79700000e+00 -1.27262000e-01\n", + " 5.21500000e-01 5.44529000e-01\n", + " 1.59600000e-01 5.80496000e-01\n", + " \n", + " 1.59600000e-01 1.00000000e+00\n", + " \n", + " \n", + "-----------------------------------------------------------------------\n", + " 4-15 P 0.0000 0.0000 0.0000 7.86827235e+00 6.89990666e-02\n", + " 1.88128854e+00 3.16423961e-01\n", + " 5.44249258e-01 7.44308291e-01\n", + " \n", + " 1.68714478e-01 1.00000000e+00\n", + " \n", + " 9.43900000e+00 3.81090000e-02\n", + " 2.00200000e+00 2.09480000e-01\n", + " 5.45600000e-01 5.08557000e-01\n", + " 1.51700000e-01 4.68842000e-01\n", + " \n", + " 1.51700000e-01 1.00000000e+00\n", + " \n", + " \n", + "-----------------------------------------------------------------------\n", + " 19-24 D 0.0000 0.0000 0.0000 5.50000000e-01 1.00000000e+00\n", + " \n", + " \n", + "-----------------------------------------------------------------------\n", + "\n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "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_g =\n", + " [ 0.220900 1.004000 3.622000 12.160000 45.870000 254.400000 ];\n", + " expo_g_inv =\n", + " [ 4.526935 0.996016 0.276091 0.082237 0.021801 0.003931 ]}}\n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val charge : int = 0\n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val simulation : Simulation.t = \n" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "let basis_filename = \"/home/scemama/qp2/data/basis/6-31g\" \n", "let aux_basis_filename = \"/home/scemama/qp2/data/basis/cc-pvdz\" \n", @@ -155,9 +520,29 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 5, + "metadata": { + "ExecuteTime": { + "end_time": "2020-03-30T10:16:00.298200Z", + "start_time": "2020-03-30T10:15:58.778Z" + }, + "hidden": true, + "init_cell": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "val n_elec_alfa : int = 3\n", + "val n_elec_beta : int = 3\n", + "val n_elec : int = 6\n" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "let n_elec_alfa, n_elec_beta, n_elec = \n", " let e = Simulation.electrons simulation in\n", @@ -166,16 +551,131 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "### Hartree-Fock" ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 6, + "metadata": { + "ExecuteTime": { + "end_time": "2020-03-30T10:16:00.474000Z", + "start_time": "2020-03-30T10:15:58.786Z" + }, + "hidden": true, + "init_cell": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "15 significant shell pairs computed in 0.013028 seconds\n", + "1\n", + "2\n", + "3\n", + "6\n", + "Computed ERIs in 0.083199 seconds\n", + "MOs =\n", + "\n", + "\n", + " -- 1 -- -- 2 -- -- 3 -- -- 4 -- -- 5 --\n", + " 1 0.995189 -0.236456 0 -0 0\n", + " 2 0.0272807 0.573253 0 2.23673E-15 0\n", + " 3 -0.00753288 0.50151 0 -2.90078E-15 0\n", + " ... ... ... ... ...\n", + " 7 -0 -0 0.028973 0.7796 -0.0935272\n", + " 8 -0 -0 0.269658 0.0144514 0.61693\n", + " 9 -0 -0 -0.260237 0.10177 0.628851\n", + " \n", + "\n", + " -- 6 -- -- 7 -- -- 8 -- -- 9 --\n", + " 1 -0 0.0649867 -0 0\n", + " 2 -4.19726E-16 -1.62901 5.23178E-16 0\n", + " 3 -0 1.63955 -1.33409E-15 0\n", + " ... ... ... ...\n", + " 7 0.0892087 0 -0.920636 0.088067\n", + " 8 0.830283 0 -0.0170658 -0.580913\n", + " 9 -0.801275 0 -0.120181 -0.592138\n", + " \n" + ] + }, + { + "data": { + "text/plain": [ + "val hf : HartreeFock.t = \n", + "======================================================================\n", + " Restricted Hartree-Fock \n", + "======================================================================\n", + "\n", + " ------------------------------------------------------------\n", + " # HF energy Convergence HOMO-LUMO\n", + " ------------------------------------------------------------\n", + " 1 -36.54848150 7.4861e-01 0.2723\n", + " ------------------------------------------------------------\n", + "\n", + "\n", + " ============================================================\n", + " One-electron energy -52.6699343629\n", + " Kinetic 42.4526846308\n", + " Potential -95.1226189937\n", + " -------------------------------------------------------- \n", + " Two-electron energy 16.1214528647\n", + " Coulomb 21.8535626014\n", + " Exchange -5.7321097367\n", + " -------------------------------------------------------- \n", + " HF HOMO 0.8237527941\n", + " HF LUMO 8.2337246898\n", + " HF LUMO-LUMO 7.4099718957\n", + " -------------------------------------------------------- \n", + " Electronic energy -36.5484814981\n", + " Nuclear repulsion 0.0000000000\n", + " Hartree-Fock energy -36.5484814981\n", + " ============================================================\n", + " \n", + "\n" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val mo_basis : MOBasis.t =\n", + " Eigenvalues: -11.127416 -0.632049 -0.269831 0.131526 0.151773 \n", + " -- 1 -- -- 2 -- -- 3 -- -- 4 -- -- 5 --\n", + " 1 0.995189 -0.236456 0 -0 0\n", + " 2 0.0272807 0.573253 0 2.23673E-15 0\n", + " 3 -0.00753288 0.50151 0 -2.90078E-15 0\n", + " ... ... ... ... ...\n", + " 7 -0 -0 0.028973 0.7796 -0.0935272\n", + " 8 -0 -0 0.269658 0.0144514 0.61693\n", + " 9 -0 -0 -0.260237 0.10177 0.628851\n", + " \n", + " Eigenvalues: 0.786913 0.828232 0.869410 0.849164 \n", + " -- 6 -- -- 7 -- -- 8 -- -- 9 --\n", + " 1 -0 0.0649867 -0 0\n", + " 2 -4.19726E-16 -1.62901 5.23178E-16 0\n", + " 3 -0 1.63955 -1.33409E-15 0\n", + " ... ... ... ...\n", + " 7 0.0892087 0 -0.920636 0.088067\n", + " 8 0.830283 0 -0.0170658 -0.580913\n", + " 9 -0.801275 0 -0.120181 -0.592138\n", + " \n", + " \n" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "let hf = HartreeFock.make ~guess:`Hcore ~max_scf:1 simulation ;;\n", "\n", @@ -184,23 +684,83 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "heading_collapsed": true + }, "source": [ "# FCI-F12" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "## Common functions" ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 7, + "metadata": { + "ExecuteTime": { + "end_time": "2020-03-30T10:16:00.507100Z", + "start_time": "2020-03-30T10:15:58.791Z" + }, + "hidden": true, + "init_cell": true + }, + "outputs": [ + { + "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_g =\n", + " [ 0.220900 1.004000 3.622000 12.160000 45.870000 254.400000 ];\n", + " expo_g_inv =\n", + " [ 4.526935 0.996016 0.276091 0.082237 0.021801 0.003931 ]}}\n" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val mo_num : int = 9\n" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val pp_spindet : Format.formatter -> Spindeterminant.t -> unit = \n" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val pp_det : Format.formatter -> Determinant.t -> unit = \n" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "let f12 = Util.of_some @@ Simulation.f12 simulation \n", "\n", @@ -219,9 +779,116 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 8, + "metadata": { + "ExecuteTime": { + "end_time": "2020-03-30T10:16:03.181400Z", + "start_time": "2020-03-30T10:15:58.799Z" + }, + "hidden": true, + "init_cell": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "val simulation_aux : Simulation.t = \n" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "66 significant shell pairs computed in 0.122791 seconds\n", + "1\n", + "2\n", + "3\n", + "6\n", + "9\n", + "10\n", + "11\n", + "12\n", + "15\n", + "18\n", + "Computed ERIs in 2.253815 seconds\n" + ] + }, + { + "data": { + "text/plain": [ + "val aux_basis : MOBasis.t =\n", + " Eigenvalues: -11.127416 -0.632049 -0.269831 0.131526 0.151773 \n", + " -- 1 -- -- 2 -- -- 3 -- -- 4 -- -- 5 --\n", + " 1 0.995189 -0.236456 8.68506E-12 4.39594E-12 -7.46805E-12\n", + " 2 0.0272807 0.573253 -4.95565E-12 -5.67405E-12 -4.47366E-13\n", + " 3 -0.00753288 0.50151 -8.79918E-13 -7.15566E-13 2.36287E-13\n", + " ... ... ... ... ...\n", + " 22 0 0 -0 -0 0\n", + " 23 0 0 -0 -0 0\n", + " 24 0 0 -0 -0 0\n", + " \n", + " Eigenvalues: 0.786913 0.828232 0.869410 0.849164 1.471540 \n", + " -- 6 -- -- 7 -- -- 8 -- -- 9 -- -- 10 --\n", + " 1 -9.38797E-12 0.0649867 -5.92171E-12 -1.01808E-11 2.71822E-13\n", + " 2 1.16966E-12 -1.62901 -8.06973E-13 6.55834E-12 -5.57185E-14\n", + " 3 6.63605E-13 1.63955 5.0525E-13 1.12933E-12 -2.61365E-15\n", + " ... ... ... ... ...\n", + " 22 0 0 0 0 -0.866025\n", + " 23 0 0 0 0 0\n", + " 24 0 0 0 0 0\n", + " \n", + " Eigenvalues: 1.031908 1.024910 1.028381 1.247024 25.947498 \n", + " -- 11 -- -- 12 -- -- 13 -- -- 14 -- -- 15 --\n", + " 1 4.00353E-10 -8.57233E-11 1.06969E-10 3.79632 -2.36312E-10\n", + " 2 1.4024E-10 1.19157E-12 9.07568E-12 -10.5897 -8.33662E-10\n", + " 3 3.79392E-12 7.24542E-12 -1.17682E-11 53.7881 -1.17393E-10\n", + " ... ... ... ... ...\n", + " 22 0 -0 0 -0 0\n", + " 23 0 -0 0 -0 1.05805E-13\n", + " 24 0 -0 0 -0 0\n", + " \n", + " Eigenvalues: 25.902263 25.896210 18.767112 377.619959 1.471980 \n", + " -- 16 -- -- 17 -- -- 18 -- -- 19 -- -- 20 --\n", + " 1 -1.81712E-09 7.6446E-09 27.2935 436.547 7.08629E-12\n", + " 2 -7.15883E-10 2.38277E-09 -127.494 -10.9671 -2.66825E-13\n", + " 3 -1.90088E-11 -1.74694E-11 -21.2025 -23.8199 -3.9535E-13\n", + " ... ... ... ... ...\n", + " 22 0 0 -0 0 -0.0208586\n", + " 23 0 0 -0 0 -0.442591\n", + " 24 0 0 -0 0 0.0417173\n", + " \n", + " Eigenvalues: 1.478279 1.478333 1.471737 \n", + " -- 21 -- -- 22 -- -- 23 --\n", + " 1 -1.36007E-11 -2.98191E-11 6.74696E-10\n", + " 2 5.12555E-13 1.12961E-12 -2.54412E-11\n", + " 3 7.58688E-13 1.66236E-12 -3.76295E-11\n", + " ... ... ...\n", + " 22 -0.0667224 0.492879 0.0467273\n", + " 23 0.890945 0.100348 0.0161496\n", + " 24 0.133445 -0.985758 -0.0934546\n", + " \n", + " \n" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": [ + "val aux_num : int = 23\n" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "let simulation_aux = \n", " let charge = Charge.to_int @@ Simulation.charge simulation \n", @@ -247,9 +914,45 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 9, + "metadata": { + "ExecuteTime": { + "end_time": "2020-03-30T10:16:41.349200Z", + "start_time": "2020-03-30T10:15:58.806Z" + }, + "hidden": true, + "init_cell": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "66 significant shell pairs computed in 1.895893 seconds\n", + "1\n", + "2\n", + "3\n", + "6\n", + "9\n", + "10\n", + "11\n", + "12\n", + "15\n", + "18\n", + "Computed ERIs in 35.957611 seconds\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "4-idx transformation \n", + "23 / 23\n", + "4-idx transformation \n", + "23 / 23\n" + ] + } + ], "source": [ "let () = ignore @@ MOBasis.f12_ints aux_basis\n", "let () = ignore @@ MOBasis.two_e_ints aux_basis" @@ -257,14 +960,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "## Integral-based functions" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "\\begin{equation}\n", "\\langle I | \\hat{H} | J \\rangle = \\begin{cases}\n", @@ -287,7 +994,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let cancel_singles = false \n", @@ -313,7 +1022,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "### H integrals" ] @@ -321,7 +1032,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let h_one =\n", @@ -346,7 +1059,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "### F12 integrals" ] @@ -354,7 +1069,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let f_two = \n", @@ -386,14 +1103,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "## Determinant-based functions" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "### Integrals" ] @@ -401,7 +1122,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let f12_integrals mo_basis =\n", @@ -437,7 +1160,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "### Determinant space" ] @@ -445,7 +1170,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let is_a_double det_space =\n", @@ -478,7 +1205,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let in_space = \n", @@ -514,6 +1243,7 @@ "cell_type": "code", "execution_count": null, "metadata": { + "hidden": true, "scrolled": false }, "outputs": [], @@ -527,14 +1257,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "Permutation operator $p_{12}$ that generates a new determinant with electrons 1 and 2 swapped." ] }, { "cell_type": "raw", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "let p12 det_space =\n", " let mo_class = DeterminantSpace.mo_class det_space in\n", @@ -578,7 +1312,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "## Matrices $\\langle I | H | \\alpha \\rangle$ and $\\langle I | F | \\alpha \\rangle$ " ] @@ -587,6 +1323,7 @@ "cell_type": "code", "execution_count": null, "metadata": { + "hidden": true, "scrolled": true }, "outputs": [], @@ -617,6 +1354,7 @@ { "cell_type": "raw", "metadata": { + "hidden": true, "raw_mimetype": "text/markdown" }, "source": [ @@ -654,7 +1392,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "Function to generate all intermediate external determinants $|\\alpha \\rangle$ between $|I\\rangle$ and $|J\\rangle$, with a positive phase factor.\n", "\n", @@ -664,7 +1404,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let generate_alphas ki kj exc =\n", @@ -792,7 +1534,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let check n integral_value exc =\n", @@ -870,7 +1614,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let array_3_init d1 d2 d3 f =\n", @@ -925,14 +1671,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "******" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "heading_collapsed": true + }, "source": [ "# 1. 0 " ] @@ -940,7 +1690,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "let m_0111_1H_1F = \n", @@ -1106,7 +1858,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "heading_collapsed": true + }, "source": [ "# 2. 1" ] @@ -1115,6 +1869,7 @@ "cell_type": "code", "execution_count": null, "metadata": { + "hidden": true, "scrolled": false }, "outputs": [], @@ -1395,14 +2150,18 @@ }, { "cell_type": "raw", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "Timing : 40.020974 0.019294" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "heading_collapsed": true + }, "source": [ "# 3. 2" ] @@ -1410,7 +2169,9 @@ { "cell_type": "code", "execution_count": 129, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [ { "data": { @@ -1918,7 +2679,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "heading_collapsed": true + }, "source": [ "# 4. 3" ] @@ -1926,7 +2689,9 @@ { "cell_type": "code", "execution_count": 130, - "metadata": {}, + "metadata": { + "hidden": true + }, "outputs": [ { "data": { @@ -2030,7 +2795,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "---------------------------" ] @@ -2039,6 +2806,7 @@ "cell_type": "code", "execution_count": null, "metadata": { + "hidden": true, "scrolled": false }, "outputs": [], @@ -2179,9 +2947,9 @@ } ], "metadata": { - "celltoolbar": "Raw Cell Format", + "celltoolbar": "Initialization Cell", "kernelspec": { - "display_name": "OCaml 4.07.1", + "display_name": "OCaml 4.10.0", "language": "OCaml", "name": "ocaml-jupyter" }, @@ -2192,7 +2960,7 @@ "name": "OCaml", "nbconverter_exporter": null, "pygments_lexer": "OCaml", - "version": "4.07.1" + "version": "4.10.0" } }, "nbformat": 4, diff --git a/Simulation.ml b/Simulation.ml index 6a010b4..c8287bd 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -1,6 +1,6 @@ type t = { charge : Charge.t; - electrons : Electrons.t ; + electrons : Electrons.t; nuclei : Nuclei.t; basis : Basis.t; ao_basis : AOBasis.t; @@ -19,12 +19,12 @@ let mu_erf t = t.mu_erf let f12 t = t.f12 let make ?cartesian:(cartesian=false) - ?multiplicity:(multiplicity=1) - ?charge:(charge=0) - ?f12 - ?mu_erf - ~nuclei - basis + ?multiplicity:(multiplicity=1) + ?charge:(charge=0) + ?f12 + ?mu_erf + ~nuclei + basis = (* Tune Garbage Collector *) @@ -35,7 +35,7 @@ let make ?cartesian:(cartesian=false) Electrons.make ~multiplicity ~charge nuclei in - let charge = + let chargre = Charge.(Nuclei.charge nuclei + Electrons.charge electrons) in From a1c29fd3f1c511803a13a3e4efe1b170f7fc0b77 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 14 Apr 2020 15:53:07 +0200 Subject: [PATCH 2/3] Fix? --- INSTALL.md | 2 +- Simulation.ml | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/INSTALL.md b/INSTALL.md index fc12576..9140d3f 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -1,7 +1,7 @@ # Generic installation ```bash -opam install ocamlbuild lacaml mpi getopt alcotest zarith +opam install ocamlbuild ocamlfind lacaml mpi getopt alcotest zarith ``` diff --git a/Simulation.ml b/Simulation.ml index c8287bd..2a1726f 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -35,15 +35,15 @@ let make ?cartesian:(cartesian=false) Electrons.make ~multiplicity ~charge nuclei in - let chargre = + let charge = Charge.(Nuclei.charge nuclei + Electrons.charge electrons) in - let ao_basis = + let ao_basis = AOBasis.make ?f12 ~basis ~cartesian nuclei in - let nuclear_repulsion = + let nuclear_repulsion = Nuclei.repulsion nuclei in From 11151c997e3d437163723d217601e2c7b5f676d1 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 16 Apr 2020 19:49:23 +0200 Subject: [PATCH 3/3] Introduced multipole integrals --- Basis/AOBasis.ml | 8 +- Basis/AOBasis.mli | 3 + Basis/Basis.ml | 6 ++ Basis/Basis.mli | 5 ++ Basis/GeneralBasis.ml | 40 ++++++--- Basis/GeneralBasis.mli | 16 ++-- Basis/Multipole.ml | 199 +++++++++++++++++++++++++++++++++++++++++ Basis/Multipole.mli | 34 +++++++ Basis/Overlap.ml | 124 ++++++++++++------------- Nuclei/Nuclei.ml | 30 +++++-- Nuclei/Nuclei.mli | 3 + 11 files changed, 383 insertions(+), 85 deletions(-) create mode 100644 Basis/Multipole.ml create mode 100644 Basis/Multipole.mli diff --git a/Basis/AOBasis.ml b/Basis/AOBasis.ml index fa2cd0b..4499393 100644 --- a/Basis/AOBasis.ml +++ b/Basis/AOBasis.ml @@ -5,6 +5,7 @@ type t = { basis : Basis.t ; overlap : Overlap.t lazy_t; + multipole : Multipole.t lazy_t; ortho : Orthonormalization.t lazy_t; eN_ints : NucInt.t lazy_t; kin_ints : KinInt.t lazy_t; @@ -17,6 +18,7 @@ type t = let basis t = t.basis let overlap t = Lazy.force t.overlap +let multipole t = Lazy.force t.multipole let ortho t = Lazy.force t.ortho let eN_ints t = Lazy.force t.eN_ints let kin_ints t = Lazy.force t.kin_ints @@ -62,7 +64,11 @@ let make ~cartesian ~basis ?f12 nuclei = ScreenedERI.of_basis basis ) in - { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; + let multipole = lazy ( + Multipole.of_basis basis + ) in + + { basis ; overlap ; multipole ; ortho ; eN_ints ; kin_ints ; ee_ints ; ee_lr_ints ; f12_ints ; f12_over_r12_ints ; cartesian ; } diff --git a/Basis/AOBasis.mli b/Basis/AOBasis.mli index 328bf4d..d00131b 100644 --- a/Basis/AOBasis.mli +++ b/Basis/AOBasis.mli @@ -10,6 +10,9 @@ val basis : t -> Basis.t val overlap : t -> Overlap.t (** Overlap matrix *) +val multipole : t -> Multipole.t +(** Multipole matrices *) + val ortho : t -> Orthonormalization.t (** Orthonormalization matrix of the overlap *) diff --git a/Basis/Basis.ml b/Basis/Basis.ml index e831d8f..9a3302f 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -93,6 +93,12 @@ let of_nuclei_and_basis_filename ~nuclei filename = in of_nuclei_and_general_basis nuclei general_basis +let of_nuclei_and_basis_string ~nuclei str = + let general_basis = + GeneralBasis.of_string str + in + of_nuclei_and_general_basis nuclei general_basis + let of_nuclei_and_basis_filenames ~nuclei filenames = let general_basis = diff --git a/Basis/Basis.mli b/Basis/Basis.mli index 83f010d..ea9297e 100644 --- a/Basis/Basis.mli +++ b/Basis/Basis.mli @@ -18,6 +18,11 @@ val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> string -> t from a file. *) +val of_nuclei_and_basis_string : nuclei:Nuclei.t -> string -> t +(** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t} + from a string. + *) + val of_nuclei_and_basis_filenames : nuclei:Nuclei.t -> string list -> t (** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t} from multiple files. diff --git a/Basis/GeneralBasis.ml b/Basis/GeneralBasis.ml index 4b8dcde..3f05319 100644 --- a/Basis/GeneralBasis.ml +++ b/Basis/GeneralBasis.ml @@ -17,23 +17,24 @@ exception No_shell exception Malformed_shell of string -let read_shell ic = +let read_shell line_stream = try let shell, n = try - let line = input_line ic in + let line = Stream.next line_stream in if String.trim line = "$END" then raise End_of_file; Scanf.sscanf line " %c %d " (fun shell n -> shell, n) with | End_of_file -> raise No_shell + | Stream.Failure -> raise No_shell | Scanf.Scan_failure m -> raise (Malformed_shell m) in let rec loop = function | 0 -> [] | i -> let contraction = - let line = (input_line ic) in + let line = Stream.next line_stream in try Scanf.sscanf line " %_d %f %f " (fun exponent coefficient -> { exponent ; coefficient }) with _ -> raise (Malformed_shell (Printf.sprintf @@ -50,36 +51,46 @@ let read_shell ic = -let read_element ic = +let read_element line_stream = try + let line = Stream.next line_stream in let element = - Scanf.sscanf (input_line ic) " %s " Element.of_string + Scanf.sscanf line " %s " Element.of_string in let rec loop () = - match read_shell ic with + match read_shell line_stream with | Some shell -> shell :: loop () | None -> [] in Some (element, Array.of_list (loop ()) ) with - | End_of_file -> None + | Stream.Failure -> None - -let read filename = - let ic = open_in filename in +let read_stream line_stream = let rec loop accu = try - match read_element ic with + match read_element line_stream with | Some e -> loop (e :: accu) | None -> accu with Element.ElementError _ -> loop accu in - loop [] + loop [] + +let read filename = + let ic = open_in filename in + let line_stream = + Stream.from (fun _ -> + try Some (input_line ic) + with End_of_file -> None ) + in + let result = read_stream line_stream in + close_in ic; + result let combine basis_list = @@ -124,3 +135,8 @@ let to_string (name, contracted_shell_array) = Printf.sprintf "%s\n%s" name (string_of_contracted_shell_array contracted_shell_array) +let of_string input_string = + String.split_on_char '\n' input_string + |> Stream.of_list + |> read_stream + diff --git a/Basis/GeneralBasis.mli b/Basis/GeneralBasis.mli index 138022a..2e0d540 100644 --- a/Basis/GeneralBasis.mli +++ b/Basis/GeneralBasis.mli @@ -46,7 +46,7 @@ exception Malformed_shell of string val read : string -> t -(** Reads a basis set file and return an association list where +(** Reads a basis-set file and return an association list where the key is an {!Element.t} and the value is the parsed basis set. *) @@ -65,8 +65,10 @@ val read_many : string list -> t -val read_element : in_channel -> element_basis option -(** Reads an element from the input channel [ic]. For example, +val read_element : string Stream.t -> element_basis option +(** Reads an element from the input [string Stream]. The [string Stream] is a +stream of lines, like a text file split on the end-of-line character. +For example, {[ HYDROGEN S 3 @@ -95,8 +97,10 @@ Some *) -val read_shell : in_channel -> general_contracted_shell option -(** Reads a shell from the input channel [ic]. For example, +val read_shell : string Stream.t -> general_contracted_shell option +(** Reads a shell from the input [string Stream]. The [string Stream] is a +stream of lines, like a text file split on the end-of-line character. +For example, {[ S 3 1 13.0100000 0.0196850 @@ -119,4 +123,6 @@ Some val to_string : string * (general_contracted_shell array) -> string (** Pretty-prints the basis set of an {Element.t}. *) +val of_string : string -> t +(** Reads a GAMESS-formatted string. *) diff --git a/Basis/Multipole.ml b/Basis/Multipole.ml new file mode 100644 index 0000000..b0289bf --- /dev/null +++ b/Basis/Multipole.ml @@ -0,0 +1,199 @@ +open Util +open Constants +open Lacaml.D + +type t = Mat.t array +(* +[| "x"; "y"; "z"; "x2"; "y2"; "z2" |] +*) + +module Am = AngularMomentum +module Bs = Basis +module Co = Coordinate +module Cs = ContractedShell +module Csp = ContractedShellPair +module Po = Powers +module Psp = PrimitiveShellPair + +let multiply a b = + let n = Mat.dim1 a in + let c = Mat.create n n in + Mat.cpab c a b; + c + +let x0 t = t.(0) +let y0 t = t.(1) +let z0 t = t.(2) +let x1 t = t.(3) +let y1 t = t.(4) +let z1 t = t.(5) +let x2 t = t.(6) +let y2 t = t.(7) +let z2 t = t.(8) + +let matrix_x t = multiply (x1 t) @@ multiply (y0 t) (z0 t) +let matrix_y t = multiply (x0 t) @@ multiply (y1 t) (z0 t) +let matrix_z t = multiply (x0 t) @@ multiply (y0 t) (z1 t) +let matrix_x2 t = multiply (x2 t) @@ multiply (y0 t) (z0 t) +let matrix_y2 t = multiply (x0 t) @@ multiply (y2 t) (z0 t) +let matrix_z2 t = multiply (x0 t) @@ multiply (y0 t) (z2 t) +let matrix_xy t = multiply (x1 t) @@ multiply (y1 t) (z0 t) +let matrix_yz t = multiply (x0 t) @@ multiply (y1 t) (z1 t) +let matrix_zx t = multiply (x1 t) @@ multiply (y0 t) (z1 t) + +let cutoff = integrals_cutoff + +let to_powers x = + let open Zkey in + match to_powers x with + | Six x -> x + | _ -> assert false + +(** Computes all the integrals of the contracted shell pair *) +let contracted_class shell_a shell_b : float Zmap.t array = + + match Csp.make shell_a shell_b with + | None -> Array.init 9 (fun _ -> Zmap.create 0) + | Some shell_p -> + begin + + (* Pre-computation of integral class indices *) + let class_indices = Csp.zkey_array shell_p in + + let contracted_class = + Array.init 9 (fun _ -> Array.make (Array.length class_indices) 0.) + in + + let a_minus_b = + Csp.a_minus_b shell_p + in + let norm_coef_scales = + Csp.norm_scales shell_p + in + + (* Compute all integrals in the shell for each pair of significant shell pairs *) + + let xyz_of_int k = + match k with + | 0 -> Co.X + | 1 -> Co.Y + | _ -> Co.Z + in + + List.iter (fun (coef_prod, psp) -> + (** Screening on the product of coefficients *) + if (abs_float coef_prod) > 1.e-6*.cutoff then + begin + let expo_inv = Psp.exponent_inv psp + and center_pa = Psp.center_minus_a psp + and xa = Co.get Co.X @@ Cs.center shell_a + and ya = Co.get Co.Y @@ Cs.center shell_a + and za = Co.get Co.Z @@ Cs.center shell_a + in + + Array.iteri (fun i key -> + let (angMomA, angMomB) = to_powers key in + (* 1D Overlap *) + let f k = + let xyz = xyz_of_int k in + Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB) + expo_inv + (Co.get xyz a_minus_b, Co.get xyz center_pa) + in + (* 1D *) + let g k = + let xyz = xyz_of_int k in + Overlap_primitives.hvrr (Po.get xyz angMomA + 1, Po.get xyz angMomB) + expo_inv + (Co.get xyz a_minus_b, Co.get xyz center_pa) + in + (* 1D *) + let h k = + let xyz = xyz_of_int k in + Overlap_primitives.hvrr (Po.get xyz angMomA + 2, Po.get xyz angMomB) + expo_inv + (Co.get xyz a_minus_b, Co.get xyz center_pa) + in + let norm = norm_coef_scales.(i) in + let f0, f1, f2, g0, g1, g2, h0, h1, h2 = + f 0, f 1, f 2, g 0, g 1, g 2, h 0, h 1, h 2 + in + let x = g0 +. f0 *. xa in + let y = g1 +. f1 *. ya in + let z = g2 +. f2 *. za in + let x2 = h0 +. xa *. (2. *. x -. xa *. f0) in + let y2 = h1 +. ya *. (2. *. y -. ya *. f1) in + let z2 = h2 +. za *. (2. *. z -. za *. f2) in + let c = contracted_class in + let d = coef_prod *. norm in + c.(0).(i) <- c.(0).(i) +. d *. f0 ; + c.(1).(i) <- c.(1).(i) +. d *. f1 ; + c.(2).(i) <- c.(2).(i) +. d *. f2 ; + c.(3).(i) <- c.(3).(i) +. d *. x ; + c.(4).(i) <- c.(4).(i) +. d *. y ; + c.(5).(i) <- c.(5).(i) +. d *. z ; + c.(6).(i) <- c.(6).(i) +. d *. x2 ; + c.(7).(i) <- c.(7).(i) +. d *. y2 ; + c.(8).(i) <- c.(8).(i) +. d *. z2 ; + ) class_indices + end + ) (Csp.coefs_and_shell_pairs shell_p); + let result = + Array.map (fun c -> Zmap.create (Array.length c) ) contracted_class + in + for j=0 to Array.length result do + let rj = result.(j) in + let cj = contracted_class.(j) in + Array.iteri (fun i key -> Zmap.add rj key cj.(i)) class_indices + done; + result + end + + +(** Create multipole matrices *) +let of_basis basis = + let to_powers x = + let open Zkey in + match to_powers x with + | Three x -> x + | _ -> assert false + in + + let n = Bs.size basis + and shell = Bs.contracted_shells basis + in + + let result = Array.init 9 (fun _ -> Mat.create n n) in + for j=0 to (Array.length shell) - 1 do + for i=0 to j do + (* Compute all the integrals of the class *) + let cls = + contracted_class shell.(i) shell.(j) + in + + for k=0 to 8 do + Array.iteri (fun j_c powers_j -> + let j_c = Cs.index shell.(j) + j_c + 1 in + let xj = to_powers powers_j in + Array.iteri (fun i_c powers_i -> + let i_c = Cs.index shell.(i) + i_c + 1 in + let xi = to_powers powers_i in + let key = + Zkey.of_powers_six xi xj + in + let value = + try Zmap.find cls.(k) key + with Not_found -> 0. + in + result.(k).{i_c,j_c} <- value; + result.(k).{j_c,i_c} <- value; + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) + done; + done; + done; + for k=0 to 8 do + Mat.detri result.(k); + done; + result + diff --git a/Basis/Multipole.mli b/Basis/Multipole.mli new file mode 100644 index 0000000..64d09de --- /dev/null +++ b/Basis/Multipole.mli @@ -0,0 +1,34 @@ +(** Multipole atomic integrals: + +{% $$ \langle \chi_i | x | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | y | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | z | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} + +*) + +open Lacaml.D + +type t + +val matrix_x : t -> Mat.t +(** {% $$ \langle \chi_i | x | \chi_j \rangle $$ %} *) + +val matrix_y : t -> Mat.t +(** {% $$ \langle \chi_i | y | \chi_j \rangle $$ %} *) + +val matrix_z : t -> Mat.t +(** {% $$ \langle \chi_i | z | \chi_j \rangle $$ %} *) + +val matrix_x2 : t -> Mat.t +(** {% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %} *) + +val matrix_y2 : t -> Mat.t +(** {% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %} *) + +val matrix_z2 : t -> Mat.t +(** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *) + +val of_basis : Basis.t -> t diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index f73dfa0..0eb858c 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -18,77 +18,77 @@ module Psp = PrimitiveShellPair let cutoff = integrals_cutoff -let to_powers x = +let to_powers x = let open Zkey in match to_powers x with | Six x -> x | _ -> assert false (** Computes all the overlap integrals of the contracted shell pair *) -let contracted_class shell_a shell_b : float Zmap.t = +let contracted_class shell_a shell_b : float Zmap.t = match Csp.make shell_a shell_b with - | Some shell_p -> + | None -> Zmap.create 0 + | Some shell_p -> begin - (* Pre-computation of integral class indices *) - let class_indices = Csp.zkey_array shell_p in + (* Pre-computation of integral class indices *) + let class_indices = Csp.zkey_array shell_p in - let contracted_class = - Array.make (Array.length class_indices) 0. - in + let contracted_class = + Array.make (Array.length class_indices) 0. + in - let a_minus_b = - Csp.a_minus_b shell_p - in - let norm_coef_scales = - Csp.norm_scales shell_p - in + let a_minus_b = + Csp.a_minus_b shell_p + in + let norm_coef_scales = + Csp.norm_scales shell_p + in - (* Compute all integrals in the shell for each pair of significant shell pairs *) + (* Compute all integrals in the shell for each pair of significant shell pairs *) - let xyz_of_int k = - match k with - | 0 -> Co.X - | 1 -> Co.Y - | _ -> Co.Z - in + let xyz_of_int k = + match k with + | 0 -> Co.X + | 1 -> Co.Y + | _ -> Co.Z + in - List.iter (fun (coef_prod, psp) -> - (** Screening on thr product of coefficients *) - if (abs_float coef_prod) > 1.e-3*.cutoff then - begin - let expo_inv = Psp.exponent_inv psp - and center_pa = Psp.center_minus_a psp - in - - Array.iteri (fun i key -> - let (angMomA,angMomB) = to_powers key in - let f k = - let xyz = xyz_of_int k in - Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB) - expo_inv - (Co.get xyz a_minus_b, - Co.get xyz center_pa) + List.iter (fun (coef_prod, psp) -> + (** Screening on the product of coefficients *) + if (abs_float coef_prod) > 1.e-6*.cutoff then + begin + let expo_inv = Psp.exponent_inv psp + and center_pa = Psp.center_minus_a psp in - let norm = norm_coef_scales.(i) in - let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in - contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral - ) class_indices - end - ) (Csp.coefs_and_shell_pairs shell_p); - let result = - Zmap.create (Array.length contracted_class) - in - Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; - result - end - | None -> Zmap.create 0 + + Array.iteri (fun i key -> + let (angMomA,angMomB) = to_powers key in + let f k = + let xyz = xyz_of_int k in + Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB) + expo_inv + (Co.get xyz a_minus_b, + Co.get xyz center_pa) + in + let norm = norm_coef_scales.(i) in + let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in + contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral + ) class_indices + end + ) (Csp.coefs_and_shell_pairs shell_p); + let result = + Zmap.create (Array.length contracted_class) + in + Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; + result + end (** Create overlap matrix *) let of_basis basis = - let to_powers x = + let to_powers x = let open Zkey in match to_powers x with | Three x -> x @@ -104,7 +104,7 @@ let of_basis basis = for i=0 to j do (* Compute all the integrals of the class *) let cls = - contracted_class shell.(i) shell.(j) + contracted_class shell.(i) shell.(j) in Array.iteri (fun j_c powers_j -> @@ -113,11 +113,11 @@ let of_basis basis = Array.iteri (fun i_c powers_i -> let i_c = Cs.index shell.(i) + i_c + 1 in let xi = to_powers powers_i in - let key = + let key = Zkey.of_powers_six xi xj in - let value = - try Zmap.find cls key + let value = + try Zmap.find cls key with Not_found -> 0. in result.{i_c,j_c} <- value; @@ -132,7 +132,7 @@ let of_basis basis = (** Create mixed overlap matrix *) let of_basis_pair first_basis second_basis = - let to_powers x = + let to_powers x = let open Zkey in match to_powers x with | Three x -> x @@ -140,7 +140,7 @@ let of_basis_pair first_basis second_basis = in let n = Bs.size first_basis - and m = Bs.size second_basis + and m = Bs.size second_basis and first = Bs.contracted_shells first_basis and second = Bs.contracted_shells second_basis in @@ -150,7 +150,7 @@ let of_basis_pair first_basis second_basis = for i=0 to (Array.length first) - 1 do (* Compute all the integrals of the class *) let cls = - contracted_class first.(i) second.(j) + contracted_class first.(i) second.(j) in Array.iteri (fun j_c powers_j -> @@ -159,11 +159,11 @@ let of_basis_pair first_basis second_basis = Array.iteri (fun i_c powers_i -> let i_c = Cs.index first.(i) + i_c + 1 in let xi = to_powers powers_i in - let key = + let key = Zkey.of_powers_six xi xj in - let value = - try Zmap.find cls key + let value = + try Zmap.find cls key with Not_found -> 0. in result.{i_c,j_c} <- value; @@ -180,7 +180,7 @@ let to_file ~filename overlap = let oc = open_out filename in let n = - Mat.dim1 overlap + Mat.dim1 overlap in for j=1 to n do diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index 31a4ebd..33da7e4 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -3,11 +3,7 @@ open Xyz_ast type t = (Element.t * Coordinate.t) array -let of_xyz_file filename = - let lexbuf = - let ic = open_in filename in - Lexing.from_channel ic - in +let of_xyz_lexbuf lexbuf = let data = Xyz_parser.input Nuclei_lexer.read_all lexbuf in @@ -24,6 +20,30 @@ let of_xyz_file filename = |> Array.of_list +let of_xyz_string buffer = + Zmatrix.of_string buffer + |> Zmatrix.to_xyz + |> Array.map (fun (e,x,y,z) -> + (e, Coordinate.(angstrom_to_bohr @@ make_angstrom { x ; y ; z} )) + ) + + +let of_xyz_string input_string = + Lexing.from_string input_string + |> of_xyz_lexbuf + + +let of_xyz_file filename = + let ic = open_in filename in + let lexbuf = + Lexing.from_channel ic + in + let result = + of_xyz_lexbuf lexbuf + in + close_in ic; + result + let of_zmt_string buffer = Zmatrix.of_string buffer diff --git a/Nuclei/Nuclei.mli b/Nuclei/Nuclei.mli index 4d216cb..9944784 100644 --- a/Nuclei/Nuclei.mli +++ b/Nuclei/Nuclei.mli @@ -5,6 +5,9 @@ of tuples ({!Element.t}, {!Coordinate.t}). type t = (Element.t * Coordinate.t) array +val of_xyz_string : string -> t +(** Create from a string, in [xyz] format. *) + val of_xyz_file : string -> t (** Create from a file, in [xyz] format. *)