Cristian Predescu

The incomplete beta function law for parallel tempering sampling of classical canonical systems

We show that the acceptance probability for swaps in the parallel tempering Monte Carlo method for classical canonical systems is given by a universal function that depends on the average statistical fluctuations of the potential and on the ratio of the temperatures. The law, called the incomplete beta function law, is valid in the limit that the two temperatures involved in swaps are close to one another. An empirical version of the law, which involves the heat capacity of the system, is developed and tested on a Lennard-Jones cluster. We argue that the best initial guess for the distribution of intermediate temperatures for parallel tempering is a geometric progression and we also propose a technique for the computation of optimal temperature schedules. Finally, we demonstrate that the swap efficiency of the parallel tempering method for condensed-phase systems decreases naturally to zero at least as fast as the inverse square root of the dimensionality of the physical system.

Thermodynamics and equilibrium structure of Ne38 cluster: quantum mechanics versus classical.

The equilibrium properties of classical Lennard-Jones (LJ38) versus quantum Ne38 Lennard-Jones clusters are investigated. The quantum simulations use both the path-integral Monte Carlo (PIMC) and the recently developed variational-Gaussian wave packet Monte Carlo (VGW-MC) methods. The PIMC and the classical MC simulations are implemented in the parallel tempering framework. The classical heat capacity Cv(T) curve agrees well with that of Neirotti et al. [J. Chem. Phys. 112, 10340 (2000)], although a much larger confining sphere is used in the present work. The classical Cv(T) shows a peak at about 6 K, interpreted as a solid-liquid transition, and a shoulder at approximately 4 K, attributed to a solid-solid transition involving structures from the global octahedral (Oh) minimum and the main icosahedral (C5v) minimum. The VGW method is used to locate and characterize the low energy states of Ne38, which are then further refined by PIMC calculations. Unlike the classical case, the ground state of Ne38 is a liquidlike structure. Among the several liquidlike states with energies below the two symmetric states (Oh and C5v), the lowest two exhibit strong delocalization over basins associated with at least two classical local minima. Because the symmetric structures do not play an essential role in the thermodynamics of Ne38, the quantum heat capacity is a featureless curve indicative of the absence of any structural transformations. Good agreement between the two methods, VGW and PIMC, is obtained. The present results are also consistent with the predictions by Calvo et al. [J. Chem. Phys. 114, 7312 (2001)] based on the quantum superposition method within the harmonic approximation. However, because of its approximate nature, the latter method leads to an incorrect assignment of the Ne38 ground state as well as to a significant underestimation of the heat capacity.

Heat capacity estimators for random series path-integral methods by finite-difference schemes

Previous heat capacity estimators used in path integral simulations either have large variances that grow to infinity with the number of path variables or require the evaluation of first- and second-order derivatives of the potential. In the present paper, we show that the evaluation of the total energy by the T-method estimator and of the heat capacity by the TT-method estimator can be implemented by a finite difference scheme in a stable fashion. As such, the variances of the resulting estimators are finite and the evaluation of the estimators requires the potential function only. By comparison with the task of computing the partition function, the evaluation of the estimators requires k+1 times more calls to the potential, where k is the order of the difference scheme employed. Quantum Monte Carlo simulations for the Ne13 cluster demonstrate that a second order central-difference scheme should suffice for most applications.

Accurate and efficient integration for molecular dynamics simulations at constant temperature and pressure

In molecular dynamics simulations, control over temperature and pressure is typically achieved by augmenting the original system with additional dynamical variables to create a thermostat and a barostat, respectively. These variables generally evolve on timescales much longer than those of particle motion, but typical integrator implementations update the additional variables along with the particle positions and momenta at each time step. We present a framework that replaces the traditional integration procedure with separate barostat, thermostat, and Newtonian particle motion updates, allowing thermostat and barostat updates to be applied infrequently. Such infrequent updates provide a particularly substantial performance advantage for simulations parallelized across many computer processors, because thermostat and barostat updates typically require communication among all processors. Infrequent updates can also improve accuracy by alleviating certain sources of error associated with limited-precision arithmetic. In addition, separating the barostat, thermostat, and particle motion update steps reduces certain truncation errors, bringing the time-average pressure closer to its target value. Finally, this framework, which we have implemented on both general-purpose and special-purpose hardware, reduces software complexity and improves software modularity.

Optimal series representations for numerical path integral simulations

By means of the Ito-Nisio theorem, we introduce and discuss a general approach to series representations of path integrals. We then argue that the optimal basis for both “primitive” and partially averaged approaches is the Wiener sine-Fourier basis. The present analysis also suggests a new approach to improving the convergence of primitive path integral methods. Current work indicates that this new technique, the “reweighted” method, converges as the cube of the number of path variables for “smooth” potentials. The technique is based on a special way of approximating the Brownian bridge which enters the Feynman-Kac formula and it does not require the Gaussian transform of the potential for its implementation.

