From b4211b8d85c8d1421d501a8e9b3a72574022acc6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 Apr 2023 19:01:42 +0200 Subject: [PATCH] Working on TREXIO basis --- DevInstructions.org | 20 +++ README.md | 11 +- ao/README.org | 1 + ao/basis.org | 97 ++++++++++---- ao/basis_gaussian.org | 32 ++--- ao/lib/basis.ml | 67 +++++++--- ao/lib/basis.mli | 20 ++- ao/lib/basis_gaussian.ml | 22 ++-- ao/lib/basis_gaussian.mli | 2 +- ao/lib/basis_poly.ml | 2 +- ao/test/dune | 1 + docs/ao.html | 137 ++++++++++++-------- docs/gaussian.html | 104 +++++++-------- docs/particles.html | 224 ++++++++++++++++++--------------- docs/top.html | 10 +- gaussian/lib/basis.ml | 103 +++++++++++++-- gaussian/lib/basis.mli | 7 ++ particles/lib/nuclei.ml | 69 +++++----- particles/lib/nuclei.mli | 17 ++- particles/lib/nuclei_lexer.mll | 4 +- particles/nuclei.org | 140 +++++++++++---------- 21 files changed, 680 insertions(+), 410 deletions(-) create mode 100644 DevInstructions.org diff --git a/DevInstructions.org b/DevInstructions.org new file mode 100644 index 0000000..8663958 --- /dev/null +++ b/DevInstructions.org @@ -0,0 +1,20 @@ +#+TITLE: Developper's instructions + +* General + +The developper writes org-mode files, and the OCaml files are produced +by running the ~make~ command. Then, the OCaml files are compiled +using ~dune~. + +In each directory, the =README.org= file is the entry point. It +contains the script to generate the dune files, which should be +created by executing ~Ctrl-C Ctrl-C~ in the appropriate code block. + +If necessary, some pure OCaml files can be added in the =lib= +directory. + +* Basis + +- =Basis_poly= is a polymorphic basis +- =Basis_gaussian= contains the basis and integrals +- =Gaussian.Basis= is the actual Gaussian basis diff --git a/README.md b/README.md index c000cf8..708ba14 100644 --- a/README.md +++ b/README.md @@ -25,15 +25,6 @@ To use the Intel MKL library: ```bash export LACAML_LIBS="-L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_rt -lpthread -lm -ldl" opam install --assume-depexts lacaml - - -odoc-ltxhtml allows to embed equations in the documentation generated by Ocamldoc. -Download the source code [here](https://github.com/scemama/odoc-ltxhtml). - -```bash -git clone https://github.com/scemama/odoc-ltxhtml -cd odoc-ltxhtml -make install ``` How to build the project @@ -42,7 +33,7 @@ How to build the project Run `make` to compile the libraries and executables that are meant to be installed. ``` -$ make +$ make -j ``` How to run tests diff --git a/ao/README.org b/ao/README.org index 2fb54c6..13a5d1b 100644 --- a/ao/README.org +++ b/ao/README.org @@ -50,6 +50,7 @@ with open(dunetest,'w') as f: (synopsis "Test for {name} library") (libraries alcotest + trexio qcaml.{name} ) ) diff --git a/ao/basis.org b/ao/basis.org index 3ea83fc..01a3a57 100644 --- a/ao/basis.org +++ b/ao/basis.org @@ -7,7 +7,7 @@ (setq ml (concat lib name ".ml")) (setq test-ml (concat testdir name ".ml")) (org-babel-tangle) -#+end_src +#+end_src * Basis :PROPERTIES: @@ -15,7 +15,7 @@ :END: Data structure for Atomic Orbitals. - + ** Dimensions :noexport: #+begin_src ocaml :tangle lib/ao_dim.mli :exports none @@ -41,7 +41,7 @@ type t = <<<~Basis.t~>>> #+begin_src ocaml :tangle (eval mli) type t -type ao = Ao_dim.t +type ao = Ao_dim.t open Common open Particles @@ -58,7 +58,7 @@ type t = type ao = Ao_dim.t open Linear_algebra -open Common +open Common #+end_src ** Conversions @@ -87,8 +87,11 @@ val of_nuclei_and_basis_filename : let b = Ao.Basis.of_nuclei_and_basis_filename ~nuclei filename;; val b : Ao.Basis.t = Gaussian Basis, spherical, 15 AOs #+end_example - + #+begin_src ocaml :tangle (eval ml) :exports none +let not_implemented () = + Util.not_implemented "Only Gaussian is implemented" + let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false) ~nuclei filename = match kind with @@ -96,14 +99,13 @@ let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false) let basis = Gaussian.Basis.of_nuclei_and_basis_filename ~nuclei filename in - let ao_basis = + let ao_basis = Basis_poly.Gaussian (Basis_gaussian.make ~basis ?operators ~cartesian nuclei ) in { ao_basis ; cartesian } - | _ -> failwith "of_nuclei_and_basis_filename needs to be called with `Gaussian" + | _ -> not_implemented () #+end_src - ** Access #+begin_src ocaml :tangle (eval mli) @@ -140,12 +142,9 @@ val values : t -> Coordinate.t -> ao Vector.t #+begin_src ocaml :tangle (eval ml) :exports none -let not_implemented () = - Util.not_implemented "Only Gaussian is implemented" - let ao_basis t = t.ao_basis - -let size t = + +let size t = match t.ao_basis with | Basis_poly.Gaussian b -> Basis_gaussian.size b | _ -> not_implemented () @@ -157,7 +156,7 @@ let overlap t = | _ -> not_implemented () end |> Matrix.relabel - + let multipole t = begin match t.ao_basis with @@ -168,7 +167,7 @@ let multipole t = |> Matrix.relabel | _ -> not_implemented () end - + let ortho t = begin match t.ao_basis with @@ -176,7 +175,7 @@ let ortho t = | _ -> not_implemented () end |> Matrix.relabel - + let eN_ints t = begin match t.ao_basis with @@ -184,7 +183,7 @@ let eN_ints t = | _ -> not_implemented () end |> Matrix.relabel - + let kin_ints t = begin match t.ao_basis with @@ -192,7 +191,7 @@ let kin_ints t = | _ -> not_implemented () end |> Matrix.relabel - + let ee_ints t = begin match t.ao_basis with @@ -200,7 +199,7 @@ let ee_ints t = | _ -> not_implemented () end |> Four_idx_storage.relabel - + let ee_lr_ints t = begin match t.ao_basis with @@ -208,7 +207,7 @@ let ee_lr_ints t = | _ -> not_implemented () end |> Four_idx_storage.relabel - + let f12_ints t = begin match t.ao_basis with @@ -216,7 +215,7 @@ let f12_ints t = | _ -> not_implemented () end |> Four_idx_storage.relabel - + let f12_over_r12_ints t = begin match t.ao_basis with @@ -224,9 +223,9 @@ let f12_over_r12_ints t = | _ -> not_implemented () end |> Four_idx_storage.relabel - + let cartesian t = t.cartesian - + let values t point = begin @@ -235,9 +234,59 @@ let values t point = | _ -> not_implemented () end |> Vector.relabel - + #+end_src +** TREXIO + +*** Read + + #+begin_src ocaml :tangle (eval mli) +(* +val of_trexio : Trexio.trexio_file -> t +*) + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +(* +let of_trexio f = + + match Trexio.read_basis_type f with + | "G" + | "Gaussian" + | "gaussian" -> + let basis = + Gaussian.Basis.of_trexio f + in + let ao_basis = + Basis_poly.Gaussian (Basis_gaussian.make ~basis ?operators ~cartesian nuclei ) + in + { ao_basis ; cartesian } + | _ -> not_implemented () + +*) + #+end_src + +*** Write + + #+begin_src ocaml :tangle (eval mli) +(* +val to_trexio : Trexio.trexio_file -> t -> unit +*) + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +(* +let to_trexio f t = + + match t.ao_basis with + | Basis_poly.Gaussian b -> Gaussian.Basis.to_trexio f (Basis_gaussian.basis b)) + | _ -> not_implemented () + +*) + + #+end_src + ** Printers #+begin_src ocaml :tangle (eval mli) diff --git a/ao/basis_gaussian.org b/ao/basis_gaussian.org index 7e2f969..26bd7bc 100644 --- a/ao/basis_gaussian.org +++ b/ao/basis_gaussian.org @@ -7,14 +7,14 @@ (setq ml (concat lib name ".ml")) (setq test-ml (concat testdir name ".ml")) (org-babel-tangle) -#+end_src +#+end_src * Gaussian basis :PROPERTIES: :header-args: :noweb yes :comments both :END: - Data structure for Gaussian Atomic Orbitals: + Data structure for Gaussian Atomic Orbitals: \[ \chi_i(\mathbf{r}) = P_i(\mathbf{r}) \sum_k c_k \exp\left( -\alpha_k (\mathbf{r-R_A})^2 \right) @@ -34,7 +34,7 @@ open Linear_algebra open Gaussian_integrals open Operators -type t +type t #+end_src #+begin_src ocaml :tangle (eval ml) :exports none @@ -42,7 +42,7 @@ open Linear_algebra open Gaussian open Gaussian_integrals open Operators - + module Basis = Gaussian.Basis type t = @@ -108,7 +108,7 @@ let overlap t = Lazy.force t.overlap let size t = Matrix.dim1 (Lazy.force t.overlap) #+end_src -** Computation +** Computation #+begin_src ocaml :tangle (eval mli) val values : t -> @@ -145,7 +145,7 @@ val make : basis:Gaussian.Basis.t -> Nuclei.t -> t #+end_src - + Creates the data structure for atomic orbitals from a Gaussian basis and the molecular geometry ~Nuclei.t~. @@ -161,28 +161,28 @@ val b : Ao.Basis_gaussian.t = Gaussian Basis, spherical, 15 AOs #+begin_src ocaml :tangle (eval ml) :exports none let make ~basis ?(operators=[]) ?(cartesian=false) nuclei = - + let overlap = lazy ( Overlap.of_basis basis ) in - + let ortho = lazy ( Lazy.force overlap - |> Orthonormalization.make ~cartesian ~basis + |> Orthonormalization.make ~cartesian ~basis ) in - + let eN_ints = lazy ( Electron_nucleus.of_basis_nuclei ~basis nuclei ) in - + let kin_ints = lazy ( Kinetic.of_basis basis ) in - + let ee_ints = lazy ( Eri.of_basis basis ) in - + let rec get_f12 = function | (Operator.F12 _ as f) :: _ -> f | [] -> failwith "Missing F12 operator" @@ -198,11 +198,11 @@ let make ~basis ?(operators=[]) ?(cartesian=false) nuclei = let ee_lr_ints = lazy ( Eri_long_range.of_basis basis~operator:(get_rs operators) ) in - + let f12_ints = lazy ( F12.of_basis basis ~operator:(get_f12 operators) ) in - + let f12_over_r12_ints = lazy ( Screened_eri.of_basis basis ~operator:(get_f12 operators) ) in @@ -227,6 +227,6 @@ let pp ppf t = let cart = if t.cartesian then "cartesian" else "spherical" in let nao = size t in Format.fprintf ppf "@[@[Gaussian Basis@], @[%s@], @[%d AOs@]@]" - cart nao + cart nao #+end_src diff --git a/ao/lib/basis.ml b/ao/lib/basis.ml index 3fddc66..8b4ad48 100644 --- a/ao/lib/basis.ml +++ b/ao/lib/basis.ml @@ -26,9 +26,12 @@ open Common * let b = Ao.Basis.of_nuclei_and_basis_filename ~nuclei filename;; * val b : Ao.Basis.t = Gaussian Basis, spherical, 15 AOs * #+end_example *) - + (* [[file:~/QCaml/ao/basis.org::*Conversions][Conversions:2]] *) +let not_implemented () = + Util.not_implemented "Only Gaussian is implemented" + let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false) ~nuclei filename = match kind with @@ -36,11 +39,11 @@ let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false) let basis = Gaussian.Basis.of_nuclei_and_basis_filename ~nuclei filename in - let ao_basis = + let ao_basis = Basis_poly.Gaussian (Basis_gaussian.make ~basis ?operators ~cartesian nuclei ) in { ao_basis ; cartesian } - | _ -> failwith "of_nuclei_and_basis_filename needs to be called with `Gaussian" + | _ -> not_implemented () (* Conversions:2 ends here *) @@ -64,12 +67,9 @@ let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false) (* [[file:~/QCaml/ao/basis.org::*Access][Access:2]] *) -let not_implemented () = - Util.not_implemented "Only Gaussian is implemented" - let ao_basis t = t.ao_basis - -let size t = + +let size t = match t.ao_basis with | Basis_poly.Gaussian b -> Basis_gaussian.size b | _ -> not_implemented () @@ -81,7 +81,7 @@ let overlap t = | _ -> not_implemented () end |> Matrix.relabel - + let multipole t = begin match t.ao_basis with @@ -92,7 +92,7 @@ let multipole t = |> Matrix.relabel | _ -> not_implemented () end - + let ortho t = begin match t.ao_basis with @@ -100,7 +100,7 @@ let ortho t = | _ -> not_implemented () end |> Matrix.relabel - + let eN_ints t = begin match t.ao_basis with @@ -108,7 +108,7 @@ let eN_ints t = | _ -> not_implemented () end |> Matrix.relabel - + let kin_ints t = begin match t.ao_basis with @@ -116,7 +116,7 @@ let kin_ints t = | _ -> not_implemented () end |> Matrix.relabel - + let ee_ints t = begin match t.ao_basis with @@ -124,7 +124,7 @@ let ee_ints t = | _ -> not_implemented () end |> Four_idx_storage.relabel - + let ee_lr_ints t = begin match t.ao_basis with @@ -132,7 +132,7 @@ let ee_lr_ints t = | _ -> not_implemented () end |> Four_idx_storage.relabel - + let f12_ints t = begin match t.ao_basis with @@ -140,7 +140,7 @@ let f12_ints t = | _ -> not_implemented () end |> Four_idx_storage.relabel - + let f12_over_r12_ints t = begin match t.ao_basis with @@ -148,9 +148,9 @@ let f12_over_r12_ints t = | _ -> not_implemented () end |> Four_idx_storage.relabel - + let cartesian t = t.cartesian - + let values t point = begin @@ -161,6 +161,37 @@ let values t point = |> Vector.relabel (* Access:2 ends here *) +(* [[file:~/QCaml/ao/basis.org::*Read][Read:2]] *) +(* +let of_trexio f = + + match Trexio.read_basis_type f with + | "G" + | "Gaussian" + | "gaussian" -> + let basis = + Gaussian.Basis.of_trexio f + in + let ao_basis = + Basis_poly.Gaussian (Basis_gaussian.make ~basis ?operators ~cartesian nuclei ) + in + { ao_basis ; cartesian } + | _ -> not_implemented () + +*) +(* Read:2 ends here *) + +(* [[file:~/QCaml/ao/basis.org::*Write][Write:2]] *) +(* +let to_trexio f t = + + match t.ao_basis with + | Basis_poly.Gaussian b -> Gaussian.Basis.to_trexio f (Basis_gaussian.basis b)) + | _ -> not_implemented () + +*) +(* Write:2 ends here *) + (* [[file:~/QCaml/ao/basis.org::*Printers][Printers:2]] *) let pp ppf t = begin diff --git a/ao/lib/basis.mli b/ao/lib/basis.mli index 00a9815..ae2c37c 100644 --- a/ao/lib/basis.mli +++ b/ao/lib/basis.mli @@ -4,7 +4,7 @@ (* [[file:~/QCaml/ao/basis.org::*Types][Types:1]] *) type t -type ao = Ao_dim.t +type ao = Ao_dim.t open Common open Particles @@ -44,6 +44,24 @@ val cartesian : t -> bool val values : t -> Coordinate.t -> ao Vector.t (* Access:1 ends here *) +(* Read *) + + +(* [[file:~/QCaml/ao/basis.org::*Read][Read:1]] *) +(* +val of_trexio : Trexio.trexio_file -> t +*) +(* Read:1 ends here *) + +(* Write *) + + +(* [[file:~/QCaml/ao/basis.org::*Write][Write:1]] *) +(* +val to_trexio : Trexio.trexio_file -> t -> unit +*) +(* Write:1 ends here *) + (* Printers *) diff --git a/ao/lib/basis_gaussian.ml b/ao/lib/basis_gaussian.ml index ba9c551..8095509 100644 --- a/ao/lib/basis_gaussian.ml +++ b/ao/lib/basis_gaussian.ml @@ -3,7 +3,7 @@ open Linear_algebra open Gaussian open Gaussian_integrals open Operators - + module Basis = Gaussian.Basis type t = @@ -79,7 +79,7 @@ let values t point = (* Computation:2 ends here *) - + (* Creates the data structure for atomic orbitals from a Gaussian basis and the * molecular geometry ~Nuclei.t~. * @@ -96,28 +96,28 @@ let values t point = (* [[file:~/QCaml/ao/basis_gaussian.org::*Creation][Creation:2]] *) let make ~basis ?(operators=[]) ?(cartesian=false) nuclei = - + let overlap = lazy ( Overlap.of_basis basis ) in - + let ortho = lazy ( Lazy.force overlap - |> Orthonormalization.make ~cartesian ~basis + |> Orthonormalization.make ~cartesian ~basis ) in - + let eN_ints = lazy ( Electron_nucleus.of_basis_nuclei ~basis nuclei ) in - + let kin_ints = lazy ( Kinetic.of_basis basis ) in - + let ee_ints = lazy ( Eri.of_basis basis ) in - + let rec get_f12 = function | (Operator.F12 _ as f) :: _ -> f | [] -> failwith "Missing F12 operator" @@ -133,11 +133,11 @@ let make ~basis ?(operators=[]) ?(cartesian=false) nuclei = let ee_lr_ints = lazy ( Eri_long_range.of_basis basis~operator:(get_rs operators) ) in - + let f12_ints = lazy ( F12.of_basis basis ~operator:(get_f12 operators) ) in - + let f12_over_r12_ints = lazy ( Screened_eri.of_basis basis ~operator:(get_f12 operators) ) in diff --git a/ao/lib/basis_gaussian.mli b/ao/lib/basis_gaussian.mli index 40ca913..7c177bd 100644 --- a/ao/lib/basis_gaussian.mli +++ b/ao/lib/basis_gaussian.mli @@ -31,7 +31,7 @@ val overlap : t -> Overlap.t val size : t -> int (* Access:1 ends here *) -(* Computation *) +(* Computation *) (* [[file:~/QCaml/ao/basis_gaussian.org::*Computation][Computation:1]] *) diff --git a/ao/lib/basis_poly.ml b/ao/lib/basis_poly.ml index 6085725..6b4b157 100644 --- a/ao/lib/basis_poly.ml +++ b/ao/lib/basis_poly.ml @@ -1,4 +1,4 @@ -(* [[file:~/QCaml/ao/basis.org::*Polymorphic%20types][Polymorphic types:2]] *) +(* [[file:~/QCaml/ao/basis.org::*Polymorphic types][Polymorphic types:2]] *) type t = | Unknown | Gaussian of Basis_gaussian.t diff --git a/ao/test/dune b/ao/test/dune index cec9ed7..f983879 100644 --- a/ao/test/dune +++ b/ao/test/dune @@ -4,6 +4,7 @@ (synopsis "Test for ao library") (libraries alcotest + trexio qcaml.ao ) ) diff --git a/docs/ao.html b/docs/ao.html index 3cadee8..aebe4ac 100644 --- a/docs/ao.html +++ b/docs/ao.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Atomic Orbitals @@ -272,31 +272,37 @@ org_html_manager.setup(); // activate after the parameters are set

