Publications

A new scheme is presented for imposing periodic boundary conditions on unit cells with arbitrary source distributions. We restrict our attention here to the Poisson, modified Helmholtz, Stokes and modified Stokes equations. The approach extends to the oscillatory equations of mathematical physics, including the Helmholtz and Maxwell equations, but we will address these in a companion paper, since the nature of the problem is somewhat different and includes the consideration of quasiperiodic boundary conditions and resonances. Unlike lattice sum-based methods, the scheme is insensitive to the unit cell's aspect ratio and is easily coupled to adaptive fast multipole methods (FMMs). Our analysis relies on classical "plane-wave" representations of the fundamental solution, and yields an explicit low-rank representation of the field due to all image sources beyond the first layer of neighboring unit cells. When the aspect ratio of the unit cell is large, our scheme can be coupled with the nonuniform fast Fourier transform (NUFFT) to accelerate the evaluation of the induced field. Its performance is illustrated with several numerial examples.

We introduce libdlr, a library implementing the recently introduced discrete Lehmann representation (DLR) of imaginary time Green's functions. The DLR basis consists of a collection of exponentials chosen by the interpolative decomposition to ensure stable and efficient recovery of Green's functions from imaginary time or Matsbuara frequency samples. The library provides subroutines to build the DLR basis and grids, and to carry out various standard operations. The simplicity of the DLR makes it straightforward to incorporate into existing codes as a replacement for less efficient representations of imaginary time Green's functions, and libdlr is intended to facilitate this process. libdlr is written in Fortran, and contains a Python module pydlr. We also introduce a stand-alone Julia implementation, Lehmann.jl.

In this paper we analyze the dynamics of single-excitation states, which model the scattering of a single photon from multiple two level atoms. For short times and weak atom-field couplings we show that the atomic amplitudes are given by a sum of decaying exponentials, where the decay rates and Lamb shifts are given by the poles of a certain analytic function. This result is a refinement of the "pole approximation" appearing in the standard Wigner-Weisskopf analysis of spontaneous emission. On the other hand, at large times, the atomic field decays like O(1/t3) with a known constant expressed in terms of the coupling parameter and the resonant frequency of the atoms. Moreover, we show that for stronger coupling, the solutions also feature a collection of oscillatory exponentials which dominate the behavior at long times. Finally, we extend the analysis to the continuum limit in which atoms are distributed according to a given density.

We consider the numerical solution of the real time equilibrium Dyson equation, which is used in calculations of the dynamical properties of quantum many-body systems. We show that this equation can be written as a system of coupled, nonlinear, convolutional Volterra integro-differential equations, for which the kernel depends self-consistently on the solution. As is typical in the numerical solution of Volterra-type equations, the computational bottleneck is the quadratic-scaling cost of history integration. However, the structure of the nonlinear Volterra integral operator precludes the use of standard fast algorithms. We propose a quasilinear-scaling FFT-based algorithm which respects the structure of the nonlinear integral operator. The resulting method can reach large propagation times, and is thus well-suited to explore quantum many-body phenomena at low energy scales. We demonstrate the solver with two standard model systems: the Bethe graph, and the Sachdev-Ye-Kitaev model.

Abstract Several antimicrobial peptides, including magainin and the human cathelicidin LL-37, act by forming pores in bacterial membranes. Bacteria such as Staphylococcus aureus modify their membrane's cardiolipin composition to resist such types of perturbations that compromise their membrane stability. Here, we used molecular dynamics simulations to quantify the role of cardiolipin on the formation of pores in simple bacterial-like membrane models composed of phosphatidylglycerol and cardiolipin mixtures. Cardiolopin modified the structure and ordering of the lipid bilayer, making it less susceptible to mechanical changes. Accordingly, the free-energy barrier for the formation of a transmembrane pore and its kinetic instability augmented by increasing the cardiolipin concentration. This is attributed to the unfavorable positioning of cardiolipin near the formed pore, due to its small polar-head and bulky hydrophobic-body. Overall, our study demonstrates how cardiolipin prevents membrane-pore formation and this constitutes a plausible mechanism used by bacteria to act against stress perturbations and, thereby, gain resistance to antimicrobial agents.

