Working n multipole integrals

Basis/Multipole.ml

@@ -141,7 +141,7 @@ let contracted_class shell_a shell_b : float Zmap.t array =
       let result =
         Array.map (fun c -> Zmap.create (Array.length c) ) contracted_class
       in
-      for j=0 to Array.length result do
+      for j=0 to Array.length result -1 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

MOBasis/MOBasis.ml

@@ -12,7 +12,7 @@ type mo_type =
   | Natural of string
   | Localized of string
 
-type t = 
+type t =
   {
     simulation      : Simulation.t;      (* Simulation which produced the MOs *)
     mo_type         : mo_type;           (* Kind of MOs (RHF, CASSCF, Localized...) *)
@@ -27,7 +27,7 @@ type t =
 
 
 let size t =
-  Mat.dim2 t.mo_coef 
+  Mat.dim2 t.mo_coef
 
 let simulation    t = t.simulation
 let mo_type       t = t.mo_type
@@ -42,9 +42,9 @@ let f12_ints      t = Lazy.force t.f12_ints
 let one_e_ints    t = Lazy.force t.one_e_ints
 
 
-let mo_energies t = 
+let mo_energies t =
   let m_C = mo_coef t in
-  let f = 
+  let f =
     let m_N = Mat.of_diag @@ mo_occupation t in
     let m_P = x_o_xt m_N m_C in
     match t.mo_type with
@@ -59,7 +59,7 @@ let mo_energies t =
   |> Mat.copy_diag
 
 
-let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = 
+let mo_matrix_of_ao_matrix ~mo_coef ao_matrix =
   xt_o_x ~x:mo_coef ~o:ao_matrix
 
 
@@ -81,16 +81,16 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef ()  =
   and kin_ints = lazy (
     AOBasis.kin_ints ao_basis
     |> KinInt.matrix
-    |> mo_matrix_of_ao_matrix ~mo_coef  
+    |> mo_matrix_of_ao_matrix ~mo_coef
     |> KinInt.of_matrix
   )
   and ee_ints = lazy (
     AOBasis.ee_ints ao_basis
-    |> ERI.four_index_transform mo_coef   
+    |> ERI.four_index_transform mo_coef
   )
   and f12_ints = lazy (
     AOBasis.f12_ints ao_basis
-    |> F12.four_index_transform mo_coef   
+    |> F12.four_index_transform mo_coef
   )
   in
   let one_e_ints = lazy (
@@ -102,8 +102,12 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef ()  =
     f12_ints }
 
 
+let values t point =
+  let c = mo_coef t in
+  let a = AOBasis.values (Simulation.ao_basis t.simulation) point in
+  gemv ~trans:`T c a
 
-let of_hartree_fock hf = 
+let of_hartree_fock hf =
   let mo_coef       = HF.eigenvectors hf in
   let simulation    = HF.simulation   hf in
   let mo_occupation = HF.occupation   hf in
@@ -118,7 +122,7 @@ let of_hartree_fock hf =
 
 let of_mo_basis simulation other =
 
-  let mo_coef = 
+  let mo_coef =
     let basis = Simulation.ao_basis simulation in
     let basis_other = ao_basis other in
     let m_S =
@@ -151,9 +155,9 @@ let of_mo_basis simulation other =
       else 0.)
   in
   make ~simulation ~mo_type:Projected ~mo_occupation ~mo_coef ()
-  
 
-  
+
+
 
 
 let pp ?(start=1) ?(finish=0) ppf t =
@@ -167,7 +171,7 @@ let pp ?(start=1) ?(finish=0) ppf t =
       | x -> x
     in
 
-    let rec aux first = 
+    let rec aux first =
 
       if (first > finish) then ()
       else
@@ -186,13 +190,13 @@ let pp ?(start=1) ?(finish=0) ppf t =
                 (Array.init rows (fun i -> Printf.sprintf "%d  " (i + 1)))
               ~col_labels:
                 (Array.init (min 5 (cols-first+1)) (fun i -> Printf.sprintf "-- %d --" (i + first) ))
-              ~print_right:false 
+              ~print_right:false
               ~print_foot:false
               () ) (lacpy ~ac:first ~n:(min 5 (cols-first+1)) (t.mo_coef)) ;
           Format.fprintf ppf "@]@;@;@]";
-          (aux [@tailcall]) (first+5) 
+          (aux [@tailcall]) (first+5)
         end
     in
-    aux start 
+    aux start
 
 

MOBasis/MOBasis.mli

@@ -55,6 +55,9 @@ val size : t -> int
 val mo_energies : t -> Vec.t
 (** Fock MO energies *)
 
+val values : t -> Coordinate.t -> Vec.t
+(** Values of the MOs evaluated at a given coordinate. *)
+
 (** {1 Creators} *)
 
 val make : simulation:Simulation.t ->