Session Z49: Computational Methods for Statistical Mechanics: Advances and Applications I

Parallel flat-histogram method for equilibrium and non-equilibrium problems

I will present a parallel implementation of the multicanonical method and its generalization to a flat-histogram method suitable for other equilibrium ensembles as well as non-equilibrium problems. After a tutorial introduction to the method, I will present several successful applications to diverse equilibrium problems, both on- and off-lattice, and new applications to spreading processes.   

Microcanonical Inflection-Point Analysis of Phases for Semiflexible Polymers

We study a generic model of semiflexible polymer with self-interactions, which is known to exhibit a multitude of structural phases. Previous studies employing canonical statistical analysis methods for the identification and characterization of these phases have been inconclusive as these approaches lead to inconsistent results for systems of finite size. For our study, we use the recently introduced generalized microcanonical inflection-point analysis method [1]. This analysis technique was developed to enable the systematic identification and classification of transitions in systems of any size. In this talk, we discuss the structural hyperphase diagram of semiflexible polymers in the range of comparatively large bending stiffness, which extends a previous case study [2]. Extensive replica-exchange Monte Carlo simulations were performed to obtain accurate estimates of the microcanonical entropies for the different values of the bending stiffness.

References

[1] K. Qi and M. Bachmann, Phys. Rev. Lett. 120, 180601 (2018).

[2] D. Aierken and M. Bachmann, Polymers 12, 3013 (2020).   

Statistical Mechanics of Chemical Disorder and Magnetic Order in MnSb2Te4 and MnBi2Te4

The magnetic ordering in the layered topological insulator materials MnSb₂Te₄ and MnBi₂Te₄ has been observed to depend strongly on the disorder and defects on the Mn and Sb/Bi sublattices. Here we will explore the effect of this chemical disorder on the magnetic interactions and we will present models of the magnetic interactions in these materials that have been extracted from large scale first principles density functional theory calculations. The first principles calculations consider unbiased disorder realizations to capture the magnetic interactions both inside the planes as well as between planes based on the chemical environment. These models form the basis for our Monte-Carlo simulations where we identify the magnetic transitions between the ferromagnetic and antiferromagetic states as a function of disorder as well as the transition temperatures for these magnetic orders.   

A graph theory-based statistical mechanics approach for nucleation of nano-porous materials

Understanding the nucleation of weak electrolytes from solution is critical for the design and synthesis of crystalline nano-porous materials such as metal-organic frameworks (MOFs) and zeolites. However, existing simulation approaches to model nucleation are often extremely limited when applied to weak electrolytes. We developed a novel graph theory-based sampling approach that overcomes limitations of existing approaches, especially, for bulk solvent systems. Our method seeks to exploit the property of materials whose crystal structure exhibit directional bonding and thus can be described as a "graph" of connected monomers. By utilizing a rigorous statistical-mechanics approach, we generate a nucleus of size N+1 from size N by performing a Monte Carlo type attachment of a solute to the surface of a nucleus followed by a thermodynamic integration step to turn on interactions. A "bootstrapping" approach in the nucleus size (N) is then used to generate an ensemble of representative nuclei and their corresponding free energies. To validate our approach, we begin with an ideal system of non interacting particles and evaluate the free energy dependence with nucleus size. The free energies show excellent agreement with reference values obtained from a GCMC (Grand Canonical Monte Carlo) approach. Subsequent work will extended our approach to treat systems with more complex interactions and geometries.   

Testing and Usage of an Automated Multi-bead Iterative Boltzmann Inversion (IBI) Code

Iterative Boltzmann Inversion (IBI) is a systematic coarse-graining (CG) method in which tabular CG potentials are derived that reproduce target distributions generated from atomistic reference simulations. We report here on a software code which automates the development of multi-bead coarse-grained potentials using IBI. Two major problems make automation difficult: 1) noisy distributions derived from sampling; 2) low sampling regions which introduces discontinuities in the sampling. Both of these make differentiation to calculate energy and forces difficult and error prone. Our code addresses these problems using an approach which combines data smoothing, fitting to functional expansions, and extrapolation schemes to handle low sampling regions without discontinuity or data distortion. A number of code features will be discussed aimed at providing guidance on usage. In particular, we discuss the importance of proper sampling and equilibration in the reference simulations and how this leads accurate target distributions in low sampling regions (as well as the suppression of erroneous features) which is a key to potential convergence. Finally, we discuss associated tools that enable robust input including automated conversion of molecular systems from all-atom (AA) to CG models.   

Fitting Hamiltonian Parameters with Bayesian Optimization of Numerical Simulations

Tuning simulation models to accurately reproduce experimental observations is a challenging but vital task for predicting material properties and informing further experiments. The search for a set of suitable Hamiltonian parameters should ideally minimize the number of dispatched simulations and be fully automated to maximize throughput. We describe a framework that utilizes Bayesian optimization[1], which fits surrogate models predicting the expected difference between target observables and results of numerical simulations, to intelligently explore and optimize in parameter space. The optimization scheme is robust in application and can function without relying on any gradient information. We apply this methodology to Monte Carlo simulations of a coarse-grained binary alloy model to successfully recover a reference interspecies coupling calculated from first principles methods[2]. The Hamiltonian fitting procedure is then demonstrated for experimental data and dynamical simulations of continuous spin models for magnetic materials.

[1] H. J. Kushner. J. Basic Eng. 86, (1964); J. Močkus. J. Glob. Optim. 4, (1975); D. Jones, M. Schonlau, W. Welch. J. Glob. Optim. 13, (1998).

[2] S. Curtarolo et al., Comp. Mat. Sci. 58 218 (2012).   

Avoiding critical slowdown in models with SALR interactions

