LCPQ
/
QCaml
mirror of https://gitlab.com/scemama/QCaml.git
10
1
Fork 0
Code Releases Activity

Compare commits

base: LCPQ:c4b022aeeecf4063f7daabff02f17fca7442e88e
Branches Tags
LCPQ:master
LCPQ:dev
LCPQ:smallarray
...
compare: LCPQ:d6ef25d55ac5254102066ff725cf94c915ea6c38
Branches Tags
LCPQ:dev
LCPQ:smallarray
LCPQ:master

3 Commits
c4b022aeee ... d6ef25d55a

Author SHA1 Message Date
Anthony Scemama d6ef25d55a Removing org-mode 2024-01-26 11:24:56 +01:00
Anthony Scemama 744f2a0552 Removing org-mode 2024-01-17 14:24:28 +01:00
Anthony Scemama e59d016ae7 Removing org-mode 2024-01-17 13:59:05 +01:00
31 changed files with 550 additions and 1740 deletions
Show Stats Expand all files Collapse all files

49
common/lib/powers.ml
View File

@ -1,30 +1,15 @@
(** Type *)
(* ~tot~ always contains ~x+y+z~. *)
(* [[file:~/QCaml/common/powers.org::*Type][Type:2]] *)
type t = {
x : int ;
y : int ;
z : int ;
tot : int ;
tot : int ; (* ~tot~ always contains ~x+y+z~. *)
}
(* Type:2 ends here *)
(** Conversions *)
(* Example:
* #+begin_example
* Powers.of_int_tuple (2,3,1);;
* - : Powers.t = x^2 + y^3 + z^1
*
* Powers.(to_int_tuple (of_int_tuple (2,3,1)));;
* - : int * int * int = (2, 3, 1)
* #+end_example *)
(* [[file:~/QCaml/common/powers.org::*Conversions][Conversions:2]] *)
let of_int_tuple t =
let result =
match t with
@ -39,29 +24,10 @@ let of_int_tuple t =
let to_int_tuple { x ; y ; z ; _ } = (x,y,z)
(* Conversions:2 ends here *)
(** Operations *)
(* | ~get~ | Returns the value of the power for $x$, $y$ or $z$
* | ~incr~ | Returns a new ~Powers.t~ with the power on the given axis incremented |
* | ~decr~ | Returns a new ~Powers.t~ with the power on the given axis decremented. As opposed to ~of_int_tuple~, the values may become negative|
*
* Example:
* #+begin_example
* Powers.get Coordinate.Y (Powers.of_int_tuple (2,3,1));;
* - : int = 3
*
* Powers.incr Coordinate.Y (Powers.of_int_tuple (2,3,1));;
* - : Powers.t = x^2 + y^4 + z^1
*
* Powers.decr Coordinate.Y (Powers.of_int_tuple (2,3,1));;
* - : Powers.t = x^2 + y^2 + z^1
* #+end_example *)
(* [[file:~/QCaml/common/powers.org::*Operations][Operations:2]] *)
let get coord t =
match coord with
| Coordinate.X -> t.x
@ -79,9 +45,10 @@ let decr coord t =
| Coordinate.X -> let r = t.x-1 in { t with x = r ; tot = t.tot-1 }
| Coordinate.Y -> let r = t.y-1 in { t with y = r ; tot = t.tot-1 }
| Coordinate.Z -> let r = t.z-1 in { t with z = r ; tot = t.tot-1 }
(* Operations:2 ends here *)
(* [[file:~/QCaml/common/powers.org::*Printers][Printers:2]] *)
(** Printers *)
let pp ppf t =
Format.fprintf ppf "@[x^%d + y^%d + z^%d@]" t.x t.y t.z
(* Printers:2 ends here *)

56
common/lib/powers.mli
View File

@ -1,36 +1,58 @@
(* Type
*
* <<<~Powers.t~>>> *)
(** Powers *)
(* Contains powers of $x$, $y$ and $z$ describing the polynomials in
* atomic basis sets.
*)
(** Type *)
(* [[file:~/QCaml/common/powers.org::*Type][Type:1]] *)
type t = private {
x : int ;
y : int ;
z : int ;
tot : int ;
tot : int ; (* ~tot~ always contains ~x+y+z~. *)
}
(* Type:1 ends here *)
(* Conversions *)
(** Conversions *)
(* [[file:~/QCaml/common/powers.org::*Conversions][Conversions:1]] *)
val of_int_tuple : int * int * int -> t
(** Powers.of_int_tuple (2,3,1);;
* - : Powers.t = x^2 + y^3 + z^1
*)
val to_int_tuple : t -> int * int * int
(* Conversions:1 ends here *)
(* Operations *)
(** Powers.(to_int_tuple (of_int_tuple (2,3,1)));;
* - : int * int * int = (2, 3, 1)
*)
(* [[file:~/QCaml/common/powers.org::*Operations][Operations:1]] *)
(** Operations *)
val get : Coordinate.axis -> t -> int
(** Returns the value of the power for $x$, $y$ or $z$.
*
* Powers.get Coordinate.Y (Powers.of_int_tuple (2,3,1));;
* - : int = 3
*
*)
val incr : Coordinate.axis -> t -> t
(** Returns a new ~Powers.t~ with the power on the given axis incremented.
*
* Powers.incr Coordinate.Y (Powers.of_int_tuple (2,3,1));;
* - : Powers.t = x^2 + y^4 + z^1
*)
val decr : Coordinate.axis -> t -> t
(* Operations:1 ends here *)
(* Printers *)
(** Returns a new ~Powers.t~ with the power on the given axis decremented.
* As opposed to ~of_int_tuple~, the values may become negative.
*
* Powers.decr Coordinate.Y (Powers.of_int_tuple (2,3,1));;
* - : Powers.t = x^2 + y^2 + z^1
*)
(* [[file:~/QCaml/common/powers.org::*Printers][Printers:1]] *)
(** Printers *)
val pp : Format.formatter -> t -> unit
(* Printers:1 ends here *)

16
common/lib/range.ml
View File

@ -1,8 +1,10 @@
(* [[file:~/QCaml/common/range.org::*Type][Type:2]] *)
type t = int list
(* Type:2 ends here *)
(** Type *)
type t = int list
(** Conversion *)
(* [[file:~/QCaml/common/range.org::*Conversion][Conversion:2]] *)
let to_int_list r = r
let expand_range r =
@ -41,9 +43,9 @@ let to_string l =
(List.map string_of_int l
|> String.concat ",") ^
"]"
(* Conversion:2 ends here *)
(* [[file:~/QCaml/common/range.org::*Printers][Printers:2]] *)
(** Printers *)
let pp ppf t =
Format.fprintf ppf "@[%s@]" (to_string t)
(* Printers:2 ends here *)

27
common/lib/range.mli
View File

@ -1,23 +1,26 @@
(* Type
(* Range *)
(* A range is a sorted list of integers in an interval.
*
* <<<~Range.t~>>> *)
* - ~"[a-b]"~ : range between a and b (included)
* - ~"[a]"~ : the list with only one integer a
* - ~"a"~ : equivalent to "[a]"
* - ~"[36-53,72-107,126-131]"~ represents the list of integers
* [ 37 ; 37 ; 38 ; ... ; 52 ; 53 ; 72 ; 73 ; ... ; 106 ; 107 ; 126 ; 127 ; ... ; 130 ; 131 ].
*)
(** Type *)
(* [[file:~/QCaml/common/range.org::*Type][Type:1]] *)
type t
(* Type:1 ends here *)
(* Conversion *)
(* [[file:~/QCaml/common/range.org::*Conversion][Conversion:1]] *)
(** Conversion *)
val of_string : string -> t
val to_string : t -> string
val to_int_list : t -> int list
(* Conversion:1 ends here *)
(* Printers *)
(* [[file:~/QCaml/common/range.org::*Printers][Printers:1]] *)
(** Printers *)
val pp : Format.formatter -> t -> unit
(* Printers:1 ends here *)

20
common/lib/spin.ml
View File

@ -1,22 +1,14 @@
(** Spin *)
(** Type *)
(* Note :
* ~Alfa~ if written with an 'f' instead of 'ph' because it has the same number of
* letters as ~Beta~, so the alignment of the code is nicer. *)
(* [[file:~/QCaml/common/spin.org::*Type][Type:2]] *)
type t = (* m_s *)
| Alfa (* {% $m_s = +1/2$ %} *)
| Beta (* {% $m_s = -1/2$ %} *)
(* Type:2 ends here *)
(** Functions *)
(* Returns the opposite spin *)
(* [[file:~/QCaml/common/spin.org::*Functions][Functions:2]] *)
let other = function
| Alfa -> Beta
| Beta -> Alfa
@ -24,9 +16,9 @@ let other = function
let to_string = function
| Alfa -> "Alpha"
| Beta -> "Beta "
(* Functions:2 ends here *)
(* [[file:~/QCaml/common/spin.org::*Printers][Printers:2]] *)
(** Printers *)
let pp ppf t =
Format.fprintf ppf "@[%s@]" (to_string t)
(* Printers:2 ends here *)

26
common/lib/spin.mli
View File

@ -1,21 +1,23 @@
(* Type
*
* <<<~Spin.t~>>> *)
(** Spin *)
(* [[file:~/QCaml/common/spin.org::*Type][Type:1]] *)
(* Electron spin *)
(** Type *)
type t = Alfa | Beta
(* Type:1 ends here *)
(* Functions *)
(** Note :
~Alfa~ if written with an 'f' instead of 'ph' because it has the same
number of letters as ~Beta~, so the alignment of the code is nicer.
*)
(* [[file:~/QCaml/common/spin.org::*Functions][Functions:1]] *)
(** Functions *)
val other : t -> t
(* Functions:1 ends here *)
(* Printers *)
(** Returns the opposite spin. *)
(* [[file:~/QCaml/common/spin.org::*Printers][Printers:1]] *)
(** Printers *)
val pp : Format.formatter -> t -> unit
(* Printers:1 ends here *)

100
common/lib/zkey.ml
View File

@ -1,5 +1,36 @@
(* [[file:~/QCaml/common/zkey.org::*Types][Types:2]] *)
type t =
(** ZKey *)
(*
* Encodes the powers of x, y, z in a compact form, suitable for being
* used as keys in a hash table.
*
* Internally, the ~Zkey.t~ is made of two integers, ~left~ and ~right~.
* The small integers x, y and z are stored compactly in this 126-bits
* space:
*
* Example:
* Left Right
* 3 [--------------------------------------------------------------] [------------------|---------------|---------------|---------------]
* x y z
*
* 6 [--------------------------------------------------------------] [---|----------|----------|----------|----------|----------|---------]
* x1 y1 z1 x2 y2 z2
*
* 9 [---------------------------------|----------|----------|---------] [---|----------|----------|----------|----------|----------|---------]
* x1 y1 z1 x2 y2 z2 x3 y3 z3
*
* 12 [---|----------|----------|----------|----------|----------|---------] [---|----------|----------|----------|----------|----------|---------]
* x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4
*
*
* The values of x,y,z should be positive and should not exceed 32767 for
* ~kind=3~. For all other kinds kinds the values should not exceed 1023.
*
*)
(** Types *)
type t =
{
mutable left : int;
mutable right : int;
@ -8,46 +39,31 @@ type t =
open Powers
type kind =
type kind =
| Three of Powers.t
| Four of (int * int * int * int)
| Six of (Powers.t * Powers.t)
| Nine of (Powers.t * Powers.t * Powers.t)
| Twelve of (Powers.t * Powers.t * Powers.t * Powers.t)
(* Types:2 ends here *)
(* | ~of_powers_three~ | Create from a ~Powers.t~ |
* | ~of_powers_six~ | Create from two ~Powers.t~ |
* | ~of_powers_nine~ | Create from three ~Powers.t~ |
* | ~of_powers_twelve~ | Create from four ~Powers.t~ |
* | ~of_powers~ | Create using the ~kind~ type |
* | ~of_int_array~ | Convert from an ~int~ array |
* | ~of_int_four~ | Create from four ~ints~ |
* | ~to_int_array~ | Convert to an ~int~ array |
* | ~to_powers~ | Convert to an ~Powers.t~ array |
* | ~to_string~ | Pretty printing | *)
(* [[file:~/QCaml/common/zkey.org::*Conversions][Conversions:2]] *)
(** Creates a Zkey. *)
let make ~kind right =
{ left = 0 ; right ; kind }
(** Move [right] to [left] and set [right = x] *)
let (<|) z x =
let (<|) z x =
z.left <- z.right;
z.right <- x;
z
(** Shift left [right] by 10 bits, and add [x]. *)
let (<<) z x =
let (<<) z x =
z.right <- (z.right lsl 10) lor x ;
z
(** Shift left [right] by 10 bits, and add [x]. *)
let (<+) z x =
let (<+) z x =
z.right <- (z.right lsl 15) lor x ;
z
@ -71,34 +87,34 @@ let of_int_four i j k l =
let of_powers_six { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } =
assert (
let alpha = a lor b lor c lor d lor e lor f in
alpha >= 0 && alpha < (1 lsl 10)
alpha >= 0 && alpha < (1 lsl 10)
);
make ~kind:6 a << b << c << d << e << f
make ~kind:6 a << b << c << d << e << f
let of_powers_nine { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ }
{ x=g ; y=h ; z=i ; _ } =
{ x=g ; y=h ; z=i ; _ } =
assert (
let alpha = a lor b lor c lor d lor e lor f lor g lor h lor i in
alpha >= 0 && alpha < (1 lsl 10)
alpha >= 0 && alpha < (1 lsl 10)
);
make ~kind:9 a << b << c << d << e << f
<| g << h << i
<| g << h << i
let of_powers_twelve { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ }
{ x=g ; y=h ; z=i ; _ } { x=j ; y=k ; z=l ; _ } =
{ x=g ; y=h ; z=i ; _ } { x=j ; y=k ; z=l ; _ } =
assert (
let alpha = a lor b lor c lor d lor e lor f
lor g lor h lor i lor j lor k lor l
in
alpha >= 0 && alpha < (1 lsl 10)
alpha >= 0 && alpha < (1 lsl 10)
);
make ~kind:12 a << b << c << d << e << f
<| g << h << i << j << k << l
let of_powers a =
let of_powers a =
match a with
| Three a -> of_powers_three a
| Six (a,b) -> of_powers_six a b
@ -117,7 +133,7 @@ let of_int_array = function
(** Transform the Zkey into an int array *)
let to_int_array { left ; right ; kind } =
let to_int_array { left ; right ; kind } =
match kind with
| 3 -> [|
mask15 land (right lsr 30) ;
@ -132,7 +148,7 @@ let to_int_array { left ; right ; kind } =
mask15 land right
|]
| 6 -> [|
| 6 -> [|
mask10 land (right lsr 50) ;
mask10 land (right lsr 40) ;
mask10 land (right lsr 30) ;
@ -172,7 +188,7 @@ let to_int_array { left ; right ; kind } =
(** Transform the Zkey into an int tuple *)
let to_powers { left ; right ; kind } =
let to_powers { left ; right ; kind } =
match kind with
| 3 -> Three (Powers.of_int_tuple (
mask15 land (right lsr 30) ,
@ -180,7 +196,7 @@ let to_powers { left ; right ; kind } =
mask15 land right
))
| 6 -> Six (Powers.of_int_tuple
| 6 -> Six (Powers.of_int_tuple
( mask10 land (right lsr 50) ,
mask10 land (right lsr 40) ,
mask10 land (right lsr 30)),
@ -190,7 +206,7 @@ let to_powers { left ; right ; kind } =
mask10 land right )
)
| 12 -> Twelve (Powers.of_int_tuple
| 12 -> Twelve (Powers.of_int_tuple
( mask10 land (left lsr 50) ,
mask10 land (left lsr 40) ,
mask10 land (left lsr 30)),
@ -208,7 +224,7 @@ let to_powers { left ; right ; kind } =
mask10 land right )
)
| 9 -> Nine (Powers.of_int_tuple
| 9 -> Nine (Powers.of_int_tuple
( mask10 land (left lsr 20) ,
mask10 land (left lsr 10) ,
mask10 land left ) ,
@ -223,17 +239,9 @@ let to_powers { left ; right ; kind } =
)
| _ -> invalid_arg (__FILE__^": to_powers")
(* Conversions:2 ends here *)
(* | ~hash~ | Associates a nonnegative integer to any Zkey |
* | ~equal~ | The equal function. True if two Zkeys are equal |
* | ~compare~ | Comparison function, used for sorting | *)
(* [[file:~/QCaml/common/zkey.org::*Functions for hash tables][Functions for hash tables:2]] *)
let hash = Hashtbl.hash
let hash = Hashtbl.hash
let equal
{ right = r1 ; left = l1 ; kind = k1 }
@ -259,9 +267,7 @@ let to_string { left ; right ; kind } =
|> Array.to_list
|> String.concat ", "
) ^ " >"
(* Functions for hash tables:2 ends here *)
(* [[file:~/QCaml/common/zkey.org::*Printers][Printers:2]] *)
(** Printers *)
let pp ppf t =
Format.fprintf ppf "@[%s@]" (to_string t)
(* Printers:2 ends here *)

51
common/lib/zkey.mli
View File

@ -1,9 +1,6 @@
(* Types
*
* <<<~Zkey.t~>>> *)
(** Types *)
(* [[file:~/QCaml/common/zkey.org::*Types][Types:1]] *)
type t
type t
type kind =
| Three of Powers.t
@ -11,36 +8,54 @@ type kind =
| Six of (Powers.t * Powers.t)
| Nine of (Powers.t * Powers.t * Powers.t)
| Twelve of (Powers.t * Powers.t * Powers.t * Powers.t)
(* Types:1 ends here *)
(* Conversions *)
(* [[file:~/QCaml/common/zkey.org::*Conversions][Conversions:1]] *)
(** Conversions *)
val of_powers_three : Powers.t -> t
(** Create from a ~Powers.t~ *)
val of_powers_six : Powers.t -> Powers.t -> t
(** Create from two ~Powers.t~ *)
val of_powers_nine : Powers.t -> Powers.t -> Powers.t -> t
(** Create from three ~Powers.t~ *)
val of_powers_twelve : Powers.t -> Powers.t -> Powers.t -> Powers.t -> t
(** Create from four ~Powers.t~ *)
val of_powers : kind -> t
(** Create using the ~kind~ type *)
val of_int_array : int array -> t
(** Convert from an ~int~ array *)
val of_int_four : int -> int -> int -> int -> t
(** Create from four ~ints~ *)
val to_int_array : t -> int array
(** Convert to an ~int~ array *)
val to_powers : t -> kind
(** Convert to an ~Powers.t~ array *)
val to_string : t -> string
(* Conversions:1 ends here *)
(* Functions for hash tables *)
(** Pretty printing *)
(* [[file:~/QCaml/common/zkey.org::*Functions for hash tables][Functions for hash tables:1]] *)
(** Functions for hash tables *)
val hash : t -> int
(** Associates a nonnegative integer to any Zkey *)
val equal : t -> t -> bool
(** The equal function. True if two Zkeys are equal *)
val compare : t -> t -> int
(* Functions for hash tables:1 ends here *)
(* Printers *)
(** Comparison function, used for sorting *)
(* [[file:~/QCaml/common/zkey.org::*Printers][Printers:1]] *)
(** Printers *)
val pp : Format.formatter -> t -> unit
(* Printers:1 ends here *)

2
common/lib/zmap.ml
View File

@ -1,4 +1,2 @@
(* [[file:~/QCaml/common/zmap.org::*Type][Type:2]] *)
module Zmap = Hashtbl.Make(Zkey)
include Zmap
(* Type:2 ends here *)

9
common/lib/zmap.mli
View File

@ -1,7 +1,6 @@
(* Type
*
* <<<~Zmap.t~>>> *)
(** ZMap *)
(* [[file:~/QCaml/common/zmap.org::*Type][Type:1]] *)
(* A hash table where the keys are ~Zkey~ *)
(** Type *)
include module type of Hashtbl.Make(Zkey)
(* Type:1 ends here *)

131
common/powers.org
View File

@ -1,131 +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
* Powers
:PROPERTIES:
:header-args: :noweb yes :comments both
:END:
Contains powers of $x$, $y$ and $z$ describing the polynomials in atomic basis sets.
** Type
<<<~Powers.t~>>>
#+begin_src ocaml :tangle (eval mli)
type t = private {
x : int ;
y : int ;
z : int ;
tot : int ;
}
#+end_src
~tot~ always contains ~x+y+z~.
#+begin_src ocaml :tangle (eval ml) :exports none
type t = {
x : int ;
y : int ;
z : int ;
tot : int ;
}
#+end_src
** Conversions
#+begin_src ocaml :tangle (eval mli)
val of_int_tuple : int * int * int -> t
val to_int_tuple : t -> int * int * int
#+end_src
Example:
#+begin_example
Powers.of_int_tuple (2,3,1);;
- : Powers.t = x^2 + y^3 + z^1
Powers.(to_int_tuple (of_int_tuple (2,3,1)));;
- : int * int * int = (2, 3, 1)
#+end_example
#+begin_src ocaml :tangle (eval ml) :exports none
let of_int_tuple t =
let result =
match t with
| (x,y,z) -> { x ; y ; z ; tot=x+y+z }
in
if result.x < 0 ||
result.y < 0 ||
result.z < 0 ||
result.tot < 0 then
invalid_arg (__FILE__^": of_int_tuple");
result
let to_int_tuple { x ; y ; z ; _ } = (x,y,z)
#+end_src
** Operations
#+begin_src ocaml :tangle (eval mli)
val get : Coordinate.axis -> t -> int
val incr : Coordinate.axis -> t -> t
val decr : Coordinate.axis -> t -> t
#+end_src
| ~get~ | Returns the value of the power for $x$, $y$ or $z$
| ~incr~ | Returns a new ~Powers.t~ with the power on the given axis incremented |
| ~decr~ | Returns a new ~Powers.t~ with the power on the given axis decremented. As opposed to ~of_int_tuple~, the values may become negative|
Example:
#+begin_example
Powers.get Coordinate.Y (Powers.of_int_tuple (2,3,1));;
- : int = 3
Powers.incr Coordinate.Y (Powers.of_int_tuple (2,3,1));;
- : Powers.t = x^2 + y^4 + z^1
Powers.decr Coordinate.Y (Powers.of_int_tuple (2,3,1));;
- : Powers.t = x^2 + y^2 + z^1
#+end_example
#+begin_src ocaml :tangle (eval ml) :exports none
let get coord t =
match coord with
| Coordinate.X -> t.x
| Coordinate.Y -> t.y
| Coordinate.Z -> t.z
let incr coord t =
match coord with
| Coordinate.X -> let r = t.x+1 in { t with x = r ; tot = t.tot+1 }
| Coordinate.Y -> let r = t.y+1 in { t with y = r ; tot = t.tot+1 }
| Coordinate.Z -> let r = t.z+1 in { t with z = r ; tot = t.tot+1 }
let decr coord t =
match coord with
| Coordinate.X -> let r = t.x-1 in { t with x = r ; tot = t.tot-1 }
| Coordinate.Y -> let r = t.y-1 in { t with y = r ; tot = t.tot-1 }
| Coordinate.Z -> let r = t.z-1 in { t with z = r ; tot = t.tot-1 }
#+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 "@[x^%d + y^%d + z^%d@]" t.x t.y t.z
#+end_src

97
common/range.org
View File

@ -1,97 +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
* Range
:PROPERTIES:
:header-args: :noweb yes :comments both
:END:
A range is a sorted list of integers in an interval.
- ~"[a-b]"~ : range between a and b (included)
- ~"[a]"~ : the list with only one integer a
- ~"a"~ : equivalent to "[a]"
- ~"[36-53,72-107,126-131]"~ represents the list of integers
[ 37 ; 37 ; 38 ; ... ; 52 ; 53 ; 72 ; 73 ; ... ; 106 ; 107 ; 126 ; 127 ; ... ; 130 ; 131 ].
** Type
<<<~Range.t~>>>
#+begin_src ocaml :tangle (eval mli)
type t
#+end_src
#+begin_src ocaml :tangle (eval ml) :exports none
type t = int list
#+end_src
** Conversion
#+begin_src ocaml :tangle (eval mli)
val of_string : string -> t
val to_string : t -> string
val to_int_list : t -> int list
#+end_src
#+begin_src ocaml :tangle (eval ml) :exports none
let to_int_list r = r
let expand_range r =
match String.split_on_char '-' r with
| s :: f :: [] ->
begin
let start = int_of_string s
and finish = int_of_string f
in
assert (start <= finish) ;
let rec do_work = function
| i when i=finish -> [ i ]
| i -> i::(do_work (i+1))
in do_work start
end
| r :: [] -> [int_of_string r]
| [] -> []
| _ -> invalid_arg "Only one range expected"
let of_string s =
match s.[0] with
| '0' .. '9' -> [ int_of_string s ]
| _ ->
assert (s.[0] = '[') ;
assert (s.[(String.length s)-1] = ']') ;
let s = String.sub s 1 ((String.length s) - 2) in
let l = String.split_on_char ',' s in
let l = List.map expand_range l in
List.concat l
|> List.sort_uniq compare
let to_string l =
"[" ^
(List.map string_of_int l
|> String.concat ",") ^
"]"
#+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 "@[%s@]" (to_string t)
#+end_src

64
common/spin.org
View File

@ -1,64 +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
* Spin
:PROPERTIES:
:header-args: :noweb yes :comments both
:END:
Electron spin
** Type
<<<~Spin.t~>>>
#+begin_src ocaml :tangle (eval mli)
type t = Alfa | Beta
#+end_src
Note :
~Alfa~ if written with an 'f' instead of 'ph' because it has the same number of
letters as ~Beta~, so the alignment of the code is nicer.
#+begin_src ocaml :tangle (eval ml) :exports none
type t = (* m_s *)
| Alfa (* {% $m_s = +1/2$ %} *)
| Beta (* {% $m_s = -1/2$ %} *)
#+end_src
** Functions
#+begin_src ocaml :tangle (eval mli)
val other : t -> t
#+end_src
Returns the opposite spin
#+begin_src ocaml :tangle (eval ml) :exports none
let other = function
| Alfa -> Beta
| Beta -> Alfa
let to_string = function
| Alfa -> "Alpha"
| Beta -> "Beta "
#+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 "@[%s@]" (to_string t)
#+end_src

348
common/zkey.org
View File

@ -1,348 +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
* Zkey
:PROPERTIES:
:header-args: :noweb yes :comments both
:END:
Encodes the powers of x, y, z in a compact form, suitable for being
used as keys in a hash table.
Internally, the ~Zkey.t~ is made of two integers, ~left~ and ~right~.
The small integers x, y and z are stored compactly in this 126-bits
space:
Example:
#+begin_example
Left Right
3 [--------------------------------------------------------------] [------------------|---------------|---------------|---------------]
x y z
6 [--------------------------------------------------------------] [---|----------|----------|----------|----------|----------|---------]
x1 y1 z1 x2 y2 z2
9 [---------------------------------|----------|----------|---------] [---|----------|----------|----------|----------|----------|---------]
x1 y1 z1 x2 y2 z2 x3 y3 z3
12 [---|----------|----------|----------|----------|----------|---------] [---|----------|----------|----------|----------|----------|---------]
x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4
#+end_example
The values of x,y,z should be positive and should not exceed 32767 for
~kind=3~. For all other kinds kinds the values should not exceed 1023.
** Types
<<<~Zkey.t~>>>
#+begin_src ocaml :tangle (eval mli)
type t
type kind =
| Three of Powers.t
| Four of (int * int * int * int)
| Six of (Powers.t * Powers.t)
| Nine of (Powers.t * Powers.t * Powers.t)
| Twelve of (Powers.t * Powers.t * Powers.t * Powers.t)
#+end_src
#+begin_src ocaml :tangle (eval ml) :exports none
type t =
{
mutable left : int;
mutable right : int;
kind : int ;
}
open Powers
type kind =
| Three of Powers.t
| Four of (int * int * int * int)
| Six of (Powers.t * Powers.t)
| Nine of (Powers.t * Powers.t * Powers.t)
| Twelve of (Powers.t * Powers.t * Powers.t * Powers.t)
#+end_src
** Conversions
#+begin_src ocaml :tangle (eval mli)
val of_powers_three : Powers.t -> t
val of_powers_six : Powers.t -> Powers.t -> t
val of_powers_nine : Powers.t -> Powers.t -> Powers.t -> t
val of_powers_twelve : Powers.t -> Powers.t -> Powers.t -> Powers.t -> t
val of_powers : kind -> t
val of_int_array : int array -> t
val of_int_four : int -> int -> int -> int -> t
val to_int_array : t -> int array
val to_powers : t -> kind
val to_string : t -> string
#+end_src
| ~of_powers_three~ | Create from a ~Powers.t~ |
| ~of_powers_six~ | Create from two ~Powers.t~ |
| ~of_powers_nine~ | Create from three ~Powers.t~ |
| ~of_powers_twelve~ | Create from four ~Powers.t~ |
| ~of_powers~ | Create using the ~kind~ type |
| ~of_int_array~ | Convert from an ~int~ array |
| ~of_int_four~ | Create from four ~ints~ |
| ~to_int_array~ | Convert to an ~int~ array |
| ~to_powers~ | Convert to an ~Powers.t~ array |
| ~to_string~ | Pretty printing |
#+begin_src ocaml :tangle (eval ml) :exports none
(** Creates a Zkey. *)
let make ~kind right =
{ left = 0 ; right ; kind }
(** Move [right] to [left] and set [right = x] *)
let (<|) z x =
z.left <- z.right;
z.right <- x;
z
(** Shift left [right] by 10 bits, and add [x]. *)
let (<<) z x =
z.right <- (z.right lsl 10) lor x ;
z
(** Shift left [right] by 10 bits, and add [x]. *)
let (<+) z x =
z.right <- (z.right lsl 15) lor x ;
z
let of_powers_three { x=a ; y=b ; z=c ; _ } =
assert (
let alpha = a lor b lor c in
alpha >= 0 && alpha < (1 lsl 15)
);
make ~kind:3 a <+ b <+ c
let of_int_four i j k l =
assert (
let alpha = i lor j lor k lor l in
alpha >= 0 && alpha < (1 lsl 15)
);
make ~kind:4 i <+ j <+ k <+ l
let of_powers_six { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } =
assert (
let alpha = a lor b lor c lor d lor e lor f in
alpha >= 0 && alpha < (1 lsl 10)
);
make ~kind:6 a << b << c << d << e << f
let of_powers_nine { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ }
{ x=g ; y=h ; z=i ; _ } =
assert (
let alpha = a lor b lor c lor d lor e lor f lor g lor h lor i in
alpha >= 0 && alpha < (1 lsl 10)
);
make ~kind:9 a << b << c << d << e << f
<| g << h << i
let of_powers_twelve { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ }
{ x=g ; y=h ; z=i ; _ } { x=j ; y=k ; z=l ; _ } =
assert (
let alpha = a lor b lor c lor d lor e lor f
lor g lor h lor i lor j lor k lor l
in
alpha >= 0 && alpha < (1 lsl 10)
);
make ~kind:12 a << b << c << d << e << f
<| g << h << i << j << k << l
let of_powers a =
match a with
| Three a -> of_powers_three a
| Six (a,b) -> of_powers_six a b
| Twelve (a,b,c,d) -> of_powers_twelve a b c d
| Nine (a,b,c) -> of_powers_nine a b c
| _ -> invalid_arg "of_powers"
let mask10 = 0x3ff
and mask15 = 0x7fff
let of_int_array = function
| [| a ; b ; c ; d |] -> of_int_four a b c d
| _ -> invalid_arg "of_int_array"
(** Transform the Zkey into an int array *)
let to_int_array { left ; right ; kind } =
match kind with
| 3 -> [|
mask15 land (right lsr 30) ;
mask15 land (right lsr 15) ;
mask15 land right
|]
| 4 -> [|
mask15 land (right lsr 45) ;
mask15 land (right lsr 30) ;
mask15 land (right lsr 15) ;
mask15 land right
|]
| 6 -> [|
mask10 land (right lsr 50) ;
mask10 land (right lsr 40) ;
mask10 land (right lsr 30) ;
mask10 land (right lsr 20) ;
mask10 land (right lsr 10) ;
mask10 land right
|]
| 12 -> [|
mask10 land (left lsr 50) ;
mask10 land (left lsr 40) ;
mask10 land (left lsr 30) ;
mask10 land (left lsr 20) ;
mask10 land (left lsr 10) ;
mask10 land left ;
mask10 land (right lsr 50) ;
mask10 land (right lsr 40) ;
mask10 land (right lsr 30) ;
mask10 land (right lsr 20) ;
mask10 land (right lsr 10) ;
mask10 land right
|]
| 9 -> [|
mask10 land (left lsr 20) ;
mask10 land (left lsr 10) ;
mask10 land left ;
mask10 land (right lsr 50) ;
mask10 land (right lsr 40) ;
mask10 land (right lsr 30) ;
mask10 land (right lsr 20) ;
mask10 land (right lsr 10) ;
mask10 land right
|]
| _ -> invalid_arg (__FILE__^": to_int_array")
(** Transform the Zkey into an int tuple *)
let to_powers { left ; right ; kind } =
match kind with
| 3 -> Three (Powers.of_int_tuple (
mask15 land (right lsr 30) ,
mask15 land (right lsr 15) ,
mask15 land right
))
| 6 -> Six (Powers.of_int_tuple
( mask10 land (right lsr 50) ,
mask10 land (right lsr 40) ,
mask10 land (right lsr 30)),
Powers.of_int_tuple
( mask10 land (right lsr 20) ,
mask10 land (right lsr 10) ,
mask10 land right )
)
| 12 -> Twelve (Powers.of_int_tuple
( mask10 land (left lsr 50) ,
mask10 land (left lsr 40) ,
mask10 land (left lsr 30)),
Powers.of_int_tuple
( mask10 land (left lsr 20) ,
mask10 land (left lsr 10) ,
mask10 land left ) ,
Powers.of_int_tuple
( mask10 land (right lsr 50) ,
mask10 land (right lsr 40) ,
mask10 land (right lsr 30)),
Powers.of_int_tuple
( mask10 land (right lsr 20) ,
mask10 land (right lsr 10) ,
mask10 land right )
)
| 9 -> Nine (Powers.of_int_tuple
( mask10 land (left lsr 20) ,
mask10 land (left lsr 10) ,
mask10 land left ) ,
Powers.of_int_tuple
( mask10 land (right lsr 50) ,
mask10 land (right lsr 40) ,
mask10 land (right lsr 30)),
Powers.of_int_tuple
( mask10 land (right lsr 20) ,
mask10 land (right lsr 10) ,
mask10 land right )
)
| _ -> invalid_arg (__FILE__^": to_powers")
#+end_src
** Functions for hash tables
#+begin_src ocaml :tangle (eval mli)
val hash : t -> int
val equal : t -> t -> bool
val compare : t -> t -> int
#+end_src
| ~hash~ | Associates a nonnegative integer to any Zkey |
| ~equal~ | The equal function. True if two Zkeys are equal |
| ~compare~ | Comparison function, used for sorting |
#+begin_src ocaml :tangle (eval ml) :exports none
let hash = Hashtbl.hash
let equal
{ right = r1 ; left = l1 ; kind = k1 }
{ right = r2 ; left = l2 ; kind = k2 } =
r1 = r2 && l1 = l2 && k1 = k2
let compare
{ right = r1 ; left = l1 ; kind = k1 }
{ right = r2 ; left = l2 ; kind = k2 } =
if k1 <> k2 then invalid_arg (__FILE__^": cmp");
if r1 < r2 then -1
else if r1 > r2 then 1
else if l1 < l2 then -1
else if l1 > l2 then 1
else 0
let to_string { left ; right ; kind } =
"< " ^ string_of_int left ^ string_of_int right ^ " | " ^ (
to_int_array { left ; right ; kind }
|> Array.map string_of_int
|> Array.to_list
|> String.concat ", "
) ^ " >"
#+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 "@[%s@]" (to_string t)
#+end_src

30
common/zmap.org
View File

@ -1,30 +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
* Zmap
:PROPERTIES:
:header-args: :noweb yes :comments both
:END:
A hash table where the keys are ~Zkey~
** Type
<<<~Zmap.t~>>>
#+begin_src ocaml :tangle (eval mli)
include module type of Hashtbl.Make(Zkey)
#+end_src
#+begin_src ocaml :tangle (eval ml) :exports none
module Zmap = Hashtbl.Make(Zkey)
include Zmap
#+end_src

175
gaussian/atomic_shell.org
View File

@ -1,175 +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
* Atomic shell
:PROPERTIES:
:header-args: :noweb yes :comments both
:END:
Set of contracted Gaussians differing only by the powers of $x$, $y$ and $z$, with a
constant ~Angular_momentum.t~, all centered on the same point.
In other words, it is the set of all contracted shells sharing the same center.
\begin{align*}
\chi_{n_x,n_y,n_z}(r) & = f(n_x,n_y,n_z) \sum_{j=1}^{n} \sum_{i=1}^{m} \mathcal{N}_{ij}\, d_{ij}\, g_{ij\,n_x,n_y,n_z}(r) \\
& = (x-X_A)^{n_x} (y-Y_A)^{n_y} (z-Z_A)^{n_z} f(n_x,n_y,n_z) \sum_{j=1}^{n} \sum_{i=1}^{m} \mathcal{N}_{ij}\, d_{ij}\, \exp \left( -\alpha_{ij} |r-R_A|^2 \right)
\end{align*}
where:
- $g_{ij\,n_x,n_y,n_z}(r)$ is the $i$-th ~PrimitiveShell.t~ of the $j$-th ~Contracted_shell.t~
- $n_x + n_y + n_z = l$, the total angular momentum
- $\alpha_{ij}$ are the exponents (tabulated) of the $j$-th ~Contracted_shell.t~
- $d_{ij}$ are the contraction coefficients of the $j$-th ~Contracted_shell.t~
- $\mathcal{N}_{ij}$ is the normalization coefficient of the $i$-th primitive shell
(~PrimitiveShell.norm_coef~) of the $j$-th ~Contracted_shell.t~
- $f(n_x,n_y,n_z)$ is a scaling factor adjusting the normalization coefficient for the
particular powers of $x,y,z$ (~PrimitiveShell.norm_coef_scale~)
** Type
#+begin_src ocaml :tangle (eval mli)
type t
open Common
#+end_src
#+begin_src ocaml :tangle (eval ml) :exports none
open Common
type t = {
expo : float array array;
coef : float array array;
norm_coef : float array array;
norm_coef_scale : float array;
contr : Contracted_shell.t array;
index : int;
center : Coordinate.t;
ang_mom : Angular_momentum.t;
}
module Am = Angular_momentum
module Co = Coordinate
module Cs = Contracted_shell
#+end_src
** Access
#+begin_src ocaml :tangle (eval mli)
val ang_mom : t -> Angular_momentum.t
val center : t -> Coordinate.t
val coefficients : t -> float array array
val contracted_shells : t -> Contracted_shell.t array
val exponents : t -> float array array
val index : t -> int
val normalizations : t -> float array array
val norm_scales : t -> float array
val size_of_shell : t -> int
val size : t -> int
#+end_src
| ~ang_mom~ | Total angular momentum : $l = n_x + n_y + n_z$. |
| ~center~ | Coordinate of the center $\mathbf{A} = (X_A,Y_A,Z_A)$. |
| ~coefficients~ | Array of contraction coefficients $d_{ij}$. The first index is the index of the contracted function, and the second index is the index of the primitive. |
| ~contracted_shells:~ | Array of contracted gaussians |
| ~exponents~ | Array of exponents $\alpha_{ij}$. The first index is the index of the contracted function, and the second index is the index of the primitive. |
| ~index~ | Index in the basis set, represented as an array of contracted shells. |
| ~normalizations~ | Normalization coefficients $\mathcal{N}_{ij}$. The first index is the index of the contracted function, and the second index is the index of the primitive. |
| ~norm_scales~ | Scaling factors $f(n_x,n_y,n_z)$, given in the same order as ~Angular_momentum.zkey_array ang_mom~. |
| ~size~ | Number of contracted functions, $n$ in the definition. |
| ~size_of_shell~ | Number of contracted functions in the shell: length of ~norm_coef_scale~. |
#+begin_example
#+end_example
#+begin_src ocaml :tangle (eval ml) :exports none
let ang_mom t = t.ang_mom
let center t = t.center
let coefficients t = t.coef
let contracted_shells t = t.contr
let exponents t = t.expo
let index t = t.index
let normalizations t = t.norm_coef
let norm_scales t = t.norm_coef_scale
let size_of_shell t = Array.length t.norm_coef_scale
let size t = Array.length t.contr
#+end_src
** Creation
#+begin_src ocaml :tangle (eval mli)
val make : ?index:int -> Contracted_shell.t array -> t
val with_index : t -> int -> t
#+end_src
| ~make~ | Creates a contracted shell from a list of coefficients and primitives. |
| ~with_index~ | Returns a copy of the contracted shell with a modified index. |
#+begin_src ocaml :tangle (eval ml) :exports none
let make ?(index=0) contr =
assert (Array.length contr > 0);
let coef = Array.map Cs.coefficients contr
and expo = Array.map Cs.exponents contr
in
let center = Cs.center contr.(0) in
let rec unique_center = function
| 0 -> true
| i -> if Cs.center contr.(i) = center then (unique_center [@tailcall]) (i-1) else false
in
if not (unique_center (Array.length contr - 1)) then
invalid_arg "ContractedAtomicShell.make Coordinate.t differ";
let ang_mom = Cs.ang_mom contr.(0) in
let rec unique_angmom = function
| 0 -> true
| i -> if Cs.ang_mom contr.(i) = ang_mom then (unique_angmom [@tailcall]) (i-1) else false
in
if not (unique_angmom (Array.length contr - 1)) then
invalid_arg "Contracted_shell.make: Angular_momentum.t differ";
let norm_coef =
Array.map Cs.normalizations contr
in
let norm_coef_scale = Cs.norm_scales contr.(0)
in
{ index ; expo ; coef ; center ; ang_mom ; norm_coef ;
norm_coef_scale ; contr }
let with_index a i =
{ a with index = i }
#+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 s =
let open Format in
fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+ (size_of_shell s)*(size s));
fprintf ppf "@[%a@ %a@] @[" Am.pp_string s.ang_mom Co.pp s.center;
Array.iter2 (fun e_arr c_arr -> fprintf ppf "@[<v>";
Array.iter2 (fun e c -> fprintf ppf "@[%16.8e %16.8e@]@;" e c)
e_arr c_arr;
fprintf ppf "@;@]") s.expo s.coef;
fprintf ppf "@]"
#+end_src

161
gaussian/atomic_shell_pair.org
View File

@ -1,161 +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
* Atomic shell pair
:PROPERTIES:
:header-args: :noweb yes :comments both
:END:
Data structure to represent pairs of atomic shells. The products of
functions in the shell pair are one-electron functions.
An atomic shell pair is an array of pairs of contracted shells.
** Type
#+begin_src ocaml :tangle (eval mli)
type t
open Common
#+end_src
#+begin_src ocaml :tangle (eval ml) :exports none
open Common
type t =
{
contracted_shell_pairs : Contracted_shell_pair.t list;
atomic_shell_a : Atomic_shell.t;
atomic_shell_b : Atomic_shell.t;
}
module Am = Angular_momentum
module As = Atomic_shell
module Co = Coordinate
module Cs = Contracted_shell
module Csp = Contracted_shell_pair
#+end_src
** Access
#+begin_src ocaml :tangle (eval mli)
val atomic_shell_a : t -> Atomic_shell.t
val atomic_shell_b : t -> Atomic_shell.t
val contracted_shell_pairs : t -> Contracted_shell_pair.t list
val ang_mom : t -> Angular_momentum.t
val monocentric : t -> bool
val norm_scales : t -> float array
val a_minus_b : t -> Coordinate.t
val a_minus_b_sq : t -> float
#+end_src
| ~atomic_shell_a~ | Returns the first ~Atomic_shell.t~ which was used to build the atomic shell pair. |
| ~atomic_shell_b~ | Returns the second ~Atomic_shell.t~ which was used to build the atomic shell pair. |
| ~contracted_shell_pairs~ | Returns an array of ~ContractedShellPair.t~, containing all the pairs of contracted functions used to build the atomic shell pair. |
| ~norm_scales~ | norm_coef.(i) / norm_coef.(0) |
| ~ang_mom~ | Total angular Momentum |
| ~monocentric~ | If true, the two atomic shells have the same center. |
| ~a_minus_b~ | Returns $A-B$ |
| ~a_minus_b_sq~ | Returns $\vert A-B \vert^2$ |
#+begin_src ocaml :tangle (eval ml) :exports none
let atomic_shell_a x = x.atomic_shell_a
let atomic_shell_b x = x.atomic_shell_b
let contracted_shell_pairs x = x.contracted_shell_pairs
let monocentric x =
Csp.monocentric @@ List.hd x.contracted_shell_pairs
let a_minus_b x =
Csp.a_minus_b @@ List.hd x.contracted_shell_pairs
let a_minus_b_sq x =
Csp.a_minus_b_sq @@ List.hd x.contracted_shell_pairs
let ang_mom x =
Csp.ang_mom @@ List.hd x.contracted_shell_pairs
let norm_scales x =
Csp.norm_scales @@ List.hd x.contracted_shell_pairs
#+end_src
** Creation
#+begin_src ocaml :tangle (eval mli)
val make : ?cutoff:float -> Atomic_shell.t -> Atomic_shell.t -> t option
#+end_src
Creates an atomic shell pair from two atomic shells.
The contracted shell pairs contains the only pairs of primitives for which
the norm is greater than ~cutoff~.
If all the contracted shell pairs are not significant, the function returns
~None~.
#+begin_src ocaml :tangle (eval mli)
val of_atomic_shell_array : ?cutoff:float -> Atomic_shell.t array -> t option array array
#+end_src
Creates all possible atomic shell pairs from an array of atomic shells.
If an atomic shell pair is not significant, sets the value to ~None~.
#+begin_src ocaml :tangle (eval ml) :exports none
let make ?(cutoff=Constants.epsilon) atomic_shell_a atomic_shell_b =
let l_a = Array.to_list (As.contracted_shells atomic_shell_a)
and l_b = Array.to_list (As.contracted_shells atomic_shell_b)
in
let contracted_shell_pairs =
List.concat_map (fun s_a ->
List.map (fun s_b ->
if Cs.index s_b <= Cs.index s_a then
Csp.make ~cutoff s_a s_b
else
None
) l_b
) l_a
|> Util.list_some
in
match contracted_shell_pairs with
| [] -> None
| _ -> Some { atomic_shell_a ; atomic_shell_b ; contracted_shell_pairs }
let of_atomic_shell_array ?(cutoff=Constants.epsilon) basis =
Array.mapi (fun i shell_a ->
Array.map (fun shell_b ->
make ~cutoff shell_a shell_b)
(Array.sub basis 0 (i+1))
) basis
#+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 s =
let open Format in
fprintf ppf "@[%a@ %a@]"
Atomic_shell.pp s.atomic_shell_a
Atomic_shell.pp s.atomic_shell_b
#+end_src

168
gaussian/atomic_shell_pair_couple.org
View File

@ -1,168 +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
* Atomic shell pair couple
:PROPERTIES:
:header-args: :noweb yes :comments both
:END:
An atomic shell pair couple is the cartesian product between two sets of functions, one
set over electron one and one set over electron two. Both sets are atomic shell
pairs.
These are usually called /shell quartets/ in the literature, but we prefer to use
/pair/ for two functions with the same electron, and /couple/ for two functions
acting on different electrons, since they will be coupled by a two-electron operator.
** Type
#+begin_src ocaml :tangle (eval mli)
type t
open Common
#+end_src
#+begin_src ocaml :tangle (eval ml) :exports none
open Common
type t =
{
contracted_shell_pair_couples : Contracted_shell_pair_couple.t list ;
atomic_shell_pair_p: Atomic_shell_pair.t ;
atomic_shell_pair_q: Atomic_shell_pair.t ;
atomic_shell_a : Atomic_shell.t ;
atomic_shell_b : Atomic_shell.t ;
atomic_shell_c : Atomic_shell.t ;
atomic_shell_d : Atomic_shell.t ;
ang_mom : Angular_momentum.t ;
}
module Am = Angular_momentum
module Co = Coordinate
module As = Atomic_shell
module Asp = Atomic_shell_pair
module Cspc = Contracted_shell_pair_couple
#+end_src
** Access
#+begin_src ocaml :tangle (eval mli)
val ang_mom : t -> Angular_momentum.t
val atomic_shell_a : t -> Atomic_shell.t
val atomic_shell_b : t -> Atomic_shell.t
val atomic_shell_c : t -> Atomic_shell.t
val atomic_shell_d : t -> Atomic_shell.t
val atomic_shell_pair_p : t -> Atomic_shell_pair.t
val atomic_shell_pair_q : t -> Atomic_shell_pair.t
val contracted_shell_pair_couples : t -> Contracted_shell_pair_couple.t list
val monocentric : t -> bool
val norm_scales : t -> float array
val zkey_array : t -> Zkey.t array
#+end_src
| ~ang_mom~ | Total angular momentum of the shell pair couple: sum of the angular momenta of all the shells. |
| ~atomic_shell_a~ | Returns the first atomic shell of the first shell pair. |
| ~atomic_shell_b~ | Returns the second atomic shell of the first shell pair. |
| ~atomic_shell_c~ | Returns the first atomic shell of the second shell pair. |
| ~atomic_shell_d~ | Returns the second atomic shell of the second shell pair. |
| ~atomic_shell_pair_p~ | Returns the first atomic shell pair that was used to build the shell pair. |
| ~atomic_shell_pair_q~ | Returns the second atomic shell pair that was used to build the shell pair. |
| ~contracted_shell_pair_couples~ | Returns the list of significant contracted shell pair couples. |
| ~monocentric~ | True if all four atomic shells have the same center. |
| ~norm_scales~ | Scaling factors of normalization coefficients inside the shell. The ordering is the same as ~zkey_array~. |
| ~zkey_array~ | Returns the array of ~Zkey.t~ relative to the four shells of the shell pair couple. |
#+begin_src ocaml :tangle (eval ml) :exports none
let contracted_shell_pair_couples t = t.contracted_shell_pair_couples
let monocentric t =
Asp.monocentric t.atomic_shell_pair_p && Asp.monocentric t.atomic_shell_pair_q &&
As.center (Asp.atomic_shell_a t.atomic_shell_pair_p) = As.center (Asp.atomic_shell_a t.atomic_shell_pair_q)
let ang_mom t = t.ang_mom
let atomic_shell_pair_p t = t.atomic_shell_pair_p
let atomic_shell_pair_q t = t.atomic_shell_pair_q
let atomic_shell_a t = t.atomic_shell_a
let atomic_shell_b t = t.atomic_shell_b
let atomic_shell_c t = t.atomic_shell_c
let atomic_shell_d t = t.atomic_shell_d
let zkey_array t =
match t.contracted_shell_pair_couples with
| f::_ -> Cspc.zkey_array f
| _ -> invalid_arg "AtomicShellPairCouple.zkey_array"
let norm_scales t =
match t.contracted_shell_pair_couples with
| f::_ -> Cspc.norm_scales f
| _ -> invalid_arg "AtomicShellPairCouple.norm_scales"
#+end_src
** Creation
#+begin_src ocaml :tangle (eval mli)
val make : ?cutoff:float -> Atomic_shell_pair.t -> Atomic_shell_pair.t -> t option
#+end_src
Default cutoff is $\epsilon$.
| ~make~ | Creates an atomic shell pair couple using two atomic shell pairs. |
#+begin_example
#+end_example
#+begin_src ocaml :tangle (eval ml) :exports none
let make ?(cutoff=Constants.epsilon) atomic_shell_pair_p atomic_shell_pair_q =
let ang_mom =
Am.(Asp.ang_mom atomic_shell_pair_p + Asp.ang_mom atomic_shell_pair_q)
in
let atomic_shell_a = Asp.atomic_shell_a atomic_shell_pair_p
and atomic_shell_b = Asp.atomic_shell_b atomic_shell_pair_p
and atomic_shell_c = Asp.atomic_shell_a atomic_shell_pair_q
and atomic_shell_d = Asp.atomic_shell_b atomic_shell_pair_q
in
let contracted_shell_pair_couples =
List.concat_map (fun ap_ab ->
List.map (fun ap_cd ->
Cspc.make ~cutoff ap_ab ap_cd
) (Asp.contracted_shell_pairs atomic_shell_pair_q)
) (Asp.contracted_shell_pairs atomic_shell_pair_p)
|> Util.list_some
in
match contracted_shell_pair_couples with
| [] -> None
| _ -> Some { atomic_shell_pair_p ; atomic_shell_pair_q ; ang_mom ;
atomic_shell_a ; atomic_shell_b ; atomic_shell_c ; atomic_shell_d ;
contracted_shell_pair_couples ;
}
#+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 "[(%d,%d),(%d,%d)]"
(Atomic_shell.index t.atomic_shell_a)
(Atomic_shell.index t.atomic_shell_b)
(Atomic_shell.index t.atomic_shell_c)
(Atomic_shell.index t.atomic_shell_d)
#+end_src

42
gaussian/lib/atomic_shell.ml
View File

@ -1,4 +1,5 @@
(* [[file:~/QCaml/gaussian/atomic_shell.org::*Type][Type:2]] *)
(** Types and modules *)
open Common
type t = {
@ -15,27 +16,11 @@ type t = {
module Am = Angular_momentum
module Co = Coordinate
module Cs = Contracted_shell
(* Type:2 ends here *)
(* | ~ang_mom~ | Total angular momentum : $l = n_x + n_y + n_z$. |
* | ~center~ | Coordinate of the center $\mathbf{A} = (X_A,Y_A,Z_A)$. |
* | ~coefficients~ | Array of contraction coefficients $d_{ij}$. The first index is the index of the contracted function, and the second index is the index of the primitive. |
* | ~contracted_shells:~ | Array of contracted gaussians |
* | ~exponents~ | Array of exponents $\alpha_{ij}$. The first index is the index of the contracted function, and the second index is the index of the primitive. |
* | ~index~ | Index in the basis set, represented as an array of contracted shells. |
* | ~normalizations~ | Normalization coefficients $\mathcal{N}_{ij}$. The first index is the index of the contracted function, and the second index is the index of the primitive. |
* | ~norm_scales~ | Scaling factors $f(n_x,n_y,n_z)$, given in the same order as ~Angular_momentum.zkey_array ang_mom~. |
* | ~size~ | Number of contracted functions, $n$ in the definition. |
* | ~size_of_shell~ | Number of contracted functions in the shell: length of ~norm_coef_scale~. |
*
* #+begin_example
*
* #+end_example *)
(** Access *)
(* [[file:~/QCaml/gaussian/atomic_shell.org::*Access][Access:2]] *)
let ang_mom t = t.ang_mom
let center t = t.center
let coefficients t = t.coef
@ -46,16 +31,9 @@ let normalizations t = t.norm_coef
let norm_scales t = t.norm_coef_scale
let size_of_shell t = Array.length t.norm_coef_scale
let size t = Array.length t.contr
(* Access:2 ends here *)
(* | ~make~ | Creates a contracted shell from a list of coefficients and primitives. |
* | ~with_index~ | Returns a copy of the contracted shell with a modified index. | *)
(* [[file:~/QCaml/gaussian/atomic_shell.org::*Creation][Creation:2]] *)
let make ?(index=0) contr =
let make ?(index=0) contr =
assert (Array.length contr > 0);
let coef = Array.map Cs.coefficients contr
@ -69,7 +47,7 @@ let make ?(index=0) contr =
in
if not (unique_center (Array.length contr - 1)) then
invalid_arg "ContractedAtomicShell.make Coordinate.t differ";
let ang_mom = Cs.ang_mom contr.(0) in
let rec unique_angmom = function
| 0 -> true
@ -77,9 +55,9 @@ let make ?(index=0) contr =
in
if not (unique_angmom (Array.length contr - 1)) then
invalid_arg "Contracted_shell.make: Angular_momentum.t differ";
let norm_coef =
Array.map Cs.normalizations contr
Array.map Cs.normalizations contr
in
let norm_coef_scale = Cs.norm_scales contr.(0)
in
@ -89,9 +67,10 @@ let make ?(index=0) contr =
let with_index a i =
{ a with index = i }
(* Creation:2 ends here *)
(* [[file:~/QCaml/gaussian/atomic_shell.org::*Printers][Printers:2]] *)
(** Printers *)
let pp ppf s =
let open Format in
fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+ (size_of_shell s)*(size s));
@ -101,4 +80,3 @@ let pp ppf s =
e_arr c_arr;
fprintf ppf "@;@]") s.expo s.coef;
fprintf ppf "@]"
(* Printers:2 ends here *)

91
gaussian/lib/atomic_shell.mli
View File

@ -1,40 +1,97 @@
(* Type *)
(** Atomic Shell *)
(**
* Set of contracted Gaussians differing only by the powers of $x$, $y$
* and $z$, with a constant ~Angular_momentum.t~, all centered on the same
* point.
*
* In other words, it is the set of all contracted shells sharing the same
* center.
*
* \begin{align*}
* \chi_{n_x,n_y,n_z}(r) & = f(n_x,n_y,n_z) \sum_{j=1}^{n} \sum_{i=1}^{m} \mathcal{N}_{ij}\, d_{ij}\, g_{ij\,n_x,n_y,n_z}(r) \\
* & = (x-X_A)^{n_x} (y-Y_A)^{n_y} (z-Z_A)^{n_z} f(n_x,n_y,n_z) \sum_{j=1}^{n} \sum_{i=1}^{m} \mathcal{N}_{ij}\, d_{ij}\, \exp \left( -\alpha_{ij} |r-R_A|^2 \right)
* \end{align*}
*
* where:
*
* - $g_{ij\,n_x,n_y,n_z}(r)$ is the $i$-th ~PrimitiveShell.t~ of the
* $j$-th ~Contracted_shell.t~
* - $n_x + n_y + n_z = l$, the total angular momentum
* - $\alpha_{ij}$ are the exponents (tabulated) of the $j$-th
* ~Contracted_shell.t~
* - $d_{ij}$ are the contraction coefficients of the $j$-th
* ~Contracted_shell.t~
* - $\mathcal{N}_{ij}$ is the normalization coefficient of the $i$-th
* primitive shell
* (~PrimitiveShell.norm_coef~) of the $j$-th ~Contracted_shell.t~
* - $f(n_x,n_y,n_z)$ is a scaling factor adjusting the normalization
* coefficient for the particular powers of $x,y,z$
* (~PrimitiveShell.norm_coef_scale~)
*)
(* [[file:~/QCaml/gaussian/atomic_shell.org::*Type][Type:1]] *)
(** Types and modules *)
type t
open Common
(* Type:1 ends here *)
(* Access *)
(* [[file:~/QCaml/gaussian/atomic_shell.org::*Access][Access:1]] *)
(** Access *)
val ang_mom : t -> Angular_momentum.t
(** Total angular momentum : $l = n_x + n_y + n_z$. *)
val center : t -> Coordinate.t
(** Coordinate of the center $\mathbf{A} = (X_A,Y_A,Z_A)$. *)
val coefficients : t -> float array array
(** Array of contraction coefficients $d_{ij}$. The first index is the index
* of the contracted function, and the second index is the index of the
* primitive.
*)
val contracted_shells : t -> Contracted_shell.t array
(** Array of contracted gaussians *)
val exponents : t -> float array array
(** Array of exponents $\alpha_{ij}$. The first index is the index of the
* contracted function, and the second index is the index of the primitive.
*)
val index : t -> int
(** Index in the basis set, represented as an array of contracted shells. *)
val normalizations : t -> float array array
(** Normalization coefficients $\mathcal{N}_{ij}$. The first index is the
* index of the contracted function, and the second index is the index of
* the primitive.
*)
val norm_scales : t -> float array
(** Scaling factors $f(n_x,n_y,n_z)$, given in the same order as
* ~Angular_momentum.zkey_array ang_mom~.
*)
val size_of_shell : t -> int
(** Number of contracted functions in the shell: length of
* ~norm_coef_scale~.
*)
val size : t -> int
(* Access:1 ends here *)
(* Creation *)
(** Number of contracted functions, $n$ in the definition. *)
(* [[file:~/QCaml/gaussian/atomic_shell.org::*Creation][Creation:1]] *)
val make : ?index:int -> Contracted_shell.t array -> t
(** Creation *)
val make : ?index:int -> Contracted_shell.t array -> t
(** Creates a contracted shell from a list of coefficients and primitives. *)
val with_index : t -> int -> t
(* Creation:1 ends here *)
(* Printers *)
(** Returns a copy of the contracted shell with a modified index. *)
(* [[file:~/QCaml/gaussian/atomic_shell.org::*Printers][Printers:1]] *)
(** Printers *)
val pp : Format.formatter -> t -> unit
(* Printers:1 ends here *)

29
gaussian/lib/atomic_shell_pair.ml
View File

@ -1,4 +1,4 @@
(* [[file:~/QCaml/gaussian/atomic_shell_pair.org::*Type][Type:2]] *)
(** Type *)
open Common
type t =
@ -14,21 +14,10 @@ module As = Atomic_shell
module Co = Coordinate
module Cs = Contracted_shell
module Csp = Contracted_shell_pair
(* Type:2 ends here *)
(** Access *)
(* | ~atomic_shell_a~ | Returns the first ~Atomic_shell.t~ which was used to build the atomic shell pair. |
* | ~atomic_shell_b~ | Returns the second ~Atomic_shell.t~ which was used to build the atomic shell pair. |
* | ~contracted_shell_pairs~ | Returns an array of ~ContractedShellPair.t~, containing all the pairs of contracted functions used to build the atomic shell pair. |
* | ~norm_scales~ | norm_coef.(i) / norm_coef.(0) |
* | ~ang_mom~ | Total angular Momentum |
* | ~monocentric~ | If true, the two atomic shells have the same center. |
* | ~a_minus_b~ | Returns $A-B$ |
* | ~a_minus_b_sq~ | Returns $\vert A-B \vert^2$ | *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair.org::*Access][Access:2]] *)
let atomic_shell_a x = x.atomic_shell_a
let atomic_shell_b x = x.atomic_shell_b
let contracted_shell_pairs x = x.contracted_shell_pairs
@ -47,15 +36,10 @@ let ang_mom x =
let norm_scales x =
Csp.norm_scales @@ List.hd x.contracted_shell_pairs
(* Access:2 ends here *)
(** Creation *)
(* Creates all possible atomic shell pairs from an array of atomic shells.
* If an atomic shell pair is not significant, sets the value to ~None~. *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair.org::*Creation][Creation:3]] *)
let make ?(cutoff=Constants.epsilon) atomic_shell_a atomic_shell_b =
let l_a = Array.to_list (As.contracted_shells atomic_shell_a)
@ -85,12 +69,13 @@ let of_atomic_shell_array ?(cutoff=Constants.epsilon) basis =
make ~cutoff shell_a shell_b)
(Array.sub basis 0 (i+1))
) basis
(* Creation:3 ends here *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair.org::*Printers][Printers:2]] *)
(** Printers *)
let pp ppf s =
let open Format in
fprintf ppf "@[%a@ %a@]"
Atomic_shell.pp s.atomic_shell_a
Atomic_shell.pp s.atomic_shell_b
(* Printers:2 ends here *)

