Fixed overflow in Boys

Notebooks/SpVec.ipynb

@@ -1746,7 +1746,7 @@
  "metadata": {
   "celltoolbar": "Raw Cell Format",
   "kernelspec": {
-   "display_name": "OCaml 4.07.1",
+   "display_name": "OCaml default",
    "language": "OCaml",
    "name": "ocaml-jupyter"
   },

Notebooks/TemplateQCaml.ipynb

@@ -0,0 +1,159 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "val png_image : string -> unit = <fun>\n"
+      ]
+     },
+     "execution_count": 4,
+     "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"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "val png_image : string -> Jupyter_notebook.display_id = <fun>\n"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "let png_image = print_endline ;;\n",
+    "\n",
+    "(* --------- *)\n",
+    "\n",
+    "#cd \"/home/scemama/QCaml\";;\n",
+    "#use \"topfind\";;\n",
+    "#require \"jupyter.notebook\";;\n",
+    "\n",
+    "let png_image name = \n",
+    "    Jupyter_notebook.display_file ~base64:true \"image/png\" (\"Notebooks/images/\"^name)\n",
+    ";;\n",
+    "\n",
+    "#require \"lacaml.top\";;\n",
+    "#require \"alcotest\";;\n",
+    "#require \"str\";;\n",
+    "#require \"bigarray\";;\n",
+    "#require \"zarith\";;\n",
+    "#require \"getopt\";;\n",
+    "#directory \"_build\";;\n",
+    "#directory \"_build/Basis\";;\n",
+    "#directory \"_build/CI\";;\n",
+    "#directory \"_build/MOBasis\";;\n",
+    "#directory \"_build/Nuclei\";;\n",
+    "#directory \"_build/Parallel\";;\n",
+    "#directory \"_build/Perturbation\";;\n",
+    "#directory \"_build/SCF\";;\n",
+    "#directory \"_build/Utils\";;\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "File Constants.cmo is not a bytecode object file.\n",
+      "File Util.cma is not a bytecode object file.\n",
+      "File Matrix.cmo is not a bytecode object file.\n",
+      "File Simulation.cmo is not a bytecode object file.\n",
+      "File Spindeterminant.cmo is not a bytecode object file.\n",
+      "File Determinant.cmo is not a bytecode object file.\n",
+      "File HartreeFock.cmo is not a bytecode object file.\n",
+      "File MOBasis.cmo is not a bytecode object file.\n",
+      "File F12CI.cmo is not a bytecode object file.\n"
+     ]
+    }
+   ],
+   "source": [
+    "#load \"Constants.cmo\";;\n",
+    "#load_rec \"Util.cma\";;\n",
+    "#load_rec \"Matrix.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\";;"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#install_printer AngularMomentum.pp_string ;;\n",
+    "#install_printer Basis.pp ;;\n",
+    "#install_printer Charge.pp ;;\n",
+    "#install_printer Coordinate.pp ;;\n",
+    "#install_printer Vector.pp;;\n",
+    "#install_printer Matrix.pp;;\n",
+    "#install_printer Util.pp_float_2darray;;\n",
+    "#install_printer Util.pp_float_array;;\n",
+    "#install_printer Util.pp_matrix;;\n",
+    "#install_printer HartreeFock.pp ;;\n",
+    "#install_printer Fock.pp ;;\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",
+    "(*\n",
+    "#install_printer DeterminantSpace.pp;;\n",
+    "*)\n",
+    "#install_printer SpindeterminantSpace.pp;;\n",
+    "#install_printer Bitstring.pp;;\n",
+    "\n",
+    "(* --------- *)\n",
+    "\n",
+    "open Lacaml.D\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "OCaml default",
+   "language": "OCaml",
+   "name": "ocaml-jupyter"
+  },
+  "language_info": {
+   "codemirror_mode": "text/x-ocaml",
+   "file_extension": ".ml",
+   "mimetype": "text/x-ocaml",
+   "name": "OCaml",
+   "nbconverter_exporter": null,
+   "pygments_lexer": "OCaml",
+   "version": "4.07.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

Utils/Util.ml

@@ -209,29 +209,26 @@ let boys_function ~maxm t =
       let result =
         Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1))
       in
-      (*
+      let power_t_inv = (maxm+maxm+1)  in
-      assert (abs_float t > 1.e-10);
+      try  
-      *)
-      if t <> 0. then
-        begin
           let fmax = 
             let t_inv = sqrt (1. /. t) in
             let n = float_of_int maxm in
             let dm = 0.5 +. n in
-            let f  = (pow t_inv (maxm+maxm+1) ) in
+            let f  = (pow t_inv power_t_inv ) in
             match classify_float f with
               | FP_normal -> (incomplete_gamma dm t) *. 0.5 *. f
               | FP_zero
               | FP_subnormal -> 0.
-              | _ -> invalid_arg "zero_m overflow"
+              | _ -> raise Exit
           in
           let emt = exp (-. t) in
           result.(maxm) <- fmax;
           for n=maxm-1 downto 0 do
             result.(n) <- ( (t+.t) *. result.(n+1) +. emt) *. result.(n)
-          done
+          done;
-        end;
+          result
-      result
+      with Exit -> result
     end