Stochastic kinetic Monte Carlo algorithms for long-range Hamiltonians

We present a higher order kinetic Monte Carlo methodology suitable to model the evolution of systems in which the transition rates are non-trivial to calculate or in which Monte Carlo moves are likely to be non-productive flicker events. The second order residence time algorithm first introduced by Athènes et al. [Phil. Mag. A 76 (1997) 565] is rederived from the n-fold way algorithm of Bortz et al. [J. Comput. Phys. 17 (1975) 10] as a fully stochastic algorithm. The second order algorithm can be dynamically called when necessary to eliminate unproductive flickering between a metastable state and its neighbours. An algorithm combining elements of the first order and second order methods is shown to be more efficient, in terms of the number of rate calculations, than the first order or second order methods alone while remaining statistically identical. This efficiency is of prime importance when dealing with computationally expensive rate functions such as those arising from long-range Hamiltonians. Our algorithm has been developed for use when considering simulations of vacancy diffusion under the influence of elastic stress fields. We demonstrate the improved efficiency of the method over that of the n-fold way in simulations of vacancy diffusion in alloys. Our algorithm is seen to be an order of magnitude more efficient than the n-fold way in these simulations. We show that when magnesium is added to an Al-2at.%Cu alloy, this has the effect of trapping vacancies. When trapping occurs, we see that our algorithm performs thousands of events for each rate calculation performed.     

Parallel Monte Carlo simulations by asynchronous domain decomposition

A simple general method for performing Metropolis Monte Carlo condensed matter simulations on parallel processors is examined. The method is based on the cyclic generation of temporary discrete domains within the system, which are separated by distances greater than the inter-particle interaction range. Particle configurations within each domain are then sampled independently by an assigned processor, whilst particles outside these domains are held fixed. Results for a simulated Lennard-Jones fluid confirm that the method rigorously satisfies the detailed balance condition, and that the efficiency of configurational sampling scales almost linearly with the number of processors. Furthermore, the number of iterations performed on a given processor can be essentially arbitrary, with very low levels of inter-process communication. Provided the CPU time per step is not state-dependent, the method can then be used to perform large calculations as unsupervised background tasks on heterogeneous networks.     

Parallel Monte Carlo Driver (PMCD)—a software package for Monte Carlo simulations in parallel

Thanks to the dramatic decrease of computer costs and the no less dramatic increase in those same computer's capabilities and also thanks to the availability of specific free software and libraries that allow the set up of small parallel computation installations the scientific community is now in a position where parallel computation is within easy reach even to moderately budgeted research groups. The software package PMCD (Parallel Monte Carlo Driver) was developed to drive the Monte Carlo simulation of a wide range of user supplied models in parallel computation environments. The typical Monte Carlo simulation involves using a software implementation of a function to repeatedly generate function values. Typically these software implementations were developed for sequential runs. Our driver was developed to enable the run in parallel of the Monte Carlo simulation, with minimum changes to the original code that implements the function of interest to the researcher. In this communication we present the main goals and characteristics of our software, together with a simple study its expected performance. Monte Carlo simulations are informally classified as “embarrassingly parallel”, meaning that the gains in parallelizing a Monte Carlo run should be close to ideal, i.e. with speed ups close to linear. In this paper our simple study shows that without compromising the easiness of use and implementation, one can get performances very close to the ideal.     

Accelerating Monte Carlo simulations with an NVIDIA graphics processor

Modern graphics cards, commonly used in desktop computers, have evolved beyond a simple interface between processor and display to incorporate sophisticated calculation engines that can be applied to general purpose computing. The Monte Carlo algorithm for modelling photon transport in turbid media has been implemented on an NVIDIA^® 8800gt graphics card using the CUDA toolkit. The Monte Carlo method relies on following the trajectory of millions of photons through the sample, often taking hours or days to complete. The graphics-processor implementation, processing roughly 110 million scattering events per second, was found to run more than 70 times faster than a similar, single-threaded implementation on a 2.67 GHz desktop computer.

Program summary

