Parallel Monte Carlo Simulations Research Articles

Parallel pseudo-random number generators: A derivative pricing perspective with the Heston stochastic volatility model

Accuracy and precision of parallel Monte Carlo (MC) simulations may be impaired by the presence of intra-thread and inter-thread correlations depending on the parallel pseudo-random number generators (PPRNGs) used. While necessary, statistical tests alone are insufficient to ensure the absence of these correlations that can manifest as bias and variance to a extent in different applications. Therefore, application-based tests designed to mimic real-life MC scenarios may uncover them in the intended applications. The results of an application-based test simulating the Heston stochastic volatility model, a widely used pricing framework, is reported in order to compare the accuracy and precision profiles among four popular libraries of scalable pseudo-random number generators implementing sequence division ( trng and RngSteam ), parameterization ( sprng ) and counter-based ( Random123 ) strategies. All pseudo-random number generators assessed demonstrate similar standard-error of mean profiles. However, the bias profiles are more varied albeit comparable. PPRNGs demonstrating the smallest bias profiles in absolute and relative terms are yarn4 from TRNG , mlfg from SPRNG , as well as Threefry2x64 from Random123 for truncated Euler scheme, and mrg5s from TRNG and lfg from SPRNG for the quadratic exponent scheme. It is recommended that, when selecting a PPRNG for parallel MC simulation from a set of competing PPRNGs with comparable bias and standard error of mean profiles in absolute terms, the PPRNG associated with the smallest parallel-sequential bias difference should be used as it reflects the smallest intra-thread correlation introduced by parallelization.

Scientific computing with Google App Engine

We investigate the potential of Google App Engine (GAE) for scientific computing. We designed a generic master–slave framework that enables implementation and integration of new algorithms by instantiating one interface and two abstract classes. Applications are dynamically scheduled, executed, and monitored by the framework on the Google infrastructure, including a fault tolerant resubmission mechanism upon certain failures. We describe the implementation of a parallel rank sort algorithm using our framework and analyse the speedup, the execution overheads, and the cost of executing it on the Google Cloud infrastructure. Experimental results show that important speedup can be obtained from using GAE, especially for embarrassingly parallel Monte Carlo simulations. The two main obstacles that hinder the performance of applications on GAE are the resource quotas upon free use (especially the 30 s sequential execution time limit) and the middleware overheads. Cost-wise, GAE offers opportunities for cheaper computation due to the CPU cycle-based payment granularity, as opposed to the hourly billing intervals of other providers such as Amazon EC2.

MCMC2 : A Monte Carlo code for multiply-charged clusters

