Showing 1–50 of 56 results for author: Stamm, B

Diagonalization of the Landau Hamiltonian with a Periodic Potential via a Galerkin Projection Method and Applications to Topological Band Properties

Authors: Rafael Antonio Lainez Reyes, Benjamin Stamm, Hans Peter Büchler

Abstract: In this paper, we study the numerical computation of eigenfunctions of the Landau Hamiltonian with a periodic potential. We propose a spectral method using the eigenfunctions of the Landau operator that allows us to numerically compute the band structure for potentials of different strengths, including the well-studied weak and strong potentials, and additionally systems in the intermediate regime… ▽ More In this paper, we study the numerical computation of eigenfunctions of the Landau Hamiltonian with a periodic potential. We propose a spectral method using the eigenfunctions of the Landau operator that allows us to numerically compute the band structure for potentials of different strengths, including the well-studied weak and strong potentials, and additionally systems in the intermediate regime. We then apply the method to simulations across the full range of coupling constants from weak to strong coupling and present an analysis of the band structure as a function of the potential strength. △ Less

Submitted 20 August, 2026; originally announced August 2026.

CP2K: An electronic structure and molecular dynamics software package - Dynamics, Transport, and Spectroscopic Response

Authors: Jan Wilhelm, Anna-Sophia Hehn, Hossam Elgabarty, Beliz Sertcan Gökmen, Maximilian Graml, Stepan Marek, Ritaj Tyagi, Frederick Stein, Johann V. Potoschnig, Augustin Bussy, Christian S. Ahart, Zehua Chen, Linnea Andersson, Zdenek Futera, Filip Ivanovic, Margherita Buraschi, Christoph Schran, Remi Pasquier, Leonard Prokisch, Bibek Samal, Jelena Schmitz, Shridhar Sanjay Shanbhag, Harald Forbert, Ole Schütt, Franz Pöschl , et al. (17 additional authors not shown)

Abstract: One of the distinguishing aspects of CP2K is its seamless integration of diverse structural and transition-state optimization techniques with advanced sampling approaches including Monte Carlo, molecular dynamics, and metadynamics, enabling the efficient exploration of complex potential- and free-energy landscapes, including rare events. These capabilities are combined with a broad hierarchy of en… ▽ More One of the distinguishing aspects of CP2K is its seamless integration of diverse structural and transition-state optimization techniques with advanced sampling approaches including Monte Carlo, molecular dynamics, and metadynamics, enabling the efficient exploration of complex potential- and free-energy landscapes, including rare events. These capabilities are combined with a broad hierarchy of energy and force evaluation methods, ranging from classical and machine-learned interaction potentials and mixed quantum-classical multiscale and semiempirical schemes, to highly accurate quantum-mechanical electronic-structure approaches. At the heart of the latter lies the Gaussian and plane-wave framework, along with its augmented all-electron generalization, which have been described in detail in our previous code review [T. D. Kühne et al., J. Chem. Phys. 152, 194103 (2020)]. Building on this foundation, the present work revisits the methods within CP2K that turn electronic structure into dynamics, transport, and spectroscopic response. Particular emphasis is placed on the coupling between static response calculations and nuclear motion: spectra may be evaluated at optimized structures, averaged over thermally sampled configurations, obtained from time-correlation functions along ab-initio or path integral molecular trajectories, or followed in real time together with electronic and nuclear dynamics. The same modular structure also enables equilibrium and biased transport simulations, from Kubo-type linear response to open-boundary approaches under external potentials, highlighting CP2K's unique capability to unify quantum chemistry with quantum and statistical mechanics within a versatile, holistic simulation environment. △ Less

Submitted 30 July, 2026; v1 submitted 24 July, 2026; originally announced July 2026.

On the Regularity and Interpolation of Coupled Cluster Amplitudes in Canonical Orbital Basis

Authors: Jonas Beck, Benjamin Stamm

Abstract: Arguably the most widely used approaches for obtaining highly accurate molecular ground-state energies are coupled cluster methods. Despite introducing two layers of approximation, a linear and a nonlinear one, coupled cluster methods remain computationally intensive, with the complexity scaling as $O(poly(N))$, where $N$ is the number of electrons. Moreover, this method must be applied over a lar… ▽ More Arguably the most widely used approaches for obtaining highly accurate molecular ground-state energies are coupled cluster methods. Despite introducing two layers of approximation, a linear and a nonlinear one, coupled cluster methods remain computationally intensive, with the complexity scaling as $O(poly(N))$, where $N$ is the number of electrons. Moreover, this method must be applied over a large set of different nuclear coordinates in order to study certain chemical phenomena. Therefore, in this work, we investigate the regularity of single-reference coupled cluster amplitudes with respect to nuclear coordinate displacements, with the aim of enabling interpolation or extrapolation approaches that rely on only a limited number of reference geometries. We show that, in theory, under certain non-degeneracy assumptions on the Hartree-Fock level of theory, and the coupled cluster level of theory the amplitudes behave real analytic. Furthermore, we analyze the artifacts that arise in practical calculations that use canonical orbitals, which hinder this high degree of regularity, and suggest strategies to mitigate these issues. Finally, we validate our findings through numerical experiments by interpolating the amplitudes and comparing the performance of the interpolants with that of the exact amplitudes. △ Less

Submitted 1 July, 2026; v1 submitted 21 May, 2026; originally announced May 2026.

Comments: 29 pages, 4 figures

MSC Class: Primary 81V55; 41A05; Secondary 65D05; 81-08

Model Order Reduction for Parametric Hermitian Eigenvalue Problems: Local Acceleration with Taylor-Reduced Basis Method

Authors: Benjamin Stamm, Zhuoyao Zeng

Abstract: This paper is concerned with the Taylor-reduced basis method (Taylor-RBM) for the efficient approximation of eigenspaces of large scale parametric Hermitian matrices. The Taylor-RBM is a local model order reduction method, which constructs an approximation space by capturing derivatives information of the spectral projector at a reference point in the parameter domain. We perform a concise error a… ▽ More This paper is concerned with the Taylor-reduced basis method (Taylor-RBM) for the efficient approximation of eigenspaces of large scale parametric Hermitian matrices. The Taylor-RBM is a local model order reduction method, which constructs an approximation space by capturing derivatives information of the spectral projector at a reference point in the parameter domain. We perform a concise error analysis to justify the Taylor-RBM for eigenvalue problems, and we present a computationally efficient procedure to assemble the Taylor-reduced basis space. Since this method is tightly connected to the classical multivariate analytic perturbation theory, we also provide a detailed analysis of the spectral approximation using the truncated power series of the eigenprojector, and compare this with the approximation obtained from the Taylor-RBM. △ Less

Submitted 27 March, 2026; v1 submitted 15 March, 2026; originally announced March 2026.

Comments: improved formulations, fixed typos, refined numerical experiments

MSC Class: Primary 41A63; 65F15; secondary 47A55; 65N25; 65F50; 65K10; 81P68

Simulation of the magnetic Ginzburg-Landau equation via vortex tracking

Authors: Thiago Carvalho Corso, Gaspard Kemlin, Christof Melcher, Benjamin Stamm

Abstract: This paper deals with the numerical simulation of the 2D magnetic time-dependent Ginzburg-Landau (TDGL) equations in the regime of small but finite (inverse) Ginzburg-Landau parameter $ε$ and constant (order $1$ in $ε$) applied magnetic field. In this regime, a well-known feature of the TDGL equation is the appearance of quantized vortices with core size of order $ε$. Moreover, in the singular lim… ▽ More This paper deals with the numerical simulation of the 2D magnetic time-dependent Ginzburg-Landau (TDGL) equations in the regime of small but finite (inverse) Ginzburg-Landau parameter $ε$ and constant (order $1$ in $ε$) applied magnetic field. In this regime, a well-known feature of the TDGL equation is the appearance of quantized vortices with core size of order $ε$. Moreover, in the singular limit $ε\searrow 0$, these vortices evolve according to an explicit ODE system. In this work, we first introduce a new numerical method for the numerical integration of this limiting ODE system, which requires to solve a linear second order PDE at each time step. We also provide a rigorous theoretical justification for this method that applies to a general class of 2D domains. We then develop and analyze a numerical strategy based on the finite-dimensional ODE system to efficiently simulate the infinite-dimensional TDGL equations in the presence of a constant external magnetic field and for small, but finite, $ε$. This method allows us to avoid resolving the $ε$-scale when solving the TDGL equations, where small values of $ε$ typically require very fine meshes and time steps. We provide numerical examples on a few test cases and justify the accuracy of the method with numerical investigations. We end the paper showing that, in the mixed flow case, the limiting ODE system is able to capture the crystallization process in which, for large times, the vortices arrange into a stable pattern. △ Less

Submitted 16 December, 2025; v1 submitted 30 October, 2025; originally announced October 2025.

Grassmanian Interpolation of Low-Pass Graph Filters: Theory and Applications

Authors: Anton Savostianov, Michael T. Schaub, Benjamin Stamm

