From b3383ac523b74d10a9db51adccdc391fcefad4c6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Mar 2020 17:58:50 +0100 Subject: [PATCH] Fixed overflow in Boys --- Notebooks/SpVec.ipynb | 2 +- Notebooks/TemplateQCaml.ipynb | 159 ++++++++++++++++++++++++++++++++++ Utils/Util.ml | 17 ++-- 3 files changed, 167 insertions(+), 11 deletions(-) create mode 100644 Notebooks/TemplateQCaml.ipynb diff --git a/Notebooks/SpVec.ipynb b/Notebooks/SpVec.ipynb index 72b1b75..4070396 100644 --- a/Notebooks/SpVec.ipynb +++ b/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" }, diff --git a/Notebooks/TemplateQCaml.ipynb b/Notebooks/TemplateQCaml.ipynb new file mode 100644 index 0000000..1f61dd1 --- /dev/null +++ b/Notebooks/TemplateQCaml.ipynb @@ -0,0 +1,159 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "val png_image : string -> unit = \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 = \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 +} diff --git a/Utils/Util.ml b/Utils/Util.ml index cdc84fc..4309fe1 100644 --- a/Utils/Util.ml +++ b/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 - (* - assert (abs_float t > 1.e-10); - *) - if t <> 0. then - begin + let power_t_inv = (maxm+maxm+1) in + try 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 - end; - result + done; + result + with Exit -> result end