Carolyn L. Phillips

Rigid body constraints realized in massively-parallel molecular dynamics on graphics processing units

Abstract Molecular dynamics (MD) methods compute the trajectory of a system of point particles in response to a potential function by numerically integrating Newtonʼs equations of motion. Extending these basic methods with rigid body constraints enables composite particles with complex shapes such as anisotropic nanoparticles, grains, molecules, and rigid proteins to be modeled. Rigid body constraints are added to the GPU-accelerated MD package, HOOMD-blue, version 0.10.0. The software can now simulate systems of particles, rigid bodies, or mixed systems in microcanonical (NVE), canonical (NVT), and isothermal-isobaric (NPT) ensembles. It can also apply the FIRE energy minimization technique to these systems. In this paper, we detail the massively parallel scheme that implements these algorithms and discuss how our design is tuned for the maximum possible performance. Two different case studies are included to demonstrate the performance attained, patchy spheres and tethered nanorods. In typical cases, HOOMD-blue on a single GTX 480 executes 2.5–3.6 times faster than LAMMPS executing the same simulation on any number of CPU cores in parallel. Simulations with rigid bodies may now be run with larger systems and for longer time scales on a single workstation than was previously even possible on large clusters.

Pseudo-random number generation for Brownian Dynamics and Dissipative Particle Dynamics simulations on GPU devices

Brownian Dynamics (BD), also known as Langevin Dynamics, and Dissipative Particle Dynamics (DPD) are implicit solvent methods commonly used in models of soft matter and biomolecular systems. The interaction of the numerous solvent particles with larger particles is coarse-grained as a Langevin thermostat is applied to individual particles or to particle pairs. The Langevin thermostat requires a pseudo-random number generator (PRNG) to generate the stochastic force applied to each particle or pair of neighboring particles during each time step in the integration of Newtons equations of motion. In a Single-Instruction-Multiple-Thread (SIMT) GPU parallel computing environment, small batches of random numbers must be generated over thousands of threads and millions of kernel calls. In this communication we introduce a one-PRNG-per-kernel-call-per-thread scheme, in which a micro-stream of pseudorandom numbers is generated in each thread and kernel call. These high quality, statistically robust micro-streams require no global memory for state storage, are more computationally efficient than other PRNG schemes in memory-bound kernels, and uniquely enable the DPD simulation method without requiring communication between threads.

Computational self-assembly of a one-component icosahedral quasicrystal

Icosahedral quasicrystals (IQCs) are a form of matter that is ordered but not periodic in any direction. All reported IQCs are intermetallic compounds and either of face-centred-icosahedral or primitive-icosahedral type, and the positions of their atoms have been resolved from diffraction data. However, unlike axially symmetric quasicrystals, IQCs have not been observed in non-atomic (that is, micellar or nanoparticle) systems, where real-space information would be directly available. Here, we show that an IQC can be assembled by means of molecular dynamics simulations from a one-component system of particles interacting via a tunable, isotropic pair potential extending only to the third-neighbour shell. The IQC is body-centred, self-assembles from a fluid phase, and in parameter space neighbours clathrates and other tetrahedrally bonded crystals. Our findings elucidate the structure and dynamics of the IQC, and suggest routes to search for it and design it in soft matter and nanoscale systems.

Stability of the double gyroid phase to nanoparticle polydispersity in polymer-tethered nanosphere systems

Recent simulations predict that aggregating nanospheres functionalized with polymer “tethers” can self-assemble to form the double gyroid (DG) phase seen in block copolymer and surfactant systems. Within the struts of the gyroid, the nanoparticles pack in icosahedral motifs, stabilizing the gyroid phase in a small region of the phase diagram. Here, we study the impact of nanoparticle size polydispersity on the stability of the double gyroid phase. We show for low amounts of polydispersity the energy of the double gyroid phase is lowered. A large amount of polydispersity raises the energy of the system, disrupts the icosahedral packing, and eventually destabilizes the gyroid. Our results show that the DG forms readily up to 10% polydispersity. Considering polydispersity as high as 30%, our results suggest no terminal polydispersity for the DG, but that higher polydispersities may kinetically inhibit the formation of phase. The inclusion of a small population of either smaller or larger nanospheres encourages low-energy icosahedral clusters and increases the gyroid stability while facilitating its formation. We also introduce a new measure for determining the volume of a component in a microphase-separated system based on the Voronoi tessellation.

An energy-conserving two-temperature model of radiation damage in single-component and binary Lennard-Jones crystals

Two-temperature models are used to represent the interaction between atoms and free electrons during thermal transients such as radiation damage, laser heating, and cascade simulations. In this paper, we introduce an energy-conserving version of an inhomogeneous finite reservoir two-temperature model using a Langevin thermostat to communicate energy between the electronic and atomic subsystems. This energy-conserving modification allows the inhomogeneous two-temperature model to be used for longer and larger simulations and simulations of small energy phenomena, without introducing nonphysical energy fluctuations that may affect simulation results. We test this model on the annealing of Frenkel defects. We find that Frenkel defect annealing is largely indifferent to the electronic subsystem, unless the electronic subsystem is very tightly coupled to the atomic subsystem. We also consider radiation damage due to local deposition of heat in two idealized systems. We first consider radiation damage in a large face-centered-cubic Lennard-Jones (LJ) single-component crystal that readily recrystallizes. Second, we consider radiation damage in a large binary glass-forming LJ crystal that retains permanent damage. We find that the electronic subsystem parameters can influence the way heat is transported through the system and have a significant impact on the number of defects after the heat deposition event. We also find that the two idealized systems have different responses to the electronic subsystem. The single-component LJ system anneals most rapidly with an intermediate electron-ion coupling and a high electronic thermal conductivity. If sufficiently damaged, the binary glass-forming LJ system retains the least permanent damage with both a high electron-ion coupling and a high electronic thermal conductivity. In general, we find that the presence of an electronic gas can affect short and long term material annealing.

