Rearrange G4. (28b648c7) · Commits · NDIP / Tool Sources / Direct-Geometry Spectroscopy / DCA / DCA Main

include/dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp

+8 −8

Original line number	Diff line number	Diff line
		@@ -178,22 +178,22 @@ void BseClusterSolver<ParametersType, DcaDataType, ScalarType>::load_G_II(
		G4.linind_2_subind(i, coor_2);

		// coordinate 0 1 2 3 4 5 6 7 8 9
		// G4 : b b b b k k k_ex w w w_ex
		// G4 : b b b b k w k w k_ex w_ex
		// G_II : b b k w b b k w

		coor_1[0] = coor_2[0];
		coor_1[1] = coor_2[1];
		coor_1[2] = coor_2[4]; // k_1
		coor_1[3] = coor_2[7]; // w_1
		coor_1[3] = coor_2[5]; // w_1
		coor_1[4] = coor_2[2];
		coor_1[5] = coor_2[3];
		coor_1[6] = coor_2[5]; // k_2
		coor_1[7] = coor_2[8]; // w_2
		coor_1[6] = coor_2[6]; // k_2
		coor_1[7] = coor_2[7]; // w_2

		// Note: only the first momentum and frequency exchange (7th and 10th index of G4) are analyzed.
		G_II(coor_1[0], coor_1[1], coor_1[2], coor_1[3], coor_1[4], coor_1[5], coor_1[6], coor_1[7]) =
		G4(coor_2[0], coor_2[1], coor_2[2], coor_2[3], coor_2[4], coor_2[5], 0, coor_2[7],
		coor_2[8], 0);
		G4(coor_2[0], coor_2[1], coor_2[2], coor_2[3], coor_2[4], coor_2[5], coor_2[6], coor_2[7],
		0, 0);
		}

		delete[] coor_1;

include/dca/phys/dca_data/dca_data.hpp

+2 −6

Original line number	Diff line number	Diff line
		@@ -95,12 +95,8 @@ public:
		func::function<std::complex<double>, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, WDmn>>;
		using TpGreensFunction =
		func::function<std::complex<TpAccumulatorScalar>,
		func::dmn_variadic<BDmn, BDmn, BDmn, BDmn, KClusterDmn, KClusterDmn,
		KExchangeDmn, WVertexDmn, WVertexDmn, WExchangeDmn>>;
		// The following typedef is for testing purposes
		using ReducedTpGreensFunction = func::function<
		std::complex<TpAccumulatorScalar>,
		func::dmn_variadic<BDmn, BDmn, BDmn, BDmn, KClusterDmn, KClusterDmn, WVertexDmn, WVertexDmn>>;
		func::dmn_variadic<BDmn, BDmn, BDmn, BDmn, KClusterDmn, WVertexDmn,
		KClusterDmn, WVertexDmn, KExchangeDmn, WExchangeDmn>>;

		DcaData(Parameters& parameters_ref);

include/dca/phys/dca_step/cluster_solver/exact_diagonalization_advanced/greens_functions/tp_greens_function.hpp

+25 −31