Table of Contents

-
-

1 Summary

+
+

1 Summary

This modules contains the data structures used to characterize the @@ -305,11 +311,11 @@ atomic basis set.

-
-

2 Gaussian basis

+
+

2 Gaussian basis

-Data structure for Gaussian Atomic Orbitals: +Data structure for Gaussian Atomic Orbitals:

@@ -324,27 +330,27 @@ nucleus \(A\).

-
-

2.1 Type

+
+

2.1 Type

-Gaussian_basis.t +Gaussian_basis.t

-
open Common
+
open Common
 open Particles
 open Linear_algebra
 open Gaussian_integrals
 open Operators
 
-type t 
+type t
-
-

2.2 Access

+
+

2.2 Access

val basis             : t -> Gaussian.Basis.t
@@ -435,8 +441,8 @@ nucleus \(A\).
-
-

2.3 Computation

+
+

2.3 Computation

val values : t ->
@@ -463,8 +469,8 @@ nucleus \(A\).
-
-

2.4 Creation

+
+

2.4 Creation

val make : basis:Gaussian.Basis.t ->
@@ -491,15 +497,15 @@ Defaults:
 

 Example:
 

-
+
 let b = Ao.Basis_gaussian.make ~basis nuclei ;;
 val b : Ao.Basis_gaussian.t = Gaussian Basis, spherical, 15 AOs
