Added simulation

ao_basis/lib/ao_basis.ml

@@ -1,8 +1,117 @@
+open Qcaml_linear_algebra
+    
 type basis =
   | Unknown
   | Gaussian of Ao_basis_gaussian.t
 
 type t =
-  { basis : basis  ;
+  { ao_basis : basis  ;
     cartesian : bool
   }
+
+let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false)
+    ~nuclei filename =
+  match kind with
+  | `Gaussian ->
+      let basis =
+        Qcaml_gaussian_basis.Basis.of_nuclei_and_basis_filename  ~nuclei filename
+      in
+      let ao_basis = 
+        Gaussian (Ao_basis_gaussian.make ~basis ?operators ~cartesian nuclei )
+      in
+      { ao_basis ; cartesian }
+  | _ -> failwith "of_nuclei_and_basis_filename needs to be called with `Gaussian"
+
+let not_implemented () =
+  failwith "Not implemented"
+  
+let ao_basis t = t.ao_basis
+                    
+let size t = 
+  match t.ao_basis with
+  | Gaussian b -> Ao_basis_gaussian.size b
+  | _ -> not_implemented ()
+
+let overlap t =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.overlap b
+    | _ -> not_implemented ()
+  end
+  |> Matrix.relabel
+       
+let multipole t =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.multipole b
+    | _ -> not_implemented ()
+  end
+  |> Array.map Matrix.relabel
+       
+let ortho t =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.ortho b
+    | _ -> not_implemented ()
+  end
+  |> Matrix.relabel
+       
+let eN_ints t =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.eN_ints b
+    | _ -> not_implemented ()
+  end
+  |> Matrix.relabel
+       
+let kin_ints t =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.kin_ints b
+    | _ -> not_implemented ()
+  end
+  |> Matrix.relabel
+       
+let ee_ints t =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.ee_ints b
+    | _ -> not_implemented ()
+  end
+  |> Four_idx_storage.relabel
+       
+let ee_lr_ints t =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.ee_lr_ints b
+    | _ -> not_implemented ()
+  end
+  |> Four_idx_storage.relabel
+       
+let f12_ints t =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.f12_ints b
+    | _ -> not_implemented ()
+  end
+  |> Four_idx_storage.relabel
+       
+let f12_over_r12_ints t =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.f12_over_r12_ints b
+    | _ -> not_implemented ()
+  end
+  |> Four_idx_storage.relabel
+       
+let cartesian t = t.cartesian
+                    
+
+let values t point =
+  begin
+    match t.ao_basis with
+    | Gaussian b -> Ao_basis_gaussian.values b point
+    | _ -> not_implemented ()
+  end
+  |> Vector.relabel
+       

ao_basis/lib/ao_basis.mli

@@ -1,10 +1,63 @@
 (** Data structure for Atomic Orbitals. *)
 
+open Qcaml_common
+open Qcaml_particles
+open Qcaml_operators
+open Qcaml_linear_algebra
+
 type basis =
   | Unknown
   | Gaussian of Ao_basis_gaussian.t
 