Particles with competing short-range attractive and long-range repulsive (SALR) interactions can form a broad array of microphase morphologies. Given that structural richness, minimal models solved by Monte Carlo methods help to hone in on the underlying physics. However, even at weak frustration, configurational sampling in the vicinity of order-disorder transition temperature, $T_c$, is particularly inefficient. Standard cluster algorithms, such as the Swendsen-Wang and Wolff schemes, then fail because they generate clusters that don't capture physical correlations and even percolate at $T>T_c$. Alternate approaches have hence long been sought out. In this presentation, we present a mean-field analysis of this challenge and use our findings to propose an algorithmic approach that sidesteps these difficulties.   

Asymptotic Error Analysis of the MBAR Equations

 The efficient estimation of high-dimensional integrals is central to answering many questions in statistical mechanics. An important class of sampling strategies, including  umbrella sampling and alchemical methods, involves sampling from multiple thermodynamic states. A statistically optimal estimator named Multistate Bennett Acceptance Ratio (MBAR) has been developed to systematically combine data from all the sampled states to estimate free energies and other ensemble averages over another probability distribution. However, the error of the MBAR estimator is not well-understood: previous error analysis of MBAR assumed independent samples and only gave the total error without tracing the error to individual sampled states. In this work, We derive a central limit theorem for the estimates from MBAR, which allows for error decomposition into individual sampled states. We demonstrate the error estimator for a two-dimensional umbrella sampling of the alanine dipeptide and the alchemical free energy calculation of the solvation of methane in water. Our error analysis suggests a close connection between error contribution of each state and the autocorrelation in the samples from that state. The preliminary results suggest that assessing the error contributions provides insight into the sources of error of the simulations and can be used to improve the accuracy of MBAR calculations.   

Some Numerical Issues in Polymer Density-Functional Theories for Tangent Hard-Sphere Chains

The two most popular polymer density-functional theories (PDFTs), one proposed by Yu and Wu (J. Chem. Phys. 117, 2368, 2002) and the other by Chapman and co-workers (J. Chem. Phys. 127, 244904, 2007), are both applicable only to model systems based on tangent hard-sphere chains. Due to the hard-sphere repulsion, however, many integrals involved in these PDFTs contain discontinuous and/or indifferentiable integrands, and cannot be accurately evaluated by the straightforward application of any quadrature or Fourier transform. Here we show how to solve this problem in both real-space and reciprocal-space calculations, using the later versions of both PDFTs (J. Chem. Phys. 118, 3835, 2003; J. Chem. Theory Comput. 8, 1393, 2012). We also use an iteration method that can converge the PDFT equations to much higher accuracy than the commonly used Picard iteration (also known as the simple mixing or relaxation method). These allow us to obtain highly accurate results (where the maximum absolute value of the residual errors is less than 10^-10) in much less iteration steps.   

Loop-Cluster Coupling and Algorithm for Classical Statistical Models

Potts spin systems play a fundamental role in statistical mechanics and quantum field theory and can be studied within the spin, the Fortuin–Kasteleyn (FK) bond or the q-flow (loop) representation.We introduce a Loop-Cluster (LC) joint model of bond-occupation variables interacting with q-flow variables and formulate an LC algorithm that is found to be in the same dynamical universality as the celebrated Swendsen–Wang algorithm. This leads to a theoretical unification for all the representations, and numerically, one can apply the most efficient algorithm in one representation and measure physical quantities in others. Moreover, by using the LC scheme, we construct a hierarchy of geometric objects that contain as special cases the q-flow clusters and the backbone of FK clusters, the exact values of whose fractal dimensions in two dimensions remain as an open question. Our work not only provides a unified framework and an efficient algorithm for the Potts model but also brings new insights into the rich geometric structures of the FK clusters.   

Funnel Hopping Monte Carlo: Efficient Monte Carlo simulations for multi-funnel systems

Monte Carlo simulations are a popular tool used in physics, chemistry and biology to study the behavior of atomic systems at finite temperatures.

However, in many cases, the potential energy landscape is divided into multiple funnels separated by high energy barriers, which can prevent a complete sampling within available simulation time.

Our newly developed method, named funnel hopping Monte Carlo (FHMC), overcomes this problem by introducing non-local Monte Carlo moves from one low energy region to another, bypassing the high energy barriers entirely.

These new moves are generated by sampling Gaussian mixtures, that have been fit to the Boltzmann distributions of each funnel in advance to the simulation.

The method was first tested on the 38 and 75 atom Lennard-Jones clusters, for which conventional Monte Carlo simulations fail to converge, due to their double funnel energy landscape.

By integrating our method into the parallel tempering scheme, we were able to reduce the number of energy and force evaluations needed for a convergence of the simulations significantly.

Recently we also extended the method to periodic systems for which we will present results from ongoing applications.   

Langevin Dynamics Simulations of Dielectric Response of Correlated Materials

We have developed a fast and flexible computational scheme to calculate the complex valued, frequency dependent dielectric function ε(ω)  of correlated materials. The premise of such a methodology is to use the atomistic crystal structure of materials and designate relatively simple bond length and bond-angle interactions, as well as internal field couplings to accommodate correlation. Once these interactions are appropriated, we simulate systems of  N x N x M unit cells with periodic boundary conditions via Langevin dynamics under an oscillating external electric field. We validate our method by recreating high-resoltuion infrared near-field experimental results of the dielectric function of perovskite oxide SmNiO₃. We show quantitative agreement with experimental data concerning the modulation of the nanoscale dielectric changes as a function of hydrogen doping and the prevalence of oxygen vacancies within the sample. We also show representative example of how our methodology can gain us nanoscale insight into the dielectric behavior of Moire patterned two-dimensional transition metal dichalcogenides and compare our results with high-resolution near-field measurements.