Merge branch 'master' of gitlab.com:scemama/QCaml

Basis/AOBasis.ml

@@ -5,6 +5,7 @@ type t =
 {
   basis              :  Basis.t                     ;
   overlap            :  Overlap.t             lazy_t;
+  multipole          :  Multipole.t           lazy_t;
   ortho              :  Orthonormalization.t  lazy_t;
   eN_ints            :  NucInt.t              lazy_t;
   kin_ints           :  KinInt.t              lazy_t;
@@ -17,6 +18,7 @@ type t =
 
 let basis        t = t.basis
 let overlap      t = Lazy.force t.overlap
+let multipole    t = Lazy.force t.multipole
 let ortho        t = Lazy.force t.ortho
 let eN_ints      t = Lazy.force t.eN_ints
 let kin_ints     t = Lazy.force t.kin_ints
@@ -62,7 +64,11 @@ let make ~cartesian ~basis ?f12 nuclei =
     ScreenedERI.of_basis basis
   ) in
 
-  { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ;
+  let multipole = lazy (
+    Multipole.of_basis basis
+  ) in
+
+  { basis ; overlap ; multipole ; ortho ; eN_ints ; kin_ints ; ee_ints ;
     ee_lr_ints ; f12_ints ; f12_over_r12_ints ; cartesian ;
   }
 

Basis/AOBasis.mli

@@ -10,6 +10,9 @@ val basis : t -> Basis.t
 val overlap : t -> Overlap.t
 (** Overlap matrix *)
 
+val multipole : t -> Multipole.t
+(** Multipole matrices *)
+
 val ortho : t -> Orthonormalization.t
 (** Orthonormalization matrix of the overlap *)
 

Basis/Basis.ml

@@ -93,6 +93,12 @@ let of_nuclei_and_basis_filename ~nuclei filename =
   in
   of_nuclei_and_general_basis nuclei general_basis
 
+let of_nuclei_and_basis_string ~nuclei str = 
+  let general_basis = 
+    GeneralBasis.of_string str
+  in
+  of_nuclei_and_general_basis nuclei general_basis
+
 
 let of_nuclei_and_basis_filenames ~nuclei filenames = 
   let general_basis = 

Basis/Basis.mli

@@ -18,6 +18,11 @@ val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> string -> t
   from a file. 
   *)
 
+val of_nuclei_and_basis_string : nuclei:Nuclei.t -> string -> t
+(** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t}
+  from a string. 
+  *)
+
 val of_nuclei_and_basis_filenames : nuclei:Nuclei.t -> string list -> t
 (** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t}
   from multiple files. 

Basis/GeneralBasis.ml

@@ -17,23 +17,24 @@ exception No_shell
 exception Malformed_shell of string
 
 
-let read_shell ic =
+let read_shell line_stream =
   try
 
     let shell, n =
       try
-        let line = input_line ic in
+        let line = Stream.next line_stream in
         if String.trim line = "$END" then raise End_of_file;
         Scanf.sscanf line " %c %d " (fun shell n -> shell, n)
       with
       | End_of_file -> raise No_shell
+      | Stream.Failure -> raise No_shell
       | Scanf.Scan_failure m -> raise (Malformed_shell m)
     in
 
     let rec loop = function
       | 0 -> []
       | i -> let contraction = 