59
gaussian/lib/atomic_shell_pair.mli
View File

@ -1,35 +1,52 @@
(* Type *)
(** Atomic shell pair *)
(** Data structure to represent pairs of atomic shells. The products of
* functions in the shell pair are one-electron functions.
*
* An atomic shell pair is an array of pairs of contracted shells.
*)
(* [[file:~/QCaml/gaussian/atomic_shell_pair.org::*Type][Type:1]] *)
(** Type *)
type t
open Common
(* Type:1 ends here *)
(* Access *)
(** Access *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair.org::*Access][Access:1]] *)
val atomic_shell_a : t -> Atomic_shell.t
(** Returns the first ~Atomic_shell.t~ which was used to build the atomic
* shell pair. *)
val atomic_shell_b : t -> Atomic_shell.t
(** Returns the second ~Atomic_shell.t~ which was used to build the atomic
* shell pair. *)
val contracted_shell_pairs : t -> Contracted_shell_pair.t list
(** Returns an array of ~ContractedShellPair.t~, containing all the pairs of
* contracted functions used to build the atomic shell pair.
*)
val ang_mom : t -> Angular_momentum.t
(** Total angular Momentum *)
val monocentric : t -> bool
(** If true, the two atomic shells have the same center. *)
val norm_scales : t -> float array
(** norm_coef.(i) / norm_coef.(0) *)
val a_minus_b : t -> Coordinate.t
(** Returns $A-B$ *)
val a_minus_b_sq : t -> float
(* Access:1 ends here *)
(* Creation *)
(** Returns $\vert A-B \vert^2$ *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair.org::*Creation][Creation:1]] *)
(** Creation *)
val make : ?cutoff:float -> Atomic_shell.t -> Atomic_shell.t -> t option
(* Creation:1 ends here *)
(* Creates an atomic shell pair from two atomic shells.
*
* The contracted shell pairs contains the only pairs of primitives for which
@ -38,14 +55,14 @@ val make : ?cutoff:float -> Atomic_shell.t -> Atomic_shell.t -> t option
* If all the contracted shell pairs are not significant, the function returns
* ~None~. *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair.org::*Creation][Creation:2]] *)
val of_atomic_shell_array : ?cutoff:float -> Atomic_shell.t array -> t option array array
(* Creation:2 ends here *)
(* Printers *)
(** Creates all possible atomic shell pairs from an array of atomic shells.
* If an atomic shell pair is not significant, sets the value to ~None~.
*)
(* [[file:~/QCaml/gaussian/atomic_shell_pair.org::*Printers][Printers:1]] *)
(** Printers *)
val pp : Format.formatter -> t -> unit
(* Printers:1 ends here *)

