Add MP2

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}. 
+*)

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. *)

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."))

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"
+
+

perturbation/test/dune

@@ -0,0 +1,13 @@
+(library
+ (name test_perturbation)
+ (libraries
+  alcotest
+  qcaml.mo
+  qcaml.perturbation
+ )
+ (synopsis "Tests for the perturbation"))
+
+
+
+
+

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)
+  ]
+    

qcaml/lib/dune

@@ -11,6 +11,7 @@
 	qcaml.mo
 	qcaml.operators
         qcaml.particles
+        qcaml.perturbation
 	qcaml.simulation
  )
  (synopsis "Main QCaml entry point"))

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

test/dune

@@ -9,6 +9,7 @@
   test_gaussian_integrals
   test_ao_basis
   test_mo_basis
+  test_perturbation
 ))
 
 (alias

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.Guess", Test_mo_basis.Guess.tests;
-  "Mo_basis.Hartree_Fock", Test_mo_basis.Hf.tests;
+  "Mo.Hartree_Fock", Test_mo_basis.Hf.tests;
+  "Perturbation.Mp2", Test_perturbation.Mp2.tests;
 ]
 
 let () = Alcotest.run "QCaml" test_suites