From 4e1a6ccdf20793f8d80c701646f4d175ac3d3aca Mon Sep 17 00:00:00 2001 From: Peter Doak <11856692+PDoakORNL@users.noreply.github.com> Date: Fri, 29 May 2026 17:33:12 +0000 Subject: [PATCH 1/7] Fix compute_tc.py HDF5 path for updated analysis output format The analysis application now writes leading-eigenvalues as a dataset directly under /analysis-functions/ (compound type with 'r'/'i' fields), rather than as a group containing a 'data' sub-dataset. Update the script to read the dataset directly and extract the real part via the named 'r' field instead of positional indexing. --- tutorials/tc/compute_tc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tutorials/tc/compute_tc.py b/tutorials/tc/compute_tc.py index e3f5bac84..80778720c 100644 --- a/tutorials/tc/compute_tc.py +++ b/tutorials/tc/compute_tc.py @@ -58,8 +58,8 @@ def computeTransitionTemp(T, eigval): # Store real part of leading eigenvalue (imaginary part = 0). # Eigenvalues are sorted w.r.t. size in decreasing order. - leading_eigenvalues = data['analysis-functions']['leading-eigenvalues']['data'][:] - eigval.append(leading_eigenvalues[0][0]) + leading_eigenvalues = data['analysis-functions']['leading-eigenvalues'][:] + eigval.append(leading_eigenvalues[0]['r']) data.close() From 74bb51d6db557f003a4d1dd7abc713e8db789b74 Mon Sep 17 00:00:00 2001 From: Peter Doak <11856692+PDoakORNL@users.noreply.github.com> Date: Fri, 29 May 2026 17:51:09 +0000 Subject: [PATCH 2/7] Fix compute_tc.py for h5py complex128 auto-conversion h5py automatically converts HDF5 compound types with 'r'/'i' fields into numpy complex128 arrays. The script was trying string-index field access (['r']) which fails on complex scalars. Use .real instead. --- tutorials/tc/compute_tc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/tc/compute_tc.py b/tutorials/tc/compute_tc.py index 80778720c..6584bf891 100644 --- a/tutorials/tc/compute_tc.py +++ b/tutorials/tc/compute_tc.py @@ -59,7 +59,7 @@ def computeTransitionTemp(T, eigval): # Store real part of leading eigenvalue (imaginary part = 0). # Eigenvalues are sorted w.r.t. size in decreasing order. leading_eigenvalues = data['analysis-functions']['leading-eigenvalues'][:] - eigval.append(leading_eigenvalues[0]['r']) + eigval.append(leading_eigenvalues[0].real) data.close() From e35845415ccc6f324e6469a0f1abcd397dcacdcf Mon Sep 17 00:00:00 2001 From: Peter Doak <11856692+PDoakORNL@users.noreply.github.com> Date: Sat, 30 May 2026 00:07:33 +0000 Subject: [PATCH 3/7] Fix C++ BSE analysis to match Python reference Fixes applied: - basis_transform.hpp: add is_initialized() to func::dmn_0 specialization, enabling per-patch interpolation matrix recomputation - interpolation_matrices.hpp: invalidate trafo cache each patch iteration - bse_cluster_solver.hpp: apply beta*Nc scaling to Gamma_cluster; load G_II_0 before symmetrization - bse_lattice_solver.hpp: remove spurious 1/(beta*N_host) renorm factor; skip PP_UP_DOWN symmetrization - dca_data.hpp: read chemical potential from DCA-loop-functions HDF5 group - analyzeDCA_multiOrbit.py: fix typos, add __main__ block --- .../basis_transform/basis_transform.hpp | 4 + .../bse_solver/bse_cluster_solver.hpp | 14 +- .../bse_solver/bse_lattice_solver.hpp | 7 +- include/dca/phys/dca_data/dca_data.hpp | 13 + .../coarsegraining/interpolation_matrices.hpp | 1 + tutorials/tc/analyzeDCA_multiOrbit.py | 733 ++++++++++++++++++ 6 files changed, 764 insertions(+), 8 deletions(-) create mode 100644 tutorials/tc/analyzeDCA_multiOrbit.py diff --git a/include/dca/math/function_transform/basis_transform/basis_transform.hpp b/include/dca/math/function_transform/basis_transform/basis_transform.hpp index 5198e102a..d1da9261e 100644 --- a/include/dca/math/function_transform/basis_transform/basis_transform.hpp +++ b/include/dca/math/function_transform/basis_transform/basis_transform.hpp @@ -62,6 +62,10 @@ class basis_transform, func::dmn_0> { typedef typename basis_transformation_type::matrix_type matrix_type; public: + static bool& is_initialized() { + return basis_transformation_type::is_initialized(); + } + static std::string& get_name() { return basis_transformation_type::get_name(); } diff --git a/include/dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp b/include/dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp index 952381e6d..ab84e8214 100644 --- a/include/dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp +++ b/include/dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp @@ -113,13 +113,14 @@ void BseClusterSolver::compute_Gamma_cl func::function, DCA_matrix_dmn_t> G_II("G_II"); func::function, DCA_matrix_dmn_t> G_II_0("G_II_0"); - apply_symmetries_sp(); - + // Load G_II and compute G_II_0 BEFORE symmetrizing Sigma/G_k_w, + // so the bare bubble matches Python's unsymmetrized calculation. load_G_II(G_II); - load_G_II_0(G_II_0); load_G_II_0_function(G_II_0); + apply_symmetries_sp(); + apply_symmetries_tp(G_II, G_II_0); solve_BSE_on_cluster(G_II, G_II_0); @@ -144,7 +145,8 @@ void BseClusterSolver::apply_symmetries if (concurrency.id() == concurrency.first()) std::cout << "\t" << __FUNCTION__ << "\n\n"; - if (parameters.symmetrize_Gamma()) { + if (parameters.symmetrize_Gamma() && + parameters.get_four_point_channels()[0] != FourPointType::PARTICLE_PARTICLE_UP_DOWN) { if (true) { if (concurrency.id() == concurrency.first()) std::cout << "symmetrize Gamma_lattice according to the symmetry-group \n" << std::endl; @@ -327,9 +329,11 @@ void BseClusterSolver::solve_BSE_on_clu dca::linalg::matrixop::copyArrayToMatrix(N, N, &G_II_0(0), N, G4_0_inv); dca::linalg::matrixop::inverse(G4_0_inv); + const ScalarType beta_Nc_factor = parameters.get_beta() * static_cast(k_DCA::dmn_size()); for (int j = 0; j < N; j++) for (int i = 0; i < N; i++) - Gamma_cluster(i, j) = G4_0_inv(i, j) - G4_inv(i, j); + Gamma_cluster(i, j) = (G4_0_inv(i, j) - G4_inv(i, j)) * beta_Nc_factor; + } } // namespace analysis diff --git a/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp b/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp index a06f0490f..faf9605a6 100644 --- a/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp +++ b/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp @@ -353,8 +353,8 @@ void BseLatticeSolver::computeChi0Lattice() coarsegraining_tp.execute(MOMS.H_HOST, MOMS.Sigma, chi_0_lattice); } - // Renormalize and set diagonal \chi_0 matrix. - const ScalarType renorm = 1. / (parameters.get_beta() * k_HOST_VERTEX::dmn_size()); + // Set diagonal \chi_0 matrix. + const ScalarType renorm = 1.; for (int w_ind = 0; w_ind < WVertexDmn::dmn_size(); w_ind++) for (int K_ind = 0; K_ind < k_HOST_VERTEX::dmn_size(); K_ind++) @@ -389,7 +389,8 @@ void BseLatticeSolver::computeGammaLattice( Gamma_lattice(i) = Gamma_cluster(i); } - if (parameters.symmetrize_Gamma()) { + if (parameters.symmetrize_Gamma() && + parameters.get_four_point_channels()[0] != FourPointType::PARTICLE_PARTICLE_UP_DOWN) { if (concurrency.id() == concurrency.first()) std::cout << "Symmetrize Gamma_lattice according to the symmetry group." << std::endl; Symmetrize::execute(Gamma_lattice, parameters.get_four_point_momentum_transfer()); diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index d1fafc22b..19428f70a 100644 --- a/include/dca/phys/dca_data/dca_data.hpp +++ b/include/dca/phys/dca_data/dca_data.hpp @@ -396,6 +396,19 @@ void DcaData::read(dca::io::Reader chemical_potentials; + bool cp_present = reader.execute("chemical-potential", chemical_potentials); + if (cp_present && !chemical_potentials.empty()) { + parameters_.get_chemical_potential() = chemical_potentials.back(); + if (concurrency_.id() == concurrency_.first()) + std::cout << "\tAdjusted chemical potential from DCA loop: " + << parameters_.get_chemical_potential() << std::endl; + } + reader.close_group(); + } + reader.open_group("functions"); reader.execute(Sigma); diff --git a/include/dca/phys/dca_step/cluster_mapping/coarsegraining/interpolation_matrices.hpp b/include/dca/phys/dca_step/cluster_mapping/coarsegraining/interpolation_matrices.hpp index 54ecc3889..ac1329d12 100644 --- a/include/dca/phys/dca_step/cluster_mapping/coarsegraining/interpolation_matrices.hpp +++ b/include/dca/phys/dca_step/cluster_mapping/coarsegraining/interpolation_matrices.hpp @@ -113,6 +113,7 @@ auto interpolation_matrices----------->---- + # l2, K l4, K' + # + self.addPhaseFactors() + + f.close() + + def addPhaseFactors(self): + if self.readG4 & (self.channel == "PARTICLE_HOLE_MAGNETIC"): + print("Adding phase factors to G4") + G4 = self.unfoldG4(self.G4) + for iq in range(self.Nq): + Q = self.Kvecs[self.Qindices[iq]] + for ik1 in range(self.Nc): + K1 = self.Kvecs[ik1] + iK1pQ = self.iKSum[ik1, self.Qindices[iq]] # index of K1+Q mapped to 1.BZ + for ik2 in range(self.Nc): + K2 = self.Kvecs[ik2] + iK2pQ = self.iKSum[ik2, self.Qindices[iq]] # index of K2+Q mapped to 1.BZ + for l1 in range(self.nOrb): + for l3 in range(self.nOrb): + pf1 = np.exp(-1j*np.dot(K1+Q-self.Kvecs[iK1pQ, :], self.aVecs[l1, :])) + pf2 = np.exp(1j*np.dot(K2+Q-self.Kvecs[iK2pQ, :], self.aVecs[l3, :])) + G4[:, iq, :, ik1, l1, :, :, ik2, l3, :] *= pf1 * pf2 + + def dispersion(self, kx, ky): + ek = np.zeros((self.nOrb, self.nOrb), dtype='complex') + if self.model == "square": # 1-orbital model + ek[0, 0] = -2.*self.t*(cos(kx)+cos(ky)) - 4.0*self.tp*cos(kx)*cos(ky) + elif self.model == "triangular": # 1-orbital model + ek[0, 0] = -2. * self.t * cos(kx) - 4. * self.t * cos(sqrt(3.) * ky / 2.) * cos(kx / 2.) + elif self.model == "Kagome": # 3-orbital model + ek[0, 1] = -2.0 * self.t * np.cos(0.5 * kx) + ek[0, 2] = -2.0 * self.t * np.cos(0.25 * kx + 0.25 * np.sqrt(3.0) * ky) + ek[1, 2] = -2.0 * self.t * np.cos(0.25 * kx - 0.25 * np.sqrt(3.0) * ky) + + ek[1, 0] = -2.0 * self.t * np.cos(0.5 * kx) + ek[2, 0] = -2.0 * self.t * np.cos(0.25 * kx + 0.25 * np.sqrt(3.0) * ky) + ek[2, 1] = -2.0 * self.t * np.cos(0.25 * kx - 0.25 * np.sqrt(3.0) * ky) + elif self.model == "La3Ni2O7": # two-orbital bilayer model for La3Ni2O7 + val11 = -2. * self.t11 * (cos(kx) + cos(ky)) + val22 = -2. * self.t22 * (cos(kx) + cos(ky)) + val12 = 2. * self.t12 * (cos(kx) - cos(ky)) + + ek[0, 0] = val11 + self.Delta + ek[1, 1] = val22 + ek[2, 2] = val11 + self.Delta + ek[3, 3] = val22 + + ek[0, 1] = val12 + ek[1, 0] = val12 + ek[2, 3] = val12 + ek[3, 2] = val12 + + ek[0, 2] = -self.tperp11 + ek[2, 0] = -self.tperp11 + ek[1, 3] = -self.tperp22 + ek[3, 1] = -self.tperp22 + + elif self.model == "two-orbital": # generic two-orbital model + val11 = -2. * self.t11 * (cos(kx) + cos(ky)) + val22 = -2. * self.t22 * (cos(kx) + cos(ky)) + val12 = -4. * self.t12 * sin(kx) * sin(ky) + + ek[0, 0] = val11 + ek[1, 1] = val22 + + ek[0, 1] = val12 + ek[1, 0] = val12 + + elif self.model == "bilayer-Hubbard-model": # generic two-orbital model + val = -2. * self.t * (cos(kx) + cos(ky)) - 4. * self.tprime * cos(kx)*cos(ky) + + ek[0, 0] = val + ek[1, 1] = val + + ek[0, 1] = -self.tperp + ek[1, 0] = -self.tperp + + else: + print("Model not implemented") + sys.exit(0) + return ek + + def calcChi0Cluster(self): + print("Now calculating chi0 on cluster") + NwG = self.NwG + + self.chic0 = np.zeros_like(self.G4) + # Load frequency and K domain + wnSet = self.wnSet + # Kset = self.shift21BZ(self.Kvecs) + Kset = self.Kvecs + + for iQ in range(self.G4.shape[0]): + (iwm, iq) = np.unravel_index(iQ, (self.Nwm, self.Nq)) + + for iwn, wn in enumerate(wnSet): # reduced tp frequencies ! + for iK, K in enumerate(Kset): # cluster K + if "PARTICLE_HOLE" in self.channel: + iKG2 = self.iKSum[iK, self.Qindices[iq]] # k+Q + iwG = int(iwn - self.iwG40 + self.iwG0) + iwPlusiwm = self.iwnPlusiwm(iwG, iwm, maxIndex=NwG) # iwn+iwm + G1 = self.Green[iwG, iK, :, :] + G2 = self.Green[iwPlusiwm, iKG2, :, :] + cc = np.outer(G2.reshape(self.nOrb*self.nOrb), G1.reshape(self.nOrb*self.nOrb)) + + for ind1 in range(self.nOrb*self.nOrb): + (l1, l2) = np.unravel_index(ind1, (self.nOrb, self.nOrb)) + iK1 = np.ravel_multi_index((iwn, iK, l1, l2), dims=(self.NwG4, self.Nc, self.nOrb, self.nOrb)) + for ind2 in range(self.nOrb*self.nOrb): + (l3, l4) = np.unravel_index(ind2, (self.nOrb, self.nOrb)) + iK2 = np.ravel_multi_index((iwn, iK, l3, l4), dims=(self.NwG4, self.Nc, self.nOrb, self.nOrb)) + ccind1 = np.ravel_multi_index((l3, l1), dims=(self.nOrb, self.nOrb)) + ccind2 = np.ravel_multi_index((l2, l4), dims=(self.nOrb, self.nOrb)) + # Add phase factor exp((K+Q-K[iK2])*r_l1-r_l2) + Q = self.Kvecs[self.Qindices[iq]] + pf = np.exp(1j*np.dot(K+Q-Kset[iKG2], self.aVecs[l3, :]-self.aVecs[l1, :])) + self.chic0[iQ, iK1, iK2] = -cc[ccind1, ccind2] * pf + + # Note that the disconnected vaccum like term has already been subtracted for the charge channel + # for the particle-hole-magnetic channel, it is in principle there for q=0 and wm=0, but the + # sum over the spins with a factor of sigma*sigma' cancels it for spin symmetric models + # For Rashba-Hubbard and Moire-Hubbard models, this term is finite since the spin degeneracy is lifted + # We therefore treat those models with a different analysis script. + + elif "PARTICLE_PARTICLE" in self.channel: + iKG2 = self.iKSum[self.iKDiff[self.iK0, iK], self.Qindices[iq]] # -k+Q + iwG = int(iwn - self.iwG40 + self.iwG0) + miwPlusiwm = self.iwnPlusiwm(NwG-iwG-1, iwm, maxIndex=NwG) # -iwn+iwm + G1 = self.Green[iwG, iK, :, :] + G2 = self.Green[miwPlusiwm, iKG2, :, :] + cc = np.outer(G1.reshape(self.nOrb*self.nOrb), G2.reshape(self.nOrb*self.nOrb)) + + for ind1 in range(self.nOrb*self.nOrb): + (l1, l2) = np.unravel_index(ind1, (self.nOrb, self.nOrb)) + iK1 = np.ravel_multi_index((iwn, iK, l1, l2), dims=(self.NwG4, self.Nc, self.nOrb, self.nOrb)) + for ind2 in range(self.nOrb*self.nOrb): + (l3, l4) = np.unravel_index(ind2, (self.nOrb, self.nOrb)) + iK2 = np.ravel_multi_index((iwn, iK, l3, l4), dims=(self.NwG4, self.Nc, self.nOrb, self.nOrb)) + + ccind1 = np.ravel_multi_index((l1, l3), dims=(self.nOrb, self.nOrb)) + ccind2 = np.ravel_multi_index((l2, l4), dims=(self.nOrb, self.nOrb)) + # ccind1 = np.ravel_multi_index((l2, l4), dims=(self.nOrb, self.nOrb)) + # ccind2 = np.ravel_multi_index((l3, l1), dims=(self.nOrb, self.nOrb)) + + # Add phase factor exp((K+Q-K[iK2])*r_l1-r_l2) + Q = self.Kvecs[self.Qindices[iq]] + pf = np.exp(1j*np.dot(K+Q-Kset[iKG2], self.aVecs[l3, :]-self.aVecs[l1, :])) + self.chic0[iQ, iK1, iK2] = cc[ccind1, ccind2] * pf + else: + print("Channel not implemented! Bailing out") + exit() + + # def calcGammaIrr(self): + # # Calculate the irr. vertex GammaIrr + # Nq=self.Nq; Nwm=self.Nwm; Nc=self.Nc; NwG4=self.NwG4; NwG=self.NwG; nt = self.nt + # self.GammaM = np.zeros((Nwm, Nq, nt, nt), dtype='complex') + # for iwm in range(Nwm): + # for iq in range(Nq): + # G4 = self.G4[iwm,iq,...].reshape(nt,nt) + # chic0 = np.diag(self.chic0[iwm,iq,...].reshape(nt)) + # G4 = linalg.inv(G4) + # chic0 = linalg.inv(chic0) + # self.GammaM[iwm, iq, ...] = chic0 - G4 + # self.GammaM[iwm, iq, ...] *= float(Nc)*self.invT + # self.Gamma = self.GammaM.reshape(Nwm, Nq, NwG4, Nc, NwG4, Nc) + + + def calcChiFromG4(self, G4, form_factor=None): + # G4 data layout: G4(Q, K, K') + # Without form factor: Sum over K, K' + if form_factor is None: + form_factor = self.swave + G4 = self.unfoldG4(G4) + + susK = np.zeros((self.Nq, self.Nwm, self.Nc, self.nOrb, self.nOrb, self.Nc, self.nOrb, self.nOrb)) + sus = np.zeros((self.Nq, self.Nwm, self.nOrb, self.nOrb, self.nOrb, self.nOrb)) + for iwm in range(self.Nwm): + for iQ in range(self.Nq): + susK[iQ, iwm, ...] = np.sum(G4[iwm, iQ, :, :, :, :, :, :, :, :].real, axis=(0, 4)) + + for iK1 in range(self.Nc): + gk1 = form_factor(self.Kvecs[iK1, 0], self.Kvecs[iK1, 1]) + for iK2 in range(self.Nc): + gk2 = form_factor(self.Kvecs[iK2, 0], self.Kvecs[iK2, 1]) + sus[:, :, ...] += gk1 * susK[:, :, iK1, :,:, iK2, :,:] * gk2 + + sus /= self.invT * self.Nc + # print("Susceptibility from G4: ", sus) + return sus[0,0, ...], form_factor + + def calcChiFromG4_Test(self, G4): + # G4 data layout: G4(Q, K, K') + # Without form factor: Sum over K, K' + G4 = self.unfoldG4(G4)[0, 0, ...] + + susK = np.zeros((self.Nc, self.nOrb, self.nOrb, self.Nc, self.nOrb, self.nOrb)) + sus = 0.0 + susD = 0.0 + susDM = np.zeros((self.nOrb, self.nOrb)) + susK = np.sum(G4[:, :, :, :, :, :, :, :].real, axis=(0, 4)) + + gks = self.swave(self.Kvecs[:, 0], self.Kvecs[:, 1]) + gkxs = self.xswave(self.Kvecs[:, 0], self.Kvecs[:, 1]) + gkspm = self.spmwave(self.Kvecs[:, 0], self.Kvecs[:, 1]) + gkd = self.dwave(self.Kvecs[:, 0], self.Kvecs[:, 1]) + + for iK1 in range(self.Nc): + for iK2 in range(self.Nc): + for l1 in range(self.nOrb): + for l2 in range(self.nOrb): + for l3 in range(self.nOrb): + for l4 in range(self.nOrb): + if (l1//2==l1/2) == (l2//2==l2/2): # both even or both odd --> same orbital + gk1 = gks[iK1]+gkxs[iK1]+gkspm[iK1] + else: + gk1 = gkd[iK1] + if (l3//2==l3/2) == (l4//2==l4/2): # both even or both odd --> same orbital + gk2 = gks[iK2]+gkxs[iK2]+gkspm[iK2] + else: + gk2 = gkd[iK2] + sus += gk1 * susK[iK1, l1, l2, iK2, l3, l4] * gk2 + if (l3==l1) & (l4==l2): + susD += gk1 * susK[iK1, l1, l2, iK2, l1, l2] * gk2 + susDM[l1, l2] += gk1 * susK[iK1, l1, l2, iK2, l1, l2] * gk2 + return susDM + + sus /= self.invT * self.Nc + print("Susceptibility from G4 Test: ", sus, susD) + # return sus[0,0, ...], form_factor + + def swave(self, kx, ky): + return np.ones_like(self.Kvecs[:,0]) + + def xswave(self, kx, ky): + return np.cos(kx) + np.cos(ky) + + def spmwave(self, kx, ky): + return np.cos(kx) * np.cos(ky) + + def dwave(self, kx, ky): + return np.cos(kx) - np.cos(ky) + + def buildFullChi0Lattice(self, nkfine=128): + Nq = self.Nq + Nwm = self.Nwm + Nc = self.Nc + NwG4 = self.NwG4 + nOrb = self.nOrb + nt = NwG4*Nc*nOrb*nOrb + self.chi0 = np.zeros_like(self.G4) + + kPatch = self.build_kGrid(cluster=True, nkfine=nkfine) + + for iQ in range(self.G4.shape[0]): + self.chi0[iQ, :, :] = self.buildChi0Lattice(iQ, kPatch) + + def buildChi0Lattice(self, iQ, kPatch): + (iwm, iq) = np.unravel_index(iQ, (self.Nwm, self.Nq)) + Qvec = self.Kvecs[self.Qindices[iq], :] + print("Now calculating chi0 on lattice for Q=", Qvec," and iwm=", iwm) + + NwG = self.NwG + + # Load frequency and K domain + wnSet = self.wnSet + # Kset = self.shift21BZ(self.Kvecs) + Kset = self.Kvecs + + # Now coarse-grain G*G to build chi0(K) = Nc/N sum_k Gc(K+k')Gc(-K-k') + nOrb = self.nOrb + chi0 = np.zeros_like(self.G4[iQ, :, :]) + + # Pre-calculate energies + ek1 = np.zeros((self.Nc, kPatch.shape[0], self.nOrb, self.nOrb), dtype='complex') + ek2 = np.zeros((self.Nc, kPatch.shape[0], self.nOrb, self.nOrb), dtype='complex') + for iK, K in enumerate(Kset): # cluster K + for ik, k in enumerate(kPatch): + kx = K[0] + k[0] + ky = K[1] + k[1] + if "PARTICLE_HOLE" in self.channel: + ek1[iK, ik, ...] = self.dispersion(kx, ky) + ek2[iK, ik, ...] = self.dispersion(kx + Qvec[0], ky + Qvec[1]) # e(k+Q) + elif "PARTICLE_PARTICLE" in self.channel: + ek1[iK, ik, ...] = self.dispersion(kx, ky) + ek2[iK, ik, ...] = self.dispersion(-kx + Qvec[0], -ky + Qvec[1]) # e(-k+Q) + + self.cG1 = np.zeros((self.NwG4, self.Nc, self.nOrb, self.nOrb), dtype='complex') + self.cG2 = np.zeros((self.NwG4, self.Nc, self.nOrb, self.nOrb), dtype='complex') + for iwn, wn in enumerate(wnSet): # reduced tp frequencies ! + for iK, K in enumerate(Kset): # cluster K + cc = np.zeros((self.nOrb*self.nOrb, self.nOrb*self.nOrb), dtype='complex') + cG1 = np.zeros((self.nOrb, self.nOrb), dtype='complex') + cG2 = np.zeros((self.nOrb, self.nOrb), dtype='complex') + + if "PARTICLE_HOLE" in self.channel: + iK2 = self.iKSum[iK, self.Qindices[iq]] # k+Q + for ik, k in enumerate(kPatch): + e1 = ek1[iK, ik, ...] + e2 = ek2[iK, ik, ...] + iwG = int(iwn - self.iwG40 + self.iwG0) + # minusiw = self.minusiwn(iwG, maxIndex=NwG) # -iwn + iwPlusiwm = self.iwnPlusiwm(iwG, iwm, maxIndex=NwG) # iwn+iwm + sigmaK1 = self.sigma[iwG, iK, :, :] + sigmaK2 = self.sigma[iwPlusiwm, iK2, :, :] + G1 = np.linalg.inv((1j*wn + self.mu)*np.identity(nOrb) - e1 - sigmaK1) + G2 = np.linalg.inv((1j*self.wnG[iwPlusiwm] + self.mu)*np.identity(nOrb) - e2 - sigmaK2) + cc += -np.outer(G2.reshape(self.nOrb*self.nOrb), G1.reshape(self.nOrb*self.nOrb)) + cG1 += G1 + cG2 += G2 + + for ind1 in range(self.nOrb*self.nOrb): + (l1, l2) = np.unravel_index(ind1, (self.nOrb, self.nOrb)) + iK1 = np.ravel_multi_index((iwn, iK, l1, l2), dims=(self.NwG4, self.Nc, self.nOrb, self.nOrb)) + self.cG1[iwn, iK, :, :] = cG1 + self.cG2[iwn, iK, :, :] = cG2 + for ind2 in range(self.nOrb*self.nOrb): + (l3, l4) = np.unravel_index(ind2, (self.nOrb, self.nOrb)) + iK2 = np.ravel_multi_index((iwn, iK, l3, l4), dims=(self.NwG4, self.Nc, self.nOrb, self.nOrb)) + ccind1 = np.ravel_multi_index((l3, l1), dims=(self.nOrb, self.nOrb)) + ccind2 = np.ravel_multi_index((l2, l4), dims=(self.nOrb, self.nOrb)) + chi0[iK1, iK2] = cc[ccind1, ccind2] + + elif "PARTICLE_PARTICLE" in self.channel: + iK2 = self.iKSum[self.iKDiff[self.iK0, iK], self.Qindices[iq]] # -k+Q + for ik, k in enumerate(kPatch): + e1 = ek1[iK, ik, ...] + e2 = ek2[iK, ik, ...] + iwG = int(iwn - self.iwG40 + self.iwG0) + miwPlusiwm = self.iwnPlusiwm(NwG-iwG-1, iwm, maxIndex=NwG) # iwn+iwm + sigmaK1 = self.sigma[iwG, iK, :, :] + sigmaK2 = self.sigma[miwPlusiwm, iK2, :, :] + G1 = np.linalg.inv((1j*wn + self.mu)*np.identity(nOrb) - e1 - sigmaK1) + G2 = np.linalg.inv((1j*self.wnG[miwPlusiwm] + self.mu)*np.identity(nOrb) - e2 - sigmaK2) + cc += np.outer(G1.reshape(self.nOrb*self.nOrb), G2.reshape(self.nOrb*self.nOrb)) + + for ind1 in range(self.nOrb*self.nOrb): + (l1, l2) = np.unravel_index(ind1, (self.nOrb, self.nOrb)) + iK1 = np.ravel_multi_index((iwn, iK, l1, l2), dims=(self.NwG4, self.Nc, self.nOrb, self.nOrb)) + self.cG1[iwn, iK, :, :] = cG1 + self.cG2[iwn, iK, :, :] = cG2 + for ind2 in range(self.nOrb*self.nOrb): + (l3, l4) = np.unravel_index(ind2, (self.nOrb, self.nOrb)) + iK2 = np.ravel_multi_index((iwn, iK, l3, l4), dims=(self.NwG4, self.Nc, self.nOrb, self.nOrb)) + ccind1 = np.ravel_multi_index((l1, l3), dims=(self.nOrb, self.nOrb)) + ccind2 = np.ravel_multi_index((l2, l4), dims=(self.nOrb, self.nOrb)) + chi0[iK1, iK2] = cc[ccind1, ccind2] + + chi0 /= kPatch.shape[0] + self.cG1 /= kPatch.shape[0] + self.cG2 /= kPatch.shape[0] + + return chi0 + + def calcG4Lattice(self): + nt = self.nt + self.G4L = np.zeros_like(self.G4) + self.GammaIrr = np.zeros_like(self.G4) + for iQ in range(self.G4.shape[0]): + G4 = self.G4[iQ, ...] + chic0 = self.chic0[iQ, ...] + chi0 = self.chi0[iQ, ...] + G4Inv = np.linalg.inv(G4) + chic0Inv = np.linalg.inv(chic0) + chi0Inv = np.linalg.inv(chi0) + self.G4L[iQ, ...] = np.linalg.inv(chi0Inv - chic0Inv + G4Inv) + self.GammaIrr[iQ, ...] = (chic0Inv - G4Inv) * self.invT * self.Nc + + def diagonalizePPKernel(self, iQ=0): + self.pm = np.dot(self.GammaIrr[iQ, ...], self.chi0[iQ, ...]) / (self.invT * self.Nc) + w, v = np.linalg.eig(self.pm) + wt = abs(w - 1) + ilead = np.argsort(wt) + self.lambdas = w[ilead] + self.evecs = v[:, ilead] + + def diagonalizeSymmPPKernel(self, iQ=0): + chiD = linalg.eig(self.chi0[iQ, ...]) + +# Helper functions + + def unfoldG4(self, G4): + return G4.reshape(self.Nwm, self.Nq, self.NwG4, self.Nc, self.nOrb, self.nOrb, self.NwG4, self.Nc, self.nOrb, self.nOrb) + + def unfold_K_Index(self, index): + return np.unravel_index(index, (self.NwG4, self.Nc, self.nOrb, self.nOrb)) + + def setupMomentumTables(self): + # build tables for K+K' and K-K' + self.iK0 = self.K_2_iK(0.0, 0.0) + self.iKDiff = zeros((self.Nc, self.Nc), dtype='int') + self.iKSum = zeros((self.Nc, self.Nc), dtype='int') + Nc = self.Nc + for iK1 in range(Nc): + Kx1 = self.Kvecs[iK1, 0] + Ky1 = self.Kvecs[iK1, 1] + for iK2 in range(0, Nc): + Kx2 = self.Kvecs[iK2, 0] + Ky2 = self.Kvecs[iK2, 1] + iKS = self.K_2_iK(Kx1+Kx2, Ky1+Ky2) + iKD = self.K_2_iK(Kx1-Kx2, Ky1-Ky2) + self.iKDiff[iK1, iK2] = iKD + self.iKSum[iK1, iK2] = iKS + + def K_2_iK(self, Kx, Ky): + delta = 1.0e-5 + K = np.array([Kx, Ky]) + M = np.column_stack((self.b0, self.b1)) + MTM = np.dot(M.T, M) + MTMinv = np.linalg.inv(MTM) + + for iK in range(0, self.Nc): + K_minus_Ktarget = K - self.Kvecs[iK] + Y = np.dot(np.dot(MTMinv, M.T), K_minus_Ktarget) + # print("iK, Y: ", iK, Y) + if (abs(Y[0] - round(Y[0])) < delta) & (abs(Y[1] - round(Y[1])) < delta): + return iK + print("No Kvec found!!!", Kx, Ky) + + def iwnPlusiwm(self, iwn, iwm, maxIndex): # iwn + iwm (wn fermionic, wm bosonic) + return int(min([max([iwn + iwm, 0]), maxIndex-1])) + + def minusiwn(self, iwn, maxIndex): # -iwn (wn fermionic) + return maxIndex-iwn-1 + + def iw1Minusiw2(self, iw1, iw2, maxIndex): # both iw1 and iw2 are fermionic + iww = min(max(iw1 - iw2, 0), maxIndex) + return iww + + def build_bStar(self, cluster=False): + # Build 8 rec. lattice vectors that surround K0=(0,0) + if cluster is False: + b0 = self.b0 # use rec. lattice vectors + b1 = self.b1 + else: + b0 = self.bc0 # use basis vectors for cluster K-points + b1 = self.bc1 + b = np.zeros((8, 2)) + b[0] = b0 + b[1] = b1 + b[2] = -b0 + b[3] = -b1 + b[4] = b0 + b1 + b[5] = b0 - b1 + b[6] = -b0 + b1 + b[7] = -b0 - b1 + if cluster is False: + self.bStar = b + else: + self.bcStar = b + + def shift21BZ(self, Kin): + b = self.bStar + # Now check if K-point is closest to K0=(0,0) and, if not, subtract from it the b-vector that it is closest to + K = Kin.copy() + for iK, Kvec in enumerate(K): + dist0 = np.linalg.norm(Kvec) + distClosest = dist0 + ibClosest = -1 # corresponds to K=(0,0) + for ib in range(8): + distb = np.linalg.norm(Kvec - b[ib]) + if distb < distClosest: + distClosest = distb + ibClosest = ib + if ibClosest > -1: + Kvec -= b[ibClosest] + return K + + def build_kGrid(self, nkfine=128, cluster=False): + # build fine k mesh for 1. BZ if cluster = False, for DCA patch if cluster = True + if cluster is False: + b = self.bStar + b0 = self.b0 + b1 = self.b1 + else: + b = self.bcStar + b0 = self.bc0 + b1 = self.bc1 + + k = [] + for ik1 in range(nkfine): + for ik2 in range(nkfine): + k.append(ik1/float(nkfine)*b0 + ik2/float(nkfine)*b1) + + # Now check if k-point is closest to K0=(0,0) and, if not, subtract the b-vector + for iK, Kvec in enumerate(k): + dist0 = np.linalg.norm(Kvec) + distClosest = dist0 + ibClosest = -1 # corresponds to K=(0,0) + for ib in range(8): + distb = np.linalg.norm(Kvec - b[ib]) + if distb < distClosest: + distClosest = distb + ibClosest = ib + if ibClosest > -1: + Kvec -= b[ibClosest] + k = np.stack(list(k)) + return k + +# Plotting functions + + def plotKvecs(self, shiftTo1BZ=True, plotBZ=True, plotPatch=True): + if shiftTo1BZ is True: + Kvecs = self.shift21BZ(self.Kvecs) + else: + Kvecs = self.Kvecs + if plotBZ is True: + k = self.build_kGrid() + if plotPatch is True: + kP = self.build_kGrid(cluster=True) + fig, ax = plt.subplots(figsize=(6, 6)) + if plotBZ is True: + ax.scatter(k[:, 0], k[:, 1], alpha=0.1, c="lightgrey") + if plotPatch is True: + ax.scatter(kP[:, 0], kP[:, 1], alpha=0.1, c="grey") + + ax.scatter(Kvecs[:, 0], Kvecs[:, 1]) + ax.annotate("", xy=(self.b0[0], self.b0[1]), xytext=(0, 0), arrowprops=dict(arrowstyle="->")) + ax.annotate("", xy=(self.b1[0], self.b1[1]), xytext=(0, 0), arrowprops=dict(arrowstyle="->")) + + for iK in range(self.Kvecs.shape[0]): + ax.text(x=Kvecs[iK, 0], y=Kvecs[iK, 1], s=str(iK)) + ax.set_aspect('equal') + + def plotEvecs(self, evecs, id=0): + + fig, ax = plt.subplots(nrows=self.nOrb, ncols=self.nOrb, figsize=(4*self.nOrb, 4*self.nOrb), sharey=True, sharex=True) + for i1 in range(self.nOrb): + for i2 in range(self.nOrb): + for iK in range(self.Nc): + ax[i1, i2].plot(self.wnSet, evecs[:, id].reshape(self.NwG4, self.Nc, self.nOrb, self.nOrb)[:, iK, i1, i2].real, label="real; K="+str(iK)) + ax[i1, i2].plot(self.wnSet, evecs[:, id].reshape(self.NwG4, self.Nc, self.nOrb, self.nOrb)[:, iK, i1, i2].imag, '--', label="imag; K="+str(iK)) + ax[i1, i2].set_title(str(i1)+str(i2)) + [ax[i1, -1].set(xlabel=r"$\omega_n$") for i1 in range(self.nOrb)] + + # ax[1, 0].legend(bbox_to_anchor=(1.2, 1.0), loc='right') + # plt.subplots_adjust(left=0.1, bottom=0.1, right=0.88, top=0.9, wspace=0.5, hspace=0.2) + fig.suptitle(r"$\lambda=$"+str(self.lambdas[id])) + + def plotEvecs3(self, id, imag=False): + if imag: + ev = self.evecs[:, id].imag.reshape(self.NwG4, self.Nc, self.nOrb, self.nOrb) + else: + ev = self.evecs[:, id].real.reshape(self.NwG4, self.Nc, self.nOrb, self.nOrb) + + dd = pd.DataFrame([[self.wnSet[j], iK, str(s1)+str(s2), ev[j, iK, s1, s2]] for s1 in range( + self.nOrb) for s2 in range(self.nOrb) for iK in range(self.Nc) for j in range(self.NwG4)]) + dd.columns = [r"$\omega_n$", "K", "orbital", r"$\phi'$"] + # theme_set(theme_minimal) + self.evplot = (ggplot(dd, aes(x=r"$\omega_n$", y=r"$\phi'$", color="factor(orbital)")) + + geom_line() + + facet_wrap('K') + + labs(color="orbital") + + theme_minimal(base_size=16, base_family="Arial") + + theme(panel_grid_major=element_line(color="Darkgrey")) + + theme(figure_size=(10, 10)) + + theme(plot_background=element_rect(fill='white', color='white')) + # + theme(strip_background=element_rect(fill="Darkgrey", size=1.4, alpha=.95),) + ) + + print(self.evplot) + +if __name__ == '__main__': + leading_eigenvalues = [] + leading_eigenvecs = [] + for filename in ['T=0.1/dca_tp.hdf5','T=0.09/dca_tp.hdf5','T=0.08/dca_tp.hdf5','T=0.07/dca_tp.hdf5']: + print('\n=== Processing', filename, '===') + T_analysis = Analyze(filename,channel='PARTICLE_PARTICLE_UP_DOWN') + T_analysis.calcChi0Cluster() + T_analysis.buildFullChi0Lattice() + T_analysis.calcG4Lattice() + T_analysis.diagonalizePPKernel() + print('Leading eigenvalues (sorted by distance from 1):') + print(T_analysis.lambdas[:5]) + print('Leading eigenvalue (real part):', T_analysis.lambdas[0].real) + leading_eigenvalues.append(T_analysis.lambdas.real) + leading_eigenvecs.append(T_analysis.evecs[:,0].reshape(T_analysis.NwG4, T_analysis.Nc).real) From 928a2d4c1e3d387de69500518b49494f269e53db Mon Sep 17 00:00:00 2001 From: Peter Doak <11856692+PDoakORNL@users.noreply.github.com> Date: Mon, 1 Jun 2026 18:11:32 -0400 Subject: [PATCH 4/7] return symmetrization for ppud anaylysis --- .../dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp | 3 +-- .../dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/include/dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp b/include/dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp index ab84e8214..6a98e8298 100644 --- a/include/dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp +++ b/include/dca/phys/dca_analysis/bse_solver/bse_cluster_solver.hpp @@ -145,8 +145,7 @@ void BseClusterSolver::apply_symmetries if (concurrency.id() == concurrency.first()) std::cout << "\t" << __FUNCTION__ << "\n\n"; - if (parameters.symmetrize_Gamma() && - parameters.get_four_point_channels()[0] != FourPointType::PARTICLE_PARTICLE_UP_DOWN) { + if (parameters.symmetrize_Gamma()) { if (true) { if (concurrency.id() == concurrency.first()) std::cout << "symmetrize Gamma_lattice according to the symmetry-group \n" << std::endl; diff --git a/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp b/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp index faf9605a6..ba97d57f8 100644 --- a/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp +++ b/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp @@ -389,8 +389,7 @@ void BseLatticeSolver::computeGammaLattice( Gamma_lattice(i) = Gamma_cluster(i); } - if (parameters.symmetrize_Gamma() && - parameters.get_four_point_channels()[0] != FourPointType::PARTICLE_PARTICLE_UP_DOWN) { + if (parameters.symmetrize_Gamma()) { if (concurrency.id() == concurrency.first()) std::cout << "Symmetrize Gamma_lattice according to the symmetry group." << std::endl; Symmetrize::execute(Gamma_lattice, parameters.get_four_point_momentum_transfer()); From ea8122a16f83fe4ca5ca3fe083a3200d4a89ce15 Mon Sep 17 00:00:00 2001 From: Peter Doak <11856692+PDoakORNL@users.noreply.github.com> Date: Wed, 17 Jun 2026 18:56:31 +0000 Subject: [PATCH 5/7] Add direct-grid chi0 lattice calculation with Python agreement Implements a direct uniform-grid chi0 lattice calculation in BseLatticeSolver that matches Python's buildFullChi0Lattice to machine precision. Changes: - Add analysis parameters direct-grid-chi0 and direct-grid-nkfine - Add evaluateH0AtK() pointwise H0 interface to all lattice models - Implement computeChi0DirectGrid() in BseLatticeSolver: * per-patch grid generation using cluster basis vectors * Voronoi-cell shift matching Python's build_kGrid(cluster=True) * correct use of fixed four-point momentum transfer Q in the bubble * correct PP channel momentum combination Q-k - Update unit tests for new parameters and square_lattice evaluateH0AtK Verified on T=0.07, 0.08, 0.09, 0.10: C++ and Python chi0 lattice agree to relative error ~5e-15 (machine precision). --- .../bse_solver/bse_lattice_solver.hpp | 307 ++++++++++++++++-- .../analytic_hamiltonians/cubic_lattice.hpp | 7 + .../fourband_lattice.hpp | 7 + .../honeycomb_lattice.hpp | 7 + .../analytic_hamiltonians/square_lattice.hpp | 8 + .../analytic_hamiltonians/twoband_lattice.hpp | 7 + .../CuO2/material_lattice_CuO2.inc | 6 + .../material_lattice_CuO2_1band.inc | 6 + .../FeSn/material_lattice_FeSn.inc | 8 + .../NiO/material_lattice_NiO.inc | 8 + .../SrVO3/material_lattice_SrVO3.inc | 6 + .../phys/parameters/analysis_parameters.hpp | 24 ++ .../square_lattice_test.cpp | 21 ++ .../analysis_parameters_test.cpp | 5 + .../analysis_parameters/input_read_all.json | 4 +- 15 files changed, 403 insertions(+), 28 deletions(-) diff --git a/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp b/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp index ba97d57f8..626030e4f 100644 --- a/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp +++ b/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp @@ -37,6 +37,7 @@ #include "dca/phys/domains/cluster/centered_cluster_domain.hpp" #include "dca/phys/domains/cluster/cluster_domain.hpp" #include "dca/phys/domains/quantum/electron_band_domain.hpp" +#include "dca/phys/domains/time_and_frequency/frequency_domain.hpp" #include "dca/phys/domains/time_and_frequency/vertex_frequency_domain.hpp" #include "dca/util/print_time.hpp" @@ -98,6 +99,7 @@ class BseLatticeSolver { void write(Writer& writer); void computeChi0Lattice(); + void computeChi0DirectGrid(); template void computeGammaLattice( /*const*/ func::function, ClusterMatrixDmn>& Gamma_cluster); @@ -312,48 +314,299 @@ void BseLatticeSolver::computeChi0Lattice() if (concurrency.id() == concurrency.first()) std::cout << "\n" << __FUNCTION__ << std::endl; - clustermapping::coarsegraining_tp coarsegraining_tp(parameters); + if (parameters.direct_grid_chi0()) { + computeChi0DirectGrid(); + } + else { + clustermapping::coarsegraining_tp coarsegraining_tp(parameters); - // DCA+/DCA with post-interpolation: Compute \chi_0 from continuous lattice self-energy. - if (parameters.do_dca_plus() || parameters.doPostInterpolation()) { - latticemapping::lattice_mapping_sp lattice_map_sp(parameters); + // DCA+/DCA with post-interpolation: Compute \chi_0 from continuous lattice self-energy. + if (parameters.do_dca_plus() || parameters.doPostInterpolation()) { + latticemapping::lattice_mapping_sp lattice_map_sp(parameters); - MOMS.Sigma_lattice_interpolated = 0.; - MOMS.Sigma_lattice_coarsegrained = 0.; - MOMS.Sigma_lattice = 0.; + MOMS.Sigma_lattice_interpolated = 0.; + MOMS.Sigma_lattice_coarsegrained = 0.; + MOMS.Sigma_lattice = 0.; - if (parameters.hts_approximation()) { - clustermapping::CoarsegrainingSp CoarsegrainingSp(parameters); + if (parameters.hts_approximation()) { + clustermapping::CoarsegrainingSp CoarsegrainingSp(parameters); - DcaDataType dca_data_hts(parameters); - dca_data_hts.H_HOST = MOMS.H_HOST; - dca_data_hts.H_interactions = MOMS.H_interactions; + DcaDataType dca_data_hts(parameters); + dca_data_hts.H_HOST = MOMS.H_HOST; + dca_data_hts.H_interactions = MOMS.H_interactions; - solver::HighTemperatureSeriesExpansionSolver hts_solver( - parameters, dca_data_hts); + solver::HighTemperatureSeriesExpansionSolver hts_solver( + parameters, dca_data_hts); - lattice_map_sp.execute_with_HTS_approximation( - dca_data_hts, hts_solver, CoarsegrainingSp, MOMS.Sigma, MOMS.Sigma_lattice_interpolated, - MOMS.Sigma_lattice_coarsegrained, MOMS.Sigma_lattice); + lattice_map_sp.execute_with_HTS_approximation( + dca_data_hts, hts_solver, CoarsegrainingSp, MOMS.Sigma, MOMS.Sigma_lattice_interpolated, + MOMS.Sigma_lattice_coarsegrained, MOMS.Sigma_lattice); + } + else { + lattice_map_sp.execute(MOMS.Sigma, MOMS.Sigma_lattice_interpolated, + MOMS.Sigma_lattice_coarsegrained, MOMS.Sigma_lattice); + } + + if (parameters.do_dca_plus()) + coarsegraining_tp.execute(MOMS.H_HOST, MOMS.Sigma_lattice, chi_0_lattice); + + else // do_post_interpolation + coarsegraining_tp.execute(MOMS.H_HOST, MOMS.Sigma_lattice_interpolated, chi_0_lattice); } + + // (Standard) DCA: Compute \chi_0 from cluster self-energy. else { - lattice_map_sp.execute(MOMS.Sigma, MOMS.Sigma_lattice_interpolated, - MOMS.Sigma_lattice_coarsegrained, MOMS.Sigma_lattice); + coarsegraining_tp.execute(MOMS.H_HOST, MOMS.Sigma, chi_0_lattice); } - if (parameters.do_dca_plus()) - coarsegraining_tp.execute(MOMS.H_HOST, MOMS.Sigma_lattice, chi_0_lattice); + // Set diagonal \chi_0 matrix. + const ScalarType renorm = 1.; - else // do_post_interpolation - coarsegraining_tp.execute(MOMS.H_HOST, MOMS.Sigma_lattice_interpolated, chi_0_lattice); + for (int w_ind = 0; w_ind < WVertexDmn::dmn_size(); w_ind++) + for (int K_ind = 0; K_ind < k_HOST_VERTEX::dmn_size(); K_ind++) + for (int m2 = 0; m2 < b::dmn_size(); m2++) + for (int n2 = 0; n2 < b::dmn_size(); n2++) + for (int m1 = 0; m1 < b::dmn_size(); m1++) + for (int n1 = 0; n1 < b::dmn_size(); n1++) { + chi_0_lattice(n1, m1, n2, m2, K_ind, w_ind) *= renorm; + chi_0_lattice_matrix(n1, m1, K_ind, w_ind, n2, m2, K_ind, w_ind) = + chi_0_lattice(n1, m1, n2, m2, K_ind, w_ind); + } } +} - // (Standard) DCA: Compute \chi_0 from cluster self-energy. - else { - coarsegraining_tp.execute(MOMS.H_HOST, MOMS.Sigma, chi_0_lattice); +template +void BseLatticeSolver::computeChi0DirectGrid() { + profiler_type prof(__FUNCTION__, "BseLatticeSolver", __LINE__); + + if (concurrency.id() == concurrency.first()) + std::cout << "\n" << __FUNCTION__ << std::endl; + + chi_0_lattice = 0.; + + const int nkfine = parameters.direct_grid_nkfine(); + constexpr int dim = Parameters::lattice_type::DIMENSION; + + using LatticeType = typename Parameters::model_type::lattice_type; + using w = func::dmn_0; + + const auto& super_basis = k_DCA::parameter_type::get_super_basis_vectors(); + const auto& basis = k_DCA::parameter_type::get_basis_vectors(); + + const int n_w = w::dmn_size(); + const int n_q = k_HOST_VERTEX::dmn_size(); + const int n_w_vertex = WVertexDmn::dmn_size(); + + const int Q_ind = domains::cluster_operations::index( + parameters.get_four_point_momentum_transfer(), k_HOST_VERTEX::get_elements(), + k_HOST_VERTEX::parameter_type::SHAPE); + (void)Q_ind; + + const auto& Q = parameters.get_four_point_momentum_transfer(); + + const int W_ind = parameters.get_four_point_frequency_transfer(); + const auto channel = parameters.get_four_point_channels()[0]; + + // Number of patch points per cluster K-point. + int points_per_patch = 1; + for (int d = 0; d < dim; ++d) + points_per_patch *= nkfine; + + // Total number of (cluster K, patch point) pairs. + int nk_total = n_q * points_per_patch; + + // Generate cluster star vectors (all nearest-neighbor reciprocal lattice vectors). + // These are used to shift each patch point to the Voronoi cell around the origin, + // matching the Python build_kGrid(cluster=True) behavior. + std::vector> bcStar; + { + std::vector coeffs = {-1, 0, 1}; + if constexpr (dim == 2) { + for (int c0 : coeffs) + for (int c1 : coeffs) { + if (c0 == 0 && c1 == 0) + continue; + std::vector v(dim, 0.); + for (int d = 0; d < dim; ++d) + v[d] = c0 * basis[0][d] + c1 * basis[1][d]; + bcStar.push_back(v); + } + } + else if constexpr (dim == 3) { + for (int c0 : coeffs) + for (int c1 : coeffs) + for (int c2 : coeffs) { + if (c0 == 0 && c1 == 0 && c2 == 0) + continue; + std::vector v(dim, 0.); + for (int d = 0; d < dim; ++d) + v[d] = c0 * basis[0][d] + c1 * basis[1][d] + c2 * basis[2][d]; + bcStar.push_back(v); + } + } + } + + // Integration factor matching coarsegraining_tp::get_integration_factor(). + double factor = -1.; + switch (channel) { + case FourPointType::PARTICLE_HOLE_TRANSVERSE: + case FourPointType::PARTICLE_HOLE_MAGNETIC: + factor = -1.; + break; + case FourPointType::PARTICLE_HOLE_CHARGE: + factor = -2.; + break; + case FourPointType::PARTICLE_PARTICLE_UP_DOWN: + factor = 1.; + break; + default: + throw std::logic_error(__FUNCTION__); + } + + // MPI round-robin over (cluster K, patch point) pairs. + for (int idx = concurrency.id(); idx < nk_total; + idx += concurrency.number_of_processors()) { + int qi = idx / points_per_patch; + int patch_idx = idx % points_per_patch; + + int ix, iy, iz = 0; + if constexpr (dim == 2) { + ix = patch_idx % nkfine; + iy = patch_idx / nkfine; + } + else if constexpr (dim == 3) { + ix = patch_idx % nkfine; + int rem = patch_idx / nkfine; + iy = rem % nkfine; + iz = rem / nkfine; + } + + const auto& K = k_HOST_VERTEX::get_elements()[qi]; + + // Build patch point using cluster basis vectors (supercell reciprocal lattice). + // This ensures each patch covers exactly 1/Nc of the BZ. + std::vector k_patch(dim, 0.); + for (int d = 0; d < dim; ++d) { + k_patch[d] += basis[0][d] * ix / nkfine; + k_patch[d] += basis[1][d] * iy / nkfine; + if constexpr (dim == 3) + k_patch[d] += basis[2][d] * iz / nkfine; + } + + // Shift k_patch to the Voronoi cell around the origin, matching Python's + // build_kGrid(cluster=True) behavior. + { + double dist0 = math::util::l2Norm2(k_patch); + int ibClosest = -1; + for (int ib = 0; ib < bcStar.size(); ++ib) { + double distb = math::util::distance2(k_patch, bcStar[ib]); + if (distb < dist0) { + dist0 = distb; + ibClosest = ib; + } + } + if (ibClosest > -1) { + for (int d = 0; d < dim; ++d) + k_patch[d] -= bcStar[ibClosest][d]; + } + } + + // Absolute k-point = cluster K + patch point. + std::vector k_abs(dim); + for (int d = 0; d < dim; ++d) + k_abs[d] = K[d] + k_patch[d]; + + // Wrap to first primitive BZ using super_basis (primitive reciprocal lattice vectors). + const auto k = domains::cluster_operations::translate_inside_cluster(k_abs, super_basis); + + // Use the patch K-point directly for the self-energy lookup. + // Patch points are generated within the DCA patch of K_qi using cluster basis vectors, + // so k is guaranteed to be in the correct patch. + const int k_cluster_ind = qi; + + // Evaluate bare dispersion at k. + const double H0_k = LatticeType::evaluateH0AtK(parameters, k); + + // Compute second momentum using the fixed momentum transfer Q. + std::vector k_second_raw(dim); + if (channel == FourPointType::PARTICLE_PARTICLE_UP_DOWN) { + for (int d = 0; d < dim; ++d) + k_second_raw[d] = Q[d] - k_abs[d]; + } + else { + for (int d = 0; d < dim; ++d) + k_second_raw[d] = k_abs[d] + Q[d]; + } + const auto kq_wrapped = + domains::cluster_operations::translate_inside_cluster(k_second_raw, super_basis); + + // Use the same patch K-point for the second G's self-energy lookup. + // For PP with Q=0, -K is equivalent to K modulo the cluster superlattice. + const int kq_cluster_ind = qi; + + // Evaluate bare dispersion at k+q. + const double H0_kq = LatticeType::evaluateH0AtK(parameters, kq_wrapped); + + // Precompute G(k, w) for all Matsubara frequencies. + std::vector> G_k(n_w); + for (int wi = 0; wi < n_w; ++wi) { + const std::complex iw(0., w::get_elements()[wi]); + // Single-band, up-spin (nu index 0 = b=0, s=0). + const std::complex Sigma_val = MOMS.Sigma(0, 0, k_cluster_ind, wi); + G_k[wi] = std::complex(1., 0.) / + (iw + parameters.get_chemical_potential() - H0_k - Sigma_val); + } + + // Precompute G(k+q, w) for all Matsubara frequencies. + std::vector> G_kq(n_w); + for (int wi = 0; wi < n_w; ++wi) { + const std::complex iw(0., w::get_elements()[wi]); + const std::complex Sigma_val = MOMS.Sigma(0, 0, kq_cluster_ind, wi); + G_kq[wi] = std::complex(1., 0.) / + (iw + parameters.get_chemical_potential() - H0_kq - Sigma_val); + } + + // Accumulate the bubble for each vertex frequency. + for (int w_ind = 0; w_ind < n_w_vertex; ++w_ind) { + // Find the w index matching the vertex frequency. + int w1 = -1; + for (int l = 0; l < n_w; ++l) { + if (std::abs(WVertexDmn::get_elements()[w_ind] - w::get_elements()[l]) < 1.e-6) { + w1 = l; + break; + } + } + if (w1 < 0) + throw std::logic_error( + "Could not find matching Matsubara frequency index in computeChi0DirectGrid."); + + int w2 = -1; + switch (channel) { + case FourPointType::PARTICLE_HOLE_CHARGE: + case FourPointType::PARTICLE_HOLE_MAGNETIC: + case FourPointType::PARTICLE_HOLE_TRANSVERSE: { + w2 = w1 + W_ind; + } break; + case FourPointType::PARTICLE_PARTICLE_UP_DOWN: { + w2 = W_ind + (n_w - 1 - w1); + } break; + default: + throw std::logic_error(__FUNCTION__); + } + + // Single-band: all band indices are 0. + chi_0_lattice(0, 0, 0, 0, qi, w_ind) += + factor * G_k[w1] * G_kq[w2]; + } } - // Set diagonal \chi_0 matrix. + // Collect contributions from all MPI ranks. + concurrency.sum(chi_0_lattice); + + // Normalize by the number of patch points per cluster K-point. + chi_0_lattice /= static_cast(points_per_patch); + + // Set diagonal chi_0 matrix (same as in computeChi0Lattice). const ScalarType renorm = 1.; for (int w_ind = 0; w_ind < WVertexDmn::dmn_size(); w_ind++) diff --git a/include/dca/phys/models/analytic_hamiltonians/cubic_lattice.hpp b/include/dca/phys/models/analytic_hamiltonians/cubic_lattice.hpp index 538d0991e..4a63472e7 100644 --- a/include/dca/phys/models/analytic_hamiltonians/cubic_lattice.hpp +++ b/include/dca/phys/models/analytic_hamiltonians/cubic_lattice.hpp @@ -57,6 +57,13 @@ class cubic_lattice { template static std::complex get_LDA_Hamiltonians(parameters_type& parameters, std::vector k, int b1, int s1, int b2, int s2); + + template + static double evaluateH0AtK(const ParametersType& parameters, const std::vector& k) { + return std::real( + get_LDA_Hamiltonians(const_cast(parameters), k, 0, 0, 0, 0)); + } + }; template diff --git a/include/dca/phys/models/analytic_hamiltonians/fourband_lattice.hpp b/include/dca/phys/models/analytic_hamiltonians/fourband_lattice.hpp index a9718a7fd..19ecf1fb7 100644 --- a/include/dca/phys/models/analytic_hamiltonians/fourband_lattice.hpp +++ b/include/dca/phys/models/analytic_hamiltonians/fourband_lattice.hpp @@ -72,6 +72,13 @@ class fourband_lattice { static std::complex get_LDA_Hamiltonians(parameters_type& parameters, std::vector /*k*/, int b1, int s1, int b2, int s2); + + template + static double evaluateH0AtK(const ParametersType& parameters, const std::vector& k) { + return std::real( + get_LDA_Hamiltonians(const_cast(parameters), k, 0, 0, 0, 0)); + } + }; template diff --git a/include/dca/phys/models/analytic_hamiltonians/honeycomb_lattice.hpp b/include/dca/phys/models/analytic_hamiltonians/honeycomb_lattice.hpp index 3446e4912..b24572bcd 100644 --- a/include/dca/phys/models/analytic_hamiltonians/honeycomb_lattice.hpp +++ b/include/dca/phys/models/analytic_hamiltonians/honeycomb_lattice.hpp @@ -64,6 +64,13 @@ class HoneycombLattice { template static std::complex get_LDA_Hamiltonians(parameters_type& parameters, std::vector k, int b1, int s1, int b2, int s2); + + template + static double evaluateH0AtK(const ParametersType& parameters, const std::vector& k) { + return std::real( + get_LDA_Hamiltonians(const_cast(parameters), k, 0, 0, 0, 0)); + } + }; template diff --git a/include/dca/phys/models/analytic_hamiltonians/square_lattice.hpp b/include/dca/phys/models/analytic_hamiltonians/square_lattice.hpp index e5a2e1c16..ccc76076d 100644 --- a/include/dca/phys/models/analytic_hamiltonians/square_lattice.hpp +++ b/include/dca/phys/models/analytic_hamiltonians/square_lattice.hpp @@ -78,6 +78,14 @@ class square_lattice { func::function, func::dmn_variadic, KDmn>>& H_0, typename KDmn::element_type& q); + template + static double evaluateH0AtK(const ParametersType& parameters, const std::vector& k) { + const auto t = parameters.get_t(); + const auto t_prime = parameters.get_t_prime(); + return -2. * t * (std::cos(k[0]) + std::cos(k[1])) - + 4. * t_prime * std::cos(k[0]) * std::cos(k[1]); + } + }; template diff --git a/include/dca/phys/models/analytic_hamiltonians/twoband_lattice.hpp b/include/dca/phys/models/analytic_hamiltonians/twoband_lattice.hpp index db48e788c..e534e3017 100644 --- a/include/dca/phys/models/analytic_hamiltonians/twoband_lattice.hpp +++ b/include/dca/phys/models/analytic_hamiltonians/twoband_lattice.hpp @@ -69,6 +69,13 @@ class twoband_lattice { static std::complex get_LDA_Hamiltonians(parameters_type& parameters, std::vector /*k*/, int b1, int s1, int b2, int s2); + + template + static double evaluateH0AtK(const ParametersType& parameters, const std::vector& k) { + return std::real( + get_LDA_Hamiltonians(const_cast(parameters), k, 0, 0, 0, 0)); + } + }; template diff --git a/include/dca/phys/models/material_hamiltonians/CuO2/material_lattice_CuO2.inc b/include/dca/phys/models/material_hamiltonians/CuO2/material_lattice_CuO2.inc index 626627b1d..595989f64 100644 --- a/include/dca/phys/models/material_hamiltonians/CuO2/material_lattice_CuO2.inc +++ b/include/dca/phys/models/material_hamiltonians/CuO2/material_lattice_CuO2.inc @@ -42,6 +42,12 @@ public: static std::complex get_LDA_Hamiltonians(parameters_type& parameters, std::vector k, int b1, int s1, int b2, int s2); + template + static double evaluateH0AtK(const ParametersType& parameters, const std::vector& k) { + return std::real( + get_LDA_Hamiltonians(const_cast(parameters), k, 0, 0, 0, 0)); + } + private: template static std::vector>& get_t_ij(parameters_type& parameters); diff --git a/include/dca/phys/models/material_hamiltonians/CuO2_1band/material_lattice_CuO2_1band.inc b/include/dca/phys/models/material_hamiltonians/CuO2_1band/material_lattice_CuO2_1band.inc index c195d2c57..f85d97fe8 100644 --- a/include/dca/phys/models/material_hamiltonians/CuO2_1band/material_lattice_CuO2_1band.inc +++ b/include/dca/phys/models/material_hamiltonians/CuO2_1band/material_lattice_CuO2_1band.inc @@ -41,6 +41,12 @@ public: static std::complex get_LDA_Hamiltonians(parameters_type& parameters, std::vector k, int b1, int s1, int b2, int s2); + template + static double evaluateH0AtK(const ParametersType& parameters, const std::vector& k) { + return std::real( + get_LDA_Hamiltonians(const_cast(parameters), k, 0, 0, 0, 0)); + } + private: template static std::vector>& get_t_ij(parameters_type& parameters); diff --git a/include/dca/phys/models/material_hamiltonians/FeSn/material_lattice_FeSn.inc b/include/dca/phys/models/material_hamiltonians/FeSn/material_lattice_FeSn.inc index 5a022ce90..39e7488cb 100644 --- a/include/dca/phys/models/material_hamiltonians/FeSn/material_lattice_FeSn.inc +++ b/include/dca/phys/models/material_hamiltonians/FeSn/material_lattice_FeSn.inc @@ -61,6 +61,14 @@ public: func::function, func::dmn_variadic, KDmn>>& H_0); + template + static double evaluateH0AtK(const ParametersType& /*parameters*/, + const std::vector& /*k*/) { + throw std::logic_error( + "evaluateH0AtK not implemented for FeSn material lattice. " + "Use get_LDA_Hamiltonians or initializeH0 instead."); + } + static double latticeConstant() { return a; } diff --git a/include/dca/phys/models/material_hamiltonians/NiO/material_lattice_NiO.inc b/include/dca/phys/models/material_hamiltonians/NiO/material_lattice_NiO.inc index 7b41ebf59..2c316b568 100644 --- a/include/dca/phys/models/material_hamiltonians/NiO/material_lattice_NiO.inc +++ b/include/dca/phys/models/material_hamiltonians/NiO/material_lattice_NiO.inc @@ -72,6 +72,14 @@ public: func::function, func::dmn_variadic, KDmn>>& H_0); + template + static double evaluateH0AtK(const ParametersType& /*parameters*/, + const std::vector& /*k*/) { + throw std::logic_error( + "evaluateH0AtK not implemented for NiO material lattice. " + "Use get_LDA_Hamiltonians or initializeH0 instead."); + } + static double latticeConstant() { return a; } diff --git a/include/dca/phys/models/material_hamiltonians/SrVO3/material_lattice_SrVO3.inc b/include/dca/phys/models/material_hamiltonians/SrVO3/material_lattice_SrVO3.inc index 0be0657c9..a34b22126 100644 --- a/include/dca/phys/models/material_hamiltonians/SrVO3/material_lattice_SrVO3.inc +++ b/include/dca/phys/models/material_hamiltonians/SrVO3/material_lattice_SrVO3.inc @@ -41,6 +41,12 @@ public: static std::complex get_LDA_Hamiltonians(parameters_type& parameters, std::vector k, int b1, int s1, int b2, int s2); + template + static double evaluateH0AtK(const ParametersType& parameters, const std::vector& k) { + return std::real( + get_LDA_Hamiltonians(const_cast(parameters), k, 0, 0, 0, 0)); + } + private: template static std::vector>& get_t_ij(parameters_type& parameters); diff --git a/include/dca/phys/parameters/analysis_parameters.hpp b/include/dca/phys/parameters/analysis_parameters.hpp index 42c034e49..766e68c5a 100644 --- a/include/dca/phys/parameters/analysis_parameters.hpp +++ b/include/dca/phys/parameters/analysis_parameters.hpp @@ -70,6 +70,12 @@ class AnalysisParameters { bool get_dump_intermediates() const { return dump_intermediates_; } + bool direct_grid_chi0() const { + return direct_grid_chi0_; + } + int direct_grid_nkfine() const { + return direct_grid_nkfine_; + } private: bool symmetrize_Gamma_ = true; @@ -80,6 +86,8 @@ class AnalysisParameters { std::vector> q_host_; std::vector> q_host_fine_; bool dump_intermediates_ = false; + bool direct_grid_chi0_ = false; + int direct_grid_nkfine_ = 32; }; template @@ -95,6 +103,8 @@ int AnalysisParameters::getBufferSize(const Concurrency& concurrency) const { buffer_size += concurrency.get_buffer_size(q_host_fine_); buffer_size += concurrency.get_buffer_size(dump_intermediates_); + buffer_size += concurrency.get_buffer_size(direct_grid_chi0_); + buffer_size += concurrency.get_buffer_size(direct_grid_nkfine_); return buffer_size; } @@ -111,6 +121,8 @@ void AnalysisParameters::pack(const Concurrency& concurrency, char* buffer, int concurrency.pack(buffer, buffer_size, position, q_host_fine_); concurrency.pack(buffer, buffer_size, position, dump_intermediates_); + concurrency.pack(buffer, buffer_size, position, direct_grid_chi0_); + concurrency.pack(buffer, buffer_size, position, direct_grid_nkfine_); } template @@ -125,6 +137,8 @@ void AnalysisParameters::unpack(const Concurrency& concurrency, char* buffer, in concurrency.unpack(buffer, buffer_size, position, q_host_fine_); concurrency.unpack(buffer, buffer_size, position, dump_intermediates_); + concurrency.unpack(buffer, buffer_size, position, direct_grid_chi0_); + concurrency.unpack(buffer, buffer_size, position, direct_grid_nkfine_); } template @@ -155,6 +169,16 @@ void AnalysisParameters::readWrite(ReaderOrWriter& reader_or_writer) { } reader_or_writer.execute("dump-intermediates", dump_intermediates_); + try { + reader_or_writer.execute("direct-grid-chi0", direct_grid_chi0_); + } + catch (const std::exception& r_e) { + } + try { + reader_or_writer.execute("direct-grid-nkfine", direct_grid_nkfine_); + } + catch (const std::exception& r_e) { + } reader_or_writer.close_group(); } } diff --git a/test/unit/phys/models/analytic_hamiltonians/square_lattice_test.cpp b/test/unit/phys/models/analytic_hamiltonians/square_lattice_test.cpp index 2ae7c7665..3991ff412 100644 --- a/test/unit/phys/models/analytic_hamiltonians/square_lattice_test.cpp +++ b/test/unit/phys/models/analytic_hamiltonians/square_lattice_test.cpp @@ -145,3 +145,24 @@ TEST(SquareLatticeTest, Initialize_H_interaction) { else EXPECT_DOUBLE_EQ(0., H_interaction(0, s1, 0, s2, r)); } + +TEST(SquareLatticeTest, EvaluateH0AtK) { + using PointGroup = phys::domains::D4; + using Lattice = phys::models::square_lattice; + + phys::params::ModelParameters> params; + params.set_t(1.); + params.set_t_prime(0.5); + + std::vector k0 = {0., 0.}; + EXPECT_DOUBLE_EQ(-6., Lattice::evaluateH0AtK(params, k0)); + + std::vector k1 = {0., M_PI}; + EXPECT_DOUBLE_EQ(2., Lattice::evaluateH0AtK(params, k1)); + + std::vector k2 = {M_PI, 0.}; + EXPECT_DOUBLE_EQ(2., Lattice::evaluateH0AtK(params, k2)); + + std::vector k3 = {M_PI, M_PI}; + EXPECT_DOUBLE_EQ(2., Lattice::evaluateH0AtK(params, k3)); +} diff --git a/test/unit/phys/parameters/analysis_parameters/analysis_parameters_test.cpp b/test/unit/phys/parameters/analysis_parameters/analysis_parameters_test.cpp index 389e30a27..6b8e5c370 100644 --- a/test/unit/phys/parameters/analysis_parameters/analysis_parameters_test.cpp +++ b/test/unit/phys/parameters/analysis_parameters/analysis_parameters_test.cpp @@ -28,6 +28,9 @@ TEST(AnalysisParametersTest, DefaultValues) { EXPECT_EQ(q_host_check, pars.get_q_host()); std::vector> q_host_fine_check{{10, 0}, {0, 10}}; EXPECT_EQ(q_host_fine_check, pars.get_q_host_fine()); + + EXPECT_FALSE(pars.direct_grid_chi0()); + EXPECT_EQ(32, pars.direct_grid_nkfine()); } TEST(AnalysisParametersTest, ReadAll) { @@ -52,4 +55,6 @@ TEST(AnalysisParametersTest, ReadAll) { std::vector> q_host_fine_check{{100, 0}, {0, 100}}; EXPECT_EQ(q_host_fine_check, pars.get_q_host_fine()); + EXPECT_TRUE(pars.direct_grid_chi0()); + EXPECT_EQ(64, pars.direct_grid_nkfine()); } diff --git a/test/unit/phys/parameters/analysis_parameters/input_read_all.json b/test/unit/phys/parameters/analysis_parameters/input_read_all.json index b248fdc80..2fadf7de5 100644 --- a/test/unit/phys/parameters/analysis_parameters/input_read_all.json +++ b/test/unit/phys/parameters/analysis_parameters/input_read_all.json @@ -5,6 +5,8 @@ "project-onto-crystal-harmonics": true, "projection-cut-off-radius": 3.0, "q-host" : [[3, 3], [3, -3]], - "q-host-fine" : [[100, 0], [0, 100]] + "q-host-fine" : [[100, 0], [0, 100]], + "direct-grid-chi0": true, + "direct-grid-nkfine": 64 } } From c7c9d630e0d09faeed3e2c56446f5fc33a33b165 Mon Sep 17 00:00:00 2001 From: Peter Doak <11856692+PDoakORNL@users.noreply.github.com> Date: Wed, 17 Jun 2026 19:27:08 +0000 Subject: [PATCH 6/7] Use full complex diagonalization for PP Q=0 direct-grid chi0 For PARTICLE_PARTICLE_UP_DOWN with Q=0 and direct-grid-chi0=true, route diagonalization through diagonalizeGammaChi0Full() instead of diagonalizeGammaChi0Symmetric(). The symmetric method took the real part of sqrt(chi0)*Gamma*sqrt(chi0) before diagonalizing, which shifted eigenvalues by 1.5-7.6% relative to Python's full complex kernel Gamma*chi0/(beta*Nc). Also apply the 1/(beta*Nc) scaling in the full diagonalization path so the C++ kernel matches Python's exactly. Verified on T=0.07, 0.08, 0.09, 0.10: C++ leading eigenvalues now agree with Python to machine precision (max relative diff ~5e-14). --- .../bse_solver/bse_lattice_solver.hpp | 34 ++++++++++++++----- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp b/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp index 626030e4f..5b5948b80 100644 --- a/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp +++ b/include/dca/phys/dca_analysis/bse_solver/bse_lattice_solver.hpp @@ -119,7 +119,7 @@ class BseLatticeSolver { void initialize(); void diagonalizeGammaChi0Symmetric(); - void diagonalizeGammaChi0Full(); + void diagonalizeGammaChi0Full(ScalarType scale = 1.); void diagonalize_folded_Gamma_chi_0(); template // Element type of eigenvalues and eigenvectors. @@ -664,14 +664,21 @@ void BseLatticeSolver::diagonalizeGammaChi0 // Diagonalize the symmetric matrix \sqrt{\chi_0}\Gamma\sqrt{\chi_0}. // The origin in momentum space has always index = 0. // TODO: loop over multiple channels. - if (parameters.get_four_point_channels()[0] == FourPointType::PARTICLE_PARTICLE_UP_DOWN && + const bool is_pp_q0 = + parameters.get_four_point_channels()[0] == FourPointType::PARTICLE_PARTICLE_UP_DOWN && parameters.get_four_point_momentum_transfer_index() == 0 && - parameters.get_four_point_frequency_transfer() == 0) { + parameters.get_four_point_frequency_transfer() == 0; + if (is_pp_q0 && !parameters.direct_grid_chi0()) { + // Old symmetric method (real( sqrt(chi0) @ Gamma @ sqrt(chi0) )). diagonalizeGammaChi0Symmetric(); } - else + else { + // Scale the kernel by 1/(beta * Nc) to match Python's diagonalizePPKernel(), + // which diagonalizes GammaIrr @ chi0 / (beta * Nc). + const double scale = is_pp_q0 ? (1.0 / (parameters.get_beta() * k_HOST_VERTEX::dmn_size())) : 1.0; + diagonalizeGammaChi0Full(scale); + } #endif // DCA_ANALYSIS_TEST_WITH_FULL_DIAGONALIZATION - diagonalizeGammaChi0Full(); } characterizeLeadingEigenvectors(); @@ -736,7 +743,8 @@ void BseLatticeSolver::diagonalizeGammaChi0 } template -void BseLatticeSolver::diagonalizeGammaChi0Full() { +void BseLatticeSolver::diagonalizeGammaChi0Full( + ScalarType scale) { profiler_type prof(__FUNCTION__, "BseLatticeSolver", __LINE__); if (concurrency.id() == concurrency.first()) @@ -744,8 +752,12 @@ void BseLatticeSolver::diagonalizeGammaChi0 const int size = LatticeEigenvectorDmn::dmn_size(); - if (concurrency.id() == concurrency.first()) - std::cout << "Compute Gamma*chi_0: " << util::print_time() << std::endl; + if (concurrency.id() == concurrency.first()) { + if (std::abs(scale - 1.) > 1.e-12) + std::cout << "Compute Gamma*chi_0 / (beta*Nc): " << util::print_time() << std::endl; + else + std::cout << "Compute Gamma*chi_0: " << util::print_time() << std::endl; + } Matrix> Gamma("Gamma", size); Matrix> chi_0("chi_0", size); @@ -763,6 +775,12 @@ void BseLatticeSolver::diagonalizeGammaChi0 // Compute \Gamma\chi_0. linalg::matrixop::gemm(Gamma, chi_0, Gamma_chi_0); + // Scale: divide by (beta * Nc) for PP Q=0 to match Python's diagonalizePPKernel(). + if (std::abs(scale - 1.) > 1.e-12) + for (int j = 0; j < size; j++) + for (int i = 0; i < size; i++) + Gamma_chi_0(i, j) *= scale; + if (concurrency.id() == concurrency.first()) { std::cout << "Finished: " << util::print_time() << std::endl; std::cout << "Diagonalize Gamma*chi_0: " << util::print_time() << std::endl; From d153622eb884f0e604ac273c1c392909f18039f7 Mon Sep 17 00:00:00 2001 From: Peter Doak <11856692+PDoakORNL@users.noreply.github.com> Date: Wed, 17 Jun 2026 19:35:58 +0000 Subject: [PATCH 7/7] Update Tc tutorial for direct-grid chi0 option Update the Tc tutorial input template and regenerate preconfigured inputs to use the new direct-grid chi0 lattice calculation. Add documentation for the direct-grid-chi0 and direct-grid-nkfine analysis parameters. - tutorials/tc/input_tp.json.in: add direct-grid-chi0: true and direct-grid-nkfine: 32 to analysis section - tutorials/tc/preconfigured/T=0.07/0.08/0.09/0.1/input_tp.json: regenerated from template via gen_temps.awk - misc/wiki/Parameters.md: document the two new parameters Reference outputs are intentionally not updated; they will be regenerated separately when the tutorial is next validated end-to-end. --- misc/wiki/Parameters.md | 6 ++++++ tutorials/tc/input_tp.json.in | 4 +++- tutorials/tc/preconfigured/T=0.07/input_tp.json | 4 +++- tutorials/tc/preconfigured/T=0.08/input_tp.json | 4 +++- tutorials/tc/preconfigured/T=0.09/input_tp.json | 4 +++- tutorials/tc/preconfigured/T=0.1/input_tp.json | 4 +++- 6 files changed, 21 insertions(+), 5 deletions(-) diff --git a/misc/wiki/Parameters.md b/misc/wiki/Parameters.md index 361323dcf..b9a128cab 100644 --- a/misc/wiki/Parameters.md +++ b/misc/wiki/Parameters.md @@ -681,6 +681,12 @@ Project Γ χ0 onto the crystal harmonic functions and diagon `"projection-cut-off-radius":` double (1.5) For the projection use crystal-harmonic functions of lattice vectors **r** with ||**r**||2 < *projection-cut-off-radius*. +`"direct-grid-chi0":` boolean (false) +If true, compute the lattice bare susceptibility χ0 and the screened interaction W on a uniform fine k-grid using direct evaluation of the Hamiltonian and cluster self-energy, instead of the default tetrahedron-mesh coarse-graining. This method is required to match the Python reference implementation for the particle-particle channel. + +`"direct-grid-nkfine":` integer (32) +Number of uniform grid points per dimension for the direct-grid χ0 integration. For a 2D system the total number of points per DCA patch is direct-grid-nkfine2. + #### Example diff --git a/tutorials/tc/input_tp.json.in b/tutorials/tc/input_tp.json.in index c0449e6c5..e51c25a15 100644 --- a/tutorials/tc/input_tp.json.in +++ b/tutorials/tc/input_tp.json.in @@ -88,6 +88,8 @@ "analysis": { "symmetrize-Gamma": true, "Gamma-deconvolution-cut-off": 0.5, - "project-onto-crystal-harmonics": false + "project-onto-crystal-harmonics": false, + "direct-grid-chi0": true, + "direct-grid-nkfine": 32 } } diff --git a/tutorials/tc/preconfigured/T=0.07/input_tp.json b/tutorials/tc/preconfigured/T=0.07/input_tp.json index 9e62aac57..03e8a9a87 100644 --- a/tutorials/tc/preconfigured/T=0.07/input_tp.json +++ b/tutorials/tc/preconfigured/T=0.07/input_tp.json @@ -88,6 +88,8 @@ "analysis": { "symmetrize-Gamma": true, "Gamma-deconvolution-cut-off": 0.5, - "project-onto-crystal-harmonics": false + "project-onto-crystal-harmonics": false, + "direct-grid-chi0": true, + "direct-grid-nkfine": 32 } } diff --git a/tutorials/tc/preconfigured/T=0.08/input_tp.json b/tutorials/tc/preconfigured/T=0.08/input_tp.json index 27eb1e01d..74f3d1a61 100644 --- a/tutorials/tc/preconfigured/T=0.08/input_tp.json +++ b/tutorials/tc/preconfigured/T=0.08/input_tp.json @@ -88,6 +88,8 @@ "analysis": { "symmetrize-Gamma": true, "Gamma-deconvolution-cut-off": 0.5, - "project-onto-crystal-harmonics": false + "project-onto-crystal-harmonics": false, + "direct-grid-chi0": true, + "direct-grid-nkfine": 32 } } diff --git a/tutorials/tc/preconfigured/T=0.09/input_tp.json b/tutorials/tc/preconfigured/T=0.09/input_tp.json index 81d2ec6d7..fe76d8325 100644 --- a/tutorials/tc/preconfigured/T=0.09/input_tp.json +++ b/tutorials/tc/preconfigured/T=0.09/input_tp.json @@ -88,6 +88,8 @@ "analysis": { "symmetrize-Gamma": true, "Gamma-deconvolution-cut-off": 0.5, - "project-onto-crystal-harmonics": false + "project-onto-crystal-harmonics": false, + "direct-grid-chi0": true, + "direct-grid-nkfine": 32 } } diff --git a/tutorials/tc/preconfigured/T=0.1/input_tp.json b/tutorials/tc/preconfigured/T=0.1/input_tp.json index 509f26989..13496e9f8 100644 --- a/tutorials/tc/preconfigured/T=0.1/input_tp.json +++ b/tutorials/tc/preconfigured/T=0.1/input_tp.json @@ -88,6 +88,8 @@ "analysis": { "symmetrize-Gamma": true, "Gamma-deconvolution-cut-off": 0.5, - "project-onto-crystal-harmonics": false + "project-onto-crystal-harmonics": false, + "direct-grid-chi0": true, + "direct-grid-nkfine": 32 } }