Added simulation

common/charge.org

@@ -65,6 +65,7 @@ val ( + ) : t -> t -> t
 val ( - ) : t -> t -> t
 val ( * ) : t -> float -> t
 val ( / ) : t -> float -> t 
+val is_null : t -> bool
    #+end_src
 
    #+begin_src ocaml :tangle (eval ml) :exports none
@@ -77,6 +78,8 @@ let ( + ) = gen_op ( +. )
 let ( - ) = gen_op ( -. )
 let ( * ) = gen_op ( *. )
 let ( / ) = gen_op ( /. )
+
+let is_null t = t == 0.
    #+end_src
 
 ** Printers
@@ -87,6 +90,6 @@ val pp : Format.formatter -> t -> unit
 
    #+begin_src ocaml :tangle (eval ml) :exports none
 let pp ppf x = 
-  Format.fprintf ppf "@[+%s@]" (to_string x)
+  Format.fprintf ppf "@[%s@]" (to_string x)
    #+end_src
 

common/lib/charge.ml

@@ -3,11 +3,11 @@
 (* This type should be used for all charges in the program (electrons, nuclei,...). *)
 
 
-(* [[file:~/QCaml/common/charge.org::*Type][Type:2]] *)
+(* [[file:../charge.org::*Type][Type:2]] *)
 type t = float
 (* Type:2 ends here *)
 
-(* [[file:~/QCaml/common/charge.org::*Conversions][Conversions:2]] *)
+(* [[file:../charge.org::*Conversions][Conversions:2]] *)
 external of_float : float -> t = "%identity"
 external to_float : t -> float = "%identity"
 
@@ -25,7 +25,7 @@ let to_string x =
     "0.0"
 (* Conversions:2 ends here *)
 
-(* [[file:~/QCaml/common/charge.org::*Simple%20operations][Simple operations:2]] *)
+(* [[file:../charge.org::*Simple operations][Simple operations:2]] *)
 let gen_op op =
   fun a b ->
   op (to_float a) (to_float b)
@@ -35,9 +35,11 @@ let ( + ) = gen_op ( +. )
 let ( - ) = gen_op ( -. )
 let ( * ) = gen_op ( *. )
 let ( / ) = gen_op ( /. )
+
+let is_null t = t == 0.
 (* Simple operations:2 ends here *)
 
-(* [[file:~/QCaml/common/charge.org::*Printers][Printers:2]] *)
+(* [[file:../charge.org::*Printers][Printers:2]] *)
 let pp ppf x = 
-  Format.fprintf ppf "@[+%s@]" (to_string x)
+  Format.fprintf ppf "@[%s@]" (to_string x)
 (* Printers:2 ends here *)

common/lib/charge.mli

@@ -1,14 +1,14 @@
 (* Type *)
 
 
-(* [[file:~/QCaml/common/charge.org::*Type][Type:1]] *)
+(* [[file:../charge.org::*Type][Type:1]] *)
 type t
 (* Type:1 ends here *)
 
 (* Conversions *)
 
 
-(* [[file:~/QCaml/common/charge.org::*Conversions][Conversions:1]] *)
+(* [[file:../charge.org::*Conversions][Conversions:1]] *)
 val of_float : float -> t
 val to_float : t -> float
 
@@ -22,16 +22,17 @@ val to_string: t -> string
 (* Simple operations *)
 
 
-(* [[file:~/QCaml/common/charge.org::*Simple%20operations][Simple operations:1]] *)
+(* [[file:../charge.org::*Simple operations][Simple operations:1]] *)
 val ( + ) : t -> t -> t
 val ( - ) : t -> t -> t
 val ( * ) : t -> float -> t
 val ( / ) : t -> float -> t 
+val is_null : t -> bool
 (* Simple operations:1 ends here *)
 
 (* Printers *)
 
 
-(* [[file:~/QCaml/common/charge.org::*Printers][Printers:1]] *)
+(* [[file:../charge.org::*Printers][Printers:1]] *)
 val pp : Format.formatter -> t -> unit
 (* Printers:1 ends here *)

particles/lib/nuclei.ml

@@ -119,13 +119,33 @@ let to_xyz_string t =
 
 
 
