Added MP2

CI/CI.ml

@@ -1,6 +1,7 @@
 open Lacaml.D
 
 module Ds = DeterminantSpace
+module Sd = Spindeterminant
 
 type t =
   {
@@ -336,6 +337,9 @@ let create_matrix_spin f det_space =
 
 
 
+
+
+
 let make ?(n_states=1) det_space = 
 
   let m_H = 
@@ -380,3 +384,65 @@ let make ?(n_states=1) det_space =
   { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
 
 
+
+(*
+let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } =
+
+  let mo_indices =
+    let cls = MOClass.mo_class_array (DeterminantSpace.mo_class det_space) in
+    Util.list_range 1 (Ds.mo_basis det_space |> MOBasis.size)
+    |> List.filter (fun i -> match cls.(i) with
+      | MOClass.Deleted  _ 
+      | MOClass.Core     _ -> false
+      | _ -> true
+      )
+  in
+
+  let psi0, e0 = Lazy.force eigensystem in
+
+  let psi0 =
+    let stream = 
+      Ds.determinant_stream det_space
+    in
+    Array.init (Ds.size det_space) (fun i ->
+      (Stream.next stream, psi0.{i,1})
+  in
+  let e0 = e0.{1} in 
+
+  let det_contribution i =
+    let psi_filtered =
+      List.filter (fun (det, _) ->
+          Determinant.degree det i < 4) psi0
+    in
+    List.fold_left (fun accu spin1 ->
+        accu +.
+        List.fold_left (fun accu particle ->
+            accu +.
+            List.fold_left (fun accu hole ->
+                let alfa =
+                  Determinant.single_excitation spin1 hole particle i
+                in
+                if Determinant.is_none alfa then
+                  accu
+                else
+                  let psi_h_alfa =
+                    Array.fold_left (fun (det, coef) ->
+                        coef *. (h_ij det alfa)
+                      ) 0. psi_filtered
+                  in
+                  let h_aa = h_ij alfa alfa in
+                  accu +. psi_h_alfa *. psi_h_alfa / (e0 - h_aa)
+              ) 0. mo_indices
+          ) 0. mo_indices
+      ) 0. [ Spin.Alfa ; Spin.Beta ]
+
+  in
+
+  match det_space with
+  | Ds.Arbitrary -> assert false
+  | Ds.Spin alfa_dets beta_dets ->
+    Array.map ( fun alfa ->
+        Array.map ( fun beta ->
+          ) beta_dets
+      ) alfa_dets
+ *) 

CI/DeterminantSpace.ml

@@ -39,6 +39,25 @@ let n_beta   t = t.n_beta
 let mo_class t = t.mo_class
 let mo_basis t = t.mo_basis
 
+let active_mos t =
+  mo_class t 
+  |> MOClass.active_mos
+
+
+let inactive_mos t =
+  mo_class t 
+  |> MOClass.inactive_mos
+
+
+let virtual_mos t =
+  mo_class t 
+  |> MOClass.inactive_mos
+
+
+let mo_class_array t =
+  mo_class t
+  |> MOClass.mo_class_array 
+
 
 let size t = 
   match t.determinants with

MOBasis/MOClass.ml

@@ -19,6 +19,9 @@ let pp_mo_occ ppf = function
 
 let of_list t = t
 
+let to_list t = t
+
+
 let core_mos t =
   List.map (fun x ->
       match x with 
@@ -26,6 +29,7 @@ let core_mos t =
       | _ -> None) t
   |> Util.list_some
 
+
 let inactive_mos t =
   List.map (fun x ->
       match x with 
@@ -33,6 +37,7 @@ let inactive_mos t =
       | _ -> None ) t
   |> Util.list_some
 
+
 let active_mos t =
   List.map (fun x ->
       match x with 
@@ -40,6 +45,7 @@ let active_mos t =
       | _ -> None ) t
   |> Util.list_some
 
+
 let virtual_mos t =
   List.map (fun x ->
       match x with 
@@ -47,6 +53,7 @@ let virtual_mos t =
       | _ -> None ) t
   |> Util.list_some
 
+
 let deleted_mos t =
   List.map (fun x ->
       match x with 
@@ -55,6 +62,20 @@ let deleted_mos t =
   |> Util.list_some
 
 
+let mo_class_array t =
+  let sze = List.length t + 1 in
+  let result = Array.make sze (Deleted 0) in
+  List.iter (fun c ->
+      match c with
+      |  Core      i  ->  result.(i)  <-  Core      i
+      |  Inactive  i  ->  result.(i)  <-  Inactive  i
+      |  Active    i  ->  result.(i)  <-  Active    i
+      |  Virtual   i  ->  result.(i)  <-  Virtual   i
+      |  Deleted   i  ->  result.(i)  <-  Deleted   i
+    ) t;
+  result
+
+
 let fci ?(frozen_core=true) mo_basis =
   let mo_num = MOBasis.size mo_basis in
   let ncore  = (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 in
@@ -71,3 +92,25 @@ let fci ?(frozen_core=true) mo_basis =
       |> List.map (fun i -> Active i)
   )
 
+
+let cas_sd mo_basis n m =
+  let mo_num = MOBasis.size mo_basis in
+  let n_alfa = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_alfa in
+  let n_beta = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_alfa in
+  let n_unpaired = n_alfa - n_beta in
+  let last_inactive = n_alfa - (n + n_unpaired)/2 in
+  let last_active  = last_inactive + m in
+  let ncore  = (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 in
+  of_list (
+    List.concat [
+        Util.list_range 1 ncore
+        |> List.map (fun i -> Core i) ;
+        Util.list_range (ncore+1) last_inactive
+        |> List.map (fun i -> Inactive i) ;
+        Util.list_range (last_inactive+1) last_active
+        |> List.map (fun i -> Active i) ;
+        Util.list_range (last_active+1) mo_num
+        |> List.map (fun i -> Virtual i) 
+      ]
+  )
+

MOBasis/MOClass.mli

@@ -12,11 +12,19 @@ type t
 (** Creation *)
 val of_list : mo_class list -> t
 
+val to_list : t -> mo_class list 
+
 val fci : ?frozen_core:bool -> MOBasis.t -> t
 (** Creates the MO classes for FCI calculations : all [Active]. The 
     [n] lowest MOs are [Core] if [frozen_core = true].
 *)
 
+val cas_sd: MOBasis.t -> int -> int -> t
+(** [cas_sd mo_basis n m ] creates the MO classes for CAS(n,m) + SD
+    calculations.  lowest MOs are [Core], then all the next MOs are [Inactive],
+    then [Active], then [Virtual].
+*)
+
 
 val core_mos : t -> int list
 (** Returns a list containing the indices of the core MOs. *)
@@ -33,6 +41,10 @@ val inactive_mos : t -> int list
 val deleted_mos : t -> int list
 (** Returns a list containing the indices of the deleted MOs. *)
 
+val mo_class_array : t -> mo_class array
+(** Returns an array [a] such that [a.(i)] returns the class of MO [i].
+    As the MO indices start from [1], the array has an extra zero entry
+    that should be ignored. *)
 
 (** {2 Printers} *)
 

Makefile.include

@@ -1,6 +1,6 @@
 .NOPARALLEL:
 
-INCLUDE_DIRS=Parallel,Nuclei,Utils,Basis,SCF,MOBasis,CI,F12
+INCLUDE_DIRS=Parallel,Nuclei,Utils,Basis,SCF,MOBasis,CI,F12,Perturbation
 LIBS=
 PKGS=
 OCAMLBUILD=ocamlbuild -j 0 -cflags $(ocamlcflags) -lflags $(ocamllflags) $(ocamldocflags) -Is $(INCLUDE_DIRS) -ocamlopt $(ocamloptflags) $(mpi)

Perturbation/MP2.ml

@@ -0,0 +1,60 @@
+type t = float
+
+let make ?(frozen_core=true) hf =
+  let mo_basis = 
+    MOBasis.of_hartree_fock hf
+  in
+  let epsilon =
+    HartreeFock.eigenvalues hf
+  in
+  let mo_class =
+    MOClass.cas_sd mo_basis 0 0
+    |> MOClass.to_list
+  in
+  let eri = 
+    MOBasis.ee_ints mo_basis
+  in
+  let inactives =
+    List.filter (fun i ->
+      match i with MOClass.Inactive _ -> true | _ -> false) mo_class
+  and virtuals =
+    List.filter (fun i ->
+      match i with MOClass.Virtual _ -> true | _ -> false) mo_class
+  in
+
+  let rmp2 () = 
+    List.fold_left (fun accu b ->
+        match b with MOClass.Virtual b ->
+          let eps = -. epsilon.{b} in
+          accu +. 
+          List.fold_left (fun accu a ->
+              match a with MOClass.Virtual a ->
+                let eps = eps -. epsilon.{a} in
+                accu +. 
+                List.fold_left (fun accu j ->
+                    match j with MOClass.Inactive j ->
+                      let eps = eps +. epsilon.{j} in
+                      accu +. 
+                      List.fold_left (fun accu i ->
+                          match i with MOClass.Inactive i ->
+                            let eps = eps +. epsilon.{i} in
+                            let ijab = ERI.get_phys eri i j a b
+                            and abji = ERI.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 HartreeFock.kind hf with
+  | HartreeFock.RHF   -> rmp2 ()
+  | _ -> failwith "Not implemented"
+
+

SCF/HartreeFock.ml

@@ -280,7 +280,7 @@ let make
       | _ -> Fock.make_rhf  ~density:m_P ao_basis
     in
 
-    let m_F, m_Hc, m_J, m_K =
+    let m_F0, m_Hc, m_J, m_K =
       let x = fock in
       Fock.(fock x, core x, coulomb x, exchange x)
     in
@@ -291,7 +291,7 @@ let make
         gemm  m_S  m_C
       in
       gemm  m_SC (gemm  m_LSmo  m_SC  ~transb:`T)
-      |> Mat.add m_F
+      |> Mat.add m_F0
     in
 
 
@@ -322,10 +322,19 @@ let make
 
 
     (* MOs in orthogonal MO basis *)
-    let m_C', eigenvalues =
+    let m_C', _ =
       diagonalize_symm  m_F_diis
     in
 
+    (* Re-compute eigenvalues to remove level-shift *)
+    let eigenvalues = 
+      let m_F_ortho =
+        xt_o_x m_F0 m_X
+      in
+      xt_o_x m_F_ortho m_C'
+      |> Mat.copy_diag 
+    in
+
     (* MOs in AO basis *)
     let m_C =
       gemm  m_X  m_C'
@@ -344,6 +353,7 @@ let make
       |> amax
       |> abs_float
     in
+
     { empty with
       iteration      = nSCF              ;
       eigenvalues    = Some eigenvalues  ;
@@ -468,7 +478,7 @@ let make
       gemm  m_S  m_C
     in
 
-    let m_F = 
+    let m_F0 = 
       x_o_xt ~x:m_SC ~o:m_F_mo
     in
 
@@ -476,7 +486,7 @@ let make
     (* Add level shift in AO basis *)
     let m_F =
       x_o_xt ~x:m_SC ~o:m_LSmo
-      |> Mat.add m_F
+      |> Mat.add m_F0
     in
 
     (* Fock matrix in orthogonal basis *)
@@ -510,6 +520,15 @@ let make
       diagonalize_symm  m_F_diis
     in
 
+    (* Re-compute eigenvalues to remove level-shift *)
+    let eigenvalues = 
+      let m_F_ortho =
+        xt_o_x m_F0 m_X
+      in
+      xt_o_x m_F_ortho m_C'
+      |> Mat.copy_diag 
+    in
+
     (* MOs in AO basis *)
     let m_C =
       gemm  m_X  m_C'

run_mp2.ml

@@ -0,0 +1,81 @@
+open Lacaml.D 
+
+let () =
+  let open Command_line in
+  begin
+    set_header_doc (Sys.argv.(0) ^ " - QCaml command");
+    set_description_doc "Runs a Hartree-Fock calculation";
+    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. Default is 0"; } ;
+
+        { short='g' ; long="guess" ; opt=Optional;
+          arg=With_arg "<string>";
+          doc="Guess with another basis set."; } ;
+      ]
+  end;
+
+  let ppf = 
+    if Parallel.master then Format.std_formatter
+    else Printing.ppf_dev_null
+  in
+
+  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 -> int_of_string x
+    | None   -> 0
+  in
+
+  let multiplicity =
+    match Command_line.get "multiplicity" with
+    | Some x -> int_of_string x
+    | None -> 1
+  in
+
+  let s = Simulation.of_filenames
+      ~charge ~multiplicity ~nuclei:nuclei_file basis_file in
+
+  let guess = 
+    match Command_line.get "guess" with
+    | None -> `Huckel
+    | Some filename ->
+      let s_guess = Simulation.of_filenames
+          ~charge ~multiplicity ~nuclei:nuclei_file filename in
+      let hf = HartreeFock.make s_guess in
+      Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf;
+      let guess_mos =
+        MOBasis.of_hartree_fock hf
+        |> MOBasis.of_mo_basis s
+      in
+      `Matrix (MOBasis.mo_coef guess_mos)
+  in
+
+  let hf = HartreeFock.make ~guess s in
+
+  Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf;
+
+  let e_hf = HartreeFock.energy hf in
+
+  let mp2 = MP2.make hf in
+  Format.fprintf ppf "@[MP2   = %15.10f@]@." mp2;
+  Format.fprintf ppf "@[E+MP2 = %15.10f@]@." (mp2 +. e_hf)
+
+
+