Session H21: Advances in Computational Methods for Statistical Physics and Their Applications III

Population Annealing for Equilibrium Sampling in Statistical Physics

Population annealing is an algorithm that is well-suited for sampling equilibrium distributions for systems such as spin glasses and configurational glasses with rough free energy landscapes. Population annealing is massively parallel and easily implemented either on large clusters or GPUs. In addition to being massively parallel, population annealing has several attractive features: (1) it gives direct access to thermodynamic potentials, (2) multiple independent runs can be combined using weighted averaging to improve both statistical and systematic errors and (3) equilibration can be assessed with an intrinsic measure. In this talk population annealing will be introduced and put in the context of other sequential Monte Carlo and annealing algorithms. Applications to spin glasses, configurational glasses and first order transitions will be discussed.   

Population annealing: A massively parallel algorithm for simulating systems with rough free energy landscapes

Population annealing is a massively parallel sequential Monte Carlo algorithm that is designed to sample equilibrium states of systems with rough free energy landscapes. Population annealing is closely related to simulated annealing. In population annealing, the relevant equilibrium ensemble is represented by a large population of replicas of the system. The population is initialized in an easy to equilibrate region of parameter space and is then "annealed" in parameter space to a more difficult, target region. Resampling the population of replicas during each annealing step ensures that the population remains near equilibrium. The entropy and thermodynamic potentials together with intrinsic estimates of systematic errors are readily accessible from the simulations. We will describe both canonical and microcanonical versions of the population annealing algorithm, and discuss applications to spin glasses, binary hard sphere fluids in the glassy regime, and large q Potts models.   

Quantum Free Energy Differences from Non-Equilibrium Path Integral Methods

In this work, we discuss how the imaginary-time path integral representation of the quantum canonical partition function and non-equilibrium work fluctuation relations can be combined to yield methods for computing free energy differences in quantum systems using non-equilibrium processes. The path integral representation is isomorphic to the configurational partition function of a classical field theory to which a natural Hamiltonian dynamics can be associated. It is then shown that both, Jarzynski nonequilibrium work relation and Crooks fluctuation relation, formally hold for this classical field theory. Since the energy diverges in canonical equilibrium, regularization methods need to be introduced in order to limit the number of degrees of freedom M to be finite. The convergence of the work distribution as M tends to infinity is demonstrated analytically for a system composed of a quantum particle trapped in a harmonic well, and numerically for a quartic double-well potential with varying asymmetry. Finally, the method is used to study the relevance of protonic quantum effects in ionic water clusters.   

Stable Recursion Relation for the Canonical Partition Function of Non-Interacting Fermions

While the thermodynamics of a system of non-interacting fermions can be straightforwardly determined in the grand canonical ensemble, results for a specific number of N particles are more difficult to obtain. In this talk we will present a general recursion scheme for the canonical partition function of free fermions with a quadratic energy level spacing as might be present for ultracold fermionic atoms confined inside a box trap. Exact results for the entropy, specific heat and one and two particle occupation probabilities are numerically obtained to arbitrary precision and compared with their corresponding grand canonical values. The numerical stability of the recursion relation allows us to quantify deviations from Wicks theorem in the canonical ensemble for a variety of temperatures and densities.   

New Approaches to Tensor Network Simulation of 2D Quantum Systems

Accessing generic many-body quantum ground states in two full dimensions remains an outstanding problem in numerical condensed matter physics. Quantum Monte Carlo suffers from the sign problem, exact diagonalization cannot access large systems, DMRG struggles to reach beyond quasi-1D ladders, and many existing 2D tensor network approaches still have much room for improvement in terms of computational costs and numerical stability. We develop a PEPS approach, inspired by successful DMRG methods in one spatial dimension, which has much better performance than DMRG while retaining many of its advantages. We discuss the benefits and some drawbacks of the method, and we present some preliminary successes in fully two dimensional simulations with an eye to working towards simulating models on the frontier of numerical and analytical understanding.   

Accelerating Monte Carlo Simulations of Two-Dimensional Spin Models using GPUs

Utilizing the vast computational power of Graphics Processing Units, we develop novel algorithms in both OpenCL and CUDA to study the statistical mechanical properties of two-dimensional spin models. Our basic technique is the general Markov chain Monte Carlo method which uses a Metropolis-Hastings step. In contrast to the popular method of checkerboard updates, we have devised an unconventional procedure that traverses rows to generate the spin configurations. This special implementation allows us to simulate large system sizes (~ 10^9 spins) and, when tested on our fastest GPU, produces an average time per spin-flip of 0.005 ns with a speed-up factor of 600 compared to an optimized single-core CPU algorithm. We test the performance characteristics of our techniques for simulating 2D spin models such as the Ising, XY, and Potts on hexagonal and square lattices. Finally, as a theoretical proof of our computational concept we present critical temperatures of the models based on finite-size scaling methods.   

Accelerated simulation of gelation

The formation of gels is a complex issue that has to be resolved to investigate manifold synthetic materials - among them: porous materials such as cement, high-quality glass fiber and geomaterials for radioactive waste sealing.

Being able to simulate the structural, mechanical and dynamical properties of gels would have far-reaching consequences to improve the synthetization of materials and to lengthen their lifetime. This requires rather demanding computational power, which is worsen by the fact that gels have weak long-distance ordering (which requires large simulation cells) and that the gelation is a slow process due to high energy barrier for bond formation.

In order to reproduce the gelation we coupled GCMC and Parallel tempering methods to enhance the formation of silicate chains and reduce the computational costs. Using such cost-efficient method, we were able to study the formation of alkali-silica gels that is key to understand concrete structures failure.   

Heuristic optimization and sampling with tensor networks for quasi-2D spin glass problems

