Molecular Dynamics on Parallel Computers

The following sections are included:

• PREFACE

• CONTENTS

The investigation of material properties of fullerenes and especially nanotubes in experiments is very difficult. Here, computer simulations are an important tool to gain further insight. To obtain realistic results we use the bond-order potential due to Brenner. This potential describes the bonding structure (the forming/breaking of bonds) and other properties of graphite and hydrocarbons quite well. Based on this potential, we implemented a parallel algorithm for nanotube simulation. It uses the well known linked cell technique for short range pair potentials. However, due to the sophisticated structure of Brenner's many-body potential, a special implementation is necessary to gain an O(N) complexity of the overall computation. For its parallelization, the domain decomposition method can be used in a straightforward way. We present the results of our numerical experiments. We considered the reactive collision of C₆₀ fullerenes with benzyne and, moreover, the bending and stretching of nanotubes with nearly 100 million atoms. To our knowledge, this is the largest nanotube simulation up to now. The parallel algorithm runs on a workstation cluster Parnass2 and a CRAY T3E (512 proc.) It shows a nearly linear speed up and exhibits a very good scaling behavior. It turns out that nanotubes, at least in our simulations, can sustain hundreds of Giga Pascal in stretching experiments before fraction occurs.

The function of many important biomolecules comes from their dynamic properties and their ability to switch between different conformations. In a conformation, the large scale geometric structure of the molecule is understood to be conserved, whereas on smaller scales the system may well rotate, oscillate or fluctuate. Most approaches to the identification of conformations neglect the dynamical aspect: they are interested only in finding clusters of molecular configurations with significantly different large scale geometric structure. Unlike these approaches, we herein advocate to directly attack the determination of conformations together with the computation of their stability time spans and the rate of transitions between them. Recently, a novel algorithmic concept has been proposed ¹, which allows to decompose the position space into such dynamically defined conformational subsets and to compute the corresponding transition probabilities. The article reviews the basic steps of this algorithm with special emphasis on the generation of appropriate sampling and transition data via a specific reweighted Hybrid Monte Carlo variant and presents the performance of the algorithm when applied to a triribonucleotide.

A new modular program system, DMMD (Distributed Memory Molecular Dynamics) has been designed which employs modern elements of the Fortran 90 language for abstraction and encapsulation. The program uses the message passing paradigm which is applicable to both distributed and shared memory systems. Communication is handled in a special interface library which efficiently translates high level Fortran 90 concepts into the more primitive addressing mechanisms of current message passing protocols. Furthermore it allows transparent switching between message passing libraries. This makes it a powerful, flexible tool for parallel Fortran 90 applications. The modular concept of DMMD allows easy modification and extensions of its capabilities towards new interaction models or integrators. The design has been tested in a serial version which also serves as a reference point for the parallel version which will be used to benchmark several strategies for data and work distribution. We report first results obtained with DMMD on a CRAY T3E system for a simple Lennard-Jones mixture and comparisons with other massively parallel molecular dynamics codes.

NWChem is the massively parallel computational chemistry software package developed by the High Performance Computational Chemistry group for the Environmental Molecular Sciences Laboratory at the U.S. Department of Energy's Pacific Northwest National Laboratory. This software includes modules for quantum-mechanical and classical computational chemistry. This contribution details the implementation and performance of the molecular dynamics simulation module. The module is based on a domain decomposition approach to make it efficient for the simulation of large molecular systems on massively parallel computer architectures, despite the added complexity to deal with the necessary periodic atomic reassignments and dynamic load balancing techniques for heterogeneous systems. The calculation of atomic forces is the component that requires most of the communication and computational resources in molecular simulations. The use of the Global Array tools for one-sided asynchronous communication between processors allows this component to be carried out very efficiently when dynamic load balancing techniques are used.

The intention of this article is two-fold: First, we want to describe IMD (ITAP Molecular Dynamics), a software package for classical molecular dynamics simulations on massively-parallel computers. The software is a general purpose package which can be used for any kind of condensed matter, but it contains a number of special features for the simulation of the mechanical properties of solids, especially of quasicrystals and silicon interfaces. Second we want to describe a representative application of IMD to the simulation of shock waves in binary icosahedral quasicrystals. These simulations have been carried out to demonstrate the feasibility of metacomputing (clustering of massively-parallel supercomputers over long distances).