-type t =
-  { basis : basis  ;
-    cartesian : bool
-  }
+type t 
+
+(** {1 Accessors} *)
+
+val size : t -> int
+(** Number of atomic orbitals in the AO basis set *)
+  
+val ao_basis : t -> basis
+(** One-electron basis set *)
+
+val overlap : t -> ('a,'a) Matrix.t
+(** Overlap matrix *)
+
+val multipole : t -> ('a,'a) Matrix.t array
+(** Multipole matrices *)
+
+val ortho : t -> ('a,'a) Matrix.t
+(** Orthonormalization matrix of the overlap *)
+
+val eN_ints : t -> ('a,'a) Matrix.t
+(** Electron-nucleus potential integrals *)
+
+val kin_ints : t -> ('a,'a) Matrix.t
+(** Kinetic energy integrals *)
+
+val ee_ints : t -> 'a Four_idx_storage.t
+(** Electron-electron potential integrals *)
+
+val ee_lr_ints : t -> 'a Four_idx_storage.t
+(** Electron-electron long-range potential integrals *)
+
+val f12_ints : t -> 'a Four_idx_storage.t
+(** Electron-electron potential integrals *)
+
+val f12_over_r12_ints : t -> 'a Four_idx_storage.t
+(** Electron-electron potential integrals *)
+
+val cartesian : t -> bool
+(** If true, use cartesian Gaussians (6d, 10f, ...) *)
+
+val values : t -> Coordinate.t -> 'a Vector.t
+(** Values of the AOs evaluated at a given point *)
+
+
+
+(** {1 Creators} *)
+
+val of_nuclei_and_basis_filename :
+  ?kind:[> `Gaussian ] -> ?operators:Operator.t list -> ?cartesian:bool ->
+  nuclei:Nuclei.t -> string -> t
+(** Creates the data structure for atomic orbitals from a {Basis.t} and the
+    molecular geometry {Nuclei.t} *)

ao_basis/lib/ao_basis_gaussian.ml

@@ -19,6 +19,7 @@ type t =
 }
 
 let basis        t = t.basis
+let size         t = Matrix.dim1 (Lazy.force t.overlap)
 let overlap      t = Lazy.force t.overlap
 let multipole    t = Lazy.force t.multipole
 let ortho        t = Lazy.force t.ortho

ao_basis/lib/ao_basis_gaussian.mli

@@ -6,24 +6,13 @@ open Qcaml_gaussian_basis
 open Qcaml_gaussian_integrals
 open Qcaml_operators
   
-type t =
-{
-  basis              :  Basis.t                     ;
-  overlap            :  Overlap.t             lazy_t;
-  multipole          :  Multipole.t           lazy_t;
-  ortho              :  Orthonormalization.t  lazy_t;
-  eN_ints            :  Electron_nucleus.t    lazy_t;
-  kin_ints           :  Kinetic.t             lazy_t;
-  ee_ints            :  Eri.t                 lazy_t;
-  ee_lr_ints         :  Eri_long_range.t      lazy_t;
-  f12_ints           :  F12.t                 lazy_t;
-  f12_over_r12_ints  :  Screened_eri.t        lazy_t;
-  cartesian          :  bool ;
-}
-
+type t 
 
 (** {1 Accessors} *)
 
+val size : t -> int
+(** Number of atomic orbitals *) 
+
 val basis : t -> Basis.t
 (** One-electron basis set *)
 

ao_basis/test/ao_basis.ml

@@ -0,0 +1,100 @@
+open Alcotest
+open Qcaml_common
+open Qcaml_ao_basis
+open Qcaml_particles
+open Qcaml_operators
+open Qcaml_linear_algebra
+open Ao_basis
+
+let wd = Qcaml.root ^ Filename.dir_sep ^ "test"
+ 
+let make_tests name t =
+
+  let check_matrix title a r =
+    let ax =
+      Matrix.as_vec_inplace a
+      |> Vector.to_bigarray_inplace
+    in
+    Matrix.as_vec_inplace r
+    |> Vector.iteri (fun i x ->
+          let message =
+            Printf.sprintf "%s line %d" title i
+          in
+          check (float 1.e-10) message ax.{i} x
+        ) 
+  in
+
+  let check_eri a r =
+    let f { Four_idx_storage.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } =
+      (i_r1, (j_r2, (k_r1, (l_r2, value))))
+    in
+    let a = Four_idx_storage.to_list a |> List.rev_map f in
+    let r = Four_idx_storage.to_list r |> List.rev_map f in
+    check (list (pair int (pair int (pair int (pair int (float 1.e-10)))))) "ERI" a r
+  in
+
+  let test_overlap () =
+    let reference =
+      Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref")
+    in
+    check_matrix "Overlap" (overlap t) reference
+  in
+
+  let test_eN_ints () =
+    let reference =
+      Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref")
+    in
+    check_matrix "eN_ints" (eN_ints t) reference
+  in
+
+  let test_kin_ints () =
+    let reference =
+      Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^  "_kin.ref")
+    in
+    check_matrix "kin_ints" (kin_ints t) reference
+  in
+
+  let test_ee_ints () =
+    let reference =
+      Four_idx_storage.of_file (wd ^ Filename.dir_sep ^ name ^ "_eri.ref")
+        ~sparsity:`Dense ~size:(size t)
+    in
+    check_eri (ee_ints t) reference
+    ;
+  in
+
+  let test_ee_lr_ints () =
+    let reference =
+      Four_idx_storage.of_file (wd ^ Filename.dir_sep ^ name ^ "_eri_lr.ref")
+        ~sparsity:`Dense ~size:(size t)
+    in
+    check_eri (ee_lr_ints t) reference
+  in
+
+  [
+    "Overlap",     `Quick,  test_overlap;
+    "eN_ints",     `Quick,  test_eN_ints;
+    "kin_ints",    `Quick,  test_kin_ints;
+    "ee_ints",     `Quick,  test_ee_ints;
+    "ee_lr_ints",  `Quick,  test_ee_lr_ints;
+  ]
+  
+let tests =
+  List.concat [ 
+    let nuclei =
+      wd ^ Filename.dir_sep ^ "water.xyz"
+    |> Nuclei.of_xyz_file 
+    in
+    let basis_filename =
+      wd ^ Filename.dir_sep ^ "cc-pvdz"
+    in
+    let rs = 0.5 in
+    let operators = [ Operator.of_range_separation rs ] in
+    let ao_basis = 
+      Ao_basis.of_nuclei_and_basis_filename ~kind:`Gaussian
+        ~operators ~cartesian:false ~nuclei basis_filename
+    in
+    make_tests "water" ao_basis ;
+
+  ]
+    

ao_basis/test/ao_basis_gaussian.ml

@@ -35,65 +35,41 @@ let make_tests name t =
     check (list (pair int (pair int (pair int (pair int (float 1.e-10)))))) "ERI" a r
   in
 
-  let check_eri_lr a r =
-    let f { Eri_long_range.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } =
-      (i_r1, (j_r2, (k_r1, (l_r2, value))))
-    in
-    let a = Eri_long_range.to_list a |> List.rev_map f in
-    let r = Eri_long_range.to_list r |> List.rev_map f in
-    check (list (pair int (pair int (pair int (pair int (float 1.e-10)))))) "ERI_lr" a r
-  in
-
   let test_overlap () =
     let reference =
       Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref")
     in
-    let overlap =
-      Lazy.force t.overlap |> Overlap.matrix
-    in
-    check_matrix "Overlap" overlap reference
+    check_matrix "Overlap" (overlap t) reference
   in
 
   let test_eN_ints () =
     let reference =
       Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref")
     in
-    let eN_ints =
-      Lazy.force t.eN_ints |> Electron_nucleus.matrix
-    in
-    check_matrix "eN_ints" eN_ints reference
+    check_matrix "eN_ints" (eN_ints t) reference
   in
 
   let test_kin_ints () =
     let reference =
       Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^  "_kin.ref")
     in
-    let kin_ints =
-      Lazy.force t.kin_ints |> Kinetic.matrix
-    in
-    check_matrix "kin_ints" kin_ints reference
+    check_matrix "kin_ints" (kin_ints t) reference
   in
 
   let test_ee_ints () =
     let reference =
-      Eri.of_file (wd ^ Filename.dir_sep ^ name ^ "_eri.ref") ~sparsity:`Dense ~size:(Basis.size t.basis)
+      Eri.of_file (wd ^ Filename.dir_sep ^ name ^ "_eri.ref") ~sparsity:`Dense ~size:(size t)
     in
-    let ee_ints =
-      Lazy.force t.ee_ints
-    in
-    check_eri ee_ints reference
+    check_eri (ee_ints t) reference
     ;
   in
 
   let test_ee_lr_ints () =
     let reference =
       Eri_long_range.of_file (wd ^ Filename.dir_sep ^ name ^ "_eri_lr.ref") ~sparsity:`Dense
-          ~size:(Basis.size t.basis)
+          ~size:(size t)
     in
-    let ee_lr_ints =
-      Lazy.force t.ee_lr_ints
-    in
-    check_eri_lr ee_lr_ints reference
+    check_eri (ee_lr_ints t) reference
   in
 
   [

gaussian_integrals/lib/electron_nucleus.ml

@@ -12,8 +12,6 @@ module Bs = Basis
 module Cs = Contracted_shell
 
 type t = (Bs.t, Bs.t) Matrix.t
-external matrix : t -> (Bs.t, Bs.t) Matrix.t = "%identity"
-external of_matrix : (Bs.t, Bs.t) Matrix.t -> t = "%identity"
 
 (** (0|0)^m : Fundamental electron-nucleus repulsion integral
     $ \int \phi_p(r1) 1/r_{C} dr_1 $

gaussian_integrals/lib/kinetic.ml

@@ -14,8 +14,6 @@ module Psp = Primitive_shell_pair
 module Ps  = Primitive_shell
 
 type t = (Basis.t, Basis.t) Matrix.t
-external matrix : t -> (Basis.t, Basis.t) Matrix.t = "%identity"
-external of_matrix : (Basis.t, Basis.t) Matrix.t -> t = "%identity"
 
 let cutoff = integrals_cutoff
 

gaussian_integrals/lib/matrix_on_basis.mli

@@ -13,9 +13,3 @@ val of_basis_pair : Basis.t -> Basis.t -> t
 
 val to_file : filename:string -> t -> unit
 (** Write the integrals in a file. *)
-
-val matrix : t -> (Basis.t, Basis.t) Matrix.t
-(** Returns the matrix. *)
-
-val of_matrix : (Basis.t, Basis.t) Matrix.t -> t 
-(** Build from a matrix. *)

gaussian_integrals/lib/multipole.mli

@@ -12,7 +12,7 @@
 open Qcaml_linear_algebra
 open Qcaml_gaussian_basis
 
-type t
+type t = (Basis.t, Basis.t) Matrix.t array
 
 val matrix_x : t -> (Basis.t, Basis.t) Matrix.t
 (** {% $$ \langle \chi_i | x | \chi_j \rangle $$ %} *)

gaussian_integrals/lib/orthonormalization.ml

@@ -13,7 +13,7 @@ type t = (Bs.t, Bs.t) Matrix.t
   
 let make_canonical ~thresh ~basis ~cartesian ~overlap =
 
-  let overlap_matrix =  Ov.matrix overlap in
+  let overlap_matrix =  overlap in
 
   let make_canonical_spherical basis =
     let ao_num = Bs.size basis in
@@ -47,9 +47,8 @@ let make_canonical ~thresh ~basis ~cartesian ~overlap =
 
 let make_lowdin ~thresh ~overlap =
   
-  let overlap_matrix = Ov.matrix overlap in
   let u_vec, u_val =
-    Matrix.diagonalize_symm overlap_matrix
+    Matrix.diagonalize_symm overlap
   in
   let u_vec_x = Matrix.to_bigarray_inplace u_vec in
 

gaussian_integrals/lib/overlap.ml

@@ -13,8 +13,6 @@ module Po  = Powers
 module Psp = Primitive_shell_pair
 
 type t = (Basis.t, Basis.t) Matrix.t
-external matrix : t -> (Basis.t, Basis.t) Matrix.t = "%identity"
-external of_matrix : (Basis.t, Basis.t) Matrix.t -> t = "%identity"
 
 let cutoff = integrals_cutoff
 

linear_algebra/lib/four_idx_storage.ml

@@ -18,6 +18,17 @@ type 'a t =
     four_index       : 'a storage_t ;
   }
 
+let relabel t =
+  { size = t.size ;
+    two_index = Matrix.relabel t.two_index ; 
+    two_index_anti = Matrix.relabel t.two_index_anti ; 
+    three_index = Matrix.relabel t.three_index ; 
+    three_index_anti = Matrix.relabel t.three_index_anti ; 
+    four_index = match t.four_index with
+      | Dense x -> Dense (Matrix.relabel x)
+      | Sparse x -> Sparse x
+  }
+
 let key_of_indices ~r1 ~r2 =
   let { first=i ; second=k } = r1 and { first=j ; second=l } = r2 in
   let f i k = 

linear_algebra/lib/four_idx_storage.mli

@@ -77,3 +77,4 @@ val of_file : size:int -> sparsity:[< `Dense | `Sparse ]
   -> Scanf.Scanning.file_name -> 'a t
 (** Read from a text file with format ["%d %d %d %d %f"]. *)
 
+val relabel : 'a t -> 'b t

simulation/lib/dune

@@ -0,0 +1,14 @@
+; name = name of the supermodule that will wrap all source files as submodules
+; public_name = name of the library for ocamlfind and opam
+(library
+ (name qcaml_simulation)
+ (public_name qcaml.simulation)
+ (libraries
+        qcaml.common
+        qcaml.particles
+        qcaml.gaussian_basis
+        qcaml.gaussian_integrals
+	qcaml.operators
+	qcaml.ao_basis
+ )
+ (synopsis "Characterizes a simulation."))

simulation/lib/simulation.ml

@@ -0,0 +1,46 @@
+open Qcaml_common
+open Qcaml_particles
+open Qcaml_ao_basis
+open Qcaml_operators
+    
+type t = {
+  charge            : Charge.t;
+  electrons         : Electrons.t;
+  nuclei            : Nuclei.t;
+  ao_basis          : Ao_basis.t;
+  operators         : Operator.t list;
+}
+
+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
+
+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
+  }
+
+
+

simulation/lib/simulation.mli

@@ -0,0 +1,33 @@
+(* Contains the state of a simulation *)
+
+open Qcaml_common
+open Qcaml_particles
+open Qcaml_ao_basis
+open Qcaml_operators
+    
+type t
+
+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) *)
+
+(** {1 Creation} *)
+
+val make : ?multiplicity:int -> ?charge:int -> 
+  ?operators:Operator.t list-> nuclei:Nuclei.t ->
+  Ao_basis.t -> t
+

test/run_tests.ml

@@ -10,6 +10,7 @@ let test_suites: unit Alcotest.test list = [
   "Particles.Electrons", Test_particles.Electrons.tests;
   "Gaussian_basis.General_basis", Test_gaussian_basis.General_basis.tests;
   "Ao_basis.Ao_basis_gaussian", Test_ao_basis.Ao_basis_gaussian.tests;
+  "Ao_basis.Ao_basis", Test_ao_basis.Ao_basis.tests;
 ]
 
 let () = Alcotest.run "QCaml" test_suites