Abstract: Low-pass graph filters are fundamental for signal processing on graphs and other non-Euclidean domains. However, the computation of such filters for parametric graph families can be prohibitively expensive as computation of the corresponding low-frequency subspaces, requires the repeated solution of an eigenvalue problem. We suggest a novel algorithm of low-pass graph filter interpolation based on… ▽ More Low-pass graph filters are fundamental for signal processing on graphs and other non-Euclidean domains. However, the computation of such filters for parametric graph families can be prohibitively expensive as computation of the corresponding low-frequency subspaces, requires the repeated solution of an eigenvalue problem. We suggest a novel algorithm of low-pass graph filter interpolation based on Riemannian interpolation in normal coordinates on the Grassmann manifold. We derive an error bound estimate for the subspace interpolation and suggest two possible applications for induced parametric graph families. First, we argue that the temporal evolution of the node features may be translated to the evolving graph topology via a similarity correction to adjust the homophily degree of the network. Second, we suggest a dot product graph family induced by a given static graph which allows to infer improved message passing scheme for node classification facilitated by the filter interpolation. △ Less

Submitted 27 October, 2025; originally announced October 2025.

Comments: 13 pages

Parallel athermal quasistatic deformation stepping of molecular systems

Authors: Maximilian Reihn, Franz Bamer, Benjamin Stamm

Abstract: The athermal quasistatic deformation method provides an elegant solution to overcome the limitation of short time spans in molecular simulations. It provides overdamped conditions, allowing for the extraction of purely structural responses in the absence of thermal vibration. However, it requires computationally expensive sequences of affine deformation followed by minimization of the potential en… ▽ More The athermal quasistatic deformation method provides an elegant solution to overcome the limitation of short time spans in molecular simulations. It provides overdamped conditions, allowing for the extraction of purely structural responses in the absence of thermal vibration. However, it requires computationally expensive sequences of affine deformation followed by minimization of the potential energy to incrementally find the path in the potential energy landscape that corresponds to the correct solution trajectory. Therefore, we propose an athermal parallel stepping scheme that significantly improves the computational time necessary to find the correct solution trajectory using a multi-thread approach. Our approach proposes stepping at two levels. Level I stepping provides a sequence of initial guesses at large increments by affine deformation of the system and land-marking anchor points on the potential energy landscape. Level II stepping performs a set of individual finely resolved athermal quasistatic deformation steps between the inherent structures of the initial level I guesses executed in parallel. The evaluated candidate trajectory is then verified by consecutively comparing the configuration of every last level II result with the corresponding inherent structure of the level I guesses at the same strain states. If the two configurations are not equivalent, the solution must be rejected and recalculated from this point. Rigorous numerical testing with $4,8,16$ and $32$ parallel threads and different values of hyper-parameters demonstrates that our method achieves computational average speed-ups of factors ranging from $2.02$ to $6.33$, while maintaining simulation accuracy, offering a powerful new tool for athermal molecular simulations. △ Less

Submitted 29 April, 2026; v1 submitted 28 July, 2025; originally announced July 2025.

Certified Model Order Reduction for parametric Hermitian eigenproblems

Authors: Mattia Manucci, Benjamin Stamm, Zhuoyao Zeng

Abstract: This article deals with the efficient and certified numerical approximation of the smallest eigenvalue and the associated eigenspace of a large-scale parametric Hermitian matrix. For this aim, we rely on projection-based model order reduction (MOR), i.e., we approximate the large-scale problem by projecting it onto a suitable subspace and reducing it to one of a much smaller dimension. Such a subs… ▽ More This article deals with the efficient and certified numerical approximation of the smallest eigenvalue and the associated eigenspace of a large-scale parametric Hermitian matrix. For this aim, we rely on projection-based model order reduction (MOR), i.e., we approximate the large-scale problem by projecting it onto a suitable subspace and reducing it to one of a much smaller dimension. Such a subspace is constructed by means of weak greedy-type strategies. After detailing the connections with the reduced basis method for source problems, we introduce a novel error estimate for the approximation error related to the eigenspace associated with the smallest eigenvalue. Since the difference between the second smallest and the smallest eigenvalue, the so-called spectral gap, is crucial for the reliability of the error estimate, we propose efficiently computable upper and lower bounds for higher eigenvalues and for the spectral gap, which enable the assembly of a subspace for the MOR approximation of the spectral gap. Based on that, a second subspace is then generated for the MOR approximation of the eigenspace associated with the smallest eigenvalue. We also provide efficiently computable conditions to ensure that the multiplicity of the smallest eigenvalue is fully captured in the reduced space. This work is motivated by a specific application: the repeated identifications of the states with minimal energy, the so-called ground states, of parametric quantum spin system models. △ Less

Submitted 13 January, 2026; v1 submitted 3 April, 2025; originally announced April 2025.

Comments: various parts substantially reviewed and rewritten. miscellaneous typos fixed. latex compilation error fixed

MSC Class: 41A63; 65F15; 65N25; 65F50; 65K10; 81P68

A symmetry-preserving and transferable representation for learning the Kohn-Sham density matrix

Authors: Liwei Zhang, Patrizia Mazzeo, Michele Nottoli, Edoardo Cignoni, Lorenzo Cupellini, Benjamin Stamm

Abstract: The Kohn-Sham (KS) density matrix is one of the most essential properties in KS density functional theory (DFT), from which many other physical properties of interest can be derived. In this work, we present a parameterized representation for learning the mapping from a molecular configuration to its corresponding density matrix using the Atomic Cluster Expansion (ACE) framework, which preserves t… ▽ More The Kohn-Sham (KS) density matrix is one of the most essential properties in KS density functional theory (DFT), from which many other physical properties of interest can be derived. In this work, we present a parameterized representation for learning the mapping from a molecular configuration to its corresponding density matrix using the Atomic Cluster Expansion (ACE) framework, which preserves the physical symmetries of the mapping, including isometric equivariance and Grassmannianity. Trained on several typical molecules, the proposed representation is shown to be systematically improvable with the increase of the model parameters and is transferable to molecules that are not part of and even more complex than those in the training set. The models generated by the proposed approach are illustrated as being able to generate reasonable predictions of the density matrix to either accelerate the DFT calculations or to provide approximations to some properties of the molecules. △ Less

Submitted 11 March, 2025; originally announced March 2025.

Comments: Including Supporting Information

Gradient Flow Finite Element Discretisations with Energy-Based $hp$-Adaptivity for the Gross-Pitaevskii Equation with Angular Momentum

Authors: Pascal Heid, Paul Houston, Benjamin Stamm, Thomas P. Wihler

Abstract: This article deals with the stationary Gross-Pitaevskii non-linear eigenvalue problem in the presence of a rotating magnetic field that is used to model macroscopic quantum effects such as Bose-Einstein condensates (BECs). In this regime, the ground-state wave-function can exhibit an a priori unknown number of quantum vortices at unknown locations, which necessitates the exploitation of adaptive n… ▽ More This article deals with the stationary Gross-Pitaevskii non-linear eigenvalue problem in the presence of a rotating magnetic field that is used to model macroscopic quantum effects such as Bose-Einstein condensates (BECs). In this regime, the ground-state wave-function can exhibit an a priori unknown number of quantum vortices at unknown locations, which necessitates the exploitation of adaptive numerical strategies. To this end, we consider the conforming finite element method in combination with a discrete Sobolev gradient descent, which is guided by the energy-topology of the problem, to address the nonlinearity. In addition, a key novelty of this work is an $hp$-adaptive strategy that is solely based on energy decay rather than a posteriori error estimators for the refinement process. Numerical results demonstrate that the $hp$-adaptive strategy is highly efficient in terms of accuracy to compute the ground-state wave function and energy for several test problems where we observe exponential convergence. △ Less

Submitted 18 December, 2025; v1 submitted 23 December, 2024; originally announced December 2024.

Approximations of the Iterative Stockholder Analysis scheme using exponential basis functions

Authors: YingXing Cheng, Benjamin Stamm

Abstract: In this work, we introduce several approximations of the Iterative Stockholder Analysis (ISA) method based on exponential basis functions. These approximations are categorized into linear and non-linear models, referred to as LISA and NLIS, respectively. By particular choices of hyperparameters in the NLIS model, both LISA and the Minimal-Basis Iterative Stockholder (MBIS) method can be reproduced… ▽ More In this work, we introduce several approximations of the Iterative Stockholder Analysis (ISA) method based on exponential basis functions. These approximations are categorized into linear and non-linear models, referred to as LISA and NLIS, respectively. By particular choices of hyperparameters in the NLIS model, both LISA and the Minimal-Basis Iterative Stockholder (MBIS) method can be reproduced. Four LISA variants are constructed using systematically generated exponential basis functions derived from the NLIS model applied to atomic systems. The performance of these LISA variants and NLIS models is benchmarked on 15 small molecules, including neutral, anionic, and cationic species. To facilitate comparison, we propose several metrics designed to highlight differences between the methods. Our results demonstrate that LISA, employing Gaussian basis functions derived from the NLIS model on isolated atomic systems, achieves an optimal balance of computational accuracy, robustness, and efficiency, particularly in minimizing the objective function. △ Less

Submitted 30 June, 2025; v1 submitted 6 December, 2024; originally announced December 2024.

Fully guaranteed and computable error bounds on the energy for periodic Kohn-Sham equations with convex density functionals

Authors: Andrea Bordignon, Geneviève Dusson, Éric Cancès, Gaspard Kemlin, Rafael Antonio Lainez Reyes, Benjamin Stamm

