XC nuclear Hessian: basis, moving-grid and partition-weight classes

include/gauxc/xc_integrator.hpp

@@ -39,6 +39,8 @@ class XCIntegrator {
   using exc_vxc_type_uks  = std::tuple< value_type, matrix_type, matrix_type >;  
   using exc_vxc_type_gks  = std::tuple< value_type, matrix_type, matrix_type, matrix_type, matrix_type >;
   using exc_grad_type = std::vector< value_type >;
+  //! Nuclear Hessian, row-major (3*natoms) x (3*natoms).
+  using exc_hess_type = std::vector< value_type >;
   using exx_type      = matrix_type;
   using fxc_contraction_type_rks = matrix_type;
   using fxc_contraction_type_uks = std::tuple< matrix_type, matrix_type >;
@@ -77,6 +79,8 @@ class XCIntegrator {
   exc_grad_type eval_exc_grad( const MatrixType&, const IntegratorSettingsXC& = IntegratorSettingsXC{} );
   exc_grad_type eval_exc_grad( const MatrixType&, const MatrixType&, const IntegratorSettingsXC& = IntegratorSettingsXC{} );
 
+  exc_hess_type eval_exc_hess( const MatrixType&, const IntegratorSettingsXC& = IntegratorSettingsXC{} );
+
   exx_type      eval_exx     ( const MatrixType&, 
                                const IntegratorSettingsEXX& = IntegratorSettingsEXX{} );
 

include/gauxc/xc_integrator/impl.hpp

@@ -92,6 +92,13 @@ typename XCIntegrator<MatrixType>::exc_grad_type
   return pimpl_->eval_exc_grad(Ps, Pz, ks_settings);
 };
 
+template <typename MatrixType>
+typename XCIntegrator<MatrixType>::exc_hess_type
+  XCIntegrator<MatrixType>::eval_exc_hess( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) {
+  if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED();
+  return pimpl_->eval_exc_hess(P, ks_settings);
+};
+
 template <typename MatrixType>
 typename XCIntegrator<MatrixType>::exx_type
   XCIntegrator<MatrixType>::eval_exx( const MatrixType&     P,

include/gauxc/xc_integrator/replicated/impl.hpp

@@ -175,6 +175,21 @@ typename ReplicatedXCIntegrator<MatrixType>::exc_grad_type
 
 }
 
+template <typename MatrixType>
+typename ReplicatedXCIntegrator<MatrixType>::exc_hess_type
+  ReplicatedXCIntegrator<MatrixType>::eval_exc_hess_( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) {
+
+  if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED();
+
+  const size_t n3 = 3*pimpl_->load_balancer().molecule().natoms();
+  std::vector<value_type> EXC_HESS( n3*n3 );
+  pimpl_->eval_exc_hess( P.rows(), P.cols(), P.data(), P.rows(),
+                         EXC_HESS.data(), ks_settings );
+
+  return EXC_HESS;
+
+}
+
 template <typename MatrixType>
 typename ReplicatedXCIntegrator<MatrixType>::exc_grad_type 
   ReplicatedXCIntegrator<MatrixType>::eval_exc_grad_( const MatrixType& Ps, const MatrixType& Pz, const IntegratorSettingsXC& ks_settings ) {

include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp

@@ -80,6 +80,11 @@ class ReplicatedXCIntegratorImpl {
 
   virtual void eval_exc_grad_( int64_t m, int64_t n, const value_type* P, int64_t ldp, 
                                value_type* EXC_GRAD, const IntegratorSettingsXC& ks_settings ) = 0;
+  virtual void eval_exc_hess_( int64_t m, int64_t n, const value_type* P, int64_t ldp,
+                              value_type* EXC_HESS, const IntegratorSettingsXC& ks_settings ) {
+    GAUXC_GENERIC_EXCEPTION("EXC Hessian Not Implemented For This Integrator");
+  }
+
   virtual void eval_exc_grad_( int64_t m, int64_t n, const value_type* P, int64_t ldps, 
                                const value_type* Pz, int64_t lpdz, value_type* EXC_GRAD, const IntegratorSettingsXC& ks_settings ) = 0;
   virtual void eval_exx_( int64_t m, int64_t n, const value_type* P,
@@ -155,6 +160,11 @@ class ReplicatedXCIntegratorImpl {
 
   void eval_exc_grad( int64_t m, int64_t n, const value_type* P, int64_t ldp, 
                       value_type* EXC_GRAD, const IntegratorSettingsXC& ks_settings );
+  void eval_exc_hess( int64_t m, int64_t n, const value_type* P, int64_t ldp,
+                     value_type* EXC_HESS, const IntegratorSettingsXC& ks_settings ) {
+    eval_exc_hess_(m, n, P, ldp, EXC_HESS, ks_settings);
+  }
+
   void eval_exc_grad( int64_t m, int64_t n, const value_type* Ps, int64_t ldps, 
                       const value_type* Pz, int64_t ldpz, value_type* EXC_GRAD, const IntegratorSettingsXC& ks_settings );
 

include/gauxc/xc_integrator/replicated_xc_integrator.hpp

@@ -36,6 +36,7 @@ class ReplicatedXCIntegrator : public XCIntegratorImpl<MatrixType> {
   using exc_vxc_type_uks   = typename XCIntegratorImpl<MatrixType>::exc_vxc_type_uks;
   using exc_vxc_type_gks   = typename XCIntegratorImpl<MatrixType>::exc_vxc_type_gks;
   using exc_grad_type  = typename XCIntegratorImpl<MatrixType>::exc_grad_type;
+  using exc_hess_type  = typename XCIntegratorImpl<MatrixType>::exc_hess_type;
   using exx_type       = typename XCIntegratorImpl<MatrixType>::exx_type;
   using fxc_contraction_type_rks   = typename XCIntegratorImpl<MatrixType>::fxc_contraction_type_rks;
   using fxc_contraction_type_uks   = typename XCIntegratorImpl<MatrixType>::fxc_contraction_type_uks;
@@ -56,6 +57,7 @@ class ReplicatedXCIntegrator : public XCIntegratorImpl<MatrixType> {
   exc_vxc_type_gks  eval_exc_vxc_ ( const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, const IntegratorSettingsXC& ) override;
   exc_grad_type eval_exc_grad_( const MatrixType&, const IntegratorSettingsXC& ) override;
   exc_grad_type eval_exc_grad_( const MatrixType&, const MatrixType&, const IntegratorSettingsXC& ) override;
+  exc_hess_type eval_exc_hess_( const MatrixType&, const IntegratorSettingsXC& ) override;
   exx_type      eval_exx_     ( const MatrixType&, const IntegratorSettingsEXX& ) override;
   fxc_contraction_type_rks  eval_fxc_contraction_ ( const MatrixType&, const MatrixType&, const IntegratorSettingsXC& ) override;
   fxc_contraction_type_uks  eval_fxc_contraction_ ( const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, const IntegratorSettingsXC&) override;

include/gauxc/xc_integrator/xc_integrator_impl.hpp

@@ -28,6 +28,7 @@ class XCIntegratorImpl {
   using exc_vxc_type_uks   = typename XCIntegrator<MatrixType>::exc_vxc_type_uks;
   using exc_vxc_type_gks   = typename XCIntegrator<MatrixType>::exc_vxc_type_gks;
   using exc_grad_type  = typename XCIntegrator<MatrixType>::exc_grad_type;
+  using exc_hess_type  = typename XCIntegrator<MatrixType>::exc_hess_type;
   using exx_type       = typename XCIntegrator<MatrixType>::exx_type;
   using fxc_contraction_type_rks   = typename XCIntegrator<MatrixType>::fxc_contraction_type_rks;
   using fxc_contraction_type_uks   = typename XCIntegrator<MatrixType>::fxc_contraction_type_uks;
@@ -48,6 +49,12 @@ class XCIntegratorImpl {
                                             const IntegratorSettingsXC& ks_settings ) = 0;
   virtual exc_grad_type eval_exc_grad_( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) = 0;
   virtual exc_grad_type eval_exc_grad_( const MatrixType& Ps, const MatrixType& Pz, const IntegratorSettingsXC& ks_settings ) = 0;
+  /** Nuclear Hessian. Defaulted to a refusal rather than pure virtual:
+   *  only the reference host integrator implements it so far, and every
+   *  other backend would otherwise need a stub. */
+  virtual exc_hess_type eval_exc_hess_( const MatrixType&, const IntegratorSettingsXC& ) {
+    GAUXC_GENERIC_EXCEPTION("EXC Hessian Not Implemented For This Integrator");
+  }
   virtual exx_type      eval_exx_     ( const MatrixType&     P, 
                                         const IntegratorSettingsEXX& settings ) = 0;
   virtual fxc_contraction_type_rks  eval_fxc_contraction_ ( const MatrixType& P,
@@ -142,6 +149,15 @@ class XCIntegratorImpl {
     return eval_exc_grad_(Ps, Pz, ks_settings);
   }
 
+  /** Integrate the EXC nuclear Hessian for RKS
+   *
+   *  @param[in] P The alpha density matrix
+   *  @returns EXC Hessian, row-major (3*natoms) x (3*natoms)
+   */
+  exc_hess_type eval_exc_hess( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) {
+    return eval_exc_hess_(P, ks_settings);
+  }
+
   /** Integrate Exact Exchange for RHF
    *
    *  @param[in] P The alpha density matrix

include/gauxc/xc_integrator_settings.hpp

@@ -29,4 +29,8 @@ struct IntegratorSettingsEXC_GRAD : public IntegratorSettingsKS {
   bool include_weight_derivatives= true; // whether to include grid weight contribution and employ translational invariance, or just use Hellmann-Feynman gradient
 };
 
+struct IntegratorSettingsEXC_HESS : public IntegratorSettingsKS {
+  bool include_weight_derivatives = true; // as for EXC_GRAD: moving grid (translational invariance) plus the partition-weight first and second derivatives; false gives the fixed-grid, fixed-weight (basis-only) Hessian
+};
+
 }

src/xc_integrator/local_work_driver/host/gauxc_hess_call_cell.inc

@@ -0,0 +1,9 @@
+// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: cell-function call site] <==
+// Arguments are bound BY NAME from each kernel's own signature;
+// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir <dir>
+if( is_becke ) {
+  xckernel::xck_gauxc_weight_cell_becke( mu, s, t, u );
+} else {
+  xckernel::xck_gauxc_weight_cell_ssf( mu, s, t, u );
+}
+// ==> END GENERATED CODE <==

src/xc_integrator/local_work_driver/host/gauxc_hess_call_egrad.inc

@@ -0,0 +1,11 @@
+// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: egrad call site] <==
+// Arguments are bound BY NAME from each kernel's own signature;
+// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir <dir>
+if( is_mgga ) {
+  xckernel::xck_gauxc_hess_egrad_mgga_tau( ROW(a,d,0,ip), ROW(a,d,1,ip), ROW(a,d,2,ip), vrho[ip], vgamma[ip], vtau[ip], de );
+} else if( is_gga ) {
+  xckernel::xck_gauxc_hess_egrad_gga( ROW(a,d,0,ip), ROW(a,d,1,ip), vrho[ip], vgamma[ip], de );
+} else {
+  xckernel::xck_gauxc_hess_egrad_lda( ROW(a,d,0,ip), vrho[ip], de );
+}
+// ==> END GENERATED CODE <==

src/xc_integrator/local_work_driver/host/gauxc_hess_call_mu.inc

@@ -0,0 +1,5 @@
+// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: mu call site] <==
+// Arguments are bound BY NAME from each kernel's own signature;
+// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir <dir>
+xckernel::xck_gauxc_weight_mu( RD[0], RD[1], RD[2], RE[0], RE[1], RE[2], rg[0], rg[1], rg[2], mu, dmu[0], dmu[1], dmu[2], dmu[3], dmu[4], dmu[5], d2mu[0], d2mu[1], d2mu[2], d2mu[3], d2mu[4], d2mu[5], d2mu[6], d2mu[7], d2mu[8], d2mu[9], d2mu[10], d2mu[11], d2mu[12], d2mu[13], d2mu[14], d2mu[15], d2mu[16], d2mu[17], d2mu[18], d2mu[19], d2mu[20], d2mu[21], d2mu[22], d2mu[23], d2mu[24], d2mu[25], d2mu[26], d2mu[27], d2mu[28], d2mu[29], d2mu[30], d2mu[31], d2mu[32], d2mu[33], d2mu[34], d2mu[35] );
+// ==> END GENERATED CODE <==

src/xc_integrator/local_work_driver/host/gauxc_hess_call_pair.inc

@@ -0,0 +1,11 @@
+// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: pair call site] <==
+// Arguments are bound BY NAME from each kernel's own signature;
+// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir <dir>
+if( is_mgga ) {
+  xckernel::xck_gauxc_hess_pair_mgga_tau( ROW(a,dx,0,ip), ROW(b,dy,0,ip), ROW(a,dx,1,ip), ROW(b,dy,1,ip), ROW(a,dx,2,ip), ROW(b,dy,2,ip), ROW(a,dx,3,ip), ROW(a,dx,4,ip), ROW(a,dx,5,ip), ROW(b,dy,3,ip), ROW(b,dy,4,ip), ROW(b,dy,5,ip), v2rho2[ip], v2rhogamma[ip], v2rhotau[ip], v2gamma2[ip], v2gammatau[ip], v2tau2[ip], vgamma[ip], h );
+} else if( is_gga ) {
+  xckernel::xck_gauxc_hess_pair_gga( ROW(a,dx,0,ip), ROW(b,dy,0,ip), ROW(a,dx,1,ip), ROW(b,dy,1,ip), ROW(a,dx,2,ip), ROW(a,dx,3,ip), ROW(a,dx,4,ip), ROW(b,dy,2,ip), ROW(b,dy,3,ip), ROW(b,dy,4,ip), v2rho2[ip], v2rhogamma[ip], v2gamma2[ip], vgamma[ip], h );
+} else {
+  xckernel::xck_gauxc_hess_pair_lda( ROW(a,dx,0,ip), ROW(b,dy,0,ip), v2rho2[ip], h );
+}
+// ==> END GENERATED CODE <==

src/xc_integrator/local_work_driver/host/gauxc_hess_call_pulayW.inc

@@ -0,0 +1,11 @@
+// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: pulayW call site] <==
+// Arguments are bound BY NAME from each kernel's own signature;
+// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir <dir>
+if( is_mgga ) {
+  xckernel::xck_gauxc_hess_pulayW_mgga_tau( dden_x[ip], dden_y[ip], dden_z[ip], vrho[ip], vgamma[ip], vtau[ip], weights[ip], W[0], W[1], W[2], W[3], W[4], W[5], W[8], W[10], W[12], W[15] );
+} else if( is_gga ) {
+  xckernel::xck_gauxc_hess_pulayW_gga( dden_x[ip], dden_y[ip], dden_z[ip], vrho[ip], vgamma[ip], weights[ip], W[0], W[1], W[2], W[3], W[4], W[8], W[12] );
+} else {
+  xckernel::xck_gauxc_hess_pulayW_lda( vrho[ip], weights[ip], W[0] );
+}
+// ==> END GENERATED CODE <==

src/xc_integrator/local_work_driver/host/gauxc_hess_call_rows.inc

@@ -0,0 +1,26 @@
+// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: rows call site] <==
+// Arguments are bound BY NAME from each kernel's own signature;
+// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir <dir>
+if( is_mgga ) {
+  xckernel::xck_gauxc_hess_rows_mgga_tau( U0, U1, U2, U3, dchi, ddchi[0], ddchi[1], ddchi[2], dden_x[ip], dden_y[ip], dden_z[ip], F_rho, F_sigma, F_tau, Gx, Gy, Gz );
+  ROW(a,d,0,ip) += F_rho;
+  ROW(a,d,1,ip) += F_sigma;
+  ROW(a,d,2,ip) += F_tau;
+  ROW(a,d,3,ip) += Gx;
+  ROW(a,d,4,ip) += Gy;
+  ROW(a,d,5,ip) += Gz;
+} else if( is_gga ) {
+  xckernel::xck_gauxc_hess_rows_gga( U0, U1, U2, U3, dchi, ddchi[0], ddchi[1], ddchi[2], dden_x[ip], dden_y[ip], dden_z[ip], F_rho, F_sigma, Gx, Gy, Gz );
+  ROW(a,d,0,ip) += F_rho;
+  ROW(a,d,1,ip) += F_sigma;
+  ROW(a,d,2,ip) += Gx;
+  ROW(a,d,3,ip) += Gy;
+  ROW(a,d,4,ip) += Gz;
+} else {
+  xckernel::xck_gauxc_hess_rows_lda( U0, U1, U2, U3, dchi, ddchi[0], ddchi[1], ddchi[2], F_rho, Gx, Gy, Gz );
+  ROW(a,d,0,ip) += F_rho;
+  ROW(a,d,1,ip) += Gx;
+  ROW(a,d,2,ip) += Gy;
+  ROW(a,d,3,ip) += Gz;
+}
+// ==> END GENERATED CODE <==

src/xc_integrator/local_work_driver/host/gauxc_hess_call_same.inc

@@ -0,0 +1,11 @@
+// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: same call site] <==
+// Arguments are bound BY NAME from each kernel's own signature;
+// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir <dir>
+if( is_mgga ) {
+  xckernel::xck_gauxc_hess_same_mgga_tau( xmat[k], xmat_x[k], xmat_y[k], xmat_z[k], d2c, d3c[0], d3c[1], d3c[2], dden_x[ip], dden_y[ip], dden_z[ip], vrho[ip], vgamma[ip], vtau[ip], weights[ip], sv );
+} else if( is_gga ) {
+  xckernel::xck_gauxc_hess_same_gga( xmat[k], xmat_x[k], xmat_y[k], xmat_z[k], d2c, d3c[0], d3c[1], d3c[2], dden_x[ip], dden_y[ip], dden_z[ip], vrho[ip], vgamma[ip], weights[ip], sv );
+} else {
+  xckernel::xck_gauxc_hess_same_lda( xmat[k], d2c, vrho[ip], weights[ip], sv );
+}
+// ==> END GENERATED CODE <==