Some first attempts to systematically link different time and length scales in computer simulations of dense polymer systems are discussed. The simulations cover the range from atomistic (force field) simulations all the way up to the semimacro-scopic regime, where the chains are represented by soft overlapping ellipsoidal particles.

Configurational Bias Monte Carlo (CBMC) simulations are often used to compute thermo-dynamic properties of flexible chain molecules. In particular for long chains, these calculations are still computationally expensive. In this article we will present several improved algorithms for CBMC and illustrate how this method can be implemented efficiently on a parallel computer.

The shape and the dynamics of tethered polymers in flow are described by bead–spring models taking into account hydrodynamic interactions and excluded volume effects. The Brownian dynamics of these bead–spring models is simulated by replacing the discretized Langevin equation with a scheme which introduces artificial inertia for the beads. With this scheme the integration time step can be chosen up to a factor of 10 larger for a prescribed error in the equilibrium distribution than in conventional Brownian dynamics simulations. Besides various applications of this scheme, we devise an efficient way to calculate the relaxation spectrum and -modes from the simulation data using the Karhunen–Loève method.

In this paper we discuss the particle mesh Ewald (PME) method and its application to the simulation of biological membranes. Since the PME method allows a fast calculation of the electrostatic interactions a simulation of a membrane patch for 10ns is performed. The length of the simulation allows the calculation of dynamical properties such as the translational diffusion coefficients. The simulation allows furthermore, the characterization of the motions of individual groups within the lipid molecules. The time scale of the movement of individual groups as well as the time scale of other physical quantities are discussed.

Molecular dynamics (MD) simulations of proteins provide descriptions of atomic motions, which allow relating observable properties of proteins to microscopic processes. In particular, the 'mechanics' at the atomic level of these 'nano-machines' can be probed by recently developed efficient MD techniques such as single molecule unbinding simulations or 'conformational flooding'. MD simulations in general are computationally very expensive and, therefore, are limited to time scales of nanoseconds, even if the latest supercomputers are used. In contrast, the elementary steps of protein function are typically in the micro- to millisecond range. Moreover, conventional MD simulations are based on force fields and thus, generally, cannot describe chemical reactions. This means, for example, that enzymatic catalysis is beyond the scope of conventional MD methods. In this contribution we describe a concerted effort to overcome these limitations. A fast multiple time step structure adapted multipole method (FAMUSAMM) has been developed to speed up the evaluation of the computationally most demanding Coulomb interactions in solvated protein models, thereby providing the basis for efficient MD simulations of any kind. On top of that, and based on a description of protein dynamics within the framework of non-equilibrium statistical mechanics, a method to predict the reaction pathway and target structure of slow conformational transitions ('conformational flooding') is described. To enable treatment of localized chemical reactions, e.g., at the active sites of enzymes, and utilizing the efficient FAMUSAMM description, an interface to the CPMD density functional package by M. Parrinello et al. (J. Hutter et al., MPI für Festkörperforschung, Stuttgart, 1995-1998) is introduced, which allows hybrid quantum mechanical / molecular dynamics (QM/MD) simulations. All these methods have been integrated within the parallel MD program EGO.

In this paper we give a short overview of our investigations on a model glass forming liquid, consisting of a binary Lennard-Jones system. The measurement of the two times autocorrelation function below the glass transition temperature T_g shows a strong dependence of the time t_w elapsed since the preparing of the system. This phenomenon is called ageing. We characterize it by means of scaling and try to find its origin.

We developed a parallelized ab initio molecular dynamics code for the investigation of materials properties and processes based on the local density approximation to the density functional theory, using separable, norm-conserving pseudo-potentials and a plane-wave representation of the wave functions. We briefly describe the methodological formulation of the many-body potential of interacting atoms generated by a quantum mechanical treatment of the valence electrons and discuss the parallelization strategies and speed-ups. Two different parallelization strategies are discussed. It is shown, that distributing the wave-function representations in real and reciprocal space to computing nodes leads to a flexible approach where many computing nodes can be efficiently used in parallel. We found that the degree of parallelization of the program is 99.5% which permits very efficient calculations on 128 computing nodes for unit cells containing up to 400 atoms. As a physical application we present results of calculations of reaction rates of ad-atoms on As-terminated Si(111).