-
-

2.5 Printers

+
+

2.5 Printers

val pp : Format.formatter -> t -> unit
@@ -509,22 +515,22 @@ val b : Ao.Basis_gaussian.t = Gaussian Basis, spherical, 15 AOs
-
-

3 Basis

+
+

3 Basis

Data structure for Atomic Orbitals.

-
-

3.1 Polymorphic types

+
+

3.1 Polymorphic types

-Basis.t +Basis.t

-
type t =
+
type t =
   | Unknown
   | Gaussian of Basis_gaussian.t
 

@@ -532,15 +538,15 @@ Data structure for Atomic Orbitals.
-
-

3.2 Types

+
+

3.2 Types

-Basis.t +Basis.t 

type t
-type ao = Ao_dim.t 
+type ao = Ao_dim.t
 
 open Common
 open Particles
@@ -551,8 +557,8 @@ Data structure for Atomic Orbitals.
-
-

3.3 Conversions

+
+

3.3 Conversions

val of_nuclei_and_basis_filename :
@@ -593,15 +599,15 @@ Defaults:
 

 Example:
 

-
+
 let b = Ao.Basis.of_nuclei_and_basis_filename ~nuclei filename;;
 val b : Ao.Basis.t = Gaussian Basis, spherical, 15 AOs
-
-

