Rephasing

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
-      eigenvectors, eigenvalues
+      (Conventions.rephase 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
-      eigenvectors, eigenvalues
+      (Conventions.rephase eigenvectors), eigenvalues
     in
 
     match algo with

CI/F12CI.ml

@@ -190,6 +190,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;
+      *)
       let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in
 
       let delta = 
@@ -197,6 +201,10 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
          dressing_vector ~frozen_core hf12_integrals psi ci
          |> Matrix.to_mat
       in
+      (*
+      Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat delta;
+      *)
+
 
       Printf.printf "Cmax : %e\n" psi.{column_idx,state};
       Printf.printf "Norm : %e\n" (sqrt (gemm ~transa:`T delta delta).{state,state});
@@ -251,10 +259,13 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
 
 
         Parallel.broadcast (lazy (
-          Davidson.make ~threshold:1.e-9 ~guess:psi ~n_states:state diagonal matrix_prod
+          Davidson.make ~threshold:1.e-10 ~guess:psi ~n_states:state diagonal matrix_prod
           ))
 
       in
+      let eigenvectors = 
+         Conventions.rephase eigenvectors
+      in
 
 
       Vec.iter (fun energy -> Printf.printf "%g\t" energy) eigenvalues;

MOBasis/MOBasis.ml

@@ -41,6 +41,7 @@ let two_e_ints    t = Lazy.force t.ee_ints
 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 m_C = mo_coef t in
   let f = 
@@ -138,6 +139,7 @@ let of_mo_basis simulation other =
     in
     (* Gram-Schmidt Orthonormalization *)
     gemm m_X @@ (Util.qr_ortho @@ gemm ~transa:`T m_X result)
+    |> Conventions.rephase
   in
 
   let mo_occupation =
@@ -152,8 +154,6 @@ let of_mo_basis simulation other =
   
 
 
-
-
 let pp_mo ?(start=1) ?finish ppf t =
     let open Lacaml.Io in
     let rows = Mat.dim1 t.mo_coef

MOBasis/MOBasis.mli

@@ -55,7 +55,6 @@ val size : t -> int
 val mo_energies : t -> Vec.t
 (** Fock MO energies *)
 
-
 (** {1 Creators} *)
 
 val make : simulation:Simulation.t ->

SCF/Guess.ml

@@ -17,6 +17,7 @@ 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 =
@@ -46,6 +47,7 @@ let huckel_guess ao_basis =
         else
           diag.{i}
     )
+    |> Conventions.rephase
 
 
 let make ?(nocc=0) ~guess ao_basis = 
@@ -66,6 +68,7 @@ let test_case ao_basis =
         Lacaml.D.Mat.add 
           (AOBasis.eN_ints  ao_basis |> NucInt.matrix)
           (AOBasis.kin_ints ao_basis |> KinInt.matrix)
+        |> Conventions.rephase
         |> Lacaml.D.Mat.to_array
       in
       Array.iteri (fun i x ->

SCF/HartreeFock.ml

@@ -188,8 +188,6 @@ let exchange_energy t =
 
 
 
-
-
 let make
     ?kind
     ?guess:(guess=`Huckel)
@@ -348,6 +346,7 @@ let make
     (* MOs in AO basis *)
     let m_C =
       gemm  m_X  m_C'
+      |> Conventions.rephase
     in
 
     (* Hartree-Fock energy *)
@@ -542,6 +541,7 @@ let make
     (* MOs in AO basis *)
     let m_C =
       gemm  m_X  m_C'
+      |> Conventions.rephase
     in
 
     (* Hartree-Fock energy *)
@@ -672,25 +672,26 @@ let pp_summary ppf t =
   let ha_to_ev = Constants.ha_to_ev in
   let e = eigenvalues t in
 
-  Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth);
-  Format.fprintf ppf "@[<v>";
-  print "One-electron energy" (kin_energy t +. eN_energy t) ;
-  print "Kinetic"   (kin_energy t) ;
-  print "Potential" (eN_energy  t) ;
-  line () ;
-  print "Two-electron energy" (coulomb_energy t +. exchange_energy t) ;
-  print "Coulomb"  (coulomb_energy  t) ;
-  print "Exchange" (exchange_energy t) ;
-  line ();
-  print "HF HOMO" (ha_to_ev *. e.{nocc t});
-  print "HF LUMO" (ha_to_ev *. e.{nocc t + 1});
-  print "HF LUMO-LUMO" (ha_to_ev *. (e.{nocc t + 1} -. e.{nocc t }));
-  line ();
-  print "Electronic   energy" (energy t -. nuclear_repulsion t) ;
-  print "Nuclear   repulsion" (nuclear_repulsion t) ;
-  print "Hartree-Fock energy" (energy t) ;
-  Format.fprintf ppf "@]";
   Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth)
+  ; Format.fprintf ppf "@[<v>"
+  ; print "One-electron energy" (kin_energy t +. eN_energy t) 
+  ; print "Kinetic"   (kin_energy t) 
+  ; print "Potential" (eN_energy  t) 
+  ; line () 
+  ; print "Two-electron energy" (coulomb_energy t +. exchange_energy t) 
+  ; print "Coulomb"  (coulomb_energy  t) 
+  ; print "Exchange" (exchange_energy t) 
+  ; line ()
+  ; print "HF HOMO" (ha_to_ev *. e.{nocc t})
+  ; print "HF LUMO" (ha_to_ev *. e.{nocc t + 1})
+  ; print "HF LUMO-LUMO" (ha_to_ev *. (e.{nocc t + 1} -. e.{nocc t }))
+  ; line ()
+  ; print "Electronic   energy" (energy t -. nuclear_repulsion t) 
+  ; print "Nuclear   repulsion" (nuclear_repulsion t) 
+  ; print "Hartree-Fock energy" (energy t) 
+  ; Format.fprintf ppf "@]"
+  ; Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth)
+
 
 
 let pp_hf ppf t =

Utils/Davidson.ml

@@ -136,11 +136,14 @@ let make
       Mat.of_col_vecs_list u_proposed
     in
     let maxu = lange u_proposed ~norm:`M in
-    let thr = maxu *. 0.001 in
+    let thr = maxu *. 0.00001 in
     let u_proposed = 
       Mat.map (fun x -> if abs_float x < thr then 0. else x) u_proposed
       |> Mat.to_col_vecs_list
     in
+    (*
+    Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat m_new_U;
+    *)
 
 
     if residual_norm > threshold && macro > 0 then

run_hartree_fock.ml

@@ -71,7 +71,7 @@ let () =
 
   Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf
   (*
-  let mos = MOBasis.of_hartree_fock hf in
+  ; let mos = MOBasis.of_hartree_fock hf in
   Format.fprintf ppf "@[%a@]@." (fun ppf x -> MOBasis.pp_mo ppf x) mos
   *)