Session N21: Computational Methods: Multiscale Modeling

The Improvement of Length Scaling in the Hyperdynamics Method

Many important physical phenomena, such as film growth, bulk diffusion, radiation damage annealing, dislocation climb and catalysis, require both long time scale and large length scale molecular dynamics, where conventional molecular dynamics methods are not applicable due to the computational costs. The hyperdynamics method has enabled us to perform molecular dynamics for a longer time scale. However, this method is limited in length scale because the overall computational speedups achieved by the current bias potential methodologies decrease rapidly with the size of the system. To overcome this disadvantage, we are designing new bias potential methodologies to maintain the overall computational speedup in larger systems. We calculate the hyperdynamics rates and the overall speedups with the current and new approaches and discuss the fundamental aspects of both approaches. The early results show that these new methods are promising for reaching greater time and length scales simultaneously.   

Finding the minimum-energy atomic configuration in large multi-atom structures: Genetic Algorithm versus the Virtual-Atom Approach

In many problems in molecular and solid state structures one needs to determine the energy-minimizing decoration of sites by different atom-types (i.~e.\emph{configuration}). The sheer size of this configurational space can be horrendous even if the underlying lattice-type is known. The ab-initio total-energy surface for different (relaxed) configurations can often be parameterized by a spin-like Hamiltonian (\emph{Cluster-Expansion}) with discrete spin -variables denoting the type of atom occupying each site. We compare two search strategies for the energy-minimizing configuration: (i) A discrete-variable genetic-algorithm approach( S. V. Dudiy and A. Zunger, PRL {\bf 97}, 046401 (2006) ) and (ii) a continuous-variable approach (M. Wang et al, J. Am. Chem. Soc. {\bf 128}, 3228 (2006) ) where the discrete-spin functional is mapped onto a continuous-spin functional (\emph{virtual atoms}) and the search is guided by local gradients with respect to each spin. We compare their efficiency at locating the ground-state configurations of fcc Au-Pd Alloy in terms of number of calls to the functional. We show that a GA approach with diversity-enhancing constraints and reciprocal-space mating easily outperforms the VA approach.   

The global space-group optimization approach to crystal structure prediction

We present the global space-group optimization (GSGO) approach to the prediction of both the lattice structure and the atomic configuration of a crystalline solid. The GSGO method is based on an evolutionary algorithm within which a population of crystal structures is evolved substituting the highest total-energy structures with new ones. The search is performed directly on the atomic positions and the unit-cell vectors, after a similarity transformation is applied to bring structures of different unit-cell shapes to a common basis. Following this transformation, we can define a crossover operation that treats on the same footing structures with different unit-cell shapes. Newly generated structures are fully relaxed to the closest local total-energy minimum. Starting from random unit-cell vectors and atomic positions, and using the VASP code, the GSGO procedure found for Si, GaAs, SiC the correct lattice structure and configuration. In the case of Au$_{8}$Pd$_{4}$, the search retrieved the correct underlying lattice type (fcc), but energetically closely spaced ($\sim 2$ meV/atom) alloy configurations were not resolved. The GSGO approach opens the way to predicting unsuspected structures, using, in the cases noted above, an order of $\leq 100$ total-energy {\em ab initio} calculations.   

Quantifying cluster expansions in multicomponent systems:\ Precise expansions from noisy databases

We have performed a systematic analysis of the ubiquitous numerical errors contained in the databases used in cluster expansions of multicomponent alloys. Our results underscore the importance of numerical noise on the effective cluster interactions and on the selection mechanisms. The relevance of the size of and the information contained in the input database is highlighted. It is shown that cross-validatory approaches by themselves can produce unphysical expansions characterized by non-negligible, long-ranged coefficients. A selection criterion that combines both forecasting ability and the physical limiting behavior for the expansion is proposed. Expansions performed under this criterion exhibit the remarkable property of noise filtering. We illustrate our findings on bcc-based Fe-Co alloys.   

Diffusion, coarsening and plasticity in alloys using the phase field crystal model

The phase field crystal model describes materials at the nanoscale on diffusive time scales, and can capture elasticity, crystallography, dislocation and grain boundary dynamics, as well as solidification processes. Here we present the extension to binary alloys, taking into account vacancies. We show how the model can be applied to technologically important phenomena, such as diffusion, grain coarsening and plasticity.   

