A scalable parallel algorithm for large-scale reactive force-field molecular dynamics simulations

A scalable parallel algorithm has been designed to perform multimillion-atom molecular dynamics (MD) simulations, in which first principles-based reactive force fields (ReaxFF) describe chemical reactions. Environment-dependent bond orders associated with atomic pairs and their derivatives are reused extensively with the aid of linked-list cells to minimize the computation associated with atomic n-tuple interactions (n?4 explicitly and ?6 due to chain-rule differentiation). These n-tuple computations are made modular, so that they can be reconfigured effectively with a multiple time-step integrator to further reduce the computation time. Atomic charges are updated dynamically with an electronegativity equalization method, by iteratively minimizing the electrostatic energy with the charge-neutrality constraint. The ReaxFF-MD simulation algorithm has been implemented on parallel computers based on a spatial decomposition scheme combined with distributed n-tuple data structures. The measured parallel efficiency of the parallel ReaxFF-MD algorithm is 0.998 on 131,072 IBM BlueGene/L processors for a 1.01 billion-atom RDX system.     

Pathfinder: A parallel search algorithm for concerted atomistic events

An algorithm has been designed to search for the escape paths with the lowest activation barriers when starting from a local minimum-energy configuration of a many-atom system. The pathfinder algorithm combines: (1) a steered eigenvector-following method that guides a constrained escape from the convex region and subsequently climbs to a transition state tangentially to the eigenvector corresponding to the lowest negative Hessian eigenvalue; (2) discrete abstraction of the atomic configuration to systematically enumerate concerted events as linear combinations of atomistic events; (3) evolutionary control of the population dynamics of low activation-barrier events; and (4) hybrid task + spatial decompositions to implement massive search for complex events on parallel computers. The program exhibits good scalability on parallel computers and has been used to study concerted bond-breaking events in the fracture of alumina.     

Embedded divide-and-conquer algorithm on hierarchical real-space grids: parallel molecular dynamics simulation based on linear-scaling density functional theory

A linear-scaling algorithm has been developed to perform large-scale molecular-dynamics (MD) simulations, in which interatomic forces are computed quantum mechanically in the framework of the density functional theory. A divide-and-conquer algorithm is used to compute the electronic structure, where non-additive contribution to the kinetic energy is included with an embedded cluster scheme. Electronic wave functions are represented on a real-space grid, which is augmented with coarse multigrids to accelerate the convergence of iterative solutions and adaptive fine grids around atoms to accurately calculate ionic pseudopotentials. Spatial decomposition is employed to implement the hierarchical-grid algorithm on massively parallel computers. A converged solution to the electronic-structure problem is obtained for a 32,768-atom amorphous CdSe system on 512 IBM POWER4 processors. The total energy is well conserved during MD simulations of liquid Rb, showing the applicability of this algorithm to first principles MD simulations. The parallel efficiency is 0.985 on 128 Intel Xeon processors for a 65,536-atom CdSe system.     

A general purpose parallel molecular dynamics simulation program

We present a general purpose parallel molecular dynamics simulation code. The code can handle NVE, NVT, and NPT ensemble molecular dynamics, Langevin dynamics, and dissipative particle dynamics. Long-range interactions are handled by using the smooth particle mesh Ewald method. The implicit solvent model using solvent-accessible surface area was also implemented. Benchmark results using molecular dynamics, Langevin dynamics, and dissipative particle dynamics are given.

Program summary

