GKS

include/gauxc/xc_integrator.hpp

@@ -33,6 +33,7 @@ class XCIntegrator {
 
   using exc_vxc_type_rks  = std::tuple< value_type, matrix_type >;
   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 >;
   using exx_type      = matrix_type;
 
@@ -53,8 +54,12 @@ class XCIntegrator {
   XCIntegrator( XCIntegrator&& ) noexcept;
 
   value_type    integrate_den( const MatrixType& );
-  exc_vxc_type_rks  eval_exc_vxc ( const MatrixType& );
-  exc_vxc_type_uks  eval_exc_vxc ( const MatrixType&, const MatrixType& );
+  exc_vxc_type_rks  eval_exc_vxc ( const MatrixType&, 
+                                   const IntegratorSettingsXC& = IntegratorSettingsXC{} );
+  exc_vxc_type_uks  eval_exc_vxc ( const MatrixType&, const MatrixType&,
+                                   const IntegratorSettingsXC& = IntegratorSettingsXC{} );
+  exc_vxc_type_gks  eval_exc_vxc ( const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&,
+                                   const IntegratorSettingsXC& = IntegratorSettingsXC{});
   exc_grad_type eval_exc_grad( const MatrixType& );
   exx_type      eval_exx     ( const MatrixType&, 
                                const IntegratorSettingsEXX& = IntegratorSettingsEXX{} );

include/gauxc/xc_integrator/impl.hpp

@@ -33,18 +33,26 @@ typename XCIntegrator<MatrixType>::value_type
 
 template <typename MatrixType>
 typename XCIntegrator<MatrixType>::exc_vxc_type_rks
-  XCIntegrator<MatrixType>::eval_exc_vxc( const MatrixType& P ) {
+  XCIntegrator<MatrixType>::eval_exc_vxc( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) {
   if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED();
-  return pimpl_->eval_exc_vxc(P);
+  return pimpl_->eval_exc_vxc(P, ks_settings);
 };
 
 template <typename MatrixType>
 typename XCIntegrator<MatrixType>::exc_vxc_type_uks
-  XCIntegrator<MatrixType>::eval_exc_vxc( const MatrixType& Pscalar, const MatrixType& Pz ) {
+  XCIntegrator<MatrixType>::eval_exc_vxc( const MatrixType& Ps, const MatrixType& Pz, const IntegratorSettingsXC& ks_settings ) {
   if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED();
-  return pimpl_->eval_exc_vxc(Pscalar, Pz);
+  return pimpl_->eval_exc_vxc(Ps, Pz, ks_settings);
 };
 
+template <typename MatrixType>
+typename XCIntegrator<MatrixType>::exc_vxc_type_gks
+  XCIntegrator<MatrixType>::eval_exc_vxc( const MatrixType& Ps, const MatrixType& Pz, const MatrixType& Py, const MatrixType& Px, 
+                                          const IntegratorSettingsXC& ks_settings ) {
+      if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED();
+        return pimpl_->eval_exc_vxc(Ps, Pz, Py, Px, ks_settings);
+  };
+
 template <typename MatrixType>
 typename XCIntegrator<MatrixType>::exc_grad_type
   XCIntegrator<MatrixType>::eval_exc_grad( const MatrixType& P ) {

include/gauxc/xc_integrator/replicated/impl.hpp

@@ -63,34 +63,59 @@ typename ReplicatedXCIntegrator<MatrixType>::value_type
 
 template <typename MatrixType>
 typename ReplicatedXCIntegrator<MatrixType>::exc_vxc_type_rks 
-  ReplicatedXCIntegrator<MatrixType>::eval_exc_vxc_( const MatrixType& P ) {
+  ReplicatedXCIntegrator<MatrixType>::eval_exc_vxc_( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) {
 
   if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED();
   matrix_type VXC( P.rows(), P.cols() );
   value_type  EXC;
 
   pimpl_->eval_exc_vxc( P.rows(), P.cols(), P.data(), P.rows(),
-                        VXC.data(), VXC.rows(), &EXC );
+                        VXC.data(), VXC.rows(), &EXC, ks_settings );
 
   return std::make_tuple( EXC, VXC );
 
 }
 
 template <typename MatrixType>
 typename ReplicatedXCIntegrator<MatrixType>::exc_vxc_type_uks
-  ReplicatedXCIntegrator<MatrixType>::eval_exc_vxc_( const MatrixType& Pscalar, const MatrixType& Pz ) {
+  ReplicatedXCIntegrator<MatrixType>::eval_exc_vxc_( const MatrixType& Ps, const MatrixType& Pz, const IntegratorSettingsXC& ks_settings ) {
 
   if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED();
-  matrix_type VXCscalar( Pscalar.rows(), Pscalar.cols() );
+  matrix_type VXCs( Ps.rows(), Ps.cols() );
   matrix_type VXCz( Pz.rows(), Pz.cols() );
   value_type  EXC;
 
-  pimpl_->eval_exc_vxc( Pscalar.rows(), Pscalar.cols(), Pscalar.data(), Pscalar.rows(),
+  pimpl_->eval_exc_vxc( Ps.rows(), Ps.cols(), Ps.data(), Ps.rows(),
                         Pz.data(), Pz.rows(),
-                        VXCscalar.data(), VXCscalar.rows(),
-                        VXCz.data(), VXCz.rows(), &EXC );
+                        VXCs.data(), VXCs.rows(),
+                        VXCz.data(), VXCz.rows(), &EXC, ks_settings );
 
-  return std::make_tuple( EXC, VXCscalar, VXCz );
+  return std::make_tuple( EXC, VXCs, VXCz );
+
+}
+
+template <typename MatrixType>
+typename ReplicatedXCIntegrator<MatrixType>::exc_vxc_type_gks
+  ReplicatedXCIntegrator<MatrixType>::eval_exc_vxc_( const MatrixType& Ps, const MatrixType& Pz, const MatrixType& Py, const MatrixType& Px,
+                                                     const IntegratorSettingsXC& ks_settings) {
+
+  if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED();
+  matrix_type VXCs( Ps.rows(), Ps.cols() );
+  matrix_type VXCz( Pz.rows(), Pz.cols() );
+  matrix_type VXCy( Py.rows(), Py.cols() );
+  matrix_type VXCx( Px.rows(), Px.cols() );
+  value_type  EXC;
+
+  pimpl_->eval_exc_vxc( Ps.rows(), Ps.cols(), Ps.data(), Ps.rows(),
+                        Pz.data(), Pz.rows(),
+                        Py.data(), Py.rows(),
+                        Px.data(), Px.rows(),
+                        VXCs.data(), VXCs.rows(),
+                        VXCz.data(), VXCz.rows(),
+                        VXCy.data(), VXCy.rows(),
+                        VXCx.data(), VXCx.rows(), &EXC, ks_settings );
+
+  return std::make_tuple( EXC, VXCs, VXCz, VXCy, VXCx);
 
 }
 

include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp

@@ -40,14 +40,28 @@ class ReplicatedXCIntegratorImpl {
                                int64_t ldp, value_type* N_EL ) = 0;
   virtual void eval_exc_vxc_( int64_t m, int64_t n, const value_type* P,
                               int64_t ldp, value_type* VXC, int64_t ldvxc,
-                              value_type* EXC ) = 0;
-  virtual void eval_exc_vxc_( int64_t m, int64_t n, const value_type* Pscalar,
-                              int64_t ldpscalar,
+                              value_type* EXC, const IntegratorSettingsXC& ks_settings ) = 0;
+  virtual void eval_exc_vxc_( int64_t m, int64_t n, const value_type* Ps,
+                              int64_t ldps,
                               const value_type* Pz,
                               int64_t ldpz,
-                              value_type* VXCscalar, int64_t ldvxcscalar,
+                              value_type* VXCs, int64_t ldvxcs,
                               value_type* VXCz, int64_t ldvxcz,
-                              value_type* EXC ) = 0;
+                              value_type* EXC, const IntegratorSettingsXC& ks_settings ) = 0;
+  virtual void eval_exc_vxc_( int64_t m, int64_t n, const value_type* Ps,
+                              int64_t ldps,
+                              const value_type* Pz,
+                              int64_t ldpz,
+                              const value_type* Py,
+                              int64_t ldpy,
+                              const value_type* Px,
+                              int64_t ldpx,
+                              value_type* VXCs, int64_t ldvxcs,
+                              value_type* VXCz, int64_t ldvxcz,
+                              value_type* VXCy, int64_t ldvxcy,
+                              value_type* VXCx, int64_t ldvxcx,
+                              value_type* EXC, const IntegratorSettingsXC& ks_settings ) = 0;
+
   virtual void eval_exc_grad_( int64_t m, int64_t n, const value_type* P,
                                int64_t ldp, value_type* EXC_GRAD ) = 0;
   virtual void eval_exx_( int64_t m, int64_t n, const value_type* P,
@@ -69,15 +83,29 @@ class ReplicatedXCIntegratorImpl {
 
   void eval_exc_vxc( int64_t m, int64_t n, const value_type* P,
                      int64_t ldp, value_type* VXC, int64_t ldvxc,
-                     value_type* EXC ); 
+                     value_type* EXC, const IntegratorSettingsXC& ks_settings ); 
 
-  void eval_exc_vxc( int64_t m, int64_t n, const value_type* Pscalar,
-                     int64_t ldpscalar,
+  void eval_exc_vxc( int64_t m, int64_t n, const value_type* Ps,
+                     int64_t ldps,
                      const value_type* Pz,
                      int64_t ldpz,
-                     value_type* VXCscalar, int64_t ldvxcscalar,
+                     value_type* VXCs, int64_t ldvxcs,
                      value_type* VXCz, int64_t ldvxcz,
-                     value_type* EXC );
+                     value_type* EXC, const IntegratorSettingsXC& ks_settings );
+  void eval_exc_vxc( int64_t m, int64_t n, const value_type* Ps,
+                     int64_t ldps,
+                     const value_type* Pz,
+                     int64_t ldpz,
+                     const value_type* Py,
+                     int64_t ldpy,
+                     const value_type* Px,
+                     int64_t ldpx,
+                     value_type* VXCs, int64_t ldvxcs,
+                     value_type* VXCz, int64_t ldvxcz,
+                     value_type* VXCy, int64_t ldvxcy,
+                     value_type* VXCx, int64_t ldvxcx,
+                     value_type* EXC, const IntegratorSettingsXC& ks_settings );
+
 
   void eval_exc_grad( int64_t m, int64_t n, const value_type* P,
                       int64_t ldp, value_type* EXC_GRAD );

include/gauxc/xc_integrator/replicated_xc_integrator.hpp

@@ -30,6 +30,7 @@ class ReplicatedXCIntegrator : public XCIntegratorImpl<MatrixType> {
   using value_type     = typename XCIntegratorImpl<MatrixType>::value_type;
   using exc_vxc_type_rks   = typename XCIntegratorImpl<MatrixType>::exc_vxc_type_rks;
   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 exx_type       = typename XCIntegratorImpl<MatrixType>::exx_type;
 
@@ -39,8 +40,9 @@ class ReplicatedXCIntegrator : public XCIntegratorImpl<MatrixType> {
   std::unique_ptr< pimpl_type > pimpl_;
 
   value_type    integrate_den_( const MatrixType& ) override;
-  exc_vxc_type_rks  eval_exc_vxc_ ( const MatrixType& ) override;
-  exc_vxc_type_uks  eval_exc_vxc_ ( const MatrixType&, const MatrixType& ) override;
+  exc_vxc_type_rks  eval_exc_vxc_ ( const MatrixType&, const IntegratorSettingsXC& ) override;
+  exc_vxc_type_uks  eval_exc_vxc_ ( const MatrixType&, const MatrixType&, const IntegratorSettingsXC&) override;
+  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& ) override;
   exx_type      eval_exx_     ( const MatrixType&, const IntegratorSettingsEXX& ) override;
   const util::Timer& get_timings_() const override;

include/gauxc/xc_integrator/xc_integrator_impl.hpp

@@ -22,14 +22,17 @@ class XCIntegratorImpl {
   using value_type     = typename matrix_type::value_type;
   using exc_vxc_type_rks   = typename XCIntegrator<MatrixType>::exc_vxc_type_rks;
   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 exx_type       = typename XCIntegrator<MatrixType>::exx_type;
 
 protected:
 
   virtual value_type    integrate_den_( const MatrixType& P ) = 0;
-  virtual exc_vxc_type_rks  eval_exc_vxc_ ( const MatrixType& P ) = 0;
-  virtual exc_vxc_type_uks  eval_exc_vxc_ ( const MatrixType& Pscalar, const MatrixType& Pz ) = 0;
+  virtual exc_vxc_type_rks  eval_exc_vxc_ ( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) = 0;
+  virtual exc_vxc_type_uks  eval_exc_vxc_ ( const MatrixType& Ps, const MatrixType& Pz, const IntegratorSettingsXC& ks_settings ) = 0;
+  virtual exc_vxc_type_gks  eval_exc_vxc_ ( const MatrixType& Ps, const MatrixType& Pz, const MatrixType& Py, const MatrixType& Px, 
+                                            const IntegratorSettingsXC& ks_settings ) = 0;
   virtual exc_grad_type eval_exc_grad_( const MatrixType& P ) = 0;
   virtual exx_type      eval_exx_     ( const MatrixType&     P, 
                                         const IntegratorSettingsEXX& settings ) = 0;
@@ -62,12 +65,16 @@ class XCIntegratorImpl {
    *  @param[in] P The alpha density matrix
    *  @returns EXC / VXC in a combined structure
    */
-  exc_vxc_type_rks eval_exc_vxc( const MatrixType& P ) {
-    return eval_exc_vxc_(P);
+  exc_vxc_type_rks eval_exc_vxc( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) {
+    return eval_exc_vxc_(P, ks_settings);
   }
 
-  exc_vxc_type_uks eval_exc_vxc( const MatrixType& Pscalar, const MatrixType& Pz ) {
-    return eval_exc_vxc_(Pscalar, Pz);
+  exc_vxc_type_uks eval_exc_vxc( const MatrixType& Ps, const MatrixType& Pz, const IntegratorSettingsXC& ks_settings ) {
+    return eval_exc_vxc_(Ps, Pz, ks_settings);
+  }
+
+  exc_vxc_type_gks eval_exc_vxc( const MatrixType& Ps, const MatrixType& Pz, const MatrixType& Py, const MatrixType& Px, const IntegratorSettingsXC& ks_settings ) {
+    return eval_exc_vxc_(Ps, Pz, Py, Px, ks_settings);
   }
 
   /** Integrate EXC gradient for RKS

include/gauxc/xc_integrator_settings.hpp

@@ -16,4 +16,9 @@ struct IntegratorSettingsSNLinK : public IntegratorSettingsEXX {
   double k_tol      = 1e-10;
 };
 
+struct IntegratorSettingsXC { virtual ~IntegratorSettingsXC() noexcept = default; };
+struct IntegratorSettingsKS : public IntegratorSettingsXC {
+  double gks_dtol = 1e-12;
+};
+
 }

src/xc_integrator/local_work_driver/device/local_device_work_driver.cxx

@@ -44,12 +44,18 @@ FWD_TO_PIMPL(eval_uvvar_gga_rks)            // U/VVar GGA (density + grad, gamma
 FWD_TO_PIMPL(eval_uvvar_lda_uks)            // U/VVar LDA (density)
 FWD_TO_PIMPL(eval_uvvar_gga_uks)            // U/VVar GGA (density + grad, gamma)
 
+FWD_TO_PIMPL(eval_uvvar_lda_gks)            // U/VVar LDA (density)
+FWD_TO_PIMPL(eval_uvvar_gga_gks)            // U/VVar GGA (density + grad, gamma)
+
 FWD_TO_PIMPL(eval_zmat_lda_vxc_rks)         // Eval Z Matrix LDA VXC
 FWD_TO_PIMPL(eval_zmat_gga_vxc_rks)         // Eval Z Matrix GGA VXC
 
 FWD_TO_PIMPL(eval_zmat_lda_vxc_uks)         // Eval Z Matrix LDA VXC
 FWD_TO_PIMPL(eval_zmat_gga_vxc_uks)         // Eval Z Matrix GGA VXC
 
+FWD_TO_PIMPL(eval_zmat_lda_vxc_gks)         // Eval Z Matrix LDA VXC
+FWD_TO_PIMPL(eval_zmat_gga_vxc_gks)         // Eval Z Matrix GGA VXC
+
 FWD_TO_PIMPL(eval_exx_fmat)             // Eval EXX F Matrix
 //FWD_TO_PIMPL(eval_exx_gmat)             // Eval EXX G Matrix
 

src/xc_integrator/local_work_driver/device/local_device_work_driver.hpp

@@ -69,6 +69,9 @@ class LocalDeviceWorkDriver : public LocalWorkDriver {
   void eval_uvvar_lda_uks( XCDeviceData* );
   void eval_uvvar_gga_uks( XCDeviceData* );
 
+  void eval_uvvar_lda_gks( XCDeviceData* );
+  void eval_uvvar_gga_gks( XCDeviceData* );
+
   void eval_kern_exc_vxc_lda( const functional_type&, XCDeviceData* );
   void eval_kern_exc_vxc_gga( const functional_type&, XCDeviceData* );
 
@@ -78,6 +81,9 @@ class LocalDeviceWorkDriver : public LocalWorkDriver {
   void eval_zmat_lda_vxc_uks( XCDeviceData* );
   void eval_zmat_gga_vxc_uks( XCDeviceData* );
 
+  void eval_zmat_lda_vxc_gks( XCDeviceData* );
+  void eval_zmat_gga_vxc_gks( XCDeviceData* );
+
   void eval_exx_fmat( XCDeviceData* );
   void eval_exx_gmat( XCDeviceData*, const BasisSetMap& );
 

src/xc_integrator/local_work_driver/device/local_device_work_driver_pimpl.hpp

@@ -41,12 +41,16 @@ struct LocalDeviceWorkDriverPIMPL {
   virtual void eval_uvvar_gga_rks( XCDeviceData* ) = 0;
   virtual void eval_uvvar_lda_uks( XCDeviceData* ) = 0;
   virtual void eval_uvvar_gga_uks( XCDeviceData* ) = 0;
+  virtual void eval_uvvar_lda_gks( XCDeviceData* ) = 0;  
+  virtual void eval_uvvar_gga_gks( XCDeviceData* ) = 0;
   virtual void eval_kern_exc_vxc_lda( const functional_type&, XCDeviceData* ) = 0;
   virtual void eval_kern_exc_vxc_gga( const functional_type&, XCDeviceData* ) = 0;
   virtual void eval_zmat_lda_vxc_rks( XCDeviceData* ) = 0;
   virtual void eval_zmat_gga_vxc_rks( XCDeviceData* ) = 0;
   virtual void eval_zmat_lda_vxc_uks( XCDeviceData* ) = 0;
   virtual void eval_zmat_gga_vxc_uks( XCDeviceData* ) = 0;
+  virtual void eval_zmat_lda_vxc_gks( XCDeviceData* ) = 0;
+  virtual void eval_zmat_gga_vxc_gks( XCDeviceData* ) = 0;
   virtual void inc_exc( XCDeviceData* ) = 0;
   virtual void inc_nel( XCDeviceData* ) = 0;
   virtual void inc_vxc( XCDeviceData* ) = 0;

src/xc_integrator/local_work_driver/device/scheme1_base.cxx

@@ -288,6 +288,18 @@ void AoSScheme1Base::eval_zmat_gga_vxc_uks( XCDeviceData* ){
 
 }
 
+void AoSScheme1Base::eval_zmat_lda_vxc_gks( XCDeviceData* ){
+
+  GAUXC_GENERIC_EXCEPTION("GKS NOT YET IMPLEMENTED FOR DEVICE");
+
+}
+
+void AoSScheme1Base::eval_zmat_gga_vxc_gks( XCDeviceData* ){
+
+  GAUXC_GENERIC_EXCEPTION("GKS NOT YET IMPLEMENTED FOR DEVICE");
+
+}
+
 void AoSScheme1Base::eval_collocation( XCDeviceData* _data ) {
 
   auto* data = dynamic_cast<Data*>(_data);
@@ -489,6 +501,18 @@ void AoSScheme1Base::eval_uvvar_gga_uks( XCDeviceData* ){
 }
 
 
+void AoSScheme1Base::eval_uvvar_lda_gks( XCDeviceData* ){
+
+    GAUXC_GENERIC_EXCEPTION("GKS NOT YET IMPLEMENTED FOR DEVICE");
+
+}
+
+void AoSScheme1Base::eval_uvvar_gga_gks( XCDeviceData* ){
+
+    GAUXC_GENERIC_EXCEPTION("GKS NOT YET IMPLEMENTED FOR DEVICE");
+
+}
+
 void AoSScheme1Base::eval_kern_exc_vxc_lda( const functional_type& func, 
   XCDeviceData* _data ) {
 

src/xc_integrator/local_work_driver/device/scheme1_base.hpp

@@ -27,6 +27,11 @@ struct AoSScheme1Base : public detail::LocalDeviceWorkDriverPIMPL {
   void eval_zmat_lda_vxc_uks( XCDeviceData* ) override final;
   void eval_zmat_gga_vxc_uks( XCDeviceData* ) override final;
 
+  void eval_uvvar_lda_gks( XCDeviceData* ) override final;
+  void eval_uvvar_gga_gks( XCDeviceData* ) override final;
+  void eval_zmat_lda_vxc_gks( XCDeviceData* ) override final;
+  void eval_zmat_gga_vxc_gks( XCDeviceData* ) override final;
+
   void eval_kern_exc_vxc_lda( const functional_type&, XCDeviceData* ) override final;
   void eval_kern_exc_vxc_gga( const functional_type&, XCDeviceData* ) override final;
   void inc_exc( XCDeviceData* ) override final;