The efficiency of Hamiltonian Monte Carlo (HMC) can suffer when sampling a distribution with a wide range of length scales, because the small step sizes needed for stability in high-curvature regions are inefficient elsewhere. To address this we present a delayed rejection variant: if an initial HMC trajectory is rejected, we make one or more subsequent proposals each using a step size geometrically smaller than the last. We extend the standard delayed rejection framework by allowing the probability of a retry to depend on the probability of accepting the previous proposal. We test the scheme in several sampling tasks, including multiscale model distributions such as Neal's funnel, and statistical applications. Delayed rejection enables up to five-fold performance gains over optimally-tuned HMC, as measured by effective sample size per gradient evaluation. Even for simpler distributions, delayed rejection provides increased robustness to step size misspecification. Along the way, we provide an accessible but rigorous review of detailed balance for HMC.

The long-term goal of this work is the development of high-fidelity simulation tools for dispersive tsunami propagation. A dispersive model is especially important for short wavelength phenomena such as an asteroid impact into the ocean, and is also important in modeling other events where the simpler shallow water equations are insufficient. Adaptive simulations are crucial to bridge the scales from deep ocean to inundation, but have difficulties with the implicit system of equations that results from dispersive models. We propose a fractional step scheme that advances the solution on separate patches with different spatial resolutions and time steps. We show a simulation with 7 levels of adaptive meshes and onshore inundation resulting from a simulated asteroid impact off the coast of Washington. Finally, we discuss a number of open research questions that need to be resolved for high quality simulations.

Stacking is a widely used model averaging technique that asymptotically yields optimal predictions among linear averages. We show that stacking is most effective when model predictive performance is heterogeneous in inputs, and we can further improve the stacked mixture with a hierarchical model. We generalize stacking to Bayesian hierarchical stacking. The model weights are varying as a function of data, partially-pooled, and inferred using Bayesian inference. We further incorporate discrete and continuous inputs, other structured priors, and time series and longitudinal data. To verify the performance gain of the proposed method, we derive theory bounds, and demonstrate on several applied problems.

Simulations with an adaptive time-dependent bias, such as metadynamics, enable an efficient exploration of the conformational space of a system. However, the dynamic information of the system is altered by the bias. With infrequent metadynamics it is possible to recover the transition rate of crossing a barrier, if the collective variables are ideal and there is no bias deposition near the transition state. Unfortunately, for simulations of complex molecules, these conditions are not always fulfilled. To overcome these limitations, and inspired by single-molecule force spectroscopy, we developed a method based on Kramers' theory for calculating the barrier-crossing rate when a time-dependent bias is added to the system. We assess the quality of the bias parameter by measuring how efficiently the bias accelerates the transitions compared to ideal behavior. We present approximate analytical expressions of the survival probability that accurately reproduce the barrier-crossing time statistics, and enable the extraction of the unbiased transition rate even for challenging cases, where previous methods fail.

Well-conditioned boundary integral methods for the solution of elliptic boundary value problems (BVPs) are powerful tools for static and dynamic physical simulations. When there are many close-to-touching boundaries (eg, in complex fluids) or when the solution is needed in the bulk, nearly-singular integrals must be evaluated at many targets. We show that precomputing a linear map from surface density to an effective source representation renders this task highly efficient, in the common case where each object is "simple", ie, its smooth boundary needs only moderately many nodes. We present a kernel-independent method needing only an upsampled smooth surface quadrature, and one dense factorization, for each distinct shape. No (near-)singular quadrature rules are needed. The resulting effective sources are drop-in compatible with fast algorithms, with no local corrections nor bookkeeping. Our extensive numerical tests include 2D FMM-based Helmholtz and Stokes BVPs with up to 1000 objects (281000 unknowns), and a 3D Laplace BVP with 10 ellipsoids separated by 1/30 of a diameter. We include a rigorous analysis for analytic data in 2D and 3D.