Debugging

Basis/GeneralBasis.ml

@@ -83,7 +83,7 @@ let combine basis_list =
   List.concat basis_list
   |> List.iter (fun (k,v) ->
       let l = Hashtbl.find_all h k in
-      Hashtbl.replace h k (Array.concat (v::l) )
+      Hashtbl.replace h k (Array.concat (l@[v]) )
       ) ;
   Hashtbl.fold (fun k v accu -> (k, v)::accu) h []
 

CI/CI.ml

@@ -676,7 +676,7 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
         Parallel.broadcast (lazy result)
       in
       let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in
-      (Conventions.rephase eigenvectors), eigenvalues
+      (Conventions.rephase @@ Util.remove_epsilons eigenvectors), eigenvalues
     in
 
     let eigensystem_direct () = 
@@ -702,7 +702,7 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
         Parallel.broadcast (lazy result)
       in
       let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in
-      (Conventions.rephase eigenvectors), eigenvalues
+      (Conventions.rephase @@ Util.remove_epsilons eigenvectors), eigenvalues
     in
 
     match algo with

CI/F12CI.ml

@@ -70,6 +70,7 @@ let single_matrices hf12_integrals density =
     done;
 *)
   
+(*---
 let hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj =
   let integrals = [
 
@@ -185,6 +186,104 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
 
   Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes)
 
+--- *)
+
+
+let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj =
+  let integrals = [
+
+    let { HF12.
+          simulation ; aux_basis ;
+          hf12 ; hf12_anti ;
+          hf12_single ; hf12_single_anti } = hf12_integrals
+    in
+
+    let kia = De.alfa ki and kib = De.beta ki in
+    let kja = De.alfa kj and kjb = De.beta kj in
+
+    let mo_a =
+      Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja)
+      |> Bitstring.to_list
+    and mo_b =
+      Bitstring.logand (Sp.bitstring kib) (Sp.bitstring kjb)
+      |> Bitstring.to_list
+    in
+
+    let one_e _ _ _ = 0. in
+
+    let two_e i j k l s s' =
+      if s = s' then
+	hf12_anti.{i,j,k,l} -. (
+	(List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +.
+	(List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) 
+        )
+      else
+	hf12.{i,j,k,l}  -. (
+	(List.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +.
+	(List.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) 
+        )
+    in
+
+    let h = MOBasis.ee_ints aux_basis in
+    let two_e_h i j k l s s' =
+      if s' <> s then
+        ERI.get_phys h l k j i
+      else
+        (ERI.get_phys h l k j i) -.  (ERI.get_phys h k l j i)
+    in
+
+    let f = MOBasis.f12_ints aux_basis in
+    let two_e_f i j k l s s' =
+      if s' <> s then
+        F12.get_phys f i j k l
+      else
+        (F12.get_phys f i j k l) -.  (F12.get_phys f i j l k)
+    in
+
+    let mo_of_s = function
+    | Spin.Alfa -> mo_a
+    | Spin.Beta -> mo_b
+    in
+
+    let three_e i j k l m n s s' s'' =
+         List.fold_left (fun accu a -> accu +. two_e_h i j l a s   s'  *. two_e_f a k m n s'  s'') 0. (mo_of_s s' )
+      -. List.fold_left (fun accu a -> accu +. two_e_h j i m a s'  s   *. two_e_f a k l n s   s'') 0. (mo_of_s s  )
+      +. List.fold_left (fun accu a -> accu +. two_e_h j k m a s'  s'' *. two_e_f a i n l s'' s  ) 0. (mo_of_s s'')
+      -. List.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s'  *. two_e_f a i m l s'  s  ) 0. (mo_of_s s' )
+      +. List.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s   *. two_e_f a j l m s   s' ) 0. (mo_of_s s  )
+      -. List.fold_left (fun accu a -> accu +. two_e_h i k l a s   s'' *. two_e_f a j n m s'' s' ) 0. (mo_of_s s'')
+    in
+
+    (one_e, two_e, Some three_e)
+  ]
+  in
+  CIMatrixElement.non_zero integrals deg_a deg_b ki kj
+  |> List.hd
+
+
+let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
+
+  if Parallel.master then
+    Printf.printf "Building matrix\n%!";
+  let det_space =
+    ci.CI.det_space
+  in
+
+  let m_HF =
+
+    let f =
+      match Ds.determinants det_space with
+      | Ds.Arbitrary _ -> CI.create_matrix_arbitrary
+      | Ds.Spin      _ -> CI.create_matrix_spin_computed ~nmax:3
+    in
+    f (fun deg_a deg_b ki kj ->
+         hf_ij_non_zero hf12_integrals deg_a deg_b ki kj
+      ) det_space
+
+  in
+
+  Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes)
+
 
 let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename ?(state=1) () =
 
@@ -210,10 +309,10 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
 
       
     let rec iteration ~state psi = 
-      (*
-      Format.printf "%a@." DeterminantSpace.pp_det_space @@ CI.det_space ci;
-      Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi;
-      *)
+(*
+Format.printf "%a@." DeterminantSpace.pp_det_space @@ CI.det_space ci;
+Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi;
+*)
       let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in
 
       let delta = 
@@ -284,7 +383,7 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
 
       in
       let eigenvectors = 
-         Conventions.rephase eigenvectors
+        Conventions.rephase @@ Util.remove_epsilons eigenvectors
       in
 
 

MOBasis/MOBasis.ml

@@ -49,6 +49,7 @@ let mo_energies t =
     let m_P = x_o_xt m_N m_C in
     match t.mo_type with
     | RHF  -> Fock.make_rhf  ~density:m_P (ao_basis t)
+    | Projected
     | ROHF -> (Mat.scal 0.5 m_P;
                Fock.make_uhf  ~density_same:m_P ~density_other:m_P (ao_basis t))
     | _ -> failwith "Not implemented"
@@ -139,6 +140,7 @@ let of_mo_basis simulation other =
     in
     (* Gram-Schmidt Orthonormalization *)
     gemm m_X @@ (Util.qr_ortho @@ gemm ~transa:`T m_X result)
+    |> Util.remove_epsilons
     |> Conventions.rephase
   in
 

SCF/Guess.ml

@@ -17,7 +17,6 @@ let hcore_guess ao_basis =
   and kin_ints = Ao.kin_ints ao_basis |> KinInt.matrix
   in
   Mat.add eN_ints kin_ints
-  |> Conventions.rephase
 
 
 let huckel_guess ao_basis =
@@ -47,7 +46,6 @@ let huckel_guess ao_basis =
         else
           diag.{i}
     )
-    |> Conventions.rephase
 
 
 let make ?(nocc=0) ~guess ao_basis = 

SCF/HartreeFock.ml

@@ -346,6 +346,7 @@ let make
     (* MOs in AO basis *)
     let m_C =
       gemm  m_X  m_C'
+      |> Util.remove_epsilons
       |> Conventions.rephase
     in
 
@@ -541,6 +542,7 @@ let make
     (* MOs in AO basis *)
     let m_C =
       gemm  m_X  m_C'
+      |> Util.remove_epsilons
       |> Conventions.rephase
     in
 
@@ -691,6 +693,8 @@ let pp_summary ppf t =
   ; print "Hartree-Fock energy" (energy t) 
   ; Format.fprintf ppf "@]"
   ; Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth)
+  ; Util.debug_matrix "MOs" (eigenvectors t)
+
 
 
 

Utils/FourIdxStorage.ml

@@ -344,7 +344,7 @@ let to_stream d =
     
 
 (** Write all integrals to a file with the <ij|kl> convention *)
-let to_file ?(cutoff=Constants.epsilon) ~filename data =
+let to_file ?(cutoff=Constants.integrals_cutoff) ~filename data =
   let oc = open_out filename in
   to_stream data 
   |> Stream.iter (fun {i_r1 ; j_r2 ; k_r1 ; l_r2 ; value} ->
@@ -496,7 +496,7 @@ let four_index_transform_dense_sparse ds coef source =
           List.iter (fun beta ->
               List.iter (fun alpha ->
                   let x = u.{alpha,beta,gamma} in
-                  if x <> 0. then
+                  if abs_float x > Constants.integrals_cutoff then
                     result := (alpha, beta, gamma, delta, x) :: !result;
                 ) (list_range 1 beta) 
             ) range_mo
@@ -541,7 +541,7 @@ let svd_of_dense t =
       scal (1. /. d.{1}) w;
       let rec aux i = 
         if i = Vec.dim w then i else
-        if w.{i} < 1.e-15 then i else
+        if w.{i} < epsilon_float then i else
           aux (i+1)
       in aux 1
     in
@@ -668,7 +668,7 @@ let four_index_transform_svd coef source =
       for c=1 to mo_num do
         for d=c to mo_num do
           let x = result.{a,b,c,d} in
-          if abs_float x > epsilon then
+          if abs_float x > Constants.integrals_cutoff then
             set_chem destination a b c d x;
         done
       done

Utils/Util.ml

@@ -302,6 +302,15 @@ let x_o_xt ~o ~x =
   gemm o x ~transb:`T 
   |> gemm x
 
+let remove_epsilons m =
+  let vecs =
+    Mat.to_col_vecs m
+  in
+  Array.map (fun v ->
+    let m = abs_float (2. *. amax v) in
+    Vec.map (fun x -> if abs_float x < m *. epsilon_float then 0. else x) v
+    ) vecs
+  |> Mat.of_col_vecs
 
 let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c =
   let d, u, _ = gesvd (lacpy overlap) in

Utils/Util.mli

@@ -130,6 +130,9 @@ val xt_o_x : o:Mat.t -> x:Mat.t -> Mat.t
 val x_o_xt : o:Mat.t -> x:Mat.t -> Mat.t
 (** Computes {% $\mathbf{X\, O\, X^\dag}$ %} *)
 
+val remove_epsilons : Mat.t -> Mat.t
+(** Remove values below machine precision *)
+
 val qr_ortho : Mat.t -> Mat.t
 (** QR orthogonalization of the input matrix *)