10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 20:33:36 +01:00

Added a test for Hartree-Fock

This commit is contained in:
Anthony Scemama 2020-10-26 11:33:51 +01:00
parent 718d39924f
commit 2a7a586c92
9 changed files with 232 additions and 18 deletions

View File

@ -1,6 +1,7 @@
(executable (executables
(name ex_integrals) (names
(libraries ex_integrals
qcaml ex_hartree_fock
)) )
(libraries qcaml))

View File

@ -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 "<string>";
doc="Name of the file containing the basis set"; } ;
{ short='x' ; long="xyz" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
{ short='m' ; long="multiplicity" ; opt=Optional;
arg=With_arg "<int>";
doc="Spin multiplicity (2S+1). Default is singlet"; } ;
{ short='c' ; long="charge" ; opt=Optional;
arg=With_arg "<int>";
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 *)

View File

@ -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 "<string>";
doc="Name of the file containing the basis set"; } ;
{ short='x' ; long="xyz" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
{ short='m' ; long="multiplicity" ; opt=Optional;
arg=With_arg "<int>";
doc="Spin multiplicity (2S+1). Default is singlet"; } ;
{ short='c' ; long="charge" ; opt=Optional;
arg=With_arg "<int>";
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

View File

@ -3,7 +3,7 @@
#+PROPERTY #+PROPERTY
In this example, we write a program that reads the geometry of a 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- format. The output is a set of files containing the one- and two-
electron integrals. electron integrals.

View File

@ -176,17 +176,17 @@ let gemm ?m ?n ?k ?beta ?c ?(transa=`N) ?(alpha=1.0) a ?(transb=`N) b =
in in
gemm ?m ?n ?k ?beta ?c ~transa ~alpha a ~transb b 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 = 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 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 = 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 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 = 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 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 = 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 gemm ?m ?n ?k ?beta ?c ~alpha a b ~transa:`T ~transb:`T
let gemm_trace ?(transa=`N) a ?(transb=`N) b = let gemm_trace ?(transa=`N) a ?(transb=`N) b =

24
mo/test/hf.ml Normal file
View File

@ -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);
]

View File

@ -4,11 +4,13 @@
(name qcaml) (name qcaml)
(public_name qcaml) (public_name qcaml)
(libraries (libraries
qcaml.ao
qcaml.common qcaml.common
qcaml.particles
qcaml.gaussian qcaml.gaussian
qcaml.gaussian_integrals qcaml.gaussian_integrals
qcaml.mo
qcaml.operators qcaml.operators
qcaml.ao qcaml.particles
qcaml.simulation
) )
(synopsis "Main QCaml entry point")) (synopsis "Main QCaml entry point"))

View File

@ -1,8 +1,9 @@
module Common = Common
module Ao = Ao module Ao = Ao
module Common = Common
module Gaussian = Gaussian module Gaussian = Gaussian
module Gaussian_integrals = Gaussian_integrals module Gaussian_integrals = Gaussian_integrals
module Linear_algebra = Linear_algebra module Linear_algebra = Linear_algebra
module Mo = Mo
module Operators = Operators module Operators = Operators
module Particles = Particles module Particles = Particles
module Simulation = Simulation

View File

@ -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_gaussian", Test_ao_basis.Ao_basis_gaussian.tests;
"Ao_basis.Ao_basis", Test_ao_basis.Ao_basis.tests; "Ao_basis.Ao_basis", Test_ao_basis.Ao_basis.tests;
"Mo_basis.Guess", Test_mo_basis.Guess.tests; "Mo_basis.Guess", Test_mo_basis.Guess.tests;
"Mo_basis.Hartree_Fock", Test_mo_basis.Hf.tests;
] ]
let () = Alcotest.run "QCaml" test_suites let () = Alcotest.run "QCaml" test_suites