Abstract: In this article, we derive fully guaranteed error bounds for the energy of convex nonlinear mean-field models. These results apply in particular to Kohn-Sham equations with convex density functionals, which includes the reduced Hartree-Fock (rHF) model, as well as the Kohn-Sham model with exact exchange-density functional (which is unfortunately not explicit and therefore not usable in practice).… ▽ More In this article, we derive fully guaranteed error bounds for the energy of convex nonlinear mean-field models. These results apply in particular to Kohn-Sham equations with convex density functionals, which includes the reduced Hartree-Fock (rHF) model, as well as the Kohn-Sham model with exact exchange-density functional (which is unfortunately not explicit and therefore not usable in practice). We then decompose the obtained bounds into two parts, one depending on the chosen discretization and one depending on the number of iterations performed in the self-consistent algorithm used to solve the nonlinear eigenvalue problem, paving the way for adaptive refinement strategies. The accuracy of the bounds is demonstrated on a series of test cases, including a Silicon crystal and an Hydrogen Fluoride molecule simulated with the rHF model and discretized with planewaves. We also show that, although not anymore guaranteed, the error bounds remain very accurate for a Silicon crystal simulated with the Kohn-Sham model using nonconvex exchangecorrelation functionals of practical interest. △ Less

Submitted 19 February, 2025; v1 submitted 18 September, 2024; originally announced September 2024.

On reduced basis methods for eigenvalue problems, and on its coupling with perturbation theory

Authors: Louis Garrigue, Benjamin Stamm

Abstract: In this article, we study eigenvalue problems associated to self-adjoint operators and their approximation obtained by subspace projection, as used in the reduced basis method for instance. We provide error bounds between the exact eigenmodes and the approximated ones and also consider degenerate cases in the analysis. When the operator depends on a parameter, we apply the bounds assuming that the… ▽ More In this article, we study eigenvalue problems associated to self-adjoint operators and their approximation obtained by subspace projection, as used in the reduced basis method for instance. We provide error bounds between the exact eigenmodes and the approximated ones and also consider degenerate cases in the analysis. When the operator depends on a parameter, we apply the bounds assuming that the reduced space contains the derivatives of the eigenfunction with respect to the parameter. Finally, we provide some numerical examples that reflect the analytical results. △ Less

Submitted 13 June, 2025; v1 submitted 21 August, 2024; originally announced August 2024.

Multi-center decomposition of molecular densities: A numerical perspective

Authors: YingXing Cheng, Eric Cancès, Virginie Ehrlacher, Alston J. Misquitta, Benjamin Stamm

Abstract: In this study, we analyze various Iterative Stockholder Analysis (ISA) methods for molecular density partitioning, focusing on the numerical performance of the recently proposed Linear approximation of Iterative Stockholder Analysis model (LISA) [J. Chem. Phys. 156, 164107 (2022)]. We first provide a systematic derivation of various iterative solvers to find the unique LISA solution. In a subseque… ▽ More In this study, we analyze various Iterative Stockholder Analysis (ISA) methods for molecular density partitioning, focusing on the numerical performance of the recently proposed Linear approximation of Iterative Stockholder Analysis model (LISA) [J. Chem. Phys. 156, 164107 (2022)]. We first provide a systematic derivation of various iterative solvers to find the unique LISA solution. In a subsequent systematic numerical study, we evaluate their performance on 48 organic and inorganic, neutral and charged molecules and also compare LISA to two other well-known ISA variants: the Gaussian Iterative Stockholder Analysis (GISA) and Minimum Basis Iterative Stockholder analysis (MBIS). The study reveals that LISA-family methods can offer a numerically more efficient approach with better accuracy compared to the two comparative methods. Moreover, the well-known issue with the MBIS method, where atomic charges obtained for negatively charged molecules are anomalously negative, is not observed in LISA-family methods. Despite the fact that LISA occasionally exhibits elevated entropy as a consequence of the absence of more diffuse basis functions, this issue can be readily mitigated by incorporating additional or integrating supplementary basis functions within the LISA framework. This research provides the foundation for future studies on the efficiency and chemical accuracy of molecular density partitioning schemes. △ Less

Submitted 24 October, 2024; v1 submitted 14 May, 2024; originally announced May 2024.

Numerical simulation of the Gross-Pitaevskii equation via vortex tracking

Authors: Thiago Carvalho Corso, Gaspard Kemlin, Christof Melcher, Benjamin Stamm

Abstract: This paper deals with the numerical simulation of the Gross-Pitaevskii (GP) equation, for which a well-known feature is the appearance of quantized vortices with core size of the order of a small parameter $\varepsilon$. Without a magnetic field and with suitable initial conditions, these vortices interact, in the singular limit $\varepsilon\to0$, through an explicit Hamiltonian dynamics. Using th… ▽ More This paper deals with the numerical simulation of the Gross-Pitaevskii (GP) equation, for which a well-known feature is the appearance of quantized vortices with core size of the order of a small parameter $\varepsilon$. Without a magnetic field and with suitable initial conditions, these vortices interact, in the singular limit $\varepsilon\to0$, through an explicit Hamiltonian dynamics. Using this analytical framework, we develop and analyze a numerical strategy based on the reduced-order Hamiltonian system to efficiently simulate the infinite-dimensional GP equation for small, but finite, $\varepsilon$. This method allows us to avoid numerical stability issues in solving the GP equation, where small values of $\varepsilon$ typically require very fine meshes and time steps. We also provide a mathematical justification of our method in terms of rigorous error estimates of the error in the supercurrent, together with numerical illustrations. △ Less

Submitted 25 October, 2024; v1 submitted 2 April, 2024; originally announced April 2024.

Trace estimates for harmonic functions along circular arcs with applications to domain decomposition on overlapping disks

Authors: Thiago Carvalho Corso, Muhammad Hassan, Abhinav Jha, Benjamin Stamm

Abstract: In this paper we derive several (and in many cases sharp) estimates for the $\mathrm{L}^2$-trace norm of harmonic functions along circular arcs. More precisely, we obtain geometry-dependent estimates on the norm, spectral radius, and numerical range of the Dirchlet-to-Dirichlet (DtD) operator sending data on the boundary of the disk to the restriction of its harmonic extension along circular arcs… ▽ More In this paper we derive several (and in many cases sharp) estimates for the $\mathrm{L}^2$-trace norm of harmonic functions along circular arcs. More precisely, we obtain geometry-dependent estimates on the norm, spectral radius, and numerical range of the Dirchlet-to-Dirichlet (DtD) operator sending data on the boundary of the disk to the restriction of its harmonic extension along circular arcs inside the disk. The estimates we derive here have applications in the convergence analysis of the Schwarz domain decomposition method for overlapping disks in two dimensions. In particular, they allow us to establish a rigorous convergence proof for the discrete parallel Schwarz method applied to the Conductor-like Screening Model (COSMO) from theoretical chemistry in the two-disk case, and to derive error estimates with respect to the discretization parameter, the number of Schwarz iterations, and the geometry of the domain. Our analysis addresses challenges beyond classical domain decomposition theory, especially the weak enforcement of boundary conditions. △ Less

Submitted 13 November, 2024; v1 submitted 29 January, 2024; originally announced January 2024.

MSC Class: 35J05; 30C40; 65N55

A Scalable Two-Level Domain Decomposition Eigensolver for Periodic Schrödinger Eigenstates in Anisotropically Expanding Domains

Authors: Lambert Theisen, Benjamin Stamm

Abstract: Accelerating iterative eigenvalue algorithms is often achieved by employing a spectral shifting strategy. Unfortunately, improved shifting typically leads to a smaller eigenvalue for the resulting shifted operator, which in turn results in a high condition number of the underlying solution matrix, posing a major challenge for iterative linear solvers. This paper introduces a two-level domain decom… ▽ More Accelerating iterative eigenvalue algorithms is often achieved by employing a spectral shifting strategy. Unfortunately, improved shifting typically leads to a smaller eigenvalue for the resulting shifted operator, which in turn results in a high condition number of the underlying solution matrix, posing a major challenge for iterative linear solvers. This paper introduces a two-level domain decomposition preconditioner that addresses this issue for the linear Schrödinger eigenvalue problem, even in the presence of a vanishing eigenvalue gap in non-uniform, expanding domains. Since the quasi-optimal shift, which is already available as the solution to a spectral cell problem, is required for the eigenvalue solver, it is logical to also use its associated eigenfunction as a generator to construct a coarse space. We analyze the resulting two-level additive Schwarz preconditioner and obtain a condition number bound that is independent of the domain's anisotropy, despite the need for only one basis function per subdomain for the coarse solver. Several numerical examples are presented to illustrate its flexibility and efficiency. △ Less

Submitted 18 June, 2024; v1 submitted 15 November, 2023; originally announced November 2023.

Comments: 26 pages, 7 figures, 2 tables

MSC Class: 65N25; 65F15; 65N30; 65F10; 65N22; 65N55; 65F08

Journal ref: SIAM Journal on Scientific Computing Vol. 46, Iss. 5 (2024)

Domain Decomposition Method for Poisson--Boltzmann Equations based on Solvent Excluded Surface

Authors: Abhinav Jha, Benjamin Stamm