Rapid First-Principles Design Estimates of Alloy Order-Disorder Temperatures

From DFT calculations, we propose a rapid, mean-field estimate for order-disorder temperatures T$_{c}$ and phase diagrams via cluster expansion Hamiltonians $H=\sum_i V_{i} \phi_{i}$, where $V_{i}$ and $\phi_{i}$ are, respectively, the $i$-th cluster interaction and correlation function. We discuss when the estimate is valid and confirm its accuracy via Monte Carlo simulation. As the cost of Monte Carlo (MC) increases with number and size of clusters, such rapid estimates are desirable both for design and to limit the T and composition range needed for MC. We show two broad classes of systems as determined by $V_{i}$ in which T$_{c}$ is given accurately by (i) $\Delta H_{d-o}/\Delta S_{d-o}$ or (ii) $\Phi \Delta H_{d-o}$, where $\Delta H_{d-o}$ and $\Delta S_{d-o}$ are the enthalpy and entropy differences between fully disordered and ordered phases, respectively, and $\Phi $ is a lattice-topology dependent constant. With no finite-T intermediate phases, phase boundaries are found analytically by T$_{c} (x-x_{s})=\eta^{2} (x-x_{s})$T$_{c} (x_{s})$, where $\eta $ ($0\le \eta \le 1$) is the long-range order parameter and ($x-x_{s}$) is deviation from stiochiometry, $x_{s}$, found rapidly by CE ground-state analysis, and T$_{c} (x_{s})$ is from (i) or (ii). We exemplify results for several alloys in each class.   

Molecular Dynamics Simulations of Interface Failure

The mechanical integrity of silicon/silicon nitride interfaces is of great importance in their applications in micro electronics and solar cells. Large-scale molecular dynamics simulations are an excellent tool to study mechanical and structural failure of interfaces subjected to externally applied stresses and strains. When pulling the system parallel to the interface, cracks in silicon nitride and slip and pit formation in silicon are typical failure mechanisms. Hypervelocity impact perpendicular to the interface plane leads to structural transformation and delamination at the interface. Influence of system temperature, strain rate, impact velocity, and system size on type and characteristics of failure will be discussed.   

Spatial Stratification of Order As Used in Failure Analysis

Silicon nitride deposited on silicon substrates has application in dielectric layers for microelectronics as well as in photovoltaics. During production and operation of components involving silicon/silicon nitride interfaces, stresses and strains can build up at various temperatures resulting in component failure. Using molecular dynamics simulations the influence of temperature and rate of externally applied strain on silicon/silicon nitride interfaces has been analyzed. The primary purpose of this research is to understand the mechanisms leading to the failure of these films. Analyses involving bond lengths and angles have been developed to gain insight into these mechanisms. Methods for stratifying bond lengths and bond angles into unique sub-populations on the basis of spatial orientation have been developed, and have given much insight to how the material behaves, particularly with regards to the Poisson effect. Possible extensions of this stratification method to primitive rings will also be examined. In combination with experimental observations, this analysis will deepen our understanding of the structural properties of silicon/silicon nitride interfaces.   

Non-equilibrium Molecular Dynamics Study of the Thermal Resistance at the Interface Between Two Materials

Two different crystalline systems comprised of atoms interacting through Lennard-Jones (LJ) potentials were set in contact. The thermal conduction through such solid-solid interface was studied as a function of temperature and relative materials parameters, where the species differ in mass, hard-core atomic diameter and well depth. The computational setup simulated a solid sample with two different materials separated at a central interface. A non-equilibrium Molecular Dynamics approach was taken to calculate the Kapitza thermal resistance across the interface and its dependence on the two species and LJ parameters. It is found that the Kapitza resistance decreases as a function of temperature for mostly all combinations of the two materials LJ parameters.   

Graded-Sequence-of-Approximations: Quantum Mechanical Forces for Molecular Dynamics