48
gaussian/lib/atomic_shell_pair_couple.ml
View File

@ -1,7 +1,10 @@
(* [[file:~/QCaml/gaussian/atomic_shell_pair_couple.org::*Type][Type:2]] *)
(** Atomic shell pair couple *)
open Common
type t =
(** Type *)
type t =
{
contracted_shell_pair_couples : Contracted_shell_pair_couple.t list ;
atomic_shell_pair_p: Atomic_shell_pair.t ;
@ -18,24 +21,10 @@ module Co = Coordinate
module As = Atomic_shell
module Asp = Atomic_shell_pair
module Cspc = Contracted_shell_pair_couple
(* Type:2 ends here *)
(** Access *)
(* | ~ang_mom~ | Total angular momentum of the shell pair couple: sum of the angular momenta of all the shells. |
* | ~atomic_shell_a~ | Returns the first atomic shell of the first shell pair. |
* | ~atomic_shell_b~ | Returns the second atomic shell of the first shell pair. |
* | ~atomic_shell_c~ | Returns the first atomic shell of the second shell pair. |
* | ~atomic_shell_d~ | Returns the second atomic shell of the second shell pair. |
* | ~atomic_shell_pair_p~ | Returns the first atomic shell pair that was used to build the shell pair. |
* | ~atomic_shell_pair_q~ | Returns the second atomic shell pair that was used to build the shell pair. |
* | ~contracted_shell_pair_couples~ | Returns the list of significant contracted shell pair couples. |
* | ~monocentric~ | True if all four atomic shells have the same center. |
* | ~norm_scales~ | Scaling factors of normalization coefficients inside the shell. The ordering is the same as ~zkey_array~. |
* | ~zkey_array~ | Returns the array of ~Zkey.t~ relative to the four shells of the shell pair couple. | *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair_couple.org::*Access][Access:2]] *)
let contracted_shell_pair_couples t = t.contracted_shell_pair_couples
let monocentric t =
@ -59,24 +48,15 @@ let zkey_array t =
| f::_ -> Cspc.zkey_array f
| _ -> invalid_arg "AtomicShellPairCouple.zkey_array"
let norm_scales t =
let norm_scales t =
match t.contracted_shell_pair_couples with
| f::_ -> Cspc.norm_scales f
| _ -> invalid_arg "AtomicShellPairCouple.norm_scales"
(* Access:2 ends here *)
(* Default cutoff is $\epsilon$.
*
* | ~make~ | Creates an atomic shell pair couple using two atomic shell pairs. |
*
* #+begin_example
*
* #+end_example *)
(** Creation *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair_couple.org::*Creation][Creation:2]] *)
let make ?(cutoff=Constants.epsilon) atomic_shell_pair_p atomic_shell_pair_q =
let make ?(cutoff=Constants.epsilon) atomic_shell_pair_p atomic_shell_pair_q =
let ang_mom =
Am.(Asp.ang_mom atomic_shell_pair_p + Asp.ang_mom atomic_shell_pair_q)
in
@ -86,26 +66,24 @@ let make ?(cutoff=Constants.epsilon) atomic_shell_pair_p atomic_shell_pair_q =
and atomic_shell_d = Asp.atomic_shell_b atomic_shell_pair_q
in
let contracted_shell_pair_couples =
List.concat_map (fun ap_ab ->
List.map (fun ap_cd ->
List.concat_map (fun ap_ab ->
List.map (fun ap_cd ->
Cspc.make ~cutoff ap_ab ap_cd
) (Asp.contracted_shell_pairs atomic_shell_pair_q)
) (Asp.contracted_shell_pairs atomic_shell_pair_p)
|> Util.list_some
in
match contracted_shell_pair_couples with
| [] -> None
| [] -> None
| _ -> Some { atomic_shell_pair_p ; atomic_shell_pair_q ; ang_mom ;
atomic_shell_a ; atomic_shell_b ; atomic_shell_c ; atomic_shell_d ;
contracted_shell_pair_couples ;
}
(* Creation:2 ends here *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair_couple.org::*Printers][Printers:2]] *)
(** Printers *)
let pp ppf t =
Format.fprintf ppf "[(%d,%d),(%d,%d)]"
(Atomic_shell.index t.atomic_shell_a)
(Atomic_shell.index t.atomic_shell_b)
(Atomic_shell.index t.atomic_shell_c)
(Atomic_shell.index t.atomic_shell_d)
(* Printers:2 ends here *)