The MD simulation program TABB is motivated by the need of long time simulations for the investigation of slow processes near the glass transition of glass forming alloys. TABB is written in C++ with a high degree of flexibility: TABB allows the use of any short ranged pair potentials or EAM potentials, by generating and using a spline representation of all functions and their derivatives. TABB supports several numerical integration algorithms like the Runge-Kotta or the modified Gear-predictor-corrector algorithm of order five. The boundary conditions can be chosen to resemble the geometry of bulk materials or films. The simulation box length or the pressure can be fixed for each dimension separately. TABB may be used in isokinetic, isoenergeric or canonic (with random forces) mode. TABB contains a simple instruction interpreter to easily control the parameters and options during the simulation. The same source code can be compiled either for workstations or for parallel computers. The main optimization goal of TABB is to allow long time simulations of medium or small sized systems. To make this possible, much attention is spent on the optimized communication between the nodes. TABB uses a domain decomposition procedure. To use many nodes with a small system, the domain size has to be small compared to the range of particle interactions. In the limit of many nodes for only few atoms, the bottle neck of communication is the latency time. TABB minimizes the number of pairs of domains containing atoms that interact between these domains. This procedure minimizes the need of communication calls between pairs of nodes. TABB decides automatically, to how many, and to which directions the decomposition shall be applied. E.g., in the case of one dimensional domain decomposition, the simulation box is only split into "slabs" along a selected direction. The three dimensional domain decomposition is best with respect to the number of interacting domains only for simulations with medium range interactions on machines with many nodes (112 and more). TABB uses a modified cell algorithm to update a neighbor list for the force evaluation where the cells are small compared to the interaction range and contain the atoms in a continuous array. The updates of the neighbor list and of the domain decomposition are carried out synchronously with a flexible frequency. A small Ni_0.5Zr_0.5 system near the glass transition can be computed over 1 μs using 2.5 fs time steps within less than one week on a T3E-900 machine. Some examples of the computational results will be presented ¹.

Molecular Dynamics simulations in the Tight-Binding approach allow the study of the ionic and electronic structures of semiconductors. The Tight-Binding codes are characterized by inhomogeneous data distribution and require the repeated diagonalization of a large sparse matrix to compute the whole body of its eigenvalues and eigenvectors. The code parallelization, by using the High Performance Fortran (HPF) environment, and the integration of optimized parallel mathematical routines is described.

For the last ten years there has been an enormous progress in the simulation of granular media like sand or powders. These simulations consist in simulating trajectories of each particle individually. Essentially one has to solve the Newton's equations including the effects of Coulomb friction and the physics occuring at a collision. But the details of the trajectories are not important for the collective behaviour. Therefore simplifications are introduced on the smallest scales.

I will introduce various methods like molecular dynamics that are used to simulate large amounts of particles (over 10⁹). Some of these medhods are based on the exploitation of parallelisation and metacomputing. Other approaches are more stochastic (DSMC Direct Simulation Monte Carlo) which simplify the calculation of collisions, positions and collision times. Very successful has been also the use of cellular automata which have been able to predict details such as the logarithmic tale of sand heaps. I will also discuss numerical techniques used for the surrounding fluid. This can be water in the case of sedimentation or air when one studies the formation of dunes in the desert. The calculation of velocity and pressure field of the fluid are done using multigrid techniques on parallel computers.

We will compare the performance of the various techniques and show some benchmarks on the dependence on the size of the system, the density of particles and the number of processors used.

When dry granular media are processed the grains tend to get charged. We show what effect electrical charge has on collisional cooling, a basic phenomenon of granular media. The free cooling of monopolar charged granular systems is studied by molecular dynamics simulations. Comparing our computer results for the dissipation rate to the known dissipation rate for uncharged granular systems, we find that the electrical charges reduce the dissipation rate by a Boltzmann factor because of the Coulomb repulsion between the particles.

Parallel threaded algorithms to speed up discrete–element simulations for granular matter are presented in this paper. Threads are a powerful tool for parallelization on shared-memory machines. First, a compact introduction to programming with threads will be given. Then an algorithm will be presented, which allows to determine neighborhood relations in polydisperse mixtures of particles of arbitrary shape, either discs, ellipses, or polygons. Another algorithm described here calculates the distance of two polygons in constant time, independent of the complexity of the shape of the polygons. This makes fast simulations of polygonal assemblies possible.