We have implemented a Monte Carlo code in reduced units providing structural and thermodynamical properties of multiply-charged Lennard-Jones droplets ANn+, composed of N individual particles among which n are charged particles, each one carrying a charge qi (qi can be positive or negative). The cluster has a total net charge Q=∑i=1,nqi (Q>0 or Q<0). The interactions between particles are modelled by a sum of pairwise Lennard-Jones potentials and electrostatic terms, including polarisation. The cluster statistical properties can be obtained from (i) parallel Monte Carlo simulations whose replicas are run at different temperatures, from configurations with same number of charged particles n and same individual charges qi (Parallel Tempering Monte Carlo), or (ii) parallel Monte Carlo simulations whose replicas are run at the same temperature but from configurations with different qi or n (Parallel Charging Monte Carlo). The code provides statistical data (evaporation rates, acceptance/rejection rates, etc.), energetic data (mean energies, heat capacities, etc.), and structural data (radial and angular distributions), for comprehensive analyses. A complete manual of the code is provided. Program summaryProgram title: MCMC2Catalogue identifier: AENZ_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AENZ_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.: 146790No. of bytes in distributed program, including test data, etc.: 1501715Distribution format: tar.gzProgramming language: Fortran 90 with MPI extensions for parallelisation.Computer: x86 and IBM platforms.Operating system:1.CentOS 5.6 Intel Xeon X5670 2.93 GHz, gfortran +MPICH2.2.CentOS 5.3 Intel Xeon E5520 2.27 GHz, gfortran/g95/pgf90 +MPICH2.3.Red Hat Enterprise 5.3 Intel Xeon X5650 2.67 GHz, gfortran +IntelMPI.4.IBM Power 6 4.7 GHz, xlf +PESS (IBM parallel library).Has the code been vectorised or parallelised?: Yes, parallelised using MPI extensions. Number of CPUs used: 2 to ∼40 but the code enables the use of 999 CPUs if desired.RAM: 10–20 MBClassification: 23.Nature of problem:We provide a general parallel code to investigate structural and thermodynamic properties of multiply charged clusters.Solution method:Parallel Monte Carlo methods are implemented for the exploration of the configuration space of multiply charged clusters. Two parallel Monte Carlo methods were found appropriate to achieve such a goal: the Parallel Tempering method, where replicas of the same cluster at different temperatures are distributed among different CPUs, and Parallel Charging where replicas (at the same temperature) having different particle charges or numbers of charged particles are distributed on different CPUs.Restrictions:The current version of the code uses Lennard-Jones interactions, as the main cohesive interaction between spherical particles, and electrostatic interactions (charge–charge, charge-induced dipole, induced dipole–induced dipole, polarisation). The Monte Carlo simulations can only be performed in the NVT ensemble in the present code.Unusual features:The Parallel Charging methods, based on the same philosophy as Parallel Tempering but with particle charges and number of charged particles as parameters instead of temperature, is an interesting new approach to explore energy landscapes. Splitting of the simulations is allowed and averages are accordingly updated.Running time:The running time depends on the number of Monte Carlo steps, cluster size, and type of interactions selected (e.g., polarisation turned on or off, and method used for calculating the induced dipoles). Typically a complete simulation can last from a few tens of minutes or a few hours for small clusters (N≤100, not including polarisation interactions), to one week for large clusters (N≥1000 not including polarisation interactions), and several weeks for large clusters (N≥1000) when including polarisation interactions. A restart procedure has been implemented that enables a splitting of the simulation accumulation phase.

Confinement Effects on Water Clusters Inside Carbon Nanotubes

The effects of confinement on water clusters inside nonmetallic carbon nanotubes with radii ranging between 4 and 7.5 Å have been computationally investigated by means of global optimization and finite temperature simulations. The water–water interaction is described by the TIP4P rigid body potential, and a Lennard-Jones potential is used for the water–carbon interaction. Water clusters containing up to 20 molecules are found to form 1D chainlike configurations for the narrow (7, 5) nanotube and 2D ladderlike structures in the (7, 6) tube. In wider tubes, 3D configurations are then formed showing helical motifs, ringlike or closed cage structures, before the most stable structure on flat graphene is eventually found. The same results are obtained by replacing the fully atomistic water–nanotube potential by its continuous approximation [Bretón, J.; González-Platas, J.; Giradet, C. J. Chem. Phys.1994, 101, 3334], indicating a negligible effect of corrugation. The effects of additional nanotubes were also considered with the adsorption energies being found to converge rather quickly already for the triple-wall tube. Parallel tempering Monte Carlo simulations of the water octamer reveal a counterintuitive decrease in the melting point relative to the free-standing case. Molecular dynamics simulations show that melting is concomitant with some axial diffusion of the water molecules, and with radial diffusion perpendicular to the tube axis remaining limited. In accordance with previous studies concerned with bulk water, the weakening of the cluster thermal stability is interpreted as being caused by the hydrophobic character of the carbon–water interaction.

Transillumination Imaging for Blood Oxygen Saturation Estimation of Skin Lesions

Detecting the early stages of melanoma can be greatly assisted by an accurate estimate of subsurface blood volume and blood oxygen saturation, indicative of angiogenesis. Visualization of this blood volume present beneath a skin lesion can be achieved through the transillumination of the skin. As the absorption of major chromophores in the skin is wavelength dependent, multispectral imaging can provide the needed information to separate out relative amounts of each chromophore. However, a critical challenge to this strategy is relating the pixel intensities observed in a given image to the wavelength-dependent total absorption existing at each spatial location. Consequently, in this paper, we develop an extension to Beer's law, estimated through a novel voxel-based, parallel processing Monte Carlo simulation of light propagation in skin which takes into account the specific geometry of our transillumination imaging apparatus. We then use this relation in a linear mixing model, solved using a multispectral image set, for chromophore separation and oxygen saturation estimation of an absorbing object located at a given depth within the medium. Validation is performed through the Monte Carlo simulation, as well as by imaging on a skin phantom. Results show that subsurface oxygen saturation can be reasonably estimated with good implications for the reconstruction of 3-D skin lesion volumes using transillumination toward early detection of malignancy.