Abstract: In this paper, we develop a domain decomposition method for the nonlinear Poisson-Boltzmann equation based on a solvent-excluded surface widely used in computational chemistry. The model relies on a nonlinear equation defined in $\mathbb{R}^3$ with a space-dependent dielectric permittivity and an ion-exclusion function that accounts for steric effects. Potential theory arguments transform the nonl… ▽ More In this paper, we develop a domain decomposition method for the nonlinear Poisson-Boltzmann equation based on a solvent-excluded surface widely used in computational chemistry. The model relies on a nonlinear equation defined in $\mathbb{R}^3$ with a space-dependent dielectric permittivity and an ion-exclusion function that accounts for steric effects. Potential theory arguments transform the nonlinear equation into two coupled equations defined in a bounded domain. Then, the Schwarz decomposition method is used to formulate local problems by decomposing the cavity into overlapping balls and only solving a set of coupled sub-equations in each ball. The main novelty of the proposed method is the introduction of a hybrid linear-nonlinear solver used to solve the equation. A series of numerical experiments are presented to test the method and show the importance of the nonlinear model. △ Less

Submitted 5 April, 2024; v1 submitted 13 September, 2023; originally announced September 2023.

MSC Class: 65N35; 65N55

A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics

Authors: Federica Pes, Ètienne Polack, Patrizia Mazzeo, Geneviève Dusson, Benjamin Stamm, Filippo Lipparini

Abstract: This article proposes a so-called Quasi Time-Reversible (QTR G-Ext) scheme based on Grassmann extrapolation of density matrices for an accurate calculation of initial guesses in Born-Oppenheimer Molecular Dynamics simulations. The method shows excellent results on four large molecular systems, ranging from 21 to 94 atoms simulated with Kohn-Sham density functional theory surrounded with a classica… ▽ More This article proposes a so-called Quasi Time-Reversible (QTR G-Ext) scheme based on Grassmann extrapolation of density matrices for an accurate calculation of initial guesses in Born-Oppenheimer Molecular Dynamics simulations. The method shows excellent results on four large molecular systems, ranging from 21 to 94 atoms simulated with Kohn-Sham density functional theory surrounded with a classical environment with 6k to 16k atoms. Namely, it clearly reduces the number of self-consistent field iterations, while keeping a similar energy drift as in the extended Lagrangian Born-Oppenheimer method. △ Less

Submitted 11 July, 2023; originally announced July 2023.

Embedded corrector problems for homogenization in linear elasticity

Authors: Virginie Ehrlacher, Frederic Legoll, Benjamin Stamm, Shuyang Xiang

Abstract: In this article, we extend the study of embedded corrector problems, that we have previously introduced in the context of the homogenization of scalar diffusive equations, to the context of homogenized elastic properties of materials. This extension is not trivial and requires mathematical arguments specific to the elasticity case. Starting from a linear elasticity model with highly-oscillatory co… ▽ More In this article, we extend the study of embedded corrector problems, that we have previously introduced in the context of the homogenization of scalar diffusive equations, to the context of homogenized elastic properties of materials. This extension is not trivial and requires mathematical arguments specific to the elasticity case. Starting from a linear elasticity model with highly-oscillatory coefficients, we introduce several effective approximations of the homogenized tensor. These approximations are based on the solution to an embedded corrector problem, where a finite-size domain made of the linear elastic heterogeneous material is embedded in a linear elastic homogeneous infinite medium, the constant elasticity tensor of which has to be appropriately determined. The approximations we provide are proven to converge to the homogenized elasticity tensor when the size of the embedded domain tends to infinity. Some particular attention is devoted to the case of isotropic materials. △ Less

Submitted 7 July, 2023; originally announced July 2023.

A multipoint perturbation formula for eigenvalue problems

Authors: Geneviève Dusson, Louis Garrigue, Benjamin Stamm

Abstract: Standard perturbation theory of eigenvalue problems consists of obtaining approximations of eigenmodes in the neighborhood of an operator where the corresponding eigenmode is known. Nevertheless, if the corresponding eigenmodes of several nearby operators are known, standard perturbation theory cannot simultaneously use all this knowledge to provide a better approximation. We derive a resolvent fo… ▽ More Standard perturbation theory of eigenvalue problems consists of obtaining approximations of eigenmodes in the neighborhood of an operator where the corresponding eigenmode is known. Nevertheless, if the corresponding eigenmodes of several nearby operators are known, standard perturbation theory cannot simultaneously use all this knowledge to provide a better approximation. We derive a resolvent formula enabling such an approximation result, and provide numerical examples for which this method is more competitive than standard perturbation theory. △ Less

Submitted 28 July, 2025; v1 submitted 14 May, 2023; originally announced May 2023.

Journal ref: ESAIM: Mathematical Modelling and Numerical Analysis 59 (4), 2081-2109

Reduced basis surrogates for quantum spin systems based on tensor networks

Authors: Paul Brehmer, Michael F. Herbst, Stefan Wessel, Matteo Rizzi, Benjamin Stamm

Abstract: Within the reduced basis methods approach, an effective low-dimensional subspace of a quantum many-body Hilbert space is constructed in order to investigate, e.g., the ground-state phase diagram. The basis of this subspace is built from solutions of snapshots, i.e., ground states corresponding to particular and well-chosen parameter values. Here, we show how a greedy strategy to assemble the reduc… ▽ More Within the reduced basis methods approach, an effective low-dimensional subspace of a quantum many-body Hilbert space is constructed in order to investigate, e.g., the ground-state phase diagram. The basis of this subspace is built from solutions of snapshots, i.e., ground states corresponding to particular and well-chosen parameter values. Here, we show how a greedy strategy to assemble the reduced basis and thus to select the parameter points can be implemented based on matrix-product-states (MPS) calculations. Once the reduced basis has been obtained, observables required for the computation of phase diagrams can be computed with a computational complexity independent of the underlying Hilbert space for any parameter value. We illustrate the efficiency and accuracy of this approach for different one-dimensional quantum spin-1 models, including anisotropic as well as biquadratic exchange interactions, leading to rich quantum phase diagrams. △ Less

Submitted 12 July, 2023; v1 submitted 26 April, 2023; originally announced April 2023.

Comments: 15 pages, 13 figures; extended conclusions, typos corrected

Numerical stability and efficiency of response property calculations in density functional theory

Authors: Eric Cancès, Michael F. Herbst, Gaspard Kemlin, Antoine Levitt, Benjamin Stamm

Abstract: Response calculations in density functional theory aim at computing the change in ground-state density induced by an external perturbation. At finite temperature these are usually performed by computing variations of orbitals, which involve the iterative solution of potentially badly-conditioned linear systems, the Sternheimer equations. Since many sets of variations of orbitals yield the same var… ▽ More Response calculations in density functional theory aim at computing the change in ground-state density induced by an external perturbation. At finite temperature these are usually performed by computing variations of orbitals, which involve the iterative solution of potentially badly-conditioned linear systems, the Sternheimer equations. Since many sets of variations of orbitals yield the same variation of density matrix this involves a choice of gauge. Taking a numerical analysis point of view we present the various gauge choices proposed in the literature in a common framework and study their stability. Beyond existing methods we propose a new approach, based on a Schur complement using extra orbitals from the self-consistent-field calculations, to improve the stability and efficiency of the iterative solution of Sternheimer equations. We show the success of this strategy on nontrivial examples of practical interest, such as Heusler transition metal alloy compounds, where savings of around 40% in the number of required cost-determining Hamiltonian applications have been achieved. △ Less

Submitted 16 February, 2023; v1 submitted 10 October, 2022; originally announced October 2022.

Journal ref: Letters in Mathematical Physics, 2023, 113 (1)

Linear scaling computation of forces for the domain-decomposition linear Poisson--Boltzmann method

Authors: Abhinav Jha, Michele Nottoli, Aleksandr Mikhalev, Chaoyu Quan, Benjamin Stamm

Abstract: The Linearized Poisson--Boltzmann (LPB) equation is a popular and widely accepted model for accounting solvent effects in computational (bio-) chemistry. In the present article we derive the analytical forces of the domain-decomposition-based ddLPB-method with vdW or SAS surface. We present an efficient strategy to compute the forces and its implementation, allowing linear scaling of the method wi… ▽ More The Linearized Poisson--Boltzmann (LPB) equation is a popular and widely accepted model for accounting solvent effects in computational (bio-) chemistry. In the present article we derive the analytical forces of the domain-decomposition-based ddLPB-method with vdW or SAS surface. We present an efficient strategy to compute the forces and its implementation, allowing linear scaling of the method with respect to the number of atoms using the fast multipole method (FMM). Numerical tests illustrates the accuracy of the computation of the analytical forces and compares efficiency with other available methods. △ Less

Submitted 9 December, 2022; v1 submitted 1 March, 2022; originally announced March 2022.

Surrogate models for quantum spin systems based on reduced order modeling

Authors: Michael F. Herbst, Stefan Wessel, Matteo Rizzi, Benjamin Stamm

