Molecular Dynamics Simulations with Long-Range Interactions

Abstract

The Wolf summation (Wolf et al., J. Chem. Phys. 110, 8254 (1999)), an order O(N) method for the calculation of long-range interactions, has been adapted successfully to the simulation of metal oxides. We present the combination of the method with the Tangney–Scandolo model for polarizable oxygen and the Streitz–Mintmire model for variable charges at metal oxide–metal interfaces. The methods have been implemented successfully in our molecular dynamics package IMD and applied to the structure of several metal oxides. The new methods allow for the simulation of cracks in oxides and the study of the flexoelectricity effect.

Introduction

Computer simulations with long-range interactions are challenging for large systems, since the basic method, Ewald summation [7], is of order O(N ^3∕2) at best. In the Ewald method the charges are smeared out with a Gaussian, thus becoming short-ranged. To get the correct interaction, the smeared out part has to be computed in Fourier space where it is also short-ranged. Adapting free parameters results in the O(N ^3∕2) complexity. The drawback is that the method only works for periodic boundary conditions. Surfaces for example cannot be simulated with the standard method.

Wolf et al. [23] showed that it is possible under certain special conditions to avoid the reciprocal part. The order of the algorithm is reduced O(N) and simulations are possible for more general boundary conditions without a reciprocal space part. The Wolf summation is especially suitable for compensated systems, i.e. structures where the number of positive and negative charges is the equal and the system is neutral as in metal oxides for example. The Wolf summation leads to a performance gain of 3 orders of magnitude in molecular dynamics (MD) simulations.

In the case of oxides it is not sufficient to treat the oxygen atoms as point charges, since they are polarizable. The Tangney–Scandolo (TS) model takes care of this property by assigning an induced dipole to each oxygen atom. In the simulations the dipole strengths are computed self-consistently between each molecular dynamics (MD) simulation step. Since the TS model was developed for the Ewald summation and SiO₂, we first had to re-optimize the interaction model parameters [5]. This was carried out with our force-matching code potfit [4]. New results could then be obtained for magnesia [1]. The importance of the polarizability of the oxygen atoms was demonstrated also in [2]. Since the Wolf summation permits open boundaries, the simulation for cracks in alumina became possible [11]. This lead to the discovery of (anti-)flexoelectricity: α-alumina Al₂O₃ cannot have a permanent polarization even if it is deformed homogeneously since it is inversion symmetric. A spontaneous polarization, however, is observed in front of a crack since the deformation is inhomogeneous. Thus the crack generates a deformation gradient which leads to polarization even in the case of inversion symmetry. Only a short account of this effect is given below since it requires a more detailed study.

The third step in our sequence of methods is the combination of metal oxides and metals. Up to now fixed charges have been preset in the process of modeling the interaction with potfit by fitting the interaction to quantum mechanical ab initio data. Since metals are conducting and metal oxides are insulating, a combination of the two across an interface will lead to charge exchange and interface layers. The atoms in the metal will be neutral, while the atoms in the metal oxide will be charged. Thus we have to introduce variable charges and we must determine their size. Streitz and Mintmire [19] have developed a model (SM) which permits the determination of the charges by minimizing the chemical potential. We have implemented the SM model in IMD and applied it successfully to aluminum-alumina interfaces. New interactions have been determined again by potfit, dynamical simulations, however, have not yet been successful since the interfaces turned out to be unstable in this model.

A further challenge is the determination of the chemical potential. Since this is a global observable, the minimization has to be carried out for the whole sample which requires the application of massively parallel minimization packages. Up to now a simple conjugate gradient algorithm is used which works well but is not easily parallelizable. However, it should be possible to precompile the charges from small samples since the chemical potential is constant in the bulk far from the surface and to apply the calculated charges in the large samples. Since it varies only slowly during simulation it should be possible to adjust it iteratively like we do for the induced dipoles. The Wolf summation and the TS model on the other hand are trivially parallelizable since they fit perfectly in the standard interaction scheme which uses cell decomposition [15, 20].

Description of the Models