3-D Volume Reconstruction of Skin Lesions for Melanin and Blood Volume Estimation and Lesion Severity Analysis

Subsurface information about skin lesions, such as the blood volume beneath the lesion, is important for the analysis of lesion severity towards early detection of skin cancer such as malignant melanoma. Depth information can be obtained from diffuse reflectance based multispectral transillumination images of the skin. An inverse volume reconstruction method is presented which uses a genetic algorithm optimization procedure with a novel population initialization routine and nudge operator based on the multispectral images to reconstruct the melanin and blood layer volume components. Forward model evaluation for fitness calculation is performed using a parallel processing voxel-based Monte Carlo simulation of light in skin. Reconstruction results for simulated lesions show excellent volume accuracy. Preliminary validation is also done using a set of 14 clinical lesions, categorized into lesion severity by an expert dermatologist. Using two features, the average blood layer thickness and the ratio of blood volume to total lesion volume, the lesions can be classified into mild and moderate/severe classes with 100% accuracy. The method therefore has excellent potential for detection and analysis of pre-malignant lesions.

Highly parallel Monte-Carlo simulations of the acousto-optic effect in heterogeneous turbid media

The development of a highly parallel simulation of the acousto-optic effect is detailed. The simulation supports optically heterogeneous simulation domains under insonification by arbitrary monochromatic ultrasound fields. An adjoint method for acousto-optics is proposed to permit point-source/point-detector simulations. The flexibility and efficiency of this simulation code is demonstrated in the development of spatial absorption sensitivity maps which are in broad agreement with current experimental investigations. The simulation code has the potential to provide guidance in the feasibility and optimization of future studies of the acousto-optic technique, and its speed may permit its use as part of an iterative inversion model.

Some Considerations in Devising Effective SCALE6/MAVRIC Models for Large Shielding Applications

*Analog Monte Carlo simulations for deep-penetration shielding problems are not viable; instead, aggressive use of variance reduction methods is used aiming to achieve acceptable accuracy for the result of interest in acceptable time. This requires developing effective variance reduction parameters, which in turn requires a systematic approach and an automated method. Frequently, a hybrid deterministic-stochastic methodology is deployed, where a deterministic transport theory method is used to generate adjoint function distribution. In the FW-CADIS method implemented in MAVRIC sequence of the SCALE6 package, deterministic forward calculation is performed as well, to enable generating variance reduction parameters for Monte Carlo simulation resulting in nearly-uniform statistical uncertainty over a large region. The choice of the deterministic model is user-controlled. However, the accuracy of the forward and adjoint deterministic solutions impacts the quality of Monte Carlo variance reduction parameters. We investigate this problem on a complex real-life shielding problem. The objective is to determine the radiation field (dose rate) throughout a nuclear power plant building. The large problem size (about 50x50x50 m 3 ), is further complicated by attenuation of more than 15 orders of magnitude. The results will serve as the initial step toward developing a general methodology for performance optimization of hybrid shielding calculations.

Thermodynamics of polymer adsorption to a flexible membrane

We analyze the structural behavior of a single polymer chain grafted to an attractive, flexible surface. Our model is composed of a coarse-grained bead-and-spring polymer and a tethered membrane. By means of extensive parallel tempering Monte Carlo simulations it is shown that the system exhibits a rich phase behavior ranging from highly ordered, compact to extended random coil structures, and from desorbed to completely adsorbed or even partially embedded conformations. These findings are summarized in a pseudophase diagram indicating the predominant class of conformations as a function of the external parameters temperature and polymer-membrane interaction strength. By comparison with adsorption to a stiff membrane surface it is shown that the flexibility of the membrane gives rise to qualitatively new behavior such as stretching of adsorbed conformations.

SU‐E‐T‐688: McCloud: Toward 10 Million Monte Carlo Primaries in 5 Minutes for Clinical Use