We devise a deterministic classical algorithm to reveal the structure of low energy spectrum for certain spin-glass systems that encode classical optimization problems. We employ tensor networks to represent probability distributions of all possible configurations. We then develop efficient techniques for approximately extract the relevant information from the networks for a class of quasi-two-dimensional Ising Hamiltonians. To this end, we apply a branch and bound approach over marginal probability distributions by approximately evaluating tensor contractions. Our approach identifies configurations with the largest Boltzmann weights corresponding to low energy states. We discover spin-glass droplet structures at finite temperatures, by exploiting local nature of the problems. This droplet finding algorithm naturally encompass sampling from high quality solutions within a given approximation ratio. It is, thus, established that tensor networks techniques can provide profound insight into the structure of large low-dimension spin-glass problems, with ramifications both for machine learning and noisy intermediate-scale quantum devices. Morever, limitations of our approach highlight alternative directions to establish quantum speed-up and possible quantum supremacy experiments.   

Studying the finite temperature properties of ferroelectrics

ABO₃ perovskite ferroelectrics, as well as their solid solutions, exhibit rich transitional behavior patterns that can be exploited, e.g., to obtain large piezoelectric and dielectric responses. Due to the complexity of the phase diagrams of these materials, mesoscale-level parameterizations capable of accurately reproducing their finite temperature properties are difficult to develop. Furthermore, obtaining such parameters from first principles calculations is complicated by the large number of available exchange correlation (XC) functionals. We investigate the influence of XC functionals on the prediction accuracy of ferroelectric phase transitions in PbTiO₃. LDA, PBE, PBEsol and vdW-DF-C09 XC functionals are evaluated utilizing constant-temperature molecular dynamics, in comparison with Wang-Landau (WL) Monte Carlo and Replica Exchange WL simulations. We find that LDA, PBEsol and vdW-C09 provide good estimates of physical properties near the phase transition, as compared with experiments, while PBE significantly overestimates the transition temperature.   

Predicting the surface phase diagram of Ag(111) using ab initio grand canonical Monte Carlo

The structure of a surface can dramatically affect its properties. For example, surface reconstructions can occur that change band alignments and/or catalytic activity. Currently, ab initio thermodynamics is the method of choice for determining the stable surfaces of a material, however, the selection of surfaces to study is done manually, which induces bias that can prevent one from finding global minima in the surface energy. We present an implementation of ab initio grand canonical Monte Carlo (GCMC) that automatically predicts surface phase diagrams and apply it to the Ag(111) system. We obtain an Ag₇O₁₀ overlayer, which is consistent with the most stable reconstruction found experimentally and computationally. We extracted structure-stability trends from our simulation data using machine learning and find that surface coordination and bond angles are important descriptors for stability. We analyzed the stochastic evolution of the surface and discovered a possible mechanism for the formation of the Ag₇O₁₀ overlayer. Ab initio GCMC therefore offers a rich set of possibilities for studying interfacial systems.   

A Lattice Boltzmann Method for Simulating Dry and Dense Active Fluids

Symmetry serves a foundational role in all areas of physics today. In classical many-body problems, by first clarifying the underlying symmetries of a system of interest, one can derive the hydrodynamic equations of motion that govern the dynamics of the system. Analysis of a hydrodynamic theory can elucidate the universal behavior exhibited by all generic systems respecting the prescribed set of symmetries; conversely, any particular many-body system defined by microscopic rules that respect the same set of symmetries can also be used to study the associated universal behavior in the hydrodynamic limit. An example of the latter is the use of lattice gas cellular automata to study the Navier-Stokes equations. Superseding the lattice gas cellular automata is the celebrated lattice Boltzmann method, which led to a drastic improvement in computational efficiency. Surprisingly, the development of a lattice Boltzmann method for dry active fluids is still lacking, which is what we accomplish here. We will demonstrate the usefulness of our approach by clarifying the phase behaviour of polar active fluids and motility-induced phase separation. In particular, we show that there are generically three distinct phases in polar active fluids separated by two discontinuous phase transitions.   

Spin-flop transition in the 3D anisotropic Heisenberg antiferromagnet: Finite size scaling for a first order transition where a continuous symmetry is broken

We use Monte Carlo simulations to explore the 3D anisotropic Heisenberg antiferromagnet in a field in order to study the finite size behavior of the first order “spin-flop” transition between the Ising-like antiferromagnetic state and the canted, XY-like state[1]. Finite size scaling for a first order phase transition where a continuous symmetry is broken is developed using an approximation of Gaussian probability distributions with a phenomenological “degeneracy” factor, q, included. Our theory yields q = π, and it predicts that for large linear dimension L the field dependence of all moments of the order parameters as well as the fourth-order cumulants exhibit universal intersections. The values of these intersections at the spin-flop transition point can be expressed in terms of the factor q. Our theory and simulation imply a heretofore unknown universality can be invoked for first order phase transitions.

[1] S. Hu, S.-H. Tsai, and D. P. Landau, Phys. Rev. E 89, 032118 (2014)   

Evaluating the Jones polynomial with tensor networks

We introduce tensor network contraction algorithms for the evaluation of the Jones polynomial of arbitrary knots. The value of the Jones polynomial of a knot maps to the partition function of a q-state Potts model defined as a planar graph with weighted edges that corresponds to the knot. For any integer q, we cast this partition function into tensor network form and employ fast tensor network contraction protocols to obtain the exact tensor trace, and thus the value of the Jones polynomial. By sampling random knots via a grid-walk procedure and computing the full tensor trace, we demonstrate numerically that the Jones polynomial can be evaluated in time that scales subexponentially with the number of crossings in the typical case. This allows us to evaluate the Jones polynomial of knots that are too complex to be treated with other available methods. Our results establish tensor network methods as a practical tool for the study of knots.