Finite-Element Anisotropic Poisson Solver (FEAPS)

• Finite-Element Anisotropic Poisson Solver (FEAPS) is a method that models dielectric responses via tensor-valued finite elements in anisotropic media.
• It employs a variational formulation coupled with CP2K and FEniCSx to accurately compute electrostatic potentials and forces at solid–liquid interfaces.
• The solver achieves scalable parallel performance using MPI and integrates mesh and boundary adaptations to handle complex anisotropic effects.

The Finite-Element Anisotropic Poisson Solver (FEAPS) is a parallel finite-element solver for anisotropic Poisson equations developed in connection with the anisotropic interfacial continuum solvation (AICS) model and implemented through the FEniCSx platform with an interface to CP2K (Chai et al., 7 Aug 2025). In this setting, FEAPS addresses electrostatics in media where the dielectric response is tensorial and spatially varying, particularly at solid–liquid interfaces where in-plane and out-of-plane dielectric components differ and vary along the surface normal. More broadly, FEAPS belongs to the class of finite-element anisotropic Poisson solvers for problems of the form $-\nabla \cdot (K(x)\nabla u)=f$, where anisotropy may arise from heterogeneous diffusion tensors, singular domain geometry, or interface-dependent dielectric structure (Chen, 2018, Li et al., 2010).

1. Definition and problem class

In the AICS formulation, the electrostatic problem solved by FEAPS is a generalized Poisson equation with a dielectric tensor field,

$$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$

where $\phi^{\mathrm{tot}}(\mathbf{r})$ is the total electrostatic potential, $\rho^{\mathrm{solute}}(\mathbf{r})$ is the solute charge density, and $\boldsymbol{\epsilon}(\mathbf{r})$ is a spatially varying dielectric tensor (Chai et al., 7 Aug 2025). In the reported implementation, the tensor is diagonal,

$$\boldsymbol{\epsilon}(\mathbf{r}) = \begin{pmatrix} \epsilon_{xx}(\mathbf{r}) & 0 & 0\ 0 & \epsilon_{yy}(\mathbf{r}) & 0\ 0 & 0 & \epsilon_{zz}(\mathbf{r}) \end{pmatrix},$$

and the off-diagonal components are set to zero (Chai et al., 7 Aug 2025).

This formulation generalizes the classical Poisson problem $-\Delta u=f$ and the standard anisotropic diffusion equation

$-\nabla \cdot (\mathbf{A}\nabla u)=f,$

with $\mathbf{A}(x,y)$ or $\mathbb{D}(\mathbf{x})$ symmetric positive definite and possibly spatially varying (Sadaka, 2012, Li, 2017). In the finite-element literature, this broader problem class includes heterogeneous anisotropic diffusion, porous-medium-type operators with anisotropic permeability, and Poisson equations on polyhedral domains with edge and vertex singularities (Li et al., 2010, Li, 2016).

FEAPS is therefore specific in implementation but not isolated in mathematical type. This suggests that it should be read simultaneously as a solver for AICS electrostatics and as a concrete realization of a more general finite-element anisotropic Poisson methodology.

2. Variational formulation and finite-element structure

The weak form used in FEAPS is written with a test function $$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$0 as

$$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$1

which is assembled into a sparse linear system (Chai et al., 7 Aug 2025). This is the direct tensor-valued analogue of the standard bilinear form

$$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$2

used for anisotropic diffusion in finite-element software such as FreeFem++, where the expanded form includes terms such as $$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$3, $$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$4, and $$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$5 (Sadaka, 2012).

A central finite-element consequence of anisotropy is that the local stiffness contribution depends on the diffusion tensor. In the MATLAB-oriented formulation for generalized Poisson problems,

$$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$6

and the local stiffness entries become

$$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$7

so anisotropy is incorporated by evaluating $$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$8 at quadrature points and inserting it into the gradient contraction (Chen, 2018). The same source emphasizes that “by modifying the subroutine localstiffness, one can easily adapt to new elements and new equations,” which directly clarifies how generalized anisotropic Poisson operators fit into standard FEM assembly pipelines (Chen, 2018).

Within AICS, the dielectric tensor is continuously differentiable with respect to both electron density and spatial coordinates, and analytical expressions were derived for electrostatic contributions to the Kohn–Sham potential and forces (Chai et al., 7 Aug 2025). The electrostatic energy is written as