69
gaussian/lib/atomic_shell_pair_couple.mli
View File

@ -1,39 +1,78 @@
(* Type *)
(** Atomic shell pair couple *)
(** An atomic shell pair couple is the cartesian product between two sets of
* functions, one set over electron one and one set over electron two. Both
* sets are atomic shell pairs.
*)
(** These are usually called /shell quartets/ in the literature, but we
* prefer to use /pair/ for two functions with the same electron, and
* /couple/ for two functions acting on different electrons, since they will
* be coupled by a two-electron operator.
*)
(* [[file:~/QCaml/gaussian/atomic_shell_pair_couple.org::*Type][Type:1]] *)
(** Type *)
type t
open Common
(* Type:1 ends here *)
(* Access *)
(* [[file:~/QCaml/gaussian/atomic_shell_pair_couple.org::*Access][Access:1]] *)
(** Access *)
val ang_mom : t -> Angular_momentum.t
(** Total angular momentum of the shell pair couple: sum of the angular
* momenta of all the shells.
*)
val atomic_shell_a : t -> Atomic_shell.t
(** Returns the first atomic shell of the first shell pair. *)
val atomic_shell_b : t -> Atomic_shell.t
(** Returns the second atomic shell of the first shell pair. *)
val atomic_shell_c : t -> Atomic_shell.t
(** Returns the first atomic shell of the second shell pair. *)
val atomic_shell_d : t -> Atomic_shell.t
(** Returns the second atomic shell of the second shell pair. *)
val atomic_shell_pair_p : t -> Atomic_shell_pair.t
(** Returns the first atomic shell pair that was used to build the shell
* pair.
*)
val atomic_shell_pair_q : t -> Atomic_shell_pair.t
(** Returns the second atomic shell pair that was used to build the shell
* pair.
*)
val contracted_shell_pair_couples : t -> Contracted_shell_pair_couple.t list
(** Returns the list of significant contracted shell pair couples. *)
val monocentric : t -> bool
(** True if all four atomic shells have the same center. *)
val norm_scales : t -> float array
(** Scaling factors of normalization coefficients inside the shell. The
* ordering is the same as ~zkey_array~.
*)
val zkey_array : t -> Zkey.t array
(* Access:1 ends here *)
(* Creation *)
(** Returns the array of ~Zkey.t~ relative to the four shells of the shell
* pair couple.
*)
(* [[file:~/QCaml/gaussian/atomic_shell_pair_couple.org::*Creation][Creation:1]] *)
(** Creation *)
val make : ?cutoff:float -> Atomic_shell_pair.t -> Atomic_shell_pair.t -> t option
(* Creation:1 ends here *)
(* Printers *)
(** Creates an atomic shell pair couple using two atomic shell pairs.
* Default cutoff is $\epsilon$.
*)
(* [[file:~/QCaml/gaussian/atomic_shell_pair_couple.org::*Printers][Printers:1]] *)
(** Printers *)
val pp : Format.formatter -> t -> unit
(* Printers:1 ends here *)