In the following a short account of the different models is given as they have been implemented in IMD. The molecular dynamics simulation package IMD itself has been described in [15, 20] and in several HLRS reports [17]. The code can be obtained from the webpage imd.itap.physik.uni-stuttgart.de, together with a detailed userguide.

The Wolf Summation for Long-Range Interactions

Ewald summation [7] is the method of choice for long-range interactions, but it is limited to small systems. For larger systems lattice- and tree-based methods exist which can lead to a complexity of nearly O(N). The drawback is that they do not fit very well into the simulation code IMD and are difficult to parallelized. For small systems they are typically slower than the Ewald summation due to a large pre-factor in the complexity. Methods which involve Fourier transform are limited to periodic boundary conditions and cannot be applied to crack simulations for example.

The Wolf summation [23] is especially suitable for ionic melts and solids since these systems are compensated: the ions screen their opposite charges across rather short distances already. Wolf et al. have evaluated the Ewald method and found conditions under which the sum in reciprocal space and thus the Fourier transform part is negligible. The drawback is a rather large cut-off radius which is about twice as large as without electrostatic contributions. Neighbor lists can be quite big, but on today’s machines that is not a problem. The gain of the method is a short-ranged potential which can be simulated together with all the other short-range interactions. It therefore leads to the full flexibility in the choice of boundary conditions.

The Tangney–Scandolo Model for Polarizable Oxygen

Tangney and Scandolo [21] devised a method for the simulation of silica SiO₂. The non-electrostatic parts of the interaction are modeled by a parameterized Morse-stretch potential. The ions get a fixed charge while the oxygen atoms possess an additional dipole moment. The polarizability of the oxygen ions is fixed and their dipole moment is calculated from the local field intensity. The field is computed in an iterative way from the monopoles and dipoles. The electrical field is determined by the ionic charges and dipole moments of the previous time step, in the original TS model by Ewald summation. The starting value of the iteration is predicted from the previous time steps to increase the efficiency.

After extending IMD for the simulation of monopoles with the Wolf summation we split up the generally non-convergent sums for energies and forces of the dipoles into a direct and a reciprocal space part. We could show that this second part can be neglected as it was the case for the monopoles. The real space part has to be cut off smoothly such that the forces vary continuously. For the conservation of energy it was found that the cutoff for monopoles and dipoles has to be chosen consistently. The final choice of the cutoff radius is a compromise between precision and efficiency.

The original TS calculations could be reproduced exactly with our implementation of the Ewald summation. With the Wolf summation the simulation time could be boosted by three orders of magnitude. The next step was to optimize the potential parameters for the different behavior of the Wolf summation. This was carried out with our potfit code and was applied to silica, alumina and magnesia MgO [1].

Details of the Implementation in IMD of the Wolf Summation Combined with the TS Model

Smooth Cut-Off

For MD simulations with limited-range interactions the potentials and their first derivatives must go to zero continuously at a cut-off radius r _c ; otherwise, atoms crossing this threshold might get unphysical kicks. For the Morse-stretch pair potential this is generally not problematic, as it decays with r _ij fast enough. Following Wolf [23] the potential $U_{\text{MS}}(r_{\mathit{ij}})$ is replaced by

$$\displaystyle{ \tilde{U}_{\text{MS}}(r_{\mathit{ij}}) = U_{\text{MS}}(r_{\mathit{ij}}) - U_{\text{MS}}(r_{c}) - (r_{\mathit{ij}} - r_{c})U^{\prime}_{\text{MS}}(r_{c}), }$$

(1)

where a prime denotes a derivative with respect to r.

The other functions used in the TS model have a general r dependence of the form r ⁻ⁿ, n ∈ { 1, 2, 3}. Especially the Coulomb energy with its r ⁻¹ dependency cannot simply be cut off without a treatment as in (1), for otherwise the energy of the system would fluctuate strongly with r _c , without convergence to the proper value. But even with a smooth cut-off (1), for which the Coulomb energy converges, a rather large cut-off radius would be required to make shifting of the potential negligible. Fortunately, the Wolf direct summation method [23] includes a weak exponential damping of the Coulomb potential by $\mathrm{erfc}\!\left (\kappa r\right )$. Such a damped potential can be cut off smoothly at a much smaller radius r _c without affecting the result. All integer powers of r ⁻¹ are treated in a way to conserve the differential relationship between the functions, i.e. the damped functions are

