From 59980ed51ca39225bae4212d40bfe0505bbd5c3e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 7 Dec 2020 13:52:56 +0100 Subject: [PATCH] Add MP2 --- ao/lib/basis_gaussian.mli | 3 +- linear_algebra/lib/matrix.mli | 6 +++- perturbation/lib/dune | 10 ++++++ perturbation/lib/mp2.ml | 62 +++++++++++++++++++++++++++++++++++ perturbation/test/dune | 13 ++++++++ perturbation/test/mp2.ml | 27 +++++++++++++++ qcaml/lib/dune | 1 + qcaml/lib/qcaml.ml | 1 + test/dune | 1 + test/run_tests.ml | 5 +-- 10 files changed, 125 insertions(+), 4 deletions(-) create mode 100644 perturbation/lib/dune create mode 100644 perturbation/lib/mp2.ml create mode 100644 perturbation/test/dune create mode 100644 perturbation/test/mp2.ml diff --git a/ao/lib/basis_gaussian.mli b/ao/lib/basis_gaussian.mli index dfc2630..e4849ef 100644 --- a/ao/lib/basis_gaussian.mli +++ b/ao/lib/basis_gaussian.mli @@ -55,4 +55,5 @@ val values : t -> Coordinate.t -> Gaussian.Basis.t Vector.t val make : basis:Gaussian.Basis.t -> ?operators:Operator.t list -> ?cartesian:bool -> Nuclei.t -> t (** Creates the data structure for atomic orbitals from a Gaussian basis and the - molecular geometry {Nuclei.t} *) + molecular geometry {Nuclei.t}. +*) diff --git a/linear_algebra/lib/matrix.mli b/linear_algebra/lib/matrix.mli index 6d4011e..c99c0f8 100644 --- a/linear_algebra/lib/matrix.mli +++ b/linear_algebra/lib/matrix.mli @@ -9,7 +9,11 @@ val dim2: ('a,'b) t -> int (** Second dimension of the matrix *) val make : int -> int -> float -> ('a,'b) t -(** Creates a matrix initialized with the given value. *) +(** Creates a matrix initialized with the given value. + @param m: First dimension + @param n: Seconfd dimension + @param f: Value used to initialize the matrix elements +*) val make0 : int -> int -> ('a,'b) t (** Creates a zero-initialized matrix. *) diff --git a/perturbation/lib/dune b/perturbation/lib/dune new file mode 100644 index 0000000..c6acb7a --- /dev/null +++ b/perturbation/lib/dune @@ -0,0 +1,10 @@ +; name = name of the supermodule that will wrap all source files as submodules +; public_name = name of the library for ocamlfind and opam +(library + (name perturbation) + (public_name qcaml.perturbation) + (libraries + qcaml.simulation + qcaml.mo + ) + (synopsis "Perturbation theory.")) diff --git a/perturbation/lib/mp2.ml b/perturbation/lib/mp2.ml new file mode 100644 index 0000000..443fe7f --- /dev/null +++ b/perturbation/lib/mp2.ml @@ -0,0 +1,62 @@ +open Linear_algebra + +type t = float + +let make ~frozen_core hf = + let mo_basis = + Mo.Basis.of_hartree_fock hf + in + let epsilon = + Mo.Basis.mo_energies mo_basis + in + let mo_class = + Mo.Class.cas_sd mo_basis ~frozen_core 0 0 + |> Mo.Class.to_list + in + let eri = + Mo.Basis.ee_ints mo_basis + in + let inactives = + List.filter (fun i -> + match i with Mo.Class.Inactive _ -> true | _ -> false) mo_class + and virtuals = + List.filter (fun i -> + match i with Mo.Class.Virtual _ -> true | _ -> false) mo_class + in + + let rmp2 () = + List.fold_left (fun accu b -> + match b with Mo.Class.Virtual b -> + let eps = -. (epsilon%.(b)) in + accu +. + List.fold_left (fun accu a -> + match a with Mo.Class.Virtual a -> + let eps = eps -. (epsilon%.(a)) in + accu +. + List.fold_left (fun accu j -> + match j with Mo.Class.Inactive j -> + let eps = eps +. epsilon%.(j) in + accu +. + List.fold_left (fun accu i -> + match i with Mo.Class.Inactive i -> + let eps = eps +. epsilon%.(i) in + let ijab = Four_idx_storage.get_phys eri i j a b + and abji = Four_idx_storage.get_phys eri a b j i in + let abij = ijab in + accu +. ijab *. ( abij +. abij -. abji) /. eps + | _ -> accu + ) 0. inactives + | _ -> accu + ) 0. inactives + | _ -> accu + ) 0. virtuals + | _ -> accu + ) 0. virtuals + in + + + match Mo.Hartree_fock.kind hf with + | Mo.Hartree_fock.RHF -> rmp2 () + | _ -> failwith "Not implemented" + + diff --git a/perturbation/test/dune b/perturbation/test/dune new file mode 100644 index 0000000..78a18e1 --- /dev/null +++ b/perturbation/test/dune @@ -0,0 +1,13 @@ +(library + (name test_perturbation) + (libraries + alcotest + qcaml.mo + qcaml.perturbation + ) + (synopsis "Tests for the perturbation")) + + + + + diff --git a/perturbation/test/mp2.ml b/perturbation/test/mp2.ml new file mode 100644 index 0000000..d2ffc80 --- /dev/null +++ b/perturbation/test/mp2.ml @@ -0,0 +1,27 @@ +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); + let e_mp2 = Perturbation.Mp2.make ~frozen_core:true hf in + Printf.printf "%s\n" (string_of_float e_mp2); + check (float 1.e-10) "MP2" (-0.201621141207) (e_mp2) + ] + diff --git a/qcaml/lib/dune b/qcaml/lib/dune index 0656da0..65b97eb 100644 --- a/qcaml/lib/dune +++ b/qcaml/lib/dune @@ -11,6 +11,7 @@ qcaml.mo qcaml.operators qcaml.particles + qcaml.perturbation qcaml.simulation ) (synopsis "Main QCaml entry point")) diff --git a/qcaml/lib/qcaml.ml b/qcaml/lib/qcaml.ml index 13dcd04..cbf0c0c 100644 --- a/qcaml/lib/qcaml.ml +++ b/qcaml/lib/qcaml.ml @@ -6,4 +6,5 @@ module Linear_algebra = Linear_algebra module Mo = Mo module Operators = Operators module Particles = Particles +module Perturbation = Perturbation module Simulation = Simulation diff --git a/test/dune b/test/dune index e0c8096..ccad83c 100644 --- a/test/dune +++ b/test/dune @@ -9,6 +9,7 @@ test_gaussian_integrals test_ao_basis test_mo_basis + test_perturbation )) (alias diff --git a/test/run_tests.ml b/test/run_tests.ml index 55b8f39..af83b15 100644 --- a/test/run_tests.ml +++ b/test/run_tests.ml @@ -11,8 +11,9 @@ let test_suites: unit Alcotest.test list = [ "Gaussian_basis.General_basis", Test_gaussian_basis.General_basis.tests; "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; + "Mo.Guess", Test_mo_basis.Guess.tests; + "Mo.Hartree_Fock", Test_mo_basis.Hf.tests; + "Perturbation.Mp2", Test_perturbation.Mp2.tests; ] let () = Alcotest.run "QCaml" test_suites