From 61152affa2a196e65cced221a16a066aa2553a29 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 19 Sep 2026 14:28:17 +0300 Subject: [PATCH] fxc_contraction: GKS (noncollinear) contraction for LDA eval_fxc_contraction( Ps, Pz, Py, Px, tPs, tPz, tPy, tPx ) returns ( FXCs, FXCz, FXCy, FXCx ), in the argument order of the GKS eval_exc_vxc, on the host replicated integrator; other integrators throw (the new virtuals have defaults, so no device stubs). The energy is the locally collinear one the GKS potential uses, n_+- = (rho_s +- |m|)/2. The kernel applied to a trial density, per noncollinear field slot, is generated by xckernel's ncwriter (the mechanical second derivative of that map) and assembled as the potential is, so FXC_X = d/dh VXC_X(P + h tP). Below gks_dtol the generated collinear limit is used: each magnetization component responds like the spin channel of a collinear perturbation about the spin-symmetric reference. Trial densities are symmetrized as in #225. H3/cc-pVDZ, LDA_X+LDA_C_PW, UltraFine, against central differences of the GKS potential (all four blocks): noncollinear reference s 1.6e-7 -> 3.9e-8, x/y/z 1e-7..3e-7 -> 2e-8..8e-8 m = 0 (limit branch) s 1.3e-7 -> 3.1e-8, x/y/z likewise as the step halves (the h^2 truncation of the reference), and sum_X Tr(X F_X[Y]) = sum_X Tr(Y F_X[X]) to 1e-14. GGA is not enabled: the Scalmani-Frisch gamma map's kernel carries 1/|(grad rho_s . grad m_J)_J|, singular at density critical points at finite |m|, and needs a regularization of its own. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_0135evJ9zgNL1y8U6T9bQ3UT --- include/gauxc/xc_integrator.hpp | 13 + include/gauxc/xc_integrator/impl.hpp | 9 + .../gauxc/xc_integrator/replicated/impl.hpp | 20 + .../replicated_xc_integrator_impl.hpp | 20 + .../replicated_xc_integrator.hpp | 3 + .../xc_integrator/xc_integrator_impl.hpp | 16 + .../host/gauxc_nc_kernel.hpp | 1283 +++++++++++++++++ ...eference_replicated_xc_host_integrator.cxx | 1 + ...eference_replicated_xc_host_integrator.hpp | 18 + ...xc_host_integrator_fxc_contraction_gks.hpp | 220 +++ .../replicated_xc_integrator_impl.cxx | 18 + 11 files changed, 1621 insertions(+) create mode 100644 src/xc_integrator/local_work_driver/host/gauxc_nc_kernel.hpp create mode 100644 src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction_gks.hpp diff --git a/include/gauxc/xc_integrator.hpp b/include/gauxc/xc_integrator.hpp index 03feaf934..88fa8a331 100644 --- a/include/gauxc/xc_integrator.hpp +++ b/include/gauxc/xc_integrator.hpp @@ -42,6 +42,7 @@ class XCIntegrator { using exx_type = matrix_type; using fxc_contraction_type_rks = matrix_type; using fxc_contraction_type_uks = std::tuple< matrix_type, matrix_type >; + using fxc_contraction_type_gks = std::tuple< matrix_type, matrix_type, matrix_type, matrix_type >; using dd_psi_type = std::vector< value_type >; using dd_psi_potential_type = matrix_type; @@ -80,10 +81,22 @@ 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 ). + * GKS: eval_fxc_contraction( Ps, Pz, Py, Px, tPs, tPz, tPy, tPx ), returning + * ( FXCs, FXCz, FXCy, FXCx ), in the argument order of the GKS + * eval_exc_vxc. LDA only. + */ 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&, const IntegratorSettingsXC& = IntegratorSettingsXC{} ); + fxc_contraction_type_gks eval_fxc_contraction ( const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, + const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, + const IntegratorSettingsXC& = IntegratorSettingsXC{} ); dd_psi_type eval_dd_psi( const MatrixType&, unsigned ); dd_psi_potential_type eval_dd_psi_potential( const MatrixType&, unsigned ); diff --git a/include/gauxc/xc_integrator/impl.hpp b/include/gauxc/xc_integrator/impl.hpp index 400afb7c7..0e97705a7 100644 --- a/include/gauxc/xc_integrator/impl.hpp +++ b/include/gauxc/xc_integrator/impl.hpp @@ -116,6 +116,15 @@ typename XCIntegrator::fxc_contraction_type_uks return pimpl_->eval_fxc_contraction(Ps, Pz, tPs, tPz, ks_settings); }; +template +typename XCIntegrator::fxc_contraction_type_gks + XCIntegrator::eval_fxc_contraction( const MatrixType& Ps, const MatrixType& Pz, const MatrixType& Py, const MatrixType& Px, + const MatrixType& tPs, const MatrixType& tPz, const MatrixType& tPy, const MatrixType& tPx, + const IntegratorSettingsXC& ks_settings ) { + if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED(); + return pimpl_->eval_fxc_contraction(Ps, Pz, Py, Px, tPs, tPz, tPy, tPx, ks_settings); +}; + template typename XCIntegrator::dd_psi_type XCIntegrator::eval_dd_psi(const MatrixType& P, unsigned max_Ylm) { diff --git a/include/gauxc/xc_integrator/replicated/impl.hpp b/include/gauxc/xc_integrator/replicated/impl.hpp index bfc95fc88..ae23b66b9 100644 --- a/include/gauxc/xc_integrator/replicated/impl.hpp +++ b/include/gauxc/xc_integrator/replicated/impl.hpp @@ -238,6 +238,26 @@ typename ReplicatedXCIntegrator::fxc_contraction_type_uks } +template +typename ReplicatedXCIntegrator::fxc_contraction_type_gks + ReplicatedXCIntegrator::eval_fxc_contraction_( const MatrixType& Ps, const MatrixType& Pz, const MatrixType& Py, const MatrixType& Px, + const MatrixType& tPs, const MatrixType& tPz, const MatrixType& tPy, const MatrixType& tPx, + const IntegratorSettingsXC& ks_settings ) { + + if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED(); + matrix_type FXCs( Ps.rows(), Ps.cols() ), FXCz( Pz.rows(), Pz.cols() ); + matrix_type FXCy( Py.rows(), Py.cols() ), FXCx( Px.rows(), Px.cols() ); + + pimpl_->eval_fxc_contraction( Ps.rows(), Ps.cols(), + Ps.data(), Ps.rows(), Pz.data(), Pz.rows(), Py.data(), Py.rows(), Px.data(), Px.rows(), + tPs.data(), tPs.rows(), tPz.data(), tPz.rows(), tPy.data(), tPy.rows(), tPx.data(), tPx.rows(), + FXCs.data(), FXCs.rows(), FXCz.data(), FXCz.rows(), FXCy.data(), FXCy.rows(), FXCx.data(), FXCx.rows(), + ks_settings ); + + return std::make_tuple( FXCs, FXCz, FXCy, FXCx ); + +} + template typename ReplicatedXCIntegrator::dd_psi_type ReplicatedXCIntegrator::eval_dd_psi_( const MatrixType& P, unsigned max_Ylm ) { diff --git a/include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp b/include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp index 457315122..783e7e529 100644 --- a/include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp +++ b/include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp @@ -98,6 +98,17 @@ class ReplicatedXCIntegratorImpl { value_type* FXCs, int64_t ldfxcs, value_type* FXCz, int64_t ldfxcz, const IntegratorSettingsXC& ks_settings )=0; + // GKS: not pure, so the device integrators need no stub + virtual void eval_fxc_contraction_( 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, + const value_type* tPs, int64_t ldtps, const value_type* tPz, int64_t ldtpz, + const value_type* tPy, int64_t ldtpy, const value_type* tPx, int64_t ldtpx, + value_type* FXCs, int64_t ldfxcs, value_type* FXCz, int64_t ldfxcz, + value_type* FXCy, int64_t ldfxcy, value_type* FXCx, int64_t ldfxcx, + const IntegratorSettingsXC& ) { + GAUXC_GENERIC_EXCEPTION("GKS FXC Contraction Not Implemented For This Integrator"); + } virtual void eval_dd_psi_( int64_t m, int64_t n, const value_type* P, int64_t ldp, unsigned max_Ylm, value_type* ddPsi, int64_t ldPsi ) = 0; virtual void eval_dd_psi_potential_( int64_t m, int64_t n, const value_type* X, unsigned max_Ylm, @@ -177,6 +188,15 @@ class ReplicatedXCIntegratorImpl { value_type* FXCz, int64_t ldfxcz, const IntegratorSettingsXC& ks_settings ); + void eval_fxc_contraction( 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, + const value_type* tPs, int64_t ldtps, const value_type* tPz, int64_t ldtpz, + const value_type* tPy, int64_t ldtpy, const value_type* tPx, int64_t ldtpx, + value_type* FXCs, int64_t ldfxcs, value_type* FXCz, int64_t ldfxcz, + value_type* FXCy, int64_t ldfxcy, value_type* FXCx, int64_t ldfxcx, + const IntegratorSettingsXC& ks_settings ); + void eval_dd_psi( int64_t m, int64_t n, const value_type* P, int64_t ldp, unsigned max_Ylm, value_type* ddPsi, int64_t ldPsi ); diff --git a/include/gauxc/xc_integrator/replicated_xc_integrator.hpp b/include/gauxc/xc_integrator/replicated_xc_integrator.hpp index 1ca53f917..22345dc51 100644 --- a/include/gauxc/xc_integrator/replicated_xc_integrator.hpp +++ b/include/gauxc/xc_integrator/replicated_xc_integrator.hpp @@ -39,6 +39,7 @@ class ReplicatedXCIntegrator : public XCIntegratorImpl { using exx_type = typename XCIntegratorImpl::exx_type; using fxc_contraction_type_rks = typename XCIntegratorImpl::fxc_contraction_type_rks; using fxc_contraction_type_uks = typename XCIntegratorImpl::fxc_contraction_type_uks; + using fxc_contraction_type_gks = typename XCIntegratorImpl::fxc_contraction_type_gks; using dd_psi_type = typename XCIntegratorImpl::dd_psi_type; using dd_psi_potential_type = typename XCIntegratorImpl::dd_psi_potential_type; @@ -59,6 +60,8 @@ class ReplicatedXCIntegrator : public XCIntegratorImpl { 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; + fxc_contraction_type_gks eval_fxc_contraction_ ( const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, + const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, const IntegratorSettingsXC&) override; dd_psi_type eval_dd_psi_( const MatrixType& , unsigned ) override; dd_psi_potential_type eval_dd_psi_potential_( const MatrixType& , unsigned ) override; const util::Timer& get_timings_() const override; diff --git a/include/gauxc/xc_integrator/xc_integrator_impl.hpp b/include/gauxc/xc_integrator/xc_integrator_impl.hpp index ba7bebebb..1843e7799 100644 --- a/include/gauxc/xc_integrator/xc_integrator_impl.hpp +++ b/include/gauxc/xc_integrator/xc_integrator_impl.hpp @@ -31,6 +31,7 @@ class XCIntegratorImpl { using exx_type = typename XCIntegrator::exx_type; using fxc_contraction_type_rks = typename XCIntegrator::fxc_contraction_type_rks; using fxc_contraction_type_uks = typename XCIntegrator::fxc_contraction_type_uks; + using fxc_contraction_type_gks = typename XCIntegrator::fxc_contraction_type_gks; using dd_psi_type = typename XCIntegrator::dd_psi_type; using dd_psi_potential_type = typename XCIntegrator::dd_psi_potential_type; @@ -54,6 +55,12 @@ class XCIntegratorImpl { const MatrixType& tP, const IntegratorSettingsXC& ks_settings ) = 0; virtual fxc_contraction_type_uks eval_fxc_contraction_ ( const MatrixType& Ps, const MatrixType& Pz, const MatrixType& tPs, const MatrixType& tPz, const IntegratorSettingsXC& ks_settings ) = 0; + // GKS: not pure, so integrators without it (the device ones) need no stub + virtual fxc_contraction_type_gks eval_fxc_contraction_ ( const MatrixType& Ps, const MatrixType& Pz, const MatrixType& Py, const MatrixType& Px, + const MatrixType& tPs, const MatrixType& tPz, const MatrixType& tPy, const MatrixType& tPx, + const IntegratorSettingsXC& ) { + GAUXC_GENERIC_EXCEPTION("GKS FXC Contraction Not Implemented For This Integrator"); + } virtual dd_psi_type eval_dd_psi_( const MatrixType& P, unsigned max_Ylm ) = 0; @@ -175,6 +182,15 @@ class XCIntegratorImpl { return eval_fxc_contraction_(Ps, Pz, tPs, tPz, ks_settings); } + /** Evaluate the GKS FXC contraction (LDA only), returning + * ( FXCs, FXCz, FXCy, FXCx ). + */ + fxc_contraction_type_gks eval_fxc_contraction( const MatrixType& Ps, const MatrixType& Pz, const MatrixType& Py, const MatrixType& Px, + const MatrixType& tPs, const MatrixType& tPz, const MatrixType& tPy, const MatrixType& tPx, + const IntegratorSettingsXC& ks_settings ) { + return eval_fxc_contraction_(Ps, Pz, Py, Px, tPs, tPz, tPy, tPx, ks_settings); + } + /** Evaluate Psi vector for ddX * * @param[in] P The density matrix diff --git a/src/xc_integrator/local_work_driver/host/gauxc_nc_kernel.hpp b/src/xc_integrator/local_work_driver/host/gauxc_nc_kernel.hpp new file mode 100644 index 000000000..e612a49bc --- /dev/null +++ b/src/xc_integrator/local_work_driver/host/gauxc_nc_kernel.hpp @@ -0,0 +1,1283 @@ +/** + * Noncollinear (relativistic) exchange-correlation potential and response + * kernels, for the locally collinear ansatz. + * + * MACHINE-GENERATED by xckernel (ncwriter); do not edit. + * Copyright (c) 2026 Susi Lehtola. + * + * The functions are pure, per-grid-point and allocation-free: the caller + * owns the grid loop. Conventions, per grid point: + * + * rho_s charge density + * rho_{x,y,z} magnetization vector m + * grad_rho_s_{x,y,z}, grad_rho_{x,y,z}_{x,y,z} their Cartesian gradients + * f_nabla sgn( grad rho_s . sum_J rho_J grad rho_J ), +-1 + * + * The collinear variables handed to the functional library are + * + * n_{+,-} = ( rho_s +- |m| ) / 2 + * gamma^{++,--} = ( |grad rho_s|^2 + sum_J |grad rho_J|^2 ) / 4 + * +- f_nabla sqrt( sum_J (grad rho_s . grad rho_J)^2 ) / 2 + * gamma^{+-} = ( |grad rho_s|^2 - sum_J |grad rho_J|^2 ) / 4 + * + * and the derivative arguments (vrho_0, ..., v2sigma2_5) follow the + * standard spin-polarized Libxc packing, so they can be passed straight + * through from a polarized functional evaluation. + * + * NOTE the |m| -> 0 limit: the map is only piecewise smooth there and the + * transverse kernel components carry 1/|m|. Hosts must apply their usual + * small-|m| regularization (a magnetization cutoff) before calling. + */ + +#pragma once + +#include + +namespace xckernel { + +/** Noncollinear LDA exchange-correlation potential: the derivative of the energy density with respect to each noncollinear field, v_X = dE/dX. */ +inline void nc_vxc_lda( [[maybe_unused]] double rho_s, double rho_x, double rho_y, double rho_z, double vrho_0, double vrho_1, [[maybe_unused]] double f_nabla, double& v_rho_s, double& v_rho_x, double& v_rho_y, double& v_rho_z ) { + const double t0 = (1.0/2.0)*(vrho_0 - vrho_1)/sqrt(pow(rho_x, 2) + pow(rho_y, 2) + pow(rho_z, 2)); + v_rho_s = (1.0/2.0)*(vrho_0 + vrho_1); + v_rho_x = rho_x*t0; + v_rho_y = rho_y*t0; + v_rho_z = rho_z*t0; +} + +/** Noncollinear LDA exchange-correlation kernel applied to one trial density: k_X = sum_Y f_xc[X,Y] tY, delivered in the same field slots as the potential. */ +inline void nc_fxc_contract_lda( [[maybe_unused]] double rho_s, double rho_x, double rho_y, double rho_z, double vrho_0, double vrho_1, double v2rho2_0, double v2rho2_1, double v2rho2_2, [[maybe_unused]] double f_nabla, double trho_s, double trho_x, double trho_y, double trho_z, double& k_rho_s, double& k_rho_x, double& k_rho_y, double& k_rho_z ) { + const double t0 = 2*v2rho2_1; + const double t1 = pow(rho_x, 2); + const double t2 = pow(rho_y, 2); + const double t3 = pow(rho_z, 2); + const double t4 = t1 + t2 + t3; + const double t5 = pow(t4, -1.0/2.0); + const double t6 = t5*(v2rho2_0 - v2rho2_2); + const double t7 = rho_x*trho_x; + const double t8 = rho_y*trho_y; + const double t9 = rho_z*trho_z; + const double t10 = t6*trho_s; + const double t11 = 1.0/t4; + const double t12 = t11*v2rho2_0; + const double t13 = t11*v2rho2_2; + const double t14 = 2/pow(t4, 3.0/2.0); + const double t15 = -t0*t11 + t12 + t13 - t14*vrho_0 + t14*vrho_1; + const double t16 = rho_x*t15; + const double t17 = t1*t11; + const double t18 = 2*t5; + const double t19 = t18*(t17 - 1); + const double t20 = rho_y*t15; + const double t21 = t11*t2; + const double t22 = t18*(t21 - 1); + const double t23 = rho_z*t15; + const double t24 = t11*t3; + const double t25 = t18*(t24 - 1); + k_rho_s = (1.0/4.0)*t6*t7 + (1.0/4.0)*t6*t8 + (1.0/4.0)*t6*t9 + (1.0/4.0)*trho_s*(t0 + v2rho2_0 + v2rho2_2); + k_rho_x = (1.0/4.0)*rho_x*t10 + (1.0/4.0)*t16*t8 + (1.0/4.0)*t16*t9 + (1.0/4.0)*trho_x*(-t0*t17 + t17*v2rho2_0 + t17*v2rho2_2 - t19*vrho_0 + t19*vrho_1); + k_rho_y = (1.0/4.0)*rho_y*t10 + (1.0/4.0)*t20*t7 + (1.0/4.0)*t20*t9 + (1.0/4.0)*trho_y*(-t0*t21 + t12*t2 + t13*t2 - t22*vrho_0 + t22*vrho_1); + k_rho_z = (1.0/4.0)*rho_z*t10 + (1.0/4.0)*t23*t7 + (1.0/4.0)*t23*t8 + (1.0/4.0)*trho_z*(-t0*t24 + t12*t3 + t13*t3 - t25*vrho_0 + t25*vrho_1); +} + +/** Collinear (|m| -> 0) limit of nc_fxc_contract_lda, for points below the host's magnetization cutoff: each magnetization component responds like the spin channel of a collinear perturbation about the spin-symmetric reference; no charge-spin coupling. Depends on the charge fields and the polarized derivatives only. */ +inline void nc_fxc_contract_limit_lda( [[maybe_unused]] double rho_s, [[maybe_unused]] double vrho_0, [[maybe_unused]] double vrho_1, double v2rho2_0, double v2rho2_1, double v2rho2_2, double trho_s, double trho_x, double trho_y, double trho_z, double& k_rho_s, double& k_rho_x, double& k_rho_y, double& k_rho_z ) { + const double t0 = 2*v2rho2_1; + const double t1 = v2rho2_0 + v2rho2_2; + const double t2 = -1.0/4.0*t0 + (1.0/4.0)*t1; + k_rho_s = (1.0/4.0)*trho_s*(t0 + t1); + k_rho_x = t2*trho_x; + k_rho_y = t2*trho_y; + k_rho_z = t2*trho_z; +} + +/** Noncollinear LDA exchange-correlation kernel matrix f_xc[X,Y] over the noncollinear field slots; upper triangle, symmetric. */ +inline void nc_fxc_matrix_lda( [[maybe_unused]] double rho_s, double rho_x, double rho_y, double rho_z, double vrho_0, double vrho_1, double v2rho2_0, double v2rho2_1, double v2rho2_2, [[maybe_unused]] double f_nabla, double& c_rho_s_rho_s, double& c_rho_s_rho_x, double& c_rho_s_rho_y, double& c_rho_s_rho_z, double& c_rho_x_rho_x, double& c_rho_x_rho_y, double& c_rho_x_rho_z, double& c_rho_y_rho_y, double& c_rho_y_rho_z, double& c_rho_z_rho_z ) { + const double t0 = 2*v2rho2_1; + const double t1 = (1.0/4.0)*rho_x; + const double t2 = pow(rho_x, 2); + const double t3 = pow(rho_y, 2); + const double t4 = pow(rho_z, 2); + const double t5 = t2 + t3 + t4; + const double t6 = pow(t5, -1.0/2.0); + const double t7 = t6*(v2rho2_0 - v2rho2_2); + const double t8 = (1.0/4.0)*t7; + const double t9 = 1.0/t5; + const double t10 = t2*t9; + const double t11 = 2*t6; + const double t12 = t11*(t10 - 1); + const double t13 = t9*v2rho2_0; + const double t14 = t9*v2rho2_2; + const double t15 = 2/pow(t5, 3.0/2.0); + const double t16 = -t0*t9 + t13 + t14 - t15*vrho_0 + t15*vrho_1; + const double t17 = t1*t16; + const double t18 = t3*t9; + const double t19 = t11*(t18 - 1); + const double t20 = t4*t9; + const double t21 = t11*(t20 - 1); + c_rho_s_rho_s = (1.0/4.0)*t0 + (1.0/4.0)*v2rho2_0 + (1.0/4.0)*v2rho2_2; + c_rho_s_rho_x = t1*t7; + c_rho_s_rho_y = rho_y*t8; + c_rho_s_rho_z = rho_z*t8; + c_rho_x_rho_x = -1.0/4.0*t0*t10 + (1.0/4.0)*t10*v2rho2_0 + (1.0/4.0)*t10*v2rho2_2 - 1.0/4.0*t12*vrho_0 + (1.0/4.0)*t12*vrho_1; + c_rho_x_rho_y = rho_y*t17; + c_rho_x_rho_z = rho_z*t17; + c_rho_y_rho_y = -1.0/4.0*t0*t18 + (1.0/4.0)*t13*t3 + (1.0/4.0)*t14*t3 - 1.0/4.0*t19*vrho_0 + (1.0/4.0)*t19*vrho_1; + c_rho_y_rho_z = (1.0/4.0)*rho_y*rho_z*t16; + c_rho_z_rho_z = -1.0/4.0*t0*t20 + (1.0/4.0)*t13*t4 + (1.0/4.0)*t14*t4 - 1.0/4.0*t21*vrho_0 + (1.0/4.0)*t21*vrho_1; +} + +/** Noncollinear GGA exchange-correlation potential: the derivative of the energy density with respect to each noncollinear field, v_X = dE/dX. */ +inline void nc_vxc_gga( [[maybe_unused]] double rho_s, double rho_x, double rho_y, double rho_z, double grad_rho_s_x, double grad_rho_s_y, double grad_rho_s_z, double grad_rho_x_x, double grad_rho_x_y, double grad_rho_x_z, double grad_rho_y_x, double grad_rho_y_y, double grad_rho_y_z, double grad_rho_z_x, double grad_rho_z_y, double grad_rho_z_z, double vrho_0, double vrho_1, double vsigma_0, double vsigma_1, double vsigma_2, double f_nabla, double& v_rho_s, double& v_rho_x, double& v_rho_y, double& v_rho_z, double& v_grad_rho_s_x, double& v_grad_rho_s_y, double& v_grad_rho_s_z, double& v_grad_rho_x_x, double& v_grad_rho_x_y, double& v_grad_rho_x_z, double& v_grad_rho_y_x, double& v_grad_rho_y_y, double& v_grad_rho_y_z, double& v_grad_rho_z_x, double& v_grad_rho_z_y, double& v_grad_rho_z_z ) { + const double t0 = (1.0/2.0)*(vrho_0 - vrho_1)/sqrt(pow(rho_x, 2) + pow(rho_y, 2) + pow(rho_z, 2)); + const double t1 = grad_rho_s_x*grad_rho_x_x + grad_rho_s_y*grad_rho_x_y + grad_rho_s_z*grad_rho_x_z; + const double t2 = grad_rho_s_x*grad_rho_y_x + grad_rho_s_y*grad_rho_y_y + grad_rho_s_z*grad_rho_y_z; + const double t3 = grad_rho_s_x*grad_rho_z_x + grad_rho_s_y*grad_rho_z_y + grad_rho_s_z*grad_rho_z_z; + const double t4 = f_nabla/sqrt(pow(t1, 2) + pow(t2, 2) + pow(t3, 2)); + const double t5 = t4*(grad_rho_x_x*t1 + grad_rho_y_x*t2 + grad_rho_z_x*t3); + const double t6 = t4*(grad_rho_x_y*t1 + grad_rho_y_y*t2 + grad_rho_z_y*t3); + const double t7 = t4*(grad_rho_x_z*t1 + grad_rho_y_z*t2 + grad_rho_z_z*t3); + const double t8 = t1*t4; + const double t9 = grad_rho_s_x*t8; + const double t10 = grad_rho_s_y*t8; + const double t11 = grad_rho_s_z*t8; + const double t12 = t2*t4; + const double t13 = grad_rho_s_x*t12; + const double t14 = grad_rho_s_y*t12; + const double t15 = grad_rho_s_z*t12; + const double t16 = t3*t4; + const double t17 = grad_rho_s_x*t16; + const double t18 = grad_rho_s_y*t16; + const double t19 = grad_rho_s_z*t16; + v_rho_s = (1.0/2.0)*(vrho_0 + vrho_1); + v_rho_x = rho_x*t0; + v_rho_y = rho_y*t0; + v_rho_z = rho_z*t0; + v_grad_rho_s_x = (1.0/2.0)*grad_rho_s_x*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_s_x + t5) - 1.0/2.0*vsigma_2*(-grad_rho_s_x + t5); + v_grad_rho_s_y = (1.0/2.0)*grad_rho_s_y*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_s_y + t6) - 1.0/2.0*vsigma_2*(-grad_rho_s_y + t6); + v_grad_rho_s_z = (1.0/2.0)*grad_rho_s_z*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_s_z + t7) - 1.0/2.0*vsigma_2*(-grad_rho_s_z + t7); + v_grad_rho_x_x = -1.0/2.0*grad_rho_x_x*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_x_x + t9) - 1.0/2.0*vsigma_2*(-grad_rho_x_x + t9); + v_grad_rho_x_y = -1.0/2.0*grad_rho_x_y*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_x_y + t10) - 1.0/2.0*vsigma_2*(-grad_rho_x_y + t10); + v_grad_rho_x_z = -1.0/2.0*grad_rho_x_z*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_x_z + t11) - 1.0/2.0*vsigma_2*(-grad_rho_x_z + t11); + v_grad_rho_y_x = -1.0/2.0*grad_rho_y_x*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_y_x + t13) - 1.0/2.0*vsigma_2*(-grad_rho_y_x + t13); + v_grad_rho_y_y = -1.0/2.0*grad_rho_y_y*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_y_y + t14) - 1.0/2.0*vsigma_2*(-grad_rho_y_y + t14); + v_grad_rho_y_z = -1.0/2.0*grad_rho_y_z*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_y_z + t15) - 1.0/2.0*vsigma_2*(-grad_rho_y_z + t15); + v_grad_rho_z_x = -1.0/2.0*grad_rho_z_x*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_z_x + t17) - 1.0/2.0*vsigma_2*(-grad_rho_z_x + t17); + v_grad_rho_z_y = -1.0/2.0*grad_rho_z_y*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_z_y + t18) - 1.0/2.0*vsigma_2*(-grad_rho_z_y + t18); + v_grad_rho_z_z = -1.0/2.0*grad_rho_z_z*vsigma_1 + (1.0/2.0)*vsigma_0*(grad_rho_z_z + t19) - 1.0/2.0*vsigma_2*(-grad_rho_z_z + t19); +} + +/** Noncollinear GGA exchange-correlation kernel applied to one trial density: k_X = sum_Y f_xc[X,Y] tY, delivered in the same field slots as the potential. */ +inline void nc_fxc_contract_gga( [[maybe_unused]] double rho_s, double rho_x, double rho_y, double rho_z, double grad_rho_s_x, double grad_rho_s_y, double grad_rho_s_z, double grad_rho_x_x, double grad_rho_x_y, double grad_rho_x_z, double grad_rho_y_x, double grad_rho_y_y, double grad_rho_y_z, double grad_rho_z_x, double grad_rho_z_y, double grad_rho_z_z, double vrho_0, double vrho_1, double vsigma_0, double vsigma_1, double vsigma_2, double v2rho2_0, double v2rho2_1, double v2rhosigma_0, double v2rhosigma_1, double v2rhosigma_2, double v2rho2_2, double v2rhosigma_3, double v2rhosigma_4, double v2rhosigma_5, double v2sigma2_0, double v2sigma2_1, double v2sigma2_2, double v2sigma2_3, double v2sigma2_4, double v2sigma2_5, double f_nabla, double trho_s, double trho_x, double trho_y, double trho_z, double tgrad_rho_s_x, double tgrad_rho_s_y, double tgrad_rho_s_z, double tgrad_rho_x_x, double tgrad_rho_x_y, double tgrad_rho_x_z, double tgrad_rho_y_x, double tgrad_rho_y_y, double tgrad_rho_y_z, double tgrad_rho_z_x, double tgrad_rho_z_y, double tgrad_rho_z_z, double& k_rho_s, double& k_rho_x, double& k_rho_y, double& k_rho_z, double& k_grad_rho_s_x, double& k_grad_rho_s_y, double& k_grad_rho_s_z, double& k_grad_rho_x_x, double& k_grad_rho_x_y, double& k_grad_rho_x_z, double& k_grad_rho_y_x, double& k_grad_rho_y_y, double& k_grad_rho_y_z, double& k_grad_rho_z_x, double& k_grad_rho_z_y, double& k_grad_rho_z_z ) { + const double t0 = 2*v2rho2_1; + const double t1 = v2rho2_0 - v2rho2_2; + const double t2 = pow(rho_x, 2); + const double t3 = pow(rho_y, 2); + const double t4 = pow(rho_z, 2); + const double t5 = t2 + t3 + t4; + const double t6 = pow(t5, -1.0/2.0); + const double t7 = grad_rho_x_x*v2rhosigma_4; + const double t8 = grad_rho_s_x*grad_rho_x_x; + const double t9 = grad_rho_s_y*grad_rho_x_y; + const double t10 = grad_rho_s_z*grad_rho_x_z; + const double t11 = t10 + t9; + const double t12 = t11 + t8; + const double t13 = pow(t12, 2); + const double t14 = grad_rho_s_x*grad_rho_y_x; + const double t15 = grad_rho_s_y*grad_rho_y_y; + const double t16 = grad_rho_s_z*grad_rho_y_z; + const double t17 = t15 + t16; + const double t18 = t14 + t17; + const double t19 = pow(t18, 2); + const double t20 = grad_rho_s_x*grad_rho_z_x; + const double t21 = grad_rho_s_y*grad_rho_z_y; + const double t22 = grad_rho_s_z*grad_rho_z_z; + const double t23 = t21 + t22; + const double t24 = t20 + t23; + const double t25 = pow(t24, 2); + const double t26 = t13 + t19 + t25; + const double t27 = f_nabla/sqrt(t26); + const double t28 = t12*t27; + const double t29 = grad_rho_s_x*t28; + const double t30 = -grad_rho_x_x + t29; + const double t31 = t30*v2rhosigma_5; + const double t32 = grad_rho_x_x + t29; + const double t33 = t32*v2rhosigma_3; + const double t34 = grad_rho_x_x*v2rhosigma_1 - t32*v2rhosigma_0; + const double t35 = t30*v2rhosigma_2 + t34; + const double t36 = t31 - t33 + t35 + t7; + const double t37 = grad_rho_x_y*v2rhosigma_4; + const double t38 = grad_rho_s_y*t28; + const double t39 = -grad_rho_x_y + t38; + const double t40 = t39*v2rhosigma_5; + const double t41 = grad_rho_x_y + t38; + const double t42 = t41*v2rhosigma_3; + const double t43 = grad_rho_x_y*v2rhosigma_1 - t41*v2rhosigma_0; + const double t44 = t39*v2rhosigma_2 + t43; + const double t45 = t37 + t40 - t42 + t44; + const double t46 = grad_rho_x_z*v2rhosigma_4; + const double t47 = grad_rho_s_z*t28; + const double t48 = -grad_rho_x_z + t47; + const double t49 = t48*v2rhosigma_5; + const double t50 = grad_rho_x_z + t47; + const double t51 = t50*v2rhosigma_3; + const double t52 = grad_rho_x_z*v2rhosigma_1 - t50*v2rhosigma_0; + const double t53 = t48*v2rhosigma_2 + t52; + const double t54 = t46 + t49 - t51 + t53; + const double t55 = grad_rho_y_x*v2rhosigma_4; + const double t56 = t18*t27; + const double t57 = grad_rho_s_x*t56; + const double t58 = -grad_rho_y_x + t57; + const double t59 = t58*v2rhosigma_5; + const double t60 = grad_rho_y_x + t57; + const double t61 = t60*v2rhosigma_3; + const double t62 = grad_rho_y_x*v2rhosigma_1 - t60*v2rhosigma_0; + const double t63 = t58*v2rhosigma_2 + t62; + const double t64 = t55 + t59 - t61 + t63; + const double t65 = grad_rho_y_y*v2rhosigma_4; + const double t66 = grad_rho_s_y*t56; + const double t67 = -grad_rho_y_y + t66; + const double t68 = t67*v2rhosigma_5; + const double t69 = grad_rho_y_y + t66; + const double t70 = t69*v2rhosigma_3; + const double t71 = grad_rho_y_y*v2rhosigma_1 - t69*v2rhosigma_0; + const double t72 = t67*v2rhosigma_2 + t71; + const double t73 = t65 + t68 - t70 + t72; + const double t74 = grad_rho_y_z*v2rhosigma_4; + const double t75 = grad_rho_s_z*t56; + const double t76 = -grad_rho_y_z + t75; + const double t77 = t76*v2rhosigma_5; + const double t78 = grad_rho_y_z + t75; + const double t79 = t78*v2rhosigma_3; + const double t80 = grad_rho_y_z*v2rhosigma_1 - t78*v2rhosigma_0; + const double t81 = t76*v2rhosigma_2 + t80; + const double t82 = t74 + t77 - t79 + t81; + const double t83 = grad_rho_z_x*v2rhosigma_4; + const double t84 = t24*t27; + const double t85 = grad_rho_s_x*t84; + const double t86 = -grad_rho_z_x + t85; + const double t87 = t86*v2rhosigma_5; + const double t88 = grad_rho_z_x + t85; + const double t89 = t88*v2rhosigma_3; + const double t90 = grad_rho_z_x*v2rhosigma_1 - t88*v2rhosigma_0; + const double t91 = t86*v2rhosigma_2 + t90; + const double t92 = t83 + t87 - t89 + t91; + const double t93 = grad_rho_z_y*v2rhosigma_4; + const double t94 = grad_rho_s_y*t84; + const double t95 = -grad_rho_z_y + t94; + const double t96 = t95*v2rhosigma_5; + const double t97 = grad_rho_z_y + t94; + const double t98 = t97*v2rhosigma_3; + const double t99 = grad_rho_z_y*v2rhosigma_1 - t97*v2rhosigma_0; + const double t100 = t95*v2rhosigma_2 + t99; + const double t101 = t100 + t93 + t96 - t98; + const double t102 = grad_rho_z_z*v2rhosigma_4; + const double t103 = grad_rho_s_z*t84; + const double t104 = -grad_rho_z_z + t103; + const double t105 = t104*v2rhosigma_5; + const double t106 = grad_rho_z_z + t103; + const double t107 = t106*v2rhosigma_3; + const double t108 = grad_rho_z_z*v2rhosigma_1 - t106*v2rhosigma_0; + const double t109 = t104*v2rhosigma_2 + t108; + const double t110 = t102 + t105 - t107 + t109; + const double t111 = grad_rho_s_x*v2rhosigma_4; + const double t112 = grad_rho_x_x*t12 + grad_rho_y_x*t18 + grad_rho_z_x*t24; + const double t113 = t112*t27; + const double t114 = grad_rho_s_x + t113; + const double t115 = t114*v2rhosigma_3; + const double t116 = -grad_rho_s_x + t113; + const double t117 = t116*v2rhosigma_5; + const double t118 = grad_rho_s_x*v2rhosigma_1 + t114*v2rhosigma_0 - t116*v2rhosigma_2; + const double t119 = t111 + t115 - t117 + t118; + const double t120 = grad_rho_s_y*v2rhosigma_4; + const double t121 = grad_rho_x_y*t12 + grad_rho_y_y*t18 + grad_rho_z_y*t24; + const double t122 = t121*t27; + const double t123 = grad_rho_s_y + t122; + const double t124 = t123*v2rhosigma_3; + const double t125 = -grad_rho_s_y + t122; + const double t126 = t125*v2rhosigma_5; + const double t127 = grad_rho_s_y*v2rhosigma_1 + t123*v2rhosigma_0 - t125*v2rhosigma_2; + const double t128 = t120 + t124 - t126 + t127; + const double t129 = grad_rho_s_z*v2rhosigma_4; + const double t130 = grad_rho_x_z*t12 + grad_rho_y_z*t18 + grad_rho_z_z*t24; + const double t131 = t130*t27; + const double t132 = grad_rho_s_z + t131; + const double t133 = t132*v2rhosigma_3; + const double t134 = -grad_rho_s_z + t131; + const double t135 = t134*v2rhosigma_5; + const double t136 = grad_rho_s_z*v2rhosigma_1 + t132*v2rhosigma_0 - t134*v2rhosigma_2; + const double t137 = t129 + t133 - t135 + t136; + const double t138 = 1.0/t5; + const double t139 = t138*v2rho2_0; + const double t140 = t138*v2rho2_2; + const double t141 = 2/pow(t5, 3.0/2.0); + const double t142 = -t0*t138 + t139 + t140 - t141*vrho_0 + t141*vrho_1; + const double t143 = t138*t2; + const double t144 = 2*t6; + const double t145 = t144*(t143 - 1); + const double t146 = rho_x*t6; + const double t147 = -t30; + const double t148 = t33 - t7; + const double t149 = tgrad_rho_x_x*(-t147*v2rhosigma_2 + t147*v2rhosigma_5 + t148 + t34); + const double t150 = -t39; + const double t151 = -t37 + t42; + const double t152 = tgrad_rho_x_y*(-t150*v2rhosigma_2 + t150*v2rhosigma_5 + t151 + t43); + const double t153 = -t48; + const double t154 = -t46 + t51; + const double t155 = tgrad_rho_x_z*(-t153*v2rhosigma_2 + t153*v2rhosigma_5 + t154 + t52); + const double t156 = -t58; + const double t157 = -t55 + t61; + const double t158 = tgrad_rho_y_x*(-t156*v2rhosigma_2 + t156*v2rhosigma_5 + t157 + t62); + const double t159 = -t67; + const double t160 = -t65 + t70; + const double t161 = tgrad_rho_y_y*(-t159*v2rhosigma_2 + t159*v2rhosigma_5 + t160 + t71); + const double t162 = -t76; + const double t163 = -t74 + t79; + const double t164 = tgrad_rho_y_z*(-t162*v2rhosigma_2 + t162*v2rhosigma_5 + t163 + t80); + const double t165 = -t86; + const double t166 = -t83 + t89; + const double t167 = tgrad_rho_z_x*(-t165*v2rhosigma_2 + t165*v2rhosigma_5 + t166 + t90); + const double t168 = -t95; + const double t169 = -t93 + t98; + const double t170 = tgrad_rho_z_y*(-t168*v2rhosigma_2 + t168*v2rhosigma_5 + t169 + t99); + const double t171 = -t104; + const double t172 = -t102 + t107; + const double t173 = tgrad_rho_z_z*(t108 - t171*v2rhosigma_2 + t171*v2rhosigma_5 + t172); + const double t174 = -t111 - t115 + t117 + t118; + const double t175 = -t120 - t124 + t126 + t127; + const double t176 = -t129 - t133 + t135 + t136; + const double t177 = t138*t3; + const double t178 = t144*(t177 - 1); + const double t179 = rho_y*t6; + const double t180 = t138*t4; + const double t181 = t144*(t180 - 1); + const double t182 = rho_z*t6; + const double t183 = 2*vsigma_1; + const double t184 = pow(grad_rho_s_x, 2); + const double t185 = grad_rho_s_x*v2sigma2_1; + const double t186 = grad_rho_s_x*v2sigma2_4; + const double t187 = 2*t116; + const double t188 = pow(grad_rho_x_x, 2); + const double t189 = pow(grad_rho_y_x, 2); + const double t190 = pow(grad_rho_z_x, 2); + const double t191 = t27*(t188 + t189 + t190); + const double t192 = f_nabla/pow(t26, 3.0/2.0); + const double t193 = pow(t112, 2)*t192; + const double t194 = 2*vsigma_0; + const double t195 = 2*vsigma_2; + const double t196 = t114*v2sigma2_2; + const double t197 = grad_rho_s_x*v2sigma2_3; + const double t198 = t114*v2sigma2_1; + const double t199 = t116*v2sigma2_4; + const double t200 = 1.0/t26; + const double t201 = t112*t200; + const double t202 = t12*t201; + const double t203 = t27*(grad_rho_x_x - t202); + const double t204 = grad_rho_s_y*t194; + const double t205 = t195*t203; + const double t206 = t114*v2sigma2_0; + const double t207 = t116*v2sigma2_2; + const double t208 = t116*v2sigma2_5; + const double t209 = grad_rho_s_y*t205 + grad_rho_x_y*t197 + grad_rho_x_y*t198 - grad_rho_x_y*t199 - t185*t41 + t186*t39 + t196*t39 - t203*t204 - t206*t41 + t207*t41 - t208*t39; + const double t210 = grad_rho_s_z*t194; + const double t211 = grad_rho_s_z*t205 + grad_rho_x_z*t197 + grad_rho_x_z*t198 - grad_rho_x_z*t199 - t185*t50 + t186*t48 + t196*t48 - t203*t210 - t206*t50 + t207*t50 - t208*t48; + const double t212 = t18*t201; + const double t213 = t27*(grad_rho_y_x - t212); + const double t214 = t195*t213; + const double t215 = grad_rho_s_y*t214 + grad_rho_y_y*t197 + grad_rho_y_y*t198 - grad_rho_y_y*t199 - t185*t69 + t186*t67 + t196*t67 - t204*t213 - t206*t69 + t207*t69 - t208*t67; + const double t216 = grad_rho_s_z*t214 + grad_rho_y_z*t197 + grad_rho_y_z*t198 - grad_rho_y_z*t199 - t185*t78 + t186*t76 + t196*t76 - t206*t78 + t207*t78 - t208*t76 - t210*t213; + const double t217 = t201*t24; + const double t218 = t27*(grad_rho_z_x - t217); + const double t219 = t195*t218; + const double t220 = grad_rho_s_y*t219 + grad_rho_z_y*t197 + grad_rho_z_y*t198 - grad_rho_z_y*t199 - t185*t97 + t186*t95 + t196*t95 - t204*t218 - t206*t97 + t207*t97 - t208*t95; + const double t221 = grad_rho_s_z*t219 + grad_rho_z_z*t197 + grad_rho_z_z*t198 - grad_rho_z_z*t199 + t104*t186 + t104*t196 - t104*t208 - t106*t185 - t106*t206 + t106*t207 - t210*t218; + const double t222 = t27*(-grad_rho_s_x*t202 + t11 + 2*t8); + const double t223 = grad_rho_x_x*t198 - grad_rho_x_x*t199 - t185*t32 + t186*t30 - t194*t222 + t195*t222 + t196*t30 - t206*t32 + t207*t32 - t208*t30 + t8*v2sigma2_3; + const double t224 = t27*(-grad_rho_s_x*t212 + 2*t14 + t17); + const double t225 = grad_rho_y_x*t198 - grad_rho_y_x*t199 + t14*v2sigma2_3 - t185*t60 + t186*t58 - t194*t224 + t195*t224 + t196*t58 - t206*t60 + t207*t60 - t208*t58; + const double t226 = t27*(-grad_rho_s_x*t217 + 2*t20 + t23); + const double t227 = grad_rho_z_x*t198 - grad_rho_z_x*t199 - t185*t88 + t186*t86 - t194*t226 + t195*t226 + t196*t86 + t20*v2sigma2_3 - t206*t88 + t207*t88 - t208*t86; + const double t228 = grad_rho_x_x*grad_rho_x_y; + const double t229 = grad_rho_y_x*grad_rho_y_y; + const double t230 = grad_rho_z_x*grad_rho_z_y; + const double t231 = t27*(-t121*t201 + t228 + t229 + t230); + const double t232 = grad_rho_s_y*t197 + grad_rho_s_y*t198 - grad_rho_s_y*t199 + t123*t185 + t123*t206 - t123*t207 - t125*t186 - t125*t196 + t125*t208 + t194*t231 - t195*t231; + const double t233 = grad_rho_x_x*grad_rho_x_z; + const double t234 = grad_rho_y_x*grad_rho_y_z; + const double t235 = grad_rho_z_x*grad_rho_z_z; + const double t236 = t27*(-t130*t201 + t233 + t234 + t235); + const double t237 = grad_rho_s_z*t197 + grad_rho_s_z*t198 - grad_rho_s_z*t199 + t132*t185 + t132*t206 - t132*t207 - t134*t186 - t134*t196 + t134*t208 + t194*t236 - t195*t236; + const double t238 = pow(grad_rho_s_y, 2); + const double t239 = grad_rho_s_y*v2sigma2_1; + const double t240 = grad_rho_s_y*v2sigma2_4; + const double t241 = 2*t125; + const double t242 = pow(grad_rho_x_y, 2); + const double t243 = pow(grad_rho_y_y, 2); + const double t244 = pow(grad_rho_z_y, 2); + const double t245 = t27*(t242 + t243 + t244); + const double t246 = pow(t121, 2)*t192; + const double t247 = t123*v2sigma2_2; + const double t248 = grad_rho_s_y*v2sigma2_3; + const double t249 = t123*v2sigma2_1; + const double t250 = t125*v2sigma2_4; + const double t251 = t121*t200; + const double t252 = t12*t251; + const double t253 = t27*(grad_rho_x_y - t252); + const double t254 = grad_rho_s_x*t253; + const double t255 = t123*v2sigma2_0; + const double t256 = t125*v2sigma2_2; + const double t257 = t125*v2sigma2_5; + const double t258 = grad_rho_x_x*t248 + grad_rho_x_x*t249 - grad_rho_x_x*t250 - t194*t254 + t195*t254 - t239*t32 + t240*t30 + t247*t30 - t255*t32 + t256*t32 - t257*t30; + const double t259 = grad_rho_s_z*t195; + const double t260 = grad_rho_x_z*t248 + grad_rho_x_z*t249 - grad_rho_x_z*t250 - t210*t253 - t239*t50 + t240*t48 + t247*t48 + t253*t259 - t255*t50 + t256*t50 - t257*t48; + const double t261 = t18*t251; + const double t262 = t27*(grad_rho_y_y - t261); + const double t263 = grad_rho_s_x*t262; + const double t264 = grad_rho_y_x*t248 + grad_rho_y_x*t249 - grad_rho_y_x*t250 - t194*t263 + t195*t263 - t239*t60 + t240*t58 + t247*t58 - t255*t60 + t256*t60 - t257*t58; + const double t265 = grad_rho_y_z*t248 + grad_rho_y_z*t249 - grad_rho_y_z*t250 - t210*t262 - t239*t78 + t240*t76 + t247*t76 - t255*t78 + t256*t78 - t257*t76 + t259*t262; + const double t266 = t24*t251; + const double t267 = t27*(grad_rho_z_y - t266); + const double t268 = grad_rho_s_x*t267; + const double t269 = grad_rho_z_x*t248 + grad_rho_z_x*t249 - grad_rho_z_x*t250 - t194*t268 + t195*t268 - t239*t88 + t240*t86 + t247*t86 - t255*t88 + t256*t88 - t257*t86; + const double t270 = grad_rho_z_z*t248 + grad_rho_z_z*t249 - grad_rho_z_z*t250 + t104*t240 + t104*t247 - t104*t257 - t106*t239 - t106*t255 + t106*t256 - t210*t267 + t259*t267; + const double t271 = t27*(-grad_rho_s_y*t252 + t10 + t8 + 2*t9); + const double t272 = grad_rho_x_y*t249 - grad_rho_x_y*t250 - t194*t271 + t195*t271 - t239*t41 + t240*t39 + t247*t39 - t255*t41 + t256*t41 - t257*t39 + t9*v2sigma2_3; + const double t273 = t27*(-grad_rho_s_y*t261 + t14 + 2*t15 + t16); + const double t274 = grad_rho_y_y*t249 - grad_rho_y_y*t250 + t15*v2sigma2_3 - t194*t273 + t195*t273 - t239*t69 + t240*t67 + t247*t67 - t255*t69 + t256*t69 - t257*t67; + const double t275 = t27*(-grad_rho_s_y*t266 + t20 + 2*t21 + t22); + const double t276 = grad_rho_z_y*t249 - grad_rho_z_y*t250 - t194*t275 + t195*t275 + t21*v2sigma2_3 - t239*t97 + t240*t95 + t247*t95 - t255*t97 + t256*t97 - t257*t95; + const double t277 = grad_rho_x_y*grad_rho_x_z; + const double t278 = grad_rho_y_y*grad_rho_y_z; + const double t279 = grad_rho_z_y*grad_rho_z_z; + const double t280 = t27*(-t130*t251 + t277 + t278 + t279); + const double t281 = grad_rho_s_z*t248 + grad_rho_s_z*t249 - grad_rho_s_z*t250 + t132*t239 + t132*t255 - t132*t256 - t134*t240 - t134*t247 + t134*t257 + t194*t280 - t195*t280; + const double t282 = pow(grad_rho_s_z, 2); + const double t283 = t132*v2sigma2_1; + const double t284 = grad_rho_s_z*v2sigma2_4; + const double t285 = 2*t134; + const double t286 = pow(grad_rho_x_z, 2); + const double t287 = pow(grad_rho_y_z, 2); + const double t288 = pow(grad_rho_z_z, 2); + const double t289 = t27*(t286 + t287 + t288); + const double t290 = pow(t130, 2)*t192; + const double t291 = t132*v2sigma2_2; + const double t292 = grad_rho_s_z*v2sigma2_3; + const double t293 = grad_rho_s_z*v2sigma2_1; + const double t294 = t134*v2sigma2_4; + const double t295 = t130*t200; + const double t296 = t12*t295; + const double t297 = t27*(grad_rho_x_z - t296); + const double t298 = grad_rho_s_x*t297; + const double t299 = t132*v2sigma2_0; + const double t300 = t134*v2sigma2_2; + const double t301 = t134*v2sigma2_5; + const double t302 = grad_rho_x_x*t283 + grad_rho_x_x*t292 - grad_rho_x_x*t294 - t194*t298 + t195*t298 + t284*t30 + t291*t30 - t293*t32 - t299*t32 - t30*t301 + t300*t32; + const double t303 = grad_rho_s_y*t195; + const double t304 = grad_rho_x_y*t283 + grad_rho_x_y*t292 - grad_rho_x_y*t294 - t204*t297 + t284*t39 + t291*t39 - t293*t41 + t297*t303 - t299*t41 + t300*t41 - t301*t39; + const double t305 = t18*t295; + const double t306 = t27*(grad_rho_y_z - t305); + const double t307 = grad_rho_s_x*t306; + const double t308 = grad_rho_y_x*t283 + grad_rho_y_x*t292 - grad_rho_y_x*t294 - t194*t307 + t195*t307 + t284*t58 + t291*t58 - t293*t60 - t299*t60 + t300*t60 - t301*t58; + const double t309 = grad_rho_y_y*t283 + grad_rho_y_y*t292 - grad_rho_y_y*t294 - t204*t306 + t284*t67 + t291*t67 - t293*t69 - t299*t69 + t300*t69 - t301*t67 + t303*t306; + const double t310 = t24*t295; + const double t311 = t27*(grad_rho_z_z - t310); + const double t312 = grad_rho_s_x*t311; + const double t313 = grad_rho_z_x*t283 + grad_rho_z_x*t292 - grad_rho_z_x*t294 - t194*t312 + t195*t312 + t284*t86 + t291*t86 - t293*t88 - t299*t88 + t300*t88 - t301*t86; + const double t314 = grad_rho_z_y*t283 + grad_rho_z_y*t292 - grad_rho_z_y*t294 - t204*t311 + t284*t95 + t291*t95 - t293*t97 - t299*t97 + t300*t97 - t301*t95 + t303*t311; + const double t315 = t27*(-grad_rho_s_z*t296 + 2*t10 + t8 + t9); + const double t316 = grad_rho_x_z*t283 - grad_rho_x_z*t294 + t10*v2sigma2_3 - t194*t315 + t195*t315 + t284*t48 + t291*t48 - t293*t50 - t299*t50 + t300*t50 - t301*t48; + const double t317 = t27*(-grad_rho_s_z*t305 + t14 + t15 + 2*t16); + const double t318 = grad_rho_y_z*t283 - grad_rho_y_z*t294 + t16*v2sigma2_3 - t194*t317 + t195*t317 + t284*t76 + t291*t76 - t293*t78 - t299*t78 + t300*t78 - t301*t76; + const double t319 = t27*(-grad_rho_s_z*t310 + t20 + t21 + 2*t22); + const double t320 = grad_rho_z_z*t283 - grad_rho_z_z*t294 + t104*t284 + t104*t291 - t104*t301 - t106*t293 - t106*t299 + t106*t300 - t194*t319 + t195*t319 + t22*v2sigma2_3; + const double t321 = rho_x*trho_x; + const double t322 = t6*(t148 - t31 + t35); + const double t323 = rho_y*trho_y; + const double t324 = rho_z*trho_z; + const double t325 = -t183; + const double t326 = grad_rho_x_x*v2sigma2_1; + const double t327 = grad_rho_x_x*v2sigma2_4; + const double t328 = 2*t30; + const double t329 = t13*t192; + const double t330 = t184*t329; + const double t331 = t184*t27; + const double t332 = t331 + 1; + const double t333 = 1 - t331; + const double t334 = t32*v2sigma2_2; + const double t335 = grad_rho_x_x*v2sigma2_3; + const double t336 = t30*v2sigma2_4; + const double t337 = t32*v2sigma2_0; + const double t338 = t30*v2sigma2_5; + const double t339 = t32*v2sigma2_1; + const double t340 = t30*v2sigma2_2; + const double t341 = grad_rho_s_x*t204; + const double t342 = t12*t192; + const double t343 = t18*t342; + const double t344 = grad_rho_s_x*t343; + const double t345 = t303*t344 - t341*t343; + const double t346 = grad_rho_y_y*t335 + grad_rho_y_y*t336 - grad_rho_y_y*t339 - t326*t69 + t327*t67 - t334*t67 + t337*t69 + t338*t67 - t340*t69 + t345; + const double t347 = -t210*t344 + t259*t344; + const double t348 = grad_rho_y_z*t335 + grad_rho_y_z*t336 - grad_rho_y_z*t339 - t326*t78 + t327*t76 - t334*t76 + t337*t78 + t338*t76 - t340*t78 + t347; + const double t349 = t24*t342; + const double t350 = grad_rho_s_x*t349; + const double t351 = t303*t350 - t341*t349; + const double t352 = grad_rho_z_y*t335 + grad_rho_z_y*t336 - grad_rho_z_y*t339 - t326*t97 + t327*t95 - t334*t95 + t337*t97 + t338*t95 - t340*t97 + t351; + const double t353 = -t210*t350 + t259*t350; + const double t354 = grad_rho_z_z*t335 + grad_rho_z_z*t336 - grad_rho_z_z*t339 + t104*t327 - t104*t334 + t104*t338 - t106*t326 + t106*t337 - t106*t340 + t353; + const double t355 = t184*t342; + const double t356 = t18*t194; + const double t357 = t195*t355; + const double t358 = grad_rho_y_x*t335 + grad_rho_y_x*t336 - grad_rho_y_x*t339 + t18*t357 - t326*t60 + t327*t58 - t334*t58 + t337*t60 + t338*t58 - t340*t60 - t355*t356; + const double t359 = t194*t24; + const double t360 = grad_rho_z_x*t335 + grad_rho_z_x*t336 - grad_rho_z_x*t339 + t24*t357 - t326*t88 + t327*t86 - t334*t86 + t337*t88 + t338*t86 - t340*t88 - t355*t359; + const double t361 = t27*(t13*t200 - 1); + const double t362 = grad_rho_s_x*t361; + const double t363 = grad_rho_x_y*t336 - grad_rho_x_y*t339 + t228*v2sigma2_3 + t303*t362 - t326*t41 + t327*t39 - t334*t39 + t337*t41 + t338*t39 - t340*t41 - t341*t361; + const double t364 = grad_rho_x_z*t336 - grad_rho_x_z*t339 - t210*t362 + t233*v2sigma2_3 + t259*t362 - t326*t50 + t327*t48 - t334*t48 + t337*t50 + t338*t48 - t340*t50; + const double t365 = t6*(t151 - t40 + t44); + const double t366 = grad_rho_x_y*v2sigma2_1; + const double t367 = grad_rho_x_y*v2sigma2_4; + const double t368 = 2*t39; + const double t369 = t238*t329; + const double t370 = t238*t27; + const double t371 = t370 + 1; + const double t372 = 1 - t370; + const double t373 = t41*v2sigma2_2; + const double t374 = grad_rho_x_y*v2sigma2_3; + const double t375 = t39*v2sigma2_4; + const double t376 = t41*v2sigma2_0; + const double t377 = t39*v2sigma2_5; + const double t378 = t41*v2sigma2_1; + const double t379 = t39*v2sigma2_2; + const double t380 = grad_rho_y_x*t374 + grad_rho_y_x*t375 - grad_rho_y_x*t378 + t345 - t366*t60 + t367*t58 - t373*t58 + t376*t60 + t377*t58 - t379*t60; + const double t381 = grad_rho_s_z*t204; + const double t382 = grad_rho_s_y*t259; + const double t383 = -t343*t381 + t343*t382; + const double t384 = grad_rho_y_z*t374 + grad_rho_y_z*t375 - grad_rho_y_z*t378 - t366*t78 + t367*t76 - t373*t76 + t376*t78 + t377*t76 - t379*t78 + t383; + const double t385 = grad_rho_z_x*t374 + grad_rho_z_x*t375 - grad_rho_z_x*t378 + t351 - t366*t88 + t367*t86 - t373*t86 + t376*t88 + t377*t86 - t379*t88; + const double t386 = -t349*t381 + t349*t382; + const double t387 = grad_rho_z_z*t374 + grad_rho_z_z*t375 - grad_rho_z_z*t378 + t104*t367 - t104*t373 + t104*t377 - t106*t366 + t106*t376 - t106*t379 + t386; + const double t388 = t238*t342; + const double t389 = t195*t388; + const double t390 = grad_rho_y_y*t374 + grad_rho_y_y*t375 - grad_rho_y_y*t378 + t18*t389 - t356*t388 - t366*t69 + t367*t67 - t373*t67 + t376*t69 + t377*t67 - t379*t69; + const double t391 = grad_rho_z_y*t374 + grad_rho_z_y*t375 - grad_rho_z_y*t378 + t24*t389 - t359*t388 - t366*t97 + t367*t95 - t373*t95 + t376*t97 + t377*t95 - t379*t97; + const double t392 = grad_rho_x_z*t375 - grad_rho_x_z*t378 + t277*v2sigma2_3 - t361*t381 + t361*t382 - t366*t50 + t367*t48 - t373*t48 + t376*t50 + t377*t48 - t379*t50; + const double t393 = t6*(t154 - t49 + t53); + const double t394 = grad_rho_x_z*v2sigma2_1; + const double t395 = grad_rho_x_z*v2sigma2_4; + const double t396 = 2*t48; + const double t397 = t282*t329; + const double t398 = t27*t282; + const double t399 = t398 + 1; + const double t400 = 1 - t398; + const double t401 = t50*v2sigma2_2; + const double t402 = grad_rho_x_z*v2sigma2_3; + const double t403 = t48*v2sigma2_4; + const double t404 = t50*v2sigma2_0; + const double t405 = t48*v2sigma2_5; + const double t406 = t50*v2sigma2_1; + const double t407 = t48*v2sigma2_2; + const double t408 = grad_rho_y_x*t402 + grad_rho_y_x*t403 - grad_rho_y_x*t406 + t347 - t394*t60 + t395*t58 - t401*t58 + t404*t60 + t405*t58 - t407*t60; + const double t409 = grad_rho_y_y*t402 + grad_rho_y_y*t403 - grad_rho_y_y*t406 + t383 - t394*t69 + t395*t67 - t401*t67 + t404*t69 + t405*t67 - t407*t69; + const double t410 = grad_rho_z_x*t402 + grad_rho_z_x*t403 - grad_rho_z_x*t406 + t353 - t394*t88 + t395*t86 - t401*t86 + t404*t88 + t405*t86 - t407*t88; + const double t411 = grad_rho_z_y*t402 + grad_rho_z_y*t403 - grad_rho_z_y*t406 + t386 - t394*t97 + t395*t95 - t401*t95 + t404*t97 + t405*t95 - t407*t97; + const double t412 = t282*t342; + const double t413 = t195*t412; + const double t414 = grad_rho_y_z*t402 + grad_rho_y_z*t403 - grad_rho_y_z*t406 + t18*t413 - t356*t412 - t394*t78 + t395*t76 - t401*t76 + t404*t78 + t405*t76 - t407*t78; + const double t415 = grad_rho_z_z*t402 + grad_rho_z_z*t403 - grad_rho_z_z*t406 + t104*t395 - t104*t401 + t104*t405 - t106*t394 + t106*t404 - t106*t407 + t24*t413 - t359*t412; + const double t416 = t6*(t157 - t59 + t63); + const double t417 = grad_rho_y_x*v2sigma2_1; + const double t418 = grad_rho_y_x*v2sigma2_4; + const double t419 = 2*t58; + const double t420 = t19*t192; + const double t421 = t184*t420; + const double t422 = t60*v2sigma2_2; + const double t423 = grad_rho_y_x*v2sigma2_3; + const double t424 = t58*v2sigma2_4; + const double t425 = t60*v2sigma2_0; + const double t426 = t58*v2sigma2_5; + const double t427 = t60*v2sigma2_1; + const double t428 = t58*v2sigma2_2; + const double t429 = t192*t24; + const double t430 = t18*t429; + const double t431 = grad_rho_s_x*t430; + const double t432 = t303*t431 - t341*t430; + const double t433 = grad_rho_z_y*t423 + grad_rho_z_y*t424 - grad_rho_z_y*t427 - t417*t97 + t418*t95 - t422*t95 + t425*t97 + t426*t95 - t428*t97 + t432; + const double t434 = -t210*t431 + t259*t431; + const double t435 = grad_rho_z_z*t423 + grad_rho_z_z*t424 - grad_rho_z_z*t427 + t104*t418 - t104*t422 + t104*t426 - t106*t417 + t106*t425 - t106*t428 + t434; + const double t436 = t184*t429; + const double t437 = t18*t195; + const double t438 = grad_rho_z_x*t423 + grad_rho_z_x*t424 - grad_rho_z_x*t427 - t356*t436 - t417*t88 + t418*t86 - t422*t86 + t425*t88 + t426*t86 - t428*t88 + t436*t437; + const double t439 = t27*(t19*t200 - 1); + const double t440 = grad_rho_s_x*t439; + const double t441 = grad_rho_y_y*t424 - grad_rho_y_y*t427 + t229*v2sigma2_3 + t303*t440 - t341*t439 - t417*t69 + t418*t67 - t422*t67 + t425*t69 + t426*t67 - t428*t69; + const double t442 = grad_rho_y_z*t424 - grad_rho_y_z*t427 - t210*t440 + t234*v2sigma2_3 + t259*t440 - t417*t78 + t418*t76 - t422*t76 + t425*t78 + t426*t76 - t428*t78; + const double t443 = t6*(t160 - t68 + t72); + const double t444 = grad_rho_y_y*v2sigma2_1; + const double t445 = grad_rho_y_y*v2sigma2_4; + const double t446 = 2*t67; + const double t447 = t238*t420; + const double t448 = t69*v2sigma2_2; + const double t449 = grad_rho_y_y*v2sigma2_3; + const double t450 = t67*v2sigma2_4; + const double t451 = t69*v2sigma2_0; + const double t452 = t67*v2sigma2_5; + const double t453 = t69*v2sigma2_1; + const double t454 = t67*v2sigma2_2; + const double t455 = grad_rho_z_x*t449 + grad_rho_z_x*t450 - grad_rho_z_x*t453 + t432 - t444*t88 + t445*t86 - t448*t86 + t451*t88 + t452*t86 - t454*t88; + const double t456 = -t381*t430 + t382*t430; + const double t457 = grad_rho_z_z*t449 + grad_rho_z_z*t450 - grad_rho_z_z*t453 + t104*t445 - t104*t448 + t104*t452 - t106*t444 + t106*t451 - t106*t454 + t456; + const double t458 = t238*t429; + const double t459 = grad_rho_z_y*t449 + grad_rho_z_y*t450 - grad_rho_z_y*t453 - t356*t458 + t437*t458 - t444*t97 + t445*t95 - t448*t95 + t451*t97 + t452*t95 - t454*t97; + const double t460 = grad_rho_y_z*t450 - grad_rho_y_z*t453 + t278*v2sigma2_3 - t381*t439 + t382*t439 - t444*t78 + t445*t76 - t448*t76 + t451*t78 + t452*t76 - t454*t78; + const double t461 = t6*(t163 - t77 + t81); + const double t462 = grad_rho_y_z*v2sigma2_1; + const double t463 = grad_rho_y_z*v2sigma2_4; + const double t464 = 2*t76; + const double t465 = t282*t420; + const double t466 = t78*v2sigma2_2; + const double t467 = grad_rho_y_z*v2sigma2_3; + const double t468 = t76*v2sigma2_4; + const double t469 = t78*v2sigma2_0; + const double t470 = t76*v2sigma2_5; + const double t471 = t78*v2sigma2_1; + const double t472 = t76*v2sigma2_2; + const double t473 = grad_rho_z_x*t467 + grad_rho_z_x*t468 - grad_rho_z_x*t471 + t434 - t462*t88 + t463*t86 - t466*t86 + t469*t88 + t470*t86 - t472*t88; + const double t474 = grad_rho_z_y*t467 + grad_rho_z_y*t468 - grad_rho_z_y*t471 + t456 - t462*t97 + t463*t95 - t466*t95 + t469*t97 + t470*t95 - t472*t97; + const double t475 = t282*t429; + const double t476 = grad_rho_z_z*t467 + grad_rho_z_z*t468 - grad_rho_z_z*t471 + t104*t463 - t104*t466 + t104*t470 - t106*t462 + t106*t469 - t106*t472 - t356*t475 + t437*t475; + const double t477 = t6*(t166 - t87 + t91); + const double t478 = grad_rho_z_x*v2sigma2_1; + const double t479 = grad_rho_z_x*v2sigma2_4; + const double t480 = 2*t86; + const double t481 = t192*t25; + const double t482 = t184*t481; + const double t483 = t88*v2sigma2_2; + const double t484 = t88*v2sigma2_1; + const double t485 = t86*v2sigma2_4; + const double t486 = t27*(t200*t25 - 1); + const double t487 = grad_rho_s_x*t486; + const double t488 = t88*v2sigma2_0; + const double t489 = t86*v2sigma2_2; + const double t490 = t86*v2sigma2_5; + const double t491 = -grad_rho_z_y*t484 + grad_rho_z_y*t485 + t230*v2sigma2_3 + t303*t487 - t341*t486 - t478*t97 + t479*t95 - t483*t95 + t488*t97 - t489*t97 + t490*t95; + const double t492 = -grad_rho_z_z*t484 + grad_rho_z_z*t485 + t104*t479 - t104*t483 + t104*t490 - t106*t478 + t106*t488 - t106*t489 - t210*t487 + t235*v2sigma2_3 + t259*t487; + const double t493 = t6*(t100 + t169 - t96); + const double t494 = grad_rho_z_y*v2sigma2_1; + const double t495 = grad_rho_z_y*v2sigma2_4; + const double t496 = 2*t95; + const double t497 = t238*t481; + const double t498 = t97*v2sigma2_2; + const double t499 = grad_rho_z_z*v2sigma2_1; + const double t500 = grad_rho_z_z*v2sigma2_4; + const double t501 = t106*v2sigma2_2; + const double t502 = t104*t495 - t104*t498 + t104*t95*v2sigma2_5 - t106*t494 + t106*t97*v2sigma2_0 + t279*v2sigma2_3 - t381*t486 + t382*t486 - t499*t97 + t500*t95 - t501*t95; + const double t503 = t6*(-t105 + t109 + t172); + const double t504 = 2*t104; + const double t505 = t282*t481; + k_rho_s = (1.0/4.0)*rho_x*t1*t6*trho_x + (1.0/4.0)*rho_y*t1*t6*trho_y + (1.0/4.0)*rho_z*t1*t6*trho_z - 1.0/4.0*t101*tgrad_rho_z_y - 1.0/4.0*t110*tgrad_rho_z_z + (1.0/4.0)*t119*tgrad_rho_s_x + (1.0/4.0)*t128*tgrad_rho_s_y + (1.0/4.0)*t137*tgrad_rho_s_z - 1.0/4.0*t36*tgrad_rho_x_x - 1.0/4.0*t45*tgrad_rho_x_y - 1.0/4.0*t54*tgrad_rho_x_z - 1.0/4.0*t64*tgrad_rho_y_x - 1.0/4.0*t73*tgrad_rho_y_y - 1.0/4.0*t82*tgrad_rho_y_z - 1.0/4.0*t92*tgrad_rho_z_x + (1.0/4.0)*trho_s*(t0 + v2rho2_0 + v2rho2_2); + k_rho_x = (1.0/4.0)*rho_x*rho_y*t142*trho_y + (1.0/4.0)*rho_x*rho_z*t142*trho_z + (1.0/4.0)*rho_x*t1*t6*trho_s + (1.0/4.0)*rho_x*t174*t6*tgrad_rho_s_x + (1.0/4.0)*rho_x*t175*t6*tgrad_rho_s_y + (1.0/4.0)*rho_x*t176*t6*tgrad_rho_s_z - 1.0/4.0*t146*t149 - 1.0/4.0*t146*t152 - 1.0/4.0*t146*t155 - 1.0/4.0*t146*t158 - 1.0/4.0*t146*t161 - 1.0/4.0*t146*t164 - 1.0/4.0*t146*t167 - 1.0/4.0*t146*t170 - 1.0/4.0*t146*t173 + (1.0/4.0)*trho_x*(-t0*t143 + t143*v2rho2_0 + t143*v2rho2_2 - t145*vrho_0 + t145*vrho_1); + k_rho_y = (1.0/4.0)*rho_x*rho_y*t142*trho_x + (1.0/4.0)*rho_y*rho_z*t142*trho_z + (1.0/4.0)*rho_y*t1*t6*trho_s + (1.0/4.0)*rho_y*t174*t6*tgrad_rho_s_x + (1.0/4.0)*rho_y*t175*t6*tgrad_rho_s_y + (1.0/4.0)*rho_y*t176*t6*tgrad_rho_s_z - 1.0/4.0*t149*t179 - 1.0/4.0*t152*t179 - 1.0/4.0*t155*t179 - 1.0/4.0*t158*t179 - 1.0/4.0*t161*t179 - 1.0/4.0*t164*t179 - 1.0/4.0*t167*t179 - 1.0/4.0*t170*t179 - 1.0/4.0*t173*t179 + (1.0/4.0)*trho_y*(-t0*t177 + t139*t3 + t140*t3 - t178*vrho_0 + t178*vrho_1); + k_rho_z = (1.0/4.0)*rho_x*rho_z*t142*trho_x + (1.0/4.0)*rho_y*rho_z*t142*trho_y + (1.0/4.0)*rho_z*t1*t6*trho_s + (1.0/4.0)*rho_z*t174*t6*tgrad_rho_s_x + (1.0/4.0)*rho_z*t175*t6*tgrad_rho_s_y + (1.0/4.0)*rho_z*t176*t6*tgrad_rho_s_z - 1.0/4.0*t149*t182 - 1.0/4.0*t152*t182 - 1.0/4.0*t155*t182 - 1.0/4.0*t158*t182 - 1.0/4.0*t161*t182 - 1.0/4.0*t164*t182 - 1.0/4.0*t167*t182 - 1.0/4.0*t170*t182 - 1.0/4.0*t173*t182 + (1.0/4.0)*trho_z*(-t0*t180 + t139*t4 + t140*t4 - t181*vrho_0 + t181*vrho_1); + k_grad_rho_s_x = (1.0/4.0)*rho_x*t174*t6*trho_x + (1.0/4.0)*rho_y*t174*t6*trho_y + (1.0/4.0)*rho_z*t174*t6*trho_z + (1.0/4.0)*t119*trho_s - 1.0/4.0*t209*tgrad_rho_x_y - 1.0/4.0*t211*tgrad_rho_x_z - 1.0/4.0*t215*tgrad_rho_y_y - 1.0/4.0*t216*tgrad_rho_y_z - 1.0/4.0*t220*tgrad_rho_z_y - 1.0/4.0*t221*tgrad_rho_z_z - 1.0/4.0*t223*tgrad_rho_x_x - 1.0/4.0*t225*tgrad_rho_y_x - 1.0/4.0*t227*tgrad_rho_z_x + (1.0/4.0)*t232*tgrad_rho_s_y + (1.0/4.0)*t237*tgrad_rho_s_z + (1.0/4.0)*tgrad_rho_s_x*(pow(t114, 2)*v2sigma2_0 + 2*t114*t185 + pow(t116, 2)*v2sigma2_5 + t183 + t184*v2sigma2_3 - t186*t187 - t187*t196 + t194*(t191 - t193 + 1) + t195*(-t191 + t193 + 1)); + k_grad_rho_s_y = (1.0/4.0)*rho_x*t175*t6*trho_x + (1.0/4.0)*rho_y*t175*t6*trho_y + (1.0/4.0)*rho_z*t175*t6*trho_z + (1.0/4.0)*t128*trho_s + (1.0/4.0)*t232*tgrad_rho_s_x - 1.0/4.0*t258*tgrad_rho_x_x - 1.0/4.0*t260*tgrad_rho_x_z - 1.0/4.0*t264*tgrad_rho_y_x - 1.0/4.0*t265*tgrad_rho_y_z - 1.0/4.0*t269*tgrad_rho_z_x - 1.0/4.0*t270*tgrad_rho_z_z - 1.0/4.0*t272*tgrad_rho_x_y - 1.0/4.0*t274*tgrad_rho_y_y - 1.0/4.0*t276*tgrad_rho_z_y + (1.0/4.0)*t281*tgrad_rho_s_z + (1.0/4.0)*tgrad_rho_s_y*(pow(t123, 2)*v2sigma2_0 + 2*t123*t239 + pow(t125, 2)*v2sigma2_5 + t183 + t194*(t245 - t246 + 1) + t195*(-t245 + t246 + 1) + t238*v2sigma2_3 - t240*t241 - t241*t247); + k_grad_rho_s_z = (1.0/4.0)*rho_x*t176*t6*trho_x + (1.0/4.0)*rho_y*t176*t6*trho_y + (1.0/4.0)*rho_z*t176*t6*trho_z + (1.0/4.0)*t137*trho_s + (1.0/4.0)*t237*tgrad_rho_s_x + (1.0/4.0)*t281*tgrad_rho_s_y - 1.0/4.0*t302*tgrad_rho_x_x - 1.0/4.0*t304*tgrad_rho_x_y - 1.0/4.0*t308*tgrad_rho_y_x - 1.0/4.0*t309*tgrad_rho_y_y - 1.0/4.0*t313*tgrad_rho_z_x - 1.0/4.0*t314*tgrad_rho_z_y - 1.0/4.0*t316*tgrad_rho_x_z - 1.0/4.0*t318*tgrad_rho_y_z - 1.0/4.0*t320*tgrad_rho_z_z + (1.0/4.0)*tgrad_rho_s_z*(2*grad_rho_s_z*t283 + pow(t132, 2)*v2sigma2_0 + pow(t134, 2)*v2sigma2_5 + t183 + t194*(t289 - t290 + 1) + t195*(-t289 + t290 + 1) + t282*v2sigma2_3 - t284*t285 - t285*t291); + k_grad_rho_x_x = -1.0/4.0*t223*tgrad_rho_s_x - 1.0/4.0*t258*tgrad_rho_s_y - 1.0/4.0*t302*tgrad_rho_s_z - 1.0/4.0*t321*t322 - 1.0/4.0*t322*t323 - 1.0/4.0*t322*t324 + (1.0/4.0)*t346*tgrad_rho_y_y + (1.0/4.0)*t348*tgrad_rho_y_z + (1.0/4.0)*t352*tgrad_rho_z_y + (1.0/4.0)*t354*tgrad_rho_z_z + (1.0/4.0)*t358*tgrad_rho_y_x - 1.0/4.0*t36*trho_s + (1.0/4.0)*t360*tgrad_rho_z_x + (1.0/4.0)*t363*tgrad_rho_x_y + (1.0/4.0)*t364*tgrad_rho_x_z + (1.0/4.0)*tgrad_rho_x_x*(t188*v2sigma2_3 + t194*(-t330 + t332) + t195*(t330 + t333) + pow(t30, 2)*v2sigma2_5 + pow(t32, 2)*v2sigma2_0 - 2*t32*t326 + t325 + t327*t328 - t328*t334); + k_grad_rho_x_y = -1.0/4.0*t209*tgrad_rho_s_x - 1.0/4.0*t272*tgrad_rho_s_y - 1.0/4.0*t304*tgrad_rho_s_z - 1.0/4.0*t321*t365 - 1.0/4.0*t323*t365 - 1.0/4.0*t324*t365 + (1.0/4.0)*t363*tgrad_rho_x_x + (1.0/4.0)*t380*tgrad_rho_y_x + (1.0/4.0)*t384*tgrad_rho_y_z + (1.0/4.0)*t385*tgrad_rho_z_x + (1.0/4.0)*t387*tgrad_rho_z_z + (1.0/4.0)*t390*tgrad_rho_y_y + (1.0/4.0)*t391*tgrad_rho_z_y + (1.0/4.0)*t392*tgrad_rho_x_z - 1.0/4.0*t45*trho_s + (1.0/4.0)*tgrad_rho_x_y*(t194*(-t369 + t371) + t195*(t369 + t372) + t242*v2sigma2_3 + t325 - 2*t366*t41 + t367*t368 - t368*t373 + pow(t39, 2)*v2sigma2_5 + pow(t41, 2)*v2sigma2_0); + k_grad_rho_x_z = -1.0/4.0*t211*tgrad_rho_s_x - 1.0/4.0*t260*tgrad_rho_s_y - 1.0/4.0*t316*tgrad_rho_s_z - 1.0/4.0*t321*t393 - 1.0/4.0*t323*t393 - 1.0/4.0*t324*t393 + (1.0/4.0)*t364*tgrad_rho_x_x + (1.0/4.0)*t392*tgrad_rho_x_y + (1.0/4.0)*t408*tgrad_rho_y_x + (1.0/4.0)*t409*tgrad_rho_y_y + (1.0/4.0)*t410*tgrad_rho_z_x + (1.0/4.0)*t411*tgrad_rho_z_y + (1.0/4.0)*t414*tgrad_rho_y_z + (1.0/4.0)*t415*tgrad_rho_z_z - 1.0/4.0*t54*trho_s + (1.0/4.0)*tgrad_rho_x_z*(t194*(-t397 + t399) + t195*(t397 + t400) + t286*v2sigma2_3 + t325 - 2*t394*t50 + t395*t396 - t396*t401 + pow(t48, 2)*v2sigma2_5 + pow(t50, 2)*v2sigma2_0); + k_grad_rho_y_x = -1.0/4.0*t225*tgrad_rho_s_x - 1.0/4.0*t264*tgrad_rho_s_y - 1.0/4.0*t308*tgrad_rho_s_z - 1.0/4.0*t321*t416 - 1.0/4.0*t323*t416 - 1.0/4.0*t324*t416 + (1.0/4.0)*t358*tgrad_rho_x_x + (1.0/4.0)*t380*tgrad_rho_x_y + (1.0/4.0)*t408*tgrad_rho_x_z + (1.0/4.0)*t433*tgrad_rho_z_y + (1.0/4.0)*t435*tgrad_rho_z_z + (1.0/4.0)*t438*tgrad_rho_z_x + (1.0/4.0)*t441*tgrad_rho_y_y + (1.0/4.0)*t442*tgrad_rho_y_z - 1.0/4.0*t64*trho_s + (1.0/4.0)*tgrad_rho_y_x*(t189*v2sigma2_3 + t194*(t332 - t421) + t195*(t333 + t421) + t325 - 2*t417*t60 + t418*t419 - t419*t422 + pow(t58, 2)*v2sigma2_5 + pow(t60, 2)*v2sigma2_0); + k_grad_rho_y_y = -1.0/4.0*t215*tgrad_rho_s_x - 1.0/4.0*t274*tgrad_rho_s_y - 1.0/4.0*t309*tgrad_rho_s_z - 1.0/4.0*t321*t443 - 1.0/4.0*t323*t443 - 1.0/4.0*t324*t443 + (1.0/4.0)*t346*tgrad_rho_x_x + (1.0/4.0)*t390*tgrad_rho_x_y + (1.0/4.0)*t409*tgrad_rho_x_z + (1.0/4.0)*t441*tgrad_rho_y_x + (1.0/4.0)*t455*tgrad_rho_z_x + (1.0/4.0)*t457*tgrad_rho_z_z + (1.0/4.0)*t459*tgrad_rho_z_y + (1.0/4.0)*t460*tgrad_rho_y_z - 1.0/4.0*t73*trho_s + (1.0/4.0)*tgrad_rho_y_y*(t194*(t371 - t447) + t195*(t372 + t447) + t243*v2sigma2_3 + t325 - 2*t444*t69 + t445*t446 - t446*t448 + pow(t67, 2)*v2sigma2_5 + pow(t69, 2)*v2sigma2_0); + k_grad_rho_y_z = -1.0/4.0*t216*tgrad_rho_s_x - 1.0/4.0*t265*tgrad_rho_s_y - 1.0/4.0*t318*tgrad_rho_s_z - 1.0/4.0*t321*t461 - 1.0/4.0*t323*t461 - 1.0/4.0*t324*t461 + (1.0/4.0)*t348*tgrad_rho_x_x + (1.0/4.0)*t384*tgrad_rho_x_y + (1.0/4.0)*t414*tgrad_rho_x_z + (1.0/4.0)*t442*tgrad_rho_y_x + (1.0/4.0)*t460*tgrad_rho_y_y + (1.0/4.0)*t473*tgrad_rho_z_x + (1.0/4.0)*t474*tgrad_rho_z_y + (1.0/4.0)*t476*tgrad_rho_z_z - 1.0/4.0*t82*trho_s + (1.0/4.0)*tgrad_rho_y_z*(t194*(t399 - t465) + t195*(t400 + t465) + t287*v2sigma2_3 + t325 - 2*t462*t78 + t463*t464 - t464*t466 + pow(t76, 2)*v2sigma2_5 + pow(t78, 2)*v2sigma2_0); + k_grad_rho_z_x = -1.0/4.0*t227*tgrad_rho_s_x - 1.0/4.0*t269*tgrad_rho_s_y - 1.0/4.0*t313*tgrad_rho_s_z - 1.0/4.0*t321*t477 - 1.0/4.0*t323*t477 - 1.0/4.0*t324*t477 + (1.0/4.0)*t360*tgrad_rho_x_x + (1.0/4.0)*t385*tgrad_rho_x_y + (1.0/4.0)*t410*tgrad_rho_x_z + (1.0/4.0)*t438*tgrad_rho_y_x + (1.0/4.0)*t455*tgrad_rho_y_y + (1.0/4.0)*t473*tgrad_rho_y_z + (1.0/4.0)*t491*tgrad_rho_z_y + (1.0/4.0)*t492*tgrad_rho_z_z - 1.0/4.0*t92*trho_s + (1.0/4.0)*tgrad_rho_z_x*(t190*v2sigma2_3 + t194*(t332 - t482) + t195*(t333 + t482) + t325 - 2*t478*t88 + t479*t480 - t480*t483 + pow(t86, 2)*v2sigma2_5 + pow(t88, 2)*v2sigma2_0); + k_grad_rho_z_y = -1.0/4.0*t101*trho_s - 1.0/4.0*t220*tgrad_rho_s_x - 1.0/4.0*t276*tgrad_rho_s_y - 1.0/4.0*t314*tgrad_rho_s_z - 1.0/4.0*t321*t493 - 1.0/4.0*t323*t493 - 1.0/4.0*t324*t493 + (1.0/4.0)*t352*tgrad_rho_x_x + (1.0/4.0)*t391*tgrad_rho_x_y + (1.0/4.0)*t411*tgrad_rho_x_z + (1.0/4.0)*t433*tgrad_rho_y_x + (1.0/4.0)*t459*tgrad_rho_y_y + (1.0/4.0)*t474*tgrad_rho_y_z + (1.0/4.0)*t491*tgrad_rho_z_x + (1.0/4.0)*t502*tgrad_rho_z_z + (1.0/4.0)*tgrad_rho_z_y*(t194*(t371 - t497) + t195*(t372 + t497) + t244*v2sigma2_3 + t325 - 2*t494*t97 + t495*t496 - t496*t498 + pow(t95, 2)*v2sigma2_5 + pow(t97, 2)*v2sigma2_0); + k_grad_rho_z_z = -1.0/4.0*t110*trho_s - 1.0/4.0*t221*tgrad_rho_s_x - 1.0/4.0*t270*tgrad_rho_s_y - 1.0/4.0*t320*tgrad_rho_s_z - 1.0/4.0*t321*t503 - 1.0/4.0*t323*t503 - 1.0/4.0*t324*t503 + (1.0/4.0)*t354*tgrad_rho_x_x + (1.0/4.0)*t387*tgrad_rho_x_y + (1.0/4.0)*t415*tgrad_rho_x_z + (1.0/4.0)*t435*tgrad_rho_y_x + (1.0/4.0)*t457*tgrad_rho_y_y + (1.0/4.0)*t476*tgrad_rho_y_z + (1.0/4.0)*t492*tgrad_rho_z_x + (1.0/4.0)*t502*tgrad_rho_z_y + (1.0/4.0)*tgrad_rho_z_z*(pow(t104, 2)*v2sigma2_5 + pow(t106, 2)*v2sigma2_0 - 2*t106*t499 + t194*(t399 - t505) + t195*(t400 + t505) + t288*v2sigma2_3 + t325 + t500*t504 - t501*t504); +} + +/** Collinear (|m| -> 0) limit of nc_fxc_contract_gga, for points below the host's magnetization cutoff: each magnetization component responds like the spin channel of a collinear perturbation about the spin-symmetric reference; no charge-spin coupling. Depends on the charge fields and the polarized derivatives only. */ +inline void nc_fxc_contract_limit_gga( [[maybe_unused]] double rho_s, double grad_rho_s_x, double grad_rho_s_y, double grad_rho_s_z, [[maybe_unused]] double vrho_0, [[maybe_unused]] double vrho_1, double vsigma_0, double vsigma_1, double vsigma_2, double v2rho2_0, double v2rho2_1, double v2rhosigma_0, double v2rhosigma_1, double v2rhosigma_2, double v2rho2_2, double v2rhosigma_3, double v2rhosigma_4, double v2rhosigma_5, double v2sigma2_0, double v2sigma2_1, double v2sigma2_2, double v2sigma2_3, double v2sigma2_4, double v2sigma2_5, double trho_s, double trho_x, double trho_y, double trho_z, double tgrad_rho_s_x, double tgrad_rho_s_y, double tgrad_rho_s_z, double tgrad_rho_x_x, double tgrad_rho_x_y, double tgrad_rho_x_z, double tgrad_rho_y_x, double tgrad_rho_y_y, double tgrad_rho_y_z, double tgrad_rho_z_x, double tgrad_rho_z_y, double tgrad_rho_z_z, double& k_rho_s, double& k_rho_x, double& k_rho_y, double& k_rho_z, double& k_grad_rho_s_x, double& k_grad_rho_s_y, double& k_grad_rho_s_z, double& k_grad_rho_x_x, double& k_grad_rho_x_y, double& k_grad_rho_x_z, double& k_grad_rho_y_x, double& k_grad_rho_y_y, double& k_grad_rho_y_z, double& k_grad_rho_z_x, double& k_grad_rho_z_y, double& k_grad_rho_z_z ) { + const double t0 = 2*v2rho2_1; + const double t1 = v2rho2_0 + v2rho2_2; + const double t2 = v2rhosigma_0 + v2rhosigma_5; + const double t3 = t2 + v2rhosigma_1 + v2rhosigma_2 + v2rhosigma_3 + v2rhosigma_4; + const double t4 = grad_rho_s_x*t3; + const double t5 = grad_rho_s_y*t3; + const double t6 = grad_rho_s_z*t3; + const double t7 = -t0 + t1; + const double t8 = t2 - v2rhosigma_2 - v2rhosigma_3; + const double t9 = grad_rho_s_x*t8; + const double t10 = grad_rho_s_y*t8; + const double t11 = grad_rho_s_z*t8; + const double t12 = 2*v2sigma2_1; + const double t13 = 2*v2sigma2_2; + const double t14 = 2*v2sigma2_4; + const double t15 = v2sigma2_0 + v2sigma2_5; + const double t16 = t12 + t13 + t14 + t15 + v2sigma2_3; + const double t17 = grad_rho_s_x*t16; + const double t18 = grad_rho_s_y*t17; + const double t19 = grad_rho_s_z*tgrad_rho_s_z; + const double t20 = pow(grad_rho_s_x, 2); + const double t21 = t13*t20; + const double t22 = t20*v2sigma2_0 + t20*v2sigma2_5; + const double t23 = 2*vsigma_1; + const double t24 = 2*vsigma_0 + 2*vsigma_2; + const double t25 = t23 + t24; + const double t26 = grad_rho_s_y*t16; + const double t27 = pow(grad_rho_s_y, 2); + const double t28 = t13*t27; + const double t29 = t27*v2sigma2_0 + t27*v2sigma2_5; + const double t30 = pow(grad_rho_s_z, 2); + const double t31 = t13*t30; + const double t32 = t30*v2sigma2_0 + t30*v2sigma2_5; + const double t33 = -t13 + t15; + const double t34 = grad_rho_s_x*t33; + const double t35 = grad_rho_s_y*t34; + const double t36 = grad_rho_s_z*tgrad_rho_x_z; + const double t37 = -t23 + t24; + const double t38 = -t21 + t22 + t37; + const double t39 = grad_rho_s_y*t33; + const double t40 = -t28 + t29 + t37; + const double t41 = grad_rho_s_z*t34; + const double t42 = grad_rho_s_z*t39; + const double t43 = -t31 + t32 + t37; + k_rho_s = (1.0/4.0)*t4*tgrad_rho_s_x + (1.0/4.0)*t5*tgrad_rho_s_y + (1.0/4.0)*t6*tgrad_rho_s_z + (1.0/4.0)*trho_s*(t0 + t1); + k_rho_x = (1.0/4.0)*t10*tgrad_rho_x_y + (1.0/4.0)*t11*tgrad_rho_x_z + (1.0/4.0)*t7*trho_x + (1.0/4.0)*t9*tgrad_rho_x_x; + k_rho_y = (1.0/4.0)*t10*tgrad_rho_y_y + (1.0/4.0)*t11*tgrad_rho_y_z + (1.0/4.0)*t7*trho_y + (1.0/4.0)*t9*tgrad_rho_y_x; + k_rho_z = (1.0/4.0)*t10*tgrad_rho_z_y + (1.0/4.0)*t11*tgrad_rho_z_z + (1.0/4.0)*t7*trho_z + (1.0/4.0)*t9*tgrad_rho_z_x; + k_grad_rho_s_x = (1.0/4.0)*t17*t19 + (1.0/4.0)*t18*tgrad_rho_s_y + (1.0/4.0)*t4*trho_s + (1.0/4.0)*tgrad_rho_s_x*(t12*t20 + t14*t20 + t20*v2sigma2_3 + t21 + t22 + t25); + k_grad_rho_s_y = (1.0/4.0)*t18*tgrad_rho_s_x + (1.0/4.0)*t19*t26 + (1.0/4.0)*t5*trho_s + (1.0/4.0)*tgrad_rho_s_y*(t12*t27 + t14*t27 + t25 + t27*v2sigma2_3 + t28 + t29); + k_grad_rho_s_z = (1.0/4.0)*grad_rho_s_z*t17*tgrad_rho_s_x + (1.0/4.0)*grad_rho_s_z*t26*tgrad_rho_s_y + (1.0/4.0)*t6*trho_s + (1.0/4.0)*tgrad_rho_s_z*(t12*t30 + t14*t30 + t25 + t30*v2sigma2_3 + t31 + t32); + k_grad_rho_x_x = (1.0/4.0)*t34*t36 + (1.0/4.0)*t35*tgrad_rho_x_y + (1.0/4.0)*t38*tgrad_rho_x_x + (1.0/4.0)*t9*trho_x; + k_grad_rho_x_y = (1.0/4.0)*t10*trho_x + (1.0/4.0)*t35*tgrad_rho_x_x + (1.0/4.0)*t36*t39 + (1.0/4.0)*t40*tgrad_rho_x_y; + k_grad_rho_x_z = (1.0/4.0)*t11*trho_x + (1.0/4.0)*t41*tgrad_rho_x_x + (1.0/4.0)*t42*tgrad_rho_x_y + (1.0/4.0)*t43*tgrad_rho_x_z; + k_grad_rho_y_x = (1.0/4.0)*t35*tgrad_rho_y_y + (1.0/4.0)*t38*tgrad_rho_y_x + (1.0/4.0)*t41*tgrad_rho_y_z + (1.0/4.0)*t9*trho_y; + k_grad_rho_y_y = (1.0/4.0)*t10*trho_y + (1.0/4.0)*t35*tgrad_rho_y_x + (1.0/4.0)*t40*tgrad_rho_y_y + (1.0/4.0)*t42*tgrad_rho_y_z; + k_grad_rho_y_z = (1.0/4.0)*t11*trho_y + (1.0/4.0)*t41*tgrad_rho_y_x + (1.0/4.0)*t42*tgrad_rho_y_y + (1.0/4.0)*t43*tgrad_rho_y_z; + k_grad_rho_z_x = (1.0/4.0)*t35*tgrad_rho_z_y + (1.0/4.0)*t38*tgrad_rho_z_x + (1.0/4.0)*t41*tgrad_rho_z_z + (1.0/4.0)*t9*trho_z; + k_grad_rho_z_y = (1.0/4.0)*t10*trho_z + (1.0/4.0)*t35*tgrad_rho_z_x + (1.0/4.0)*t40*tgrad_rho_z_y + (1.0/4.0)*t42*tgrad_rho_z_z; + k_grad_rho_z_z = (1.0/4.0)*t11*trho_z + (1.0/4.0)*t41*tgrad_rho_z_x + (1.0/4.0)*t42*tgrad_rho_z_y + (1.0/4.0)*t43*tgrad_rho_z_z; +} + +/** Noncollinear GGA exchange-correlation kernel matrix f_xc[X,Y] over the noncollinear field slots; upper triangle, symmetric. */ +inline void nc_fxc_matrix_gga( [[maybe_unused]] double rho_s, double rho_x, double rho_y, double rho_z, double grad_rho_s_x, double grad_rho_s_y, double grad_rho_s_z, double grad_rho_x_x, double grad_rho_x_y, double grad_rho_x_z, double grad_rho_y_x, double grad_rho_y_y, double grad_rho_y_z, double grad_rho_z_x, double grad_rho_z_y, double grad_rho_z_z, double vrho_0, double vrho_1, double vsigma_0, double vsigma_1, double vsigma_2, double v2rho2_0, double v2rho2_1, double v2rhosigma_0, double v2rhosigma_1, double v2rhosigma_2, double v2rho2_2, double v2rhosigma_3, double v2rhosigma_4, double v2rhosigma_5, double v2sigma2_0, double v2sigma2_1, double v2sigma2_2, double v2sigma2_3, double v2sigma2_4, double v2sigma2_5, double f_nabla, double& c_rho_s_rho_s, double& c_rho_s_rho_x, double& c_rho_s_rho_y, double& c_rho_s_rho_z, double& c_rho_s_grad_rho_s_x, double& c_rho_s_grad_rho_s_y, double& c_rho_s_grad_rho_s_z, double& c_rho_s_grad_rho_x_x, double& c_rho_s_grad_rho_x_y, double& c_rho_s_grad_rho_x_z, double& c_rho_s_grad_rho_y_x, double& c_rho_s_grad_rho_y_y, double& c_rho_s_grad_rho_y_z, double& c_rho_s_grad_rho_z_x, double& c_rho_s_grad_rho_z_y, double& c_rho_s_grad_rho_z_z, double& c_rho_x_rho_x, double& c_rho_x_rho_y, double& c_rho_x_rho_z, double& c_rho_x_grad_rho_s_x, double& c_rho_x_grad_rho_s_y, double& c_rho_x_grad_rho_s_z, double& c_rho_x_grad_rho_x_x, double& c_rho_x_grad_rho_x_y, double& c_rho_x_grad_rho_x_z, double& c_rho_x_grad_rho_y_x, double& c_rho_x_grad_rho_y_y, double& c_rho_x_grad_rho_y_z, double& c_rho_x_grad_rho_z_x, double& c_rho_x_grad_rho_z_y, double& c_rho_x_grad_rho_z_z, double& c_rho_y_rho_y, double& c_rho_y_rho_z, double& c_rho_y_grad_rho_s_x, double& c_rho_y_grad_rho_s_y, double& c_rho_y_grad_rho_s_z, double& c_rho_y_grad_rho_x_x, double& c_rho_y_grad_rho_x_y, double& c_rho_y_grad_rho_x_z, double& c_rho_y_grad_rho_y_x, double& c_rho_y_grad_rho_y_y, double& c_rho_y_grad_rho_y_z, double& c_rho_y_grad_rho_z_x, double& c_rho_y_grad_rho_z_y, double& c_rho_y_grad_rho_z_z, double& c_rho_z_rho_z, double& c_rho_z_grad_rho_s_x, double& c_rho_z_grad_rho_s_y, double& c_rho_z_grad_rho_s_z, double& c_rho_z_grad_rho_x_x, double& c_rho_z_grad_rho_x_y, double& c_rho_z_grad_rho_x_z, double& c_rho_z_grad_rho_y_x, double& c_rho_z_grad_rho_y_y, double& c_rho_z_grad_rho_y_z, double& c_rho_z_grad_rho_z_x, double& c_rho_z_grad_rho_z_y, double& c_rho_z_grad_rho_z_z, double& c_grad_rho_s_x_grad_rho_s_x, double& c_grad_rho_s_x_grad_rho_s_y, double& c_grad_rho_s_x_grad_rho_s_z, double& c_grad_rho_s_x_grad_rho_x_x, double& c_grad_rho_s_x_grad_rho_x_y, double& c_grad_rho_s_x_grad_rho_x_z, double& c_grad_rho_s_x_grad_rho_y_x, double& c_grad_rho_s_x_grad_rho_y_y, double& c_grad_rho_s_x_grad_rho_y_z, double& c_grad_rho_s_x_grad_rho_z_x, double& c_grad_rho_s_x_grad_rho_z_y, double& c_grad_rho_s_x_grad_rho_z_z, double& c_grad_rho_s_y_grad_rho_s_y, double& c_grad_rho_s_y_grad_rho_s_z, double& c_grad_rho_s_y_grad_rho_x_x, double& c_grad_rho_s_y_grad_rho_x_y, double& c_grad_rho_s_y_grad_rho_x_z, double& c_grad_rho_s_y_grad_rho_y_x, double& c_grad_rho_s_y_grad_rho_y_y, double& c_grad_rho_s_y_grad_rho_y_z, double& c_grad_rho_s_y_grad_rho_z_x, double& c_grad_rho_s_y_grad_rho_z_y, double& c_grad_rho_s_y_grad_rho_z_z, double& c_grad_rho_s_z_grad_rho_s_z, double& c_grad_rho_s_z_grad_rho_x_x, double& c_grad_rho_s_z_grad_rho_x_y, double& c_grad_rho_s_z_grad_rho_x_z, double& c_grad_rho_s_z_grad_rho_y_x, double& c_grad_rho_s_z_grad_rho_y_y, double& c_grad_rho_s_z_grad_rho_y_z, double& c_grad_rho_s_z_grad_rho_z_x, double& c_grad_rho_s_z_grad_rho_z_y, double& c_grad_rho_s_z_grad_rho_z_z, double& c_grad_rho_x_x_grad_rho_x_x, double& c_grad_rho_x_x_grad_rho_x_y, double& c_grad_rho_x_x_grad_rho_x_z, double& c_grad_rho_x_x_grad_rho_y_x, double& c_grad_rho_x_x_grad_rho_y_y, double& c_grad_rho_x_x_grad_rho_y_z, double& c_grad_rho_x_x_grad_rho_z_x, double& c_grad_rho_x_x_grad_rho_z_y, double& c_grad_rho_x_x_grad_rho_z_z, double& c_grad_rho_x_y_grad_rho_x_y, double& c_grad_rho_x_y_grad_rho_x_z, double& c_grad_rho_x_y_grad_rho_y_x, double& c_grad_rho_x_y_grad_rho_y_y, double& c_grad_rho_x_y_grad_rho_y_z, double& c_grad_rho_x_y_grad_rho_z_x, double& c_grad_rho_x_y_grad_rho_z_y, double& c_grad_rho_x_y_grad_rho_z_z, double& c_grad_rho_x_z_grad_rho_x_z, double& c_grad_rho_x_z_grad_rho_y_x, double& c_grad_rho_x_z_grad_rho_y_y, double& c_grad_rho_x_z_grad_rho_y_z, double& c_grad_rho_x_z_grad_rho_z_x, double& c_grad_rho_x_z_grad_rho_z_y, double& c_grad_rho_x_z_grad_rho_z_z, double& c_grad_rho_y_x_grad_rho_y_x, double& c_grad_rho_y_x_grad_rho_y_y, double& c_grad_rho_y_x_grad_rho_y_z, double& c_grad_rho_y_x_grad_rho_z_x, double& c_grad_rho_y_x_grad_rho_z_y, double& c_grad_rho_y_x_grad_rho_z_z, double& c_grad_rho_y_y_grad_rho_y_y, double& c_grad_rho_y_y_grad_rho_y_z, double& c_grad_rho_y_y_grad_rho_z_x, double& c_grad_rho_y_y_grad_rho_z_y, double& c_grad_rho_y_y_grad_rho_z_z, double& c_grad_rho_y_z_grad_rho_y_z, double& c_grad_rho_y_z_grad_rho_z_x, double& c_grad_rho_y_z_grad_rho_z_y, double& c_grad_rho_y_z_grad_rho_z_z, double& c_grad_rho_z_x_grad_rho_z_x, double& c_grad_rho_z_x_grad_rho_z_y, double& c_grad_rho_z_x_grad_rho_z_z, double& c_grad_rho_z_y_grad_rho_z_y, double& c_grad_rho_z_y_grad_rho_z_z, double& c_grad_rho_z_z_grad_rho_z_z ) { + const double t0 = 2*v2rho2_1; + const double t1 = (1.0/4.0)*rho_x; + const double t2 = pow(rho_x, 2); + const double t3 = pow(rho_y, 2); + const double t4 = pow(rho_z, 2); + const double t5 = t2 + t3 + t4; + const double t6 = pow(t5, -1.0/2.0); + const double t7 = t6*(v2rho2_0 - v2rho2_2); + const double t8 = (1.0/4.0)*t7; + const double t9 = grad_rho_s_x*v2rhosigma_4; + const double t10 = grad_rho_s_x*grad_rho_x_x; + const double t11 = grad_rho_s_y*grad_rho_x_y; + const double t12 = grad_rho_s_z*grad_rho_x_z; + const double t13 = t11 + t12; + const double t14 = t10 + t13; + const double t15 = grad_rho_s_x*grad_rho_y_x; + const double t16 = grad_rho_s_y*grad_rho_y_y; + const double t17 = grad_rho_s_z*grad_rho_y_z; + const double t18 = t16 + t17; + const double t19 = t15 + t18; + const double t20 = grad_rho_s_x*grad_rho_z_x; + const double t21 = grad_rho_s_y*grad_rho_z_y; + const double t22 = grad_rho_s_z*grad_rho_z_z; + const double t23 = t21 + t22; + const double t24 = t20 + t23; + const double t25 = grad_rho_x_x*t14 + grad_rho_y_x*t19 + grad_rho_z_x*t24; + const double t26 = pow(t14, 2); + const double t27 = pow(t19, 2); + const double t28 = pow(t24, 2); + const double t29 = t26 + t27 + t28; + const double t30 = pow(t29, -1.0/2.0); + const double t31 = f_nabla*t30; + const double t32 = t25*t31; + const double t33 = grad_rho_s_x + t32; + const double t34 = t33*v2rhosigma_3; + const double t35 = -grad_rho_s_x + t32; + const double t36 = t35*v2rhosigma_5; + const double t37 = grad_rho_s_x*v2rhosigma_1 + t33*v2rhosigma_0 - t35*v2rhosigma_2; + const double t38 = grad_rho_s_y*v2rhosigma_4; + const double t39 = grad_rho_x_y*t14 + grad_rho_y_y*t19 + grad_rho_z_y*t24; + const double t40 = t31*t39; + const double t41 = grad_rho_s_y + t40; + const double t42 = t41*v2rhosigma_3; + const double t43 = -grad_rho_s_y + t40; + const double t44 = t43*v2rhosigma_5; + const double t45 = grad_rho_s_y*v2rhosigma_1 + t41*v2rhosigma_0 - t43*v2rhosigma_2; + const double t46 = grad_rho_s_z*v2rhosigma_4; + const double t47 = grad_rho_x_z*t14 + grad_rho_y_z*t19 + grad_rho_z_z*t24; + const double t48 = t31*t47; + const double t49 = grad_rho_s_z + t48; + const double t50 = t49*v2rhosigma_3; + const double t51 = -grad_rho_s_z + t48; + const double t52 = t51*v2rhosigma_5; + const double t53 = grad_rho_s_z*v2rhosigma_1 + t49*v2rhosigma_0 - t51*v2rhosigma_2; + const double t54 = grad_rho_x_x*v2rhosigma_4; + const double t55 = t14*t31; + const double t56 = grad_rho_s_x*t55; + const double t57 = -grad_rho_x_x + t56; + const double t58 = t57*v2rhosigma_5; + const double t59 = grad_rho_x_x + t56; + const double t60 = t59*v2rhosigma_3; + const double t61 = grad_rho_x_x*v2rhosigma_1 + t57*v2rhosigma_2 - t59*v2rhosigma_0; + const double t62 = grad_rho_x_y*v2rhosigma_4; + const double t63 = grad_rho_s_y*t55; + const double t64 = -grad_rho_x_y + t63; + const double t65 = t64*v2rhosigma_5; + const double t66 = grad_rho_x_y + t63; + const double t67 = t66*v2rhosigma_3; + const double t68 = grad_rho_x_y*v2rhosigma_1 + t64*v2rhosigma_2 - t66*v2rhosigma_0; + const double t69 = grad_rho_x_z*v2rhosigma_4; + const double t70 = grad_rho_s_z*t55; + const double t71 = -grad_rho_x_z + t70; + const double t72 = t71*v2rhosigma_5; + const double t73 = grad_rho_x_z + t70; + const double t74 = t73*v2rhosigma_3; + const double t75 = grad_rho_x_z*v2rhosigma_1 + t71*v2rhosigma_2 - t73*v2rhosigma_0; + const double t76 = grad_rho_y_x*v2rhosigma_4; + const double t77 = t19*t31; + const double t78 = grad_rho_s_x*t77; + const double t79 = -grad_rho_y_x + t78; + const double t80 = t79*v2rhosigma_5; + const double t81 = grad_rho_y_x + t78; + const double t82 = t81*v2rhosigma_3; + const double t83 = grad_rho_y_x*v2rhosigma_1 + t79*v2rhosigma_2 - t81*v2rhosigma_0; + const double t84 = grad_rho_y_y*v2rhosigma_4; + const double t85 = grad_rho_s_y*t77; + const double t86 = -grad_rho_y_y + t85; + const double t87 = t86*v2rhosigma_5; + const double t88 = grad_rho_y_y + t85; + const double t89 = t88*v2rhosigma_3; + const double t90 = grad_rho_y_y*v2rhosigma_1 + t86*v2rhosigma_2 - t88*v2rhosigma_0; + const double t91 = grad_rho_y_z*v2rhosigma_4; + const double t92 = grad_rho_s_z*t77; + const double t93 = -grad_rho_y_z + t92; + const double t94 = t93*v2rhosigma_5; + const double t95 = grad_rho_y_z + t92; + const double t96 = t95*v2rhosigma_3; + const double t97 = grad_rho_y_z*v2rhosigma_1 + t93*v2rhosigma_2 - t95*v2rhosigma_0; + const double t98 = grad_rho_z_x*v2rhosigma_4; + const double t99 = t24*t31; + const double t100 = grad_rho_s_x*t99; + const double t101 = -grad_rho_z_x + t100; + const double t102 = t101*v2rhosigma_5; + const double t103 = grad_rho_z_x + t100; + const double t104 = t103*v2rhosigma_3; + const double t105 = grad_rho_z_x*v2rhosigma_1 + t101*v2rhosigma_2 - t103*v2rhosigma_0; + const double t106 = grad_rho_z_y*v2rhosigma_4; + const double t107 = grad_rho_s_y*t99; + const double t108 = -grad_rho_z_y + t107; + const double t109 = t108*v2rhosigma_5; + const double t110 = grad_rho_z_y + t107; + const double t111 = t110*v2rhosigma_3; + const double t112 = grad_rho_z_y*v2rhosigma_1 + t108*v2rhosigma_2 - t110*v2rhosigma_0; + const double t113 = grad_rho_z_z*v2rhosigma_4; + const double t114 = grad_rho_s_z*t99; + const double t115 = -grad_rho_z_z + t114; + const double t116 = t115*v2rhosigma_5; + const double t117 = grad_rho_z_z + t114; + const double t118 = t117*v2rhosigma_3; + const double t119 = grad_rho_z_z*v2rhosigma_1 + t115*v2rhosigma_2 - t117*v2rhosigma_0; + const double t120 = 1.0/t5; + const double t121 = t120*t2; + const double t122 = 2*t6; + const double t123 = t122*(t121 - 1); + const double t124 = t120*v2rho2_0; + const double t125 = t120*v2rho2_2; + const double t126 = 2/pow(t5, 3.0/2.0); + const double t127 = -t0*t120 + t124 + t125 - t126*vrho_0 + t126*vrho_1; + const double t128 = t1*t127; + const double t129 = -t34 + t36 + t37 - t9; + const double t130 = t1*t6; + const double t131 = -t38 - t42 + t44 + t45; + const double t132 = -t46 - t50 + t52 + t53; + const double t133 = t54 + t58 - t60 - t61; + const double t134 = t62 + t65 - t67 - t68; + const double t135 = t69 + t72 - t74 - t75; + const double t136 = t76 + t80 - t82 - t83; + const double t137 = t84 + t87 - t89 - t90; + const double t138 = t91 + t94 - t96 - t97; + const double t139 = t102 - t104 - t105 + t98; + const double t140 = t106 + t109 - t111 - t112; + const double t141 = t113 + t116 - t118 - t119; + const double t142 = t120*t3; + const double t143 = t122*(t142 - 1); + const double t144 = (1.0/4.0)*rho_y; + const double t145 = t144*t6; + const double t146 = t120*t4; + const double t147 = t122*(t146 - 1); + const double t148 = (1.0/4.0)*rho_z*t6; + const double t149 = 2*vsigma_1; + const double t150 = pow(grad_rho_s_x, 2); + const double t151 = grad_rho_s_x*v2sigma2_1; + const double t152 = grad_rho_s_x*v2sigma2_4; + const double t153 = 2*t35; + const double t154 = pow(grad_rho_x_x, 2); + const double t155 = pow(grad_rho_y_x, 2); + const double t156 = pow(grad_rho_z_x, 2); + const double t157 = t31*(t154 + t155 + t156); + const double t158 = f_nabla/pow(t29, 3.0/2.0); + const double t159 = t158*pow(t25, 2); + const double t160 = 2*vsigma_0; + const double t161 = 2*vsigma_2; + const double t162 = t33*v2sigma2_2; + const double t163 = grad_rho_s_x*v2sigma2_3; + const double t164 = t33*v2sigma2_1; + const double t165 = t35*v2sigma2_4; + const double t166 = t33*v2sigma2_0; + const double t167 = t35*v2sigma2_2; + const double t168 = t35*v2sigma2_5; + const double t169 = grad_rho_x_x*grad_rho_x_y; + const double t170 = grad_rho_y_x*grad_rho_y_y; + const double t171 = grad_rho_z_x*grad_rho_z_y; + const double t172 = 1.0/t29; + const double t173 = t172*t25; + const double t174 = t31*(t169 + t170 + t171 - t173*t39); + const double t175 = grad_rho_x_x*grad_rho_x_z; + const double t176 = grad_rho_y_x*grad_rho_y_z; + const double t177 = grad_rho_z_x*grad_rho_z_z; + const double t178 = t31*(-t173*t47 + t175 + t176 + t177); + const double t179 = t14*t173; + const double t180 = -grad_rho_s_x*t179 + 2*t10 + t13; + const double t181 = grad_rho_x_x - t179; + const double t182 = t181*t31; + const double t183 = grad_rho_s_y*t161; + const double t184 = grad_rho_s_z*t161; + const double t185 = t173*t19; + const double t186 = -grad_rho_s_x*t185 + 2*t15 + t18; + const double t187 = grad_rho_y_x - t185; + const double t188 = t187*t31; + const double t189 = t173*t24; + const double t190 = -grad_rho_s_x*t189 + 2*t20 + t23; + const double t191 = grad_rho_z_x - t189; + const double t192 = t191*t31; + const double t193 = pow(grad_rho_s_y, 2); + const double t194 = grad_rho_s_y*v2sigma2_1; + const double t195 = grad_rho_s_y*v2sigma2_4; + const double t196 = 2*t43; + const double t197 = pow(grad_rho_x_y, 2); + const double t198 = pow(grad_rho_y_y, 2); + const double t199 = pow(grad_rho_z_y, 2); + const double t200 = t31*(t197 + t198 + t199); + const double t201 = t158*pow(t39, 2); + const double t202 = t41*v2sigma2_2; + const double t203 = grad_rho_s_y*v2sigma2_3; + const double t204 = t41*v2sigma2_1; + const double t205 = t43*v2sigma2_2; + const double t206 = grad_rho_x_y*grad_rho_x_z; + const double t207 = grad_rho_y_y*grad_rho_y_z; + const double t208 = grad_rho_z_y*grad_rho_z_z; + const double t209 = t172*t39; + const double t210 = t31*(t206 + t207 + t208 - t209*t47); + const double t211 = t14*t209; + const double t212 = grad_rho_x_y - t211; + const double t213 = t212*t31; + const double t214 = -grad_rho_s_y*t211 + t10 + 2*t11 + t12; + const double t215 = t19*t209; + const double t216 = grad_rho_y_y - t215; + const double t217 = t216*t31; + const double t218 = -grad_rho_s_y*t215 + t15 + 2*t16 + t17; + const double t219 = t209*t24; + const double t220 = grad_rho_z_y - t219; + const double t221 = t220*t31; + const double t222 = -grad_rho_s_y*t219 + t20 + 2*t21 + t22; + const double t223 = pow(grad_rho_s_z, 2); + const double t224 = pow(grad_rho_x_z, 2); + const double t225 = pow(grad_rho_y_z, 2); + const double t226 = pow(grad_rho_z_z, 2); + const double t227 = t31*(t224 + t225 + t226); + const double t228 = t158*pow(t47, 2); + const double t229 = t49*v2sigma2_2; + const double t230 = grad_rho_s_z*v2sigma2_3; + const double t231 = grad_rho_s_z*v2sigma2_4; + const double t232 = t49*v2sigma2_1; + const double t233 = t172*t47; + const double t234 = t14*t233; + const double t235 = grad_rho_x_z - t234; + const double t236 = t235*t31; + const double t237 = t51*v2sigma2_2; + const double t238 = -grad_rho_s_z*t234 + t10 + t11 + 2*t12; + const double t239 = t19*t233; + const double t240 = grad_rho_y_z - t239; + const double t241 = t240*t31; + const double t242 = -grad_rho_s_z*t239 + t15 + t16 + 2*t17; + const double t243 = t233*t24; + const double t244 = grad_rho_z_z - t243; + const double t245 = t244*t31; + const double t246 = -grad_rho_s_z*t243 + t20 + t21 + 2*t22; + const double t247 = -t149; + const double t248 = grad_rho_x_x*v2sigma2_1; + const double t249 = grad_rho_x_x*v2sigma2_4; + const double t250 = 2*t57; + const double t251 = t158*t26; + const double t252 = t150*t251; + const double t253 = t150*t31; + const double t254 = t253 + 1; + const double t255 = 1 - t253; + const double t256 = t59*v2sigma2_2; + const double t257 = t59*v2sigma2_1; + const double t258 = t57*v2sigma2_4; + const double t259 = t31*(t172*t26 - 1); + const double t260 = grad_rho_s_x*grad_rho_s_y*t160; + const double t261 = grad_rho_s_x*t259; + const double t262 = t59*v2sigma2_0; + const double t263 = t57*v2sigma2_2; + const double t264 = t57*v2sigma2_5; + const double t265 = grad_rho_s_z*t160; + const double t266 = grad_rho_x_x*v2sigma2_3; + const double t267 = t14*t158; + const double t268 = t150*t267; + const double t269 = t160*t19; + const double t270 = t161*t268; + const double t271 = t19*t267; + const double t272 = grad_rho_s_x*t271; + const double t273 = t183*t272 - t260*t271; + const double t274 = t184*t272 - t265*t272; + const double t275 = t160*t24; + const double t276 = t24*t267; + const double t277 = grad_rho_s_x*t276; + const double t278 = t183*t277 - t260*t276; + const double t279 = t184*t277 - t265*t277; + const double t280 = grad_rho_x_y*v2sigma2_1; + const double t281 = grad_rho_x_y*v2sigma2_4; + const double t282 = 2*t64; + const double t283 = t193*t251; + const double t284 = t193*t31; + const double t285 = t284 + 1; + const double t286 = 1 - t284; + const double t287 = t66*v2sigma2_2; + const double t288 = t66*v2sigma2_1; + const double t289 = t64*v2sigma2_4; + const double t290 = grad_rho_s_y*t265; + const double t291 = grad_rho_s_z*t183; + const double t292 = t66*v2sigma2_0; + const double t293 = t64*v2sigma2_2; + const double t294 = t64*v2sigma2_5; + const double t295 = grad_rho_x_y*v2sigma2_3; + const double t296 = t193*t267; + const double t297 = t161*t296; + const double t298 = -t271*t290 + t271*t291; + const double t299 = -t276*t290 + t276*t291; + const double t300 = grad_rho_x_z*v2sigma2_1; + const double t301 = grad_rho_x_z*v2sigma2_4; + const double t302 = 2*t71; + const double t303 = t223*t251; + const double t304 = t223*t31; + const double t305 = t304 + 1; + const double t306 = 1 - t304; + const double t307 = t73*v2sigma2_2; + const double t308 = grad_rho_x_z*v2sigma2_3; + const double t309 = t71*v2sigma2_4; + const double t310 = t73*v2sigma2_0; + const double t311 = t71*v2sigma2_5; + const double t312 = t73*v2sigma2_1; + const double t313 = t71*v2sigma2_2; + const double t314 = t223*t267; + const double t315 = t161*t314; + const double t316 = grad_rho_y_x*v2sigma2_1; + const double t317 = grad_rho_y_x*v2sigma2_4; + const double t318 = 2*t79; + const double t319 = t158*t27; + const double t320 = t150*t319; + const double t321 = t81*v2sigma2_2; + const double t322 = t81*v2sigma2_1; + const double t323 = t79*v2sigma2_4; + const double t324 = t31*(t172*t27 - 1); + const double t325 = grad_rho_s_x*t324; + const double t326 = t81*v2sigma2_0; + const double t327 = t79*v2sigma2_2; + const double t328 = t79*v2sigma2_5; + const double t329 = grad_rho_y_x*v2sigma2_3; + const double t330 = t158*t24; + const double t331 = t150*t330; + const double t332 = t161*t19; + const double t333 = t19*t330; + const double t334 = grad_rho_s_x*t333; + const double t335 = t183*t334 - t260*t333; + const double t336 = t184*t334 - t265*t334; + const double t337 = grad_rho_y_y*v2sigma2_1; + const double t338 = grad_rho_y_y*v2sigma2_4; + const double t339 = 2*t86; + const double t340 = t193*t319; + const double t341 = t88*v2sigma2_2; + const double t342 = t88*v2sigma2_1; + const double t343 = t86*v2sigma2_4; + const double t344 = t88*v2sigma2_0; + const double t345 = t86*v2sigma2_2; + const double t346 = t86*v2sigma2_5; + const double t347 = grad_rho_y_y*v2sigma2_3; + const double t348 = t193*t330; + const double t349 = -t290*t333 + t291*t333; + const double t350 = grad_rho_y_z*v2sigma2_1; + const double t351 = grad_rho_y_z*v2sigma2_4; + const double t352 = 2*t93; + const double t353 = t223*t319; + const double t354 = t95*v2sigma2_2; + const double t355 = grad_rho_y_z*v2sigma2_3; + const double t356 = t93*v2sigma2_4; + const double t357 = t95*v2sigma2_0; + const double t358 = t93*v2sigma2_5; + const double t359 = t95*v2sigma2_1; + const double t360 = t93*v2sigma2_2; + const double t361 = t223*t330; + const double t362 = grad_rho_z_x*v2sigma2_1; + const double t363 = grad_rho_z_x*v2sigma2_4; + const double t364 = 2*t101; + const double t365 = t158*t28; + const double t366 = t150*t365; + const double t367 = t103*v2sigma2_2; + const double t368 = t103*v2sigma2_1; + const double t369 = t101*v2sigma2_4; + const double t370 = t31*(t172*t28 - 1); + const double t371 = grad_rho_s_x*t370; + const double t372 = t103*v2sigma2_0; + const double t373 = t101*v2sigma2_2; + const double t374 = t101*v2sigma2_5; + const double t375 = grad_rho_z_y*v2sigma2_1; + const double t376 = grad_rho_z_y*v2sigma2_4; + const double t377 = 2*t108; + const double t378 = t193*t365; + const double t379 = t110*v2sigma2_2; + const double t380 = grad_rho_z_z*v2sigma2_1; + const double t381 = grad_rho_z_z*v2sigma2_4; + const double t382 = t117*v2sigma2_2; + const double t383 = 2*t115; + const double t384 = t223*t365; + c_rho_s_rho_s = (1.0/4.0)*t0 + (1.0/4.0)*v2rho2_0 + (1.0/4.0)*v2rho2_2; + c_rho_s_rho_x = t1*t7; + c_rho_s_rho_y = rho_y*t8; + c_rho_s_rho_z = rho_z*t8; + c_rho_s_grad_rho_s_x = (1.0/4.0)*t34 - 1.0/4.0*t36 + (1.0/4.0)*t37 + (1.0/4.0)*t9; + c_rho_s_grad_rho_s_y = (1.0/4.0)*t38 + (1.0/4.0)*t42 - 1.0/4.0*t44 + (1.0/4.0)*t45; + c_rho_s_grad_rho_s_z = (1.0/4.0)*t46 + (1.0/4.0)*t50 - 1.0/4.0*t52 + (1.0/4.0)*t53; + c_rho_s_grad_rho_x_x = -1.0/4.0*t54 - 1.0/4.0*t58 + (1.0/4.0)*t60 - 1.0/4.0*t61; + c_rho_s_grad_rho_x_y = -1.0/4.0*t62 - 1.0/4.0*t65 + (1.0/4.0)*t67 - 1.0/4.0*t68; + c_rho_s_grad_rho_x_z = -1.0/4.0*t69 - 1.0/4.0*t72 + (1.0/4.0)*t74 - 1.0/4.0*t75; + c_rho_s_grad_rho_y_x = -1.0/4.0*t76 - 1.0/4.0*t80 + (1.0/4.0)*t82 - 1.0/4.0*t83; + c_rho_s_grad_rho_y_y = -1.0/4.0*t84 - 1.0/4.0*t87 + (1.0/4.0)*t89 - 1.0/4.0*t90; + c_rho_s_grad_rho_y_z = -1.0/4.0*t91 - 1.0/4.0*t94 + (1.0/4.0)*t96 - 1.0/4.0*t97; + c_rho_s_grad_rho_z_x = -1.0/4.0*t102 + (1.0/4.0)*t104 - 1.0/4.0*t105 - 1.0/4.0*t98; + c_rho_s_grad_rho_z_y = -1.0/4.0*t106 - 1.0/4.0*t109 + (1.0/4.0)*t111 - 1.0/4.0*t112; + c_rho_s_grad_rho_z_z = -1.0/4.0*t113 - 1.0/4.0*t116 + (1.0/4.0)*t118 - 1.0/4.0*t119; + c_rho_x_rho_x = -1.0/4.0*t0*t121 + (1.0/4.0)*t121*v2rho2_0 + (1.0/4.0)*t121*v2rho2_2 - 1.0/4.0*t123*vrho_0 + (1.0/4.0)*t123*vrho_1; + c_rho_x_rho_y = rho_y*t128; + c_rho_x_rho_z = rho_z*t128; + c_rho_x_grad_rho_s_x = t129*t130; + c_rho_x_grad_rho_s_y = t130*t131; + c_rho_x_grad_rho_s_z = t130*t132; + c_rho_x_grad_rho_x_x = t130*t133; + c_rho_x_grad_rho_x_y = t130*t134; + c_rho_x_grad_rho_x_z = t130*t135; + c_rho_x_grad_rho_y_x = t130*t136; + c_rho_x_grad_rho_y_y = t130*t137; + c_rho_x_grad_rho_y_z = t130*t138; + c_rho_x_grad_rho_z_x = t130*t139; + c_rho_x_grad_rho_z_y = t130*t140; + c_rho_x_grad_rho_z_z = t130*t141; + c_rho_y_rho_y = -1.0/4.0*t0*t142 + (1.0/4.0)*t124*t3 + (1.0/4.0)*t125*t3 - 1.0/4.0*t143*vrho_0 + (1.0/4.0)*t143*vrho_1; + c_rho_y_rho_z = rho_z*t127*t144; + c_rho_y_grad_rho_s_x = t129*t145; + c_rho_y_grad_rho_s_y = t131*t145; + c_rho_y_grad_rho_s_z = t132*t145; + c_rho_y_grad_rho_x_x = t133*t145; + c_rho_y_grad_rho_x_y = t134*t145; + c_rho_y_grad_rho_x_z = t135*t145; + c_rho_y_grad_rho_y_x = t136*t145; + c_rho_y_grad_rho_y_y = t137*t145; + c_rho_y_grad_rho_y_z = t138*t145; + c_rho_y_grad_rho_z_x = t139*t145; + c_rho_y_grad_rho_z_y = t140*t145; + c_rho_y_grad_rho_z_z = t141*t145; + c_rho_z_rho_z = -1.0/4.0*t0*t146 + (1.0/4.0)*t124*t4 + (1.0/4.0)*t125*t4 - 1.0/4.0*t147*vrho_0 + (1.0/4.0)*t147*vrho_1; + c_rho_z_grad_rho_s_x = t129*t148; + c_rho_z_grad_rho_s_y = t131*t148; + c_rho_z_grad_rho_s_z = t132*t148; + c_rho_z_grad_rho_x_x = t133*t148; + c_rho_z_grad_rho_x_y = t134*t148; + c_rho_z_grad_rho_x_z = t135*t148; + c_rho_z_grad_rho_y_x = t136*t148; + c_rho_z_grad_rho_y_y = t137*t148; + c_rho_z_grad_rho_y_z = t138*t148; + c_rho_z_grad_rho_z_x = t139*t148; + c_rho_z_grad_rho_z_y = t140*t148; + c_rho_z_grad_rho_z_z = t141*t148; + c_grad_rho_s_x_grad_rho_s_x = (1.0/4.0)*t149 + (1.0/4.0)*t150*v2sigma2_3 + (1.0/2.0)*t151*t33 - 1.0/4.0*t152*t153 - 1.0/4.0*t153*t162 + (1.0/4.0)*t160*(t157 - t159 + 1) + (1.0/4.0)*t161*(-t157 + t159 + 1) + (1.0/4.0)*pow(t33, 2)*v2sigma2_0 + (1.0/4.0)*pow(t35, 2)*v2sigma2_5; + c_grad_rho_s_x_grad_rho_s_y = (1.0/4.0)*grad_rho_s_y*t163 + (1.0/4.0)*grad_rho_s_y*t164 - 1.0/4.0*grad_rho_s_y*t165 + (1.0/4.0)*t151*t41 - 1.0/4.0*t152*t43 + (1.0/4.0)*t160*t174 - 1.0/4.0*t161*t174 - 1.0/4.0*t162*t43 + (1.0/4.0)*t166*t41 - 1.0/4.0*t167*t41 + (1.0/4.0)*t168*t43; + c_grad_rho_s_x_grad_rho_s_z = (1.0/4.0)*grad_rho_s_z*t163 + (1.0/4.0)*grad_rho_s_z*t164 - 1.0/4.0*grad_rho_s_z*t165 + (1.0/4.0)*t151*t49 - 1.0/4.0*t152*t51 + (1.0/4.0)*t160*t178 - 1.0/4.0*t161*t178 - 1.0/4.0*t162*t51 + (1.0/4.0)*t166*t49 - 1.0/4.0*t167*t49 + (1.0/4.0)*t168*t51; + c_grad_rho_s_x_grad_rho_x_x = (1.0/2.0)*f_nabla*t180*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_x*t59*v2sigma2_1 - 1.0/4.0*grad_rho_x_x*t164 + (1.0/4.0)*grad_rho_x_x*t35*v2sigma2_4 - 1.0/4.0*t10*v2sigma2_3 - 1.0/4.0*t152*t57 - 1.0/4.0*t161*t180*t31 - 1.0/4.0*t162*t57 - 1.0/4.0*t167*t59 + (1.0/4.0)*t33*t59*v2sigma2_0 + (1.0/4.0)*t35*t57*v2sigma2_5; + c_grad_rho_s_x_grad_rho_x_y = (1.0/2.0)*f_nabla*grad_rho_s_y*t181*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_x*t66*v2sigma2_1 - 1.0/4.0*grad_rho_x_y*t163 - 1.0/4.0*grad_rho_x_y*t164 + (1.0/4.0)*grad_rho_x_y*t35*v2sigma2_4 - 1.0/4.0*t152*t64 - 1.0/4.0*t162*t64 - 1.0/4.0*t167*t66 - 1.0/4.0*t182*t183 + (1.0/4.0)*t33*t66*v2sigma2_0 + (1.0/4.0)*t35*t64*v2sigma2_5; + c_grad_rho_s_x_grad_rho_x_z = (1.0/2.0)*f_nabla*grad_rho_s_z*t181*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_x*t73*v2sigma2_1 - 1.0/4.0*grad_rho_x_z*t163 - 1.0/4.0*grad_rho_x_z*t164 + (1.0/4.0)*grad_rho_x_z*t35*v2sigma2_4 - 1.0/4.0*t152*t71 - 1.0/4.0*t162*t71 - 1.0/4.0*t167*t73 - 1.0/4.0*t182*t184 + (1.0/4.0)*t33*t73*v2sigma2_0 + (1.0/4.0)*t35*t71*v2sigma2_5; + c_grad_rho_s_x_grad_rho_y_x = (1.0/2.0)*f_nabla*t186*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_x*t81*v2sigma2_1 - 1.0/4.0*grad_rho_y_x*t164 + (1.0/4.0)*grad_rho_y_x*t35*v2sigma2_4 - 1.0/4.0*t15*v2sigma2_3 - 1.0/4.0*t152*t79 - 1.0/4.0*t161*t186*t31 - 1.0/4.0*t162*t79 - 1.0/4.0*t167*t81 + (1.0/4.0)*t33*t81*v2sigma2_0 + (1.0/4.0)*t35*t79*v2sigma2_5; + c_grad_rho_s_x_grad_rho_y_y = (1.0/2.0)*f_nabla*grad_rho_s_y*t187*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_x*t88*v2sigma2_1 - 1.0/4.0*grad_rho_y_y*t163 - 1.0/4.0*grad_rho_y_y*t164 + (1.0/4.0)*grad_rho_y_y*t35*v2sigma2_4 - 1.0/4.0*t152*t86 - 1.0/4.0*t162*t86 - 1.0/4.0*t167*t88 - 1.0/4.0*t183*t188 + (1.0/4.0)*t33*t88*v2sigma2_0 + (1.0/4.0)*t35*t86*v2sigma2_5; + c_grad_rho_s_x_grad_rho_y_z = (1.0/2.0)*f_nabla*grad_rho_s_z*t187*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_x*t95*v2sigma2_1 - 1.0/4.0*grad_rho_y_z*t163 - 1.0/4.0*grad_rho_y_z*t164 + (1.0/4.0)*grad_rho_y_z*t35*v2sigma2_4 - 1.0/4.0*t152*t93 - 1.0/4.0*t162*t93 - 1.0/4.0*t167*t95 - 1.0/4.0*t184*t188 + (1.0/4.0)*t33*t95*v2sigma2_0 + (1.0/4.0)*t35*t93*v2sigma2_5; + c_grad_rho_s_x_grad_rho_z_x = (1.0/2.0)*f_nabla*t190*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_x*t103*v2sigma2_1 - 1.0/4.0*grad_rho_z_x*t164 + (1.0/4.0)*grad_rho_z_x*t35*v2sigma2_4 - 1.0/4.0*t101*t152 - 1.0/4.0*t101*t162 + (1.0/4.0)*t101*t35*v2sigma2_5 - 1.0/4.0*t103*t167 + (1.0/4.0)*t103*t33*v2sigma2_0 - 1.0/4.0*t161*t190*t31 - 1.0/4.0*t20*v2sigma2_3; + c_grad_rho_s_x_grad_rho_z_y = (1.0/2.0)*f_nabla*grad_rho_s_y*t191*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_x*t110*v2sigma2_1 - 1.0/4.0*grad_rho_z_y*t163 - 1.0/4.0*grad_rho_z_y*t164 + (1.0/4.0)*grad_rho_z_y*t35*v2sigma2_4 - 1.0/4.0*t108*t152 - 1.0/4.0*t108*t162 + (1.0/4.0)*t108*t35*v2sigma2_5 - 1.0/4.0*t110*t167 + (1.0/4.0)*t110*t33*v2sigma2_0 - 1.0/4.0*t183*t192; + c_grad_rho_s_x_grad_rho_z_z = (1.0/2.0)*f_nabla*grad_rho_s_z*t191*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_x*t117*v2sigma2_1 - 1.0/4.0*grad_rho_z_z*t163 - 1.0/4.0*grad_rho_z_z*t164 + (1.0/4.0)*grad_rho_z_z*t35*v2sigma2_4 - 1.0/4.0*t115*t152 - 1.0/4.0*t115*t162 + (1.0/4.0)*t115*t35*v2sigma2_5 - 1.0/4.0*t117*t167 + (1.0/4.0)*t117*t33*v2sigma2_0 - 1.0/4.0*t184*t192; + c_grad_rho_s_y_grad_rho_s_y = (1.0/4.0)*t149 + (1.0/4.0)*t160*(t200 - t201 + 1) + (1.0/4.0)*t161*(-t200 + t201 + 1) + (1.0/4.0)*t193*v2sigma2_3 + (1.0/2.0)*t194*t41 - 1.0/4.0*t195*t196 - 1.0/4.0*t196*t202 + (1.0/4.0)*pow(t41, 2)*v2sigma2_0 + (1.0/4.0)*pow(t43, 2)*v2sigma2_5; + c_grad_rho_s_y_grad_rho_s_z = (1.0/4.0)*grad_rho_s_z*t203 + (1.0/4.0)*grad_rho_s_z*t204 - 1.0/4.0*grad_rho_s_z*t43*v2sigma2_4 + (1.0/4.0)*t160*t210 - 1.0/4.0*t161*t210 + (1.0/4.0)*t194*t49 - 1.0/4.0*t195*t51 - 1.0/4.0*t202*t51 - 1.0/4.0*t205*t49 + (1.0/4.0)*t41*t49*v2sigma2_0 + (1.0/4.0)*t43*t51*v2sigma2_5; + c_grad_rho_s_y_grad_rho_x_x = (1.0/2.0)*f_nabla*grad_rho_s_x*t212*t30*vsigma_0 - 1.0/4.0*grad_rho_s_x*t161*t213 + (1.0/4.0)*grad_rho_s_y*t59*v2sigma2_1 - 1.0/4.0*grad_rho_x_x*t203 - 1.0/4.0*grad_rho_x_x*t204 + (1.0/4.0)*grad_rho_x_x*t43*v2sigma2_4 - 1.0/4.0*t195*t57 - 1.0/4.0*t202*t57 - 1.0/4.0*t205*t59 + (1.0/4.0)*t41*t59*v2sigma2_0 + (1.0/4.0)*t43*t57*v2sigma2_5; + c_grad_rho_s_y_grad_rho_x_y = (1.0/2.0)*f_nabla*t214*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_y*t66*v2sigma2_1 - 1.0/4.0*grad_rho_x_y*t204 + (1.0/4.0)*grad_rho_x_y*t43*v2sigma2_4 - 1.0/4.0*t11*v2sigma2_3 - 1.0/4.0*t161*t214*t31 - 1.0/4.0*t195*t64 - 1.0/4.0*t202*t64 - 1.0/4.0*t205*t66 + (1.0/4.0)*t41*t66*v2sigma2_0 + (1.0/4.0)*t43*t64*v2sigma2_5; + c_grad_rho_s_y_grad_rho_x_z = (1.0/2.0)*f_nabla*grad_rho_s_z*t212*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_y*t73*v2sigma2_1 - 1.0/4.0*grad_rho_x_z*t203 - 1.0/4.0*grad_rho_x_z*t204 + (1.0/4.0)*grad_rho_x_z*t43*v2sigma2_4 - 1.0/4.0*t184*t213 - 1.0/4.0*t195*t71 - 1.0/4.0*t202*t71 - 1.0/4.0*t205*t73 + (1.0/4.0)*t41*t73*v2sigma2_0 + (1.0/4.0)*t43*t71*v2sigma2_5; + c_grad_rho_s_y_grad_rho_y_x = (1.0/2.0)*f_nabla*grad_rho_s_x*t216*t30*vsigma_0 - 1.0/4.0*grad_rho_s_x*t161*t217 + (1.0/4.0)*grad_rho_s_y*t81*v2sigma2_1 - 1.0/4.0*grad_rho_y_x*t203 - 1.0/4.0*grad_rho_y_x*t204 + (1.0/4.0)*grad_rho_y_x*t43*v2sigma2_4 - 1.0/4.0*t195*t79 - 1.0/4.0*t202*t79 - 1.0/4.0*t205*t81 + (1.0/4.0)*t41*t81*v2sigma2_0 + (1.0/4.0)*t43*t79*v2sigma2_5; + c_grad_rho_s_y_grad_rho_y_y = (1.0/2.0)*f_nabla*t218*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_y*t88*v2sigma2_1 - 1.0/4.0*grad_rho_y_y*t204 + (1.0/4.0)*grad_rho_y_y*t43*v2sigma2_4 - 1.0/4.0*t16*v2sigma2_3 - 1.0/4.0*t161*t218*t31 - 1.0/4.0*t195*t86 - 1.0/4.0*t202*t86 - 1.0/4.0*t205*t88 + (1.0/4.0)*t41*t88*v2sigma2_0 + (1.0/4.0)*t43*t86*v2sigma2_5; + c_grad_rho_s_y_grad_rho_y_z = (1.0/2.0)*f_nabla*grad_rho_s_z*t216*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_y*t95*v2sigma2_1 - 1.0/4.0*grad_rho_y_z*t203 - 1.0/4.0*grad_rho_y_z*t204 + (1.0/4.0)*grad_rho_y_z*t43*v2sigma2_4 - 1.0/4.0*t184*t217 - 1.0/4.0*t195*t93 - 1.0/4.0*t202*t93 - 1.0/4.0*t205*t95 + (1.0/4.0)*t41*t95*v2sigma2_0 + (1.0/4.0)*t43*t93*v2sigma2_5; + c_grad_rho_s_y_grad_rho_z_x = (1.0/2.0)*f_nabla*grad_rho_s_x*t220*t30*vsigma_0 - 1.0/4.0*grad_rho_s_x*t161*t221 + (1.0/4.0)*grad_rho_s_y*t103*v2sigma2_1 - 1.0/4.0*grad_rho_z_x*t203 - 1.0/4.0*grad_rho_z_x*t204 + (1.0/4.0)*grad_rho_z_x*t43*v2sigma2_4 - 1.0/4.0*t101*t195 - 1.0/4.0*t101*t202 + (1.0/4.0)*t101*t43*v2sigma2_5 - 1.0/4.0*t103*t205 + (1.0/4.0)*t103*t41*v2sigma2_0; + c_grad_rho_s_y_grad_rho_z_y = (1.0/2.0)*f_nabla*t222*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_y*t110*v2sigma2_1 - 1.0/4.0*grad_rho_z_y*t204 + (1.0/4.0)*grad_rho_z_y*t43*v2sigma2_4 - 1.0/4.0*t108*t195 - 1.0/4.0*t108*t202 + (1.0/4.0)*t108*t43*v2sigma2_5 - 1.0/4.0*t110*t205 + (1.0/4.0)*t110*t41*v2sigma2_0 - 1.0/4.0*t161*t222*t31 - 1.0/4.0*t21*v2sigma2_3; + c_grad_rho_s_y_grad_rho_z_z = (1.0/2.0)*f_nabla*grad_rho_s_z*t220*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_y*t117*v2sigma2_1 - 1.0/4.0*grad_rho_z_z*t203 - 1.0/4.0*grad_rho_z_z*t204 + (1.0/4.0)*grad_rho_z_z*t43*v2sigma2_4 - 1.0/4.0*t115*t195 - 1.0/4.0*t115*t202 + (1.0/4.0)*t115*t43*v2sigma2_5 - 1.0/4.0*t117*t205 + (1.0/4.0)*t117*t41*v2sigma2_0 - 1.0/4.0*t184*t221; + c_grad_rho_s_z_grad_rho_s_z = (1.0/2.0)*grad_rho_s_z*t49*v2sigma2_1 - 1.0/2.0*grad_rho_s_z*t51*v2sigma2_4 + (1.0/4.0)*t149 + (1.0/4.0)*t160*(t227 - t228 + 1) + (1.0/4.0)*t161*(-t227 + t228 + 1) + (1.0/4.0)*t223*v2sigma2_3 - 1.0/2.0*t229*t51 + (1.0/4.0)*pow(t49, 2)*v2sigma2_0 + (1.0/4.0)*pow(t51, 2)*v2sigma2_5; + c_grad_rho_s_z_grad_rho_x_x = (1.0/2.0)*f_nabla*grad_rho_s_x*t235*t30*vsigma_0 - 1.0/4.0*grad_rho_s_x*t161*t236 + (1.0/4.0)*grad_rho_s_z*t59*v2sigma2_1 - 1.0/4.0*grad_rho_x_x*t230 - 1.0/4.0*grad_rho_x_x*t232 + (1.0/4.0)*grad_rho_x_x*t51*v2sigma2_4 - 1.0/4.0*t229*t57 - 1.0/4.0*t231*t57 - 1.0/4.0*t237*t59 + (1.0/4.0)*t49*t59*v2sigma2_0 + (1.0/4.0)*t51*t57*v2sigma2_5; + c_grad_rho_s_z_grad_rho_x_y = (1.0/2.0)*f_nabla*grad_rho_s_y*t235*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_z*t66*v2sigma2_1 - 1.0/4.0*grad_rho_x_y*t230 - 1.0/4.0*grad_rho_x_y*t232 + (1.0/4.0)*grad_rho_x_y*t51*v2sigma2_4 - 1.0/4.0*t183*t236 - 1.0/4.0*t229*t64 - 1.0/4.0*t231*t64 - 1.0/4.0*t237*t66 + (1.0/4.0)*t49*t66*v2sigma2_0 + (1.0/4.0)*t51*t64*v2sigma2_5; + c_grad_rho_s_z_grad_rho_x_z = (1.0/2.0)*f_nabla*t238*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_z*t73*v2sigma2_1 - 1.0/4.0*grad_rho_x_z*t232 + (1.0/4.0)*grad_rho_x_z*t51*v2sigma2_4 - 1.0/4.0*t12*v2sigma2_3 - 1.0/4.0*t161*t238*t31 - 1.0/4.0*t229*t71 - 1.0/4.0*t231*t71 - 1.0/4.0*t237*t73 + (1.0/4.0)*t49*t73*v2sigma2_0 + (1.0/4.0)*t51*t71*v2sigma2_5; + c_grad_rho_s_z_grad_rho_y_x = (1.0/2.0)*f_nabla*grad_rho_s_x*t240*t30*vsigma_0 - 1.0/4.0*grad_rho_s_x*t161*t241 + (1.0/4.0)*grad_rho_s_z*t81*v2sigma2_1 - 1.0/4.0*grad_rho_y_x*t230 - 1.0/4.0*grad_rho_y_x*t232 + (1.0/4.0)*grad_rho_y_x*t51*v2sigma2_4 - 1.0/4.0*t229*t79 - 1.0/4.0*t231*t79 - 1.0/4.0*t237*t81 + (1.0/4.0)*t49*t81*v2sigma2_0 + (1.0/4.0)*t51*t79*v2sigma2_5; + c_grad_rho_s_z_grad_rho_y_y = (1.0/2.0)*f_nabla*grad_rho_s_y*t240*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_z*t88*v2sigma2_1 - 1.0/4.0*grad_rho_y_y*t230 - 1.0/4.0*grad_rho_y_y*t232 + (1.0/4.0)*grad_rho_y_y*t51*v2sigma2_4 - 1.0/4.0*t183*t241 - 1.0/4.0*t229*t86 - 1.0/4.0*t231*t86 - 1.0/4.0*t237*t88 + (1.0/4.0)*t49*t88*v2sigma2_0 + (1.0/4.0)*t51*t86*v2sigma2_5; + c_grad_rho_s_z_grad_rho_y_z = (1.0/2.0)*f_nabla*t242*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_z*t95*v2sigma2_1 - 1.0/4.0*grad_rho_y_z*t232 + (1.0/4.0)*grad_rho_y_z*t51*v2sigma2_4 - 1.0/4.0*t161*t242*t31 - 1.0/4.0*t17*v2sigma2_3 - 1.0/4.0*t229*t93 - 1.0/4.0*t231*t93 - 1.0/4.0*t237*t95 + (1.0/4.0)*t49*t95*v2sigma2_0 + (1.0/4.0)*t51*t93*v2sigma2_5; + c_grad_rho_s_z_grad_rho_z_x = (1.0/2.0)*f_nabla*grad_rho_s_x*t244*t30*vsigma_0 - 1.0/4.0*grad_rho_s_x*t161*t245 + (1.0/4.0)*grad_rho_s_z*t103*v2sigma2_1 - 1.0/4.0*grad_rho_z_x*t230 - 1.0/4.0*grad_rho_z_x*t232 + (1.0/4.0)*grad_rho_z_x*t51*v2sigma2_4 - 1.0/4.0*t101*t229 - 1.0/4.0*t101*t231 + (1.0/4.0)*t101*t51*v2sigma2_5 - 1.0/4.0*t103*t237 + (1.0/4.0)*t103*t49*v2sigma2_0; + c_grad_rho_s_z_grad_rho_z_y = (1.0/2.0)*f_nabla*grad_rho_s_y*t244*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_z*t110*v2sigma2_1 - 1.0/4.0*grad_rho_z_y*t230 - 1.0/4.0*grad_rho_z_y*t232 + (1.0/4.0)*grad_rho_z_y*t51*v2sigma2_4 - 1.0/4.0*t108*t229 - 1.0/4.0*t108*t231 + (1.0/4.0)*t108*t51*v2sigma2_5 - 1.0/4.0*t110*t237 + (1.0/4.0)*t110*t49*v2sigma2_0 - 1.0/4.0*t183*t245; + c_grad_rho_s_z_grad_rho_z_z = (1.0/2.0)*f_nabla*t246*t30*vsigma_0 + (1.0/4.0)*grad_rho_s_z*t117*v2sigma2_1 - 1.0/4.0*grad_rho_z_z*t232 + (1.0/4.0)*grad_rho_z_z*t51*v2sigma2_4 - 1.0/4.0*t115*t229 - 1.0/4.0*t115*t231 + (1.0/4.0)*t115*t51*v2sigma2_5 - 1.0/4.0*t117*t237 + (1.0/4.0)*t117*t49*v2sigma2_0 - 1.0/4.0*t161*t246*t31 - 1.0/4.0*t22*v2sigma2_3; + c_grad_rho_x_x_grad_rho_x_x = (1.0/4.0)*t154*v2sigma2_3 + (1.0/4.0)*t160*(-t252 + t254) + (1.0/4.0)*t161*(t252 + t255) + (1.0/4.0)*t247 - 1.0/2.0*t248*t59 + (1.0/4.0)*t249*t250 - 1.0/4.0*t250*t256 + (1.0/4.0)*pow(t57, 2)*v2sigma2_5 + (1.0/4.0)*pow(t59, 2)*v2sigma2_0; + c_grad_rho_x_x_grad_rho_x_y = -1.0/4.0*grad_rho_x_y*t257 + (1.0/4.0)*grad_rho_x_y*t258 + (1.0/4.0)*t169*v2sigma2_3 + (1.0/4.0)*t183*t261 - 1.0/4.0*t248*t66 + (1.0/4.0)*t249*t64 - 1.0/4.0*t256*t64 - 1.0/4.0*t259*t260 + (1.0/4.0)*t262*t66 - 1.0/4.0*t263*t66 + (1.0/4.0)*t264*t64; + c_grad_rho_x_x_grad_rho_x_z = -1.0/4.0*grad_rho_x_z*t257 + (1.0/4.0)*grad_rho_x_z*t258 + (1.0/4.0)*t175*v2sigma2_3 + (1.0/4.0)*t184*t261 - 1.0/4.0*t248*t73 + (1.0/4.0)*t249*t71 - 1.0/4.0*t256*t71 - 1.0/4.0*t261*t265 + (1.0/4.0)*t262*t73 - 1.0/4.0*t263*t73 + (1.0/4.0)*t264*t71; + c_grad_rho_x_x_grad_rho_y_x = -1.0/4.0*grad_rho_y_x*t257 + (1.0/4.0)*grad_rho_y_x*t258 + (1.0/4.0)*grad_rho_y_x*t266 + (1.0/4.0)*t19*t270 - 1.0/4.0*t248*t81 + (1.0/4.0)*t249*t79 - 1.0/4.0*t256*t79 + (1.0/4.0)*t262*t81 - 1.0/4.0*t263*t81 + (1.0/4.0)*t264*t79 - 1.0/4.0*t268*t269; + c_grad_rho_x_x_grad_rho_y_y = -1.0/4.0*grad_rho_y_y*t257 + (1.0/4.0)*grad_rho_y_y*t258 + (1.0/4.0)*grad_rho_y_y*t266 - 1.0/4.0*t248*t88 + (1.0/4.0)*t249*t86 - 1.0/4.0*t256*t86 + (1.0/4.0)*t262*t88 - 1.0/4.0*t263*t88 + (1.0/4.0)*t264*t86 + (1.0/4.0)*t273; + c_grad_rho_x_x_grad_rho_y_z = -1.0/4.0*grad_rho_y_z*t257 + (1.0/4.0)*grad_rho_y_z*t258 + (1.0/4.0)*grad_rho_y_z*t266 - 1.0/4.0*t248*t95 + (1.0/4.0)*t249*t93 - 1.0/4.0*t256*t93 + (1.0/4.0)*t262*t95 - 1.0/4.0*t263*t95 + (1.0/4.0)*t264*t93 + (1.0/4.0)*t274; + c_grad_rho_x_x_grad_rho_z_x = -1.0/4.0*grad_rho_z_x*t257 + (1.0/4.0)*grad_rho_z_x*t258 + (1.0/4.0)*grad_rho_z_x*t266 + (1.0/4.0)*t101*t249 - 1.0/4.0*t101*t256 + (1.0/4.0)*t101*t264 - 1.0/4.0*t103*t248 + (1.0/4.0)*t103*t262 - 1.0/4.0*t103*t263 + (1.0/4.0)*t24*t270 - 1.0/4.0*t268*t275; + c_grad_rho_x_x_grad_rho_z_y = -1.0/4.0*grad_rho_z_y*t257 + (1.0/4.0)*grad_rho_z_y*t258 + (1.0/4.0)*grad_rho_z_y*t266 + (1.0/4.0)*t108*t249 - 1.0/4.0*t108*t256 + (1.0/4.0)*t108*t264 - 1.0/4.0*t110*t248 + (1.0/4.0)*t110*t262 - 1.0/4.0*t110*t263 + (1.0/4.0)*t278; + c_grad_rho_x_x_grad_rho_z_z = -1.0/4.0*grad_rho_z_z*t257 + (1.0/4.0)*grad_rho_z_z*t258 + (1.0/4.0)*grad_rho_z_z*t266 + (1.0/4.0)*t115*t249 - 1.0/4.0*t115*t256 + (1.0/4.0)*t115*t264 - 1.0/4.0*t117*t248 + (1.0/4.0)*t117*t262 - 1.0/4.0*t117*t263 + (1.0/4.0)*t279; + c_grad_rho_x_y_grad_rho_x_y = (1.0/4.0)*t160*(-t283 + t285) + (1.0/4.0)*t161*(t283 + t286) + (1.0/4.0)*t197*v2sigma2_3 + (1.0/4.0)*t247 - 1.0/2.0*t280*t66 + (1.0/4.0)*t281*t282 - 1.0/4.0*t282*t287 + (1.0/4.0)*pow(t64, 2)*v2sigma2_5 + (1.0/4.0)*pow(t66, 2)*v2sigma2_0; + c_grad_rho_x_y_grad_rho_x_z = -1.0/4.0*grad_rho_x_z*t288 + (1.0/4.0)*grad_rho_x_z*t289 + (1.0/4.0)*t206*v2sigma2_3 - 1.0/4.0*t259*t290 + (1.0/4.0)*t259*t291 - 1.0/4.0*t280*t73 + (1.0/4.0)*t281*t71 - 1.0/4.0*t287*t71 + (1.0/4.0)*t292*t73 - 1.0/4.0*t293*t73 + (1.0/4.0)*t294*t71; + c_grad_rho_x_y_grad_rho_y_x = -1.0/4.0*grad_rho_y_x*t288 + (1.0/4.0)*grad_rho_y_x*t289 + (1.0/4.0)*grad_rho_y_x*t295 + (1.0/4.0)*t273 - 1.0/4.0*t280*t81 + (1.0/4.0)*t281*t79 - 1.0/4.0*t287*t79 + (1.0/4.0)*t292*t81 - 1.0/4.0*t293*t81 + (1.0/4.0)*t294*t79; + c_grad_rho_x_y_grad_rho_y_y = -1.0/4.0*grad_rho_y_y*t288 + (1.0/4.0)*grad_rho_y_y*t289 + (1.0/4.0)*grad_rho_y_y*t295 + (1.0/4.0)*t19*t297 - 1.0/4.0*t269*t296 - 1.0/4.0*t280*t88 + (1.0/4.0)*t281*t86 - 1.0/4.0*t287*t86 + (1.0/4.0)*t292*t88 - 1.0/4.0*t293*t88 + (1.0/4.0)*t294*t86; + c_grad_rho_x_y_grad_rho_y_z = -1.0/4.0*grad_rho_y_z*t288 + (1.0/4.0)*grad_rho_y_z*t289 + (1.0/4.0)*grad_rho_y_z*t295 - 1.0/4.0*t280*t95 + (1.0/4.0)*t281*t93 - 1.0/4.0*t287*t93 + (1.0/4.0)*t292*t95 - 1.0/4.0*t293*t95 + (1.0/4.0)*t294*t93 + (1.0/4.0)*t298; + c_grad_rho_x_y_grad_rho_z_x = -1.0/4.0*grad_rho_z_x*t288 + (1.0/4.0)*grad_rho_z_x*t289 + (1.0/4.0)*grad_rho_z_x*t295 + (1.0/4.0)*t101*t281 - 1.0/4.0*t101*t287 + (1.0/4.0)*t101*t294 - 1.0/4.0*t103*t280 + (1.0/4.0)*t103*t292 - 1.0/4.0*t103*t293 + (1.0/4.0)*t278; + c_grad_rho_x_y_grad_rho_z_y = -1.0/4.0*grad_rho_z_y*t288 + (1.0/4.0)*grad_rho_z_y*t289 + (1.0/4.0)*grad_rho_z_y*t295 + (1.0/4.0)*t108*t281 - 1.0/4.0*t108*t287 + (1.0/4.0)*t108*t294 - 1.0/4.0*t110*t280 + (1.0/4.0)*t110*t292 - 1.0/4.0*t110*t293 + (1.0/4.0)*t24*t297 - 1.0/4.0*t275*t296; + c_grad_rho_x_y_grad_rho_z_z = -1.0/4.0*grad_rho_z_z*t288 + (1.0/4.0)*grad_rho_z_z*t289 + (1.0/4.0)*grad_rho_z_z*t295 + (1.0/4.0)*t115*t281 - 1.0/4.0*t115*t287 + (1.0/4.0)*t115*t294 - 1.0/4.0*t117*t280 + (1.0/4.0)*t117*t292 - 1.0/4.0*t117*t293 + (1.0/4.0)*t299; + c_grad_rho_x_z_grad_rho_x_z = (1.0/4.0)*t160*(-t303 + t305) + (1.0/4.0)*t161*(t303 + t306) + (1.0/4.0)*t224*v2sigma2_3 + (1.0/4.0)*t247 - 1.0/2.0*t300*t73 + (1.0/4.0)*t301*t302 - 1.0/4.0*t302*t307 + (1.0/4.0)*pow(t71, 2)*v2sigma2_5 + (1.0/4.0)*pow(t73, 2)*v2sigma2_0; + c_grad_rho_x_z_grad_rho_y_x = (1.0/4.0)*grad_rho_y_x*t308 + (1.0/4.0)*grad_rho_y_x*t309 - 1.0/4.0*grad_rho_y_x*t312 + (1.0/4.0)*t274 - 1.0/4.0*t300*t81 + (1.0/4.0)*t301*t79 - 1.0/4.0*t307*t79 + (1.0/4.0)*t310*t81 + (1.0/4.0)*t311*t79 - 1.0/4.0*t313*t81; + c_grad_rho_x_z_grad_rho_y_y = (1.0/4.0)*grad_rho_y_y*t308 + (1.0/4.0)*grad_rho_y_y*t309 - 1.0/4.0*grad_rho_y_y*t312 + (1.0/4.0)*t298 - 1.0/4.0*t300*t88 + (1.0/4.0)*t301*t86 - 1.0/4.0*t307*t86 + (1.0/4.0)*t310*t88 + (1.0/4.0)*t311*t86 - 1.0/4.0*t313*t88; + c_grad_rho_x_z_grad_rho_y_z = (1.0/4.0)*grad_rho_y_z*t308 + (1.0/4.0)*grad_rho_y_z*t309 - 1.0/4.0*grad_rho_y_z*t312 + (1.0/4.0)*t19*t315 - 1.0/4.0*t269*t314 - 1.0/4.0*t300*t95 + (1.0/4.0)*t301*t93 - 1.0/4.0*t307*t93 + (1.0/4.0)*t310*t95 + (1.0/4.0)*t311*t93 - 1.0/4.0*t313*t95; + c_grad_rho_x_z_grad_rho_z_x = (1.0/4.0)*grad_rho_z_x*t308 + (1.0/4.0)*grad_rho_z_x*t309 - 1.0/4.0*grad_rho_z_x*t312 + (1.0/4.0)*t101*t301 - 1.0/4.0*t101*t307 + (1.0/4.0)*t101*t311 - 1.0/4.0*t103*t300 + (1.0/4.0)*t103*t310 - 1.0/4.0*t103*t313 + (1.0/4.0)*t279; + c_grad_rho_x_z_grad_rho_z_y = (1.0/4.0)*grad_rho_z_y*t308 + (1.0/4.0)*grad_rho_z_y*t309 - 1.0/4.0*grad_rho_z_y*t312 + (1.0/4.0)*t108*t301 - 1.0/4.0*t108*t307 + (1.0/4.0)*t108*t311 - 1.0/4.0*t110*t300 + (1.0/4.0)*t110*t310 - 1.0/4.0*t110*t313 + (1.0/4.0)*t299; + c_grad_rho_x_z_grad_rho_z_z = (1.0/4.0)*grad_rho_z_z*t308 + (1.0/4.0)*grad_rho_z_z*t309 - 1.0/4.0*grad_rho_z_z*t312 + (1.0/4.0)*t115*t301 - 1.0/4.0*t115*t307 + (1.0/4.0)*t115*t311 - 1.0/4.0*t117*t300 + (1.0/4.0)*t117*t310 - 1.0/4.0*t117*t313 + (1.0/4.0)*t24*t315 - 1.0/4.0*t275*t314; + c_grad_rho_y_x_grad_rho_y_x = (1.0/4.0)*t155*v2sigma2_3 + (1.0/4.0)*t160*(t254 - t320) + (1.0/4.0)*t161*(t255 + t320) + (1.0/4.0)*t247 - 1.0/2.0*t316*t81 + (1.0/4.0)*t317*t318 - 1.0/4.0*t318*t321 + (1.0/4.0)*pow(t79, 2)*v2sigma2_5 + (1.0/4.0)*pow(t81, 2)*v2sigma2_0; + c_grad_rho_y_x_grad_rho_y_y = -1.0/4.0*grad_rho_y_y*t322 + (1.0/4.0)*grad_rho_y_y*t323 + (1.0/4.0)*t170*v2sigma2_3 + (1.0/4.0)*t183*t325 - 1.0/4.0*t260*t324 - 1.0/4.0*t316*t88 + (1.0/4.0)*t317*t86 - 1.0/4.0*t321*t86 + (1.0/4.0)*t326*t88 - 1.0/4.0*t327*t88 + (1.0/4.0)*t328*t86; + c_grad_rho_y_x_grad_rho_y_z = -1.0/4.0*grad_rho_y_z*t322 + (1.0/4.0)*grad_rho_y_z*t323 + (1.0/4.0)*t176*v2sigma2_3 + (1.0/4.0)*t184*t325 - 1.0/4.0*t265*t325 - 1.0/4.0*t316*t95 + (1.0/4.0)*t317*t93 - 1.0/4.0*t321*t93 + (1.0/4.0)*t326*t95 - 1.0/4.0*t327*t95 + (1.0/4.0)*t328*t93; + c_grad_rho_y_x_grad_rho_z_x = -1.0/4.0*grad_rho_z_x*t322 + (1.0/4.0)*grad_rho_z_x*t323 + (1.0/4.0)*grad_rho_z_x*t329 + (1.0/4.0)*t101*t317 - 1.0/4.0*t101*t321 + (1.0/4.0)*t101*t328 - 1.0/4.0*t103*t316 + (1.0/4.0)*t103*t326 - 1.0/4.0*t103*t327 - 1.0/4.0*t269*t331 + (1.0/4.0)*t331*t332; + c_grad_rho_y_x_grad_rho_z_y = -1.0/4.0*grad_rho_z_y*t322 + (1.0/4.0)*grad_rho_z_y*t323 + (1.0/4.0)*grad_rho_z_y*t329 + (1.0/4.0)*t108*t317 - 1.0/4.0*t108*t321 + (1.0/4.0)*t108*t328 - 1.0/4.0*t110*t316 + (1.0/4.0)*t110*t326 - 1.0/4.0*t110*t327 + (1.0/4.0)*t335; + c_grad_rho_y_x_grad_rho_z_z = -1.0/4.0*grad_rho_z_z*t322 + (1.0/4.0)*grad_rho_z_z*t323 + (1.0/4.0)*grad_rho_z_z*t329 + (1.0/4.0)*t115*t317 - 1.0/4.0*t115*t321 + (1.0/4.0)*t115*t328 - 1.0/4.0*t117*t316 + (1.0/4.0)*t117*t326 - 1.0/4.0*t117*t327 + (1.0/4.0)*t336; + c_grad_rho_y_y_grad_rho_y_y = (1.0/4.0)*t160*(t285 - t340) + (1.0/4.0)*t161*(t286 + t340) + (1.0/4.0)*t198*v2sigma2_3 + (1.0/4.0)*t247 - 1.0/2.0*t337*t88 + (1.0/4.0)*t338*t339 - 1.0/4.0*t339*t341 + (1.0/4.0)*pow(t86, 2)*v2sigma2_5 + (1.0/4.0)*pow(t88, 2)*v2sigma2_0; + c_grad_rho_y_y_grad_rho_y_z = -1.0/4.0*grad_rho_y_z*t342 + (1.0/4.0)*grad_rho_y_z*t343 + (1.0/4.0)*t207*v2sigma2_3 - 1.0/4.0*t290*t324 + (1.0/4.0)*t291*t324 - 1.0/4.0*t337*t95 + (1.0/4.0)*t338*t93 - 1.0/4.0*t341*t93 + (1.0/4.0)*t344*t95 - 1.0/4.0*t345*t95 + (1.0/4.0)*t346*t93; + c_grad_rho_y_y_grad_rho_z_x = -1.0/4.0*grad_rho_z_x*t342 + (1.0/4.0)*grad_rho_z_x*t343 + (1.0/4.0)*grad_rho_z_x*t347 + (1.0/4.0)*t101*t338 - 1.0/4.0*t101*t341 + (1.0/4.0)*t101*t346 - 1.0/4.0*t103*t337 + (1.0/4.0)*t103*t344 - 1.0/4.0*t103*t345 + (1.0/4.0)*t335; + c_grad_rho_y_y_grad_rho_z_y = -1.0/4.0*grad_rho_z_y*t342 + (1.0/4.0)*grad_rho_z_y*t343 + (1.0/4.0)*grad_rho_z_y*t347 + (1.0/4.0)*t108*t338 - 1.0/4.0*t108*t341 + (1.0/4.0)*t108*t346 - 1.0/4.0*t110*t337 + (1.0/4.0)*t110*t344 - 1.0/4.0*t110*t345 - 1.0/4.0*t269*t348 + (1.0/4.0)*t332*t348; + c_grad_rho_y_y_grad_rho_z_z = -1.0/4.0*grad_rho_z_z*t342 + (1.0/4.0)*grad_rho_z_z*t343 + (1.0/4.0)*grad_rho_z_z*t347 + (1.0/4.0)*t115*t338 - 1.0/4.0*t115*t341 + (1.0/4.0)*t115*t346 - 1.0/4.0*t117*t337 + (1.0/4.0)*t117*t344 - 1.0/4.0*t117*t345 + (1.0/4.0)*t349; + c_grad_rho_y_z_grad_rho_y_z = (1.0/4.0)*t160*(t305 - t353) + (1.0/4.0)*t161*(t306 + t353) + (1.0/4.0)*t225*v2sigma2_3 + (1.0/4.0)*t247 - 1.0/2.0*t350*t95 + (1.0/4.0)*t351*t352 - 1.0/4.0*t352*t354 + (1.0/4.0)*pow(t93, 2)*v2sigma2_5 + (1.0/4.0)*pow(t95, 2)*v2sigma2_0; + c_grad_rho_y_z_grad_rho_z_x = (1.0/4.0)*grad_rho_z_x*t355 + (1.0/4.0)*grad_rho_z_x*t356 - 1.0/4.0*grad_rho_z_x*t359 + (1.0/4.0)*t101*t351 - 1.0/4.0*t101*t354 + (1.0/4.0)*t101*t358 - 1.0/4.0*t103*t350 + (1.0/4.0)*t103*t357 - 1.0/4.0*t103*t360 + (1.0/4.0)*t336; + c_grad_rho_y_z_grad_rho_z_y = (1.0/4.0)*grad_rho_z_y*t355 + (1.0/4.0)*grad_rho_z_y*t356 - 1.0/4.0*grad_rho_z_y*t359 + (1.0/4.0)*t108*t351 - 1.0/4.0*t108*t354 + (1.0/4.0)*t108*t358 - 1.0/4.0*t110*t350 + (1.0/4.0)*t110*t357 - 1.0/4.0*t110*t360 + (1.0/4.0)*t349; + c_grad_rho_y_z_grad_rho_z_z = (1.0/4.0)*grad_rho_z_z*t355 + (1.0/4.0)*grad_rho_z_z*t356 - 1.0/4.0*grad_rho_z_z*t359 + (1.0/4.0)*t115*t351 - 1.0/4.0*t115*t354 + (1.0/4.0)*t115*t358 - 1.0/4.0*t117*t350 + (1.0/4.0)*t117*t357 - 1.0/4.0*t117*t360 - 1.0/4.0*t269*t361 + (1.0/4.0)*t332*t361; + c_grad_rho_z_x_grad_rho_z_x = (1.0/4.0)*pow(t101, 2)*v2sigma2_5 + (1.0/4.0)*pow(t103, 2)*v2sigma2_0 - 1.0/2.0*t103*t362 + (1.0/4.0)*t156*v2sigma2_3 + (1.0/4.0)*t160*(t254 - t366) + (1.0/4.0)*t161*(t255 + t366) + (1.0/4.0)*t247 + (1.0/4.0)*t363*t364 - 1.0/4.0*t364*t367; + c_grad_rho_z_x_grad_rho_z_y = -1.0/4.0*grad_rho_z_y*t368 + (1.0/4.0)*grad_rho_z_y*t369 + (1.0/4.0)*t108*t363 - 1.0/4.0*t108*t367 + (1.0/4.0)*t108*t374 - 1.0/4.0*t110*t362 + (1.0/4.0)*t110*t372 - 1.0/4.0*t110*t373 + (1.0/4.0)*t171*v2sigma2_3 + (1.0/4.0)*t183*t371 - 1.0/4.0*t260*t370; + c_grad_rho_z_x_grad_rho_z_z = -1.0/4.0*grad_rho_z_z*t368 + (1.0/4.0)*grad_rho_z_z*t369 + (1.0/4.0)*t115*t363 - 1.0/4.0*t115*t367 + (1.0/4.0)*t115*t374 - 1.0/4.0*t117*t362 + (1.0/4.0)*t117*t372 - 1.0/4.0*t117*t373 + (1.0/4.0)*t177*v2sigma2_3 + (1.0/4.0)*t184*t371 - 1.0/4.0*t265*t371; + c_grad_rho_z_y_grad_rho_z_y = (1.0/4.0)*pow(t108, 2)*v2sigma2_5 + (1.0/4.0)*pow(t110, 2)*v2sigma2_0 - 1.0/2.0*t110*t375 + (1.0/4.0)*t160*(t285 - t378) + (1.0/4.0)*t161*(t286 + t378) + (1.0/4.0)*t199*v2sigma2_3 + (1.0/4.0)*t247 + (1.0/4.0)*t376*t377 - 1.0/4.0*t377*t379; + c_grad_rho_z_y_grad_rho_z_z = (1.0/4.0)*t108*t115*v2sigma2_5 + (1.0/4.0)*t108*t381 - 1.0/4.0*t108*t382 + (1.0/4.0)*t110*t117*v2sigma2_0 - 1.0/4.0*t110*t380 + (1.0/4.0)*t115*t376 - 1.0/4.0*t115*t379 - 1.0/4.0*t117*t375 + (1.0/4.0)*t208*v2sigma2_3 - 1.0/4.0*t290*t370 + (1.0/4.0)*t291*t370; + c_grad_rho_z_z_grad_rho_z_z = (1.0/4.0)*pow(t115, 2)*v2sigma2_5 + (1.0/4.0)*pow(t117, 2)*v2sigma2_0 - 1.0/2.0*t117*t380 + (1.0/4.0)*t160*(t305 - t384) + (1.0/4.0)*t161*(t306 + t384) + (1.0/4.0)*t226*v2sigma2_3 + (1.0/4.0)*t247 + (1.0/4.0)*t381*t383 - 1.0/4.0*t382*t383; +} + +} // namespace xckernel diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.cxx b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.cxx index 6695d9121..d78388fdb 100644 --- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.cxx +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.cxx @@ -15,6 +15,7 @@ #include "reference_replicated_xc_host_integrator_exc_grad.hpp" #include "reference_replicated_xc_host_integrator_exx.hpp" #include "reference_replicated_xc_host_integrator_fxc_contraction.hpp" +#include "reference_replicated_xc_host_integrator_fxc_contraction_gks.hpp" #include "reference_replicated_xc_host_integrator_dd_psi.hpp" #include "reference_replicated_xc_host_integrator_dd_psi_potential.hpp" diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.hpp index a32748eb5..466751f85 100644 --- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.hpp +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.hpp @@ -104,6 +104,16 @@ class ReferenceReplicatedXCHostIntegrator : value_type* FXCz, int64_t ldfxcz, const IntegratorSettingsXC& ks_settings ) override; + // GKS FXC contraction (LDA) + void eval_fxc_contraction_( 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, + const value_type* tPs, int64_t ldtps, const value_type* tPz, int64_t ldtpz, + const value_type* tPy, int64_t ldtpy, const value_type* tPx, int64_t ldtpx, + value_type* FXCs, int64_t ldfxcs, value_type* FXCz, int64_t ldfxcz, + value_type* FXCy, int64_t ldfxcy, value_type* FXCx, int64_t ldfxcx, + const IntegratorSettingsXC& ks_settings ) override; + /// ddX PSi void eval_dd_psi_( int64_t m, int64_t n, const value_type* P, int64_t ldp, unsigned max_Ylm, value_type* ddPsi, int64_t ldPsi ) override; @@ -145,6 +155,14 @@ class ReferenceReplicatedXCHostIntegrator : value_type *N_EL, const IntegratorSettingsXC& ks_settings, task_iterator task_begin, task_iterator task_end ); + // Implementation details of GKS FXC contraction + void fxc_contraction_gks_local_work_( const basis_type& basis, + const value_type* const P[4], const int64_t ldP[4], + const value_type* const tP[4], const int64_t ldtP[4], + value_type* const FXC[4], const int64_t ldFXC[4], + const IntegratorSettingsXC& ks_settings, + task_iterator task_begin, task_iterator task_end ); + // Implementation details of ddX Psi void dd_psi_local_work_( const value_type* P, int64_t ldp, unsigned max_Ylm, value_type* ddPsi, int64_t ldPsi ); diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction_gks.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction_gks.hpp new file mode 100644 index 000000000..f5ea20891 --- /dev/null +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction_gks.hpp @@ -0,0 +1,220 @@ +/** + * GauXC Copyright (c) 2020-2024, The Regents of the University of California, + * through Lawrence Berkeley National Laboratory (subject to receipt of + * any required approvals from the U.S. Dept. of Energy). + * + * (c) 2024-2025, Microsoft Corporation + * + * All rights reserved. + * + * See LICENSE.txt for details + */ +#pragma once + +#include "reference_replicated_xc_host_integrator.hpp" +#include "integrator_util/integrator_common.hpp" +#include "host/local_host_work_driver.hpp" +#include "host/blas.hpp" +#include "host/gauxc_nc_kernel.hpp" +#include +#include + +namespace GauXC::detail { + +/** + * GKS (two-component, noncollinear) FXC contraction. + * + * The energy is the locally collinear one exc_vxc evaluates for GKS: + * n_+- = (rho_s +- |m|)/2 fed to the spin-polarized functional. The + * kernel applied to a trial density, per noncollinear field slot, is + * GENERATED (xckernel ncwriter: the mechanical second derivative of + * that map), and assembled exactly as the potential is -- so + * FXC_X = d/dh VXC_X(P + h tP) for X = s, z, y, x. + * + * Below gks_dtol the map is not twice differentiable (the transverse + * kernel carries 1/|m|); there the generated collinear limit is used: + * every magnetization component responds like the spin channel of a + * collinear perturbation about the spin-symmetric reference. + * + * LDA only for now. The GGA map of Scalmani and Frisch has a second + * singularity, 1/|(grad rho_s . grad m_J)_J|, which vanishes at density + * critical points at finite |m|; its regularization is still open. + * + * Argument order follows the GKS eval_exc_vxc: (s, z, y, x). + */ +template +void ReferenceReplicatedXCHostIntegrator:: + eval_fxc_contraction_( 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, + const value_type* tPs, int64_t ldtps, const value_type* tPz, int64_t ldtpz, + const value_type* tPy, int64_t ldtpy, const value_type* tPx, int64_t ldtpx, + value_type* FXCs, int64_t ldfxcs, value_type* FXCz, int64_t ldfxcz, + value_type* FXCy, int64_t ldfxcy, value_type* FXCx, int64_t ldfxcx, + const IntegratorSettingsXC& ks_settings ) { + + const auto& basis = this->load_balancer_->basis(); + const int64_t nbf = basis.nbf(); + if( m != n ) GAUXC_GENERIC_EXCEPTION("P/FXC Must Be Square"); + if( m != nbf ) GAUXC_GENERIC_EXCEPTION("P/FXC Must Have Same Dimension as Basis"); + for( int64_t ld : { ldps, ldpz, ldpy, ldpx, ldtps, ldtpz, ldtpy, ldtpx, + ldfxcs, ldfxcz, ldfxcy, ldfxcx } ) + if( ld < nbf ) GAUXC_GENERIC_EXCEPTION("Invalid Leading Dimension"); + + // Symmetrize the trial densities: Exc sees only the symmetric part of a + // density matrix, but the gradient channel of a noncollinear GGA would + // not (see #225); done here so the contraction never depends on it. + std::vector tsym[4]; + const value_type* tPin[4] = { tPs, tPz, tPy, tPx }; + const int64_t ldtin[4] = { ldtps, ldtpz, ldtpy, ldtpx }; + const value_type* tP[4]; + int64_t ldtP[4]; + for( int k = 0; k < 4; ++k ) { + tsym[k].resize( nbf*nbf ); + for( int64_t j = 0; j < nbf; ++j ) + for( int64_t i = 0; i < nbf; ++i ) + tsym[k][i + j*nbf] = 0.5 * ( tPin[k][i + j*ldtin[k]] + tPin[k][j + i*ldtin[k]] ); + tP[k] = tsym[k].data(); ldtP[k] = nbf; + } + const value_type* P[4] = { Ps, Pz, Py, Px }; + const int64_t ldP[4] = { ldps, ldpz, ldpy, ldpx }; + value_type* FXC[4] = { FXCs, FXCz, FXCy, FXCx }; + const int64_t ldFXC[4] = { ldfxcs, ldfxcz, ldfxcy, ldfxcx }; + + auto& tasks = this->load_balancer_->get_tasks(); + this->timer_.time_op("XCIntegrator.LocalWork", [&](){ + fxc_contraction_gks_local_work_( basis, P, ldP, tP, ldtP, FXC, ldFXC, + ks_settings, tasks.begin(), tasks.end() ); + }); + + this->timer_.time_op("XCIntegrator.Allreduce", [&](){ + if( not this->reduction_driver_->takes_host_memory() ) + GAUXC_GENERIC_EXCEPTION("This Module Only Works With Host Reductions"); + for( int k = 0; k < 4; ++k ) + this->reduction_driver_->allreduce_inplace( FXC[k], nbf*nbf, ReductionOp::Sum ); + }); +} + + +template +void ReferenceReplicatedXCHostIntegrator:: + fxc_contraction_gks_local_work_( const basis_type& basis, + const value_type* const P[4], const int64_t ldP[4], + const value_type* const tP[4], const int64_t ldtP[4], + value_type* const FXC[4], const int64_t ldFXC[4], + const IntegratorSettingsXC& settings, + task_iterator task_begin, task_iterator task_end ) { + + IntegratorSettingsKS ks_settings; + if( auto* tmp = dynamic_cast(&settings) ) + ks_settings = *tmp; + const double dtol = ks_settings.gks_dtol; + + auto* lwd = dynamic_cast(this->local_work_driver_.get()); + const auto& func = *this->func_; + const auto& mol = this->load_balancer_->molecule(); + + if( not func.is_polarized() ) + GAUXC_GENERIC_EXCEPTION("GKS FXC Contraction Requires A Polarized Functional"); + if( not func.is_lda() ) + GAUXC_GENERIC_EXCEPTION("GKS FXC Contraction Only Implemented For LDA"); + + auto& lb_state = this->load_balancer_->state(); + if( not lb_state.modified_weights_are_stored ) + GAUXC_GENERIC_EXCEPTION("Weights Have Not Been Modified"); + + BasisSetMap basis_map(basis, mol); + const int32_t nbf = basis.nbf(); + + for( int k = 0; k < 4; ++k ) + for( int32_t j = 0; j < nbf; ++j ) + for( int32_t i = 0; i < nbf; ++i ) + FXC[k][i + j*ldFXC[k]] = 0.; + + const size_t ntasks = std::distance(task_begin, task_end); + + #pragma omp parallel + { + std::vector basis_eval, nbe_scr, X, tX, Z, fields, tfields, + npm, vrho, v2rho2; + + #pragma omp for schedule(dynamic) + for( size_t iT = 0; iT < ntasks; ++iT ) { + const auto& task = *(task_begin + iT); + const int32_t npts = static_cast(task.points.size()); + const int32_t nbe = task.bfn_screening.nbe; + const int32_t nshells = static_cast(task.bfn_screening.shell_list.size()); + const auto* points = task.points.data()->data(); + const auto* weights = task.weights.data(); + const int32_t* shell_list = task.bfn_screening.shell_list.data(); + const size_t blk = size_t(npts)*nbe; + + basis_eval.resize( blk ); nbe_scr.resize( size_t(nbe)*nbe ); + X.resize( 4*blk ); tX.resize( 4*blk ); Z.resize( 4*blk ); + fields.resize( 4*size_t(npts) ); tfields.resize( 4*size_t(npts) ); + npm.resize( 2*size_t(npts) ); vrho.resize( 2*size_t(npts) ); v2rho2.resize( 3*size_t(npts) ); + + std::vector< std::array > submat_map; + std::tie(submat_map, std::ignore) = + gen_compressed_submat_map(basis_map, task.bfn_screening.shell_list, nbf, nbf); + + lwd->eval_collocation( npts, nshells, nbe, points, basis, shell_list, basis_eval.data() ); + + // fields and trial fields, slots (s, z, y, x); GKS uses xmat_fac = 1 + for( int k = 0; k < 4; ++k ) { + lwd->eval_xmat( npts, nbf, nbe, submat_map, 1.0, P[k], ldP[k], basis_eval.data(), nbe, + X.data() + k*blk, nbe, nbe_scr.data() ); + lwd->eval_xmat( npts, nbf, nbe, submat_map, 1.0, tP[k], ldtP[k], basis_eval.data(), nbe, + tX.data() + k*blk, nbe, nbe_scr.data() ); + for( int32_t ip = 0; ip < npts; ++ip ) { + const auto* b = basis_eval.data() + size_t(ip)*nbe; + fields [k*npts + ip] = blas::dot( nbe, b, 1, X .data() + k*blk + size_t(ip)*nbe, 1 ); + tfields[k*npts + ip] = blas::dot( nbe, b, 1, tX.data() + k*blk + size_t(ip)*nbe, 1 ); + } + } + + // the locally collinear variables at the TRUE |m| (not the vxc fallback's + // component average): the limit kernel below needs n+ = n- at m = 0 + for( int32_t ip = 0; ip < npts; ++ip ) { + const double mz = fields[1*npts+ip], my = fields[2*npts+ip], mx = fields[3*npts+ip]; + const double mn = std::sqrt( mx*mx + my*my + mz*mz ); + npm[2*ip] = 0.5*( fields[ip] + mn ); + npm[2*ip+1] = 0.5*( fields[ip] - mn ); + } + func.eval_vxc_fxc( npts, npm.data(), vrho.data(), v2rho2.data() ); + + // kernel per point, in the potential's field slots; Z_X = 1/2 k_X chi + for( int32_t ip = 0; ip < npts; ++ip ) { + const double rs = fields[ip], rz = fields[npts+ip], ry = fields[2*npts+ip], rx = fields[3*npts+ip]; + const double ts = tfields[ip], tz = tfields[npts+ip], ty = tfields[2*npts+ip], tx = tfields[3*npts+ip]; + const double v0 = vrho[2*ip], v1 = vrho[2*ip+1]; + const double f0 = v2rho2[3*ip], f1 = v2rho2[3*ip+1], f2 = v2rho2[3*ip+2]; + double ks, kx, ky, kz; + if( std::sqrt( rx*rx + ry*ry + rz*rz ) > dtol ) + xckernel::nc_fxc_contract_lda( rs, rx, ry, rz, v0, v1, f0, f1, f2, 1.0, + ts, tx, ty, tz, ks, kx, ky, kz ); + else + xckernel::nc_fxc_contract_limit_lda( rs, v0, v1, f0, f1, f2, + ts, tx, ty, tz, ks, kx, ky, kz ); + const double kslot[4] = { ks, kz, ky, kx }; + for( int k = 0; k < 4; ++k ) { + const double c = 0.5 * weights[ip] * kslot[k]; + const auto* b = basis_eval.data() + size_t(ip)*nbe; + auto* z = Z.data() + k*blk + size_t(ip)*nbe; + for( int32_t mu = 0; mu < nbe; ++mu ) z[mu] = c * b[mu]; + } + } + + for( int k = 0; k < 4; ++k ) + lwd->inc_vxc( npts, nbf, nbe, basis_eval.data(), submat_map, Z.data() + k*blk, nbe, + FXC[k], ldFXC[k], nbe_scr.data() ); + } // tasks + } // omp parallel + + for( int k = 0; k < 4; ++k ) + for( int32_t j = 0; j < nbf; ++j ) + for( int32_t i = j+1; i < nbf; ++i ) + FXC[k][ j + i*ldFXC[k] ] = FXC[k][ i + j*ldFXC[k] ]; +} + +} // namespace GauXC::detail diff --git a/src/xc_integrator/replicated/replicated_xc_integrator_impl.cxx b/src/xc_integrator/replicated/replicated_xc_integrator_impl.cxx index 071afe312..92d2f937a 100644 --- a/src/xc_integrator/replicated/replicated_xc_integrator_impl.cxx +++ b/src/xc_integrator/replicated/replicated_xc_integrator_impl.cxx @@ -192,6 +192,24 @@ eval_fxc_contraction( int64_t m, int64_t n, const value_type* Ps, } +template +void ReplicatedXCIntegratorImpl:: +eval_fxc_contraction( 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, + const value_type* tPs, int64_t ldtps, const value_type* tPz, int64_t ldtpz, + const value_type* tPy, int64_t ldtpy, const value_type* tPx, int64_t ldtpx, + value_type* FXCs, int64_t ldfxcs, value_type* FXCz, int64_t ldfxcz, + value_type* FXCy, int64_t ldfxcy, value_type* FXCx, int64_t ldfxcx, + const IntegratorSettingsXC& ks_settings ) { + + eval_fxc_contraction_( m, n, Ps, ldps, Pz, ldpz, Py, ldpy, Px, ldpx, + tPs, ldtps, tPz, ldtpz, tPy, ldtpy, tPx, ldtpx, + FXCs, ldfxcs, FXCz, ldfxcz, FXCy, ldfxcy, FXCx, ldfxcx, + ks_settings ); + +} + template void ReplicatedXCIntegratorImpl:: eval_dd_psi( int64_t m, int64_t n, const value_type* P,