$$\displaystyle\begin{array}{rcl}{ r}^{-1} \rightarrow & {r}^{-1}\mathrm{erfc}(\kappa r) =: f_{ -1}(r),&{}\end{array}$$

(2)

$$\displaystyle\begin{array}{rcl}{ r}^{-2} = -\frac{d({r}^{-1})} {\mathit{dr}} \rightarrow & -& \frac{d({r}^{-1}\mathrm{erfc}(\kappa r))} {\mathit{dr}} \\ & =& {r}^{-2}\mathrm{erfc}(\kappa r) -\frac{2\kappa \exp ({-\kappa }^{2}{r}^{2})} {\sqrt{\pi }r} =: f_{-2}(r).{}\end{array}$$

(3)

This procedure is also required to conserve the energy during an MD simulation, as discussed in more detail in Sect. 2.3.

The damped potentials are then shifted to zero and zero derivative at the cut-off radius, as in (1). This allows for limited-range MD simulations with a standard MD code. The computational effort of such a simulation scales linearly in the number of particles (as the number of interactions that need to be evaluated per particle does not increase with the number of particles), but scales roughly with $O(r_{c}^{3})$.

Energy Conservation

In MD simulations the energy is conserved, if the forces on the particles are exactly equal to the gradient of the potential energy with respect to the atomic coordinates. Otherwise, the energy might oscillate or even drift off if not controlled by a thermostat. In standard MD simulations, the requirement is usually automatically fulfilled: The forces are calculated as the derivative of the potential, which depends directly on the atomic positions. In the TS model, there is also an indirect dependence, as the potential is also a function of the dipole moments:

$$\displaystyle{ U = U(\{\boldsymbol{r}_{i}\},\{\boldsymbol{p}_{i}(\{\boldsymbol{r}_{j}\})\}). }$$

(4)

This would in principle lead to an extra contribution to the derivative of the potential,

$$\displaystyle{ \frac{dU} {d\{\boldsymbol{r}_{i}\}} = \frac{\partial U} {\partial \{\boldsymbol{r}_{i}\}} + \frac{\partial U} {\partial \{\boldsymbol{p}_{i}\}} \frac{\partial \{\boldsymbol{p}_{i}\}} {\partial \{\boldsymbol{r}_{j}\}}, }$$

(5)

which would be practically impossible to determined effectively. Luckily, if the dipole moments are iterated until convergence is reached, we are at an extremal value of the potential energy, with $\partial U/\partial \{\boldsymbol{p}_{i}\} = 0$, and so this part need not be evaluated. Imperfections in convergence may lead to a drift in the energy, however, as was already observed by Tangney and Scandolo [21].

When applying the Wolf formalism to the TS potential, another issue arises concerning the conservation of energy. It can most easily be explained with a simple one-dimensional example. Given are two oppositely charged point charges ± q at a mutual distance r. If the negatively charged one is polarizable with polarizability α, it will get a dipole moment $p =\alpha q/(k{r}^{2})$, with k = 4π ε ₀. This leads to a total interaction energy

$$\displaystyle\begin{array}{rcl} U& =& \underbrace{\mathop{-2 \cdot \frac{1} {2} \frac{1} {k} \frac{{q}^{2}} {r} }}\limits _{q-q} -\underbrace{\mathop{ 2 \cdot \frac{1} {2} \frac{q} {k} \frac{p} {{r}^{2}} }}\limits _{q-p} +\underbrace{\mathop{ \frac{1} {2} \frac{{p}^{2}} {\alpha } }}\limits _{\text{dipole}},{}\end{array}$$

(6)

from which it follows that

$$\displaystyle\begin{array}{rcl} \frac{\partial U} {\partial p} & =& -\frac{1} {k} \frac{q} {{r}^{2}} +\underbrace{\mathop{ \frac{p} {\alpha } }}\limits _{= \frac{1} {k} \frac{q} {{r}^{2}} } = 0.{}\end{array}$$