We investigated numerically the dust growth driven by Brownian motion in a proto-planetary disc around a solar-type young stellar object. This process is considered as the first stage in the transformation of the initially micron-sized solid particles to a planetary system. In contrast to earlier studies the growth was investigated at the small particle number densities typical for the conditions in a proto-planetary disc. Under such circumstances, the mean particle distance exceeds the typical aggregate diameter by orders of magnitude, and a collision will be a very rare event. We derived a criterion which allows an efficient detection of candidates for imminent collisions. The N-particle-method we used is based upon an adaptive time step scheme respecting the individual dynamical states of the aggregates. Its basic concept is to perform on average constant "length steps", instead of using constant time steps. The numerical cost of the algorithm scales with the particle number better than N log N. In order to minimise the influence of the decreasing number of particles within the simulation box, a new rescaling method is used throughout the aggregation process. Our numerical results indicate that at very low number densities, the growth process is influenced by spatial number density fluctuations.

Astrophysical challenges which demand the solution of the one million (or more) gravitating body problem are briefly discussed for the fields of cosmology, galactic nuclei and globular star clusters. Results from the classical three-body problem are to be combined with methods of statistical mechanics or thermodynamics in order to provide a reasonable model for astrophysically relevant particle numbers of the order 10⁴ – 10⁸. For high-precision modelling of star clusters and galactic nuclei over the required hundreds and thousands of crossing times specially tailored "brute-force" N-body integration methods are used together with special purpose hardware. A recent implementation of such a code on a general purpose parallel computer is presented. After a comparison of the relevance of different N-body algorithms a new concept for a more flexible customized special purpose computer based on a combination of GRAPE and FPGA is proposed. It is an ideal machine for all kinds of N-body simulations using neighbour schemes, as the Ahmad-Cohen direct N-body codes and smoothed particle hydrodynamics (SPH) for systems including interstellar gas dynamics.

We present a brief introduction to Smoothed Particle Hydrodynamics (SPH), where we focus on an approach to treat the physical viscosity. Additionally we describe in detail the basic principles of our parallel implementation of the SPH method. In the second part we summarize some results of test simulations performed to achieve general features of the method, and we end discussing some astrophysical applications, namely the simulation of accretion disks in Cataclysmic Variables.

During the last years, immense progress has been achieved in the stochastic treatment of quantum chromodynamics (QCD) on the lattice. Algorithms based on a molecular dynamics like treatment of gauge and quark fields have opened the door to exact ab initio simulations. The Hybrid Monte Carlo (HMC) algorithm, the most prominent representative of full QCD algorithms, is the favorite scheme to include virtual dynamical fermion creation and annihilation processes, which are the main obstacle on the way towards realistic QCD simulations. In this talk, I want to discuss methodical and practical aspects of the HMC for full QCD simulations, with emphasis on molecular dynamics. I will touch basics, merits and shortcomings of the HMC, recent algorithmic improvements and current state-of-the-art simulations, and I will dare to give a prognosis for the costs of future full QCD simulations deep in the chiral regime of vanishing quark masses.

The following sections are included:

• Introduction

• Multi-bosonic actions

• Two-step multi-bosonic algorithm

• Least-squares optimized polynomials

• Outlook

• References

References

The ability to determine the free energy of solvation for a number of small organic molecules of varying size and properties from the coordinate trajectory of a single simulation of a given reference state is investigated. The relative free energies were estimated from a single step perturbation using the perturbation formula. The reference state consisted of a cavity surrounded by solvent. To enhance sampling a soft-core interaction was used for the cavity. The effect of the size of the cavity, the effective barrier height and the length of simulation on the ability to reproduce results obtained from thermodynamic integration calculations is considered. It is shown that results using a single step perturbation from an appropriately chosen initial state are comparable to results from thermodynamic integration calculations for a wide range of compounds. Using a large number of compounds the computational efficiency is potentially increased by 2 to 3 orders of magnitude over traditional free energy approaches. Factors determining the efficiency of the approach are discussed.

The following sections are included:

• Acknowledgments

• References

