Unverified Commit ed9406d2 authored by Peter Doak's avatar Peter Doak Committed by GitHub
Browse files

Merge pull request #128 from CompFUSE/fix_ctaux_initialization 

Correct anitoperiodic interpolation for G0, vertex_pair.delta_r no longer  dropped
parents 0074aa1b 4657512b

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

+22 −12
Original line number Diff line number Diff line
@@ -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

+10 −6
Original line number Diff line number Diff line
@@ -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

+2 −2
Original line number Diff line number Diff line
@@ -17,8 +17,8 @@ 


    "domains": {

        "real-space-grids": {

            "cluster": [[2, 0],

                [0, 2]]

            "cluster": [[2, 2],

                [-4, 2]]

        }

    },