118
mo/frozen_core.org
View File

@ -1,118 +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
* Frozen core
:PROPERTIES:
:header-args: :noweb yes :comments both
:END:
Defines how the core electrons are frozen, for each atom.
** Type
#+NAME: types
#+begin_src ocaml :tangle (eval mli)
type kind =
| All_electron
| Small
| Large
#+end_src
#+begin_src ocaml :tangle (eval mli)
type t
#+end_src
#+begin_src ocaml :tangle (eval ml) :exports none
<<types>>
type t = int array
#+end_src
** Creation
#+begin_src ocaml :tangle (eval mli)
val make : kind -> Particles.Nuclei.t -> t
val of_int_list : int list -> t
val of_int_array : int array -> t
#+end_src
| ~make~ | Creates a ~Frozen_core.t~ with the same kind for all atoms |
| ~of_int_array~ | Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom |
| ~of_int_list~ | Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom |
#+begin_example
let f = Frozen_core.(make Small nuclei) ;;
val f : Frozen_core.t = [|0; 2; 2; 0|]
let f = Frozen_core.(of_int_list [0; 2; 2; 0])
val f : Frozen_core.t = [|0; 2; 2; 0|]
#+end_example
#+begin_src ocaml :tangle (eval ml) :exports none
let make_ae nuclei =
Array.map (fun _ -> 0) nuclei
let make_small nuclei =
Array.map (fun (e,_) -> Particles.Element.small_core e) nuclei
let make_large nuclei =
Array.map (fun (e,_) -> Particles.Element.large_core e) nuclei
let make = function
| All_electron -> make_ae
| Small -> make_small
| Large -> make_large
external of_int_array : int array -> t = "%identity"
let of_int_list = Array.of_list
#+end_src
** Access
#+begin_src ocaml :tangle (eval mli)
val num_elec : t -> int
val num_mos : t -> int
#+end_src
| ~num_elec~ | Number of frozen electrons |
| ~num_mos~ | Number of frozen molecular orbitals |
#+begin_example
Frozen_core.num_elec f ;;
- : int = 4
Frozen_core.num_mos f ;;
- : int = 2
#+end_example
#+begin_src ocaml :tangle (eval ml) :exports none
let num_elec t =
Array.fold_left ( + ) 0 t
let num_mos t =
(num_elec t) / 2
#+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 "@[[|";
Array.iter (fun x -> Format.fprintf ppf "@,@[%d@]" x) t;
Format.fprintf ppf "|]@]";
#+end_src