Purpose: At AAPM 2010, we presented the concept of using Cloud Computing to speed up clinical Monte Carlo calculations. Cloud Computing refers to a set of technologies that offers on‐demand computing resources using a pay‐as‐you‐go pricing model. We presented a feasibility study that included a variety of clinical relevant calculations using protons, electrons, and photons in water and a voxel patient geometry (Zubal phantom). While our initial study was designed for proof‐of‐concept, this research presents a more refined platform with an emphasis on particle therapy and demonstrates its potential to enable Monte Carlo dose calculation for routine clinical use.Methods: A Cloud Computing based distributed Monte Carlo dose calculation system called McCloud was developed. The system is based on Amazon EC2, Fluka, and a variety of “glue” technologies. The calculation involves these key steps: (1) A cluster is allocated in the Cloud. (2) A Monte Carlo task is launched on each node in the cluster. Extra care is used for random seed generations to ensure the correctness of the simulations. (3) The progress of the computation is dynamically monitored. (4) Once the computation is finished in the cluster, the results are aggregated using a linear model or a tree‐based distributive model, and a final dose distribution is returned to the user. Results: We found that for 14 Million 75MeV protons, an analog Monte Carlo simulation can be completed in 5 minutes 36 seconds. For this particular computation, the cost is less than 2 USD (based on official pricing at 8 cents per node per hour). Studies using various cluster sizes have shown the 1/n speed up expected from analytical Monte Carlo simulations Conclusions: Our implementations demonstrate the power of the existing cloud computing model, an extremely promising new paradigm for clinical and research computing. Further research is underway.Supported in part by NSF CBET‐0853157 and NSF CBET‐0755054.

Charge flow model for atomic ordering in nonisovalent alloys

Nonisovalent alloys, also known as aliovalent alloys, are formed by mixing two semiconductors of different valences for both cations and anions. These alloys exhibit many interesting properties and have found applications in optoelectronics, refractory materials, photovoltaics, photocatalytic splitting of water, etc. For example, the alloy of GaN and ZnO has a surprisingly large band gap bowing, enabling visible light absorption and overall water splitting at a record quantum efficiency. The understanding of the properties of nonisovalent alloys, however, is hindered by the lack of knowledge of the detailed atomic structures. We recently developed a charge flow model which can predict the total energy of different atomic configurations of nonisovalent alloys. In this work, we extend this model and apply it to a number of alloy systems---GaN/ZnO, GaAs/ZnSe, InP/CdS, AlN/ZnO, AlN/MgO, and AlN/SiC in wurtzite, zinc blende, and/or rocksalt structures. Good agreement between the model-predicted and ab initio results is found. We also employ the charge flow model in parallel tempering Monte Carlo simulations to calculate the thermodynamic properties of nonisovalent alloys of wurtzite structure. A phase transition between the phase separated and alloying regions is found and a general phase diagram is obtained. This model could be used to guide the design and synthesis of new nonisovalent alloy materials.

Molecular Simulations of Sequence-Specific Association of Transmembrane Proteins in Lipid Bilayers

Association of membrane proteins is central in material and information flow across the cellular membranes. Amino-acid sequence and the membrane environment are two critical factors controlling association, however, quantitative knowledge on such contributions is limited. In this work, we study the dimerization of helices in lipid bilayers using extensive parallel Monte Carlo simulations with recently developed algorithms [1]. The dimerization of Glycophorin A is examined employing a coarse-grain model that retains a level of amino-acid specificity, in three different phospholipid bilayers. Association is driven by a balance of protein-protein and lipid-induced interactions with the latter playing a major role at short separations. In all bilayers, sequence-specificity is evident by the formation of a clear interface between the helices that is modulated by the lipid environment. Extracted estimates of the dimerization affinity are in excellent agreement with experimental data [2]. Following a different approach, the effect of amino-acid sequence is studied using the four transmembrane domains of the epidermal growth factor receptor family in identical lipid environments. Detailed characterization of dimer formation and estimates of the free energy of association reveal that these helices present significant affinity to self-associate with certain dimers forming non-specific interfaces. We present results that support the role of lipid-mediated contributions to such effects with major implications on protein function. [1] Janosi L. and Doxastakis M., “Accelerating flat-histogram methods for potential of mean force calculations”, J. Chem. Phys., 131, 054105 (2009) [2] Janosi L., Prakash A. and Doxastakis M., “Lipid-Modulated sequence-specific association of Glycophorin A in membranes”, Biophys. J., 99, 284-292 (2010)