The rate-limiting step in the multiscale simulation of materials, biomolecular, and other complex systems is quite generally the generation of the quantum mechanical (QM) forces in the chemically active region. A sequence of approximations involving both QM and classical approximations is used to reduce the computational intensity of the problem. More computationally intensive approximations of greater accuracy, are used at infrequent simulation steps to recalibrate forces from the less intensive calculations, which have lesser accuracy. The graded-sequence-of-approximations technique is illustrated in a particularly demanding case in which we have used a published classical potential for silica with QM forces generated by a quantum chemical technique independently trained to reproduce relevant coupled-cluster forces.   

Incoporating Existing Large Applications in the PUPIL System: Amber

PUPIL (Program for User Package Interfacing and Linking)$^1$ inter-operates existing codes for multi-threaded, multi-scale quantum and classical mechanical simulations via JAVA, XML, JAVA, a C++ library, and minimally intrusive wrappers for each code. An architectural challenge for PUPIL is support of modules from a multi-scale QM-MD suite with much internal coupling. We have succeeded with the AMBER suite MD module (Sander), with Gaussian03 for QM. Our demonstration study is the decomposition of Angelis' salt with explict water. A variable quantum zone (solute and first solvation cell) was used, with the remaining waters via TIP3P. Sander calculated the Potential of Mean Force for the reaction through umbrella sampling, with the QM forces from Gaussian. We summarize PUPIL architecture and implementation aspects, report efficiency and overhead measures, and discuss the computed results. $^1$J.~Torras, E.~Deumens and S.B.~Trickey, J. Computer Aided Mat. Des. {\bf 13}, 201 (2006); J.~Torras et al. Comp. Phys. Comm. 2006 [accepted]   

Concepts of Multi-Scale Modeling

The approximate representation of a quantum solid as an equivalent composite semi-classical solid is considered. In the classical bulk domain this potential energy is represented by potentials constructed to give the same structure and elastic properties as the underlying quantum solid. In a small local quantum domain the potential is determined from a detailed quantum calculation of the electronic structure.The features of this problem are the representation of the classical domain by potentials focused on reproducing the specific quantum response being studied, development of `pseudo-atoms' for a realistic treatment of charge, and inclusion of polarization effects on the quantum domain due to its distant bulk environment   

Quasi-continuum orbital-free density-functional theory (QC-OFDFT)

Density-functional theory has provided insights into various materials properties in the recent decade. However, its computational complexity has made other aspects, especially those involving defects, beyond reach. Here, we present a seamless coarse-graining scheme for orbital-free density-functional theory (OFDFT), that enables the study of multi-million atom clusters with no spurious physics and at no significant loss of accuracy. The key ideas are (i) a real-space formulation, (ii) a nested finite-element implementation of the formulation and (iii) a systematic means of adaptive coarse-graining retaining full resolution where necessary and coarsening elsewhere with no patches, assumptions or structure. Fully-resolved OFDFT and finite lattice-elasticity are obtained as special limits of this scheme. This methodological development has enabled OFDFT calculations on large systems, which have revealed interesting physics and phenomena that have not been observed to date.   

First-passage Monte Carlo for simulations of alloy microstructure

We unveil a principally new Monte Carlo algorithm for simulations of multiple diffusing particles of finite dimensions that coalesce or annihilate on collisions. The algorithm is derived from the theory of first-passage processes and a time-dependent Green's function formalism. The new method circumvents the need for long and tedious diffusion hops by which the particles find each other in space. At the same time, the algorithm is exact and its computational efficiency is astonishing. The new algorithm is generally applicable in 1d, 2d, 3d, ... and to a wide variety of important physical situations, including nucleation, growth and coarsening of alloy particles, interstitial and vacancy clusters after quench or under irradiation. We will present simulation of multi-million particle ensembles covering over 10 decades of time of microstructural evolution.   

Adaptive Resolution in Molecular Dynamics Simulations

For the study of complex synthetic and biological molecular systems by computer simulations one is still restricted to simple model systems or to by far too small time scales. To overcome this problem multiscale techniques are being developed. However in almost all cases, the regions treated at different level of resolution are kept fixed and do not allow for a free exchange. We here give a basic theoretical framework for an efficient and flexible coupling of the different regimes. The approach leads to a concept, which can be seen as a geometry induced phase transition and to a counterpart of the equipartition theorem for fractional degrees of freedom. The efficiency of the presented approach is illustrated on two numerical examples, i.e., the molecular dynamics simulations of bulk water and a generic polymer in a solvent.