diff --git a/perturbation/lib/mp2.ml b/perturbation/lib/mp2.ml index 6586cb6..380296f 100644 --- a/perturbation/lib/mp2.ml +++ b/perturbation/lib/mp2.ml @@ -1,24 +1,13 @@ -(* [[file:~/QCaml/perturbation/mp2.org::*Type][Type:2]] *) +(** Type *) + type t = { energy : float ; mo_basis : Mo.Basis.t ; frozen_core : Mo.Frozen_core.t ; } -(* Type:2 ends here *) +(** Creation *) - -(* | ~make~ | Creates an MP2 data structure | - * - * #+begin_example - * let mp2 = - * Perturbation.Mp2.make ~frozen_core:(Mo.Frozen_core.(make Small nuclei)) mo_basis - * ;; - * val mp2 : Perturbation.Mp2.t = E(MP2)=-0.185523 - * #+end_example *) - - -(* [[file:~/QCaml/perturbation/mp2.org::*Creation][Creation:2]] *) open Linear_algebra let make_rmp2 mo_basis mo_class = @@ -37,22 +26,22 @@ let make_rmp2 mo_basis mo_class = List.fold_left (fun accu b -> match b with Mo.Class.Virtual b -> let eps = -. (epsilon%.(b)) in - accu +. + accu +. List.fold_left (fun accu a -> match a with Mo.Class.Virtual a -> let eps = eps -. (epsilon%.(a)) in - accu +. + accu +. List.fold_left (fun accu j -> match j with Mo.Class.Inactive j -> let eps = eps +. epsilon%.(j) in - accu +. + 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 +. ijab *. ( abij +. abij -. abji) /. eps | _ -> accu ) 0. inactives | _ -> accu @@ -70,7 +59,7 @@ let make ~frozen_core mo_basis = |> Mo.Class.to_list in - let energy = + let energy = match Mo.Basis.mo_type mo_basis with | RHF -> make_rmp2 mo_basis mo_class | ROHF -> Common.Util.not_implemented "ROHF MP2" @@ -78,26 +67,17 @@ let make ~frozen_core mo_basis = | _ -> invalid_arg "MP2 needs RHF or ROHF MOs" in { energy ; mo_basis ; frozen_core } -(* Creation:2 ends here *) +(** Access *) -(* | ~energy~ | Returns the MP2 energy | - * | ~mo_basis~ | Returns the MO basis on which the MP2 energy was computed | - * | ~frozen_core~ | Returns the frozen_core scheme used to compute the MP2 energy | - * - * #+begin_example - * - * #+end_example *) - - -(* [[file:~/QCaml/perturbation/mp2.org::*Access][Access:2]] *) let energy t = t.energy let mo_basis t = t.mo_basis let frozen_core t = t.frozen_core -(* Access:2 ends here *) -(* [[file:~/QCaml/perturbation/mp2.org::*Printers][Printers:2]] *) + +(** Printers *) + let pp ppf t = Format.fprintf ppf "@[E(MP2)=%f@]" t.energy -(* Printers:2 ends here *) + diff --git a/perturbation/lib/mp2.mli b/perturbation/lib/mp2.mli index 5017a83..f9d294d 100644 --- a/perturbation/lib/mp2.mli +++ b/perturbation/lib/mp2.mli @@ -1,29 +1,36 @@ -(* Type *) +(** Data structure for an MP2 calculation *) + +(** + * let mp2 = + * Perturbation.Mp2.make ~frozen_core:(Mo.Frozen_core.(make Small nuclei)) mo_basis + * ;; + * val mp2 : Perturbation.Mp2.t = E(MP2)=-0.185523 + * + *) -(* [[file:~/QCaml/perturbation/mp2.org::*Type][Type:1]] *) +(** Type *) + type t -(* Type:1 ends here *) - -(* Creation *) -(* [[file:~/QCaml/perturbation/mp2.org::*Creation][Creation:1]] *) +(** Creation *) + val make : frozen_core:Mo.Frozen_core.t -> Mo.Basis.t -> t -(* Creation:1 ends here *) - -(* Access *) -(* [[file:~/QCaml/perturbation/mp2.org::*Access][Access:1]] *) +(** Access *) + val energy : t -> float +(* Returns the MP2 energy *) + val mo_basis : t -> Mo.Basis.t +(* Returns the MO basis on which the MP2 energy was computed *) + val frozen_core : t -> Mo.Frozen_core.t -(* Access:1 ends here *) - -(* Printers *) +(* Returns the frozen_core scheme used to compute the MP2 energy *) -(* [[file:~/QCaml/perturbation/mp2.org::*Printers][Printers:1]] *) +(** Printers *) + val pp : Format.formatter -> t -> unit -(* Printers:1 ends here *) diff --git a/perturbation/mp2.org b/perturbation/mp2.org deleted file mode 100644 index 773de2f..0000000 --- a/perturbation/mp2.org +++ /dev/null @@ -1,180 +0,0 @@ -#+begin_src elisp tangle: no :results none :exports none -(setq pwd (file-name-directory buffer-file-name)) -(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) -(setq lib (concat pwd "lib/")) -(setq testdir (concat pwd "test/")) -(setq mli (concat lib name ".mli")) -(setq ml (concat lib name ".ml")) -(setq test-ml (concat testdir name ".ml")) -(org-babel-tangle) -#+end_src - -* MP2 - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - -** Type - - #+begin_src ocaml :tangle (eval mli) -type t - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -type t = { - energy : float ; - mo_basis : Mo.Basis.t ; - frozen_core : Mo.Frozen_core.t ; -} - #+end_src - -** Creation - - #+begin_src ocaml :tangle (eval mli) -val make : frozen_core:Mo.Frozen_core.t -> Mo.Basis.t -> t - #+end_src - - | ~make~ | Creates an MP2 data structure | - - #+begin_example -let mp2 = - Perturbation.Mp2.make ~frozen_core:(Mo.Frozen_core.(make Small nuclei)) mo_basis - ;; -val mp2 : Perturbation.Mp2.t = E(MP2)=-0.185523 - #+end_example - - #+begin_src ocaml :tangle (eval ml) :exports none -open Linear_algebra - -let make_rmp2 mo_basis mo_class = - - let epsilon = Mo.Basis.mo_energies mo_basis 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 - - 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 - - -let make ~frozen_core mo_basis = - - let mo_class = - Mo.Class.cas_sd mo_basis ~frozen_core 0 0 - |> Mo.Class.to_list - in - - let energy = - match Mo.Basis.mo_type mo_basis with - | RHF -> make_rmp2 mo_basis mo_class - | ROHF -> Common.Util.not_implemented "ROHF MP2" - | UHF -> Common.Util.not_implemented "UHF MP2" - | _ -> invalid_arg "MP2 needs RHF or ROHF MOs" - in - { energy ; mo_basis ; frozen_core } - #+end_src - -** Access - - #+begin_src ocaml :tangle (eval mli) -val energy : t -> float -val mo_basis : t -> Mo.Basis.t -val frozen_core : t -> Mo.Frozen_core.t - #+end_src - - | ~energy~ | Returns the MP2 energy | - | ~mo_basis~ | Returns the MO basis on which the MP2 energy was computed | - | ~frozen_core~ | Returns the frozen_core scheme used to compute the MP2 energy | - - #+begin_example - - #+end_example - - #+begin_src ocaml :tangle (eval ml) :exports none -let energy t = t.energy -let mo_basis t = t.mo_basis -let frozen_core t = t.frozen_core - #+end_src - -** Printers - - #+begin_src ocaml :tangle (eval mli) -val pp : Format.formatter -> t -> unit - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -let pp ppf t = - Format.fprintf ppf "@[E(MP2)=%f@]" t.energy - #+end_src - -** Tests - - #+begin_src ocaml :tangle (eval test-ml) :exports none -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 2.e-10) "Energy" (-76.0267987005) (Mo.Hartree_fock.energy hf); - - let mo_basis = Mo.Basis.of_hartree_fock hf in - - let frozen_core = Mo.Frozen_core.(make Small nuclei) in - - let mp2 = Perturbation.Mp2.make ~frozen_core mo_basis in - - let e_mp2 = Perturbation.Mp2.energy mp2 in - - check (float 1.e-9) "MP2" (-0.2016211415) (e_mp2) - ] - #+end_src -