Title of program:MM_PARCatalogue identifier:ADXP_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADXP_v1_0Program obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandComputer for which the program is designed and others on which it has been tested:any UNIX machine. The code has been tested on Linux cluster and IBM p690Operating systems or monitors under which the program has been tested:Linux, AIXProgramming language used:CMemory required to execute with typical data:∼60 MB for a system of atoms Has the code been vectorized or parallelized? parallelized with MPI using atom decomposition and domain decompositionNo. of lines in distributed program, including test data, etc.:171 427No. of bytes in distributed program, including test data, etc.:4 558 773Distribution format:tar.gzExternal routines/libraries used:FFTW free software (http://www.fftw.org)Nature of physical problem:Structural, thermodynamic, and dynamical properties of fluids and solids from microscopic scales to mesoscopic scales.Method of solution:Molecular dynamics simulation in NVE, NVT, and NPT ensemble, Langevin dynamics simulation, dissipative particle dynamics simulation.Typical running time:Table below shows the typical run times for the four test programs.

Benchmark results. The values in the parenthesis are the number of processors used
System	Method	Timing for 100 steps in seconds
256 TIP3P	MD	23.8 (1)
64 DMPC + 1645 TIP3P	MD	890 (1)	528 (2)	326 (4)	209 (8)
8 Aβ16-22	LD	1.02 (1)
23760 Groot-Warren particles	DPD	22.16 (1)

Full-size table

     

A cell multipole based domain decomposition algorithm for molecular dynamics simulation of systems of arbitrary shape

A domain decomposition algorithm for molecular dynamics simulation of atomic and molecular systems with arbitrary shape and non-periodic boundary conditions is described. The molecular dynamics program uses cell multipole method for efficient calculation of long range electrostatic interactions and a multiple time step method to facilitate bigger time steps. The system is enclosed in a cube and the cube is divided into a hierarchy of cells. The deepest level cells are assigned to processors such that each processor has contiguous cells and static load balancing is achieved by redistributing the cells so that each processor has approximately same number of atoms. The resulting domains have irregular shape and may have more than 26 neighbors. Atoms constituting bond angles and torsion angles may straddle more than two processors. An efficient strategy is devised for initial assignment and subsequent reassignment of such multiple-atom potentials to processors. At each step, computation is overlapped with communication greatly reducing the effect of communication overhead on parallel performance. The algorithm is tested on a spherical cluster of water molecules, a hexasaccharide and an enzyme both solvated by a spherical cluster of water molecules. In each case a spherical boundary containing oxygen atoms with only repulsive interactions is used to prevent evaporation of water molecules. The algorithm shows excellent parallel efficiency even for small number of cells/atoms per processor.     

A general-purpose coarse-grained molecular dynamics program

In this article, we describe a general-purpose coarse-grained molecular dynamics program COGNAC (COarse Grained molecular dynamics program by NAgoya Cooperation). COGNAC has been developed for general molecular dynamics simulation, especially for coarse-grained polymer chain models. COGNAC can deal with general molecular models, in which each molecule consists of coarse-grained atomic units connected by chemical bonds. The chemical bonds are specified by bonding potentials for the stretching, bending and twisting of the bonds, each of which are the functions of the position coordinates of the two, three and four atomic units. COGNAC can deal with both isotropic and anisotropic interactions between the non-bonded atomic units. As an example, the Gay-Berne potential is implemented. New potential functions can be added to the list of existing potential functions by users. COGNAC can do simulations for various situations such as under constant temperature, under constant pressure, under shear and elongational deformation, etc. Some new methods are implemented in COGNAC for modeling multiphase structures of polymer blends and block copolymers. A density biased Monte Carlo method and a density biased potential method can generate equilibrium chain configurations from the results of the self-consistent field calculations. Staggered reflective boundary conditions can generate interfacial structures with smaller system size compared with those of periodic boundary conditions.     

Partitioning 3D space for parallel many-particle simulations

In a common approach for parallel processing applied to simulations of many-particle systems with short-ranged interactions and uniform density, the cubic simulation box is partitioned into domains of equal shape and size, each of which is assigned to one processor. We compare the commonly used simple-cubic (SC) domain shape to domain shapes chosen as the Voronoi cells of BCC, FCC, and HCP sphere packings. The latter three are found to result in superior partitionings with respect to communication overhead. Scaling of the domain shape is used to extend the range of applicability of these partitionings to a large set of processor numbers. The higher efficiency with BCC and FCC partitionings is demonstrated in simulations of the sillium model for amorphous silicon.     

Optimization techniques for parallel force-decomposition algorithm in molecular dynamic simulations

The efficiency and scalability of traditional parallel force-decomposition (FD) algorithms are not good because of high communication cost which is introduced when skew-symmetric character of force matrix is applied. This paper proposed a new parallel algorithm called UTFBD (Under Triangle Force Block Decomposition), which is based on a new efficient force matrix decomposition strategy. This strategy decomposes only the under triangle force matrix and greatly reduces parallel communication cost, e.g., the communication cost of UTFBD algorithm is only one third of Taylor's FD algorithm. UTFBD algorithm is implemented on Cluster system and applied to solve a physical nucleation problem with 500,000 particles. Numerical results are analyzed and compared among three algorithms, namely, FRI, Taylor's FD and UTFBD. The efficiency of UTFBD on 105 processors is 41.3%, and the efficiencies of FRI and Taylor's FD on 100 processors are 4.3 and 35.2%, respectively. In another words, the efficiency of UTFBD on about 100 processors is 37.0 and 6.1% higher than that of FRI and Taylor's FD, respectively. Results show that UTFBD can increase the efficiency of parallel MD (Molecular Dynamics) simulation to a higher degree and has a better scalability.     

Multibillion-atom molecular dynamics simulation: Design considerations for vector-parallel processing

Progress in adapting molecular dynamics algorithms for systems with short-range interactions to utilize the features of modern supercomputers is described. Efficient utilization of the latest generation of processor architectures requires algorithms that can be both vectorized and parallelized. The approach adopted for vectorization involves combining the layer and neighbor-list methods, while parallelization employs spatial subdivision with explicit communication. The techniques presented here have been used in performance tests on the Cray X1 vector-parallel supercomputer with systems containing over 12 billion atoms.     

A finite-difference eigenvalue algorithm for calculating the band structure of a photonic crystal

A new method based on a finite difference of the governing differential equation for the eigenvalue problem is introduced to calculate the band structure of a two-dimensional photonic crystal. The effective medium technique is also used in the method. The problem is reduced to a standard matrix eigenvalue problem. Compared to the conventional plane wave expansion method, the present method improves the convergence of the solution and thus is a fast and accurate algorithm for calculating the band structure of a photonic crystal.     

An efficient parallel implementation of the smooth particle mesh Ewald method for molecular dynamics simulations

This paper focuses on the implementation and the performance analysis of a smooth particle mesh Ewald method on several parallel computers. We present the details of the algorithms and our implementation that are used to optimize parallel efficiency on such parallel computers.     

A fine grained parallel smooth particle mesh Ewald algorithm for biophysical simulation studies: Application to the 6-D torus QCDOC supercomputer

In order to model complex heterogeneous biophysical macrostructures with non-trivial charge distributions such as globular proteins in water, it is important to evaluate the long range forces present in these systems accurately and efficiently. The Smooth Particle Mesh Ewald summation technique (SPME) is commonly used to determine the long range part of electrostatic energy in large scale molecular simulations. While the SPME technique does not give rise to a performance bottleneck on a single processor, current implementations of SPME on massively parallel, supercomputers become problematic at large processor numbers, limiting the time and length scales that can be reached. Here, a synergistic investigation involving method improvement, parallel programming and novel architectures is employed to address this difficulty. A relatively simple modification of the SPME technique is described which gives rise to both improved accuracy and efficiency on both massively parallel and scalar computing platforms. Our fine grained parallel implementation of the modified SPME method for the novel QCDOC supercomputer with its 6D-torus architecture is then given. Numerical tests of algorithm performance on up to 1024 processors of the QCDOC machine at BNL are presented for two systems of interest, a β-hairpin solvated in explicit water, a system which consists of 1142 water molecules and a 20 residue protein for a total of 3579 atoms, and the HIV-1 protease solvated in explicit water, a system which consists of 9331 water molecules and a 198 residue protein for a total of 29508 atoms.     

A parallel implementation of the Wang-Landau algorithm

The Wang-Landau algorithm is a flat-histogram Monte Carlo method that performs random walks in the configuration space of a system to obtain a close estimation of the density of states iteratively. It has been applied successfully to many research fields. In this paper, we propose a parallel implementation of the Wang-Landau algorithm on computers of shared memory architectures by utilizing the OpenMP API for distributed computing. This implementation is applied to Ising model systems with promising speedups. We also examine the effects on the running speed when using different strategies in accessing the shared memory space during the updating procedure. The allowance of data race is recommended in consideration of the simulation efficiency. Such treatment does not affect the accuracy of the final density of states obtained.     

A hybrid multi-loop genetic-algorithm/simplex/spatial-grid method for locating the optimum orientation of an adsorbed protein on a solid surface

Atomistic simulation of protein adsorption on a solid surface in aqueous environment is computationally demanding, therefore the determination of preferred protein orientations on the solid surface usually serves as an initial step in simulation studies. We have developed a hybrid multi-loop genetic-algorithm/simplex/spatial-grid method to search for low adsorption-energy orientations of a protein molecule on a solid surface. In this method, the surface and the protein molecule are treated as rigid bodies, whereas the bulk fluid is represented by spatial grids. For each grid point, an effective interaction region in the surface is defined by a cutoff distance, and the possible interaction energy between an atom at the grid point and the surface is calculated and recorded in a database. In searching for the optimum position and orientation, the protein molecule is translated and rotated as a rigid body with the configuration obtained from a previous Molecular Dynamic simulation. The orientation-dependent protein-surface interaction energy is obtained using the generated database of grid energies. The hybrid search procedure consists of two interlinked loops. In the first loop A, a genetic algorithm (GA) is applied to identify promising regions for the global energy minimum and a local optimizer with the derivative-free Nelder-Mead simplex method is used to search for the lowest-energy orientation within the identified regions. In the second loop B, a new population for GA is generated and competitive solution from loop A is improved. Switching between the two loops is adaptively controlled by the use of similarity analysis. We test the method for lysozyme adsorption on a hydrophobic hydrogen-terminated silicon (110) surface in implicit water (i.e., a continuum distance-dependent dielectric constant). The results show that the hybrid search method has faster convergence and better solution accuracy compared with the conventional genetic algorithm.     

Oedometric test, Bauer's law and the micro-macro connection for a dry sand

What is the relationship between the macroscopic parameters of the constitutive equation for a granular soil and the microscopic forces between grains? In order to investigate this connection, we have simulated by molecular dynamics the oedometric compression of a granular soil (a dry and bad-graded sand) and computed the hypoplastic parameters hs (the granular skeleton hardness) and η (the exponent in the compression law) by following the same procedure than in experiments, that is by fitting the Bauer's law e/e₀=exp(−n(3p/hs)), where p is the pressure and e₀ and e are the initial and present void ratios. The micro-mechanical simulation includes elastic and dissipative normal forces plus slip, rolling and static friction between grains. By this way we have explored how the macroscopic parameters change by modifying the grains stiffness, V; the dissipation coefficient, γn; the static friction coefficient, μs; and the dynamic friction coefficient, μk. Cumulating all simulations, we obtained an unexpected result: the two macroscopic parameters seems to be related by a power law, hs=0.068(4)η9.88(3). Moreover, the experimental result for a Guamo sand with the same granulometry fits perfectly into this power law. Is this relation real? What is the final ground of the Bauer's Law? We conclude by exploring some hypothesis.     

A domain decomposition molecular dynamics program for the simulation of flexible molecules of spherically-symmetrical and nonspherical sites. II. Extension to NVT and NPT ensembles

We describe the algorithms for NVT and NPT-ensemble simulations developed within the parallel molecular dynamics program GBMOLDD. This program uses the domain decomposition algorithm and is targeted at large-scale simulations of molecular systems (particularly polymers and liquid crystals) composed of both spherically-symmetric and nonspherical sites. The nonspherical sites can be described either by a Gay-Berne potential or by soft repulsive spherocylinders. The molecules can be of arbitrary topology and the intramolecular forces are described via standard force fields. We tested the stability of both leap-frog and velocity-Verlet integrators on two “real-life” systems—a nematic liquid crystal phase of 1944 one-site Gay-Berne molecules and on 512 flexible liquid-crystalline dimers. In both cases the algorithm demonstrates good stability over the typical simulation times required for new phase formation and/or molecular relaxation processes.     

On the density correlation of the spontaneous fluctuation in a real-coded lattice gas

We investigate the spontaneous fluctuation in a real-coded lattice gas (RLG) model by studying the density correlation function. In particular, the dynamic structure factor obtained from RLG is in agreement with the Rayleigh-Brillouin spectrum, a spectrum which can also be measured from a real fluid at the continuum limit. We also work out the analytic form of the static structure factor for the RLG model, which is supported by the numerical results.     

Massively parallel quantum computer simulator

We describe portable software to simulate universal quantum computers on massive parallel computers. We illustrate the use of the simulation software by running various quantum algorithms on different computer architectures, such as a IBM BlueGene/L, a IBM Regatta p690+, a Hitachi SR11000/J1, a Cray X1E, a SGI Altix 3700 and clusters of PCs running Windows XP. We study the performance of the software by simulating quantum computers containing up to 36 qubits, using up to 4096 processors and up to 1 TB of memory. Our results demonstrate that the simulator exhibits nearly ideal scaling as a function of the number of processors and suggest that the simulation software described in this paper may also serve as benchmark for testing high-end parallel computers.     

A recursive algorithm for finding HDMR terms for sensitivity analysis

High Dimensional Model Representation (HDMR) is an efficient technique which decomposes a multivariate function into a constant, univariate, bivariate functions and so on. These functions are forced to be mutually orthogonal by means of an orthogonality condition. The technique which is generally used for high-dimensional input-output systems can be applied to various disciplines including sensitivity analysis, differential equations, inversion of data and so on. In this article we present a computer program that computes individual components of HDMR resolution of a given multivariate function. The program also calculates the global sensitivity indices. Lastly the results of the numerical experiments for different set of functions are introduced.