diff --git a/include/gauxc/c/xc_integrator.h b/include/gauxc/c/xc_integrator.h index 76a0095f8..ec15d90da 100644 --- a/include/gauxc/c/xc_integrator.h +++ b/include/gauxc/c/xc_integrator.h @@ -343,10 +343,15 @@ extern void gauxc_integrator_eval_exx_rks( * @param n Number of columns in the density matrix. * @param density_matrix Pointer to the density matrix data. * @param ldp Leading dimension of the density matrix. - * @param t_density_matrix Pointer to the density matrix data. - * @param ldtp Leading dimension of the density matrix. + * @param t_density_matrix Pointer to the trial density matrix data. + * @param ldtp Leading dimension of the trial density matrix. * @param fxc Matrix container to store the FXC contraction result. * @param ldfxc Leading dimension of the FXC contraction matrix. + * @note Exc depends only on the symmetric part of a density matrix, so the + * kernel is contracted with the symmetric part (tP + tP^T)/2 of each + * trial density tP; any antisymmetric part contributes exactly zero. + * For linear response, pass the perturbed density P1 = T + T^T: + * passing a one-sided transition density T yields half that response. */ extern void gauxc_integrator_eval_fxc_contraction_rks( GauXCStatus* status, @@ -371,14 +376,19 @@ extern void gauxc_integrator_eval_fxc_contraction_rks( * @param ldp_s Leading dimension of the total density matrix. * @param density_matrix_z Pointer to the density matrix data for spin density. * @param ldp_z Leading dimension of the spin density matrix. - * @param t_density_matrix_s Pointer to the density matrix data for total density. - * @param ldtp_s Leading dimension of the total density matrix. - * @param t_density_matrix_z Pointer to the density matrix data for spin density. - * @param ldtp_z Leading dimension of the spin density matrix. + * @param t_density_matrix_s Pointer to the trial density matrix data for total density. + * @param ldtp_s Leading dimension of the total trial density matrix. + * @param t_density_matrix_z Pointer to the trial density matrix data for spin density. + * @param ldtp_z Leading dimension of the spin trial density matrix. * @param fxc_s Matrix container to store the FXC contraction result for total density. * @param ldfxc_s Leading dimension of the FXC contraction matrix for total density. * @param fxc_z Matrix container to store the FXC contraction result for spin density. * @param ldfxc_z Leading dimension of the FXC contraction matrix for spin density. + * @note Exc depends only on the symmetric part of a density matrix, so the + * kernel is contracted with the symmetric part (tP + tP^T)/2 of each + * trial density tP; any antisymmetric part contributes exactly zero. + * For linear response, pass the perturbed density P1 = T + T^T: + * passing a one-sided transition density T yields half that response. */ extern void gauxc_integrator_eval_fxc_contraction_uks( GauXCStatus* status, diff --git a/include/gauxc/xc_integrator.hpp b/include/gauxc/xc_integrator.hpp index 03feaf934..8dd8eed65 100644 --- a/include/gauxc/xc_integrator.hpp +++ b/include/gauxc/xc_integrator.hpp @@ -80,6 +80,20 @@ class XCIntegrator { exx_type eval_exx ( const MatrixType&, const IntegratorSettingsEXX& = IntegratorSettingsEXX{} ); + /** + * @brief Contract the XC kernel with a trial density matrix. + * + * RKS: eval_fxc_contraction( P, tP ). + * UKS: eval_fxc_contraction( Ps, Pz, tPs, tPz ), returning ( FXCs, FXCz ). + * + * Exc depends only on the symmetric part of a density matrix, so the + * kernel is contracted with the symmetric part (tP + tP^T)/2 of each + * trial density; any antisymmetric part contributes exactly zero. + * + * In particular, for linear response the trial density should be the + * perturbed density P1 = T + T^T. Passing a one-sided transition density + * T instead yields half of that response. + */ fxc_contraction_type_rks eval_fxc_contraction ( const MatrixType&, const MatrixType&, const IntegratorSettingsXC& = IntegratorSettingsXC{} ); fxc_contraction_type_uks eval_fxc_contraction ( const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, diff --git a/src/xc_integrator/local_work_driver/host/local_host_work_driver.cxx b/src/xc_integrator/local_work_driver/host/local_host_work_driver.cxx index 0fa970bad..57c04e2c3 100644 --- a/src/xc_integrator/local_work_driver/host/local_host_work_driver.cxx +++ b/src/xc_integrator/local_work_driver/host/local_host_work_driver.cxx @@ -439,14 +439,14 @@ void LocalHostWorkDriver::eval_tmat_mgga_vxc_rks( size_t npts, const double* vga const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, const double* v2gamma2, const double* v2gammalapl, const double* v2gammatau, const double* v2lapl2, const double* v2lapltau, const double* v2tau2, - const double* tden_eval, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C) { + const double* tden_eval, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D) { throw_if_invalid_pimpl(pimpl_); pimpl_->eval_tmat_mgga_vxc_rks(npts, vgamma, v2rho2, v2rhogamma, v2rholapl, v2rhotau, v2gamma2, v2gammalapl, v2gammatau, v2lapl2, v2lapltau, v2tau2, - tden_eval, tdden_x_eval, tdden_y_eval, tdden_z_eval, ttau, dden_x_eval, - dden_y_eval, dden_z_eval, A, B, C); + tden_eval, tdden_x_eval, tdden_y_eval, tdden_z_eval, ttau, tlapl, dden_x_eval, + dden_y_eval, dden_z_eval, A, B, C, D); } @@ -471,14 +471,14 @@ void LocalHostWorkDriver::eval_tmat_mgga_vxc_uks( size_t npts, const double* vga const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, const double* v2gamma2, const double* v2gammalapl, const double* v2gamma_tau, const double* v2lapl2, const double* v2tau_lapl, const double* v2tau2, - const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C) { + const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D) { throw_if_invalid_pimpl(pimpl_); pimpl_->eval_tmat_mgga_vxc_uks(npts, vgamma, v2rho2, v2rhogamma, v2rholapl, v2rhotau, v2gamma2, v2gammalapl, v2gamma_tau, v2lapl2, v2tau_lapl, v2tau2, - trho, tdden_x_eval, tdden_y_eval, tdden_z_eval, ttau, dden_x_eval, - dden_y_eval, dden_z_eval, A, B, C); + trho, tdden_x_eval, tdden_y_eval, tdden_z_eval, ttau, tlapl, dden_x_eval, + dden_y_eval, dden_z_eval, A, B, C, D); } diff --git a/src/xc_integrator/local_work_driver/host/local_host_work_driver.hpp b/src/xc_integrator/local_work_driver/host/local_host_work_driver.hpp index 41cf430ec..228c6ec3f 100644 --- a/src/xc_integrator/local_work_driver/host/local_host_work_driver.hpp +++ b/src/xc_integrator/local_work_driver/host/local_host_work_driver.hpp @@ -554,14 +554,14 @@ class LocalHostWorkDriver : public LocalWorkDriver { const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, const double* v2gamma2, const double* v2gammalapl, const double* v2gammatau, const double* v2lapl2, const double* v2lapltau, const double* v2tau2, - const double* tden_eval, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C); + const double* tden_eval, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D); void eval_tmat_mgga_vxc_uks( size_t npts, const double* vgamma, const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, const double* v2gamma2, const double* v2gammalapl, const double* v2gamma_tau, const double* v2lapl2, const double* v2tau_lapl, const double* v2tau2, - const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C); + const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D); void eval_zmat_lda_vxc_uks_ts( size_t npts, size_t nbe, const double* vrho, diff --git a/src/xc_integrator/local_work_driver/host/local_host_work_driver_pimpl.hpp b/src/xc_integrator/local_work_driver/host/local_host_work_driver_pimpl.hpp index c5e4182d7..556dbf15f 100644 --- a/src/xc_integrator/local_work_driver/host/local_host_work_driver_pimpl.hpp +++ b/src/xc_integrator/local_work_driver/host/local_host_work_driver_pimpl.hpp @@ -193,14 +193,14 @@ struct LocalHostWorkDriverPIMPL { const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, const double* v2gamma2, const double* v2gammalapl, const double* v2gammatau, const double* v2lapl2, const double* v2lapltau, const double* v2tau2, - const double* tden_eval, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C) = 0; + const double* tden_eval, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D) = 0; virtual void eval_tmat_mgga_vxc_uks( size_t npts, const double* vgamma, const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, const double* v2gamma2, const double* v2gammalapl, const double* v2gamma_tau, const double* v2lapl2, const double* v2tau_lapl, const double* v2tau2, - const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C) = 0; + const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D) = 0; virtual void eval_zmat_lda_vxc_uks_ts( size_t npts, size_t nbe, const double* vrho, const double* basis_eval, double* Za, size_t ldza, double* Zb, size_t ldzb ) = 0; diff --git a/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.cxx b/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.cxx index 30429432c..82f37cd90 100644 --- a/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.cxx +++ b/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.cxx @@ -1222,11 +1222,11 @@ void ReferenceLocalHostWorkDriver::eval_tmat_gga_vxc_uks( size_t npts, const dou void ReferenceLocalHostWorkDriver::eval_tmat_mgga_vxc_rks( size_t npts, const double* vgamma, - const double* v2rho2, const double* v2rhogamma, [[maybe_unused]] const double* v2rholapl, const double* v2rhotau, - const double* v2gamma2, [[maybe_unused]] const double* v2gammalapl, const double* v2gammatau, - [[maybe_unused]] const double* v2lapl2, [[maybe_unused]] const double* v2lapltau, const double* v2tau2, - const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C){ + const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, + const double* v2gamma2, const double* v2gammalapl, const double* v2gammatau, + const double* v2lapl2, const double* v2lapltau, const double* v2tau2, + const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D){ for( int32_t i = 0; i < (int32_t)npts; ++i ) { @@ -1237,6 +1237,17 @@ void ReferenceLocalHostWorkDriver::eval_tmat_mgga_vxc_rks( size_t npts, const do C[i] = v2rhotau[i] * trho[i] + 2 * v2gammatau[i] * tgamma + v2tau2[i] * ttau[i]; auto B_coef = v2rhogamma[i] * trho[i] + 2 * v2gamma2[i] * tgamma + v2gammatau[i] * ttau[i]; + + // Laplacian-dependent functionals: the trial Laplacian enters every + // channel, and D = delta v_lapl is a fourth one + if( v2rholapl != nullptr ) { + const auto tl = tlapl[i]; + A[i] += v2rholapl[i] * tl; + C[i] += v2lapltau[i] * tl; + B_coef += v2gammalapl[i] * tl; + D[i] = v2rholapl[i] * trho[i] + 2 * v2gammalapl[i] * tgamma + v2lapl2[i] * tl + + v2lapltau[i] * ttau[i]; + } B[i * 3] = 2 * B_coef * dden_x_eval[i] + 2 * vgamma[i] * tdden_x_eval[i]; B[i * 3 + 1] = 2 * B_coef * dden_y_eval[i] + 2 * vgamma[i] * tdden_y_eval[i]; @@ -1251,12 +1262,10 @@ void ReferenceLocalHostWorkDriver::eval_tmat_mgga_vxc_uks( size_t npts, const do const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, const double* v2gamma2, const double* v2gammalapl, const double* v2gammatau, const double* v2lapl2, const double* v2lapltau, const double* v2tau2, - const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C){ + const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D){ - // Laplacian is not supported now - if( v2rholapl != nullptr || v2gammalapl != nullptr || v2lapltau != nullptr || v2lapl2 != nullptr ) - GAUXC_GENERIC_EXCEPTION(std::string("Laplacian not supported")); + const bool lapl = v2rholapl != nullptr; for( int32_t i = 0; i < (int32_t)npts; ++i ) { @@ -1328,19 +1337,55 @@ void ReferenceLocalHostWorkDriver::eval_tmat_mgga_vxc_uks( size_t npts, const do C[2 * i + 1] = v2rhotau_b_b * trho_b + 2 * v2gammatau_bb_b * tgamma_bb + v2gammatau_ab_b * tgamma_ab + v2tau2_b_b * ttau_b + v2rhotau_a_b * trho_a + 2 * v2gammatau_aa_b * tgamma_aa + v2tau2_a_b * ttau_a; + // Laplacian-dependent functionals. Libxc orders the mixed blocks with + // the non-Laplacian variable first (rho_x lapl_y, sigma_xy lapl_z) + // and puts lapl before tau (lapl_x tau_y). + double tl_a = 0., tl_b = 0.; + double v2gammalapl_aa_a = 0., v2gammalapl_aa_b = 0., v2gammalapl_ab_a = 0., + v2gammalapl_ab_b = 0., v2gammalapl_bb_a = 0., v2gammalapl_bb_b = 0.; + if( lapl ) { + tl_a = tlapl[2*i]; tl_b = tlapl[2*i+1]; + const auto v2rholapl_a_a = v2rholapl[4*i], v2rholapl_a_b = v2rholapl[4*i+1]; + const auto v2rholapl_b_a = v2rholapl[4*i+2], v2rholapl_b_b = v2rholapl[4*i+3]; + v2gammalapl_aa_a = v2gammalapl[6*i]; v2gammalapl_aa_b = v2gammalapl[6*i+1]; + v2gammalapl_ab_a = v2gammalapl[6*i+2]; v2gammalapl_ab_b = v2gammalapl[6*i+3]; + v2gammalapl_bb_a = v2gammalapl[6*i+4]; v2gammalapl_bb_b = v2gammalapl[6*i+5]; + const auto v2lapl2_a_a = v2lapl2[3*i], v2lapl2_a_b = v2lapl2[3*i+1], v2lapl2_b_b = v2lapl2[3*i+2]; + const auto v2lapltau_a_a = v2lapltau[4*i], v2lapltau_a_b = v2lapltau[4*i+1]; + const auto v2lapltau_b_a = v2lapltau[4*i+2], v2lapltau_b_b = v2lapltau[4*i+3]; + + A[2*i] += v2rholapl_a_a * tl_a + v2rholapl_a_b * tl_b; + A[2*i+1] += v2rholapl_b_a * tl_a + v2rholapl_b_b * tl_b; + C[2*i] += v2lapltau_a_a * tl_a + v2lapltau_b_a * tl_b; + C[2*i+1] += v2lapltau_a_b * tl_a + v2lapltau_b_b * tl_b; + + D[2*i] = v2rholapl_a_a * trho_a + v2rholapl_b_a * trho_b + + 2 * v2gammalapl_aa_a * tgamma_aa + v2gammalapl_ab_a * tgamma_ab + 2 * v2gammalapl_bb_a * tgamma_bb + + v2lapl2_a_a * tl_a + v2lapl2_a_b * tl_b + + v2lapltau_a_a * ttau_a + v2lapltau_a_b * ttau_b; + D[2*i+1] = v2rholapl_a_b * trho_a + v2rholapl_b_b * trho_b + + 2 * v2gammalapl_aa_b * tgamma_aa + v2gammalapl_ab_b * tgamma_ab + 2 * v2gammalapl_bb_b * tgamma_bb + + v2lapl2_a_b * tl_a + v2lapl2_b_b * tl_b + + v2lapltau_b_a * ttau_a + v2lapltau_b_b * ttau_b; + } + auto B_coef1 = v2rhogamma_a_aa * trho_a + 2 * v2gamma2_aa_aa * tgamma_aa + v2gamma2_aa_ab * tgamma_ab + v2gammatau_aa_a * ttau_a - + v2rhogamma_b_aa * trho_b + 2 * v2gamma2_aa_bb * tgamma_bb + v2gammatau_aa_b * ttau_b; + + v2rhogamma_b_aa * trho_b + 2 * v2gamma2_aa_bb * tgamma_bb + v2gammatau_aa_b * ttau_b + + v2gammalapl_aa_a * tl_a + v2gammalapl_aa_b * tl_b; auto B_coef2 = v2rhogamma_a_ab * trho_a + 2 * v2gamma2_aa_ab * tgamma_aa + v2gamma2_ab_ab * tgamma_ab + v2gammatau_ab_a * ttau_a - + v2rhogamma_b_ab * trho_b + 2 * v2gamma2_ab_bb * tgamma_bb + v2gammatau_ab_b * ttau_b; + + v2rhogamma_b_ab * trho_b + 2 * v2gamma2_ab_bb * tgamma_bb + v2gammatau_ab_b * ttau_b + + v2gammalapl_ab_a * tl_a + v2gammalapl_ab_b * tl_b; B[i * 6] = 2 * B_coef1 * dden_x_eval_a + B_coef2 * dden_x_eval_b + 2 * vgamma_aa * tdden_x_eval_a + vgamma_ab * tdden_x_eval_b; B[i * 6 + 1] = 2 * B_coef1 * dden_y_eval_a + B_coef2 * dden_y_eval_b + 2 * vgamma_aa * tdden_y_eval_a + vgamma_ab * tdden_y_eval_b; B[i * 6 + 2] = 2 * B_coef1 * dden_z_eval_a + B_coef2 * dden_z_eval_b + 2 * vgamma_aa * tdden_z_eval_a + vgamma_ab * tdden_z_eval_b; B_coef1 = v2rhogamma_b_bb * trho_b + 2 * v2gamma2_bb_bb * tgamma_bb + v2gamma2_ab_bb * tgamma_ab + v2gammatau_bb_b * ttau_b - + v2rhogamma_a_bb * trho_a + 2 * v2gamma2_aa_bb * tgamma_aa + v2gammatau_bb_a * ttau_a; + + v2rhogamma_a_bb * trho_a + 2 * v2gamma2_aa_bb * tgamma_aa + v2gammatau_bb_a * ttau_a + + v2gammalapl_bb_a * tl_a + v2gammalapl_bb_b * tl_b; B_coef2 = v2rhogamma_b_ab * trho_b + 2 * v2gamma2_ab_bb * tgamma_bb + v2gamma2_ab_ab * tgamma_ab + v2gammatau_ab_b * ttau_b - + v2rhogamma_a_ab * trho_a + 2 * v2gamma2_aa_ab * tgamma_aa + v2gammatau_ab_a * ttau_a; + + v2rhogamma_a_ab * trho_a + 2 * v2gamma2_aa_ab * tgamma_aa + v2gammatau_ab_a * ttau_a + + v2gammalapl_ab_a * tl_a + v2gammalapl_ab_b * tl_b; B[i * 6 + 3] = 2 * B_coef1 * dden_x_eval_b + B_coef2 * dden_x_eval_a + 2 * vgamma_bb * tdden_x_eval_b + vgamma_ab * tdden_x_eval_a; B[i * 6 + 4] = 2 * B_coef1 * dden_y_eval_b + B_coef2 * dden_y_eval_a + 2 * vgamma_bb * tdden_y_eval_b + vgamma_ab * tdden_y_eval_a; diff --git a/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.hpp b/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.hpp index 3560b85b9..08ef11f61 100644 --- a/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.hpp +++ b/src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.hpp @@ -198,14 +198,14 @@ struct ReferenceLocalHostWorkDriver : public detail::LocalHostWorkDriverPIMPL { const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, const double* v2gamma2, const double* v2gammalapl, const double* v2gammatau, const double* v2lapl2, const double* v2lapltau, const double* v2tau2, - const double* tden_eval, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C) override; + const double* tden_eval, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D) override; void eval_tmat_mgga_vxc_uks( size_t npts, const double* vgamma, const double* v2rho2, const double* v2rhogamma, const double* v2rholapl, const double* v2rhotau, const double* v2gamma2, const double* v2gammalapl, const double* v2gamma_tau, const double* v2lapl2, const double* v2tau_lapl, const double* v2tau2, - const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, - const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C) override; + const double* trho, const double* tdden_x_eval, const double* tdden_y_eval, const double* tdden_z_eval, const double* ttau, const double* tlapl, + const double* dden_x_eval, const double* dden_y_eval, const double* dden_z_eval, double* A, double* B, double* C, double* D) override; void eval_zmat_lda_vxc_uks_ts( size_t npts, size_t nbe, const double* vrho, const double* basis_eval, double* Za, size_t ldza, double* Zb, size_t ldzb ) override; diff --git a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp index ffc0ca417..f1c1cd70d 100644 --- a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp +++ b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp @@ -13,6 +13,7 @@ #include "incore_replicated_xc_device_integrator.hpp" #include #include +#include namespace GauXC::detail { @@ -66,6 +67,26 @@ namespace GauXC::detail { GAUXC_GENERIC_EXCEPTION("Invalid LDFXCz"); } + // Symmetrize the TRIAL densities -- see the host integrator for why: + // Exc sees only the symmetric part of a density matrix, but the + // gradient channel of the uvvar kernels assumes symmetry rather than + // enforcing it, so an antisymmetric trial direction (a TDDFT + // B-matrix or stability block) produced a nonzero kernel and broke + // the Hermiticity of the Kohn-Sham Hessian for every GGA and + // meta-GGA. https://github.com/wavefunction91/GauXC/issues/225 + std::vector tPs_sym, tPz_sym; + auto symmetrize = [nbf]( const value_type* A, int64_t lda, + std::vector& out ) -> const value_type* { + if( !A ) return nullptr; + out.resize( nbf*nbf ); + for( int64_t j = 0; j < nbf; ++j ) + for( int64_t i = 0; i < nbf; ++i ) + out[i + j*nbf] = 0.5 * ( A[i + j*lda] + A[j + i*lda] ); + return out.data(); + }; + tPs = symmetrize( tPs, ldtps, tPs_sym ); if( tPs ) ldtps = nbf; + tPz = symmetrize( tPz, ldtpz, tPz_sym ); if( tPz ) ldtpz = nbf; + // Get Tasks auto& tasks = this->load_balancer_->get_tasks(); diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp index 7320039c7..6d9d5081a 100644 --- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp @@ -16,6 +16,7 @@ #include "host/local_host_work_driver.hpp" #include "host/blas.hpp" #include +#include namespace GauXC::detail { @@ -57,6 +58,43 @@ void ReferenceReplicatedXCHostIntegrator:: GAUXC_GENERIC_EXCEPTION("Invalid LDFXCZ"); + // Symmetrize the TRIAL densities. + // + // The XC energy is a functional of n = Tr[P B B], which sees only the + // symmetric part of P: the derivative of Exc along an antisymmetric + // direction is identically zero. The density channel gets that for + // free -- Tr[P B B] kills an antisymmetric P by itself -- but the + // GRADIENT channel does not. eval_uvvar_* builds it as + // d_x n = 2 X . d_x B, X = fac * P * B, + // which equals Tr[(P + P^T) B d_x B] only when P is symmetric. That + // holds for the reference density, so exc_vxc and exc_grad are + // unaffected; it does NOT hold for a trial density, which a caller has + // no reason to hand over symmetric -- a TDDFT B-matrix or stability + // block is contracted against the ANTISYMMETRIC combination on + // purpose. The result was an XC kernel that answered nonzero to a + // direction along which Exc does not vary, breaking the Hermiticity + // of the Kohn-Sham Hessian for every GGA and meta-GGA while leaving + // LDA (which never touches the gradient channel) looking correct. + // See https://github.com/wavefunction91/GauXC/issues/225. + // + // Symmetrizing here rather than in the uvvar kernels keeps the fix in + // one place for every backend and every rung, and costs one nbf^2 + // copy against a grid loop. + std::vector tPs_sym, tPz_sym; + auto symmetrize = [nbf]( const value_type* A, int64_t lda, + std::vector& out ) -> const value_type* { + if( !A ) return nullptr; + out.resize( nbf*nbf ); + for( int64_t j = 0; j < nbf; ++j ) + for( int64_t i = 0; i < nbf; ++i ) + out[i + j*nbf] = 0.5 * ( A[i + j*lda] + A[j + i*lda] ); + return out.data(); + }; + const value_type* tPs_use = symmetrize( tPs, ldtps, tPs_sym ); + const value_type* tPz_use = symmetrize( tPz, ldtpz, tPz_sym ); + const int64_t ldtps_use = tPs_use ? nbf : ldtps; + const int64_t ldtpz_use = tPz_use ? nbf : ldtpz; + // Get Tasks auto& tasks = this->load_balancer_->get_tasks(); @@ -66,7 +104,7 @@ void ReferenceReplicatedXCHostIntegrator:: // Compute Local contributions to FXC contraction this->timer_.time_op("XCIntegrator.LocalWork", [&](){ fxc_contraction_local_work_( basis, Ps, ldps, Pz, ldpz, - tPs, ldtps, tPz, ldtpz, + tPs_use, ldtps_use, tPz_use, ldtpz_use, FXCs, ldfxcs, FXCz, ldfxcz, &N_EL, ks_settings, tasks.begin(), tasks.end() ); @@ -117,10 +155,6 @@ void ReferenceReplicatedXCHostIntegrator:: const auto& mol = this->load_balancer_->molecule(); const bool needs_laplacian = func.needs_laplacian(); - // not suppport laplacian yet - if( needs_laplacian ) { - GAUXC_GENERIC_EXCEPTION("Laplacian Not Supported Yet for FXC Contraction"); - } // Get basis map BasisSetMap basis_map(basis,mol); @@ -254,6 +288,7 @@ void ReferenceReplicatedXCHostIntegrator:: host_data.v2gammalapl.resize(npts * spin_dim_rhogamma); host_data.v2lapltau .resize(npts * spin_dim_rhotau); host_data.tlapl .resize(npts * spin_dim_scal); + host_data.FXC_D .resize(npts * spin_dim_scal); } else { host_data.basis_eval .resize( 4 * npts * nbe ); // basis + grad (3) @@ -298,6 +333,13 @@ void ReferenceReplicatedXCHostIntegrator:: auto* FXC_A = host_data.FXC_A.data(); auto* FXC_B = host_data.FXC_B.data(); auto* FXC_C = host_data.FXC_C.data(); + // delta v_lapl; null for functionals without a Laplacian, which is + // what switches the Laplacian terms off downstream + auto* FXC_D = needs_laplacian ? host_data.FXC_D.data() : nullptr; + auto* v2rholapl_t = needs_laplacian ? v2rholapl : nullptr; + auto* v2gammalapl_t = needs_laplacian ? v2gammalapl : nullptr; + auto* v2lapl2_t = needs_laplacian ? v2lapl2 : nullptr; + auto* v2lapltau_t = needs_laplacian ? v2lapltau : nullptr; value_type* dbasis_x_eval = nullptr; @@ -453,9 +495,9 @@ void ReferenceReplicatedXCHostIntegrator:: lwd->eval_uvvar_mgga_rks( npts, nbe, basis_eval, dbasis_x_eval, dbasis_y_eval, dbasis_z_eval, lbasis_eval, zmat, nbe, mmat_x, mmat_y, mmat_z, nbe, tden_eval, tdden_x_eval, tdden_y_eval, tdden_z_eval, gamma, ttau, tlapl); - lwd->eval_tmat_mgga_vxc_rks( npts, vgamma, v2rho2, v2rhogamma, v2rholapl, v2rhotau, v2gamma2, - v2gammalapl, v2gammatau, v2lapl2, v2lapltau, v2tau2, tden_eval, tdden_x_eval, - tdden_y_eval, tdden_z_eval, ttau, dden_x_eval, dden_y_eval, dden_z_eval, FXC_A, FXC_B, FXC_C ); + lwd->eval_tmat_mgga_vxc_rks( npts, vgamma, v2rho2, v2rhogamma, v2rholapl_t, v2rhotau, v2gamma2, + v2gammalapl_t, v2gammatau, v2lapl2_t, v2lapltau_t, v2tau2, tden_eval, tdden_x_eval, + tdden_y_eval, tdden_z_eval, ttau, tlapl, dden_x_eval, dden_y_eval, dden_z_eval, FXC_A, FXC_B, FXC_C, FXC_D ); } else if (is_uks) { // tgamma is not needed since it has different definitions than gamma // gamma = nabla rho * nabla rho, but tgamma = nabla trho * nabla rho, not both trho @@ -463,9 +505,9 @@ void ReferenceReplicatedXCHostIntegrator:: dbasis_z_eval, lbasis_eval, zmat, nbe, zmat_z, nbe, mmat_x, mmat_y, mmat_z, nbe, mmat_x_z, mmat_y_z, mmat_z_z, nbe, tden_eval, tdden_x_eval, tdden_y_eval, tdden_z_eval, gamma, ttau, tlapl); - lwd->eval_tmat_mgga_vxc_uks( npts, vgamma, v2rho2, v2rhogamma, v2rholapl, v2rhotau, v2gamma2, - v2gammalapl, v2gammatau, v2lapl2, v2lapltau, v2tau2, tden_eval, tdden_x_eval, - tdden_y_eval, tdden_z_eval, ttau, dden_x_eval, dden_y_eval, dden_z_eval, FXC_A, FXC_B, FXC_C ); + lwd->eval_tmat_mgga_vxc_uks( npts, vgamma, v2rho2, v2rhogamma, v2rholapl_t, v2rhotau, v2gamma2, + v2gammalapl_t, v2gammatau, v2lapl2_t, v2lapltau_t, v2tau2, tden_eval, tdden_x_eval, + tdden_y_eval, tdden_z_eval, ttau, tlapl, dden_x_eval, dden_y_eval, dden_z_eval, FXC_A, FXC_B, FXC_C, FXC_D ); } } else if ( func.is_gga() ) { if(is_rks) { @@ -517,6 +559,12 @@ void ReferenceReplicatedXCHostIntegrator:: if(not is_rks) FXC_C[sds*i+1] *= weights[i]; } } + if( FXC_D ){ + for( int32_t i = 0; i < npts; ++i) { + FXC_D[sds*i] *= weights[i]; + if(not is_rks) FXC_D[sds*i+1] *= weights[i]; + } + } // Scalar integrations double NEL_local = 0.0; @@ -532,16 +580,25 @@ void ReferenceReplicatedXCHostIntegrator:: // Evaluate Z matrix for VXC if( func.is_mgga() ) { if(is_rks) { - // Because we do not support Laplacian, so mgga will do the same operation as GGA + // The GGA part, then delta v_lapl exactly where the ground-state + // potential puts v_lapl: D lapl(chi) in Z, D grad(chi) in M lwd->eval_zmat_gga_vxc_rks_ts( npts, nbe, FXC_A, FXC_B, basis_eval, dbasis_x_eval, dbasis_y_eval, dbasis_z_eval, zmat, nbe); - lwd->eval_mmat_mgga_vxc_rks( npts, nbe, FXC_C, vlapl, dbasis_x_eval, dbasis_y_eval, dbasis_z_eval, + if( FXC_D ) + for( int32_t i = 0; i < npts; ++i ) + blas::axpy( nbe, FXC_D[i], lbasis_eval + size_t(i)*nbe, 1, zmat + size_t(i)*nbe, 1 ); + lwd->eval_mmat_mgga_vxc_rks( npts, nbe, FXC_C, FXC_D, dbasis_x_eval, dbasis_y_eval, dbasis_z_eval, mmat_x, mmat_y, mmat_z, nbe); } else if (is_uks) { - // Because we do not support Laplacian, so mgga will do the same operation as GGA + // as for RKS, per spin (alpha in Z, beta in Z_z at this stage) lwd->eval_zmat_gga_vxc_uks_ts( npts, nbe, FXC_A, FXC_B, basis_eval, dbasis_x_eval, dbasis_y_eval, dbasis_z_eval, zmat, nbe, zmat_z, nbe); - lwd->eval_mmat_mgga_vxc_uks_ts( npts, nbe, FXC_C, vlapl, dbasis_x_eval, dbasis_y_eval, dbasis_z_eval, + if( FXC_D ) + for( int32_t i = 0; i < npts; ++i ) { + blas::axpy( nbe, FXC_D[2*i], lbasis_eval + size_t(i)*nbe, 1, zmat + size_t(i)*nbe, 1 ); + blas::axpy( nbe, FXC_D[2*i+1], lbasis_eval + size_t(i)*nbe, 1, zmat_z + size_t(i)*nbe, 1 ); + } + lwd->eval_mmat_mgga_vxc_uks_ts( npts, nbe, FXC_C, FXC_D, dbasis_x_eval, dbasis_y_eval, dbasis_z_eval, mmat_x, mmat_y, mmat_z, nbe, mmat_x_z, mmat_y_z, mmat_z_z, nbe); } } diff --git a/src/xc_integrator/replicated/host/xc_host_data.hpp b/src/xc_integrator/replicated/host/xc_host_data.hpp index 1c7fc9a21..58e50c1ad 100644 --- a/src/xc_integrator/replicated/host/xc_host_data.hpp +++ b/src/xc_integrator/replicated/host/xc_host_data.hpp @@ -51,6 +51,7 @@ struct XCHostData { std::vector FXC_A; std::vector FXC_B; std::vector FXC_C; + std::vector FXC_D; std::vector tden_scr; std::vector ttau; std::vector tlapl; diff --git a/tests/2nd_derivative_test.cxx b/tests/2nd_derivative_test.cxx index eedcb27f7..f7bd8ddcc 100644 --- a/tests/2nd_derivative_test.cxx +++ b/tests/2nd_derivative_test.cxx @@ -134,6 +134,53 @@ void test_fxc_contractioin(ExecutionSpace ex, const RuntimeEnvironment& rt, CHECK(FXCz_diff_nrm / basis.nbf() < 1e-10); } + + // The stored trial densities are symmetric, so the comparison above + // cannot see how a non-symmetric trial density (e.g. a one-sided + // transition density) is handled. Exc depends only on the symmetric + // part of a density matrix, so the contraction must act on (M + M^T)/2 + // and vanish identically along (M - M^T)/2, at every rung. The gradient + // channel of GGAs and meta-GGAs used to violate both, see + // https://github.com/wavefunction91/GauXC/issues/225. + // + // M is built so that (M + M^T)/2 equals the stored trial density exactly + // (strict upper triangle doubled, lower triangle zero), so F[M] can be + // checked against the same reference data. + auto make_nonsymmetric = []( const matrix_type& A ) -> matrix_type { + matrix_type M = 2.0 * A.triangularView().toDenseMatrix(); + M.diagonal() = A.diagonal(); + return M; + }; + auto antisymmetric_part = []( const matrix_type& M ) -> matrix_type { + return 0.5 * ( M - M.transpose() ); + }; + + if (rks) { + matrix_type tM = make_nonsymmetric(tP); + matrix_type tA = antisymmetric_part(tM); + REQUIRE( tA.norm() > 0.1 * tP.norm() ); + + auto FXC_M = integrator.eval_fxc_contraction(P, tM); + CHECK((FXC_M - FXC_ref).norm() / basis.nbf() < 1e-10); + + auto FXC_A = integrator.eval_fxc_contraction(P, tA); + CHECK(FXC_A.norm() / basis.nbf() < 1e-10); + } else if (uks) { + matrix_type tMs = make_nonsymmetric(tP); + matrix_type tMz = make_nonsymmetric(tPz); + matrix_type tAs = antisymmetric_part(tMs); + matrix_type tAz = antisymmetric_part(tMz); + REQUIRE( tAs.norm() > 0.1 * tP.norm() ); + REQUIRE( tAz.norm() > 0.1 * tPz.norm() ); + + auto [FXCs_M, FXCz_M] = integrator.eval_fxc_contraction(P, Pz, tMs, tMz); + CHECK((FXCs_M - FXC_ref ).norm() / basis.nbf() < 1e-10); + CHECK((FXCz_M - FXCz_ref).norm() / basis.nbf() < 1e-10); + + auto [FXCs_A, FXCz_A] = integrator.eval_fxc_contraction(P, Pz, tAs, tAz); + CHECK(FXCs_A.norm() / basis.nbf() < 1e-10); + CHECK(FXCz_A.norm() / basis.nbf() < 1e-10); + } } void test_integrator_2nd(std::string reference_file, functional_type& func, PruningScheme pruning_scheme) {