From 2a7a586c9212b48e5e65ba7af5d8e8a556957f05 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 26 Oct 2020 11:33:51 +0100 Subject: [PATCH] Added a test for Hartree-Fock --- examples/dune | 11 ++-- examples/ex_hartree_fock.ml | 76 ++++++++++++++++++++++++ examples/ex_hartree_fock.org | 109 +++++++++++++++++++++++++++++++++++ examples/ex_integrals.org | 2 +- linear_algebra/lib/matrix.ml | 16 ++--- mo/test/hf.ml | 24 ++++++++ qcaml/lib/dune | 6 +- qcaml/lib/qcaml.ml | 5 +- test/run_tests.ml | 1 + 9 files changed, 232 insertions(+), 18 deletions(-) create mode 100644 examples/ex_hartree_fock.ml create mode 100644 examples/ex_hartree_fock.org create mode 100644 mo/test/hf.ml diff --git a/examples/dune b/examples/dune index accdf4d..3e5d340 100644 --- a/examples/dune +++ b/examples/dune @@ -1,6 +1,7 @@ -(executable - (name ex_integrals) - (libraries - qcaml -)) +(executables + (names + ex_integrals + ex_hartree_fock +) +(libraries qcaml)) diff --git a/examples/ex_hartree_fock.ml b/examples/ex_hartree_fock.ml new file mode 100644 index 0000000..76fc8cf --- /dev/null +++ b/examples/ex_hartree_fock.ml @@ -0,0 +1,76 @@ +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Header][Header:1]] *) +module Command_line = Qcaml.Common.Command_line +module Util = Qcaml.Common.Util + +let () = +(* Header:1 ends here *) + +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Definition][Definition:1]] *) +let open Command_line in +begin + set_header_doc (Sys.argv.(0) ^ " - QuAcK command"); + set_description_doc "Computes the one- and two-electron hartree_fock on the Gaussian atomic basis set."; + set_specs + [ { short='b' ; long="basis" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the basis set"; } ; + + { short='x' ; long="xyz" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the nuclear coordinates in xyz format"; } ; + + { short='m' ; long="multiplicity" ; opt=Optional; + arg=With_arg ""; + doc="Spin multiplicity (2S+1). Default is singlet"; } ; + + { short='c' ; long="charge" ; opt=Optional; + arg=With_arg ""; + doc="Total charge of the molecule. Specify negative charges with 'm' instead of the minus sign, for example m1 instead of -1. Default is 0"; } ; + + ] +end; +(* Definition:1 ends here *) + +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Interpretation][Interpretation:1]] *) +let basis_file = Util.of_some @@ Command_line.get "basis" in +let nuclei_file = Util.of_some @@ Command_line.get "xyz" in + +let charge = + match Command_line.get "charge" with + | Some x -> ( if x.[0] = 'm' then + ~- (int_of_string (String.sub x 1 (String.length x - 1))) + else + int_of_string x ) + | None -> 0 +in + +let multiplicity = + match Command_line.get "multiplicity" with + | Some x -> int_of_string x + | None -> 1 +in +(* Interpretation:1 ends here *) + +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:1]] *) +let nuclei = + Qcaml.Particles.Nuclei.of_xyz_file nuclei_file +in +(* Computation:1 ends here *) + +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:2]] *) +let ao_basis = + Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file +in +(* Computation:2 ends here *) + +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:3]] *) +let simulation = Qcaml.Simulation.make ~charge ~multiplicity ~nuclei ao_basis in +(* Computation:3 ends here *) + +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:4]] *) +let hf = Qcaml.Mo.Hartree_fock.make ~guess:`Huckel simulation in +(* Computation:4 ends here *) + +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Output][Output:1]] *) +Format.printf "@[%a@]" (Mo.Hartree_fock.pp) hf +(* Output:1 ends here *) diff --git a/examples/ex_hartree_fock.org b/examples/ex_hartree_fock.org new file mode 100644 index 0000000..3887973 --- /dev/null +++ b/examples/ex_hartree_fock.org @@ -0,0 +1,109 @@ +#+TITLE: Hartree-Fock + +#+PROPERTY + +In this example, we write a program that makes a Hartree-Fock +caculation. The molecule is read in =xyz= format and a Gaussian +atomic basis set in GAMESS format. + +* Header + +#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml +module Command_line = Qcaml.Common.Command_line +module Util = Qcaml.Common.Util + +let () = +#+END_SRC + +* Command-line arguments + + We use the =Command_line= module to define the following possible + arguments: + - =-b --basis= : The name of the file containing the basis set + - =-x --xyz= : The name of the file containing the atomic coordinates + - =-c --charge= : The charge of the molecule + - =-m --multiplicity= : The spin multiplicity + +** Definition + + #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml +let open Command_line in +begin + set_header_doc (Sys.argv.(0) ^ " - QuAcK command"); + set_description_doc "Computes the one- and two-electron hartree_fock on the Gaussian atomic basis set."; + set_specs + [ { short='b' ; long="basis" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the basis set"; } ; + + { short='x' ; long="xyz" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the nuclear coordinates in xyz format"; } ; + + { short='m' ; long="multiplicity" ; opt=Optional; + arg=With_arg ""; + doc="Spin multiplicity (2S+1). Default is singlet"; } ; + + { short='c' ; long="charge" ; opt=Optional; + arg=With_arg ""; + doc="Total charge of the molecule. Specify negative charges with 'm' instead of the minus sign, for example m1 instead of -1. Default is 0"; } ; + + ] +end; + #+END_SRC + +** Interpretation + + #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml +let basis_file = Util.of_some @@ Command_line.get "basis" in +let nuclei_file = Util.of_some @@ Command_line.get "xyz" in + +let charge = + match Command_line.get "charge" with + | Some x -> ( if x.[0] = 'm' then + ~- (int_of_string (String.sub x 1 (String.length x - 1))) + else + int_of_string x ) + | None -> 0 +in + +let multiplicity = + match Command_line.get "multiplicity" with + | Some x -> int_of_string x + | None -> 1 +in + #+END_SRC + +* Computation + We first read the =xyz= file to create a molecule: + #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml +let nuclei = + Qcaml.Particles.Nuclei.of_xyz_file nuclei_file +in + #+END_SRC + + Then we create an Gaussian AO basis using the atomic coordinates: + #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml +let ao_basis = + Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file +in + #+END_SRC + + We create a simulation from the nuclei and the basis set: + #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml +let simulation = Qcaml.Simulation.make ~nuclei ao_basis in + #+END_SRC + + and we can make the Hartree-Fock computation: + #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml +let hf = Qcaml.Mo.Hartree_fock.make ~guess:`Huckel simulation in + #+END_SRC + +* Output + + We print the convergence of the calculation: + #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml +Format.printf "@[%a@]" (Mo.Hartree_fock.pp) hf + #+END_SRC + + diff --git a/examples/ex_integrals.org b/examples/ex_integrals.org index 0ce91cf..241c3b6 100644 --- a/examples/ex_integrals.org +++ b/examples/ex_integrals.org @@ -3,7 +3,7 @@ #+PROPERTY In this example, we write a program that reads the geometry of a -molecule in in =xyz= format and a Gaussian atomic basis set in GAMESS +molecule in =xyz= format and a Gaussian atomic basis set in GAMESS format. The output is a set of files containing the one- and two- electron integrals. diff --git a/linear_algebra/lib/matrix.ml b/linear_algebra/lib/matrix.ml index f01f700..b30c853 100644 --- a/linear_algebra/lib/matrix.ml +++ b/linear_algebra/lib/matrix.ml @@ -176,17 +176,17 @@ let gemm ?m ?n ?k ?beta ?c ?(transa=`N) ?(alpha=1.0) a ?(transb=`N) b = in gemm ?m ?n ?k ?beta ?c ~transa ~alpha a ~transb b -let gemm_nn ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b = - gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`N ~transb:`N +let gemm_nn ?m ?n ?k ?beta ?c ?(alpha=1.0) a b = + gemm ?m ?n ?k ?beta ?c ~alpha a b ~transa:`N ~transb:`N -let gemm_nt ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b = - gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`N ~transb:`T +let gemm_nt ?m ?n ?k ?beta ?c ?(alpha=1.0) a b = + gemm ?m ?n ?k ?beta ?c ~alpha a b ~transa:`N ~transb:`T -let gemm_tn ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b = - gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`T ~transb:`N +let gemm_tn ?m ?n ?k ?beta ?c ?(alpha=1.0) a b = + gemm ?m ?n ?k ?beta ?c ~alpha a b ~transa:`T ~transb:`N -let gemm_tt ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b = - gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`T ~transb:`T +let gemm_tt ?m ?n ?k ?beta ?c ?(alpha=1.0) a b = + gemm ?m ?n ?k ?beta ?c ~alpha a b ~transa:`T ~transb:`T let gemm_trace ?(transa=`N) a ?(transb=`N) b = diff --git a/mo/test/hf.ml b/mo/test/hf.ml new file mode 100644 index 0000000..72925b0 --- /dev/null +++ b/mo/test/hf.ml @@ -0,0 +1,24 @@ +open Alcotest +open Particles + +let wd = Common.Qcaml.root ^ Filename.dir_sep ^ "test" + +let tests = + [ "HF Water", `Quick, fun () -> + let nuclei = + wd ^ Filename.dir_sep ^ "water.xyz" + |> Nuclei.of_xyz_file + in + let basis_filename = + wd ^ Filename.dir_sep ^ "cc-pvdz" + in + let ao_basis = + Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian + ~cartesian:false ~nuclei basis_filename + in + let simulation = Simulation.make ~nuclei ao_basis in + let hf = Mo.Hartree_fock.make ~guess:`Huckel simulation in + Format.printf "%a" (Mo.Hartree_fock.pp) hf; + check (float 1.e-10) "Energy" (-76.0267987006) (Mo.Hartree_fock.energy hf); + ] + diff --git a/qcaml/lib/dune b/qcaml/lib/dune index 9453e6b..0656da0 100644 --- a/qcaml/lib/dune +++ b/qcaml/lib/dune @@ -4,11 +4,13 @@ (name qcaml) (public_name qcaml) (libraries + qcaml.ao qcaml.common - qcaml.particles qcaml.gaussian qcaml.gaussian_integrals + qcaml.mo qcaml.operators - qcaml.ao + qcaml.particles + qcaml.simulation ) (synopsis "Main QCaml entry point")) diff --git a/qcaml/lib/qcaml.ml b/qcaml/lib/qcaml.ml index 8784bd8..13dcd04 100644 --- a/qcaml/lib/qcaml.ml +++ b/qcaml/lib/qcaml.ml @@ -1,8 +1,9 @@ -module Common = Common module Ao = Ao +module Common = Common module Gaussian = Gaussian module Gaussian_integrals = Gaussian_integrals module Linear_algebra = Linear_algebra +module Mo = Mo module Operators = Operators module Particles = Particles - +module Simulation = Simulation diff --git a/test/run_tests.ml b/test/run_tests.ml index fecebd5..55b8f39 100644 --- a/test/run_tests.ml +++ b/test/run_tests.ml @@ -12,6 +12,7 @@ let test_suites: unit Alcotest.test list = [ "Ao_basis.Ao_basis_gaussian", Test_ao_basis.Ao_basis_gaussian.tests; "Ao_basis.Ao_basis", Test_ao_basis.Ao_basis.tests; "Mo_basis.Guess", Test_mo_basis.Guess.tests; + "Mo_basis.Hartree_Fock", Test_mo_basis.Hf.tests; ] let () = Alcotest.run "QCaml" test_suites