Original line number	Diff line number	Diff line
		@@ -21,6 +21,7 @@

		#include "dca/function/domains.hpp"
		#include "dca/function/function.hpp"
		#include "dca/phys/dca_data/dca_data.hpp"
		#include "dca/phys/dca_step/cluster_solver/exact_diagonalization_advanced/fermionic_overlap_matrices.hpp"
		#include "dca/phys/dca_step/cluster_solver/exact_diagonalization_advanced/fock_space.hpp"
		#include "dca/phys/dca_step/cluster_solver/exact_diagonalization_advanced/greens_functions/c_operator.hpp"
		@@ -91,6 +92,9 @@ public:
		using w_VERTEX = func::dmn_0<domains::vertex_frequency_domain<domains::COMPACT>>;
		using w_VERTEX_EXTENDED = func::dmn_0<domains::vertex_frequency_domain<domains::EXTENDED>>;

		using Data = DcaData<parameters_type>;
		using TpGreensFunctionData = typename Data::TpGreensFunction;

		public:
		TpGreensFunction(parameters_type& parameters_ref, fermionic_Hamiltonian_type& Hamiltonian_ref,
		fermionic_overlap_type& overlap_ref);
		@@ -100,14 +104,8 @@ public:

		void compute_two_particle_Greens_function(bool interacting);

		void compute_particle_particle_superconducting_A(
		func::function<complex_type,
		func::dmn_variadic<b_dmn, b_dmn, b_dmn, b_dmn, KClusterDmn, KClusterDmn,
		KExchangeDmn, w_VERTEX, w_VERTEX, WExchangeDmn>>& G4);
		void compute_particle_particle_superconducting_B(
		func::function<complex_type,
		func::dmn_variadic<b_dmn, b_dmn, b_dmn, b_dmn, KClusterDmn, KClusterDmn,
		KExchangeDmn, w_VERTEX, w_VERTEX, WExchangeDmn>>& G4);
		void compute_particle_particle_superconducting_A(TpGreensFunctionData& G4);
		void compute_particle_particle_superconducting_B(TpGreensFunctionData& G4);

		void compute_two_particle_Greens_function(
		func::function<complex_type, func::dmn_variadic<w_VERTEX_EXTENDED, w_VERTEX_EXTENDED, w_VERTEX_EXTENDED,
		@@ -182,10 +180,10 @@ private:
		func::function<vector_type, fermionic_Fock_dmn_type>& eigen_energies;
		func::function<matrix_type, fermionic_Fock_dmn_type>& eigen_states;

		func::function<int, func::dmn_variadic<fermionic_Fock_dmn_type, fermionic_Fock_dmn_type,
		b_s_r_dmn_type>>& creation_set_all;
		func::function<int, func::dmn_variadic<fermionic_Fock_dmn_type, fermionic_Fock_dmn_type,
		b_s_r_dmn_type>>& annihilation_set_all;
		func::function<int, func::dmn_variadic<fermionic_Fock_dmn_type, fermionic_Fock_dmn_type, b_s_r_dmn_type>>&
		creation_set_all;
		func::function<int, func::dmn_variadic<fermionic_Fock_dmn_type, fermionic_Fock_dmn_type, b_s_r_dmn_type>>&
		annihilation_set_all;

		func::function<int, KClusterDmn> min_k_dmn_t;
		func::function<int, KClusterDmn> q_plus_;
		@@ -326,9 +324,7 @@ void TpGreensFunction<parameters_type, ed_options>::write(Writer& writer) {
		*/
		template <typename parameters_type, typename ed_options>
		void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_superconducting_A(
		func::function<complex_type,
		func::dmn_variadic<b_dmn, b_dmn, b_dmn, b_dmn, KClusterDmn, KClusterDmn,
		KExchangeDmn, w_VERTEX, w_VERTEX, WExchangeDmn>>& G4) {
		TpGreensFunctionData& G4) {
		if (concurrency.id() == concurrency.first())
		std::cout << "\t" << __FUNCTION__ << std::endl;

		@@ -364,7 +360,7 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		// int w2 = min_w_vertex_ext(wn_ext);
		// int w3 = wm_ext;

		G4(b_0, b_1, b_2, b_3, 0, 0, 0, wn, wm, w_nu_idx) +=
		G4(b_0, b_1, b_2, b_3, 0, wn, 0, wm, w_nu_idx, 0, 0) +=
		G_tp_int(w1, w2, w3, b_0, 0, b_1, 1, b_2, 1, b_3, 0, 0, 0, 0);
		}
		}
		@@ -386,10 +382,10 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		for (int wn = 0; wn < w_VERTEX::dmn_size(); wn++) {
		std::cout << w_VERTEX::get_elements()[wn] << "\t\t";
		for (int wm = 0; wm < w_VERTEX::dmn_size(); wm++)
		if (std::abs(real(G4(0, 0, 0, 0, 0, 0, 0, wn, wm, w_nu_idx))) < 1.e-10)
		if (std::abs(real(G4(0, 0, 0, 0, 0, wn, 0, wm, w_nu_idx, 0, 0))) < 1.e-10)
		std::cout << 0.0 << "\t";
		else
		std::cout << real(G4(0, 0, 0, 0, 0, 0, 0, wn, wm, w_nu_idx)) << "\t";
		std::cout << real(G4(0, 0, 0, 0, 0, wn, 0, wm, w_nu_idx, 0, 0)) << "\t";
		std::cout << "\n";
		}
		std::cout << "\n";
		@@ -406,7 +402,7 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		if (std::abs(imag(G4(0, 0, 0, 0, 0, 0, 0, wn, wm, w_nu_idx))) < 1.e-10)
		std::cout << 0.0 << "\t";
		else
		std::cout << imag(G4(0, 0, 0, 0, 0, 0, 0, wn, wm, w_nu_idx)) << "\t";
		std::cout << imag(G4(0, 0, 0, 0, 0, wn, 0, wm, w_nu_idx, 0, 0)) << "\t";
		std::cout << "\n";
		}
		std::cout << "\n";
		@@ -417,7 +413,7 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		std::vector<double> x, y;
		for (int wn = 0; wn < w_VERTEX::dmn_size(); wn++) {
		x.push_back(w_VERTEX::get_elements()[wn]);
		y.push_back(real(G4(0, 0, 0, 0, 0, 0, 0, wn, wn, 0)));
		y.push_back(real(G4(0, 0, 0, 0, 0, wn, 0, wn, 0, 0)));
		}

		util::Plot::plotPoints(x, y);
		@@ -427,7 +423,7 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		std::vector<double> x, y;
		for (int wn = 0; wn < w_VERTEX::dmn_size(); wn++) {
		x.push_back(w_VERTEX::get_elements()[wn]);
		y.push_back(real(G4(0, 0, 0, 0, 0, 0, 0, wn, w_VERTEX::dmn_size() - 1 - wn, 0)));
		y.push_back(real(G4(0, 0, 0, 0, 0, wn, 0, w_VERTEX::dmn_size() - 1 - wn, 0, 0)));
		}

		util::Plot::plotPoints(x, y);
		@@ -454,9 +450,7 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		*/
		template <typename parameters_type, typename ed_options>
		void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_superconducting_B(
		func::function<complex_type,
		func::dmn_variadic<b_dmn, b_dmn, b_dmn, b_dmn, KClusterDmn, KClusterDmn,
		KExchangeDmn, w_VERTEX, w_VERTEX, WExchangeDmn>>& G4) {
		TpGreensFunctionData& G4) {
		if (concurrency.id() == concurrency.first())
		std::cout << "\n\n\t" << __FUNCTION__ << "\n\n";

		@@ -493,7 +487,7 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		// int w3 = wm_ext;

		// TODO check if ignoring the momentum is correct.
		G4(b_0, b_1, b_2, b_3, 0, 0, 0, wn, wm, w_nu_idx) +=
		G4(b_0, b_1, b_2, b_3, 0, wn, 0, wm, 0, w_nu_idx) +=
		G_tp_int(w1, w2, w3, b_0, 0, b_1, 1, b_2, 1, b_3, 0, 0, 0, 0);
		}
		}
		@@ -514,10 +508,10 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		for (int wn = 0; wn < w_VERTEX::dmn_size(); wn++) {
		std::cout << w_VERTEX::get_elements()[wn] << "\t\t";
		for (int wm = 0; wm < w_VERTEX::dmn_size(); wm++)
		if (abs(real(G4(0, 0, 0, 0, 0, 0, 0, wn, wm, w_nu_idx))) < 1.e-10)
		if (abs(real(G4(0, 0, 0, 0, wn, 0, wm, 0, w_nu_idx))) < 1.e-10)
		std::cout << 0.0 << "\t";
		else
		std::cout << real(G4(0, 0, 0, 0, 0, 0, 0, wn, wm, w_nu_idx)) << "\t";
		std::cout << real(G4(0, 0, 0, 0, 0, wn, 0, wm, 0, w_nu_idx)) << "\t";
		std::cout << "\n";
		}
		std::cout << "\n";
		@@ -531,10 +525,10 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		for (int wn = 0; wn < w_VERTEX::dmn_size(); wn++) {
		std::cout << w_VERTEX::get_elements()[wn] << "\t\t";
		for (int wm = 0; wm < w_VERTEX::dmn_size(); wm++)
		if (abs(imag(G4(0, 0, 0, 0, 0, 0, 0, wn, wm, w_nu_idx))) < 1.e-10)
		if (abs(imag(G4(0, 0, 0, 0, 0, wn, 0, wm, 0, w_nu_idx))) < 1.e-10)
		std::cout << 0.0 << "\t";
		else
		std::cout << imag(G4(0, 0, 0, 0, 0, 0, 0, wn, wm, w_nu_idx)) << "\t";
		std::cout << imag(G4(0, 0, 0, 0, 0, wn, 0, wm, 0, w_nu_idx)) << "\t";
		std::cout << "\n";
		}
		std::cout << "\n";
		@@ -545,7 +539,7 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		std::vector<double> x, y;
		for (int wn = 0; wn < w_VERTEX::dmn_size(); wn++) {
		x.push_back(w_VERTEX::get_elements()[wn]);
		y.push_back(real(G4(0, 0, 0, 0, 0, 0, 0, wn, wn, 0)));
		y.push_back(real(G4(0, 0, 0, 0, 0, wn, 0, wn, 0, 0)));
		}

		util::Plot::plotPoints(x, y);
		@@ -555,7 +549,7 @@ void TpGreensFunction<parameters_type, ed_options>::compute_particle_particle_su
		std::vector<double> x, y;
		for (int wn = 0; wn < w_VERTEX::dmn_size(); wn++) {
		x.push_back(w_VERTEX::get_elements()[wn]);
		y.push_back(real(G4(0, 0, 0, 0, 0, 0, wn, w_VERTEX::dmn_size() - 1 - wn)));
		y.push_back(real(G4(0, 0, 0, 0, 0, wn, 0, w_VERTEX::dmn_size() - 1 - wn, 0, 0)));
		}

		util::Plot::plotPoints(x, y);