3.4 Access

+
+

3.4 Access

val size              : t -> int
@@ -699,9 +705,38 @@ val b : Ao.Basis.t = Gaussian Basis, spherical, 15 AOs
-
-

3.5 Printers

+
+

3.5 TREXIO

+
+
+

3.5.1 Read

+
+
+
(*
+val of_trexio : Trexio.trexio_file -> t
+*)
+
+
+
+
+ +
+

3.5.2 Write

+
+
+
(*
+val to_trexio : Trexio.trexio_file -> t -> unit
+*)
+
+
+
+
+
+ +
+

3.6 Printers

+
val pp : Format.formatter -> t -> unit
@@ -712,7 +747,7 @@ val b : Ao.Basis.t = Gaussian Basis, spherical, 15 AOs

Author: Anthony Scemama

-

Created: 2021-01-28 Thu 00:21

+

Created: 2023-04-24 Mon 18:06

Validate

diff --git a/docs/gaussian.html b/docs/gaussian.html index 711cf85..350d18b 100644 --- a/docs/gaussian.html +++ b/docs/gaussian.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Gaussian @@ -272,41 +272,41 @@ org_html_manager.setup(); // activate after the parameters are set

Table of Contents

-
-

1 Summary

+
+

1 Summary

-
-

2 Atomic shell

+
+

2 Atomic shell

Set of contracted Gaussians differing only by the powers of \(x\), \(y\) and \(z\), with a @@ -339,8 +339,8 @@ particular powers of \(x,y,z\) (PrimitiveShell.norm_coef_scale) -

-

2.1 Type

+
+

2.1 Type

type t
@@ -351,8 +351,8 @@ particular powers of \(x,y,z\) (PrimitiveShell.norm_coef_scale)
-
-

2.2 Access

+
+

2.2 Access

val ang_mom           : t -> Angular_momentum.t
@@ -429,14 +429,14 @@ particular powers of \(x,y,z\) (PrimitiveShell.norm_coef_scale)
 
 
-
+
-
-

2.3 Creation

+
+

2.3 Creation

val make : ?index:int -> Contracted_shell.t array -> t 
@@ -468,8 +468,8 @@ particular powers of \(x,y,z\) (PrimitiveShell.norm_coef_scale)
-
-

2.4 Printers

+
+

2.4 Printers

val pp : Format.formatter -> t -> unit
@@ -479,8 +479,8 @@ particular powers of \(x,y,z\) (PrimitiveShell.norm_coef_scale)
-
-

3 Atomic shell pair couple

+
+

3 Atomic shell pair couple

An atomic shell pair couple is the cartesian product between two sets of functions, one @@ -496,8 +496,8 @@ acting on different electrons, since they will be coupled by a two-electron oper

-
-

3.1 Type

+
+

3.1 Type

type t
@@ -508,8 +508,8 @@ acting on different electrons, since they will be coupled by a two-electron oper
-
-

3.2 Access

+
+

3.2 Access

val ang_mom                       : t -> Angular_momentum.t
@@ -594,8 +594,8 @@ acting on different electrons, since they will be coupled by a two-electron oper
-
-

3.3 Creation

+
+

3.3 Creation

val make : ?cutoff:float -> Atomic_shell_pair.t -> Atomic_shell_pair.t -> t option
@@ -621,14 +621,14 @@ Default cutoff is \(\epsilon\).
 
 
 
-
+
-
-

3.4 Printers

+
+

3.4 Printers

val pp : Format.formatter -> t -> unit
@@ -638,8 +638,8 @@ Default cutoff is \(\epsilon\).
-
-

4 Atomic shell pair

+
+

4 Atomic shell pair

Data structure to represent pairs of atomic shells. The products of @@ -651,8 +651,8 @@ An atomic shell pair is an array of pairs of contracted shells.

-
-

4.1 Type

+
+

4.1 Type

type t
@@ -663,8 +663,8 @@ An atomic shell pair is an array of pairs of contracted shells.
-
-

4.2 Access

+
+

4.2 Access

val atomic_shell_a         : t -> Atomic_shell.t
@@ -731,8 +731,8 @@ An atomic shell pair is an array of pairs of contracted shells.
-
-

4.3 Creation

+
+

4.3 Creation