Abstract: We present a methodology to investigate phase-diagrams of quantum models based on the principle of the reduced basis method (RBM). The RBM is built from a few ground-state snapshots, i.e., lowest eigenvectors of the full system Hamiltonian computed at well-chosen points in the parameter space of interest. We put forward a greedy-strategy to assemble such small-dimensional basis, i.e., to select wh… ▽ More We present a methodology to investigate phase-diagrams of quantum models based on the principle of the reduced basis method (RBM). The RBM is built from a few ground-state snapshots, i.e., lowest eigenvectors of the full system Hamiltonian computed at well-chosen points in the parameter space of interest. We put forward a greedy-strategy to assemble such small-dimensional basis, i.e., to select where to spend the numerical effort needed for the snapshots. Once the RBM is assembled, physical observables required for mapping out the phase-diagram (e.g., structure factors) can be computed for any parameter value with a modest computational complexity, considerably lower than the one associated to the underlying Hilbert space dimension. We benchmark the method in two test cases, a chain of excited Rydberg atoms and a geometrically frustrated antiferromagnetic two-dimensional lattice model, and illustrate the accuracy of the approach. In particular, we find that the ground-manifold can be approximated to sufficient accuracy with a moderate number of basis functions, which increases very mildly when the number of microscopic constituents grows - in stark contrast to the exponential growth of the Hilbert space needed to describe each of the few snapshots. A combination of the presented RBM approach with other numerical techniques circumventing even the latter big cost, e.g., Tensor Network methods, is a tantalising outlook of this work. △ Less

Submitted 11 March, 2022; v1 submitted 29 October, 2021; originally announced October 2021.

A Quasi-Optimal Factorization Preconditioner for Periodic Schrödinger Eigenstates in Anisotropically Expanding Domains

Authors: Benjamin Stamm, Lambert Theisen

Abstract: This paper provides a provably quasi-optimal preconditioning strategy of the linear Schrödinger eigenvalue problem with periodic potentials for a possibly non-uniform spatial expansion of the domain. The quasi-optimality is achieved by having the iterative eigenvalue algorithms converge in a constant number of iterations for different domain sizes. In the analysis, we derive an analytic factorizat… ▽ More This paper provides a provably quasi-optimal preconditioning strategy of the linear Schrödinger eigenvalue problem with periodic potentials for a possibly non-uniform spatial expansion of the domain. The quasi-optimality is achieved by having the iterative eigenvalue algorithms converge in a constant number of iterations for different domain sizes. In the analysis, we derive an analytic factorization of the spectrum and asymptotically describe it using concepts from the homogenization theory. This decomposition allows us to express the eigenpair as an easy-to-calculate cell problem solution combined with an asymptotically vanishing remainder. We then prove that the easy-to-calculate limit eigenvalue can be used in a shift-and-invert preconditioning strategy to bound the number of eigensolver iterations uniformly. Several numerical examples illustrate the effectiveness of this quasi-optimal preconditioning strategy. △ Less

Submitted 24 May, 2022; v1 submitted 28 October, 2021; originally announced October 2021.

Comments: 29 pages, 9 figures, 2 tables

MSC Class: 65N25; 65F15; 65N30; 35B27; 35B40

Journal ref: SIAM Journal on Numerical Analysis, Vol. 60, Iss. 5 (2022)

Multi-center decomposition of molecular densities: a mathematical perspective

Authors: Robert Benda, Eric Cancès, Virginie Ehrlacher, Benjamin Stamm

Abstract: The aim of this paper is to analyze from a mathematical perspective some existing schemes to partition a molecular density into several atomic contributions, with a specific focus on Iterative Stockholder Atom (ISA) methods. We provide a unified mathematical framework to describe the latter family of methods and propose a new scheme, named L-ISA (for linear approximation of ISA). We prove several… ▽ More The aim of this paper is to analyze from a mathematical perspective some existing schemes to partition a molecular density into several atomic contributions, with a specific focus on Iterative Stockholder Atom (ISA) methods. We provide a unified mathematical framework to describe the latter family of methods and propose a new scheme, named L-ISA (for linear approximation of ISA). We prove several important mathematical properties of the ISA and L-ISA minimization problems and show that the so-called ISA algorithms can be viewed as alternating minimization schemes, which in turn enables us to obtain new convergence results for these numerical methods. Specific mathematical properties of the ISA decomposition for diatomic systems are also presented. We also review the basis-space oriented Distributed Multipole Analysis method, the mathematical formulation of which is also clarified. Different schemes are numerically compared on different molecules and we discuss the advantages and drawbacks of each approach. △ Less

Submitted 7 September, 2021; originally announced September 2021.

MSC Class: 81Q99

Grassmann extrapolation of density matrices for Born-Oppenheimer molecular dynamics

Authors: Etienne Polack, Geneviève Dusson, Benjamin Stamm, Filippo Lipparini

Abstract: Born-Oppenheimer Molecular Dynamics (BOMD) is a powerful but expensive technique. The main bottleneck in a density functional theory bomd calculation is the solution to the Kohn-Sham (KS) equations, that requires an iterative procedure that starts from a guess for the density matrix. Converged densities from previous points in the trajectory can be used to extrapolate a new guess, however, the non… ▽ More Born-Oppenheimer Molecular Dynamics (BOMD) is a powerful but expensive technique. The main bottleneck in a density functional theory bomd calculation is the solution to the Kohn-Sham (KS) equations, that requires an iterative procedure that starts from a guess for the density matrix. Converged densities from previous points in the trajectory can be used to extrapolate a new guess, however, the non-linear constraint that an idempotent density needs to satisfy make the direct use of standard linear extrapolation techniques not possible. In this contribution, we introduce a locally bijective map between the manifold where the density is defined and its tangent space, so that linear extrapolation can be performed in a vector space while, at the same time, retaining the correct physical properties of the extrapolated density using molecular descriptors. We apply the method to real-life, multiscale polarizable QM/MM. △ Less

Submitted 22 September, 2021; v1 submitted 28 July, 2021; originally announced July 2021.

Continuity estimates for Riesz potentials on polygonal boundaries

Authors: Xavier Claeys, Muhammad Hassan, Benjamin Stamm

Abstract: Riesz potentials are well known objects of study in the theory of singular integrals that have been the subject of recent, increased interest from the numerical analysis community due to their connections with fractional Laplace problems and proposed use in certain domain decomposition methods. While the L$^p$-mapping properties of Riesz potentials on flat geometries are well-established, comparab… ▽ More Riesz potentials are well known objects of study in the theory of singular integrals that have been the subject of recent, increased interest from the numerical analysis community due to their connections with fractional Laplace problems and proposed use in certain domain decomposition methods. While the L$^p$-mapping properties of Riesz potentials on flat geometries are well-established, comparable results on rougher geometries for Sobolev spaces are very scarce. In this article, we study the continuity properties of the surface Riesz potential generated by the $1/\sqrt{x}$ singular kernel on a polygonal domain $Ω\subset \mathbb{R}^2$. We prove that this surface Riesz potential maps L$^{2}(\partialΩ)$ into H$^{+1/2}(\partialΩ)$. Our proof is based on a careful analysis of the Riesz potential in the neighbourhood of corners of the domain $Ω$. The main tool we use for this corner analysis is the Mellin transform which can be seen as a counterpart of the Fourier transform that is adapted to corner geometries. △ Less

Submitted 22 July, 2021; originally announced July 2021.

Comments: 19 pages, 2 figures

MSC Class: 45P05; 47G10; 47G30; 65R99

The Feshbach-Schur map and perturbation theory

Authors: Geneviève Dusson, Israel Sigal, Benjamin Stamm

Abstract: This paper deals with perturbation theory for discrete spectra of linear operators. To simplify exposition we consider here self-adjoint operators. This theory is based on the Feshbach-Schur map and it has advantages with respect to the standard perturbation theory in three aspects: (a) it readily produces rigorous estimates on eigenvalues and eigenfunctions with explicit constants; (b) it is comp… ▽ More This paper deals with perturbation theory for discrete spectra of linear operators. To simplify exposition we consider here self-adjoint operators. This theory is based on the Feshbach-Schur map and it has advantages with respect to the standard perturbation theory in three aspects: (a) it readily produces rigorous estimates on eigenvalues and eigenfunctions with explicit constants; (b) it is compact and elementary (it uses properties of norms and the fundamental theorem of algebra about solutions of polynomial equations); and (c) it is based on a self-contained formulation of a fixed point problem for the eigenvalues and eigenfunctions, allowing for easy iterations. We apply our abstract results to obtain rigorous bounds on the ground states of Helium-type ions. △ Less

Submitted 5 May, 2021; originally announced May 2021.

Subdiffusive-Brownian crossover in membrane proteins: a Generalized Langevin Equation-based approach

Authors: Loris Di Cairano, Benjamin Stamm, Vania Calandrini

Abstract: In this paper, we propose a Generalized Langevin Equation (GLE)-based model to describe the lateral diffusion of a protein in a lipid bilayer. The memory kernel is represented in terms of a viscous (instantaneous) and an elastic (non instantaneous) component modeled respectively through a Dirac delta function and a three-parameter Mittag-Leffler type function. By imposing a specific relationship b… ▽ More In this paper, we propose a Generalized Langevin Equation (GLE)-based model to describe the lateral diffusion of a protein in a lipid bilayer. The memory kernel is represented in terms of a viscous (instantaneous) and an elastic (non instantaneous) component modeled respectively through a Dirac delta function and a three-parameter Mittag-Leffler type function. By imposing a specific relationship between the parameters of the three-parameters Mittag-Leffler function, the different dynamical regimes, namely ballistic, subdiffusive and Brownian, as well as the crossover from one regime to another, are retrieved. Within this approach, the transition time from the ballistic to the subdiffusive regime and the distribution of relaxation times underlying the transition from the subdiffusive to the Brownian regime are given. The reliability of the model is tested by comparing the Mean Squared Displacement (MSD) derived in the framework of this model and the MSD of a protein diffusing in a membrane calculated through molecular dynamics (MD) simulations. △ Less