$$\nabla \cdot \left( \boldsymbol{\epsilon}(\mathbf{r}) \nabla \phi^{\mathrm{tot}}(\mathbf{r}) \right) = -4\pi \rho^{\mathrm{solute}}(\mathbf{r}),$$9

and the contribution to the KS potential includes dielectric-derivative terms,

$\phi^{\mathrm{tot}}(\mathbf{r})$0

with analytical force expressions given in the paper’s equations (15)–(23) (Chai et al., 7 Aug 2025).

3. Software architecture and parallel implementation

FEAPS is implemented on top of FEniCSx / DOLFINx and coupled to CP2K through a Fortran-to-C-to-Python interface (Chai et al., 7 Aug 2025). The solver is designed for distributed-memory execution via MPI, and the simulation cell is subdivided into a uniform 3D hexahedral mesh with linear Lagrange (CG1) basis functions defined at vertices (Chai et al., 7 Aug 2025). The mesh topology, shape, and subdivision mirror CP2K’s real-space grid for consistent data mapping (Chai et al., 7 Aug 2025).

Periodic boundary conditions are imposed in all three Cartesian directions using mesh techniques identical with those used in CP2K wavefunction grids (Chai et al., 7 Aug 2025). To realize periodicity in the finite-element setting, the outermost mesh layer is assigned to process 0, inner layers are distributed across processes, and specialized routines collapse periodic vertices on opposite faces, edges, and corners to unique indices (Chai et al., 7 Aug 2025).

The linear algebra layer uses PETSc with the conjugate gradient Krylov method; the default preconditioner is Geometric-Algebraic Multigrid (GAMG), with JACOBI as fallback (Chai et al., 7 Aug 2025). The reported tolerances are relative $\phi^{\mathrm{tot}}(\mathbf{r})$1, absolute $\phi^{\mathrm{tot}}(\mathbf{r})$2, and maximum iterations typically $\phi^{\mathrm{tot}}(\mathbf{r})$3 (Chai et al., 7 Aug 2025). Data exchange between CP2K and FEAPS is performed by direct memory pointer exchange for electron density, charge density, and dielectric tensor fields, with the stated purpose of minimizing overhead (Chai et al., 7 Aug 2025).

The self-consistent workflow is correspondingly explicit. At each SCF step, CP2K computes a new charge density including electron density and pseudocharge on its grid; FEAPS interpolates that data and the spatially varying dielectric tensor onto the finite-element mesh; the anisotropic Poisson equation is solved; and the resulting potential is mapped back to the CP2K grid for energy, force, and subsequent SCF updates (Chai et al., 7 Aug 2025). Mesh and boundary structures are initialized at the first SCF step, whereas later iterations update only the field data unless remeshing becomes necessary (Chai et al., 7 Aug 2025).

This implementation strategy closely matches a recurring pattern in finite-element PDE platforms. FreeFem++, for example, was developed at the Laboratoire Jacques-Louis Lions, Université Pierre et Marie Curie, Paris by Frédéric Hecht in collaboration with Olivier Pironneau, Jacques Morice, Antoine Le Hyaric and Kohji Ohtsuka as an open-source platform for solving PDEs numerically through finite-element methods, with high-level variational specification and interfaces to direct sparse solvers such as UMFPACK and SuperLU (Sadaka, 2012). FEAPS differs in application target and parallel coupling, but the underlying design principle—direct coding of variational forms combined with external sparse linear algebra—is consistent.

4. Anisotropy, meshing, and numerical analysis context

Anisotropic Poisson solvers are not defined only by tensor coefficients; they are also shaped by mesh design. In the finite-element literature, “metric-based anisotropic mesh adaptation” is a standard mechanism for aligning element size and orientation with anisotropic physics or anisotropic solution features (Sadaka, 2012). For the anisotropic porous medium equation, mesh adaptation is formulated through a metric tensor $\phi^{\mathrm{tot}}(\mathbf{r})$4 under equidistribution and alignment conditions, with diffusion-aligned, solution-adaptive, and combined metrics denoted $\phi^{\mathrm{tot}}(\mathbf{r})$5, $\phi^{\mathrm{tot}}(\mathbf{r})$6, and $\phi^{\mathrm{tot}}(\mathbf{r})$7 (Li, 2017). In heterogeneous anisotropic diffusion, the DMP-oriented metric

$\phi^{\mathrm{tot}}(\mathbf{r})$8