-(* | ~repulsion~  | Nuclear repulsion energy, in atomic units        |
+(* | ~formula~    | Returns the chemical formula                     |
+ * | ~repulsion~  | Nuclear repulsion energy, in atomic units        |
  * | ~charge~     | Sum of the charges of the nuclei                 |
  * | ~small_core~ | Number of core electrons in the small core model |
  * | ~large_core~ | Number of core electrons in the large core model | *)
 
 
 (* [[file:../nuclei.org::*Query][Query:2]] *)
+let formula t =
+  let dict = Hashtbl.create 67 in
+  Array.iter (fun (e,_) ->
+    let e = Element.to_string e in
+    let value =
+      try (Hashtbl.find dict e) + 1
+      with Not_found -> 1
+    in
+    Hashtbl.replace dict e value
+  ) t;
+  Hashtbl.to_seq_keys dict
+  |> List.of_seq
+  |> List.sort String.compare
+  |> List.fold_left (fun accu key ->
+    let x = Hashtbl.find dict key in
+    accu ^ key ^ "_{" ^ (string_of_int x) ^ "}") ""
+
+
+  
 let repulsion nuclei =
   let get_charge e = 
     Element.to_charge e

particles/lib/nuclei.mli

@@ -24,10 +24,11 @@ val to_string      : t -> string
 val of_filename : string -> t
 (* Conversion:1 ends here *)
 
-(* Query *)
+(* TODO Query *)
 
 
 (* [[file:../nuclei.org::*Query][Query:1]] *)
+val formula    : t -> string
 val repulsion  : t -> float
 val charge     : t -> Charge.t
 val small_core : t -> int

particles/nuclei.org

@@ -24,7 +24,9 @@ type t = (Element.t * Coordinate.t) array
    #+end_src
 
    #+begin_src ocaml :tangle (eval ml) :exports none
-<<types>>
+open Common
+
+type t = (Element.t * Coordinate.t) array
 open Xyz_ast
    #+end_src
 
@@ -312,21 +314,42 @@ let to_xyz_string t =
   |> String.concat "\n"
    #+end_src
 
-** Query
+** TODO Query
 
    #+begin_src ocaml :tangle (eval mli)
+val formula    : t -> string
 val repulsion  : t -> float
 val charge     : t -> Charge.t
 val small_core : t -> int
 val large_core : t -> int
    #+end_src
 
+   | ~formula~    | Returns the chemical formula                     |
    | ~repulsion~  | Nuclear repulsion energy, in atomic units        |
    | ~charge~     | Sum of the charges of the nuclei                 |
    | ~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
+let formula t =
+  let dict = Hashtbl.create 67 in
+  Array.iter (fun (e,_) ->
+    let e = Element.to_string e in
+    let value =
+      try (Hashtbl.find dict e) + 1
+      with Not_found -> 1
+    in
+    Hashtbl.replace dict e value
+  ) t;
+  Hashtbl.to_seq_keys dict
+  |> List.of_seq
+  |> List.sort String.compare
+  |> List.fold_left (fun accu key ->
+    let x = Hashtbl.find dict key in
+    accu ^ key ^ "_{" ^ (string_of_int x) ^ "}") ""
+
+
+  
 let repulsion nuclei =
   let get_charge e = 
     Element.to_charge e

qcaml/README.org

@@ -10,14 +10,15 @@
 : Main QCaml entry point
 
 
-* Dune files                                                         :noexport:
+* File generation                                                  :noexport:
   
-
-   Use [C-c C-c] on the code below to create the output for the dune files
+   Use [C-c C-c] on the code below. This will create the dune files
+   and the =qcaml.ml= file.
 
    #+header: :noweb yes
    #+header: :var name=(file-name-directory buffer-file-name) 
    #+header: :var dune="lib/dune"
+   #+header: :var mlfile="lib/qcaml.ml"
    #+header: :var dunetest="test/dune"
    #+header: :var dependencies=dependencies
    #+begin_src python :exports none :results output none
@@ -26,13 +27,15 @@ synopsis = """
 <<synopsis>>
 """
 
-dependencies = '\n '.join(map(lambda x: x[0], dependencies))
+excluded_modules = [ "top" ]
+
+dependencies = '\n '.join([x[0] for x in dependencies if x[0].replace("qcaml.","") not in excluded_modules])
 
 with open(dune,'w') as f:
     f.write(f"""
 (library
- (name {name})
- (public_name qcaml.{name})
+ (name qcaml)
+ (public_name qcaml)
  (synopsis {synopsis}  )
  (libraries
  {dependencies}
@@ -41,19 +44,29 @@ with open(dune,'w') as f:
  <<c-files>>
 )
 <<lex-yacc>>
-""".replace("qcaml.qcaml","qcaml"))
+""")
 
 with open(dunetest,'w') as f:
     f.write(f"""
 (library
- (name test_{name})
- (synopsis "Test for {name} library")
+ (name test_qcaml)
+ (synopsis "Test for qcaml library")
  (libraries
   alcotest
-  qcaml.{name}
+  qcaml
  )
 )
-""".replace("qcaml.qcaml","qcaml"))
+""")
+
+with open(mlfile,'w') as f:
+    s = dependencies.split()
+    s = [x.replace("qcaml.","") for x in s ]
+    s = [x.capitalize() for x in s]
+    s = [ f"module {x} = {x}\n" for x in s ]
+    s = ''.join(s)
+    f.write("(* Auto-generated by qcaml/README.org *)\n")
+    f.write(s)
+
    #+end_src
 
 