val make : ?cutoff:float -> Atomic_shell.t -> Atomic_shell.t -> t option
@@ -765,8 +765,8 @@ If an atomic shell pair is not significant, sets the value to None.
-
-

4.4 Printers

+
+

4.4 Printers

val pp : Format.formatter -> t -> unit
@@ -778,7 +778,7 @@ If an atomic shell pair is not significant, sets the value to None.

Author: Anthony Scemama

-

Created: 2023-04-20 Thu 12:08

+

Created: 2023-04-24 Mon 19:28

Validate

diff --git a/docs/particles.html b/docs/particles.html index f4e96e4..0edc3db 100644 --- a/docs/particles.html +++ b/docs/particles.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Common @@ -272,75 +272,81 @@ org_html_manager.setup(); // activate after the parameters are set

Table of Contents

-
-

1 Summary

+
+

1 Summary

-
-

2 Electrons

+
+

2 Electrons

Data structure which contains the number of α and β electrons.

-
-

2.1 Type

+
+

2.1 Type

-
type t 
+
type t
-
-

2.2 Creation

+
+

2.2 Creation

open Common
@@ -372,15 +378,15 @@ Data structure which contains the number of α and β electrons.
 
 
 of_atoms
-Creates the data relative to electrons for a molecular system described by Nuclei.t for a given total charge and spin multiplicity.
+Creates the data relative to electrons for a molecular system described by Nuclei.t for a given total charge and spin multiplicity.
-
-

2.3 Access

+
+

2.3 Access

val charge       : t -> Charge.t
@@ -429,8 +435,8 @@ Data structure which contains the number of α and β electrons.
-
-

2.4 Printers

+
+

2.4 Printers

val pp : Format.formatter -> t -> unit
@@ -439,24 +445,24 @@ Data structure which contains the number of α and β electrons.
-
-

2.5 Tests

+
+

2.5 Tests

-
-

3 Element

+
+

3 Element

Chemical elements.

-
-

3.1 Type

+
+

3.1 Type

-
type t =
+
type t =
   |X
   |H                                                 |He
   |Li|Be                              |B |C |N |O |F |Ne
@@ -473,8 +479,8 @@ Chemical elements.
-
-

3.2 Conversion

+
+

3.2 Conversion

val of_string      : string -> t
@@ -535,7 +541,7 @@ Chemical elements.
 
 
 
-
+
 Element.of_string "Fe" ;;
 - : Element.t = Particles.Element.Fe
 
@@ -554,8 +560,8 @@ Element.(to_string Fe);;
-
-

3.3 Database information

+
+

3.3 Database information

val covalent_radius : t -> Non_negative_float.t
@@ -604,8 +610,8 @@ Element.(to_string Fe);;
-
-

3.4 Printers

+
+

3.4 Printers

val pp      : Format.formatter -> t -> unit
@@ -616,34 +622,34 @@ Element.(to_string Fe);;
-
-

4 Atomic mass

+
+

4 Atomic mass

Atomic mass, a non-negative float.

-
include module type of Common.Non_negative_float
+
include module type of Common.Non_negative_float
-
-

5 Nuclei

+
+

5 Nuclei

-
-

5.1 Type

+
+

5.1 Type

-Nuclei.t +Nuclei.t

-
open Common
+
open Common
 
 type t = (Element.t * Coordinate.t) array
 

@@ -651,12 +657,12 @@ Atomic mass, a non-negative float.
-
-

5.2 xyz file lexer/parser

+
+

5.2 xyz file lexer/parser

-
-

5.2.1 Lexer

+
+

5.2.1 Lexer

nuclei_lexer.mll contains the description of the lexemes used in @@ -709,11 +715,11 @@ rule read_all = parse

-
-

5.2.2 Parser

+
+

5.2.2 Parser

-xyz_parser.mly parses nuclear coordinates in xyz format. +xyz_parser.mly parses nuclear coordinates in xyz format. 