Brownian dynamics simulations of a tethered polymer in uniform flow show that polymers are neither impenetrable nor free-draining. According to these results a f-shell blob model has been introduced which contains the well known non-draining blob model as well as the free-draining blob model as limiting cases. In addition, this model describes the transition from a nearly impenetrable polymer coil at small flow velocities to a nearly free-draining uncoiled polymer at large flow velocities.

The growth dynamics of convection patterns in incompressible one-component and binary fluid layers heated from below are discussed. The investigated system is two-dimensional with periodic lateral and no-slip impermeable horizontal boundary conditions. Starting from different initial states the transients are investigated. In pure fluids the further temporal dynamics is independent of the structure of the initial state and depends only on one critical mode. In binary mixtures, however, the whole transient dynamics is strongly influenced by the initial condition due to two critical modes that finally interact nonlinearly. The behavior close to and far from the boundaries is also compared.

The self-diffusion coefficient tensor of chiral liquid crystals was investigated by molecular dynamics (MD) simulations in the NVE ensemble. The bulk systems were composed of uniaxial chiral molecules with interactions described by the Gay-Berne potential and an additive term for the energy of chiral interaction. A cholesteric phase (N*) and a blue phase (BP II) were observed along an isotherm in dependence on the chirality parameter.

The following sections are included:

• Acknowledgments

• References

Molecular dynamics simulation of melittin at high dilution in close proximity to a fully hydrated l-palmitol-2-oleoyl-3-phosphatidylcholine (POPC) bilayer membrane has been performed at the atomic resolution. Our simulation shows that a single melittin molecule could not insert into the bilayer under ambient conditions. Furthermore, melittin exhibits various orientations at the bilayer/water interface.

A pore formed by four melittin helices has been studied by the atomic level computer simulation. A combined Monte Carlo(MC)/Molecular Dynamics(MD) technique was employed in this simulation to avoid the spurious effects caused by unequilibrated water and lipids, and thereby we can examine the stability of the pore. The pore was found to expand from its initial closed-packed structure, and the lipids and water molecules were found to fill the gaps between the helices. The permeation of water molecules through such an expanding pore has also been investigated.

Using Larson's model and its later generalization we have developed a microscopic model of gemini surfactants in aqueous medium. Gemini surfactants are two-tailed surfactants whose heads are connected by a hydrocarbon spacer chain; the spacer chain can be either hydrophobicor hydrophilic. By carrying out Monte Carlo (MC) simulations, we have studied the variation of the Critical Micellar Concentration (CMC) with the variation of a) the length of the spacer chain, b) the length of the hydrophobic tail and c) the bending rigidity of the chain forming the spacer and the tail. These results are in agreement with the available experimental results. We have also looked at the morphology of the micellar aggregates formed by the gemini surfactants. In contrast to the spherical micelles formed by the single tail surfactants, gemini surfactants with short spacer form thread-like micelles and those with long spacer (in particular when the spacer length is comparable to the tail length) form rod-like micelles. In an ongoing investigation, this approach is being used for a mixture of gemini surfactants and single tail surfactants. Preliminary results indicate that this system might exhibit micelle crosslinking as well as micelle swelling, depending on the parameters.

To predict how proteins or peptides fold based on sequence information alone is one of the grand challenges in computational biology. In principle, molecular dynamics simulation techniques may be used to simulate the process of folding directly. This is true given an appropriate description of the inter-atomic interactions and an ability to simulate on an appropriate time scale …

We use molecular dynamics simulations to study polyelectrolytes end-grafted to a surface. Our model includes counterions explicitly and treats the full long-range Coulomb interaction. For strongly charged brushes with no added salt, the counterions are confined inside the brush layer. We find that the brush thickness is linearly proportional to the chain length and the grafting density. The counterion distribution is strongly inhomogeneous, and counterion condensation can be observed although the Bjerrum length is smaller than the average bond length. The osmotic counterion pressure is much higher than the ideal gas value. Counterion diffusion is anisotropic and is enhanced at higher grafting densities.

We have performed ab initio total energy calculations and molecular dynamics simulations of dye molecules on NaCl(100) surfaces and in water. The flat dye molecule trimethine, [C₁₉H₁₇N₂O₂]⁺, which forms sandwich-like structures, and the cya-nine molecule monomethine, [C₂₁H₂₃N₂]⁺, which shows a typical stereochemical deformation due to two repulsive methyl groups, have been simulated. The molecular dynamics simulations are able to reproduce the experimentally observed configuration of the charged molecules. The behavior of trimethine, [C₇H₁₅N₂]⁺, in water has also been simulated.

