LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 04:13:33 +01:00
Code Releases Activity

Remove org-mode in mp2

Browse Source
This commit is contained in:
Anthony Scemama 2024-02-28 11:07:24 +01:00
parent ac44492ffe
commit 50698ce0ff
3 changed files with 35 additions and 228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
perturbation/lib/mp2.ml
View File

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

37
perturbation/lib/mp2.mli
View File

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

180
perturbation/mp2.org
View File

@ -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