and the combined DMP-plus-adaptivity metric

$\phi^{\mathrm{tot}}(\mathbf{r})$9

were derived to enforce the anisotropic non-obtuse angle condition while adapting to solution variation (Li et al., 2010).

A distinct but related line of work studies anisotropy caused by singular geometry rather than by coefficient tensors alone. For the Poisson equation on three-dimensional polyhedral domains, anisotropic tetrahedral refinement algorithms classify tetrahedra by proximity to singular vertices and edges and use grading parameters $\rho^{\mathrm{solute}}(\mathbf{r})$0 to compress the mesh toward singular sets (Li, 2016). In the weighted-space analysis developed there, if the solution lies in the appropriate anisotropic weighted space, then after $\rho^{\mathrm{solute}}(\mathbf{r})$1 refinements with $\rho^{\mathrm{solute}}(\mathbf{r})$2,

$\rho^{\mathrm{solute}}(\mathbf{r})$3

and for the finite-element solution $\rho^{\mathrm{solute}}(\mathbf{r})$4,

$\rho^{\mathrm{solute}}(\mathbf{r})$5

yielding optimal convergence with respect to degrees of freedom when grading matches singularity structure (Li, 2016).

Another misconception in anisotropic FEM is that meaningful error estimates require shape-regular or maximum-angle meshes. For three-dimensional anisotropic meshes, error analyses for piecewise linear nonconforming Crouzeix–Raviart and lowest-order Raviart–Thomas methods explicitly avoid both the shape-regularity condition and the maximum-angle condition, instead using the anisotropy parameter

$\rho^{\mathrm{solute}}(\mathbf{r})$6

with the assumption $\rho^{\mathrm{solute}}(\mathbf{r})$7 (Ishizaka et al., 2020). Under that assumption, the paper proves, for example,

$\rho^{\mathrm{solute}}(\mathbf{r})$8

and

$\rho^{\mathrm{solute}}(\mathbf{r})$9

on meshes that may contain very flat or sliver tetrahedra (Ishizaka et al., 2020).

These results do not describe FEAPS directly, because the FEAPS implementation reported in AICS uses a uniform 3D hexahedral mesh with CG1 basis functions (Chai et al., 7 Aug 2025). They do, however, place FEAPS in a larger theoretical landscape where anisotropy may be handled through coefficient tensors, mesh metrics, or singularity-adapted grading, and where angle restrictions are not universally necessary.

5. Validation, numerical behavior, and computational performance

The reported FEAPS validation covers forces, electrostatic potentials, and parallel timing (Chai et al., 7 Aug 2025). Analytical forces, including dielectric derivative contributions, were compared with finite-difference calculations for a five-layer Ag(111) slab in implicit anisotropic water, yielding a maximum absolute force difference of $\boldsymbol{\epsilon}(\mathbf{r})$0 Hartree/Bohr, reported as less than $\boldsymbol{\epsilon}(\mathbf{r})$1 (Chai et al., 7 Aug 2025). This establishes agreement between the analytical derivatives and numerical reference forces within the tested setting.

For electrostatic potentials under vacuum conditions, FEAPS with $\boldsymbol{\epsilon}(\mathbf{r})$2 was benchmarked against the conventional CP2K FFT-based Poisson solver (Chai et al., 7 Aug 2025). The reported agreement was within $\boldsymbol{\epsilon}(\mathbf{r})$3 eV for work functions, and the mean absolute difference in planar or grid-based comparison was $\boldsymbol{\epsilon}(\mathbf{r})$4 Hartree (Chai et al., 7 Aug 2025). The cited explanation for residual discrepancies is the difference in functional representation: near nuclei, the finite-element solution uses linear interpolation per element, whereas the FFT solver assumes a uniform value per grid cell (Chai et al., 7 Aug 2025). Under isotropic, uniform dielectric conditions with $\boldsymbol{\epsilon}(\mathbf{r})$5, FEAPS and the SCCS model agreed within numerical noise, again quantified as $\boldsymbol{\epsilon}(\mathbf{r})$6 eV in work function comparison (Chai et al., 7 Aug 2025).