A GPU approach to parallel replica-exchange polymer simulations

We investigate new programming techniques for parallel tempering Monte Carlo simulations of an elementary bead-spring homopolymer model using graphics processing units (GPUs). For a precise estimation of statistical quantities, like the peak structure of the specific heat, a large number of conformations with substantial statistical data is needed. Therefore the advantage of gathering this data faster by improving the performance of Monte Carlo simulations cannot be overrated. With the huge computational capability of the large number of cores on GPUs, that can be exploited by means of multithreaded programming, we find significant increases in e_ciency compared to CPU-only simulations.

Availability of a general k-out-of-n:G system with non-identical components considering shut-off rules using quasi-birth–death process

A general repairable k-out-of-n:G system with non-identical components, which is a common form of redundancy, is considered in this paper. The number of repairmen is assumed to be r (1≤ r≤ n− k+1) while components can have similar or different repair priorities. The objective of this work is to address the problem of efficient evaluation of the system's availability in a way that steady state solutions can be obtained systematically in a reasonable computational time. This problem is modeled as a finite state-dependent non-homogeneous quasi-birth–death (QBD) process. An algorithm is introduced to systematically generate the system state vectors and transition rate matrix and then an iterative method based on the Block Gauss–Seidel method is employed to determine the steady state probabilities. These are the novel contributions made in this paper. An analog Monte Carlo simulation is presented to demonstrate the correctness and the efficiency of the proposed method.

Self-Association of Models of Transmembrane Domains of ErbB Receptors in a Lipid Bilayer

Association of transmembrane (TM) helices is facilitated by the close packing of small residues present along the amino-acid sequence. Extensive studies have established the role of such small residue motifs (GxxxG) in the dimerization of Glycophorin A (GpA) and helped to elucidate the association of TM domains in the epidermal growth factor family of receptors (ErbBs). Although membrane-mediated interactions are known to contribute under certain conditions to the dimerization of proteins, their effect is often considered nonspecific, and any potential dependence on protein sequence has not been thoroughly investigated. We recently reported that the association of GpA is significantly assisted by membrane-induced contributions as quantified in different lipid bilayers. Herein we extend our studies to explore the origin of these effects and quantify their magnitude using different amino-acid sequences in the same lipid environment. Using a coarse-grained model that accounts for amino-acid specificity, we perform extensive parallel Monte Carlo simulations of ErbB homodimerization in dipalmitoyl-phosphatidylcholine lipid bilayers. A detailed characterization of dimer formation and estimates of the free energy of association reveal that the TM domains show a significant affinity to self-associate in lipid bilayers, in qualitative agreement with experimental findings. The presence of GxxxG motifs enhances favorable protein-protein interactions at short separations. However, the lipid-induced attraction presents a more complex character than anticipated. Depending on the interfacial residues, lipid-entropic contributions support a decrease of separation or a parallel orientation to the membrane normal, with important implications for protein function.

Parallel Monte Carlo simulation for reliability and cost evaluation of equipment and systems

An algorithm for reliability simulation of equipment and systems using a parallel computing environment is developed. A sequential Monte Carlo simulation-based method is applied for generating reliability and maintainability history charts for equipment, from which the reliability indices, as well as their probability distributions, are calculated. A parallel methodology is developed to discretize the simulation into periods, and the simulation of individual periods is dispatched to different processors for efficient computing. Case studies utilizing the proposed method for reliability evaluation of equipment and simple systems are presented. Simulation results are compared with the analytical approach results for the same cases. The study provides the basis for using a parallel computing environment in power system reliability and cost evaluations. The environment used in this paper is Rock-131 from the San Diego Supercomputer Center in La Jolla, CA.

Structural transitions and dipole moment of water clusters (H2O)n=4–100