-               let line = (input_line ic) in
+               let line = Stream.next line_stream in
                try Scanf.sscanf line " %_d %f %f "
                      (fun exponent coefficient -> { exponent ; coefficient })
                with _ -> raise (Malformed_shell (Printf.sprintf
@@ -50,30 +51,29 @@ let read_shell ic =
 
 
 
-let read_element ic =
+let read_element line_stream =
   try
 
+    let line = Stream.next line_stream in
     let element =
-      Scanf.sscanf (input_line ic) " %s " Element.of_string
+      Scanf.sscanf line " %s " Element.of_string
     in
 
     let rec loop () =
-      match read_shell ic with
+      match read_shell line_stream with
       | Some shell -> shell :: loop ()
       | None -> []
     in
     Some (element, Array.of_list (loop ()) )
 
   with
-  | End_of_file -> None
+  | Stream.Failure -> None
   
 
-
-let read filename =
-  let ic = open_in filename in
+let read_stream line_stream = 
   let rec loop accu =
     try
-      match read_element ic with
+      match read_element line_stream with
       | Some e -> loop (e :: accu)
       | None -> accu
     with
@@ -81,6 +81,17 @@ let read filename =
   in
   loop [] 
 
+let read filename =
+  let ic = open_in filename in
+  let line_stream =
+    Stream.from (fun _ -> 
+      try Some (input_line ic)
+      with End_of_file -> None )
+  in
+  let result = read_stream line_stream in
+  close_in ic;
+  result
+
 
 let combine basis_list =
   let h = Hashtbl.create 63 in
@@ -124,3 +135,8 @@ let to_string (name, contracted_shell_array) =
   Printf.sprintf "%s\n%s" name (string_of_contracted_shell_array contracted_shell_array)
 
 
+let of_string input_string = 
+  String.split_on_char '\n' input_string
+  |> Stream.of_list
+  |> read_stream
+

Basis/GeneralBasis.mli

@@ -46,7 +46,7 @@ exception Malformed_shell of string
 
 
 val read : string -> t
-(** Reads a basis set file and return an association list where
+(** Reads a basis-set file and return an association list where
     the key is an {!Element.t} and the value is the parsed basis set.
 *)
 
@@ -65,8 +65,10 @@ val read_many : string list -> t
 
 
 
-val read_element : in_channel -> element_basis option
-(** Reads an element from the input channel [ic]. For example,
+val read_element : string Stream.t -> element_basis option
+(** Reads an element from the input [string Stream]. The [string Stream] is a
+stream of lines, like a text file split on the end-of-line character.
+For example,
 {[
 HYDROGEN
 S   3
@@ -95,8 +97,10 @@ Some
 *)
   
 
-val read_shell : in_channel -> general_contracted_shell option
-(** Reads a shell from the input channel [ic]. For example,
+val read_shell : string Stream.t -> general_contracted_shell option
+(** Reads a shell from the input [string Stream]. The [string Stream] is a 
+stream of lines, like a text file split on the end-of-line character.
+For example,
 {[
 S   3
   1     13.0100000              0.0196850
@@ -119,4 +123,6 @@ Some
 val to_string : string * (general_contracted_shell array) -> string
 (** Pretty-prints the basis set of an {Element.t}. *)
 
+val of_string : string -> t
+(** Reads a GAMESS-formatted string. *)
 

Basis/Multipole.ml

@@ -0,0 +1,199 @@
+open Util
+open Constants
+open Lacaml.D
+
+type t = Mat.t array
+(*
+[| "x"; "y"; "z"; "x2"; "y2"; "z2" |]
+*)
+
+module Am  = AngularMomentum
+module Bs  = Basis
+module Co  = Coordinate
+module Cs  = ContractedShell
+module Csp = ContractedShellPair
+module Po  = Powers
+module Psp = PrimitiveShellPair
+
+let multiply a b =
+  let n = Mat.dim1 a in
+  let c = Mat.create n n in
+  Mat.cpab c a b;
+  c
+
+let x0 t = t.(0)
+let y0 t = t.(1)
+let z0 t = t.(2)
+let x1 t = t.(3)
+let y1 t = t.(4)
+let z1 t = t.(5)
+let x2 t = t.(6)
+let y2 t = t.(7)
+let z2 t = t.(8)
+
+let matrix_x  t = multiply (x1 t) @@ multiply (y0 t) (z0 t)
+let matrix_y  t = multiply (x0 t) @@ multiply (y1 t) (z0 t)
+let matrix_z  t = multiply (x0 t) @@ multiply (y0 t) (z1 t)
+let matrix_x2 t = multiply (x2 t) @@ multiply (y0 t) (z0 t)
+let matrix_y2 t = multiply (x0 t) @@ multiply (y2 t) (z0 t)
+let matrix_z2 t = multiply (x0 t) @@ multiply (y0 t) (z2 t)
+let matrix_xy t = multiply (x1 t) @@ multiply (y1 t) (z0 t)
+let matrix_yz t = multiply (x0 t) @@ multiply (y1 t) (z1 t)
+let matrix_zx t = multiply (x1 t) @@ multiply (y0 t) (z1 t)
+
+let cutoff = integrals_cutoff
+
+let to_powers x =
+  let open Zkey in
+  match to_powers x with
+  | Six x -> x
+  | _ -> assert false
+
+(** Computes all the integrals of the contracted shell pair *)
+let contracted_class shell_a shell_b : float Zmap.t array =
+
+  match Csp.make shell_a shell_b with
+  | None -> Array.init 9 (fun _ -> Zmap.create 0)
+  | Some shell_p ->
+    begin
+
+      (* Pre-computation of integral class indices *)
+      let class_indices = Csp.zkey_array shell_p in
+
+      let contracted_class =
+        Array.init 9 (fun _ -> Array.make (Array.length class_indices) 0.)
+      in
+
+      let a_minus_b =
+        Csp.a_minus_b shell_p
+      in
+      let norm_coef_scales =
+        Csp.norm_scales shell_p
+      in
+
+      (* Compute all integrals in the shell for each pair of significant shell pairs *)
+
+      let xyz_of_int k =
+        match k with
+        | 0 -> Co.X
+        | 1 -> Co.Y
+        | _ -> Co.Z
+      in
+
+      List.iter (fun (coef_prod, psp) ->
+          (** Screening on the product of coefficients *)
+          if (abs_float coef_prod) > 1.e-6*.cutoff then
+            begin
+              let expo_inv = Psp.exponent_inv psp
+              and center_pa = Psp.center_minus_a psp
+              and xa = Co.get Co.X @@ Cs.center shell_a
+              and ya = Co.get Co.Y @@ Cs.center shell_a
+              and za = Co.get Co.Z @@ Cs.center shell_a
+              in
+
+              Array.iteri (fun i key ->
+                  let (angMomA, angMomB) = to_powers key in
+                  (* 1D Overlap <i|j> *)
+                  let f k =
+                    let xyz = xyz_of_int k in
+                    Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB)
+                      expo_inv
+                      (Co.get xyz a_minus_b, Co.get xyz center_pa)
+                  in
+                  (* 1D <i|x-Xa|j> *)
+                  let g k =
+                    let xyz = xyz_of_int k in
+                    Overlap_primitives.hvrr (Po.get xyz angMomA + 1, Po.get xyz angMomB)
+                      expo_inv
+                      (Co.get xyz a_minus_b, Co.get xyz center_pa)
+                  in
+                  (* 1D <i|(x-Xa)^2|j> *)
+                  let h k =
+                    let xyz = xyz_of_int k in
+                    Overlap_primitives.hvrr (Po.get xyz angMomA + 2, Po.get xyz angMomB)
+                      expo_inv
+                      (Co.get xyz a_minus_b, Co.get xyz center_pa)
+                  in
+                  let norm =  norm_coef_scales.(i) in
+                  let f0, f1, f2, g0, g1, g2, h0, h1, h2 =
+                    f 0, f 1, f 2, g 0, g 1, g 2, h 0, h 1, h 2
+                  in
+                  let  x = g0 +. f0 *. xa in
+                  let  y = g1 +. f1 *. ya in
+                  let  z = g2 +. f2 *. za in
+                  let x2 = h0 +. xa *. (2. *. x -. xa *. f0) in
+                  let y2 = h1 +. ya *. (2. *. y -. ya *. f1) in
+                  let z2 = h2 +. za *. (2. *. z -. za *. f2) in
+                  let c = contracted_class in
+                  let d = coef_prod *. norm in
+                  c.(0).(i) <- c.(0).(i) +. d *. f0 ;
+                  c.(1).(i) <- c.(1).(i) +. d *. f1 ;
+                  c.(2).(i) <- c.(2).(i) +. d *. f2 ;
+                  c.(3).(i) <- c.(3).(i) +. d *.  x ;
+                  c.(4).(i) <- c.(4).(i) +. d *.  y ;
+                  c.(5).(i) <- c.(5).(i) +. d *.  z ;
+                  c.(6).(i) <- c.(6).(i) +. d *. x2 ;
+                  c.(7).(i) <- c.(7).(i) +. d *. y2 ;
+                  c.(8).(i) <- c.(8).(i) +. d *. z2 ;
+                ) class_indices
+            end
+        ) (Csp.coefs_and_shell_pairs shell_p);
+      let result =
+        Array.map (fun c -> Zmap.create (Array.length c) ) contracted_class
+      in
+      for j=0 to Array.length result do
+        let rj = result.(j) in
+        let cj = contracted_class.(j) in
+        Array.iteri (fun i key -> Zmap.add rj key cj.(i)) class_indices
+      done;
+      result
+    end
+
+
+(** Create multipole matrices *)
+let of_basis basis =
+  let to_powers x =
+    let open Zkey in
+    match to_powers x with
+    | Three x -> x
+    | _ -> assert false
+  in
+
+  let n     = Bs.size basis
+  and shell = Bs.contracted_shells basis
+  in
+
+  let result = Array.init 9 (fun _ -> Mat.create n n) in
+  for j=0 to (Array.length shell) - 1 do
+    for i=0 to j do
+      (* Compute all the integrals of the class *)
+      let cls =
+        contracted_class shell.(i) shell.(j)
+      in
+
+      for k=0 to 8 do
+        Array.iteri (fun j_c powers_j ->
+            let j_c = Cs.index shell.(j) + j_c + 1 in
+            let xj = to_powers powers_j in
+            Array.iteri (fun i_c powers_i ->
+                let i_c = Cs.index shell.(i) + i_c + 1 in
+                let xi = to_powers powers_i in
+                let key =
+                  Zkey.of_powers_six xi xj
+                in
+                let value =
+                  try Zmap.find cls.(k) key
+                  with Not_found -> 0.
+                in
+                result.(k).{i_c,j_c} <- value;
+                result.(k).{j_c,i_c} <- value;
+              ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i))))
+          ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j))))
+      done;
+    done;
+  done;
+  for k=0 to 8 do
+    Mat.detri result.(k);
+  done;
+  result
+

Basis/Multipole.mli

@@ -0,0 +1,34 @@
+(** Multipole atomic integrals:
+
+{% $$ \langle \chi_i | x | \chi_j \rangle $$ %}
+{% $$ \langle \chi_i | y | \chi_j \rangle $$ %}
+{% $$ \langle \chi_i | z | \chi_j \rangle $$ %}
+{% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %}
+{% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %}
+{% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %}
+
+*)
+
+open Lacaml.D
+
+type t
+
+val matrix_x : t -> Mat.t
+(** {% $$ \langle \chi_i | x | \chi_j \rangle $$ %} *)
+
+val matrix_y : t -> Mat.t
+(** {% $$ \langle \chi_i | y | \chi_j \rangle $$ %} *)
+
+val matrix_z : t -> Mat.t
+(** {% $$ \langle \chi_i | z | \chi_j \rangle $$ %} *)
+
+val matrix_x2 : t -> Mat.t
+(** {% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %} *)
+
+val matrix_y2 : t -> Mat.t
+(** {% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %} *)
+
+val matrix_z2 : t -> Mat.t
+(** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *)
+
+val of_basis : Basis.t -> t

Basis/Overlap.ml

@@ -28,62 +28,62 @@ let to_powers x =
 let contracted_class shell_a shell_b : float Zmap.t =
 
   match Csp.make shell_a shell_b with
+  | None -> Zmap.create 0
   | Some shell_p ->
     begin
 
-  (* Pre-computation of integral class indices *)
-  let class_indices = Csp.zkey_array shell_p in
+      (* Pre-computation of integral class indices *)
+      let class_indices = Csp.zkey_array shell_p in
 
-  let contracted_class = 
-    Array.make (Array.length class_indices) 0.
-  in
+      let contracted_class =
+        Array.make (Array.length class_indices) 0.
+      in
 
-  let a_minus_b = 
-    Csp.a_minus_b shell_p
-  in
-  let norm_coef_scales =
-    Csp.norm_scales shell_p
-  in
+      let a_minus_b =
+        Csp.a_minus_b shell_p
+      in
+      let norm_coef_scales =
+        Csp.norm_scales shell_p
+      in
 
-  (* Compute all integrals in the shell for each pair of significant shell pairs *)
+      (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
-  let xyz_of_int k =
-    match k with
-    | 0 -> Co.X
-    | 1 -> Co.Y
-    | _ -> Co.Z
-  in
+      let xyz_of_int k =
+        match k with
+        | 0 -> Co.X
+        | 1 -> Co.Y
+        | _ -> Co.Z
+      in
 
-  List.iter (fun (coef_prod, psp) ->
-      (** Screening on thr product of coefficients *)
-      if (abs_float coef_prod) > 1.e-3*.cutoff then
-        begin
-          let expo_inv = Psp.exponent_inv psp
-          and center_pa = Psp.center_minus_a psp
-          in
-
-          Array.iteri (fun i key ->
-              let (angMomA,angMomB) = to_powers key in
-              let f k = 
-                let xyz = xyz_of_int k in
-                Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB)  
-                  expo_inv
-                  (Co.get xyz a_minus_b,
-                  Co.get xyz center_pa) 
+      List.iter (fun (coef_prod, psp) ->
+          (** Screening on the product of coefficients *)
+          if (abs_float coef_prod) > 1.e-6*.cutoff then
+            begin
+              let expo_inv = Psp.exponent_inv psp
+              and center_pa = Psp.center_minus_a psp
               in
-              let norm =  norm_coef_scales.(i) in
-              let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in 
-              contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
-            ) class_indices
-        end
-    ) (Csp.coefs_and_shell_pairs shell_p);
-  let result = 
-    Zmap.create (Array.length contracted_class)
-  in
-  Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
-  result
-  end
-  | None -> Zmap.create 0
+
+              Array.iteri (fun i key ->
+                  let (angMomA,angMomB) = to_powers key in
+                  let f k =
+                    let xyz = xyz_of_int k in
+                    Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB)
+                      expo_inv
+                      (Co.get xyz a_minus_b,
+                      Co.get xyz center_pa)
+                  in
+                  let norm =  norm_coef_scales.(i) in
+                  let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in
+                  contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+                ) class_indices
+            end
+        ) (Csp.coefs_and_shell_pairs shell_p);
+      let result =
+        Zmap.create (Array.length contracted_class)
+      in
+      Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
+      result
+    end
 
 
 (** Create overlap matrix *)

INSTALL.md

@@ -1,7 +1,7 @@
 # Generic installation
 
 ```bash
-opam install ocamlbuild lacaml mpi getopt alcotest zarith
+opam install ocamlbuild ocamlfind lacaml mpi getopt alcotest zarith
 ```
 
 

Nuclei/Nuclei.ml

@@ -3,11 +3,7 @@ open Xyz_ast
 
 type t = (Element.t * Coordinate.t) array
 
-let of_xyz_file filename =
-  let lexbuf =
-    let ic = open_in filename in
-    Lexing.from_channel ic
-  in
+let of_xyz_lexbuf lexbuf =
   let data = 
     Xyz_parser.input Nuclei_lexer.read_all lexbuf
   in
@@ -24,6 +20,30 @@ let of_xyz_file filename =
   |> Array.of_list
 
 
+let of_xyz_string buffer =
+  Zmatrix.of_string buffer
+  |> Zmatrix.to_xyz
+  |> Array.map (fun (e,x,y,z) ->
+        (e, Coordinate.(angstrom_to_bohr @@ make_angstrom { x ; y ; z} ))
+      )
+
+
+let of_xyz_string input_string = 
+  Lexing.from_string input_string
+  |> of_xyz_lexbuf
+
+
+let of_xyz_file filename =
+  let ic = open_in filename in
+  let lexbuf =
+    Lexing.from_channel ic
+  in
+  let result = 
+    of_xyz_lexbuf lexbuf
+  in
+  close_in ic;
+  result
+
 
 let of_zmt_string buffer =
   Zmatrix.of_string buffer

Nuclei/Nuclei.mli

@@ -5,6 +5,9 @@ of tuples ({!Element.t}, {!Coordinate.t}).
 type t = (Element.t * Coordinate.t) array
 
 
+val of_xyz_string : string -> t
+(** Create from a string, in [xyz] format. *)
+
 val of_xyz_file : string -> t
 (** Create from a file, in [xyz] format. *)
 

Simulation.ml

@@ -1,6 +1,6 @@
 type t = {
   charge            : Charge.t;
-  electrons         : Electrons.t ;
+  electrons         : Electrons.t;
   nuclei            : Nuclei.t;
   basis             : Basis.t;
   ao_basis          : AOBasis.t;
@@ -19,12 +19,12 @@ let mu_erf            t = t.mu_erf
 let f12               t = t.f12
 
 let make ?cartesian:(cartesian=false)
-         ?multiplicity:(multiplicity=1)
-         ?charge:(charge=0)
-         ?f12
-         ?mu_erf
-         ~nuclei
-         basis
+      ?multiplicity:(multiplicity=1)
+      ?charge:(charge=0)
+      ?f12
+      ?mu_erf
+      ~nuclei
+      basis
   =
 
   (* Tune Garbage Collector *)