fxc contraction: Laplacian-dependent meta-GGAs (RKS and UKS)

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,

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&,

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);
 
 }
 

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,

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;

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;