45
mo/lib/frozen_core.ml
View File

@ -1,27 +1,14 @@
(* [[file:~/QCaml/mo/frozen_core.org::*Type][Type:3]] *)
(** Type *)
type kind =
| All_electron
| Small
| Large
type t = int array
(* Type:3 ends here *)
(** Creation *)
(* | ~make~ | Creates a ~Frozen_core.t~ with the same kind for all atoms |
* | ~of_int_array~ | Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom |
* | ~of_int_list~ | Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom |
*
* #+begin_example
* let f = Frozen_core.(make Small nuclei) ;;
* val f : Frozen_core.t = [|0; 2; 2; 0|]
*
* let f = Frozen_core.(of_int_list [0; 2; 2; 0])
* val f : Frozen_core.t = [|0; 2; 2; 0|]
* #+end_example *)
(* [[file:~/QCaml/mo/frozen_core.org::*Creation][Creation:2]] *)
let make_ae nuclei =
Array.map (fun _ -> 0) nuclei
@ -32,41 +19,29 @@ let make_large nuclei =
Array.map (fun (e,_) -> Particles.Element.large_core e) nuclei
let make = function
| All_electron -> make_ae
| Small -> make_small
| All_electron -> make_ae
| Small -> make_small
| Large -> make_large
external of_int_array : int array -> t = "%identity"
let of_int_list = Array.of_list
(* Creation:2 ends here *)
(** Access *)
(* | ~num_elec~ | Number of frozen electrons |
* | ~num_mos~ | Number of frozen molecular orbitals |
*
* #+begin_example
* Frozen_core.num_elec f ;;
* - : int = 4
*
* Frozen_core.num_mos f ;;
* - : int = 2
* #+end_example *)
(* [[file:~/QCaml/mo/frozen_core.org::*Access][Access:2]] *)
let num_elec t =
Array.fold_left ( + ) 0 t
let num_mos t =
(num_elec t) / 2
(* Access:2 ends here *)
(* [[file:~/QCaml/mo/frozen_core.org::*Printers][Printers:2]] *)
(** Printers *)
let pp ppf t =
Format.fprintf ppf "@[[|";
Array.iter (fun x -> Format.fprintf ppf "@,@[%d@]" x) t;
Format.fprintf ppf "|]@]";
(* Printers:2 ends here *)

54
mo/lib/frozen_core.mli
View File

@ -1,39 +1,55 @@
(* Type
*
* #+NAME: types *)
(** Type *)
(** Defines how the core electrons are frozen, for each atom. *)
(* [[file:~/QCaml/mo/frozen_core.org::types][types]] *)
type kind =
| All_electron
| Small
| Large
(* types ends here *)
(* [[file:~/QCaml/mo/frozen_core.org::*Type][Type:2]] *)
type t
(* Type:2 ends here *)
(* Creation *)
(* [[file:~/QCaml/mo/frozen_core.org::*Creation][Creation:1]] *)
(** Creation *)
(** Example
*
* let f = Frozen_core.(make Small nuclei) ;;
* val f : Frozen_core.t = [|0; 2; 2; 0|]
*
* let f = Frozen_core.(of_int_list [0; 2; 2; 0])
* val f : Frozen_core.t = [|0; 2; 2; 0|]
*
*)
val make : kind -> Particles.Nuclei.t -> t
(** Creates a ~Frozen_core.t~ with the same kind for all atoms *)
val of_int_list : int list -> t
(** Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom *)
val of_int_array : int array -> t
(* Creation:1 ends here *)
(* Access *)
(** Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom *)
(* [[file:~/QCaml/mo/frozen_core.org::*Access][Access:1]] *)
(** Access *)
(** Example
*
* Frozen_core.num_elec f ;;
* - : int = 4
*
* Frozen_core.num_mos f ;;
*
*)
val num_elec : t -> int
(** Number of frozen electrons *)
val num_mos : t -> int
(* Access:1 ends here *)
(** Number of frozen molecular orbitals *)
(* Printers *)
(** Printers *)
(* [[file:~/QCaml/mo/frozen_core.org::*Printers][Printers:1]] *)
val pp : Format.formatter -> t -> unit
(* Printers:1 ends here *)

86
mo/lib/localization.ml
View File

@ -1,6 +1,8 @@
(* [[file:~/QCaml/mo/localization.org::*Type][Type:3]] *)
open Linear_algebra
open Common
(** Types *)
type localization_kind =
| Edmiston
| Boys
@ -18,7 +20,7 @@ type localization_data =
convergence : float ;
iteration : int ;
}
type t =
{
kind : localization_kind ;
@ -26,14 +28,10 @@ type t =
data : localization_data option lazy_t array ;
selected_mos : int list ;
}
open Common
(* Type:3 ends here *)
(* Edmiston-Rudenberg *)
(* [[file:~/QCaml/mo/localization.org::*Edmiston-Rudenberg][Edmiston-Rudenberg:1]] *)
(** Edmiston-Rudenberg *)
let kappa_edmiston in_basis m_C =
let ao_basis =
Basis.simulation in_basis
@ -117,7 +115,7 @@ let kappa_edmiston in_basis m_C =
) (Util.array_range 1 n_ao);
let f i j =
let f i j =
if i=j then
0.
else
@ -133,12 +131,10 @@ let kappa_edmiston in_basis m_C =
Matrix.init_cols n_mo n_mo ( fun i j -> if i<=j then f i j else -. (f j i) ),
Vector.sum (Vector.of_bigarray_inplace v_d)
)
(* Edmiston-Rudenberg:1 ends here *)
(* Boys *)
(* [[file:~/QCaml/mo/localization.org::*Boys][Boys:1]] *)
(** Boys *)
let kappa_boys in_basis =
let ao_basis =
Basis.simulation in_basis
@ -147,11 +143,11 @@ let kappa_boys in_basis =
let multipole = Ao.Basis.multipole ao_basis in
fun m_C ->
let n_mo = Matrix.dim2 m_C in
let phi_x_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "x") in
let phi_y_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "y") in
let phi_z_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "z") in
let m_b12 =
let g x i j =
let x_ii = x%:(i,i) in
@ -160,13 +156,13 @@ let kappa_boys in_basis =
(x_ii -. x_jj) *. x_ij
in
Matrix.init_cols n_mo n_mo (fun i j ->
let x =
(g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j)
let x =
(g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j)
in
if (abs_float x > 1.e-15) then x else 0.
)
)
in
let m_a12 =
let g x i j =
let x_ii = x%:(i,i) in
@ -176,18 +172,18 @@ let kappa_boys in_basis =
(x_ij *. x_ij) -. 0.25 *. y *. y
in
Matrix.init_cols n_mo n_mo (fun i j ->
let x =
(g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j)
let x =
(g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j)
in
if (abs_float x > 1.e-15) then x else 0.
)
)
in
let f i j =
let f i j =
let pi = Constants.pi in
if i=j then
0.
else
else
let x = atan2 (m_b12%:(i,j)) (m_a12%:(i,j)) in
if x >= 0. then
0.25 *. (pi -. x)
@ -205,11 +201,9 @@ let kappa_boys in_basis =
(* | ~kappa~ | Returns the $\kappa$ antisymmetric matrix used for the rotation matrix and the value of the localization function |
* | ~make~ | Performs the orbital localization | *)
(** Access *)
(* [[file:~/QCaml/mo/localization.org::*Access][Access:2]] *)
let kind t = t.kind
let simulation t = Basis.simulation t.mo_basis
let selected_mos t = t.selected_mos
@ -218,7 +212,7 @@ let kappa ~kind =
match kind with
| Edmiston -> kappa_edmiston
| Boys -> kappa_boys
let n_iterations t =
Array.fold_left (fun accu x ->
@ -238,20 +232,20 @@ let ao_basis t = Simulation.ao_basis (simulation t)
let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos =
let kappa_loc = kappa ~kind mo_basis in
let mo_array = Matrix.to_col_vecs (Basis.mo_coef mo_basis) in
let mos_vec_list = List.map (fun i -> mo_array.(i-1)) selected_mos in
let x: (ao,loc) Matrix.t =
Matrix.of_col_vecs_list mos_vec_list
Matrix.of_col_vecs_list mos_vec_list
in
let next_coef kappa x =
let next_coef kappa x =
let r = Matrix.exponential_antisymmetric kappa in
let m_C = Matrix.gemm_nt x r in
m_C
in
let data_initial =
let data_initial =
let iteration = 0
and scaling = 0.1
in
@ -262,7 +256,7 @@ let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos =
{ coefficients ; kappa ; scaling ; convergence ; loc_value ; iteration }
in
let iteration data =
let iteration data =
let iteration = data.iteration+1 in
let new_kappa, new_loc_value = kappa_loc data.coefficients in
let new_convergence = abs_float (Matrix.amax new_kappa) in
@ -318,33 +312,34 @@ let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos =
Array.init max_iter (fun i -> lazy (loop i))
in
{ kind ; mo_basis ; data = array_data ; selected_mos }
let to_basis t =
let to_basis t =
let mo_basis = t.mo_basis in
let simulation = Basis.simulation mo_basis in
let mo_occupation = Basis.mo_occupation mo_basis in
let data = last_iteration t in
let mo_coef_array = Matrix.to_col_vecs (Basis.mo_coef mo_basis) in
let new_mos =
Matrix.to_col_vecs data.coefficients
Matrix.to_col_vecs data.coefficients
in
selected_mos t
|> List.iteri (fun i j -> mo_coef_array.(j-1) <- new_mos.(i)) ;
let mo_coef = Matrix.of_col_vecs mo_coef_array in
Basis.make ~simulation ~mo_type:(Localized "Boys") ~mo_occupation ~mo_coef ()
(* Access:2 ends here *)
(* [[file:~/QCaml/mo/localization.org::*Printers][Printers:2]] *)
(** Printers *)
let linewidth = 60
let pp_iterations ppf t =
let line = (String.make linewidth '-') in
Format.fprintf ppf "@[%4s%s@]@." "" line;
Format.fprintf ppf "@[%4s@[%5s@]@,@[%16s@]@,@[%16s@]@,@[%11s@]@]@."
Format.fprintf ppf "@[%4s@[%5s@]@,@[%16s@]@,@[%16s@]@,@[%11s@]@]@."
"" "#" "Localization " "Convergence" "Scaling";
Format.fprintf ppf "@[%4s%s@]@." "" line;
Array.iter (fun data ->
@ -372,4 +367,3 @@ let pp ppf t =
) "MO Localization";
Format.fprintf ppf "@[%s@]@.@." (String.make 70 '=');
Format.fprintf ppf "@[%a@]@." pp_iterations t;
(* Printers:2 ends here *)

33
mo/lib/localization.mli
View File

@ -1,10 +1,10 @@
(* Type
*
* #+NAME: types *)
(** Orbital localization *)
(** Types *)
(* [[file:~/QCaml/mo/localization.org::types][types]] *)
open Linear_algebra
type localization_kind =
| Edmiston
| Boys
@ -12,17 +12,13 @@ type localization_kind =
type mo = Mo_dim.t
type ao = Ao.Ao_dim.t
type loc
(* types ends here *)
(* [[file:~/QCaml/mo/localization.org::*Type][Type:2]] *)
type localization_data
type localization_data
type t
(* Type:2 ends here *)
(* Access *)
(* [[file:~/QCaml/mo/localization.org::*Access][Access:1]] *)
(** Access *)
val kind : t -> localization_kind
val simulation : t -> Simulation.t
val selected_mos : t -> int list
@ -32,21 +28,20 @@ val kappa :
Basis.t ->
( ao,loc) Matrix.t ->
(loc,loc) Matrix.t * float
(** Returns the $\kappa$ antisymmetric matrix used for the rotation matrix and the value of the localization function *)
val make :
kind:localization_kind ->
?max_iter:int ->
?convergence:float ->
?max_iter:int ->
?convergence:float ->
Basis.t ->
int list ->
t
(** Performs the orbital localization *)
val to_basis : t -> Basis.t
(* Access:1 ends here *)
(* Printers *)
(* [[file:~/QCaml/mo/localization.org::*Printers][Printers:1]] *)
(** Printers *)
val pp : Format.formatter -> t -> unit
(* Printers:1 ends here *)

72
mo/test/localization.ml Normal file
View File

@ -0,0 +1,72 @@
(* Tests *)
let test_localization =
let nuclei =
Particles.Nuclei.of_xyz_string
" 10
Hydrogen chain, d=1.8 Angstrom
H -4.286335 0.000000 0.000000
H -3.333816 0.000000 0.000000
H -2.381297 0.000000 0.000000
H -1.428778 0.000000 0.000000
H -0.476259 0.000000 0.000000
H 0.476259 0.000000 0.000000
H 1.428778 0.000000 0.000000
H 2.381297 0.000000 0.000000
H 3.333816 0.000000 0.000000
H 4.286335 0.000000 0.000000
" in
let basis_file = "/home/scemama/qp2/data/basis/sto-6g" in
let ao_basis =
Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file
in
let charge = 0 in
let multiplicity = 1 in
let simulation =
Simulation.make ~charge ~multiplicity ~nuclei ao_basis
in
let hf =
Mo.Hartree_fock.make ~guess:`Hcore simulation
in
let mo_basis =
Mo.Basis.of_hartree_fock hf
in
let localized_mo_basis =
Mo.Localization.make
~kind:Mo.Localization.Boys
mo_basis
[4;5;6;7;8]
|> Mo.Localization.to_basis
in
Format.printf "%a" (Mo.Basis.pp ~start:1 ~finish:10) localized_mo_basis
(*
open Common
open Alcotest
let wd = Qcaml.root ^ Filename.dir_sep ^ "test" in
let test_xyz molecule length repulsion charge core =
let xyz = Nuclei.of_xyz_file (wd^Filename.dir_sep^molecule^".xyz") in
check int "length" length (Array.length xyz);
check (float 1.e-4) "repulsion" repulsion (Nuclei.repulsion xyz);
check int "charge" charge (Charge.to_int @@ Nuclei.charge xyz);
check int "small_core" core (Nuclei.small_core xyz);
()
let tests = [
"caffeine", `Quick, (fun () -> test_xyz "caffeine" 24 917.0684 102 28);
"water", `Quick, (fun () -> test_xyz "water" 3 9.19497 10 2);
]
*)

14
top/lib/install_printers.ml
View File

@ -1,15 +1,6 @@
(* [[file:~/QCaml/top/install_printers.org::*Intall printers][Intall printers:3]] *)
(** Intall printers printers:3]] *)
let printers =
[
"Common.Charge.pp" ;
"Common.Coordinate.pp" ;
"Common.Powers.pp" ;
"Common.Range.pp" ;
"Common.Spin.pp" ;
"Common.Zkey.pp" ;
"Gaussian.Atomic_shell.pp" ;
"Gaussian.Atomic_shell_pair.pp" ;
"Gaussian.Atomic_shell_pair_couple.pp" ;
"Mo.Frozen_core.pp" ;
"Mo.Localization.pp" ;
"Particles.Electrons.pp" ;
@ -18,7 +9,7 @@ let printers =
"Particles.Zmatrix.pp" ;
"Perturbation.Mp2.pp" ;
"Simulation.pp" ;
]
let eval_exn str =
@ -36,4 +27,3 @@ let rec install_printers = function
let () =
if not (install_printers printers) then
Format.eprintf "Problem installing QCaml-printers@."
(* Intall printers:3 ends here *)