include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator.hpp

+33 −34

Original line number	Diff line number	Diff line
		@@ -22,6 +22,7 @@
		#include "dca/linalg/matrix_view.hpp"
		#include "dca/linalg/matrixop.hpp"
		#include "dca/math/function_transform/special_transforms/space_transform_2D.hpp"
		#include "dca/phys/dca_data/dca_data.hpp"
		#include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/ndft/cached_ndft_cpu.hpp"
		#include "dca/phys/domains/cluster/momentum_exchange_domain.hpp"
		#include "dca/phys/domains/quantum/electron_band_domain.hpp"
		@@ -99,6 +100,7 @@ public:
		}

		protected:
		using Data = DcaData<Parameters>;
		using Profiler = typename Parameters::profiler_type;

		using WTpDmn = func::dmn_0<domains::vertex_frequency_domain<domains::COMPACT>>;
		@@ -111,9 +113,7 @@ protected:
		using SpGreenFunction =
		func::function<Complex, func::dmn_variadic<BDmn, BDmn, SDmn, KDmn, KDmn, WTpExtPosDmn, WTpExtDmn>>;

		using TpDomain =
		func::dmn_variadic<BDmn, BDmn, BDmn, BDmn, KDmn, KDmn, KExchangeDmn, WTpDmn, WTpDmn, WExchangeDmn>;
		using TpGreenFunction = func::function<Complex, TpDomain>;
		using TpGreenFunction = typename Data::TpGreensFunction;
		using Matrix = linalg::Matrix<Complex, linalg::CPU>;

		void initializeG0();
		@@ -402,13 +402,13 @@ double TpAccumulator<Parameters, linalg::CPU>::updateG4() {
		// (s1 == s2) G(k2 + k_ex, k1 + k_ex, s1) G(k1, k2, s1)>.
		for (int w_ex_idx = 0; w_ex_idx < exchange_frq.size(); ++w_ex_idx) {
		const int w_ex = exchange_frq[w_ex_idx];
		for (int w2 = 0; w2 < WTpDmn::dmn_size(); ++w2)
		for (int w1 = 0; w1 < WTpDmn::dmn_size(); ++w1)
		for (int k_ex_idx = 0; k_ex_idx < exchange_mom.size(); ++k_ex_idx) {
		const int k_ex = exchange_mom[k_ex_idx];
		for (int w2 = 0; w2 < WTpDmn::dmn_size(); ++w2)
		for (int k2 = 0; k2 < KDmn::dmn_size(); ++k2)
		for (int w1 = 0; w1 < WTpDmn::dmn_size(); ++w1)
		for (int k1 = 0; k1 < KDmn::dmn_size(); ++k1) {
		Complex* const G4_ptr = &(*G4_)(0, 0, 0, 0, k1, k2, k_ex_idx, w1, w2, w_ex_idx);
		Complex* const G4_ptr = &(*G4_)(0, 0, 0, 0, k1, w1, k2, w2, k_ex_idx, w_ex_idx);
		updateG4SpinDifference(G4_ptr, -1, k1, momentum_sum(k1, k_ex), w1,
		w_plus_w_ex(w1, w_ex), momentum_sum(k2, k_ex), k2,
		w_plus_w_ex(w2, w_ex), w2, sign_over_2, false);
		@@ -428,13 +428,13 @@ double TpAccumulator<Parameters, linalg::CPU>::updateG4() {
		// (s1 == s2) G(k2 + k_ex, k1 + k_ex, s1) G(k1, k2, s1)>.
		for (int w_ex_idx = 0; w_ex_idx < exchange_frq.size(); ++w_ex_idx) {
		const int w_ex = exchange_frq[w_ex_idx];
		for (int w2 = 0; w2 < WTpDmn::dmn_size(); ++w2)
		for (int w1 = 0; w1 < WTpDmn::dmn_size(); ++w1)
		for (int k_ex_idx = 0; k_ex_idx < exchange_mom.size(); ++k_ex_idx) {
		const int k_ex = exchange_mom[k_ex_idx];
		for (int w2 = 0; w2 < WTpDmn::dmn_size(); ++w2)
		for (int k2 = 0; k2 < KDmn::dmn_size(); ++k2)
		for (int w1 = 0; w1 < WTpDmn::dmn_size(); ++w1)
		for (int k1 = 0; k1 < KDmn::dmn_size(); ++k1) {
		Complex* const G4_ptr = &(*G4_)(0, 0, 0, 0, k1, k2, k_ex_idx, w1, w2, w_ex_idx);
		Complex* const G4_ptr = &(*G4_)(0, 0, 0, 0, k1, w1, k2, w2, k_ex_idx, w_ex_idx);
		updateG4SpinDifference(G4_ptr, 1, k1, momentum_sum(k1, k_ex), w1,
		w_plus_w_ex(w1, w_ex), momentum_sum(k2, k_ex), k2,
		w_plus_w_ex(w2, w_ex), w2, sign_over_2, false);
		@@ -454,13 +454,13 @@ double TpAccumulator<Parameters, linalg::CPU>::updateG4() {
		// = -1/2 G(k2 + k_ex, k1 + k_ex, s) G(k1, k2, -s).
		for (int w_ex_idx = 0; w_ex_idx < exchange_frq.size(); ++w_ex_idx) {
		const int w_ex = exchange_frq[w_ex_idx];
		for (int w2 = 0; w2 < WTpDmn::dmn_size(); ++w2)
		for (int w1 = 0; w1 < WTpDmn::dmn_size(); ++w1)
		for (int k_ex_idx = 0; k_ex_idx < exchange_mom.size(); ++k_ex_idx) {
		const int k_ex = exchange_mom[k_ex_idx];
		for (int w2 = 0; w2 < WTpDmn::dmn_size(); ++w2)
		for (int k2 = 0; k2 < KDmn::dmn_size(); ++k2)
		for (int w1 = 0; w1 < WTpDmn::dmn_size(); ++w1)
		for (int k1 = 0; k1 < KDmn::dmn_size(); ++k1) {
		Complex* const G4_ptr = &(*G4_)(0, 0, 0, 0, k1, k2, k_ex_idx, w1, w2, w_ex_idx);
		Complex* const G4_ptr = &(*G4_)(0, 0, 0, 0, k1, w1, k2, w2, k_ex_idx, w_ex_idx);
		for (int s = 0; s < 2; ++s)
		updateG4Atomic(G4_ptr, s, k1, k2, w1, w2, not s, momentum_sum(k2, k_ex),
		momentum_sum(k1, k_ex), w_plus_w_ex(w2, w_ex),
		@@ -477,13 +477,13 @@ double TpAccumulator<Parameters, linalg::CPU>::updateG4() {
		// = 1/2 G(k_ex-k2, k_ex-k1, s) G(k2, k1, -s).
		for (int w_ex_idx = 0; w_ex_idx < exchange_frq.size(); ++w_ex_idx) {
		const int w_ex = exchange_frq[w_ex_idx];
		for (int w2 = 0; w2 < WTpDmn::dmn_size(); ++w2)
		for (int w1 = 0; w1 < WTpDmn::dmn_size(); ++w1)
		for (int k_ex_idx = 0; k_ex_idx < exchange_mom.size(); ++k_ex_idx) {
		const int k_ex = exchange_mom[k_ex_idx];
		for (int w2 = 0; w2 < WTpDmn::dmn_size(); ++w2)
		for (int k2 = 0; k2 < KDmn::dmn_size(); ++k2)
		for (int w1 = 0; w1 < WTpDmn::dmn_size(); ++w1)
		for (int k1 = 0; k1 < KDmn::dmn_size(); ++k1) {
		Complex* const G4_ptr = &(*G4_)(0, 0, 0, 0, k1, k2, k_ex_idx, w1, w2, w_ex_idx);
		Complex* const G4_ptr = &(*G4_)(0, 0, 0, 0, k1, w1, k2, w2, k_ex_idx, w_ex_idx);
		for (int s = 0; s < 2; ++s)
		updateG4Atomic(G4_ptr, s, k1, k2, w1, w2, !s, q_minus_k(k1, k_ex),
		q_minus_k(k2, k_ex), w_ex_minus_w(w1, w_ex),
		@@ -548,8 +548,7 @@ void TpAccumulator<Parameters, linalg::CPU>::updateG4SpinDifference(
		// + sign * G(down,k1_a, k2_a, w1_a, w2_a)) *
		// (G(up,k1_b,k2_b,w1_b,w2_b) + sign * G(down,k1_b,k2_b,w1_b,w2_b))
		if (n_bands_ == 1) {
		G4_ptr += alpha
		(getGSingleband(0, k1_a, k2_a, w1_a, w2_a) +
		G4_ptr += alpha (getGSingleband(0, k1_a, k2_a, w1_a, w2_a) +
		Complex(sign) * getGSingleband(1, k1_a, k2_a, w1_a, w2_a)) *
		(getGSingleband(0, k1_b, k2_b, w1_b, w2_b) +
		Complex(sign) * getGSingleband(1, k1_b, k2_b, w1_b, w2_b));

include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu.hpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -249,7 +249,7 @@ void TpAccumulator<Parameters, linalg::GPU>::resetG4() {
		// Note: this method is not thread safe by itself.
		auto& G4 = get_G4();
		try {
		typename BaseClass::TpDomain tp_dmn;
		typename BaseClass::TpGreenFunction::this_domain_type tp_dmn;
		if (!multiple_accumulators_) {
		G4.setStream(streams_[0]);
		}