%{
@@ -829,8 +835,8 @@ an xyz_file data structure.
-
-

5.3 Conversion

+
+

5.3 Conversion

val of_xyz_string : string -> t
@@ -843,9 +849,6 @@ an xyz_file data structure.
 val to_string      : t -> string
 
 val of_filename : string -> t
-
-val of_trexio : Trexio.trexio_file -> t
-val to_trexio : Trexio.trexio_file -> t -> unit
@@ -887,23 +890,13 @@ an xyz_file data structure. of_filename Detects the type of file (xyz, z-matrix) and reads the file - - -of_trexio -Reads the geometry from a TREXIO file - - - -to_trexio -Writes the geometry in a TREXIO file -
-
-

5.4 Query

+
+

5.4 Query

val formula    : t -> string
@@ -952,9 +945,34 @@ an xyz_file data structure.
-
-

5.5 Printers

+
+

5.5 TREXIO

+
+
+

5.5.1 Read

+
+
+
val of_trexio : Trexio.trexio_file -> t
+
+
+
+
+ +
+

5.5.2 Write

+
+
+
val to_trexio : Trexio.trexio_file -> t -> unit
+
+
+
+
+
+ +
+

5.6 Printers

+
val pp : Format.formatter -> t -> unit
@@ -962,31 +980,31 @@ an xyz_file data structure.
-
-

5.6 Tests

+
+

5.7 Tests

-
-

6 Z-matrix

+
+

6 Z-matrix

Z-matrix representation of nuclear coordinates.

-
-

6.1 Type

+
+

6.1 Type

-
type t 
+
type t
-
-

6.2 Conversion

+
+

6.2 Conversion

val of_string      : string -> t
@@ -1021,7 +1039,7 @@ Z-matrix representation of nuclear coordinates.
 
 
 
-
+
 let zmt = Zmatrix.of_string "
  n
  n    1 nn
@@ -1066,8 +1084,8 @@ H -0.568803 -0.793910 1.726048"
-
-

6.3 Printers

+
+

6.3 Printers

val pp : Format.formatter -> t -> unit
@@ -1079,7 +1097,7 @@ H -0.568803 -0.793910 1.726048"

Author: Anthony Scemama

-

Created: 2023-04-24 Mon 14:13

+

Created: 2023-04-24 Mon 18:07

Validate

diff --git a/docs/top.html b/docs/top.html index 54a8721..bdfa85a 100644 --- a/docs/top.html +++ b/docs/top.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Top-level @@ -250,18 +250,18 @@ org_html_manager.setup(); // activate after the parameters are set

Table of Contents

-
-

1 Summary

+
+

1 Summary

Author: Anthony Scemama

-

Created: 2023-04-24 Mon 14:13

+

Created: 2023-04-24 Mon 19:28

Validate

diff --git a/gaussian/lib/basis.ml b/gaussian/lib/basis.ml index 6667248..9df7650 100644 --- a/gaussian/lib/basis.ml +++ b/gaussian/lib/basis.ml @@ -1,4 +1,4 @@ -type t = +type t = { size : int; contracted_shells : Contracted_shell.t array ; @@ -16,11 +16,11 @@ let general_basis t = t.general_basis (** Returns an array of the basis set per atom *) let of_nuclei_and_general_basis nucl bas = let index_ = ref 0 in - let contracted_shells = + let contracted_shells = Array.map (fun (e, center) -> List.assoc e bas |> Array.map (fun (ang_mom, shell) -> - let lc = + let lc = Array.map (fun Gb.{exponent ; coefficient} -> coefficient, Ps.make ang_mom center exponent) shell in @@ -38,11 +38,11 @@ let of_nuclei_and_general_basis nucl bas = |> Array.to_list |> List.sort_uniq compare in - let csl = + let csl = Array.to_list contracted_shells in List.map (fun (center, ang_mom) -> - let a = + let a = List.filter (fun x -> Cs.center x = center && Cs.ang_mom x = ang_mom) csl |> Array.of_list in @@ -65,7 +65,7 @@ let to_string b = let b = atomic_shells b in let line =" ----------------------------------------------------------------------- -" in +" in " Atomic Basis set ---------------- @@ -85,26 +85,105 @@ let to_string b = -let of_nuclei_and_basis_filename ~nuclei filename = - let general_basis = +let of_nuclei_and_basis_filename ~nuclei filename = + let general_basis = General_basis.read filename in of_nuclei_and_general_basis nuclei general_basis -let of_nuclei_and_basis_string ~nuclei str = - let general_basis = +let of_nuclei_and_basis_string ~nuclei str = + let general_basis = General_basis.of_string str in of_nuclei_and_general_basis nuclei general_basis -let of_nuclei_and_basis_filenames ~nuclei filenames = - let general_basis = +let of_nuclei_and_basis_filenames ~nuclei filenames = + let general_basis = General_basis.read_many filenames in of_nuclei_and_general_basis nuclei general_basis +let of_trexio f = + let nuclei = Particles.Nuclei.of_trexio f in + (* + let nucl_num = Array.length nuclei in + let prim_num = Trexio.read_basis_prim_num f in + let shell_num = Trexio.read_basis_shell_num f in + let shell_factor = Trexio.read_basis_shell_factor f in + let prim_factor = Trexio.read_basis_prim_factor f in + *) + let shell_ang_mom = Trexio.read_basis_shell_ang_mom f + |> Array.map Common.Angular_momentum.of_int + in + let shell_index = Trexio.read_basis_shell_index f + |> Array.to_list + in + let nucleus_index = Trexio.read_basis_nucleus_index f + |> Array.to_list + in + let exponent = Trexio.read_basis_exponent f in + let coefficient = Trexio.read_basis_coefficient f in + + let filter_indices ref_index full_array = + List.mapi (fun j idx -> + if idx = ref_index then Some j + else None) full_array + |> List.filter (function + | None -> false + | _ -> true) + |> List.map (function + | Some x -> x + | None -> assert false) + in + + let index_ = ref 0 in + let contracted_shells = + Array.mapi (fun i (_,center) -> + let shell_indices = filter_indices i nucleus_index in + List.map (fun ishell -> + let prim_indices = filter_indices ishell shell_index in + let lc = List.map (fun k -> + coefficient.(k), + Ps.make shell_ang_mom.(ishell) center exponent.(k)) prim_indices + |> Array.of_list + in + let result = Cs.make ~index:!index_ lc in + index_ := !index_ + Cs.size_of_shell result; + result + ) shell_indices + |> Array.of_list + ) nuclei + |> Array.to_list + |> Array.concat + in + let atomic_shells = lazy( + let uniq_center_angmom = + Array.map (fun x -> Cs.center x, Cs.ang_mom x) contracted_shells + |> Array.to_list + |> List.sort_uniq compare + in + let csl = + Array.to_list contracted_shells + in + List.map (fun (center, ang_mom) -> + let a = + List.filter (fun x -> Cs.center x = center && Cs.ang_mom x = ang_mom) csl + |> Array.of_list + in + As.make ~index:(Cs.index a.(0)) a + ) uniq_center_angmom + |> List.sort (fun x y -> compare (As.index x) (As.index y)) + |> Array.of_list + ) in + { contracted_shells ; atomic_shells ; size = !index_; + general_basis = []; + } + +let to_trexio _ _ = + () + let pp ppf t = Format.fprintf ppf "@[%s@]" (to_string t) diff --git a/gaussian/lib/basis.mli b/gaussian/lib/basis.mli index f258f8d..469d203 100644 --- a/gaussian/lib/basis.mli +++ b/gaussian/lib/basis.mli @@ -42,6 +42,13 @@ val contracted_shells : t -> Contracted_shell.t array val general_basis : t -> General_basis.t (** Returns the [!GeneralBasis] that was used to build the current basis. *) + +(** {TREXIO} *) + +val of_trexio : Trexio.trexio_file -> t +val to_trexio : Trexio.trexio_file -> t -> unit + + (** {2 Printers} *) val to_string : t -> string diff --git a/particles/lib/nuclei.ml b/particles/lib/nuclei.ml index b62f966..508cfcf 100644 --- a/particles/lib/nuclei.ml +++ b/particles/lib/nuclei.ml @@ -12,20 +12,18 @@ open Xyz_ast * | ~of_zmt_string~ | Create from a string, in z-matrix format | * | ~of_zmt_file~ | Create from a file, in z-matrix format | * | ~to_string~ | Transform to a string, for printing | - * | ~of_filename~ | Detects the type of file (xyz, z-matrix) and reads the file | - * | ~of_trexio~ | Reads the geometry from a TREXIO file | - * | ~to_trexio~ | Writes the geometry in a TREXIO file | *) + * | ~of_filename~ | Detects the type of file (xyz, z-matrix) and reads the file | *) (* [[file:~/QCaml/particles/nuclei.org::*Conversion][Conversion:2]] *) let of_xyz_lexbuf lexbuf = - let data = + let data = Xyz_parser.input Nuclei_lexer.read_all lexbuf in let len = List.length data.nuclei in if len <> data.number_of_atoms then - Printf.sprintf "Error: expected %d atoms but %d read" + Printf.sprintf "Error: expected %d atoms but %d read" data.number_of_atoms len |> failwith; @@ -35,7 +33,7 @@ let of_xyz_lexbuf lexbuf = |> Array.of_list -let of_xyz_string input_string = +let of_xyz_string input_string = Lexing.from_string input_string |> of_xyz_lexbuf @@ -45,7 +43,7 @@ let of_xyz_file filename = let lexbuf = Lexing.from_channel ic in - let result = + let result = of_xyz_lexbuf lexbuf in close_in ic; @@ -90,7 +88,7 @@ let to_string atoms = bohr_to_angstrom coord in Printf.sprintf " %5d %5d %5s %12.6f %12.6f %12.6f" - (i+1) (Element.to_int e) (Element.to_string e) + (i+1) (Element.to_int e) (Element.to_string e) coord.x coord.y coord.z ) atoms |> Array.to_list @@ -114,22 +112,9 @@ let to_xyz_string t = in Printf.sprintf " %5s %12.6f %12.6f %12.6f" (Element.to_string e) coord.x coord.y coord.z - ) t + ) t |> Array.to_list ) |> String.concat "\n" - - -let of_trexio f = - let num = Trexio.read_nucleus_num f in - let charge = Trexio.read_nucleus_charge f - |> Array.map Charge.of_float in - let coord = Trexio.read_nucleus_coord f in - Array.init num (fun i -> - let coord = Coordinate.{ x = coord.(3*i) ; - y = coord.(3*i+1) ; - z = coord.(3*i+2) } in - (Element.of_charge charge.(i), Coordinate.make coord) - ) (* Conversion:2 ends here *) @@ -160,14 +145,14 @@ let formula t = accu ^ key ^ "_{" ^ (string_of_int x) ^ "}") "" - + let repulsion nuclei = - let get_charge e = + let get_charge e = Element.to_charge e |> Charge.to_float in - Array.fold_left ( fun accu (e1, coord1) -> - accu +. + Array.fold_left ( fun accu (e1, coord1) -> + accu +. Array.fold_left (fun accu (e2, coord2) -> let r = Coordinate.(norm (coord1 |- coord2)) in if r > 0. then @@ -177,26 +162,42 @@ let repulsion nuclei = ) 0. nuclei -let charge nuclei = +let charge nuclei = Array.fold_left (fun accu (e, _) -> accu + Charge.to_int (Element.to_charge e) ) - 0 nuclei + 0 nuclei |> Charge.of_int -let small_core a = +let small_core a = Array.fold_left (fun accu (e,_) -> accu + (Element.small_core e)) 0 a -let large_core a = +let large_core a = Array.fold_left (fun accu (e,_) -> accu + (Element.large_core e)) 0 a +(* Query:2 ends here *) +(* [[file:~/QCaml/particles/nuclei.org::*Read][Read:2]] *) +let of_trexio f = + let num = Trexio.read_nucleus_num f in + let charge = Trexio.read_nucleus_charge f + |> Array.map Charge.of_float in + let coord = Trexio.read_nucleus_coord f in + Array.init num (fun i -> + let coord = Coordinate.{ x = coord.(3*i) ; + y = coord.(3*i+1) ; + z = coord.(3*i+2) } in + (Element.of_charge charge.(i), Coordinate.make coord) + ) +(* Read:2 ends here *) + +(* [[file:~/QCaml/particles/nuclei.org::*Write][Write:2]] *) let to_trexio f t = let num = Array.length t in Trexio.write_nucleus_num f num; - Array.map (fun (e, _) -> Element.to_charge e |> Charge.to_float) t + Array.map (fun (e, _) -> Element.to_charge e |> Charge.to_float) t |> Trexio.write_nucleus_charge f; - + Array.map (fun (e, _) -> Element.to_string e) t |> Trexio.write_nucleus_label f; @@ -206,10 +207,10 @@ let to_trexio f t = coord.(3*i+1) <- Coordinate.(get Y xyz) ; coord.(3*i+2) <- Coordinate.(get Z xyz) ) t; Trexio.write_nucleus_coord f coord; - + repulsion t |> Trexio.write_nucleus_repulsion f -(* Query:2 ends here *) +(* Write:2 ends here *) (* [[file:~/QCaml/particles/nuclei.org::*Printers][Printers:2]] *) let pp ppf t = diff --git a/particles/lib/nuclei.mli b/particles/lib/nuclei.mli index 3c2e9a0..775ce49 100644 --- a/particles/lib/nuclei.mli +++ b/particles/lib/nuclei.mli @@ -23,9 +23,6 @@ val of_zmt_file : string -> t val to_string : t -> string val of_filename : string -> t - -val of_trexio : Trexio.trexio_file -> t -val to_trexio : Trexio.trexio_file -> t -> unit (* Conversion:1 ends here *) (* Query *) @@ -39,6 +36,20 @@ val small_core : t -> int val large_core : t -> int (* Query:1 ends here *) +(* Read *) + + +(* [[file:~/QCaml/particles/nuclei.org::*Read][Read:1]] *) +val of_trexio : Trexio.trexio_file -> t +(* Read:1 ends here *) + +(* Write *) + + +(* [[file:~/QCaml/particles/nuclei.org::*Write][Write:1]] *) +val to_trexio : Trexio.trexio_file -> t -> unit +(* Write:1 ends here *) + (* Printers *) diff --git a/particles/lib/nuclei_lexer.mll b/particles/lib/nuclei_lexer.mll index 128a4f7..ac36dad 100644 --- a/particles/lib/nuclei_lexer.mll +++ b/particles/lib/nuclei_lexer.mll @@ -1,8 +1,8 @@ (* Lexer - * + * * =nuclei_lexer.mll= contains the description of the lexemes used in * an xyz file. *) - + (* [[file:~/QCaml/particles/nuclei.org::*Lexer][Lexer:1]] *) { diff --git a/particles/nuclei.org b/particles/nuclei.org index 6aaa19c..669413a 100644 --- a/particles/nuclei.org +++ b/particles/nuclei.org @@ -7,7 +7,7 @@ (setq ml (concat lib name ".ml")) (setq test-ml (concat testdir name ".ml")) (org-babel-tangle) -#+end_src +#+end_src * Nuclei :PROPERTIES: @@ -34,10 +34,10 @@ open Xyz_ast ** xyz file lexer/parser *** Lexer - + =nuclei_lexer.mll= contains the description of the lexemes used in an xyz file. - + #+begin_src ocaml :tangle lib/nuclei_lexer.mll :export none { open Xyz_parser @@ -79,11 +79,11 @@ rule read_all = parse done; *) } - #+end_src + #+end_src *** Parser - - =xyz_parser.mly= parses nuclear coordinates in xyz format. + + =xyz_parser.mly= parses nuclear coordinates in xyz format. #+begin_src ocaml :tangle lib/xyz_parser.mly :export none :comments none %{ open Common @@ -169,12 +169,12 @@ atoms_list: | atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_int $5 $7 $9 $3 :: $1 } | atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_int $5 $7 $9 $3 :: $1 } ; - #+end_src + #+end_src When an xyz file is read by =xyz_parser.mly=, it is converted into an ~xyz_file~ data structure. - #+begin_src ocaml :tangle lib/xyz_ast.mli + #+begin_src ocaml :tangle lib/xyz_ast.mli open Common type nucleus = @@ -204,9 +204,6 @@ val of_zmt_file : string -> t val to_string : t -> string val of_filename : string -> t - -val of_trexio : Trexio.trexio_file -> t -val to_trexio : Trexio.trexio_file -> t -> unit #+end_src | ~of_xyz_string~ | Create from a string, in xyz format | @@ -215,18 +212,16 @@ val to_trexio : Trexio.trexio_file -> t -> unit | ~of_zmt_file~ | Create from a file, in z-matrix format | | ~to_string~ | Transform to a string, for printing | | ~of_filename~ | Detects the type of file (xyz, z-matrix) and reads the file | - | ~of_trexio~ | Reads the geometry from a TREXIO file | - | ~to_trexio~ | Writes the geometry in a TREXIO file | #+begin_src ocaml :tangle (eval ml) :exports none let of_xyz_lexbuf lexbuf = - let data = + let data = Xyz_parser.input Nuclei_lexer.read_all lexbuf in let len = List.length data.nuclei in if len <> data.number_of_atoms then - Printf.sprintf "Error: expected %d atoms but %d read" + Printf.sprintf "Error: expected %d atoms but %d read" data.number_of_atoms len |> failwith; @@ -236,7 +231,7 @@ let of_xyz_lexbuf lexbuf = |> Array.of_list -let of_xyz_string input_string = +let of_xyz_string input_string = Lexing.from_string input_string |> of_xyz_lexbuf @@ -246,7 +241,7 @@ let of_xyz_file filename = let lexbuf = Lexing.from_channel ic in - let result = + let result = of_xyz_lexbuf lexbuf in close_in ic; @@ -291,7 +286,7 @@ let to_string atoms = bohr_to_angstrom coord in Printf.sprintf " %5d %5d %5s %12.6f %12.6f %12.6f" - (i+1) (Element.to_int e) (Element.to_string e) + (i+1) (Element.to_int e) (Element.to_string e) coord.x coord.y coord.z ) atoms |> Array.to_list @@ -315,48 +310,11 @@ let to_xyz_string t = in Printf.sprintf " %5s %12.6f %12.6f %12.6f" (Element.to_string e) coord.x coord.y coord.z - ) t + ) t |> Array.to_list ) |> String.concat "\n" -let of_trexio f = - let num = Trexio.read_nucleus_num f in - let charge = Trexio.read_nucleus_charge f - |> Array.map Charge.of_float in - let coord = Trexio.read_nucleus_coord f in - Array.init num (fun i -> - let coord = Coordinate.{ x = coord.(3*i) ; - y = coord.(3*i+1) ; - z = coord.(3*i+2) } in - (Element.of_charge charge.(i), Coordinate.make coord) - ) - - - #+end_src - - #+NAME: trexio_export - #+begin_src ocaml :tangle no :exports none -let to_trexio f t = - let num = Array.length t in - Trexio.write_nucleus_num f num; - - Array.map (fun (e, _) -> Element.to_charge e |> Charge.to_float) t - |> Trexio.write_nucleus_charge f; - - Array.map (fun (e, _) -> Element.to_string e) t - |> Trexio.write_nucleus_label f; - - let coord = Array.init (num*3) (fun _ -> 0.) in - Array.iteri (fun i (_, xyz) -> - coord.(3*i) <- Coordinate.(get X xyz) ; - coord.(3*i+1) <- Coordinate.(get Y xyz) ; - coord.(3*i+2) <- Coordinate.(get Z xyz) ) t; - Trexio.write_nucleus_coord f coord; - - repulsion t - |> Trexio.write_nucleus_repulsion f - #+end_src ** Query @@ -375,7 +333,7 @@ val large_core : t -> int | ~small_core~ | Number of core electrons in the small core model | | ~large_core~ | Number of core electrons in the large core model | - #+begin_src ocaml :tangle (eval ml) :exports none :noweb yes + #+begin_src ocaml :tangle (eval ml) :exports none let formula t = let dict = Hashtbl.create 67 in Array.iter (fun (e,_) -> @@ -394,14 +352,14 @@ let formula t = accu ^ key ^ "_{" ^ (string_of_int x) ^ "}") "" - + let repulsion nuclei = - let get_charge e = + let get_charge e = Element.to_charge e |> Charge.to_float in - Array.fold_left ( fun accu (e1, coord1) -> - accu +. + Array.fold_left ( fun accu (e1, coord1) -> + accu +. Array.fold_left (fun accu (e2, coord2) -> let r = Coordinate.(norm (coord1 |- coord2)) in if r > 0. then @@ -411,22 +369,72 @@ let repulsion nuclei = ) 0. nuclei -let charge nuclei = +let charge nuclei = Array.fold_left (fun accu (e, _) -> accu + Charge.to_int (Element.to_charge e) ) - 0 nuclei + 0 nuclei |> Charge.of_int -let small_core a = +let small_core a = Array.fold_left (fun accu (e,_) -> accu + (Element.small_core e)) 0 a -let large_core a = +let large_core a = Array.fold_left (fun accu (e,_) -> accu + (Element.large_core e)) 0 a -<> #+end_src +** TREXIO + +*** Read + + #+begin_src ocaml :tangle (eval mli) +val of_trexio : Trexio.trexio_file -> t + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +let of_trexio f = + let num = Trexio.read_nucleus_num f in + let charge = Trexio.read_nucleus_charge f + |> Array.map Charge.of_float in + let coord = Trexio.read_nucleus_coord f in + Array.init num (fun i -> + let coord = Coordinate.{ x = coord.(3*i) ; + y = coord.(3*i+1) ; + z = coord.(3*i+2) } in + (Element.of_charge charge.(i), Coordinate.make coord) + ) + #+end_src + +*** Write + + #+begin_src ocaml :tangle (eval mli) +val to_trexio : Trexio.trexio_file -> t -> unit + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +let to_trexio f t = + let num = Array.length t in + Trexio.write_nucleus_num f num; + + Array.map (fun (e, _) -> Element.to_charge e |> Charge.to_float) t + |> Trexio.write_nucleus_charge f; + + Array.map (fun (e, _) -> Element.to_string e) t + |> Trexio.write_nucleus_label f; + + let coord = Array.init (num*3) (fun _ -> 0.) in + Array.iteri (fun i (_, xyz) -> + coord.(3*i) <- Coordinate.(get X xyz) ; + coord.(3*i+1) <- Coordinate.(get Y xyz) ; + coord.(3*i+2) <- Coordinate.(get Z xyz) ) t; + Trexio.write_nucleus_coord f coord; + + repulsion t + |> Trexio.write_nucleus_repulsion f + + #+end_src + ** Printers #+begin_src ocaml :tangle (eval mli)