Program title: Phoogle-C/Phoogle-GCatalogue identifier: AEEB_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEEB_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 51 264No. of bytes in distributed program, including test data, etc.: 2 238 805Distribution format: tar.gzProgramming language: C++Computer: Designed for Intel PCs. Phoogle-G requires a NVIDIA graphics card with support for CUDA 1.1Operating system: Windows XPHas the code been vectorised or parallelized?: Phoogle-G is written for SIMD architecturesRAM: 1 GBClassification: 21.1External routines: Charles Karney Random number library. Microsoft Foundation Class library. NVIDA CUDA library [1].Nature of problem: The Monte Carlo technique is an effective algorithm for exploring the propagation of light in turbid media. However, accurate results require tracing the path of many photons within the media. The independence of photons naturally lends the Monte Carlo technique to implementation on parallel architectures. Generally, parallel computing can be expensive, but recent advances in consumer grade graphics cards have opened the possibility of high-performance desktop parallel-computing.Solution method: In this pair of programmes we have implemented the Monte Carlo algorithm described by Prahl et al. [2] for photon transport in infinite scattering media to compare the performance of two readily accessible architectures: a standard desktop PC and a consumer grade graphics card from NVIDIA.Restrictions: The graphics card implementation uses single precision floating point numbers for all calculations. Only photon transport from an isotropic point-source is supported. The graphics-card version has no user interface. The simulation parameters must be set in the source code. The desktop version has a simple user interface; however some properties can only be accessed through an ActiveX client (such as Matlab).Additional comments: The random number library used has a LGPL (http://www.gnu.org/copyleft/lesser.html) licence.Running time: Runtime can range from minutes to months depending on the number of photons simulated and the optical properties of the medium.References:

[1]

http://www.nvidia.com/object/cuda_home.html.

[2]

S. Prahl, M. Keijzer, Sl. Jacques, A. Welch, SPIE Institute Series 5 (1989) 102.

     

Reduction of the self-forces in Monte Carlo simulations of semiconductor devices on unstructured meshes

When using an unstructured mesh for device geometry, the ensemble Monte Carlo simulations of semiconductor devices may be affected by unwanted self-forces resulting from the particle-mesh coupling. We report on the progress in minimisation of the self-forces on arbitrary meshes by showing that they can be greatly reduced on a finite element mesh with proper interpolation functions. The developed methodology is included into a self-consistent finite element 3D Monte Carlo device simulator. Minimising of the self-forces using the proper interpolation functions is tested by simulating the electron transport in a 10 nm gate length, 6.1 nm body thick, double gate metal-oxide-semiconductor field-effect transistor (MOSFET). We demonstrate the reduction in the self-force and illustrate the practical distinction by showing I-V characteristics for the device.     

An optimized algorithm for ionized impurity scattering in Monte Carlo simulations

We present a new optimized model of Brookes-Herring ionized impurity scattering for use in Monte Carlo simulations of semiconductors. When implemented, it greatly decreases the execution time needed for simulations (typically by a factor of the order of 100), and also properly incorporates the great proportion of small angle scatterings that are neglected in the standard algorithm. It achieves this performance by using an anisotropic choice of scattering angle which accurately mimics the true angular distribution of ionized impurity scattering.     

Non-reversible Monte Carlo simulations of spin models

Monte Carlo simulations are used to study simple systems where the underlying Markov chain satisfies the necessary condition of global balance but does not obey the more restrictive condition of detailed balance. Here, we show that non-reversible Markov chains can be set up that generate correct stationary distributions, but reduce or eliminate the diffusive motion in phase space typical of the usual Monte Carlo dynamics. Our approach is based on splitting the dynamics into a set of replicas with each replica representing a biased movement in reaction-coordinate space. This introduction of an additional bias in a given replica is compensated for by choosing an appropriate dynamics on the other replicas such as to ensure the validity of global balance. First, we apply this method to a mean-field Ising model, splitting the system into two replicas: one trying to increase magnetization and the other trying to decrease it. For this simple test system, our results show that the altered dynamics is able to reduce the dynamical critical exponent. Generalizations of this scheme to simulations of the Ising model in two dimensions are discussed.     

Large scale atomistic polymer simulations using Monte Carlo methods for parallel vector processors

In this paper we discuss the implementation of advanced variable connectivity Monte Carlo (MC) simulation methods for studying large (>10⁵ atom) polymer systems at the atomic level. Such codes are intrinsically difficult to optimize since they involve a mixture of many different elementary MC steps, such as reptation, flip, end rotation, concerted rotation and volume fluctuation moves. In particular, connectivity altering MC moves, such as the recently developed directed end bridging (DEB) algorithm, are required in order to vigorously sample the configuration space. Techniques for effective vector implementation of such moves are described. We also show how a simple domain decomposition method can provide a general and efficient means of parallelizing these complex MC protocols. Benchmarks are reported for a 192,000 atom simulation of polydisperse linear polyethylene with an average chain length C₆₀₀₀, for simulations using 1 to 8 processors and a variety of MC protocols.     

Monte Carlo simulations of the site-diluted three-dimensional XY model

We present extensive Monte Carlo simulations of the classical XY model in three dimensions on a simple cubic lattice with periodic boundary conditions subject to quenched, random site dilution. By using a hybrid algorithm and finite-size scaling techniques we estimate the transition temperatures and the critical exponents associated with the diluted system. Our results for the critical exponents and universal cumulants are, within statistical errors, the same as for the zero dilution case and are independent of the amount of dilution, in agreement with the Harris criterion. The initial reduction of the critical temperature with dilution is also evaluated and compared to recent experimental results.     

Distributed implementation of the adaptive kinetic Monte Carlo method

The program EON2 is a distributed implementation of the adaptive kinetic Monte Carlo method for long time scale simulations of atomistic systems. The method is based on the transition state theory approach within the harmonic approximation and the key step is the identification of relevant saddle points on the potential energy rim surrounding the energy minimum corresponding to a state of the system. The saddle point searches are carried out in a distributed fashion starting with random initial displacements of the atoms in regions where atoms have less than optimal coordination. The main priorities of this implementation have been to (1) make the code transparent, (2) decouple the master and slaves, and (3) have a well defined interface to the energy and force evaluation. The computationally intensive parts are implemented in C++, whereas the less compute intensive server-side software is written in Python. The platform for distributed computing is BOINC. A simulation of the annealing of a twist and tilt grain boundary in a copper crystal is described as an example application.     

Error in Monte Carlo, quasi-error in Quasi-Monte Carlo

While the Quasi-Monte Carlo method of numerical integration achieves smaller integration error than standard Monte Carlo, its use in particle physics phenomenology has been hindered by the absence of a reliable way to estimate that error. The standard Monte Carlo error estimator relies on the assumption that the points are generated independently of each other and, therefore, fails to account for the error improvement advertised by the Quasi-Monte Carlo method. We advocate the construction of an estimator of stochastic nature, based on the ensemble of pointsets with a particular discrepancy value. We investigate the consequences of this choice and give some first empirical results on the suggested estimators.     

A grid procedure applied to the determination of parameters of a kinetic process

A fitting procedure for one trap and one recombination centre kinetic model is described here. The procedure makes use of a grid in the parameters space obtained by changing each parameter back and forth and calculating robust cost functions on the surfaces of this grid. The lengths of the changes are determined empirically. The best set of parameters is calculated by the projection on the grid surface with smallest cost function. The fitting procedure applied to the fit of one, two and three parameters of the kinetic model is analyzed. In all cases the optimization procedure shows reliable fitting within a feasible interval of processing time.     

Static and dynamic glass transitions in the 10-state Potts glass: What can Monte Carlo simulations contribute?

The p-state Potts glass with infinite range Gaussian interactions can be solved exactly in the thermodynamic limit and exhibits an unconventional phase behavior if p>4: A dynamical transition from ergodic to non-ergodic behavior at a temperature T_D is followed by a first order transition at T₀<T_D, where a glass order parameter appears discontinuously, although the latent heat is zero. If one assumes that a similar scenario occurs for the structural glass transition as well (though with the singular behavior at T_D rounded off), the p-state Potts glass should be a good test case to develop methods to deal with finite size effects for the static as well as the dynamic transition, and to check what remnants of these unconventional transitions are left in finite sized systems, as they are used in simulations. While it is shown that a sensible extrapolation N→∞ of the simulation results are compatible with the exact results, we find that it would be rather difficult to obtain a correct understanding of the behavior of the system in the thermodynamic limit if only the numerical data would be available.     

Monte Carlo simulations of systems with complex energy landscapes

Non-traditional Monte Carlo simulations are a powerful approach to the study of systems with complex energy landscapes. After reviewing several of these specialized algorithms we shall describe the behavior of typical systems including spin glasses, lattice proteins, and models for “real” proteins. In the Edwards-Anderson spin glass it is now possible to produce probability distributions in the canonical ensemble and thermodynamic results of high numerical quality. In the hydrophobic-polar (HP) lattice protein model Wang-Landau sampling with an improved move set (pull-moves) produces results of very high quality. These can be compared with the results of other methods of statistical physics. A more realistic membrane protein model for Glycophorin A is also examined. Wang-Landau sampling allows the study of the dimerization process including an elucidation of the nature of the process.     

Tests of the multi-spin-coding technique in Monte Carlo simulations of statistical systems

We present the multi-spin-coding technique applied to the Monte Carlo simulation of the Ising model in some detail, including a FORTRAN program, and confirm its efficiency.     

Mixed Monte Carlo/Molecular Dynamics simulations of the prion protein

In this paper we present the results of mixed Monte Carlo/Molecular Dynamics (MC/MD) simulations of the D178N mutant of the human prion protein. We have used the MC moves for polypeptide sampling known as Concerted Rotations with Angles (CRA) to selectively sample the region of the prion protein comprising the β-sheet and one of the α-helices. The results indicate that the MC/MD simulations sample the phase space substantially faster than regular Molecular Dynamics simulations starting with the same initial conditions. This work further indicates the MC/MD technique as a potentially powerful simulation tool, allowing the selective sampling of a region of a physical system that is deemed important.     

Monte Carlo simulations of a single polystyrene chain in spherical confinement

We report on Monte Carlo simulations of a single coarse-grained polystyrene chain in spherical confinement. To this end we employ a variant of the freely rotating chain model, the parameters of which are chosen to mimic polystyrene in good solvent conditions. Entanglements are analyzed as a function of molecular weight and capsid radius to provide an educated guess about the structure of a single polystyrene chain in a miniemulsion droplet. We also show that significant knotting occurs first when the radius of the confining sphere falls below the chain?s radius of gyration.