Merge pull request #128 from CompFUSE/fix_ctaux_initialization

Correct anitoperiodic interpolation for G0, vertex_pair.delta_r no longer  dropped

include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp

@@ -866,25 +866,35 @@ inline double TpEqualTimeAccumulator<parameters_type, MOMS_type>::interpolate_ak
   const static double beta = parameters.get_beta();
   const static double N_div_beta = parameters.get_sp_time_intervals() / beta;
 
-  // make sure that new_tau is positive !!
-  double new_tau = tau + beta;
+  int sign = 1;
+  // Map tau to [0, beta).
+  if (tau < 0) {
+    tau += beta;
+    sign = -1;
+  }
+  assert(0 <= tau && tau < beta);
 
-  double scaled_tau = new_tau * N_div_beta;
+  const double scaled_tau = tau * N_div_beta;
+  // Find interpolation index of on the left of tau.
+  const int t_ind = static_cast<int>(scaled_tau);
 
-  int t_ind = scaled_tau;
-  assert(shifted_t::get_elements()[t_ind] <= tau &&
-         tau < shifted_t::get_elements()[t_ind] + 1. / N_div_beta);
+#ifndef NDEBUG
+  const double* positive_times =
+      shifted_t::get_elements().data() + shifted_t::get_elements().size() / 2;
+  assert(positive_times[t_ind] <= tau && tau < positive_times[t_ind] + 1. / N_div_beta);
+#endif  // NDEBUG
 
-  double delta_tau = scaled_tau - t_ind;
-  assert(delta_tau > -1.e-16 && delta_tau <= 1 + 1.e-16);
+  const double delta_tau = scaled_tau - t_ind;
+  assert(delta_tau >= 0 && delta_tau < 1);
 
-  int linind = 4 * nu_nu_r_dmn_t_t_shifted_dmn(b_i, s_i, b_j, s_j, delta_r, t_ind);
+  const int linind = 4 * nu_nu_r_dmn_t_t_shifted_dmn(b_i, s_i, b_j, s_j, delta_r, t_ind);
 
-  double* a_ptr = &akima_coefficients(linind);
+  const double* a_ptr = &akima_coefficients(linind);
 
-  double result = (a_ptr[0] + delta_tau * (a_ptr[1] + delta_tau * (a_ptr[2] + delta_tau * a_ptr[3])));
+  const double result =
+      (a_ptr[0] + delta_tau * (a_ptr[1] + delta_tau * (a_ptr[2] + delta_tau * a_ptr[3])));
 
-  return result;
+  return sign * result;
 }
 
 template <class parameters_type, class MOMS_type>

include/dca/phys/dca_step/cluster_solver/ctaux/structs/read_write_config.hpp

@@ -24,16 +24,20 @@ namespace ctaux {
 
 template <class Parameters>
 io::Buffer& operator<<(io::Buffer& buff, const vertex_pair<Parameters>& v) {
-  return buff << v.bands << v.e_spins << v.spin_orbitals << v.r_sites << v.HS_spin << v.tau;
+  return buff << v.bands << v.e_spins << v.spin_orbitals << v.r_sites << v.delta_r << v.HS_spin
+              << v.tau;
 }
 
 template <class Parameters>
 io::Buffer& operator>>(io::Buffer& buff, vertex_pair<Parameters>& v) {
   v.creatable = false;
   v.annihilatable = true;
+  v.successfully_flipped = false;
+  v.Bennett = false;
   v.shuffled = true;
 
-  return buff >> v.bands >> v.e_spins >> v.spin_orbitals >> v.r_sites >> v.HS_spin >> v.tau;
+  return buff >> v.bands >> v.e_spins >> v.spin_orbitals >> v.r_sites >> v.delta_r >> v.HS_spin >>
+         v.tau;
 }
 
 template <class Parameters>
@@ -74,9 +78,9 @@ io::Buffer& operator>>(io::Buffer& buff, CT_AUX_HS_configuration<Parameters>& co
   return buff;
 }
 
-}  // ctaux
-}  // solver
-}  // phys
-}  // dca
+}  // namespace ctaux
+}  // namespace solver
+}  // namespace phys
+}  // namespace dca
 
 #endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTAUX_STRUCTS_READ_WRITE_CONFIG_HPP

test/unit/phys/dca_step/cluster_solver/ctaux/structs/input.json

@@ -17,8 +17,8 @@
 
     "domains": {
         "real-space-grids": {
-            "cluster": [[2, 0],
-                [0, 2]]
+            "cluster": [[2, 2],
+                [-4, 2]]
         }
     },