The atomic displacements of molten Cu₃₃Zr₆₇ are analyzed by molecular dynamics to temperatures below the glass transition temperature T_g. We observe hopping processes and an underlying flow. With decreasing temperature the jumps become increasingly reversible and their contribution to diffusion diminishes.

Calculations were made on different platforms. On parallel computers atom decomposition is used and a high degree of parallelization is achieved.

References

We study nuclear quantum effects on the structural properties of hydrogen at 50 K and beyond 150 GPa by a method that takes into account the quantum nature of both electrons and protons ¹. The structure of solid hydrogen at such high compression is found to be unexpectedly diffuse, sharing characteristics with the high–T fluid or plasma–like phases. This suggests the concept of a "quantum fluxional solid" which cannot be understood in terms of a single underlying classical structure.

We study hydrogen at very high densities using semi-classical "Wave-Packet-Molecular-Dynamics" simulations, in which the electrons are represented by Gaussian wave-packets described by parameters which follow a pseudo-Hamiltonian dynamics. Exchange effects are considered by a pairwise anti-symmetrization of the many-body Hartree wave-function. This method allows to simulate hydrogen at densities of the order of the solid-state density and at various temperatures, covering the wide range from a molecular solid state, over a liquid up to a fully ionized plasma. From these simulations we are able to extract pressure, pair correlation functions, autocorrelation functions of current and velocity as well as the resulting transport coefficients. Our results are in quite good agreement with both other simulation methods (e.g. PIMC), and analytical approaches.

An ab initio molecular dynamics procedure implemented for parallel processing is presented. The method is based on gradient corrected density functionals for exchange and correlation employing Gaussian atomic basis. Sufficiently long trajectories can be calculated at an acceptable computational cost, allowing for analysis of dynamical behaviour of small metallic clusters. This is illustrated on examples of Li₈ and clusters. It can be shown that different isomers, even those close in energy, can exhibit distinct dynamics, dependent on their particular type of structure.

Observations of the Milky Way and nearby galaxies show that stars are grouped in clusters of up to millions of stars ('globular' clusters). As their velocity dispersion measures up to their integrated mass, they are largely self-gravitating Newtonian systems. Until recently the sphericity of old clusters in our galaxy has motivated theorists to understand star clusters as non-rotating, spherical entities.⁶

However, a careful study of 100 Milky Way clusters, and 173 clusters in M31, has shown small but significant departures from spherical symmetry: fits to their projected isophotes yielded ellipticities < ε > ≡ 1- < a/b >≈ 0.08 ± 0.01, when <ε> = 0 is expected.⁷ Furthermore, young clusters in the Large Magellanic Cloud show ellipticities as large as ε = 0.3. This raises the possibility that clusters are formed as rotating structures, which then evolve towards rounder configurations by releasing angular momentum.¹ Gravitational cold collapse simulations of non-rotating but aspherical galaxies already hint to long-lived rotation profiles in young self-gravitating systems (CMB, unpublished). Our objective is to quantify this process by developing computer models of star clusters with up to 50,000 member stars, using star-by-star parallel integration methods on CRAY computers.

We present our approach to develop a more general and realistic model of spherical and non-rotating globular clusters using an anisotropic gaseous model with unequal stellar masses (mass spectrum). We describe the numerical method and focus briefly on the advantages of a parallel version of our code.

The implementation of a hybrid "self consistent field" (SCF)¹ and direct Aarseth N–body integrator (NBODY6)² which synthesises the advantages of the direct force calculation with the efficiency of the field method is described. The resulting code is aimed for use on parallel architectures and is therefore applicable for collisional N–body integrations with extraordinarily large particle numbers. It opens the perspective to simulate the dynamics of globular clusters with realistic relaxation as well as stellar systems surrounding a supermassive black hole in galactic nuclei.

A technique for modelling collisional plasmas is demonstrated which uses a hierarchical tree-code – an accelerated molecular dynamics algorithm with an N log N computation time – adapted to model periodic, non-equilibrium two-component plasmas.

The following sections are included:

• Introduction

• Smoothed particle hydrodynamics

• Hierarchial treecodes

• Application to accretion discs

• References

LIST OF PARTICIPANTS