Phase changes in selected Lennard-Jones X13−nYn clusters

Detailed studies of the thermodynamic properties of selected binary Lennard-Jones clusters of the type X13-nYn (where n=1, 2, 3) are presented. The total energy, heat capacity, and first derivative of the heat capacity as a function of temperature are calculated by using the classical and path integral Monte Carlo methods combined with the parallel tempering technique. A modification in the phase change phenomena from the presence of impurity atoms and quantum effects is investigated.

Numerical implementation of some reweighted path integral methods

The reweighted random series techniques provide finite-dimensional approximations to the quantum density matrix of a physical system that have fast asymptotic convergence. We study two special reweighted techniques that are based upon the Levy–Ciesielski and Wiener–Fourier series, respectively. In agreement with the theoretical predictions, we demonstrate by numerical examples that the asymptotic convergence of the two reweighted methods is cubic for smooth enough potentials. For each reweighted technique, we propose some minimalist quadrature techniques for the computation of the path averages. These quadrature techniques are designed to preserve the asymptotic convergence of the original methods.

Equipartition and the Calculation of Temperature in Biomolecular Simulations

Since the behavior of biomolecules can be sensitive to temperature, the ability to accurately calculate and control the temperature in molecular dynamics (MD) simulations is important. Standard analysis of equilibrium MD simulations-even constant-energy simulations with negligible long-term energy drift-often yields different calculated temperatures for different motions, however, in apparent violation of the statistical mechanical principle of equipartition of energy. Although such analysis provides a valuable warning that other simulation artifacts may exist, it leaves the actual value of the temperature uncertain. We observe that Tolmans generalized equipartition theorem should hold for long stable simulations performed using velocity-Verlet or other symplectic integrators, because the simulated trajectory is thought to sample almost exactly from a continuous trajectory generated by a shadow Hamiltonian. From this we conclude that all motions should share a single simulation temperature, and we provide a new temperature estimator that we test numerically in simulations of a diatomic fluid and of a solvated protein. Apparent temperature variations between different motions observed using standard estimators do indeed disappear when using the new estimator. We use our estimator to better understand how thermostats and barostats can exacerbate integration errors. In particular, we find that with large (albeit widely used) time steps, the common practice of using two thermostats to remedy so-called hot solvent-cold solute problems can have the counterintuitive effect of causing temperature imbalances. Our results, moreover, highlight the utility of multiple-time step integrators for accurate and efficient simulation.

Energy estimators for random series path-integral methods

We perform a thorough analysis on the choice of estimators for random series path integral methods. In particular, we show that both the thermodynamic (T-method) and the direct (H-method) energy estimators have finite variances and are straightforward to implement. It is demonstrated that the agreement between the T-method and the H-method estimators provides an important consistency check on the quality of the path integral simulations. We illustrate the behavior of the various estimators by computing the total, kinetic, and potential energies of a molecular hydrogen cluster using three different path integral techniques. Statistical tests are employed to validate the sampling strategy adopted as well as to measure the performance of the parallel random number generator utilized in the Monte Carlo simulation. Some issues raised by previous simulations of the hydrogen cluster are clarified.

Reconstruction of silicon surfaces: A stochastic optimization problem

Over the last two decades, scanning tunneling microscopy (STM) has become one of the most important ways to investigate the structure of crystal surfaces. STM has helped achieve remarkable successes in surface science such as finding the atomic structure of Si(111) and Si(001). For high-index Si surfaces the information about the local density of states obtained by scanning does not translate directly into knowledge about the positions of atoms at the surface. A commonly accepted strategy for identifying the atomic structure is to propose several possible models and analyze their corresponding simulated STM images for a match with the experimental ones. However, the number of good candidates for the lowest-energy structure is very large for high-index surfaces, and heuristic approaches are not likely to cover all the relevant structural models. In this paper, we take the view that finding the atomic structure of a surface is a problem of stochastic optimization, and we address it as such. We design a general technique for predicting the reconstruction of silicon surfaces with arbitrary orientation, which is based on parallel-tempering Monte Carlo simulations combined with an exponential cooling. The advantages of the method are illustrated using the Si(105) surface as an example, with two main results: (a) the correct single-step rebonded structure [e.g., Fujikawa, Akiyama, Nagao, Sakurai, Lagally, Hashimoto, Morikawa, and Terakura, Phys. Rev. Lett. 88, 176101 (2002)] is obtained even when starting from the paired-dimer model [Mo, Savage, Swartzentruber, and Lagally, Phys. Rev. Lett. 65, 1020 (1990)] that was assumed to be correct for many years, and (b) we have found several double-step reconstructions that have lower surface energies than any previously proposed double-step models.