The properties of water clusters (H(2)O)(n) over a broad range of sizes (n=4-100) were studied by microcanonical parallel tempering Monte Carlo and replica exchange molecular dynamics simulations at temperatures between 20 and 300 K, with special emphasis in the understanding of relation between the structural transitions and dipole behavior. The effect of the water interaction potential was analyzed using six nonpolarizable models, but more extensive calculations were performed using the TIP4P-ice water model. We find that, in general, the dipole moment of the cluster increases significantly as the cluster melts, suggesting that it could be used to discriminate between the solidlike and liquidlike phases. The effect of a moderate electric field on the cluster heat capacity and total dipole moment was found to be negligible.

Lipid-Modulated Sequence-Specific Association of Glycophorin A in Membranes

Protein association in lipid membranes is a complex process with thermodynamics directed by a multitude of different factors. Amino-acid sequence is a molecular parameter that affects dimerization as shown by limited directed mutations along the transmembrane domains. Membrane-mediated interactions are also important although details of such contributions remain largely unclear. In this study, we probe directly the free energy of association of Glycophorin A by means of extensive parallel Monte Carlo simulations with recently developed methods and a model that accounts for sequence-specificity while representing lipid membranes faithfully. We find that lipid-induced interactions are significant both at short and intermediate separations. The ability of molecules to tilt in a specific hydrophobic environment extends their accessible interfaces, leading to intermittent contacts during protein recognition. The dimer with the lowest free energy is largely determined by the favorable lipid-induced attractive interactions at the closest distance. Finally, the coarse-grained model employed herein, together with the extensive sampling performed, provides estimates of the free energy of association that are in excellent agreement with existing data.

SU-GG-I-39: Evaluation of Different Benchmark Measurements for Validating Monte Carlo MDCT Source Models Used in Estimating Radiation Dose

Purpose:Monte Carlo simulation has become a popular method of estimating dose from multidetector computed tomography (MDCT) exams. Typically, MDCT modeling techniques and scanner‐specific input information are validated using physical dose measurements made under very simple conditions, such as homogenous cylindrical phantoms (i.e. CTDI). The purpose of this study is to investigate the sensitivity of such measurements in order to assess their utility in validating scanner‐specific x‐ray source information used in MDCT simulation packages. Method and Materials: Two 64‐slice MDCT scanners were used to obtain three types of benchmark measurements: 1) half and quarter value layers (HVL/QVL), 2) bowtie profiles (exposure values across the fan‐beam), and 3) CTDI100 values. Analogous Monte Carlo simulations were performed using two types of x‐ray source models; (A) one based on physical measurements and (B) one based on information provided from other sources. Results: HVL/QVL simulations using source model A resulted in average disagreements of 8.1% and 10.5% for Scanners 1 and 2, respectively, while source model B had errors of 2.2% and 39.8%. Similarly, bowtie profile simulations had average errors of 0.7% and 4.0% for source model A and 5.3% and 19.4% for source model B. Finally, CTDI100 simulation errors were 1.3% and 3.5% for source model A and 4.1% and 34.9% for source model B. Conclusion: The results show that the errors between measurements and simulations for the HVL/QVL and bowtie profile measurements reveal more fundamental and specific information about the accuracy of source models than does CTDI100 . This implies that CTDI100 may not be as sensitive as more elemental measurements for verifying the accuracy of more complicated dose simulations. This study shows that while CTDI100 may be a necessary benchmark it is not sufficient for MDCT Monte Carlodosimetry validation. Further testing in heterogeneous phantoms is needed to extend these sensitivity evaluations.

Abnormal grain growth of Goss grains in Fe–3% Si steel driven by sub-boundary-enhanced solid-state wetting: Analysis by Monte Carlo simulation

The selective abnormal grain growth (AGG) of Goss grains in Fe–3% Si steel was investigated using a parallel Monte Carlo (MC) simulation based on the concept of sub-boundary-enhanced solid-state wetting. Goss grains that contain low angle sub-boundaries will induce solid-state wetting against grain boundaries in the matrix that exhibit a moderate variation in grain boundary energy, as expected from the variation in boundary type. AGG generates a sharp Goss texture provided that only Goss-oriented grains have the required sub-grain structure to grow selectively by sub-boundary-induced wetting and that other orientations lack the required content of low angle boundaries. This behavior is shown in three-dimensional MC simulations of microstructure evolution with textures and grain boundary distributions matched to experimental data.