(7)

Here, q − q denotes the Coulomb interaction between charges, q − p the interactions between charge and dipole, and the last term is the dipole energy. When we now damp and cut off the interactions, we replace the ${r}^{-1},{r}^{-2}$ functions by their damped and smoothed counterparts $\tilde{f}_{-1}(r),\tilde{f}_{-2}(r)$. If energy conservation is to be maintained, the differential relation between the $\tilde{f}_{-n}$ must be the same as for the r ⁻ⁿ:

$$\displaystyle\begin{array}{rcl} \frac{d\tilde{f}_{-1}(r)} {\mathit{dr}} = -\tilde{f}_{-2}(r).& &{}\end{array}$$

(8)

As a consequence, the first two derivatives of the smoothed damped Coulomb potential must be zero at r _c .

In MD simulation it is computationally advantageous to represent pair potential functions internally as functions of r ², and their derivative as ${f\hat{}^{{\prime}}}:= {r}^{-1}\mathit{df }/\mathit{dr}$. In this way, Wolf summation can be applied to dipolar interactions in the TS potential model.

The Streitz–Mintmire Model for Variable Charges

At interfaces of metals and metal oxides charges have to be variable at least in an transition region. The charge of the metal atoms will depend on the local oxygen concentration and vary from neutral to the maximal charge in the oxide. There are nearly no models for large scale simulations. It is not clear if mirror charges or interface diffusion has to be taken into account. In the model of Streitz and Mintmire [19] long-range interactions are taken into account adequately. The total energy is composed of an embedded atom potential and an electrostatic contribution which contains the charges to second order. The charge of each atom is obtained by minimizing the electrostatic energy under the condition of charge neutrality. The optimized charges minimized a quadratic expression which leads to a linear system of equations. Streitz and Mintmire [19] have inverted the matrices directly which is not suitable for large systems and massively parallel computations. We have implemented a conjugate gradient method which solves the equations iteratively. The efficiency can be increased by taken the values of the previous step as starting conditions since the optimized charges change only little from time step to time step.

Again we apply the Wolf summation at all steps: for the determination and minimization of the charge matrix as well as for the calculation of the long-range energies and forces. Figure 1 shows the variable charges across a alumina aluminum surface.

Fig. 1

Variable charge across an alumina-aluminum surface computed with the SM model, the TS model and Wolf summations. Abscissa: position in \r{A}, ordinate: charge in elemental charge units

A weak point of the SM model is the necessity to minimize the global chemical potential which involves global communication. An improvement has been studied by Elsener et al. [6] for the case of single oxygen ions in a metal matrix. There the global chemical potential can be replaced by a local chemical potential. A similar approach could be advisable for interfaces, since in the bulk far from the interface the charges should attain the values of the pure oxide or metal. Consequently, the proper values can be precompiled from smaller samples and the actual values will change only little during simulation. Thus the can be adjusted iteratively.

Results

A short overview is given of the validation of the new models in simulations with IMD for a number of structures by studying structural and thermodynamic properties. Crack simulations and first results for (anti-)flexoelectricity will also be presented.

Simulation of Oxides

Validation of the Wolf Summation

A comparisons of the Ewald and Wolf summation of silica shows that the differences are marginal. As a first step we have assumed the original TS parameters for SiO₂. Figures 2–4 clearly show a good agreement with ab-initio data and the original TS model with Ewald summation. For comparison simulations with the pair potential of Beest et al. [22] are given which is still widely used.

Fig. 2

Equation of state of liquid SiO₂ at 3,100 K for the damped and smoothly cut TS potential compared to the original TS potential [21], ab-initio calculations, simulations with the BKS potential and experiment [8]

Fig. 3

Si–O–Si angular distribution in liquid SiO₂ at 3,100 K for the new potential in comparison to the original TS potential [21], ab-initio calculations, and simulations with the BKS potential

The original potential has also been applied to other structure variants of SiO₂, namely quartz, coesite, and cristobalite [5], where it also lead to good agreement with experimental data and the results of Tangney and Scandolo. As a rule of thumb, the real space cut-off radius required for the Wolf summation should be about five times the largest nearest neighbor distance of opposite charges in the system. Still we obtain a speedup of more than two orders of magnitude compared to the original TS model with Ewald summation.