Submitted 10 August, 2021; v1 submitted 14 November, 2020; originally announced November 2020.

Analysis of the Feshbach-Schur method for the Fourier Spectral discretizations of Schr{ö}dinger operators

Authors: Geneviève Dusson, Israel Sigal, Benjamin Stamm

Abstract: In this article, we propose a new numerical method and its analysis to solve eigenvalue problems for self-adjoint Schr{ö}dinger operators, by combining the Feshbach-Schur perturbation theory with the spectral Fourier discretization. In order to analyze the method, we establish an abstract framework of Feshbach-Schur perturbation theory with minimal regularity assumptions on the potential that is… ▽ More In this article, we propose a new numerical method and its analysis to solve eigenvalue problems for self-adjoint Schr{ö}dinger operators, by combining the Feshbach-Schur perturbation theory with the spectral Fourier discretization. In order to analyze the method, we establish an abstract framework of Feshbach-Schur perturbation theory with minimal regularity assumptions on the potential that is then applied to the setting of the new spectral Fourier discretization method. Finally, we present some numerical results that underline the theoretical findings. △ Less

Submitted 8 December, 2021; v1 submitted 25 August, 2020; originally announced August 2020.

Guaranteed a posteriori bounds for eigenvalues and eigenvectors: multiplicities and clusters

Authors: Eric Cancès, Geneviève Dusson, Yvon Maday, Benjamin Stamm, Martin Vohralík

Abstract: This paper presents a posteriori error estimates for conforming numerical approximations of eigenvalue clusters of second-order self-adjoint elliptic linear operators with compact resolvent. Given a cluster of eigenvalues, we estimate the error in the sum of the eigenvalues, as well as the error in the eigenvectors represented through the density matrix, i.e., the orthogonal projector on the assoc… ▽ More This paper presents a posteriori error estimates for conforming numerical approximations of eigenvalue clusters of second-order self-adjoint elliptic linear operators with compact resolvent. Given a cluster of eigenvalues, we estimate the error in the sum of the eigenvalues, as well as the error in the eigenvectors represented through the density matrix, i.e., the orthogonal projector on the associated eigenspace. This allows us to deal with degenerate (multiple) eigenvalues within the framework. All the bounds are valid under the only assumption that the cluster is separated from the surrounding smaller and larger eigenvalues; we show how this assumption can be numerically checked. Our bounds are guaranteed and converge with the same speed as the exact errors. They can be turned into fully computable bounds as soon as an estimate on the dual norm of the residual is available, which is presented in two particular cases: the Laplace eigenvalue problem discretized with conforming finite elements, and a Schr{ö}dinger operator with periodic boundary conditions of the form $--$Δ$ + V$ discretized with planewaves. For these two cases, numerical illustrations are provided on a set of test problems. △ Less

Submitted 10 August, 2020; originally announced August 2020.

Journal ref: Mathematics of Computation, American Mathematical Society, In press

Analysis of the Schwarz domain decomposition method for the conductor-like screening continuum model

Authors: Arnold Reusken, Benjamin Stamm

Abstract: We study the Schwarz overlapping domain decomposition method applied to the Poisson problem on a special family of domains, which by construction consist of a union of a large number of fixed-size subdomains. These domains are motivated by applications in computational chemistry where the subdomains consist of van der Waals balls. As is usual in the theory of domain decomposition methods, the rate… ▽ More We study the Schwarz overlapping domain decomposition method applied to the Poisson problem on a special family of domains, which by construction consist of a union of a large number of fixed-size subdomains. These domains are motivated by applications in computational chemistry where the subdomains consist of van der Waals balls. As is usual in the theory of domain decomposition methods, the rate of convergence of the Schwarz method is related to a stable subspace decomposition. We derive such a stable decomposition for this family of domains and analyze how the stability "constant" depends on relevant geometric properties of the domain. For this, we introduce new descriptors that are used to formalize the geometry for the family of domains. We show how, for an increasing number of subdomains, the rate of convergence of the Schwarz method depends on specific local geometry descriptors and on one global geometry descriptor. The analysis also naturally provides lower bounds in terms of the descriptors for the smallest eigenvalue of the Laplace eigenvalue problem for this family of domains. △ Less

Submitted 30 June, 2021; v1 submitted 2 June, 2020; originally announced June 2020.

Comments: published in SIAM Journal on Numerical Analysis

MSC Class: 65N55

Journal ref: SIAM J. Numer. Anal., 59(2), 769-796, 2021

An approximation strategy to compute accurate initial density matrices for repeated self-consistent field calculations at different geometries

Authors: E. Polack, A. Mikhalev, Geneviève Dusson, B. Stamm, F. Lipparini

Abstract: Repeated computations on the same molecular system, but with different geometries, are often performed in quantum chemistry, for instance, in ab-initio molecular dynamics simulations or geometry optimizations. While many efficient strategies exist to provide a good guess for the self-consistent field procedure, which is usually the main computational task to be performed, little is known on how to… ▽ More Repeated computations on the same molecular system, but with different geometries, are often performed in quantum chemistry, for instance, in ab-initio molecular dynamics simulations or geometry optimizations. While many efficient strategies exist to provide a good guess for the self-consistent field procedure, which is usually the main computational task to be performed, little is known on how to efficiently exploit in this direction the abundance of information generated during the many computations. In this article, we present a strategy to provide an accurate initial guess for the density matrix, expanded in a set of localized basis functions, within the self-consistent field iterations for parametrized Hartree-Fock problems where the nuclear coordinates are changed along a few user-specified collective variables, such as the molecule's normal modes. Our approach is based on an offline-stage where the Hartree-Fock eigenvalue problem is solved for some particular parameter values and an online-stage where the initial guess is computed very efficiently for any new parameter value.The method allows non-linear approximations of density matrices, which belong to a non-linear manifold that is isomorphic to the Grassmann manifold.The so-called Grassmann exponential and logarithm map the manifold onto the tangent space and thus provides the correct geometrical setting accounting for the manifold structure when working with subspaces rather than functions itself.Numerical tests on different amino acids show promising initial results. △ Less

Submitted 5 June, 2020; v1 submitted 26 March, 2020; originally announced March 2020.

A Linear Scaling in Accuracy Numerical Method for Computing the Electrostatic Forces in the $N$-Body Dielectric Spheres Problem

Authors: Muhammad Hassan, Benjamin Stamm

Abstract: This article deals with the efficient and accurate computation of the electrostatic forces between charged, spherical dielectric particles undergoing mutual polarisation. We use the spectral Galerkin boundary integral equation framework developed by Lindgren et al. (J. Comput. Phys. 371 (2018): 712-731) and subsequently analysed in two earlier contributions of the authors to propose a linear scali… ▽ More This article deals with the efficient and accurate computation of the electrostatic forces between charged, spherical dielectric particles undergoing mutual polarisation. We use the spectral Galerkin boundary integral equation framework developed by Lindgren et al. (J. Comput. Phys. 371 (2018): 712-731) and subsequently analysed in two earlier contributions of the authors to propose a linear scaling in cost algorithm for the computation of the approximate forces. We establish exponential convergence of the method and derive error estimates for the approximate forces that do not explicitly depend on the number of dielectric particles $N$. Consequently, the proposed method requires only $\mathcal{O}(N)$ operations to compute the electrostatic forces acting on $N$ dielectric particles up to any given and fixed relative error. △ Less

Submitted 13 July, 2020; v1 submitted 4 February, 2020; originally announced February 2020.

Comments: 41 pages including Appendix and References; 10 figures

MSC Class: 65N12; 65N15; 65N35; 65R20

On the Scalability of the Parallel Schwarz Method in One-Dimension

Authors: Gabriele Ciaramella, Muhammad Hassan, Benjamin Stamm

Abstract: In contrast with classical Schwarz theory, recent results in computational chemistry have shown that for special domain geometries, the one-level parallel Schwarz method can be scalable. This property is not true in general, and the issue of quantifying the lack of scalability remains an open problem. Even though heuristic explanations are given in the literature, a rigorous and systematic analysi… ▽ More In contrast with classical Schwarz theory, recent results in computational chemistry have shown that for special domain geometries, the one-level parallel Schwarz method can be scalable. This property is not true in general, and the issue of quantifying the lack of scalability remains an open problem. Even though heuristic explanations are given in the literature, a rigorous and systematic analysis is still missing. In this short manuscript, we provide a first rigorous result that precisely quantifies the lack of scalability of the classical one-level parallel Schwarz method for the solution to the one-dimensional Laplace equation. Our analysis technique provides a possible roadmap for a systematic extension to more realistic problems in higher dimensions. △ Less

Submitted 23 December, 2019; originally announced December 2019.

Comments: To appear in the Proceedings of the 25th International Conference on Domain Decomposition Methods, St. John's, Canada, 2019

An Integral Equation Formulation of the $N$-body Dielectric Spheres Problem. Part II: Complexity Analysis

Authors: Bérenger Bramas, Muhammad Hassan, Benjamin Stamm

