ATK-ForceField: A New Generation Molecular Dynamics Software Package

TL;DR

ATK-ForceField is a new Python-based software package for atomistic simulations using classical potentials, enabling customizable molecular dynamics, geometry optimization, and advanced modeling with demonstrated applications.

Contribution

It introduces a comprehensive, Python-integrated simulation tool with detailed implementation and performance analysis, expanding capabilities for atomistic modeling.

Findings

Modeling thermal transport in silicon germanium

Simulating vapor deposition of selenium molecules

Creep simulation in copper polycrystals

Abstract

ATK-ForceField is a software package for atomistic simulations using classical interatomic potentials. It is implemented as a part of the Atomistix ToolKit (ATK), which is a Python programming environment that makes it easy to create and analyze both standard and highly customized simulations. This paper will focus on the atomic interaction potentials, molecular dynamics, and geometry optimization features of the software, however, many more advanced modeling features are available. The implementation details of these algorithms and their computational performance will be shown. We present three illustrative examples of the types of calculations that are possible with ATK-ForceField: modeling thermal transport properties in a silicon germanium crystal, vapor deposition of selenium molecules on a selenium surface, and a simulation of creep in a copper polycrystal.

Figures

Tables

Table 1. Table 1 : Selected potential models included in ATK-ForceField.

Potential model	Special properties	Reference
Stillinger-Weber	three-body	[38]
Embedded atom model (EAM)	many-body	[58]
Modified embedded atom model (MEAM)	many-body directional bonding	[59]
Tersoff	bond-order	[44]
ReaxFF	bond-order dynamical charges	[50]
COMB/COMB3	bond-order dynamical charges induced dipoles	[48]
Core-shell	dynamical charge fluctuations	[54]
Tangney-Scandolo	induced dipoles	[60]
Aspherical ion model	induced dipoles and quadrupoles dynamical ion distortion	[61]
biomolecular and valence force fields	static bonds	[56, 57]

Table 2. Table 2 : Absolute timings for the MD runs using a single core with different potentials, normalized to one atom and a single MD step.

	LJ	Tersoff	SW	BKS	EAM	ReaxFF	COMB	TS
T (s)	$1.5\cdot {10}^{-6}$	$7.7\cdot {10}^{-7}$	$3.3\cdot {10}^{-7}$	$2.8\cdot {10}^{-5}$	$2.3\cdot {10}^{-6}$	$1.9\cdot {10}^{-4}$	$1.7\cdot {10}^{-4}$	$3.6\cdot {10}^{-4}$

Equations