Fig. 4

Percentage of N-valued Si atoms in liquid SiO₂ at 3,100 K as a function of the pressure for the new potential in comparison to the original TS potential [21] and simulations with the BKS potential

Silica and Magnesia

After validating the TS model with Wolf summation for the original interactions new interactions have been fitted with potfit for silica and magnesia [1]. For both structures the parameters of the Morse-stretch potential, the polarizabilities and the charges were obtained. The potential was validated by determining thermodynamic, microstructural and vibrational properties of silica. They are in excellent agreement with ab initio data. As a test, the new parameters were applied to amorphous silica which was not used for the fitting procedure. The vibrational density of state could be reproduced very well. A further test was the application to low-pressure α-quartz which yields a hard case for the transferability of the interaction. We get reasonable results although the interaction was not optimized for such crystal structures.

A further test was the application of the TS model to magnesia MgO. Details and interaction parameters can be found in [1]. The interaction parameters were fitted with potfit to a database of small reference structures for which the forces, stresses, and energies were computed with VASP [13, 14]. The new interactions are validated in the same way for liquid magnesia. Again an excellent agreement was found if compared to a study with ab initio data [12]. Although the potential was fitted only up to 15 GPa, it yielded good values up to 160 GPa. The vibrational density of state is also in good agreement with experiment [3] and ab initio calculations [9].

Not surprisingly the scaling with system size is linear for both materials, silica and magnesia.

The Necessity of Polarizability

The question is how important the polarizability is for the simulation of the oxides. A systematic comparison has been carried out for silica, magnesia and alumina [2] for two types of force fields which only differ in the polarizable term; the non-electrostatic and Coulomb term have the same form for both force fields. The influence on the radial distribution functions is negligible, the influence on the bond angle distribution function is stronger. The interactions without polarizability overestimate the lower angles and underestimate higher angles. Interactions with polarizability yield the correct behavior for increasing temperatures and are in better agreement with ab-initio data. For alumina the cohesive energy was also calculated. It was found that it is strongly overestimated without polarizability. The equation of state for silica and magnesia is underestimated without polarizability as compared to ab initio calculations. In summary we have clearly shown that the interactions with polarizable oxygen atoms yield much better results than without polarization.

Crack Simulations of Alumina

Interactions for α-alumina Al₂O₃ have been determined in pretty much the same way as for silica and magnesia described in Sect. 3.1. Special care was taken for a correct behavior of surface relaxation and surface energies since crack simulations were intended. A good agreement with ab initio calculations and experimental data was obtained. The new interaction underestimates the energy of different twins, but qualitative agreement is achieved. Details about the crack studies can be found in [11]. Cracks do not move in close packed planes. Cracks which are not placed not in the cleavage plane tend to be deflected into The most important result is the observation of spontaneous polarization of the sample in front of the crack. This is especially remarkable since α-alumina is inversion symmetrical which forbids polarization. Since effects from the boundary conditions can be excluded by applying periodic boundary conditions the only remaining influence is an inhomogeneous deformation which leads to (anti-)flexoelectricity as described in the next section.

Flexoelectricity

If an ionic crystalline structure without inversion symmetry is deformed it will show piezoelectricity since the deformation tensor $\epsilon _{\mathit{jk}}$ will couple to the polarization vector P _i via a third order tensor d _ijk :

$$\displaystyle{ P_{i} = d_{\mathit{ijk}}\epsilon _{\mathit{jk}} }$$

(9)

This leads to a free energy f _piezo equal to

$$\displaystyle{ f_{\mathrm{piezo}} = \frac{1} {2}E_{i}P_{i} = \frac{1} {2}E_{i}d_{\mathit{ijk}}\epsilon _{\mathit{jk}} }$$

(10)

where summation over equal indices is assumed. Since E _i is an polar vector while ε _jk and d _ijk are inversion invariant, the energy must vanish for polar structures since f _piezo would change its sign. This is the case for periclase and α-alumina.