Abstract: This article is the second in a series of two papers concerning the mathematical study of a boundary integral equation of the second kind that describes the interaction of $N$ dielectric spherical particles undergoing mutual polarisation. The first article presented the numerical analysis of the Galerkin method used to solve this boundary integral equation and derived $N$-independent convergence r… ▽ More This article is the second in a series of two papers concerning the mathematical study of a boundary integral equation of the second kind that describes the interaction of $N$ dielectric spherical particles undergoing mutual polarisation. The first article presented the numerical analysis of the Galerkin method used to solve this boundary integral equation and derived $N$-independent convergence rates for the induced surface charges and total electrostatic energy. The current article will focus on computational aspects of the algorithm. We provide a convergence analysis of the iterative method used to solve the underlying linear system and show that the number of liner solver iterations required to obtain a solution is independent of $N$. Additionally, we present two linear scaling solution strategies for the computation of the approximate induced surface charges. Finally, we consider a series of numerical experiments designed to validate our theoretical results and explore the dependence of the numerical errors and computational cost of solving the underlying linear system on different system parameters. △ Less

Submitted 9 August, 2020; v1 submitted 17 November, 2019; originally announced November 2019.

Comments: 29 Page article

MSC Class: 65F35; 65N12; 65N15; 65N35; 65R20;

A Discontinuous Galerkin method for Shock Capturing using a mixed high-order and sub-grid low-order approximation space

Authors: Per-Olof Persson, Benjamin Stamm

Abstract: This article considers a new discretization scheme for conservation laws. The discretization setting is based on a discontinuous Galerkin scheme in combination with an approximation space that contains high-order polynomial modes as well as piece-wise constant modes on a sub-grid. The high-order modes can continuously be suppressed with a penalty function that is based on a sensor which is intertw… ▽ More This article considers a new discretization scheme for conservation laws. The discretization setting is based on a discontinuous Galerkin scheme in combination with an approximation space that contains high-order polynomial modes as well as piece-wise constant modes on a sub-grid. The high-order modes can continuously be suppressed with a penalty function that is based on a sensor which is intertwined with the approximation space. Numerical tests finally illustrate the performance of this scheme. △ Less

Submitted 13 October, 2021; v1 submitted 4 November, 2019; originally announced November 2019.

Gradient Flow Finite Element Discretizations with Energy-Based Adaptivity for the Gross-Pitaevskii Equation

Authors: Pascal Heid, Benjamin Stamm, Thomas P. Wihler

Abstract: We present an effective adaptive procedure for the numerical approximation of the steady-state Gross-Pitaevskii equation. Our approach is solely based on energy minimization, and consists of a combination of gradient flow iterations and adaptive finite element mesh refinements. Numerical tests show that this strategy is able to provide highly accurate results, with optimal convergence rates with r… ▽ More We present an effective adaptive procedure for the numerical approximation of the steady-state Gross-Pitaevskii equation. Our approach is solely based on energy minimization, and consists of a combination of gradient flow iterations and adaptive finite element mesh refinements. Numerical tests show that this strategy is able to provide highly accurate results, with optimal convergence rates with respect to the number of freedom. △ Less

Submitted 15 January, 2020; v1 submitted 17 June, 2019; originally announced June 2019.

MSC Class: 35P30; 47J25; 49M25; 49R05; 65N25; 65N30; 65N50

On the Scalability of the Schwarz Method

Authors: Gabriele Ciaramella, Muhammad Hassan, Benjamin Stamm

Abstract: In this article, we analyse the convergence behaviour and scalability properties of the one-level Parallel Schwarz method (PSM) for domain decomposition problems in which the boundaries of many subdomains lie in the interior of the global domain. Such problems arise, for instance, in solvation models in computational chemistry. Existing results on the scalability of the one-level PSM are limited t… ▽ More In this article, we analyse the convergence behaviour and scalability properties of the one-level Parallel Schwarz method (PSM) for domain decomposition problems in which the boundaries of many subdomains lie in the interior of the global domain. Such problems arise, for instance, in solvation models in computational chemistry. Existing results on the scalability of the one-level PSM are limited to situations where each subdomain has access to the external boundary, and at most only two subdomains have a common overlap. We develop a systematic framework that allows us to bound the norm of the Schwarz iteration operator for domain decomposition problems in which subdomains may be completely embedded in the interior of the global domain and an arbitrary number of subdomains may have a common overlap. △ Less

Submitted 18 October, 2019; v1 submitted 10 February, 2019; originally announced February 2019.

Comments: Revised Version of a previously uploaded article; 36 pages;

MSC Class: 65N55; 65F10; 65N22; 35J05; 35J57

Boundary integral equations for isotropic linear elasticity

Authors: Benjamin Stamm, Shuyang Xiang

Abstract: This articles first investigates boundary integral operators for the three-dimensional isotropic linear elasticity of a biphasic model with piecewise constant Lamé coefficients in the form of a bounded domain of arbitrary shape surrounded by a background material. In the simple case of a spherical inclusion, the vector spherical harmonics consist of eigenfunctions of the single and double layer bo… ▽ More This articles first investigates boundary integral operators for the three-dimensional isotropic linear elasticity of a biphasic model with piecewise constant Lamé coefficients in the form of a bounded domain of arbitrary shape surrounded by a background material. In the simple case of a spherical inclusion, the vector spherical harmonics consist of eigenfunctions of the single and double layer boundary operators and we provide their spectra. Further, in the case of many spherical inclusions with isotropic materials, each with its own set of Lamé parameters, we propose an integral equation and a subsequent Galerkin discretization using the vector spherical harmonics and apply the discretization to several numerical test cases. △ Less

Submitted 5 March, 2021; v1 submitted 6 February, 2019; originally announced February 2019.

Comments: 36 pages

An Integral Equation Formulation of the $N$-Body Dielectric Spheres Problem. Part I: Numerical Analysis

Authors: Muhammad Hassan, Benjamin Stamm

Abstract: In this article, we analyse an integral equation of the second kind that represents the solution of $N$ interacting dielectric spherical particles undergoing mutual polarisation. A traditional analysis can not quantify the scaling of the stability constants -- and thus the approximation error -- with respect to the number $N$ of involved dielectric spheres. We develop a new a priori error analysis… ▽ More In this article, we analyse an integral equation of the second kind that represents the solution of $N$ interacting dielectric spherical particles undergoing mutual polarisation. A traditional analysis can not quantify the scaling of the stability constants -- and thus the approximation error -- with respect to the number $N$ of involved dielectric spheres. We develop a new a priori error analysis that demonstrates $N$-independent stability of the continuous and discrete formulations of the integral equation. Consequently, we obtain convergence rates that are independent of $N$. △ Less

Submitted 13 July, 2020; v1 submitted 4 February, 2019; originally announced February 2019.

Comments: Revised version of the previously uploaded manuscript. 36 page article (38 including Appendix)

MSC Class: 65N12 (Primary) 65N15; 65N35; 65R20 (Secondary)

Theoretical analysis of screened many-body electrostatic interactions between charged polarizable particles

Authors: Eric B. Lindgren, Chaoyu Quan, Benjamin Stamm

Abstract: This paper builds on two previous works, Lindgren et al. J. Comp. Phys. 371, 712-731 (2018) and Quan et al. arXiv:1807.05384 (2018), to devise a new method to solve the problem of calculating electrostatic interactions in a system composed by many dielectric particles, embedded in a homogeneous dielectric medium, which in turn can also be permeated by charge carriers. The system is defined by the… ▽ More This paper builds on two previous works, Lindgren et al. J. Comp. Phys. 371, 712-731 (2018) and Quan et al. arXiv:1807.05384 (2018), to devise a new method to solve the problem of calculating electrostatic interactions in a system composed by many dielectric particles, embedded in a homogeneous dielectric medium, which in turn can also be permeated by charge carriers. The system is defined by the charge, size, position and dielectric constant of each particle, as well as the dielectric constant and Debye length of the medium. The effects of taking into account the dielectric nature of the particles is explored in selected scenarios where the presence of electrolytes in the medium can significantly influence the total undergoing interactions. Description of the mutual interactions between all particles in the system as being truly of many-body nature reveals how such effects can effectively influence the magnitudes and even directions of the resulting forces, especially those acting on particles that have a null net charge. Particular attention is given to a situation that can be related to colloidal particles in an electrolyte solution, where it's shown that polarization effects alone can substantially raise or lower---depending on the dielectric contrast between the particles and the medium---the energy barrier that divides particle coagulation and flocculation regions, when an interplay between electrostatic and additional van der Waals forces is considered. Overall, the results suggest that for an accurate description of the type of system in question, it is essential to consider particle polarization if the separation between the interacting particles are comparable to or smaller than the Debye length of the medium. △ Less

Submitted 31 October, 2018; originally announced October 2018.

Comments: 26 pages, 6 figures

An embedded corrector problem for homogenization. Part II: Algorithms and discretization

Authors: Eric Cancès, Virginie Ehrlacher, Frederic Legoll, Benjamin Stamm, Shuyang Xiang