For parallel performance, the paper reports that on a single node with 288 cores, the wall time for a full SCF step including FEAPS can be drastically reduced by increasing the number of MPI processes up to 12, beyond which speedup saturates because mesh generation is less parallel (Chai et al., 7 Aug 2025). The GAMG preconditioner is identified as the most efficient in wall-clock time, whereas at high cutoff or mesh densities memory becomes the limiting factor (Chai et al., 7 Aug 2025). The stated communication overhead is minimal because of direct pointer exchange between CP2K and FEAPS (Chai et al., 7 Aug 2025).

The broader anisotropic-mesh literature provides an important counterpoint concerning conditioning. In adaptive finite elements with anisotropic meshes, a common concern is that stretched elements necessarily make the linear system impractically ill-conditioned. The reported numerical study on a Poisson test problem with a corner singularity, a peak, a boundary layer, and a wavefront found that the condition number of the anisotropic system was about one order of magnitude worse than the isotropic system, but also concluded that “the conditioning with adaptive anisotropic meshes is not as bad as generally assumed” (Huang et al., 2012). When anisotropy is local, diagonal scaling can reduce the condition number to values comparable to isotropic or quasi-uniform cases (Huang et al., 2012). This suggests that FEAPS-like solvers should distinguish between conditioning growth and practical solver failure rather than treating anisotropy itself as pathological.

6. Applications, significance, and limits

The immediate application of FEAPS is the AICS model for solid–liquid interfaces, where the dielectric response near surfaces is not adequately represented by a single isotropic scalar (Chai et al., 7 Aug 2025). In low-electron-density regions, each dielectric function in the diagonal tensor components varies monotonically with distance from the solid surface along the surface normal; in high-electron-density regions near the surface, each dielectric function adopts the electron-density-based formulation proposed by Andreussi et al. (Chai et al., 7 Aug 2025). Within this framework, FEAPS was used to compute work functions, electrostatic potentials, and adsorption geometries for OH on Ag(111), and the anisotropic solvent environment was characterized by enhanced in-plane and reduced out-of-plane dielectric functions near the surface (Chai et al., 7 Aug 2025). Compared with the isotropic case, the reported outcomes included more pronounced work-function shifts, spatially modulated electrostatic profiles across different charge states, and an OH configuration tilted more toward the plane parallel to the surface under anisotropic dielectric conditions (Chai et al., 7 Aug 2025).

A broader implication is that FEAPS exemplifies the use of finite elements where tensor-valued, spatially varying coefficients make conventional isotropic Poisson solvers difficult or inapplicable (Chai et al., 7 Aug 2025). Earlier diffusion-oriented work had already shown that standard numerical methods may produce spurious oscillations for heterogeneous anisotropic diffusion problems and that mesh conditions tied to the diffusion tensor, including the anisotropic non-obtuse angle condition, are relevant for preserving the discrete maximum principle (Li et al., 2010). Likewise, recent work on Vlasov–Poisson discretization argues that isotropic stabilization may fail on stretched or anisotropic meshes and introduces direction-dependent artificial viscosity precisely because anisotropic resolution changes the stability requirements of the coupled Poisson setting (Wen et al., 10 Mar 2025). This suggests that FEAPS should be understood not merely as a replacement linear solver, but as part of a wider methodological shift toward tensor-aware discretization and anisotropy-aware numerical control.

At the same time, the reported FEAPS implementation has explicit limits. The dielectric tensor is diagonal in the current implementation, so only principal-axis anisotropy is represented directly (Chai et al., 7 Aug 2025). The boundary treatment is periodic in all three Cartesian directions, reflecting the intended first-principles simulation setting rather than a universal boundary-condition framework (Chai et al., 7 Aug 2025). The speedup saturation beyond 12 MPI processes on a single 288-core node and the onset of memory limits at high cutoff or mesh density indicate that scalability is not determined by the Krylov solve alone (Chai et al., 7 Aug 2025). More generally, the anisotropic finite-element literature shows that robustness may depend on metric quality, grading parameters, or mesh anisotropy parameters such as $\boldsymbol{\epsilon}(\mathbf{r})$7, even when classical angle conditions are dropped (Li, 2016, Ishizaka et al., 2020).

FEAPS thus occupies a specific and technically important position: it is a parallel finite-element realization of anisotropic Poisson electrostatics for continuum solvation and electronic-structure coupling, while also fitting into a longer research trajectory on anisotropic diffusion, mesh adaptation, discrete maximum principles, singularity resolution, and solver robustness in finite-element PDE computation (Chai et al., 7 Aug 2025, Li et al., 2010, Huang et al., 2012).