If the deformation is inhomogeneous, however, it could couple to P _i through its gradient ∂ _i n via a fourth order tensor μ _ijkl :

$$\displaystyle{ P_{i} =\mu _{\mathit{ijkl}}\partial _{j}\epsilon _{\mathit{kl}} }$$

(11)

This is called flexoelectricity. The free energy f _flexo is equal to

$$\displaystyle{ f_{\mathrm{flexo}} = \frac{1} {2}E_{i}P_{i} = \frac{1} {2}E_{i}\mu _{\mathit{ijk}}\partial _{j}\epsilon _{\mathit{kl}} }$$

(12)

Since the gradient ∂ _j is also a polar vector, f _flexo will not change its sign on inversion and therefore the flexoelectric effect will be present for all crystal structures. In the case of polar crystals it will be a second order correction to piezoelectricity.

The coupling of the inhomogeneous deformation to polarization was the reason for the spontaneous polarization of the sample in the crack simulations on α-alumina. Since the crack geometry is difficult to describe and analyse we invented a geometry where flexoelectricity is more easy to study: A rectangular plate was deformed in such a way that it formed a trapezoidal wedge. With this geometry two of the three independent constants of μ _ijkl in the cubic case could be determined for periclase MgO. The constants are μ ₁₁=2× 10⁻¹³ C/cm and $\mu _{12} = -12\times$ 10⁻¹³ C/cm.

If two wedges with opposite deformation are fit together, a Néel wall is formed as shown in Fig. 5.

Fig. 5

Two cubic MgO periclase blocks deformed as indicated by the boundaries. The polarization in the green part points to the upper right corner, the polarization in the lower part points to the lower right corner. The different polarizations are adjusted in the center by an Néel wall. Visualization using MegaMol [10, 18]

Things become much more complicated for α-alumina. First of all, ten independent tensor components μ _ijkl are expected since this structure is hexagonal. Second, as can bee seen in Fig. 6, the material becomes antiflexoelectric! The domains are not as clearly separated as for periclase. Up to now we have not fully understood what is happening.

Fig. 6

Two hexagonal α-alumina Al₂O₃ blocks deformed as indicated by the boundaries. The arrows indicate the local polarization. As nearby arrows point in opposite directions, we find that the sample splits into antiflexoelectric domains. Visualization using MegaMol [10, 18]

In summary we have found that it is possible to simulate flexo- and antiflexoelectricity with classical MD simulations if the potential allows for induced dipoles. The investigations of the observed structures and the determination of materials parameters through adequate deformation geometries are still at the very beginning.

Performance and Benchmarks

General benchmark data for IMD have been given by Stadler et al. [20]. The data demonstrate that IMD scales almost linearly in weak scaling (same number of atoms per processor) and fairly well for strong scaling (total number of atoms constant, thus communication load growing). This behavior is still valid as a systematic study on the Blue Gene/L clearly shows (see the previous HLRS report [17]). New data especially for the simulation of laser ablation have been reported last year (see [16]). Further benchmarks especially for long-range interactions depend on the algorithms only and are therefore given in the original publications (See e.g. [5]). Results for production simulations on the HLRS Hermite are not yet available since our project did not include access to this machine up to now. Other benchmark data for XE6, XK6 and XT5 machines obtained by Cray have been reported last year [16].

Summary

We have reported the extension of our molecular dynamics simulation package IMD to the effective simulation of ionic materials with long-range interactions including a model of polarizable oxygen. This part is fully parallelized with MPI. IMD has been extended further to variable charges and to the computation of these charges at metal oxide–metal interfaces by the minimization of the chemical potential. The latter part still has to be parallelized. As indicated there are good reasons that the chemical potential can be precalculated and iteratively adjusted during simulations. The new program capacities have been applied successfully to the study of several crystal structures of silica SiO₂, magnesia MgO and alumina Al₂O₃. The study of cracks in crystalline α-alumina lead to the discovery of flexoelectricity in cubic periclase MgO and antiflexoelectricity in hexagonal α-alumina. The new effects are still under investigation and final results can not yet been given.