From 2c8a303e404b7db8fe2c03d49f626e74ddde09a8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 18 Mar 2019 12:41:32 +0100 Subject: [PATCH] Added MP2 --- CI/CI.ml | 66 ++++++++++++++++++++++++++++++++++ CI/DeterminantSpace.ml | 19 ++++++++++ MOBasis/MOClass.ml | 43 ++++++++++++++++++++++ MOBasis/MOClass.mli | 12 +++++++ Makefile.include | 2 +- Perturbation/MP2.ml | 60 +++++++++++++++++++++++++++++++ SCF/HartreeFock.ml | 29 ++++++++++++--- run_mp2.ml | 81 ++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 306 insertions(+), 6 deletions(-) create mode 100644 Perturbation/MP2.ml create mode 100644 run_mp2.ml diff --git a/CI/CI.ml b/CI/CI.ml index 189c4d0..f77450a 100644 --- a/CI/CI.ml +++ b/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 + *) diff --git a/CI/DeterminantSpace.ml b/CI/DeterminantSpace.ml index bbc2586..d329404 100644 --- a/CI/DeterminantSpace.ml +++ b/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 diff --git a/MOBasis/MOClass.ml b/MOBasis/MOClass.ml index 7f23a98..682708a 100644 --- a/MOBasis/MOClass.ml +++ b/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) + ] + ) + diff --git a/MOBasis/MOClass.mli b/MOBasis/MOClass.mli index b43ef33..c1c43cf 100644 --- a/MOBasis/MOClass.mli +++ b/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} *) diff --git a/Makefile.include b/Makefile.include index 235218a..fec5e89 100644 --- a/Makefile.include +++ b/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) diff --git a/Perturbation/MP2.ml b/Perturbation/MP2.ml new file mode 100644 index 0000000..eb3ef7e --- /dev/null +++ b/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" + + diff --git a/SCF/HartreeFock.ml b/SCF/HartreeFock.ml index 273466b..1c8cfe7 100644 --- a/SCF/HartreeFock.ml +++ b/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' diff --git a/run_mp2.ml b/run_mp2.ml new file mode 100644 index 0000000..6be4bf8 --- /dev/null +++ b/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 ""; + 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. Default is 0"; } ; + + { short='g' ; long="guess" ; opt=Optional; + arg=With_arg ""; + 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) + + +