Self Assembled Clusters of Spheres Related to Spherical Codes

We consider the thermodynamically driven self-assembly of spheres onto the surface of a central sphere. This assembly process forms self-limiting, or terminal, anisotropic clusters (N-clusters) with well-defined structures. We use Brownian dynamics to model the assembly of N-clusters varying in size from two to twelve outer spheres and free energy calculations to predict the expected cluster sizes and shapes as a function of temperature and inner particle diameter. We show that the arrangements of outer spheres at finite temperatures are related to spherical codes, an ideal mathematical sequence of points corresponding to the densest possible sphere packings. We demonstrate that temperature and the ratio of the diameters of the inner and outer spheres dictate cluster morphology. We present a surprising result for the equilibrium structure of a 5-cluster, for which the square pyramid arrangement is preferred over a more symmetric structure. We show this result using Brownian dynamics, a Monte Carlo simulation, and a free energy approximation. Our results suggest a promising way to assemble anisotropic building blocks from constituent colloidal spheres.

Stable large-scale solver for Ginzburg-Landau equations for superconductors

Understanding the interaction of vortices with inclusions in type-II superconductors is a major outstanding challenge both for fundamental science and energy applications. At application-relevant scales, the long-range interactions between a dense configuration of vortices and the dependence of their behavior on external parameters, such as temperature and an applied magnetic field, are all important to the net response of the superconductor. Capturing these features, in general, precludes analytical description of vortex dynamics and has also made numerical simulation prohibitively expensive. Here we report on a highly optimized iterative implicit solver for the time-dependent Ginzburg-Landau equations suitable for investigations of type-II superconductors on massively parallel architectures. Its main purpose is to study vortex dynamics in disordered or geometrically confined mesoscopic systems. In this work, we present the discretization and time integration scheme in detail for two types of boundary conditions. We describe the necessary conditions for a stable and physically accurate integration of the equations of motion. Using an inclusion pattern generator, we can simulate complex pinning landscapes and the effect of geometric confinement. We show that our algorithm, implemented on a GPU, can provide static and dynamic solutions of the Ginzburg-Landau equations for mesoscopically large systems over thousands of time steps in a matter of hours. Using our formulation, studying scientifically-relevant problems is a computationally reasonable task.

Optimal Filling of Shapes

We present filling as a type of spatial subdivision problem similar to covering and packing. Filling addresses the optimal placement of overlapping objects lying entirely inside an arbitrary shape so as to cover the most interior volume. In n-dimensional space, if the objects are polydisperse n-balls, we show that solutions correspond to sets of maximal n-balls. For polygons, we provide a heuristic for finding solutions of maximal disks. We consider the properties of ideal distributions of N disks as N→∞. We note an analogy with energy landscapes.

Discovering crystals using shape matching and machine learning

As the rate at which data can be amassed through computational simulation or experimental techniques accelerates, the pace of discovery becomes limited by the rate at which data can be analyzed. In this paper, a method is introduced by which new types of crystalline structures can be automatically identified from large data sets of coordinates. By deploying a hierarchy of pattern analysis techniques using shape matching and machine learning algorithms, local structures are extracted, classified, and then used to partition a data set into groups of common crystals. This method requires no a priori knowledge of what might be present in the data set. The method is applied to two data sets that contain both simple and complex crystals, including quasicrystals. We show how phase diagrams can be automatically generated and identify a crystal phase missed in prior analyses. By integrating shape matching and machine learning techniques to analyze rapidly produced databases, the discovery of new crystal structures and materials can be accelerated. This method is especially applicable to soft matter systems, where particle interactions can be exquisitely tuned and designed to drive the self-assembly of mesoscale materials with exotic structures.

A two-temperature model of radiation damage in α-quartz

Two-temperature models are used to represent the physics of the interaction between atoms and electrons during thermal transients such as radiation damage, laser heating, and cascade simulations. We introduce a two-temperature model applied to an insulator, α-quartz, to model heat deposition in a SiO(2) lattice. Our model of the SiO(2) electronic subsystem is based on quantum simulations of the electronic response in a SiO(2) repeat cell. We observe how the parametrization of the electronic subsystem impacts the degree of permanent amorphization of the lattice, especially compared to a metallic electronic subsystem. The parametrization of the insulator electronic subsystem has a significant effect on the amount of residual defects in the crystal after 10 ps. While recognizing that more development in the application of two-temperature models to insulators is needed, we argue that the inclusion of a simple electronic subsystem substantially improves the realism of such radiation damage simulations.