$$V(\vec{x})=V_{\varnothing }+\sum _i{V}_{\{i\}}({\mathbf{r}}_i)+\sum _{i<j}{V}_{\{i,j\}}({\mathbf{r}}_i,{\mathbf{r}}_j)+\dots +V_{\{1,\dots ,N\}}({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N).$$

$$V_2(\vec{x})=\sum _{i=1}^N\sum _{j\ne i}{v}_{i,j}({\mathbf{r}}_i,{\mathbf{r}}_j),$$

$$V_2(\vec{x},L)=\sum _{i=1}^N\sum _{\mathbf{n}\in {\mathbb{Z}}^3}\sum _{\genfrac{}{}{0ex}{}{j=1}{{\mathbf{r}}_i\ne L\mathbf{n}+{\mathbf{r}}_j}}^N{v}_{i,j}({\mathbf{r}}_i,L\mathbf{n}+{\mathbf{r}}_j).$$

$$V_{ij}(r)=\frac{q_i{q}_j}{r}$$

$$V({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)=\underset{q_1,\dots ,q_N}{\min }\tilde{V}({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N,q_1,\dots ,q_N).$$

$${\mathbf{F}}_{\mathrm{simultaneous}}=(-\nabla V(\mathbf{r}),-({\sigma }_{\mathrm{cauchy}}-{\sigma }_{\mathrm{external}})\frac{\Omega }{L\sqrt{n}}),$$

$$\kappa =\frac{\dot{Q}}{A\Delta T},$$

Full text

1 Introduction

Molecular dynamics (MD) simulations have become a versatile and widely used tool in scientific computer simulations [1, 2]. MD can be used to gain insight in atomic-scale processes when experimental techniques are unable to provide sufficient resolution. MD is used to sample microscopic ensembles, to study microscopic dynamical and transport properties, such as diffusion or thermal transport, or to simulate physical processes in order to access atomic structures for which little experimental information is available, e.g. structures of complex interfaces, amorphous materials, or vapor-deposited films.

MD simulations are often used to interpret and understand the results of experimental studies, but they are also capable of predicting or screening new materials or processes which can significantly reduce the number and cost of experiments [3, 4].

Numerical optimization algorithms are another, closely related, core tool in atomic scale modeling, as locating stable structures and reaction mechanisms can be formulated as optimization problems. Minimum energy structures are obtained by performing a local minimization of the potential energy (or enthalpy for solids under isotropic pressure). Reaction mechanisms can be located by minimizing the projected forces of a nudged elastic band (NEB) [5].

A wide variety of software packages are available for MD and optimization. Some packages like LAMMPS [6], DL_POLY [7], or MBN Explorer [8], are designed to work efficiently with a wide spectrum of materials, however, most are optimized for certain classes of systems, such as crystals (e.g. VASP [9] and GULP [10]), biomolecules (e.g. AMBER [11], GROMACS [12], CHARMM [13], and NAMD [14]), or polymers (e.g. ESPResSo [15]). Almost all of these packages only support either empirical potentials or ab-initio methods. For many tasks it may be desirable to quickly and seamlessly switch between different levels of accuracy, e.g. when validating or parameterizing empirical potentials against higher-level calculations, or when using classical potentials to obtain mechanical and thermal properties and DFT for electronic properties of the same system [16, 17]. Yet, very few (e.g. CP2K [18] or ADF [19]) simulation codes provide forces and energies from both of ends of the atomic accuracy spectrum. The atomic simulation environment (ASE) [20] addresses this problem by integrating a broad range of atomic-scale simulation methods into a common Python framework which implements various tools for atomic scale simulations, including molecular dynamics and optimization functionality. However, the ASE is not optimized towards highly efficient classical MD and is therefore mainly used in combination with DFT methods.

Some physical problems require very specific techniques which reach beyond the standard algorithms commonly implemented in most MD codes or may even need to be custom-tailored to the application at hand. In this case, a flexible way of extending and combining the standard simulation algorithms is desired. Ideally this task should be done as seamlessly as possible, without having to manually interface different programs, programming languages, or converting input and output files from one format to the another. This concerns not only the simulation technique itself, but, equally important, the subsequent analysis and visualization of the simulation trajectory and other raw data. Although most major software packages come with a more or less extensive suite of analysis and visualization tools, these are often included as small standalone programs, which have to be compiled and invoked individually. This often leads to each user developing a custom workflow by patching together different tools, which significantly impedes reusability and reproducability.

As MD becomes a more important and commonly used technique, it is important to make it accessible to a broader spectrum of scientists. Thus there is a need for user-friendly ways of setting up and analyzing simulations. This becomes particularly important in industrial settings, where the time a researcher needs to become an efficient user of a new code impacts the business’ bottom line.

In this paper, we describe how Virtual NanoLab (VNL) and the Atomistix ToolKit (ATK) developed by QuantumWise are designed to meet these challenges. Generally, the main goal of VNL and ATK is to make atomistic simulation techniques easily accessible to an increasing number of researchers. To this aim, the VNL provides a graphical user interface, in which all steps involved in the workflow of atomistic simulations, i.e designing the configuration to be simulated, setting up and running the actual calculation, as well as analyzing and visualizing the results, can be carried out.

ATK provides a Python scripting interface to all the calculators and simulation methods. ATKPython is inspired by ASE [20], however, it is more extended and rich in functionality and specifically with respect to MD simulations there is a negligible overhead from the Python layer. It is designed in a modular way, meaning that individual components can easily be exchanged without having to change the entire workflow. For example, switching from an empirical potential to DFT, only requires a single line change in the input script. Most importantly, by using Python, it is possible to for the user to implement custom simulation techniques without having to re-compile the code.

For the scope of this paper we will specifically focus on the ATK-ForceField package and the dynamics and optimization module in ATK, although, in principle all scripts and examples could as well be done with DFT methods. ATK-ForceField is the calculation engine for empirical potentials. It mainly consists of the TremoloX-calculator, which is developed by the Fraunhofer Institute for Algorithms and Scientific Calculations (SCAI). The TremoloX- calculator is based on SCAI’s Tremolo-X software package [21] for numerical simulation in molecular dynamics and mechanics [2]. The TremoloX-calculator is optimized for large scale simulations, provides state-of-the-art scalable fast methods for long-range interactions [22] and has successfully performed in several applications, ranging from nanoparticles and nanocomposites over cementitious materials to electrolytes and biomolecules [23, 24, 25, 26, 27, 28, 29, 30, 31, 32]. TremoloX itself is implemented in parallel, although the ATK-ForceField currently only supports the shared-memory OpenMP version. The combination of the ATK-ForceField calculator module and the dynamics module in ATK provides state-of-the-art MD functionality and optimization algorithms, as well as more advanced techniques such as nudged elastic band (NEB) and adaptive kinetic Monte Carlo (AKMC) simulations.

Although in principle VNL and ATK are commercial software packages, VNL and ATK-ForceField are freely available for academic groups. We hope that in this way we can develop an academic community of end users which will develop new application examples for the ATKPython framework and new plugins for VNL. In fact, all examples presented in section 9 of this paper can be run under this free license.

In the first few sections, we describe the basic concepts of modeling atomic geometries and how interatomic potentials are composed and calculated inside ATK-ForceField, as well as some performance benchmarks. The next section reviews the basic concepts of MD and optimization, and their implementation in ATK to meet the specifications defined above. Subsequently, the simulation and analysis workflow in VNL and ATK is presented. We then describe the integration of the calculation engine into the Python platform, including some small code samples. In the following section, we present some example simulations to demonstrate the capability of MD using ATK. Finally, we give an overview for the future development of ATK-ForceField.

2 Atomic Configurations

At the beginning of each calculation the system under investigation must be specified. ATK offers several configuration types for different applications. The simplest case is a molecule configuration which describes a system containing an isolated collection of atoms with their positions $\mathbf{r}_{i}$ and elements $Z_{i}$ (see figure 1 (d)). For systems like crystals which are are inherently periodic and require a different treatment, ATK offers the so-called bulk configuration class. A bulk configuration contains the central cell and repeats it periodically in all directions during simulations (see figure 1 (b)). This simulation cell can be a unit cell, e.g. of a crystal (figure 1 (b)), but it can also be used as a supercell, containing systems without a pronounced periodic character, such as e.g. amorphous, liquid, or interface systems, as well as slab geometries (figure 1 (c)). In ATK, the three lattice vectors $\mathbf{l}_{i}$ spanning the cell can either be given explicitly, as a cell matrix $\mathbf{L}=[\mathbf{l}_{1},\mathbf{l}_{2},\mathbf{l}_{3}]$, or by using one of the predefined classes that represent the 14 different Bravais lattice systems, which require only the independent lattice parameters to be specified.

Apart from bulk and molecule configurations, ATK also offers so-called device (two-probe) configurations, which consist of a central region, as well as a left and right electrode. A device configuration is periodic in two directions, whereas the third dimension is bounded by two infinite electrodes (see Fig. 1 (e)). The central region and electrode parts are each defined as bulk configurations, and, to be a valid definition, the left and right part of the central region must be identical to the corresponding adjacent electrode. These configurations can be used to simulate electron or phonon transport via the non-equilibrium Green’s functions (NEGF) technique [33]. Finally, to accurately study surface chemistry, ATK provides surface (one-probe) configurations, which are device configurations with only one electrode (see Fig. 1 (f)), accurately describing a semi-infinite surface geometry beyond the common slab supercell approximation.

3 Interatomic Potentials

In general, the underlying model to describe the dynamics of an atomic systems is the full Schrödinger equation for the electrons and nuclei under consideration. Its high-dimensionality, however, makes a direct numerical treatment impossible and thus one has to resort to model approximations. A hierarchy of these approximations is typically necessary to calculate the properties of nano-scale systems. Based on the Born-Oppenheimer formulation, which disentangles the motions of the nuclei and the electrons into a classical and a quantum mechanical part [34, 2], the solution to the electronic Schrödinger equation is commonly approximated by density functional theory (DFT) or tight-binding (TB) techniques, and a variety of such methods is readily available in the ATK program package.

Although these QM methods have become very powerful, the accessible systems sizes are still not sufficient for a realistic simulation of many nano-structured systems of practical interest. To this aim, the idea is to determine an approximation $V$ to the exact potential energy surface in an efficient analytical functional form, a so-called force field or classical potential.

To overcome the curse of dimensionality of $V$ being a function of all atomic positions, classical potentials are typically based on an appropriate many-body expansion, i.e.

[TABLE]

If the size of the terms decays fast with e.g. the order $k=|u|$ of the $k$-body contributions $V_{u}$, then a proper truncation of the expansion results in a substantial reduction in computational complexity.

A very simple class of classical potentials are two-body potentials which have the form

[TABLE]

where $v_{i,j}$ is a function that depends only on the properties of two particles. A variety of well-known potential classes, such as the Lennard-Jones potential [35] or the Buckingham potential [36] can be written in the form of equation (1). In the periodic case, interactions between the particles in the unit cell and all other particles, including their periodic copies, are taken into account:

[TABLE]

In a similar fashion, three-body potentials and straightforwardly higher-order $k$-body potentials can be defined, where the periodic case can be treated analogously to two-body potentials [2, 37]. An example for a widely used three-body potential – in combination with two-body terms – is the Stillinger-Weber potential [38].

If the interaction terms $v_{u_{1},\ldots,u_{k}}$ are of short range, meaning that the error that results from discarding terms involving particles further apart than a cutoff distance $r_{cut}$ can be made arbitrarily small by choosing $r_{cut}$ sufficiently large [39], the total potential function at any fixed $\vec{x}$, involves only $O(N)$ terms. Most interaction potentials applied within an approximation potential of the electronic energy $V_{el}$ are of short range. However, in particular Coulomb interactions, i.e.

[TABLE]

or dipolar interactions do not satisfy the short range condition. Hence, the simple cutoff approach can not be applied and the interactions with all periodic images need to be taken into account [39]. To this end, a long range potential is typically decomposed in a short range but non-smooth part and a long range but smooth part which leads to methods with costs of order $O(N^{\frac{3}{2}})$ and even of order $O(N\log(N)^{\alpha})$, where $\alpha\leq 0$ depends on the specific method [2, 40, 41, 42, 22].

Besides a variety of standard two- and three-body potentials, ATK-ForceField implements a number of bond-order potentials or cluster functionals [34], like Abel and Tersoff [43, 44], embedded atom method (EAM) and modified EAM (MEAM) potential types [45, 46] or the Sutton-Chen potential [47]. Here, the expansion contributions $V_{u}$ do not depend on the particles $\mathbf{r}_{u_{1}},\ldots,\mathbf{r}_{u_{k}}$ alone, but in addition also on their neighboring environment.

A further class of more advanced potential models implemented in ATK-ForceField, are potentials in which the partial charges of the particles are not fixed but are calculated by minimizing the potential energy with respect to these charges:

[TABLE]

This includes all charge-optimized many-body (COMB) potentials [48] and ReaxFF potentials [49, 50], as well as ReaxFF+, a variant of ReaxFF, which combines charge equilibration with the bond-order principle to better distinguish between covalent and ionic bonds [51]. The equilibrium charges are obtained by using a nonlinear optimization algorithm. A similar ansatz is chosen for potentials with induced dipoles where the dipoles are once again calculated by minimizing the potential energy with respect to these dipoles. Induced dipoles are used in the third-generation COMB potentials [52] and in the aspherical ion model potential [53].

Moreover, to describe ionic polarizability in ionic materials ATK-ForceField includes also the widely used core-shell potential models [54]. Here, in contrast to the minimization problem (4), the potential energy is minimized with respect to the position of charged shells.

For a wide range of applications, it is sufficient to consider molecular systems which retain a static bond topology during simulation. A common example for these kind of force fields are biomolecular forces fields, which, apart from basic Lennard-Jones and Coulomb interactions, typically involve potentials to model bond stretching, bond angles, dihedral torsions, and in some cases improper torsions (see e.g. [55, 56]). In contrast to the non-bonded potentials, these bonded potentials only act between atoms which are connected via the bond topology. The same holds for valence force fields (VFF), which use similar bonded potential terms and which are typically used to describe elastic and vibrational properties in covalent solids [57].

A summary of the main included potential model types is given in Table 1. Several hundred pre-defined parameter sets for these potentials are currently provided in ATK-ForceField 111For a list of most provided parameter sets e.g.

http://docs.quantumwise.com/manuals/atkforcefield.html\#tremolox-potential-parameter-sets

or http//www.tremolo-x.com/en/tremoloxpotentials.html.

In some applications combinations of different potential models are required to provide an appropriate description of the considered atomic systems. Apart from a few exceptions, ATK-ForceField allows arbitrary combinations of different potentials. Moreover, if a specialized pair potential form is not implemented in ATK-ForceField, one can always resort to a tabulated representation of the potential.

Since ATK originates from an ab-initio based code, a major field of application is the simulation of small crystal unit cells. Thus, special emphasis has been placed on implementing ATK-ForceField such that it correctly calculates the interactions in cells much smaller than the interaction range of the potential, by including as many periodic images as necessary [37], but at the same time efficiently handles simulations of large systems without loss of performance.

4 Performance

In this section we will present the performance of MD simulations with different potentials in ATK-ForceField. Let us first review some implementation details of the ATK-ForceField back-engine Tremolo-X. Similar to other performance-intensive parts in ATK, the Tremolo-X back-engine is implemented in C and C++. All short range potential terms are computed with linear cost complexity $O(N)$ based on the combination of the well-known Verlet list and cell linked list algorithms [62, 1]. In this hybrid approach, the cell linked list technique is used to setup and update (when necessary) the Verlet pair lists, which are then used to efficiently compute all short range potentials. To improve the memory hierarchy performance a Hilbert-Peano curve data ordering can be applied [63, 2]. To speed up the evaluation all short range two-body potential terms can be automatically aggregated and tabulated for all pairs of particle types, respectively. To compute long-range Coulomb interactions in the periodic case (i.e. BulkConfiguration), we implemented e.g. the conventional Ewald summation [64], the smooth particle mesh Ewald [41] and the DSF approach [42], with cost complexities of order $O(N^{\frac{3}{2}})$, $O(N\log(N))$ and $O(N)$, respectively. In the case of non-periodic boundary conditions (i.e. MolecularConfiguration), we implemented a fast multipole method (FMM) [65] with linear cost complexity $O(N)$.

Tremolo-X itself is implemented in parallel using a hybrid distributed (MPI) and shared memory (OpenMP) approach. Currently, ATK-ForceField only provides the shared memory OpenMP version and in the following we limit our performance study to computations on a single node using OpenMP. The simulations were run on an Intel E5-2687W 3.1GHz CPU with 8 cores.

Six different potentials were benchmarked. We start with a simple Lennard-Jones (LJ) system, built by Argon atoms on a slightly randomized lattice with a particle density of 0.84, given in units of the LJ-minimum distance $\sigma$. Such a system is often used as benchmark system for the performance of force-field MD simulations [6]. We examine a bulk configuration of a fcc-silicon crystal, which is simulated using a Stillinger-Weber (SW) potential [38] and a Tersoff potential [44]. As an example for a system with electrostatic interactions we simulate a SiO2 cristobalite crystal, using the potential proposed by van Beest, Kramer, and van Santen (BKS) et al. [66]. The electrostatic interactions are calculated via the SPME algorithm using a cutoff of 9 Å, whereas the dispersion interactions are truncated at 10 Å. To benchmark the performance of potentials based on the embedded-atom-method (EAM), we simulate a copper crystal, using the potential of Mishin et al. [58]. Finally, we consider an even more complex, reactive ReaxFF-potential. Here, we simulated a triaminotrinitrobenzene (TATB) crystal, and we employed the potential of Strachan et al. [67], which was designed for such systems. All system comprise around 32 000 atoms. The simulations were run for 100 MD steps and the total wallclock time $T$ for the MD loop is measured. All simulations use a time step of 1 fs, except the ReaxFF simulation which used a time step of 0.2 fs.

Figure 2 shows the acceleration factor $T(1)/T(N)$ with respect to the simulation time on a single CPU. The ideal scaling behavior is depicted by the dashed line. Almost all potentials show good scaling on up to 8 threads. Even the complicated ReaxFF-potential scales very well in this range. Only the EAM potential exhibits a scaling behaviour which deviates somewhat more from the ideal speedup.

The absolute times of these potentials are given in Tab. 2, extended by the some more elaborate potentials, such as the COMB [48] and Tangney-Scandolo (TS) [60] potential. For both potentials the SiO2 cristobalite system was used with a time step of 0.2 fs. These simulations were run on a single core of a laptop computer with an Intel i7-3740QM 2.7 GHz processor. The time values are normalized to one atom and a single MD step. Although a direct comparison absolute timings is always difficult, these values are very similar to other state of the art simulation codes such as LAMMPS [6]222http://lammps.sandia.gov.

5 Molecular Dynamics

The ATK-ForceField package is designed and implemented to work highly efficiently with ATK’s molecular dynamics and optimization module. In this section we briefly review the main concepts and implementation details behind molecular dynamics in ATK, while in the next section we will discuss static approaches.

Essentially, molecular dynamics is a numerical integration to solve Newton’s equations of motion, which is typically performed by discretizing the time into small time steps $\Delta t$. The most common way to achieve this, is the velocity-Verlet algorithm [68]. All MD methods in ATK are based on this integration scheme. The resulting trajectory $\{\mathbf{r},\mathbf{v}\}(t=0,...,T_{\rm final})$ represents the motion of the particles in the microcanonical (NVE) ensemble, starting from the given initial conditions. To simulate physical processes and experiments as precisely as possible, one can additionally include external influences, such as coupling to a heat bath or applying pressure or stress, which results in different ensembles.

ATK implements several algorithms for such thermostats and barostats. The preferred types for weak coupling are the chained Nosé-Hoover thermostat [69], as well as the barostat proposed by Martyna et al. [70] for isotropic and anisotropic pressure coupling. Both integrators have been implemented based on the algorithm proposed by Tuckerman et al. [71]. These methods typically yield stable trajectories for a wide range of conditions and systems, and they exactly reproduce the canonical ensemble. For completeness, ATK also provides a Berendsen thermostat and barostat [72], although these methods only approximately reproduce the canonical ensemble. For tight temperature coupling ATK offers an impulsive version of the Langevin thermostat as described in Ref. [73]. All thermostats in ATK can optionally be used to heat and cool the system by specifying a linear increase or decrease of the reservoir temperature over the simulation time. In contrast to the force calculation, the integrators are implemented in Python, which can be efficiently accomplished in a vectorized manner using the Numpy package [74].

The essential functional blocks in a typical MD loop in ATK are schematically depicted in Fig. 3. One of the powerful features of ATK is that the MD loop works in the same way with any calculator that is attached to the given configuration. That means, instead of the ATK-ForceField calculator one could for instance employ a DFT or a semi-empirical calculator, to achieve more accurate results at higher computational cost, without having to change the MD block of the simulation.

To achieve more flexibility with respect to specialized simulation techniques, the MD loop in ATK provides the possibility of defining hook functions, which interface with the simulation at specified points during each integration loop (cf. Fig. 3). This way one can easily employ pre-defined or user-defined custom operations implementing more advanced simulation techniques. To this aim, ATK provides two kinds of hook functions, which are called either before or after the force calculation. The former type, called pre_step_hook, may modify e.g. positions, cell vectors, or velocities of the current configuration, whereas the latter type, called post_step_hook, can be used to modify the forces and stress. Thus, a pre_step_hook can for instance be used to implement custom constraints, apply strain to the simulation cell, or rigidly move atoms. The post_step_hook can be used to add external contributions or to implement bias potentials to the regular interaction forces. For example the system can be biased to explore configuration space more rapidly, as in metadynamics. In ATK-2017, the metadynamics technique is included via a post_step_hook interface to the PLUMED package [75]. Besides modifying the current state of the system, hook functions can be employed to perform and store custom on-the-fly measurements at arbitrary intervals. One advantage of this hook function framework is that the user can directly implement such hook functions in the Python simulation script without having to re-compile the code. Moreover, due to the standardized API there is no the need for the user to learn specific implementation details of the MD loop. ATK is shipped with some pre-defined hook functions, implementing e.g. thermal transport via reverse non-equilibrium molecular dynamics (RNEMD) [76].

It is often useful to define constraints on the atomic or other generalized coordinates during an MD simulation. Two of the most commonly encountered constraints are implemented in ATK: fixing the Cartesian coordinates of an atom and allowing a group of atoms to only move rigidly. It is also possible for the user to define their own custom constraints. Both the predefined and custom constraints are evaluated in the same part of the MD loop (cf. Fig. 3).

6 Optimization

6.1 Optimization Algorithms

There are two optimization algorithms implemented in ATK: the fast inertial relaxation engine (FIRE) [77] and the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) quasi-Newton optimizer [78].

The fast inertial relaxation engine (FIRE) uses a modified version of molecular dynamics to locate a minimum energy configuration. The velocities at each step are a weighted combination of the velocity at the previous step and the current forces. The time step and the weighting factor, on the previous step’s velocities, increase as long as the velocities continue to point in a descent direction (i.e. the angle between the velocity and the forces does not exceed 90 degrees). If the velocities are not a descent direction, then the timestep and weights are reset to an initial safe value. For a complete description of the algorithm see Ref. [77].

L-BFGS is a quasi-Newton optimization algorithm, which unlike other quasi-Newton methods, is able to treat problems of many variables because the memory usage and computational cost at each step are controllable [79, 78]. The implementation of L-BFGS in ATK makes use of an automatic preconditioning scheme to determine the initial inverse Hessian at each step. In the first step, the second derivative along the force is estimated using finite differences. In all subsequent steps, it is chosen according to Eqn. 2.4 in Ref. [80]. The preconditioning scheme typically ensures that the proposed L-BFGS step reduces the function value without the use of a line search. If the function value increases, then a backtracking line search is used to ensure that progress is made.

L-BFGS is superior to FIRE for most optimization problems. In one benchmark study on Lennard-Jones clusters, FIRE required 3.5 times more steps on average than L-BFGS to reach convergence [81]. In another benchmark comparing optimization algorithms for finding minimum energy pathways, L-BFGS required 1.5–3.9 times fewer iterations than FIRE for NEB calculations on surfaces [82].

6.2 Geometry Optimization

ATK features geometry optimization for isolated and periodic systems. For molecules and clusters, geometry optimization involves minimizing the atomic forces. While, for crystals, the shape and size of the unit cell as well as any external pressure (possibly anisotropic) must also be taken into account.

The approach to optimizing crystals in ATK is based upon Ref. [83]. The main idea is to express changes to the system as a combined vector of atomic and strain coordinates. The strain expresses changes to the lattice vectors relative to the initial configuration. By combining the coordinates together into a single vector it is possible to simultaneously optimize the atomic positions and the lattice vectors. The optimization algorithm is given a combined force vector which will bring the system into hydrostatic equilibrium,

[TABLE]

where $-\nabla V(\mathbf{r})$ is the usual atomic force vector, $\mathbf{\sigma}_{\mathrm{cauchy}}$ is the internal stress tensor (calculated from the atomic interactions), $\mathbf{\sigma}_{\mathrm{external}}$ is the external stress tensor, $\Omega$ is the volume of the initial unit cell, $L$ is the average distance between atoms in the initial unit cell, and $n$ is the total number of atoms. The factor of $\Omega/L\sqrt{n}$ is included as a rescaling coefficient to ensure that the stress, which has units of pressure, scales similarly to the atomic forces. The rescaling coefficient is calculated during the first optimization step and kept constant during subsequent steps.

6.3 Reaction path optimization

ATK has a state of the art NEB [5] and climbing image NEB [84] implementation that includes a newly developed algorithm for generating initial guesses [85], rapidly convergent L-BFGS optimization, and per-image parallelization for scaling to large system sizes.

The initial set of replicas can be obtained in one of two ways: linear interpolation between the end points or the image dependent pair potential (IDPP) [85]. As The IDPP avoids unphysical starting guesses, this method often leads to an initial band that is closer to the minimum energy path (MEP) than a linear interpolation and it can typically reduce the required number of optimization steps by a factor 2.

In some implementations the projected NEB forces for each image are optimized independently, however, in that case the L-BFGS algorithm is known to behave poorly [82]. In ATK, the NEB forces for each image are combined into a single vector: $\mathbf{F}_{\mathrm{NEB}}\in\mathbb{R}^{3mn}$, where $m$ is the number of images and $n$ is the number of atoms. This combined approach is much more efficient when used with L-BFGS, and has been referred to as the global L-BFGS (GL-BFGS) method [82].

6.4 Adaptive Kinetic Monte Carlo

Adaptive kinetic Monte Carlo (AKMC) is an algorithm for modeling the long timescale kinetics of solid-state materials [86, 87]. The goal is to construct a kinetic Monte Carlo (KMC) simulation by automatically locating the relevant reaction mechanisms and calculating the reaction rates using harmonic transition state theory (HTST) [88].

For a given configuration AKMC involves 3 steps: (1) locate all kinetically relevant product states; (2) determine the saddle point between the reactant and product states; (3) select a reaction using Monte Carlo.

In ATK step 1 is performed using high temperature MD. At regular intervals the MD simulation is stopped and a geometry optimization is performed. If the geometry optimization converges to a new geometry then a reaction has occured. This procedure is repeated until all relevant reactions are found (within an user specified confidence). For details of the confidence calculation see refs. [87, 89]. The saddle point geometry for each reaction is then determined by performing a NEB calculation for each reaction. A reaction is then selected using MC and the system evolves to the corresponding product configuration and the entire procedure repeated. Further details of the AKMC implementation will be presented in a separate publication [90].

7 The ATKPython Platform

In order to facilitate the integration of physical simulation methods in Python, ATK provides ATKPython, an extension of the regular Python package. Beyond the standard Python functionality, ATKPython automatically imports the relevant modules which implement all ATK-functionality. As mentioned in the previous section, the back-engines of ATK-ForceField and other performance intensive parts, such as DFT calculators, are efficiently implemented in C and C++. For a seamless integration, ATKPython uses Python front-end classes, which automatically call interface modules to these back-engines. Furthermore, ATKPython provides a built-in PhysicalQuantity module to conveniently handle and convert physical quantities including their units. The implementation of physical quantities is based on the Numpy-array structure, and ATKPython provides support for almost all native Numpy-array operations with physical quantities. Configurations are intuitively set up by specifying the coordinates and elements of the atoms, as well as details of the simulation cell, for bulk, surface and device configurations. When using bonded force fields, the bond connectivity of the configuration can be explicitly specified or automatically detected based on interatomic distances and covalent radii.

Calculators are defined by creating an instance of the desired calculator class, e.g. a LCAOCalculator() for DFT with linear combination of atomic orbitals (LCAO). When constructing the objects, parameters can be given to specify various calculator properties. A simulation using ATK-ForceField can be set up by specifying a pre-defined potential set object and initializing the TremoloXCalculator object with this potential set, e.g.

⬇

potential_set = StillingerWeber_Si_1985()

calculator = TremoloXCalculator(parameters=potential_set)

which invokes the Stillinger-Weber potential for silicon [38]. If the desired force field is not listed among the pre-defined potentials, one can set it up manually from the provided classes of interaction functions in TremoloX. To this aim, one first specifies an empty TremoloXPotentialSet object, which acts as container for the various components of the potential.

⬇

potential_set = TremoloXPotentialSet(’New potential’)

Then, one needs to add the particle type objects, one for each element contained in the system, e.g.

⬇

potential_set.addParticleType(ParticleType(

symbol=’O’,

charge=-0.5*elementary_charge

))

One can choose to add different properties to the each type, e.g. partial charges or the atomic Lennard-Jones coefficients. To distinguish between different particle types of the same element, as required e.g. in biomolecular force fields [56], particle types can be created such that they represent only a specified group of atoms, marked by a given selection tag. In the next step, one needs to set up objects for each potential function that should be included in the potential. The objects are typically initialized with the element or particle type symbols between which the potential should act, the potential parameters themselves, and some technical parameters, such as inner and outer cut-off radii. An ATKPython object defining a Morse potential between silicon and oxygen, for instance, could be set up as follows.

⬇

morse_potential = MorsePotential(

’Si’,

’O’,

r_0=2.1*Ang,

k=2.0067*1/Ang,

E_0=0.340554*eV,

r_i=6.0*Ang,

r_cut=7.5*Ang),

)

potential_set.addPotential(morse_potential)

Each of these potential function objects needs to be added to the TremoloXPotentialSet object. If electrostatic interactions should be taken into account, a coulomb solver can be added, which by default uses the partial charges specified with the particle types, although atom-specific charges can be given as well. To be used in a simulation, the calculator object must be attached to the configuration object.

Various properties, such as the total energy, the forces, or the stress, can be calculated via analysis objects, which return the corresponding properties as a physical quantity in the appropriate units independent of which type of calculator has been selected. Beyond the basic properties, ATK also provides more complex analysis objects to obtain, e.g. the phonon bandstructure, elastic constants, or local atomic stress of a given configuration with attached calculator. Some of the analysis objects, such as electron bandstructure or electron transmission spectrum, depend on the electronic structure, meaning that they are only compatible with quantum mechanical or semi-empirical calculators.

An MD simulation is set up by first defining the integration method, i.e. the ensemble and the algorithm that is used to simulate that ensemble, as an object of the desired integrator class. Here, depending on the chosen method, various thermostat parameters, as well as the initial velocities, can be invoked when setting up the integrator object.

⬇

method = NVTNoseHoover(

time_step=1*femtoSecond,

reservoir_temperature=300*Kelvin,

thermostat_timescale=100*femtoSecond,

heating_rate=0*Kelvin/picoSecond,

initial_velocity=MaxwellBoltzmannDistribution(300*Kelvin),

)

If no parameters are specified, default values are used, which generally work well for a large range of simulation purposes. Instead of a global temperature value, all NVT methods also accept a list of (tag_name, temperature)-tuples, in order to apply local thermostats to the group of atoms marked by the corresponding tag.

Finally, the MD simulation over the desired number of integration steps is started by calling the MolecularDynamics() function. Snapshots of the system are stored every log_interval steps and written to a trajectory file in HDF5, NetCDF, or XYZ format. The entire MD block in the Python script may look like the following:

⬇

md_trajectory = MolecularDynamics(

configuration,

constraints=[FixAtomConstraints([0, 1])],

trajectory_filename=’Silicon_NVT.hdf5’,

steps=5000,

log_interval=100,

pre_step_hook=None,

post_step_hook=None,

method=method

)

When all MD steps are done, the resulting md_trajectory object will contain all saved images and its associated properties, such as potential energies. Thus, the simulation script can be arbitrarily extended, by analyzing the md_trajectory object, using built-in or custom MD analysis objects. Note that in the HDF5-format the MD trajectory object is not stored in memory but as a file on disk. When quantities are extracted from the final trajectory, only the queried data is loaded into memory, which means that even large trajectories, whose size exceeds the total amount of available memory, can be analyzed.

To implement advanced functionality via hook functions as discussed in section 5, one has to define the hook, either as a function, or as a callable Python class. The parameters handed to such a function or method must be the step number, the current simulation time, the configuration, and the atomic forces and the stress. The configuration object, as well as the arrays containing forces and stress can be modified inside the function and the simulation will continue with the new values. Simple hook functions could look as in the following examples. Here, the pre_step_hook shifts the z-component of the position of the first atom, while the post_step_hook tethers the same position to a fixed point using a harmonic bias potential.

⬇

def pre_hook_example(step, time, conf, forces, stress):

positions = conf.cartesianCoordinates()

# Add an offset to the z-position of the first atom.

positions[0, 2] += 0.1*Ang

conf.setCartesianCoordinates(positions)

def post_hook_example(step, time, conf, forces, stress):

# Add a bias force to the force acting on the first atom

# in the z-direction.

positions = conf.cartesianCoordinates()

bias_force = -0.5*eV/Ang**2*(positions[0, 2] - 2.0*Ang)

forces[0, 2] += bias_force

Further examples for hook functions can be found in Sec. 9.

8 Workflow and Analysis

The simulation workflow of typical MD simulation can be set up entirely in the VNL graphical user interface. VNL is based on the ATKPython language, and each tool in VNL can interpret and generate Python scripts, thus, it is possible to seamlessly shift from the GUI to the scripting language. Most of the main functionality in ATKPython is also available in VNL. VNL is developed around a plugin concept which makes it easy to extend and add new functionality. Plugins can be downloaded and installed from an add-on server and the majority are available as source code, thus, they are easy to modify or extend with new user defined functionality.

A typical workflow to set up an MD simulation in VNL is presented in Fig. 4 with screenshots of the corresponding tools. It starts with setting up the atomic configuration in the Builder tool, using the Builder databases, combined with Builder plugins to e.g. repeat, cleave, or strain the structure, add, remove, or shift atoms, or build interfaces between different structures. To set up the ATKPython simulation protocol, the configuration is sent to the Script Generator, where a calculator, an optimization or one or more MD simulations, as well as various analysis objects can be added and their parameters set. The simulation script can either be sent to the Job Manager, from which the calculation is run on the local machine or on a remote cluster, or to the Editor. The Editor allows for modifying or extending the Python script before the calculation is run.

VNL provides various graphical analysis tools, which can be used to visualize and analyze the simulation with respect to a broad range of properties. At first, VNL provides the Movie Tool, which shows a movie of the trajectory and simultaneously displays the corresponding graph of potential, kinetic, and total energy, as well as the system temperature over simulation time. High quality graphics and more advanced visual modifications, such as changing the graphical properties of the atoms, e.g. by coloring the atoms by forces or velocities, can be carried out in the Viewer tool. Quantitative analysis of the simulation can be carried out using the various analysis objects in VNL, based on the final structure after the MD or optimization run. Here, the user can for instance calculate the local crystal structure or the local stress of each atom, and visualize the results by coloring the atoms in the Viewer. To extract averaged or time-dependent properties from the MD trajectory itself, the MD-Analyzer can be used. This tool provides a broad range of analysis quantities, such as radial distribution function, neutron scattering function, angle and coordination number distribution, vibrational density of state, mean-square displacement, void-size distribution and many more.

9 Simulation Examples

In this section we provide examples of some MD simulations using ATK-ForceField. All these examples can be run under the free academic license version of ATK.

9.1 Interfacial thermal conductance

Thermal conductance through an interface can be simulated in different ways. Two common approaches are the non-equilibrium Greens functions (NEGF) technique to calculate the phonon transmission spectrum [16] and non-equilibrium molecular dynamics (NEMD) based simulations in which a thermal flux is imposed in the system [76, 91]. Both approaches have different advantages and shortcomings. NEGF based phonon transmission only accounts for the elastic scattering of phonons and neglects contributions from phonon-phonon scattering. While NEMD simulations include both elastic and inelastic scattering, the technique is often limited by small simulations cells, which not only neglects the effects of long-wavelength phonons [92] but also introduces temperature gradients that are larger than the experimental conditions. Moreover, quantum effects, which become important at low temperatures, are neglected. ATK provides a unique possibility to use both techniques with the same calculators, and thus to exactly compare the results.

In this example, we compare both techniques by considering a simple model interface system: a silicon crystal with a single layer of germanium atoms. The atomic interactions are modeled using a Tersoff potential [93]. For the NEGF-based phonon transmission, we use a device configuration as displayed in Fig. 5 (a). The phonon transmission is calculated as described in Ref. [16] and evaluated to obtain the thermal conductance as function of temperature. The result is shown as the solid line in Fig. 5 (d).

The NEMD simulations are performed on a bulk configuration containing a slab geometry with free surfaces (cf. Fig. 5 (b)). The terminating unit cell layers on each side are used the heat source and heat sink, similar to the approach used in Ref. [94].

Note the larger silicon crystal leads in the NEMD simulation compared to the NEGF simulation. This is because in MD simulations the cell has to explicitly accommodate all phonon modes which should be included in the thermal conductance calculations, whereas in the NEGF approach the bulk character of the semi-infinite leads is taken into account via the the q-points.

After equilibrating the system in the NVT ensemble, the NEMD simulations are carried out in the NVE ensemble using a time step of 0.5 fs. To impose a thermal flux, we employ the reverse NEMD algorithm [76, 95], which works via exchanging the momenta of the hottest atom in the heat sink and the coldest atom in the heat source. This NEMD functionality is readily available in ATK as a pre-defined hook class, leaving only few parameters such as the exchange interval and the tag-names of the heat source and sink regions to be specified.

The ATKPython script for an NEMD simulation is shown in the A.1.

The temperature profile is calculated over a period of 300 ps, after discarding the initial 200 ps and plotted using the MD-Analyzer tool in VNL (see Fig. 5 (c)). At the location of the germanium layer a pronounced temperature drop $\Delta T$ is visible. The interfacial thermal conductance is obtained, based on Fourier’s law, via

[TABLE]

where $\dot{Q}$ is the average thermal flux, which is printed at the end of the simulation, and $A$ is the cross-sectional area of the interface.

The results for our model interface system are displayed in Fig. 5. The thermal conductance obtained by NEMD is slightly smaller than the results of the NEGF calculation. This can to a some extent be attributed to finite size effects in the NEMD configuration, which suppresses the contribution of some long-wavelength modes. In fact, some test simulations with very long systems (up to 5 times longer leads on both sides) showed that the NEMD-conductance increases to values slightly larger than the NEGF results for all considered temperatures. However, for the scope of this example, we did not perform an exhaustive convergence study with respect to the system size. Importantly, we note that both methods yield values that are consistent with each other.

9.2 Vapor Deposition

Vapor deposition is another area where MD simulations can be used to understand the underlying microscopic processes and the effect of experimental parameters.

For technical reasons, most MD codes do not allow the number of atoms to change during a single MD simulation. This makes depositions simulations challenging. A user must run an initial MD simulation, then add new atoms or molecules in a script outside the MD-code, prepare and write the updated input files, start a new MD simulation, and repeat until the required timescale is reached. In ATK, however, this can be achieved in a single Python script, by adding an outer loop around the MD function, in which newly deposited atoms or molecules are added to the current configuration at each cycle. Another advantage of this approach, compared to a built-in deposition functionality which is offered e.g. by the LAMMPS package, is that the user can use Python scripting to exactly specify how the deposition is supposed to take place. In this example, we demonstrate a vapor deposition simulation of selenium molecules onto a crystalline substrate using a slightly modified Stillinger-Weber potential [96]. To account for the polydispersity of selenium vapor, the deposited molecules are drawn from a distribution of ring-like selenium molecules of different sizes. In a recent study, we have shown that such simulations can help elucidate the structural differences in vapor-deposited films compared to melt-quenched amorphous selenium [97].

A schematic picture of the simulation setup is shown in Fig. 6. The starting configuration consists of a slab model of the (100) surface of a trigonal selenium crystal. In this substrate configuration, two regions are selected and marked via tags: The bottom layer, which will be fixed in the entire deposition simulation to mimic an infinitely extended perfect crystal substrate, and the substrate region above this fixed layer. In the MD simulation, a thermostat will be applied only to that region to remove excess thermal energy from the substrate, while leaving the ballistic dynamics of the incoming molecules unaltered.

Schematically this looks as shown in the code snippet in the A.2. First, a function (only the signature is shown) is defined, that draws a new molecule from a given distribution of molecules. Then the deposition loop is started. In each loop cycle, a new molecule is drawn and its center of mass is shifted to a random lateral position at a fixed height above the substrate. Thermal velocities, corresponding to the desired vapor temperature, are assigned making sure that the center-of-mass velocity is directed towards the surface. The new molecule is inserted into the last MD configuration and a new MD simulation over 20 ps is started from the updated configuration.

The simulation is run for 200 deposition cycles producing a 25 Å thick amorphous film. Snapshots of the intermediate stages of the deposition simulation are shown in Fig. 7. Upon adsorption at the substrate surface, a fraction of the initially ring-like molecules begin polymerizing to form longer chain-like structures. In contrast to experiments, where a precise determination of the exact molecular topology is nearly impossible, the simulation results provide atomistic insight into the initial structure immediately after deposition. In this case, it is possible to determine the ratio between atoms in ring- or chain-like topologies. The simulation results suggest that the vapor-deposited film exhibits predominantly ring-like structures [97] and that the ratio is temperature dependent. This is in contrast to the results obtained from the simulation of melt-quenched amorphous selenium where the dominating molecular structure is chain-like.

9.3 Creep simulation of a copper polycrystal

Even in microscopic devices, many components are not composed of perfect single crystals, but instead contain small crystalline grains, which form a polycrystalline material. Thus, the properties of such a material are not only dominated by the properties of the corresponding crystal itself but by an interplay with the grain boundaries.

Here, we present the simulation of a creep experiment of a copper polycrystal as an example of how such systems can be simulated in VNL and ATK. The initial configuration is set up in the polycrystal builder within VNL (see Fig. 8)333The Polycrystal builder plugin is part of the SCAITools package which can be installed in VNL using the AddOn-manager. Besides generating the grain centers using random numbers, ATK has an autoseed tool, which can generate grain seeds according to a given grain size distribution, using a genetic algorithm [98]. The generated seed positions can be loaded into the polycrystal builder and the polycrystal builder fills the grains with the given crystal materials.. It is possible to specify the size, the estimated number of grains, and the algorithm by which the grain seeds are determined. In this example, we use a cell volume of $12\times 12\times 12\,{\rm nm}^{3}$. Eight initial seed positions are generated based on a uniform random distribution, and the grains are filled by copper FCC crystals with random orientations. The resulting grains have an average volume of 216 nm3 corresponding to an average diameter of 7.4 nm.

To describe the interactions in the system we use an EAM potential [99]. Figure 9 (a) and (b) show the grain structure of the initial polycrystal and after the creep simulation. The grain boundaries are visualized via the Local Structure analysis in VNL, which assigns each atom a local crystal structure by analyzing the arrangement of the neighboring atoms [100]. Before the actual creep simulation, the configuration is equilibrated at a temperature of 500K and zero pressure. For the creep simulation, a tensile stress of 1.5 GPa is applied in the z-direction, while the lateral directions remain coupled to zero pressure. The creep simulation is run for 800 ps using a time step of 2 fs. The code snippet in the A.3 shows the ATKPython script for this simulation block. Since the configuration contains around 150 000 atoms, saving snapshots of the entire system at a high frequency will produce huge trajectory files. On the other hand, it may be necessary to analyze some observables, such as the strain, at a much higher time resolution. As shown in the example script, this can be achieved via a hook function, which measures and stores the desired quantities, in this case the strain, every 100 MD steps, while regular snapshots of the system are stored only every 25 000 steps.

The strain relaxation as function of the simulation time is shown in Fig. 9 (c). Almost immediately after applying the tensile stress, the system responds by elongating by about 2%. This deformation is primarily elastic. In the following simulation time, the strain increases more slowly, which is characteristic of creep behavior. Fig. 9 (c) shows the corresponding changes in potential energy per atom. The initial, mainly elastic deformation is accompanied by an increase in energy, while the creep relaxation results in a continuous decrease in potential energy. However, due to the small number of grains, the strain does not exhibit a nearly linear increase, as typically expected for polycrystalline materials, but a more step-like behavior. The average creep rate $\dot{\varepsilon}=5\cdot 10^{7}$ s*-1* is depicted by the dotted line in Fig. 9 (c).

10 Outlook

In this paper, we have presented the dynamics and optimization functionality in ATK, including the ATK-ForceField module for empirical potentials. One of the main strengths of the ATK suite is its full integration with the Python scripting language which makes it very flexible and allows for defining complex workflows without the need for recompiling any code. Due to the implementation of the performance-intensive parts in C/C++, this flexibility is obtained without sacrificing computational efficiency and the performance of ATK-ForceField is on par, or in some cases faster, than other state-of-the-art classical MD codes. The ATK is seamlessly integrated with the graphical user interface VNL, which makes it easy to set up and interpret ATKPython scripts. In this respect, the usability of the VNL-ATK suite is unparalleled and it can both lower the barrier for beginners in MD simulations as well as make it easier for experts to develop and implement new MD functionality. A recent example, highlighting the advantage of integrating calculators with different levels of accuracy, shows how ATK can be used to combine MD simulations using ATK-ForceField with advanced electronic transport calculations to efficiently obtain the electronic mobility and conductivity of various bulk systems and nano-wires [17].

ATK and VNL are under active development by QuantumWise. ATK-ForceField is developed in collaboration with Fraunhofer SCAI. A number of features not available in ATK-2017, in particular supporting TremoloX’s native hybrid MPI/OpenMP-parallelization in ATK, are planned for ATK-2018. Furthermore, a toolbox for fitting potential parameters to ab-initio calculations is under development. We particularly expect such workflows for fitting parameters to benefit from synergies present in the ATKPython framework, e.g. calculators with increasing levels of accuracy are readily available to provide reference observables, and a broad range of optimization libraries are natively included in Python modules, such as scipy.

11 Acknowledgements

This work was funded in parts by the German Federal Ministry for Education and Research and the Innovation Fund Denmark under the Eurostars projects E!6935 ATOMMODEL and E!9389 MULTIMODEL.

Appendix A Code Examples

Here, we provide some snippets containing the ATKPython code for the simulation examples in section 9.

A.1 NEMD Simulations of Thermal Transport

The following snippet shows how to run a non-equilibrium MD simulation to obtain the thermal conductance through an impurity layer of germanium in a silicon crystal. The exchange of momenta according to the RNEMD method is invoked as a hook function object in the MD simulation, which is supplied as a pre-defined class in ATKPython. It is initialized by passing the two tag names referring to the groups of atoms in the heat source and sink, as well as the interval at which momenta are exchanged.

⬇

Add tags to atoms in the heat source and heat sink

bulk_configuration.addTags(’heat_source’, heat_source_indices)

bulk_configuration.addTags(’heat_sink’, heat_sink_indices)

Define the hook function which invokes the exchange of momenta.

momentum_exchange_hook = NonEquilibriumMomentumExchange(

configuration=bulk_configuration,

exchange_interval=200,

heat_source=’heat_source’,

heat_sink=’heat_sink’

)

Set up the MD integrator method.

method = NVEVelocityVerlet(

time_step=0.5*femtoSecond,

initial_velocity=ConfigurationVelocities()

)

Run the RNEMD simulation for 1 000 000 steps.

md_trajectory = MolecularDynamics(

bulk_configuration,

constraints=[],

trajectory_filename=’Si_w_Ge_NEMD.hdf5’,

steps=1000000,

log_interval=1000,

post_step_hook=momentum_exchange_hook,

method=method,

)

Print the average thermal flux in the NEMD simulation

nlprint(momentum_exchange_hook)

A.2 Vapor Deposition Simulation of Selenium

The following snippet shows how a vapor-deposition simulation of selenium molecules onto a crystalline substrate can be set up in ATKPython. The function to select a new molecule from a given distribution is only shown schematically and can be implemented as desired for the given application. The central part of the script is a loop over all deposition events. In each iteration a new molecule with thermal velocities directed towards the surface is added to the configuration. This is followed by an MD simulation in which the deposited molecule can adsorb and the configuration can equilibrate for a short time.

⬇

Define helper function to randomly draw a molecule from

a molecular distribution with given probabilities.

def selectFromDistribution():

...

return coordinates

Loop over the deposition events.

for i_deposit in range(deposition_cycles):

# Generate the coordinates for a new vapor molecule.

new_positions = selectFromDistribution()

mol_size = len(new_positions)

# Shift the molecule to a random lateral position at fixed height.

new_positions[:,0] += numpy.random.uniform()*Lx

new_positions[:,1] += numpy.random.uniform()*Ly

new_positions[:,2] += deposition_height

# Set random velocities at vapor temperature.

thermal_vel = (T_vapor*boltzmann_constant/Selenium.atomicMass())**0.5

new_velocities = numpy.random.normal(size=(mol_size, 3))*thermal_vel

com_velocity = new_velocities[:,2].sum()/mol_size

# The z-component of the velocity should point towards the substrate.

if com_velocity > 0.0*Ang/fs:

  new_velocities[:, 2] -= 2.0*com_velocity

new_elements = [Selenium]*mol_size

# The new molecule is added to the current configuration.

bulk_configuration._insertAtoms(

    elements=new_elements,

    positions=new_positions,

    velocities=new_velocities,

)

# Run an MD simulation depositing the atom towards the surface.

md_trajectory = MolecularDynamics(...)

# Extract the final snapshot as starting configuration

# for the next deposition cycle.

bulk_configuration = md_trajectory.lastImage()

A.3 Creep Simulation of a Copper Polycrystal

The following snippet shows how creep simulation of a copper polycrystal can be run using ATKPython. Here, a hook function is used to calculate and store the current strain in the system.

⬇

Define a hook class to store the time-dependent strain.

class StrainHook(object):

def __init__(self, interval):

    self._interval = interval

    self._cells = []

def __call__(self, step, time, configuration, forces, stress):

    """ Apply the hook in the MD simulation. """

    if (step % self._interval == 0):

        cell = configuration.primitiveVectors().inUnitsOf(Ang)

        self._cells.append(cell)

def strain(self):

    """ Return the strain in z-direction. """

    cells = numpy.array(self._cells)

    strain_zz = cells[:, 2, 2]/cells[0, 2, 2] - 1.0

    return strain_zz

Set up the barostat with anisotropic coupling to tensile stress of 1.5 GPa.

The lateral target stress is set to zero.

method = NPTMartynaTobiasKlein(

time_step=2*femtoSecond,

reservoir_temperature=500*Kelvin,

reservoir_pressure=[0, 0, -1.5]*GPa,

initial_velocity=ConfigurationVelocities(),

)

Calculate and store the strain every 100 steps.

strain_hook = StrainHook(interval=100)

Run the creep simulation.

md_trajectory = MolecularDynamics(

bulk_configuration,

constraints=[],

trajectory_filename=’Cu_polycrystal_creep.hdf5’,

steps=400000,

log_interval=25000,

pre_step_hook=strain_hook,

method=method

)

bulk_configuration = md_trajectory.lastImage()

Obtain the time-dependent strain stored in the hook.

strain_zz = strain_hook.strain()