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

.merlin

@@ -1,12 +1,8 @@
-PKG str unix bigarray lacaml zarith mpi alcotest
+EXCLUDE_QUERY_DIR
+B _build/default/.run_hartree_fock.eobjs/byte
+B _build/default/Basis/.Basis.objs/byte
+B _build/default/Nuclei/.Nuclei.objs/byte
 S .
-S Test
-S Nuclei
-S Parallel
-S CI
-S Utils
 S Basis
-S SCF
+S Nuclei
-S MOBasis
+FLG -w @a-4-29-40-41-42-44-45-48-58-59-60-40 -strict-sequence -strict-formats -short-paths -keep-locs
-S CI
-B _build/**

Basis/TwoElectronRRVectorized.ml

@@ -756,14 +756,14 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
             in
             let empty = Array.make (maxm+1) 0. in
             let center_qc_tmp = Array.init nq (fun cd ->
-                Bohr.make { Point.
+                Coordinate.make { Coordinate.
                   x = (center_qc Co.X).(cd) ;
                   y = (center_qc Co.Y).(cd) ;
                   z = (center_qc Co.Z).(cd) ;
                 })
             in
             Array.iteri (fun ab shell_ab ->
-                let center_pa = Bohr.make { Point.
+                let center_pa = Coordinate.make { Coordinate.
                     x = (center_pa Co.X).(ab) ;
                     y = (center_pa Co.Y).(ab) ;
                     z = (center_pa Co.Z).(ab) ;
@@ -789,7 +789,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                         in
                         zero_m Zp.{
                           maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
-                          center_pq = Bohr.make Point.{x ; y ; z} ;
+                          center_pq = Coordinate.make Coordinate.{x ; y ; z} ;
                           center_pa ; center_qc = center_qc_tmp.(cd) ;
                           normalization ;
                         }

Nuclei/Nuclei.ml

@@ -28,7 +28,9 @@ let of_xyz_file filename =
 let of_zmt_string buffer =
   Zmatrix.of_string buffer
   |> Zmatrix.to_xyz
-  |> Array.map (fun (e,x,y,z) -> (e, Coordinate.angstrom_to_bohr (Angstrom.make {Point.x ; y ; z} )))
+  |> Array.map (fun (e,x,y,z) ->
+        (e, Coordinate.(angstrom_to_bohr @@ make_angstrom { x ; y ; z} ))
+      )
     
 
 let of_zmt_file filename =
@@ -56,12 +58,13 @@ let to_string atoms =
 -----------------------------------------------------------------------
 " ^
     (Array.mapi (fun i (e, coord) ->
+         let open Coordinate in
          let coord =
-           Coordinate.bohr_to_angstrom coord
+           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) 
-            coord.Angstrom.x coord.Angstrom.y coord.Angstrom.z
+            coord.x coord.y coord.z
        ) atoms
      |> Array.to_list
      |> String.concat "\n" ) ^
@@ -100,11 +103,12 @@ let charge nuclei =
 let to_xyz_string t =
   [ string_of_int (Array.length t) ; "" ] @
   ( Array.mapi (fun i (e, coord) ->
+         let open Coordinate in
          let coord =
-           Coordinate.bohr_to_angstrom coord
+           bohr_to_angstrom coord
          in
          Printf.sprintf " %5s   %12.6f  %12.6f  %12.6f"
-            (Element.to_string e) coord.Angstrom.x coord.Angstrom.y coord.Angstrom.z
+            (Element.to_string e) coord.x coord.y coord.z
        ) t 
     |> Array.to_list )
   |> String.concat "\n"
@@ -117,8 +121,9 @@ let to_t2_string t =
     (2 * Array.length t); 
     "# Znuc   x            y           z" ]
   @ (Array.mapi (fun i (e, coord) ->
+         let open Coordinate in
          Printf.sprintf " %5f    %12.6f  %12.6f  %12.6f"
-            (Element.to_int e |> float_of_int)  coord.Bohr.x coord.Bohr.y coord.Bohr.z
+            (Element.to_int e |> float_of_int)  coord.x coord.y coord.z
        ) t 
     |> Array.to_list )
   |> String.concat "\n"

Nuclei/Nuclei.mli

@@ -2,7 +2,7 @@
 of tuples ({!Element.t}, {!Coordinate.t}).
 *)
 
-type t = (Element.t * Bohr.t) array
+type t = (Element.t * Coordinate.t) array
 
 
 val of_xyz_file : string -> t

Nuclei/Xyz_ast.mli

@@ -4,7 +4,7 @@
 type nucleus =
   {
     element: Element.t ;
-    coord  : Angstrom.t;
+    coord  : Coordinate.angstrom Coordinate.point;
   }
 
 type xyz_file =

Nuclei/Xyz_parser.mly

@@ -3,10 +3,9 @@
 %{
 exception InputError of string
 let make_angstrom x y z =
-  Angstrom.make {
+  Coordinate.(make_angstrom {
-    Point.
     x ; y ; z
-  }
+  })
 
 let output_of f x y z =
   let a = make_angstrom x y z in

Utils/Angstrom.ml

@@ -1 +0,0 @@
-include Bohr

Utils/Angstrom.mli

@@ -1,4 +0,0 @@
-(** 3D point in cartesian coordinates, where units are Angstroms. *)
-
-include module type of Bohr
-

Utils/Bohr.ml

@@ -1,11 +0,0 @@
-type t = {
-  x : float ;
-  y : float ;
-  z : float ;
-}
-
-
-external make : Point.t -> t = "%identity"
-  
-
-

Utils/Bohr.mli

@@ -1,12 +0,0 @@
-(** 3D point in cartesian coordinates, where units are atomic units (Bohr). *)
-
-type t = private {
-  x : float ;
-  y : float ;
-  z : float ;
-}
-
-
-val make : Point.t -> t
-
-

Utils/Cholesky.ml

@@ -3,13 +3,16 @@ open Lacaml.D
 let full_ldl m_A =
 
   let n = Mat.dim1 m_A in
+
   assert (Mat.dim2 m_A = n);
 
   let v_D  = Vec.make0 n in
+
   let m_Lt = Mat.identity n in
 
   let v_D_inv  = Vec.make0 n in
 
+
   let compute_d j =
     let l_jk = 
       Mat.col m_Lt j
@@ -20,6 +23,7 @@ let full_ldl m_A =
     m_A.{j,j} -. dot ~n:(j-1) l_jk  l_jk__d_k
   in
 
+
   let compute_l i =
     let l_ik__d_k = 
       Mat.col m_Lt i
@@ -33,6 +37,7 @@ let full_ldl m_A =
       v_D_inv.{j} *. ( m_A.{j,i} -. dot ~n:(j-1) l_ik__d_k l_jk )
   in
 
+
   for i=1 to n do
     for j=1 to (i-1) do
       m_Lt.{j,i} <- compute_l i j;
@@ -85,6 +90,7 @@ let pivoted_ldl threshold m_A =
       v_D_inv.{j} *. ( m_A.{pi.(j),pi.(i)} -. dot ~n:(j-1) l_ik__d_k l_jk )
   in
 
+
   let maxloc v i =
     let rec aux pos value i =
       if i > n then
@@ -106,6 +112,8 @@ let pivoted_ldl threshold m_A =
   in                                                                                          
 
 
+
+
   let () = 
     try
       for i=1 to n do
@@ -126,6 +134,10 @@ let pivoted_ldl threshold m_A =
 
 
 
+let make_ldl ?(threshold=Constants.epsilon) m_A =
+  pivoted_ldl m_A
+
+
 
 
 

Utils/Coordinate.ml

@@ -1,7 +1,19 @@
-type bohr = Bohr.t
+type bohr 
-type angstrom = Angstrom.t
+type angstrom 
+
+type xyz = {
+  x : float ;
+  y : float ;
+  z : float ;
+}
+
+type 'a point = xyz
+
+type t = bohr point
+
+external make : 'a point -> t = "%identity"
+external make_angstrom : 'a point -> angstrom point = "%identity"
   
-type t = bohr
 
 type axis = X | Y | Z
 
@@ -9,19 +21,17 @@ type axis = X | Y | Z
 let a_to_b a = Constants.a0_inv *. a
 let b_to_a b = Constants.a0     *. b
 
-let bohr_to_angstrom { Bohr.x ; y ; z } = 
+let bohr_to_angstrom { x ; y ; z } =
-  Angstrom.make 
+  make 
   {
-    Point.
     x = b_to_a x ;
     y = b_to_a y ;
     z = b_to_a z ;
   }
 
-let angstrom_to_bohr { Angstrom.x ; y ; z } =  
+let angstrom_to_bohr { x ; y ; z } =  
-  Bohr.make 
+  make 
   {
-    Point.
     x = a_to_b x ;
     y = a_to_b y ;
     z = a_to_b z ;
@@ -29,30 +39,27 @@ let angstrom_to_bohr { Angstrom.x ; y ; z } =
 
 
 let zero =
-  Bohr.make { Point.x = 0. ; y = 0. ; z = 0. } 
+  make { x = 0. ; y = 0. ; z = 0. } 
 
 
 (** Linear algebra *)
-let ( |. )  s { Bohr.x ; y ; z } =
+let ( |. )  s { x ; y ; z } =
-    Bohr.make {
+    make {
-      Point.
       x = s *. x ;
       y = s *. y ;
       z = s *. z ;
     }
 
-let ( |+ ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
+let ( |+ ) { x = x1 ; y = y1 ; z = z1 } { x = x2 ; y = y2 ; z = z2 } =
-  Bohr.make {
+  make {
-    Point.
     x = x1 +. x2 ;
     y = y1 +. y2 ;
     z = z1 +. z2 ;
   }
 
 
-let ( |- ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
+let ( |- ) { x = x1 ; y = y1 ; z = z1 } { x = x2 ; y = y2 ; z = z2 } =
-  Bohr.make {
+  make {
-    Point.
     x = x1 -. x2 ;
     y = y1 -. y2 ;
     z = z1 -. z2 ;
@@ -62,7 +69,7 @@ let ( |- ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
 let neg a = -1. |. a
 
 
-let dot { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
+let dot { x = x1 ; y = y1 ; z = z1 } { x = x2 ; y = y2 ; z = z2 } =
   x1 *. x2 +. y1 *. y2 +. z1 *. z2
 
 
@@ -70,7 +77,7 @@ let norm u =
   sqrt ( dot u u )
 
 
-let get axis { Bohr.x ; y ; z } = 
+let get axis { x ; y ; z } = 
   match axis with
   | X -> x 
   | Y -> y 
@@ -79,14 +86,12 @@ let get axis { Bohr.x ; y ; z } =
 
 open Format
 let pp ppf c = 
-  let open Bohr in
   fprintf ppf "@[@[%8.4f@] @[%8.4f@] @[%8.4f@]@]" c.x c.y c.z
 
 let pp_bohr ppf c = 
-  let open Bohr in
   fprintf ppf "@[(@[%10f@], @[%10f@], @[%10f@] Bohr)@]" c.x c.y c.z
 
 let pp_angstrom ppf c = 
   let c = bohr_to_angstrom c in
-  let open Angstrom in
   fprintf ppf "@[(@[%10f@], @[%10f@], @[%10f@] Angs)@]" c.x c.y c.z
+

Utils/Coordinate.mli

@@ -1,92 +1,108 @@
 (** Type for coordinates in 3D space.
 
   All operations on points are done in atomic units. Therefore, 
-  all coordinates are given in {!Bohr.t} and manipulated with this
+  all coordinates are given in bohr and manipulated with this
   module.
  *)
 
-type t = Bohr.t
+type bohr 
+type angstrom 
+
+type xyz = {
+  x : float ;
+  y : float ;
+  z : float ;
+}
+
+type 'a point = xyz
+
+type t = bohr point
+
+val make : 'a point -> t
+(** Creates a point in atomic units *)
+
+val make_angstrom : 'a point -> angstrom point
+(** Creates a point in Angstrom *)
 
 type axis = X | Y | Z
 
-val bohr_to_angstrom : Bohr.t -> Angstrom.t 
+val bohr_to_angstrom : bohr point -> angstrom point
-(** Converts a point in Bohr to Angstrom. *)
+(** Converts a point in bohr to angstrom. *)
 
 
-val angstrom_to_bohr : Angstrom.t -> Bohr.t
+val angstrom_to_bohr : angstrom point -> bohr point
-(** Converts a point in Angstrom to Bohr. *)
+(** Converts a point in Angstrom to bohr. *)
 
 
-val zero : Bohr.t
+val zero : bohr point
-(** [zero = { Bohr.x = 0. ; y=0. ; z=0. }] *)
+(** [zero = { x = 0. ; y=0. ; z=0. }] *)
 
 
-val get : axis -> Bohr.t -> float
+val get : axis -> bohr point -> float
 (** Extracts the projection of the coordinate on an axis.
 
   Example:
   
-  [Coordinate.(get Y) { Bohr.x=1. ; y=2. ; z=3. } -> 2.]
+  [Coordinate.(get Y) { x=1. ; y=2. ; z=3. } -> 2.]
 
 *)
 
-
 (** {1 Vector operations} *)
 
 
-val ( |. ) : float -> Bohr.t -> Bohr.t
+val ( |. ) : float -> t -> t
 (** Scale by a float.
 
   Example:
 
-  [ 2. |. { Bohr.x=1. ; y=2. ; z=3. } -> { Bohr.x=2. ; y=4. ; z=6. } ]
+  [ 2. |. { x=1. ; y=2. ; z=3. } -> { x=2. ; y=4. ; z=6. } ]
 
  *)
 
 
-val ( |+ ) : Bohr.t -> Bohr.t -> Bohr.t
+val ( |+ ) : t -> t -> t
 (** Add two vectors.
 
   Example:
 
-  {[{ Bohr.x=1. ; y=2. ; z=3. } |+ { Bohr.x=2. ; y=3. ; z=1. } ->
+  {[{ x=1. ; y=2. ; z=3. } |+ { x=2. ; y=3. ; z=1. } ->
-    { Bohr.x=3. ; y=5. ; z=4. }]}
+    { x=3. ; y=5. ; z=4. }]}
 
  *)
 
 
-val ( |- ) : Bohr.t -> Bohr.t -> Bohr.t
+val ( |- ) : t -> t -> t
 (** Subtract two vectors.
 
   Example:
 
-  {[{ Bohr.x=1. ; y=2. ; z=3. } |- { Bohr.x=2. ; y=3. ; z=1. } ->
+  {[{ x=1. ; y=2. ; z=3. } |- { x=2. ; y=3. ; z=1. } ->
-    { Bohr.x=-1. ; y=-1. ; z=2. }]}
+    { x=-1. ; y=-1. ; z=2. }]}
 
  *)
 
 
-val neg : Bohr.t -> Bohr.t
+val neg : t -> t
 (** Example:
 
-{[Coordinate.neg { Bohr.x=1. ; y=2. ; z=-3. } ->
+{[Coordinate.neg { x=1. ; y=2. ; z=-3. } ->
-    { Bohr.x=-1. ; y=-2. ; z=3. }]}
+    { x=-1. ; y=-2. ; z=3. }]}
 
  *)
 
-val dot : Bohr.t -> Bohr.t -> float
+val dot : t -> t -> float
 (** Dot product. *)
 
 
-val norm : Bohr.t -> float
+val norm : t -> float
 (** L{^2} norm of the vector. *)
 
 
 (** {2 Printers} *)
 
-val pp: Format.formatter -> Bohr.t -> unit
+val pp: Format.formatter -> t -> unit
 
-val pp_bohr: Format.formatter -> Bohr.t -> unit
+val pp_bohr: Format.formatter -> t -> unit
 
-val pp_angstrom : Format.formatter -> Bohr.t -> unit
+val pp_angstrom : Format.formatter -> t -> unit