@@ -63,7 +76,6 @@ grep "public_name" ../*/lib/dune | grep -v "qcaml)" | cut -d ' ' -f 3 | tr ')' '
    #+end_src
 
    #+RESULTS: dependencies
-   | (libraries               |
    | qcaml.ao                 |
    | qcaml.common             |
    | qcaml.gaussian           |
@@ -74,4 +86,4 @@ grep "public_name" ../*/lib/dune | grep -v "qcaml)" | cut -d ' ' -f 3 | tr ')' '
    | qcaml.particles          |
    | qcaml.perturbation       |
    | qcaml.simulation         |
-   | )                        |
+   | qcaml.top                |

qcaml/lib/qcaml.ml

@@ -1,3 +1,4 @@
+(* Auto-generated by qcaml/README.org *)
 module Ao = Ao
 module Common = Common
 module Gaussian = Gaussian

simulation/lib/simulation.ml

@@ -1,7 +1,10 @@
+(* [[file:../simulation.org::*Simulation][Simulation:2]] *)
 open Common
 open Particles
 open Operators
+(* Simulation:2 ends here *)
 
+(* [[file:../simulation.org::*Type][Type:2]] *)
 type t = {
   charge            : Charge.t;
   electrons         : Electrons.t;
@@ -9,15 +12,38 @@ type t = {
   ao_basis          : Ao.Basis.t;
   operators         : Operator.t list;
 }
+(* Type:2 ends here *)
 
+
+
+(* | ~nuclei~            | Nuclear coordinates used in the smiulation           |
+ * | ~charge~            | Total charge (electrons + nuclei)                    |
+ * | ~electrons~         | Electrons used in the simulation                     |
+ * | ~ao_basis~          | Atomic basis set                                     |
+ * | ~nuclear_repulsion~ | Nuclear repulsion energy                             |
+ * | ~operators~         | List of extra operators (range-separation, f12, etc) | *)
+
+
+(* [[file:../simulation.org::*Access][Access:2]] *)
 let nuclei            t = t.nuclei
 let charge            t = t.charge
 let electrons         t = t.electrons
 let ao_basis          t = t.ao_basis
 let nuclear_repulsion t = Nuclei.repulsion @@ nuclei t
 let operators         t = t.operators
+(* Access:2 ends here *)
 
-let make ?(multiplicity=1)
+
+
+(* Defaults:
+ * - multiplicity : ~1~
+ * - charge : ~0~
+ * - operators : ~[]~ *)
+
+
+(* [[file:../simulation.org::*Creation][Creation:2]] *)
+let make
+    ?(multiplicity=1)
     ?(charge=0)
     ?(operators=[])
     ~nuclei
@@ -36,10 +62,14 @@ let make ?(multiplicity=1)
     Charge.(Nuclei.charge nuclei + Electrons.charge electrons)
   in
 
-  {
-    charge ; nuclei ; electrons ; ao_basis ; 
-    operators
-  }
-
-
+  { charge ; nuclei ; electrons ; ao_basis ; operators}
+(* Creation:2 ends here *)
 
+(* [[file:../simulation.org::*Printers][Printers:2]] *)
+let pp ppf t =
+  let formula = Nuclei.formula t.nuclei in
+  let n_aos = Ao.Basis.size t.ao_basis in
+  let n_ops = List.length t.operators in
+  Format.fprintf ppf "@[@[%s@], @[%a@], @[%d AOs@], @[%d operators@]@]"
+    formula  Electrons.pp t.electrons n_aos n_ops
+(* Printers:2 ends here *)

simulation/lib/simulation.mli

@@ -1,32 +1,50 @@
-(* Contains the state of a simulation *)
+(* Simulation
+ *   :PROPERTIES:
+ *   :header-args: :noweb yes :comments both
+ *   :END:
+ * 
+ *   Contains the state of the simulation.
+ *   
+ *   #+NAME: open *)
 
+(* [[file:../simulation.org::open][open]] *)
 open Common
 open Particles
 open Operators
+(* open ends here *)
 
+(* Type
+ * 
+ *    #+NAME: types *)
+
+(* [[file:../simulation.org::types][types]] *)
 type t
+(* types ends here *)
 
+(* Access *)
+
+
+(* [[file:../simulation.org::*Access][Access:1]] *)
 val nuclei            : t -> Nuclei.t
-(** Nuclear coordinates used in the smiulation *)
-
 val charge            : t -> Charge.t
-(** Total charge (electrons + nuclei) *)
-
 val electrons         : t -> Electrons.t
-(** Electrons used in the simulation *)
-
 val ao_basis          : t -> Ao.Basis.t
-(** Atomic basis set *)
-
 val nuclear_repulsion : t -> float
-(** Nuclear repulsion energy *)
-
 val operators         : t -> Operator.t list
-(** List of extra operators (range-separation, f12, etc) *)
+(* Access:1 ends here *)
 
-(** {1 Creation} *)
+(* Creation *)
 
+
+(* [[file:../simulation.org::*Creation][Creation:1]] *)
 val make : ?multiplicity:int -> ?charge:int -> 
   ?operators:Operator.t list-> nuclei:Nuclei.t ->
   Ao.Basis.t -> t
+(* Creation:1 ends here *)
 
+(* Printers *)
+
+
+(* [[file:../simulation.org::*Printers][Printers:1]] *)
+val pp : Format.formatter -> t -> unit
+(* Printers:1 ends here *)

simulation/simulation.org

@@ -0,0 +1,131 @@
+#+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 
+
+* Simulation
+  :PROPERTIES:
+  :header-args: :noweb yes :comments both
+  :END:
+
+  Contains the state of the simulation.
+  
+  #+NAME: open
+  #+begin_src ocaml :tangle (eval mli)
+open Common
+open Particles
+open Operators
+  #+end_src
+
+   #+begin_src ocaml :tangle (eval ml) :exports none
+<<open>>
+   #+end_src
+
+** Type
+
+   #+NAME: types
+   #+begin_src ocaml :tangle (eval mli)
+type t
+   #+end_src
+
+   #+begin_src ocaml :tangle (eval ml) :exports none
+type t = {
+  charge            : Charge.t;
+  electrons         : Electrons.t;
+  nuclei            : Nuclei.t;
+  ao_basis          : Ao.Basis.t;
+  operators         : Operator.t list;
+}
+   #+end_src
+
+** Access
+
+   #+begin_src ocaml :tangle (eval mli)
+val nuclei            : t -> Nuclei.t
+val charge            : t -> Charge.t
+val electrons         : t -> Electrons.t
+val ao_basis          : t -> Ao.Basis.t
+val nuclear_repulsion : t -> float
+val operators         : t -> Operator.t list
+   #+end_src
+
+ | ~nuclei~            | Nuclear coordinates used in the smiulation           |
+ | ~charge~            | Total charge (electrons + nuclei)                    |
+ | ~electrons~         | Electrons used in the simulation                     |
+ | ~ao_basis~          | Atomic basis set                                     |
+ | ~nuclear_repulsion~ | Nuclear repulsion energy                             |
+ | ~operators~         | List of extra operators (range-separation, f12, etc) |
+
+   #+begin_src ocaml :tangle (eval ml) :exports none
+let nuclei            t = t.nuclei
+let charge            t = t.charge
+let electrons         t = t.electrons
+let ao_basis          t = t.ao_basis
+let nuclear_repulsion t = Nuclei.repulsion @@ nuclei t
+let operators         t = t.operators
+   #+end_src
+
+** Creation
+
+   #+begin_src ocaml :tangle (eval mli)
+val make : ?multiplicity:int -> ?charge:int -> 
+  ?operators:Operator.t list-> nuclei:Nuclei.t ->
+  Ao.Basis.t -> t
+   #+end_src
+
+   Defaults:
+   - multiplicity : ~1~
+   - charge : ~0~
+   - operators : ~[]~
+     
+   #+begin_src ocaml :tangle (eval ml) :exports none
+let make
+    ?(multiplicity=1)
+    ?(charge=0)
+    ?(operators=[])
+    ~nuclei
+    ao_basis
+  =
+
+  (* Tune Garbage Collector *)
+  let gc = Gc.get () in
+  Gc.set { gc with space_overhead = 1000  };
+
+  let electrons =
+    Electrons.of_atoms ~multiplicity ~charge nuclei
+  in
+
+  let charge = 
+    Charge.(Nuclei.charge nuclei + Electrons.charge electrons)
+  in
+
+  { charge ; nuclei ; electrons ; ao_basis ; operators}
+   #+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 =
+  let formula = Nuclei.formula t.nuclei in
+  let n_aos = Ao.Basis.size t.ao_basis in
+  let n_ops = List.length t.operators in
+  Format.fprintf ppf "@[@[%s@], @[%a@], @[%d AOs@], @[%d operators@]@]"
+    formula  Electrons.pp t.electrons n_aos n_ops
+   #+end_src
+
+   #+RESULTS:
+   : Line 2, characters 16-30:
+   : 2 |   let formula = Nuclei.formula t.nuclei in
+   :                     ^^^^^^^^^^^^^^
+   : Error: Unbound module Nuclei
+

top/lib/install_printers.ml

@@ -11,7 +11,9 @@ let printers =
     "Common.Zkey.pp" ;
     "Particles.Electrons.pp" ;
     "Particles.Element.pp" ;
+    "Particles.Nuclei.pp" ;
     "Particles.Zmatrix.pp" ;
+    "Simulation.Simulation.pp" ;
     
   ]