Procedures

@brief Compute A = A - A^T of a square matrix @param[inout] a square NxN matrix @param[in] n matrix dimension

@brief Compute overlap matrix between two basis sets @details Overlap integrals are computed using Gauss-Hermite quadrature formula @author Igor S. Gerasimov @date Oct, 2022 Initial release

@brief Convert raw null-terminated C-string to Fortran string

@brief Compute primitive block of multipole integrals up to an order MXMOM @param[in] cp shell pair data @param[in] id current pair of primitives @param[in] r point in space to compute integrals @param[in] mxmom multiplole moment order (1-dipole, 2-quadrupole, 3-octopole) @param[inout] blk block of 1e multipole moment integrals @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Subtract damping function term from ESP block @details Compute one-electron Coulomb integrals with the damping function: \f$ |r-r_C|^{-1} (1 - \beta e^{-\alpha(r-r_C)^2}) \f$ Only the part \f$ - |r-r_C|^{-1} \beta e^{-\alpha(r-r_C)^2}) \f$ is computed here; the other part is regular Coulomb potential computed elsewhere @param[in] cp shell pair data @param[in] id current pair of primitives @param[in] alpha dumping exponent @param[in] beta dumping function scaling factor @param[in] c coordinates of the charged particle @param[in] znuc particle charge @param[inout] vblk block of 1e Coulomb integrals @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute 1e Coulomb contribution to the gradient (v.r.t. shifts of shell's centers) @param[in] nroots roots for GaussRys @param[in] cp shell pair data @param[in] c coordinates of the charged particle @param[in] znuc particle charge @param[in] dij density matrix block @param[inout] dernuc dimension(3), contribution to gradient @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute 1e Hellmann-Feynman contribution to the gradient @param[in] nroots roots for GaussRys @param[in] cp shell pair data @param[in] c coordinates of the charged particle @param[in] znuc particle charge @param[in] dij density matrix block @param[inout] derhf dimension(3), contribution to gradient @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute primitive block of 1e Coulomb atraction integrals @param[in] cp shell pair data @param[in] id current pair of primitives @param[in] c coordinates of the charged particle @param[in] znuc particle charge @param[inout] vblk block of 1e Coulomb integrals @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute sum of 1e Coulomb integrals over primitive shell pair @param[in] cp shell pair data @param[in] id current pair of primitives @param[in] c coordinates of the charged particle @param[in] den normalized density matrix block @param[inout] vsum sum of Coulomb integrals over pair of primitives @author Vladimir Mironov @date Oct, 2018 Initial release

This routine calculates the determinate of a square matrix. Gauss Elimination Method array the matrix of order norder which is to be evaluated. this subprogram destroys the matrix array norder the order of the square matrix to be evaluated.

@brief Compute 1e Ewald long-range contribution to the gradient (v.r.t. shifts of shell's centers) @param[in] nroots roots for GaussRys @param[in] cp shell pair data @param[in] c coordinates of the charged particle @param[in] znuc particle charge @param[in] dij density matrix block @param[in] omega Ewald splitting parameter @param[inout] dernuc dimension(3), contribution to gradient @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute Ewald long-range 1e Hellmann-Feynman contribution to the gradient @param[in] nroots roots for GaussRys @param[in] cp shell pair data @param[in] c coordinates of the charged particle @param[in] znuc particle charge @param[in] dij density matrix block @param[in] omega Ewald splitting parameter @param[inout] derhf dimension(3), contribution to gradient @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute primitive block of 1e Coulomb atraction integrals for Ewald summation, long-range part @details 1e integrals using modified Coulomb potential: \f$ \frac{Erf(\omega^{1/2}|r-r_C|)}{|r-r_C|} \f$ @param[in] cp shell pair data @param[in] id current pair of primitives @param[in] c coordinates of the charged particle @param[in] znuc particle charge @param[in] omega Ewald splitting parameter @param[inout] vblk block of 1e Coulomb integrals @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute primitive block of overlap and kinetic energy 1e integrals @param[in] cp shell pair data @param[in] id current pair of primitives @param[in] dokinetic if .FALSE. compute only overlap integrals @param[inout] sblk block of 1e overlap integrals @param[inout] tblk block of 1e kinetic energy integrals @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute 1e kinetic contribution to the gradient @param[in] cp shell pair data @param[in] dij density matrix block @param[inout] de dimension(3), contribution to gradient @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute primitive block of 1e Coulomb ESP integrals in FMO method @param[in] cp shell pair data @param[in] id current pair of primitives @param[inout] zblk block of 1e Lz-integrals @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute primitive block of multipole integrals of order MOM @param[in] cp shell pair data @param[in] id current pair of primitives @param[in] r point in space to compute integrals @param[in] mom multiplole moment order (1-dipole, 2-quadrupole, 3-octopole) @param[inout] blk block of 1e multipole moment integrals @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute 1e overlap contribution to the gradient @param[in] cp shell pair data @param[in] dij density matrix block @param[inout] de dimension(3), contribution to gradient @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute XC contributions to the gradient from a grid point, GGA part @param[inout] bfGrad array of gradient contributions per basis function @param[in] fgrad XC gradient @param[in] moV MO-like orbital values @param[in] moG1 MO-like orbital gradients @author Vladimir Mironov

@brief Compute XC contributions to the gradient from a grid point, LDA part @param[in] iPt index of a grid point @param[inout] bfGrad array of gradient contributions per basis function @param[in] fgrad XC gradient @param[in] moV MO-like orbital values @author Vladimir Mironov

@brief Compute XC contributions to the gradient from a grid point, mGGA part @param[in] iPt index of a grid point @param[inout] bfGrad array of gradient contributinos per basis function @param[in] dedta XC energy, mGGA contribution, alpha-spin @param[in] dedtb XC energy, mGGA contribution, beta-spin @author Vladimir Mironov

@brief Compute overlap integrals between MRSF response states at different MD time steps

@brief Corresponding orbital projection @param[in] nproj number of orbitals from vb to projected onto @param[in] l0 number of the first orbitals from the va space to project vb to @param[in,out] va a orbitals on entry, replaced by corresponding orbitals on exit @param[in] vb b orbitals on entry, unchanged on exit @param[in] sba overlap integrals between the two bases,

@brief This function return count of substring in string @author Igor S. Gerasimov @date Sep, 2019 --Initial release-- @date May, 2021 Moved to strings @param substring - (in) @param string - (in)

@brief Copy density block from the triangular density matrix @details This subroutine assumes shi%shid>=shj%shid @param[in] shi first shell data @param[in] shj second shell data @param[in] dij density matrix in packed triangular form @param[out] denab density matrix block for shells shi and shj

@brief Copy density block from the triangular density matrix @details This subroutine assumes arbitrary order of shell IDs @param[in] shi first shell data @param[in] shj second shell data @param[in] dij density matrix in packed triangular form @param[out] denab density matrix block for shells shi and shj

@brief Compute grid XC contribution to the nuclear gradient

@brief Assemble numerical atomic DFT grids to a molecular grid @param[in] atmxvec array of atomic X coordinates @param[in] atmyvec array of atomic Y coordinates @param[in] atmzvec array of atomic Z coordinates @param[in] rij interatomic distances @param[in] nat number of atoms @param[in] curAt index of current atom @param[in] rad effective (e.g. Bragg-Slater) radius of current atom @param[inout] wtab normalized cell function values for LRD @param[in] aij surface shifting factors for Becke's method @author Vladimir Mironov

@brief Analytical DFT gradient

@brief Find eigenvalues and eigenvectors of symmetric matrix in full format @param[in] mode algorithm of diagonalization (not used now) @param[in] n matrix dimension @param[in,out] a matrix to be diagonalized, overwritten by the eigenvectors on the exit @param[in] lda leading dimension of the matrix @param[out] eig eigenvalues @param[out] ierr status

@brief Find eigenvalues and eigenvectors of symmetric matrix in packed format @param[in] mode algorithm of diagonalization (not used now) @param[in] n matrix dimension @param[in] ldvect leading dimension of eigenvector matrix @param[in] nvect required number of eigenvectors @param[in,out] h matrix to be diagonalized @param[out] eig eigenvalues @param[out] vector eigenvectors @param[out] ierr status

@brief Compute grid XC contribution to the Kohn-Sham matrix @param[in] coeffa MO coefficients, alpha-spin @param[in] coeffb MO coefficients, beta-spin @param[inout] fa KS matrix, alpha-spin @param[inout] fb KS matrix, beta-spin @param[out] exc XC energy @param[out] totele electronic denisty integral @param[out] totkin kinetic energy integral @param[in] mxAngMom max. needed ang. mom. value (incl. derivatives) @param[in] nbf basis set size @param[in] urohf .TRUE. if open-shell calculation @author Vladimir Mironov

@brief Gauss-Hermite quadrature using minimum point formula @details Compute: xint = sum( w(1:npts,npts) * (h(1:npts,npts)t+dxi)(ni-1) * (h(1:npts,npts)t+dxj)(nj-1) ) yint = sum( w(1:npts,npts) * (h(1:npts,npts)*t+dyi)(ni-1) * (h(1:npts,npts)t+dyj)(nj-1) ) zint = sum( w(1:npts,npts) * (h(1:npts,npts)t+dzi)(ni-1) * (h(1:npts,npts)*t+dzj)(nj-1) )

@brief Compute "energy weighted density matrix"

@brief Compute electronic contribution to electrostatic potential on a grid @author Vladimir Mironov @date Sep, 2023 Initial release

@breif Add (E_a-E_i)*Z_ai to Pmo

@brief Convert Fortran string to a raw null-terminated C-string

@brief Form the ROHF Fock matrix in MO basis.

@brief This function return fortran allocatable string @author Igor S. Gerasimov @date April, 2022 --Initial release-- @param string - (in) C-like string

@brief this will calculatte the density matrix

@brief Solve \f$ F C = \eps S C \f$

@brief Compute derivative coupling vectors (DCV) using the finite difference method, typically denoted as d_IJ between states I and J.

@brief Set the radial points and weights @param ptrad [inout] quadrature points @param wtrad [inout] quadrature weights @param rad_grid_type [in] type of radial grid @param rad_grid_type [in] which quadrature to map on the selected grid

@brief C-interoperable wrapper for get_states_overlap

@brief Main subroutine for calculating state overlaps and derivative coupling matrix elements

@brief Calculate AO overlap between two different geometries @param[in,out] infos Information structure containing all necessary data

@brief C-interoperable wrapper for get_states_overlap

@brief Effective core potential gradient

@brief Basis function derivative contributions to gradient @details Compute derivative integrals of type =

@brief Compute overlap energy derivative contribution to gradient @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Hellmann-Feynman force @details Compute derivative contributions due to the Hamiltonian operator change w.r.t. shifts of nuclei. The contribution of the form is evaluated by Gauss-Rys quadrature. This version handles spdfg and L shells.

@brief Basis function derivative contributions to gradient @details Compute derivative integrals of type =

@brief Gradient of nuclear repulsion energy

@brief The driver for the two electron gradient

@brief Compute density factors for \Gamma term of 2-electron contribution to TD-DFT gradients

@brief This function return index of i'th substring in string @details if ind is negative, backward search will be result is equal 0 if search was failed @author Igor S. Gerasimov @date May, 2021 --Initial release-- @param substring - (in) @param string - (in) @param ith - (in) index of needed substring

@brief Calculate the basic H, S, and T 1e-integrals

@brief Compute H^-[V] over a set of integrals

@brief Compute H^+[V] over a set of integrals

@brief Destroy internal variables of functional @author Igor S. Gerasimov @date Dec, 2020 - Initial release - @date Dec, 2022 Pass local functional instead of global

@brief setting up of using libxc functionals (Analog of INPGDFT) @detail @author Igor S. Gerasimov @date July, 2019 - Initial release - @date July, 2021 Using messages module Adding optional arguments @date Dec, 2022 Pass functional instead of using global @params functional_name (in) functional's name @param infos (inout) info datatype @param functional (inout) constructing functional

@brief Compute matrix inverse square root using SVD and removing linear dependency @detail This subroutine is used to obtain set of canonical orbitals by diagonalization of the basis set overlap matrix Q = S^{-1/2}, Q^T * S * Q = I @param[in] s Overlap matrix, symmetric packed format @param[out] q Matrix inverse square root, square matrix @param[in] nbf Dimeension of matrices S and Q, basis set size @param[out] qrnk Rank of matrix Q @param[in] tol optional, tolerance to remove linear dependency, default = 1.0e-8

@brief Generate transformation to spherical harmonics basis

@brief Calculate Molecular Orbital (MO) overlap between two geometries

@brief Back-transform a symmetric operator Fmo expressed in the MO basis V to the AO basis @detail compute the transformation: Fao = S*V * Fmo * (SV)^T

@brief Compute overlap integrals between CSFs of MRSF-TDDFT using TLF approximation

@brief Spin-pairing parts of singlet and triplet MRSF Lagrangian

@details This subroutine transforms Multi-Reference Spin-Flip (MRSF) response vectors from a compressed representation to an expanded form. It handles both singlet (mrst=1) and triplet (mrst=3) cases.

@brief Gauss-Hermite quadrature using minimum point formula @details Compute: xint = sum( w(1:npts,npts) * (h(1:npts,npts)t+dxi)(ni-1) * (h(1:npts,npts)t+dxj)(nj-1) ) yint = sum( w(1:npts,npts) * (h(1:npts,npts)*t+dyi)(ni-1) * (h(1:npts,npts)t+dyj)(nj-1) ) zint = sum( w(1:npts,npts) * (h(1:npts,npts)t+dzi)(ni-1) * (h(1:npts,npts)*t+dzj)(nj-1) )

@brief Driver for multipole integrals @details Compute one electron multipole integrals Integrals are evaluated by Gauss-Hermite quadrature, @author Vladimir Mironov @date Feb, 2023 Initial release

@brief Driver for conventional h, S, and T integrals @details Compute one electron integrals and core Hamiltonian, - S is evaluated by Gauss-Hermite quadrature, - T is an overlap with -2,0,+2 angular momentum shifts, - V is evaluated by Gauss-Rys quadrature, then \f$ h = T+V \f$ Also, do \f$ L_z \f$ integrals if requested

Compute the inverse of a real matrix using LAPACK routine dgetri.

@brief Compute ESP charges fitted by modified Merz-Kollman method

@brief Compute density matrix from a set of orbitals and respective occupation numbers @detail Compute the transformation: D = V * X * V^T @param[out] d density matrix @param[in] v matrix of orbitals @param[in] x vector of occupation numbers @param[in] m number of columns in V @param[in] n dimension of D @param[in] ldv leading dimension of V

@brief Compute orthogonal transformation of a square marix @param[in] trans If trans='n' compute B = U^T * A * U If trans='t' compute B = U * A * U^T @param[in] ld Dimension of matrices @param[in] u Square orthogonal matrix @param[inout] a Matrix to transform, optionally output matrix @param[out] b Result, can be absent for in-place transform of matrix A @param[inout] wrk Scratch space, optional @author Vladimir Mironov

@brief Compute orthogonal transformation of a square marix @param[in] trans If trans='n' compute B = U^T * A * U If trans='t' compute B = U * A * U^T @param[in] ld Dimension of matrices @param[in] u Square orthogonal matrix @param[in] a Matrix to transform @param[out] b Result @param[inout] wrk Scratch space @author Vladimir Mironov

@brief Compute orthogonal transformation of a symmetric marix A in packed format: B = U^T * A * U @param[in] a Matrix to transform @param[in] u Orthogonal matrix U(ldu,m) @param[in] n dimension of matrix A @param[in] m dimension of matrix B @param[in] ldu leading dimension of matrix U @param[out] b Result @param[inout] wrk Scratch space @author Vladimir Mironov

@brief Compute overlap and integrals @details Overlap integrals are computed using Gauss-Hermite quadrature formula @author Vladimir Mironov @date Mar, 2023 Initial release

@brief Fortran-90 routine for packing symmetric matrix to 1D array

@brief pFON Implementation in SCF Module Author: Alireza Lashkaripour Date: January 2025 Reference paper: https://doi.org/10.1063/1.478177 This subroutine incorporates the Partial Fractional Occupation Number (pFON) method into SCF calculations, ensuring smooth occupation numbers using Fermi-Dirac distribution. It dynamically adjusts temperature and beta factors to enhance SCF convergence, particularly for near-degenerate states.

@brief print eigenvector/values, with MO symmetry labels

@brief Print the solution of the eigenvalue problem @param[in] v matrix of eigenvectors, v(ldv,m) @param[in] e array of eigenvectors, e(m) @param[in] m dimension of column space @param[in] n dimension of row space @param[in] ldv leading dimension of V

@brief Print geometry information

@brief Print energy gradient vector

@brief Printing out MOs

@brief Print MODULE information @detail Printout the information of each MODULES of OQP

@brief Print results of MO overlap between two geometries

@brief Print out a square matrix @param[in] v rectanbular matrix @param[in] m number of columns in V @param[in] n number of rows in V @param[in] ndim leading dimension of V @param[in] tag optional, will be printed at the beginning of each line @param[in] maxcolumns optional, number of columns to wrap printing

@brief Print symmetric matrix D in square format @param[in] d matrix to print, only lower triangle is referenced @param[in] n rank of matrix D @param[in] ld leading dimension of matrix D

@brief Print symmetric packed matrix d of dimension n @detail The rows will be labeled with basis function tags

@brief Print out a symmetric matrix in packed format @param[in] d symmetric matrix in packed format @param[in] n matric dimension

@brief Find shells and primitives which are significant in a given set of 3D coordinates @author Vladimir Mironov

@brief This routine remove not needed spaces from section line @detail This routine used revert reading of lines. Example, that this routine do: > DFTTYP = PBE0 BASNAM=APC4 , ACC5 < DFTTYP=PBE0 BASNAM=APC4,ACC5 @author Igor S. Gerasimov @date Sep, 2019 --Initial release-- @date May, 2021 Moved to strings @param line - (inout) worked line

@brief This routine reorders orbitals to maximum overlap.

@brief Compute ESP charges fitted by modified Merz-Kollman method @details This is a C interface for a Fortran subroutine

@brief Remove negative eigenvalues

@brief Diagonalize small reduced RPA matrix

@brief Expand reduced vectors to real size space

@brief Orthonormalize q(xvec_dim,ndsr*2) and append to bvec

@brief Print current excitation energies and errors

@brief Construct residual vectors and check convergence

@brief Normalize V1 and V2 by biorthogonality condition

@brief Run population analysis

@brief Compute the solution to a real system of linear equations A * X = B @detai Wrapper for DSYSV from Lapack @param[in,out] A general matrix, destroyed on exit. @param[in,out] B RHS on entry. The solution on exit. @param[in] n size of the problem. @param[in] nrhs number of RHS vectors @param[in] lda leading dimension of matrix A @param[out] ierr Error flag, ierr=0 if no errors. Read DSYEV manual for details

@brief Compute A = A + A^T of a square matrix @param[inout] a square NxN matrix @param[in] n matrix dimension

@brief Compute derivative XC contribution to the TD-DFT KS-like matrices @param[in] basis basis set @param[in] isVecs .true. if orbitals are provided instead of density matrix @param[in] wf density matrix/orbitals @param[inout] fx fock-like matrices @param[inout] dx densities @param[in] nMtx number of density/Fock-like matrices @param[in] threshold tolerance @param[in] infos OQP metadata @author Vladimir Mironov

@brief Compute derivative XC contribution to the TD-DFT KS-like matrices @param[in] basis basis set @param[in] isVecs .true. if orbitals are provided instead of density matrix @param[in] wf density matrix/orbitals @param[inout] fx fock-like matrices @param[inout] dx densities @param[in] nMtx number of density/Fock-like matrices @param[in] threshold tolerance @param[in] infos OQP metadata @author Vladimir Mironov

@brief Compute derivative XC contribution to the TD-DFT KS-like matrices @param[in] basis basis set @param[in] wf density matrix/orbitals @param[inout] fx fock-like matrices @param[inout] dx densities @param[in] nMtx number of density/Fock-like matrices @param[in] threshold tolerance @param[in] isGGA .TRUE. if GGA/mGGA functional used @param[in] infos OQP metadata @author Vladimir Mironov

@brief The driver for the two electron gradient

@brief Compute unrelaxed difference density matrix for R-TDDFT

@brief return string in lower case @author Igor S. Gerasimov @date Sep, 2021 --Initial release-- @param line - (inout)

@brief return string in upper case @author Igor S. Gerasimov @date Sep, 2019 --Initial release-- @date May, 2021 Moved to strings @date Sep, 2021 subroutine -> function @param line - (in)

@brief Compute the trace of the product of two symmetric matrices in packed format @detail The trace is actually an inner product of matrices, assuming they are vectors @param[in] a first matrix @param[in] b second matrix @param[in] n dimension of matrices A and B

@brief Fill the upper/lower triangle of the symmetric matrix @param[inout] a square NxN matrix in triangular form @param[in] n matrix dimension @param[in] uplo U if A is upper triangular, L if lower triangular

@brief Fortran-90 routine for unpacking 1D array to symmetric matrix

@brief Add contribution of the 1e-integral block to the rectangular matrix @param[in] shi first shell data @param[in] shj second shell data @param[in] mblk square block of 1e integrals passed as 1D array @param[inout] m rectangular matrix of 1-e integral contribution @author Igor S. Gerasimov @date Oct, 2022 Initial release

@brief Add contribution of the 1e-integral block to the triangular matrix @param[in] shi first shell data @param[in] shj second shell data @param[in] mblk square block of 1e integrals passed as 1D array @param[inout] m packed triangular matrix of 1e integral contribution @author Vladimir Mironov @date Sep, 2018 Initial release

@brief Compute derivative XC contribution to the TD-DFT KS-like matrices @param[in] basis basis set @param[in] isVecs .true. if orbitals are provided instead of density matrix @param[in] wf density matrix/orbitals @param[inout] fx fock-like matrices @param[inout] dx densities @param[in] nMtx number of density/Fock-like matrices @param[in] threshold tolerance @param[in] isGGA .TRUE. if GGA/mGGA functional used @param[in] infos OQP metadata @author Vladimir Mironov

@brief Compute derivative XC contribution to the TD-DFT KS-like matrices @param[in] basis basis set @param[in] wf density matrix/orbitals @param[inout] fx fock-like matrices @param[inout] dx densities @param[in] nMtx number of density/Fock-like matrices @param[in] threshold tolerance @param[in] isGGA .TRUE. if GGA/mGGA functional used @param[in] infos OQP metadata @author Vladimir Mironov

@brief Get first derivative of the XC functional @param[in] xce XC engine @param[in] beta Whether to return spin-polarized quantities @param[out] d_r dE_xc / d_rho (alpha, beta) @param[out] d_s dE_xc / d_sigma (alpha-alpha, beta-beta, alpha-beta) @param[out] d_t dE_xc / d_tau (alpha, beta)

@brief Get second derivative of the XC functional contracted with response densities @param[in] xce XC engine @param[in] beta Whether to return spin-polarized quantities @param[in] d_r \delta_rho (alpha, beta) @param[in] d_s \delta_sigma (alpha-alpha, beta-beta, alpha-beta) @param[in] d_t \delta_tau (alpha, beta) @param[out] f_r \sum_i d2E_xc / (d_rho * d_zeta_i) (alpha, beta) @param[out] f_s \sum_i d2E_xc / (d_sigma * d_zeta_i) (alpha-alpha, beta-beta, alpha-beta) @param[out] f_t \sum_i d2E_xc / (d_tau * d_zeta_i) (alpha, beta)

@brief Get third derivative of the XC functional contracted with response densities, spin-polarized version @param[in] xce XC engine @param[in] d_r \delta_rho (alpha, beta) @param[in] d_s \delta_sigma (alpha-alpha, beta-beta, alpha-beta) @param[in] d_t \delta_tau (alpha, beta) @param[in] ss (\nabla\rho(T)\nabla\rho(T)) (alpha-alpha, beta-beta, alpha-beta) @param[out] g_r \sum_i,j d2E_xc / (d_rho * d_zeta_id_zeta_j) (alpha, beta) @param[out] g_s \sum_i,j d2E_xc / (d_sigma * d_zeta_id_zeta_j) (alpha-alpha, beta-beta, alpha-beta) @param[out] g_t \sum_i,j d2E_xc / (d_tau * d_zeta_id_zeta_j) (alpha, beta)