Abstract: This contribution is the numerically oriented companion article of the work [E. Cancès, V. Ehrlacher, F. Legoll, B. Stamm and S. Xiang, arxiv preprint 1807.05131]. We focus here on the numerical resolution of the embedded corrector problem introduced in [E. Cancès, V. Ehrlacher, F. Legoll and B. Stamm, CRAS 2015; E. Cancès, V. Ehrlacher, F. Legoll, B. Stamm and S. Xiang, arxiv preprint 1807.05131]… ▽ More This contribution is the numerically oriented companion article of the work [E. Cancès, V. Ehrlacher, F. Legoll, B. Stamm and S. Xiang, arxiv preprint 1807.05131]. We focus here on the numerical resolution of the embedded corrector problem introduced in [E. Cancès, V. Ehrlacher, F. Legoll and B. Stamm, CRAS 2015; E. Cancès, V. Ehrlacher, F. Legoll, B. Stamm and S. Xiang, arxiv preprint 1807.05131] in the context of homogenization of diffusion equations. Our approach consists in considering a corrector-type problem, posed on the whole space, but with a diffusion matrix which is constant outside some bounded domain. In [E. Cancès, V. Ehrlacher, F. Legoll, B. Stamm and S. Xiang, arxiv preprint 1807.05131], we have shown how to define three approximate homogenized diffusion coefficients on the basis of the embedded corrector problems. We have also proved that these approximations all converge to the exact homogenized coefficients when the size of the bounded domain increases. We show here that, under the assumption that the diffusion matrix is piecewise constant, the corrector problem to solve can be recast as an integral equation. In case of spherical inclusions with isotropic materials, we explain how to efficiently discretize this integral equation using spherical harmonics, and how to use the fast multipole method (FMM) to compute the resulting matrix-vector products at a cost which scales only linearly with respect to the number of inclusions. Numerical tests illustrate the performance of our approach in various settings. △ Less

Submitted 23 October, 2018; originally announced October 2018.

A Domain Decomposition Method for the Poisson-Boltzmann Solvation Models

Authors: Chaoyu Quan, Benjamin Stamm, Yvon Maday

Abstract: In this paper, a domain decomposition method for the Poisson-Boltzmann (PB) solvation model that is widely used in computational chemistry is proposed. This method, called ddLPB for short, solves the linear Poisson-Boltzmann (LPB) equation defined in $\mathbb R^3$ using the van der Waals cavity as the solute cavity. The Schwarz domain decomposition method is used to formulate local problems by dec… ▽ More In this paper, a domain decomposition method for the Poisson-Boltzmann (PB) solvation model that is widely used in computational chemistry is proposed. This method, called ddLPB for short, solves the linear Poisson-Boltzmann (LPB) equation defined in $\mathbb R^3$ using the van der Waals cavity as the solute cavity. The Schwarz domain decomposition method is used to formulate local problems by decomposing the cavity into overlapping balls and only solving a set of coupled sub-equations in balls. A series of numerical experiments is presented to test the robustness and the efficiency of this method including the comparisons with some existing methods. We observe exponential convergence of the solvation energy with respect to the number of degrees of freedom which allows this method to reach the required level of accuracy when coupling with quantum mechanical descriptions of the solute. △ Less

Submitted 19 December, 2018; v1 submitted 14 July, 2018; originally announced July 2018.

An embedded corrector problem for homogenization. Part I: Theory

Authors: Eric Cancès, Virginie Ehrlacher, Frederic Legoll, Benjamin Stamm, Shuyang Xiang

Abstract: This article is the first part of a two-fold study, the objective of which is the theoretical analysis and numerical investigation of new approximate corrector problems in the context of stochastic homogenization. We present here three new alternatives for the approximation of the homogenized matrix for diffusion problems with highly-oscillatory coefficients. These different approximations all rel… ▽ More This article is the first part of a two-fold study, the objective of which is the theoretical analysis and numerical investigation of new approximate corrector problems in the context of stochastic homogenization. We present here three new alternatives for the approximation of the homogenized matrix for diffusion problems with highly-oscillatory coefficients. These different approximations all rely on the use of an embedded corrector problem (that we previously introduced in [Cancès, Ehrlacher, Legoll and Stamm, C. R. Acad. Sci. Paris, 2015]), where a finite-size domain made of the highly oscillatory material is embedded in a homogeneous infinite medium whose diffusion coefficients have to be appropriately determined. The motivation for considering such embedded corrector problems is made clear in the companion article [Cancès, Ehrlacher, Legoll, Stamm and Xiang, in preparation], where a very efficient algorithm is presented for the resolution of such problems for particular heterogeneous materials. In the present article, we prove that the three different approximations we introduce converge to the homogenized matrix of the medium when the size of the embedded domain goes to infinity. △ Less

Submitted 13 July, 2018; originally announced July 2018.

A coherent derivation of the Ewald summation for arbitrary orders of multipoles: The self-terms

Authors: Benjamin Stamm, Louis Lagardère, Étienne Polack, Yvon Maday, Jean-Philip Piquemal

Abstract: In this work, we provide the mathematical elements we think essential for a proper understanding of the calculus of the electrostatic energy of point-multipoles of arbitrary order under periodic boundary conditions. The emphasis is put on the expressions of the so-called self parts of the \Ewald\, summation where different expressions can be found in literature. Indeed, such expressions are of pri… ▽ More In this work, we provide the mathematical elements we think essential for a proper understanding of the calculus of the electrostatic energy of point-multipoles of arbitrary order under periodic boundary conditions. The emphasis is put on the expressions of the so-called self parts of the \Ewald\, summation where different expressions can be found in literature. Indeed, such expressions are of prime importance in the context of new generation polarizable force field where the self field appears in the polarization equations. We provide a general framework, where the idea of the \Ewald\ splitting is applied to the electric potential and subsequently, all other quantities such as the electric field, the energy and the forces are derived consistently thereof. Mathematical well-posedness is shown for all these contributions for any order of multipolar distribution. △ Less

Submitted 6 September, 2018; v1 submitted 25 May, 2018; originally announced May 2018.

Journal ref: The Journal of Chemical Physics,149, 124103 (2018)

A posteriori error estimates for discontinuous Galerkin methods using non-polynomial basis functions. Part II: Eigenvalue problems

Authors: Lin Lin, Benjamin Stamm

Abstract: We present the first systematic work for deriving a posteriori error estimates for general non-polynomial basis functions in an interior penalty discontinuous Galerkin (DG) formulation for solving eigenvalue problems associated with second order linear operators. Eigenvalue problems of such types play important roles in scientific and engineering applications, particularly in theoretical chemistry… ▽ More We present the first systematic work for deriving a posteriori error estimates for general non-polynomial basis functions in an interior penalty discontinuous Galerkin (DG) formulation for solving eigenvalue problems associated with second order linear operators. Eigenvalue problems of such types play important roles in scientific and engineering applications, particularly in theoretical chemistry, solid state physics and material science. Based on the framework developed in [{\it L. Lin, B. Stamm, http://dx.doi.org/10.1051/m2an/2015069}] for second order PDEs, we develop residual type upper and lower bound error estimates for measuring the a posteriori error for eigenvalue problems. The main merit of our method is that the method is parameter-free, in the sense that all but one solution-dependent constants appearing in the upper and lower bound estimates are explicitly computable by solving local and independent eigenvalue problems, and the only non-computable constant can be reasonably approximated by a computable one without affecting the overall effectiveness of the estimates in practice. Compared to the PDE case, we find that a posteriori error estimators for eigenvalue problems must neglect certain terms, which involves explicitly the exact eigenvalues or eigenfunctions that are not accessible in numerical simulations. We define such terms carefully, and justify numerically that the neglected terms are indeed numerically high order terms compared to the computable estimators. Numerical results for a variety of problems in 1D and 2D demonstrate that both the upper bound and lower bound are effective for measuring the error of eigenvalues and eigenfunctions. △ Less

Submitted 14 March, 2016; originally announced March 2016.

Comments: arXiv admin note: text overlap with arXiv:1502.01738

MSC Class: 65J10; 65N15; 65N30

A posteriori error estimates for discontinuous Galerkin methods using non-polynomial basis functions. Part I: Second order linear PDE

Authors: Lin Lin, Benjamin Stamm

Abstract: We present the first systematic work for deriving a posteriori error estimates for general non-polynomial basis functions in an interior penalty discontinuous Galerkin (DG) formulation for solving second order linear PDEs. Our residual type upper and lower bound error estimates measure the error in the energy norm. The main merit of our method is that the method is parameter-free, in the sense tha… ▽ More We present the first systematic work for deriving a posteriori error estimates for general non-polynomial basis functions in an interior penalty discontinuous Galerkin (DG) formulation for solving second order linear PDEs. Our residual type upper and lower bound error estimates measure the error in the energy norm. The main merit of our method is that the method is parameter-free, in the sense that all but one solution-dependent constants appearing in the upper and lower bound estimates are explicitly computable by solving local eigenvalue problems, and the only non-computable constant can be reasonably approximated by a computable one without affecting the overall effectiveness of the estimates in practice. As a side product of our formulation, the penalty parameter in the interior penalty formulation can be automatically determined as well. We develop an efficient numerical procedure to compute the error estimators. Numerical results for a variety of problems in 1D and 2D demonstrate that both the upper bound and lower bound are effective. △ Less

Submitted 17 June, 2015; v1 submitted 5 February, 2015; originally announced February 2015.

MSC Class: 65J10; 65N15; 65N30