# The Environment-Dependent Interatomic Potential (EDIP) This directory contains the environment-dependent interatomic potential (EDIP) model driver [[1,2,3,4,5](#references)]. ## EDIP model driver This model driver provides a C++ implementation of the Environment Dependent Interatomic Potential (EDIP) [[1,2,3](#references)] and is based on the source code of the ['edip'](https://docs.lammps.org/pair_edip.html) and ['edip/multi'](https://docs.lammps.org/pair_edip.html) pair styles found in the LAMMPS software package [[4,5](#references)]. It additionally features an `edip/c` mode not found in any pair styles in LAMMPS, which is an extension of EDIP to carbon that accounts for pi-bonding effects [[2,3](#references)]. The three modes of the EDIP driver are, - ['edip'](#edip-style) - ['edip/multi'](#edipmulti-style) - ['edip/c'](#edipc-style) The mode of the driver (`edip` for single-species models, `edip/multi` for multispecies models, or `edip/c` for carbon models with pi-bonding effects) is automatically detected based on the input files provided, which must be ASCII text files. For the `edip` and `edip/multi` modes of operation, the driver expects an [**element**](#edip-element-file-format) file and a [**parameter**](#edip-parameter-file-format) file. The element file contains a unique list of chemical elements. In the case of a single species, it automatically chooses the `edip` mode, whereas if multiple species are present, it automatically selects the `edip/multi` mode. However, this driver only expects a [**parameter**](#edipc-parameter-file-format) file for the `edip/c` mode, since only pure carbon is supported in that case. **NOTE:** We adapted the LAMMPS format from the [LAMMPS documentation](https://docs.lammps.org) as of `July 30, 2021`. ## EDIP style The current model is an empirical potential popular for modeling silicon materials [[1](#references)]. It consists of two and three-body interactions with theoretically motivated functional forms that capture chemical and physical trends. The two and three-body terms have environment dependence controlled by the atomic coordination Z. The total energy is written as a sum of on-site energies, ```tex U_i = \sum_j U_{2}(r_{ij}, Z_{i}) + \sum_{j < k} U_{3}(r_{ij}, r_{ik}, Z_{i}), ``` The pair potential `U2` is an interaction between atoms `i` and `j` representing pairwise bonds and `U3` is the interaction between atoms `i`, `j`, and `k` centered at atom `i` representing angular forces. ```tex \begin{align*} &U_{2}(r, Z) = A\left[\left(\frac{B}{r}\right)^{\rho} - e^{-\beta Z^2}\right]exp{\left(\frac{\sigma}{r-a}\right)} \\ &U_{3}(r_{ij}, r_{ik}, Z_i) = exp{\left(\frac{\gamma}{r_{ij}-a}\right)}exp{\left(\frac{\gamma}{r_{ik}-a}\right)}h(cos\theta_{ijk},Z_i). \end{align*} ``` The two-body term, `U2` includes repulsive and attractive interactions, which go to zero at a distance equals `a`. The three-body term, `U3` contains radial and angular factors. Both types of interactions depend on the local environment of atom `i` through its effective coordination number, defined by, ```tex Z_i = \sum_{j \ne i} f(R_{ij}), ``` where `f` is a cutoff function that measures the contribution of neighbor `j` to the coordination of atom `i` in terms of the separation `r`. The neighbor function is unity for `r` less than `c`, with a gentle drop from `1` to `0` between `c` and `a`, and it is `0` for `r` greater than `a`. ```tex f(r) = \begin{cases} 1 & \quad ra \end{cases} ``` All neighbors within the distance `c` to the particle `i` have full contribution, while others between `c` and `a` have a partial contribution to the coordination number `Z`. The angular function, ```tex h(l,Z)=\lambda [(1-e^{-Q(Z)(l+\tau(Z))^2}) + \eta Q(Z)(l+\tau(Z))^2], ``` is strongly dependent on the local coordination through two functions, ```tex \tau(Z)=u_1 + u_2 (u_3 e^{-u_4 Z} - e^{-2u_4 Z}), ``` and ```tex Q(Z)=Q_0 e^{-\mu Z}. ``` These two functions control the equilibrium angle and the strength of the interaction, respectively. The various parameters in the EDIP formulas are listed in a file. The EDIP parameter file. The name usually contains the element name and ends in the `.edip` extension. For example, for a system containing `Si`, the `Si.edip` or `Si.parameter.edip` file contains the specific parameter settings. ### EDIP element file format The element file contains a unique list of chemical element names separated by spaces: ```bash Elem1 Elem2 Elem3 ... ``` For example, for a Si system, the element file may look like this: ```bash Si ``` For a Si and C alloy system, the element file may look like this: ```bash Si C ``` For an alloy system containing Al, Si, Mg, Cu, and Fe, it could look like this: ```bash Al Si Mg Cu Fe ``` ### EDIP parameter file format The EDIP parameter file defines parameters for a triplet of elements. The parameters in a single entry corresponding to the two-body and three-body coefficients in the formula above, ```txt element1 element2 element3 A B cutoffA cutoffC alpha beta eta gamma lambda mu rho sigma Q0 u1 u2 u3 u4 ``` `element1` is the center atom in a 3-body interaction. The `A`, `B`, `beta`, and `sigma` parameters are used only for two-body interactions. The `eta`, `gamma`, `lambda`, `mu`, `Q0` and all `u1` to `u4` parameters are used for three-body interactions. The `alpha` and `cutoffC` parameters are used for the coordination environment function. `A` and `lambda` have energy units, while `B`, `cutoffA`, `cutoffC`, `gamma`, and `sigma` have distance units, and the rest of the parameters are pure numbers [[1](#references)]. ### EDIP style example As an example, a KIM Portable Model (PM) for pure Si systems that uses this model driver with the edip style contains the following two input files, ```bash elems.edip Si.edip ``` where `elems.edip` element file contains, ```bash # list of elements # Elem1 Si ``` and the `Si.edip` parameter file contains, ```bash # DATE: 2011-09-15 CONTRIBUTOR: Unknown CITATION: Justo, Bazant, Kaxiras, Bulatov and Yip, Phys Rev B, 58, 2539 (1998) # EDIP parameters for various elements and mixtures # multiple entries can be added to this file, LAMMPS reads the ones it needs # these entries are in LAMMPS "metal" units # format of a single entry (one or more lines) # # element 1, element 2, element 3, # A B cutoffA cutoffC alpha beta eta # gamma lambda mu rho sigma Q0 # u1 u2 u3 u4 # # units for each parameter: # A , lambda are in eV # B, cutoffA, cutoffC, gamma, sigma are in Angstrom # alpha, beta, eta, mu, rho, Q0, u1-u4 are pure numbers # Here are the original parameters in metal units, for silicon from: # J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov, S. Yip # Phys. Rev. B 58, 2539 (1998) # Si Si Si 7.9821730 1.5075463 3.1213820 2.5609104 3.1083847 0.0070975 0.2523244 1.1247945 1.4533108 0.6966326 1.2085196 0.5774108 312.1341346 -0.165799 32.557 0.286198 0.66 ``` ## EDIP/MULTI style This style supports multi-element systems, where the element file contains a unique list of elements. In addition, its parameter file must contain the number of unique elements to the power of 3 entries. For example, for a two-element system, the file must contain 8 entries and for a three-element system, the file must have 27 entries. ### EDIP/MULTI style example For example, a KIM Portable Model (PM) for modeling silicon carbide is a two-element system. It contains the following two input files, ```bash elems.edip SiC.edip ``` where `elems.edip` element file contains, ```bash # list of elements # Elem1 Elem2 Si C ``` and the `SiC.edip` parameter file contains, ```bash # DATE: 2017-05-16 CONTRIBUTOR: Laurent Pizzagalli CITATION: G. Lucas, M. Bertolus, and L. Pizzagalli, J. Phys. : Condens. Matter 22, 035802 (2010) # # element 1, element 2, element 3, # A B cutoffA cutoffC alpha beta eta # gamma lambda mu rho sigma Q0 # u1 u2 u3 u4 # # units for each parameter: # A , lambda are in eV # B, cutoffA, cutoffC, gamma, sigma are in Angstrom # alpha, beta, eta, mu, rho, Q0, u1-u4 are pure numbers # Si Si Si 5.488043 1.446435 2.941586 2.540193 3.066580 0.008593 0.589390 1.135256 2.417497 0.629131 1.343679 0.298443 208.924548 -0.165799 32.557 0.286198 0.66 C C C 10.222599 0.959814 2.212263 1.741598 1.962090 0.025661 0.275605 1.084183 3.633621 0.594236 2.827634 0.536561 289.305617 -0.165799 32.557 0.286198 0.66 C Si Si 7.535967 1.177019 2.534972 1.973974 2.507738 0.015347 0.432497 1.191567 3.025559 0.611684 2.061835 0.423863 249.115082 -0.165799 32.557000 0.286198 0.660000 Si C C 7.535967 1.177019 2.534972 1.973974 2.507738 0.015347 0.432497 1.191567 3.025559 0.611684 2.061835 0.423863 249.115082 -0.165799 32.557000 0.286198 0.660000 Si Si C 5.488043 1.446435 2.941586 2.540193 3.066580 0.008593 0.510944 1.135256 2.721528 0.620407 1.343679 0.298443 229.019815 -0.165799 32.557000 0.286198 0.660000 Si C Si 7.535967 1.177019 2.534972 1.973974 2.507738 0.015347 0.510944 1.191567 2.721528 0.620407 2.061835 0.423863 229.019815 -0.165799 32.557000 0.286198 0.660000 C C Si 10.222599 0.959814 2.212263 1.741598 1.962090 0.025661 0.354051 1.084183 3.329590 0.602960 2.827634 0.536561 269.210350 -0.165799 32.557000 0.286198 0.660000 C Si C 7.535967 1.177019 2.534972 1.973974 2.507738 0.015347 0.354051 1.191567 3.329590 0.602960 2.061835 0.423863 269.210350 -0.165799 32.557000 0.286198 0.660000 ``` The file must contain 8 entries for `SiSiSi`, `SiSiC`, `SiCSi`, `SiCC`, `CSiSi`, `CSiC`, `CCSi`, `CCC`, that specify EDIP parameters for all permutations of the two elements interacting in three-body configurations. ## EDIP/C style The `edip/c`[[2](#references)] style is a transferable empirical potential developed for carbon by extending the environment-dependent interaction potential proposed for silicon [[1](#references)]. This implementation extends EDIP [[1](#references)] to describe carbon, and in doing so, addresses the most significant weakness of silicon EDIP, namely, the absence of -bonding effects. With this improvement important phenomena such as dihedral rotation penalties and -repulsion are described. The functional form of carbon EDIP consists of three components: - Two-body pair energy, - Three-body angular penalty, - Generalized coordination. Following silicon EDIP [[1](#references)], the two and three-body terms have environment dependence controlled by the atomic coordination Z. Within this framework, the total energy is written as a sum of on-site energies. ```tex U_i=\sum_j U_{2}(r_{ij}, Z_{i})+\sum_{j ```tex U_{2}(r,Z)=\epsilon\left[\left(\frac{B}{r}\right)^4-e^{-\beta Z^2}\right]\exp{\left(\frac{\sigma}{r-a-a'Z}\right)}. ``` The describes the bond-order. The three-body term follows silicon EDIP [[1](#references)] in using SW-like radial and angular terms with environmental dependence, having the form, ```tex U_{3}(r_{ij},r_{ik},\theta,Z)=\lambda(Z)g(r_{ij},Z)g(r_{ik},Z)h(\theta,Z). ``` The components of `U3` are defined by, ```tex \begin{align*} &\lambda(Z)=\lambda_0\exp[-\lambda'(Z-Z_0)^2],\\ &g(r,Z)=\exp[\gamma/(r-a-a'Z)],\\ &h(\theta,Z)=1-\exp\left\{-q[\cos\theta+\tau(Z)]^2\right\}. \end{align*} ``` describes the variation in ideal bonding angle and is defined by,
```tex \tau(Z)=1-\frac{Z}{12}\tanh[t_1(Z-t_2)]. ``` The spherical coordination contribution is determined by the sum , where is a three-parameter function defined by, ```tex f(r) = \begin{cases} 1 & \quad rf_{high} \end{cases} ``` The generalized coordination augments , ```tex Z_i=z_i+\pi_3(z_i)X_i^{dih}+\pi_3(z_i)X_i^{rep3}+\pi_2(z_i)X_i^{rep2}. ``` and are switching functions identifying 2-coordinated and 3-coordinated atoms, and describes dihedral rotation, -repulsion at a threefold site, and -repulsion at a twofold site. Vector products capture the appropriate symmetries via the functions, ```tex \begin{align*} &X_i^{dih}=Z_{dih}\sum \pi_3(z_j)(\mathbf{\hat{R}_{jm}.\hat{R}_{ik}\times\hat{R}_{il}})^2C^{dih}_{ijklm},\\ &X_i^{rep3}=Z_{rep}\sum \pi(z_j)(\mathbf{\hat{R}_{ij}.\hat{R}_{ik}\times\hat{R}_{il}})^2C^{rep3}_{ijkl},\\ &X_i^{rep2}=Z_{rep}\sum \pi(z_j)[1-(\mathbf{\hat{R}_{ij}.\hat{R}_{ik}})^2]C^{rep2}_{ijk}. \end{align*} ``` where `jkl` are neighbors of `i`, and `m` is a neighbor of `j`, and selects for twofold or threesfold sites. The distance-based cutoff functions `C` are, ```tex \begin{align*} &C^{dih}_{ijklm}=p(R_{ij})p(R_{ik})p(R_{il})p(R_{jm}),\\ &C^{rep3}_{ijkl}=(R_{ij}-c_0)^2[1-p(R_{ij})]p(R_{ik})p(R_{il}),\\ &C^{rep2}_{ijk}=(R_{ij}-c_0)^2[1-p(R_{ij})]p(R_{ik}). \end{align*} ``` with the function `p(r)` equivalent to `f(r)` except for different end points. The various parameters in the EDIP/C formulas are listed in a parameter file. Usually, it is called the `C.edip` parameter file. **NOTE:** To accurately model close-ranged interactions between atoms, the pair potential within the EDIP for carbon (EDIP/C) formalism is smoothly switched to the Ziegler-Biersacke-Littermark (ZBL) pair potential at small separations. However, due to the environmental dependence, this transition is less straightforward than pair potentials, where one can use interpolation functions to connect the potentials. Here [[3](#references)], the approach is to use two Fermi-type scaling functions, ```tex \begin{align*} f_{zbl}(r)=\left[1+\exp{\left(\frac{r-r_c+\delta}{s}\right)}\right]^{-1},\\ f_{edip}(r)=1-\left[1+\exp{\left(\frac{r-r_c-\delta}{s}\right)}\right]^{-1}, \end{align*} ``` where is the cutoff distance, is the shift and `s` is the skin (controls the sharpness of the transition). ### EDIP/C parameter file format The `EDIP/C` parameter file defines parameters coefficients, ```txt eps BB beta sigma a aprime Z0 lambda0 lambdaprime gamma q flow fhigh alpha plow phigh Zdih Zrep c0 zbl cutoff skin shift ``` The `eps`, `BB`, `beta`, `sigma`, `a`, and `aprime` parameters are used for two-body interactions. The `Z0`, `lambda0`, `lambdaprime`, `gamma`, and `q` parameters are used for three-body interactions. The `flow`, `fhigh`, `alpha`, `plow`, `phigh`, `Zdih`, and `Zrep` parameters are used for the coordination environment function. The `zbl` parameter indicates whether the pair potential within the EDIP formalism is smoothly switched to the ZBL pair potential at small separations. Due to environmental dependence, this transition is not straightforward. The approach is to use two Fermi-type scaling functions which are defined using three parameters `cutoff`, `skin`, and `shift`. `eps`, `lambda0`, `Zdih`, `Zrep` parameters have `energy units`, and `sigma`, `a`, `aprime`, `gamma`, `flow`, `fhigh`, `plow`, `phigh`, `c0`, `cutoff`, `skin`, and `shift` have `distance units`. The `BB` parameter has the power of four `distance unit`. And the rest of the parameters are pure numbers. **NOTE:** Setting the `zbl` parameter to `1` indicates using ZBL pair potential and setting it to `0` indicates not using it. If `zbl` parameter is set to `0` any other parameter after that is optional and is ignored. ### EDIP/C style example As an example, a KIM Portable Model (PM) for modeling carbon contains the following parameter file, ```bash C.edip ``` where the `C.edip` parameter file contains, ```bash # DATE: 2021-06-06 Reference: Marks, N.A., Phys. Rev.B, 63(3):035401, Dec 2000. # format of a single entry (one or more lines) # # eps BB beta sigma a aprime # Z0 lambda0 lambdaprime gamma q # flow fhigh alpha plow phigh Zdih Zrep c0 # zbl cutoff skin shift # # if zbl is 0, then \c cutoff, \c skin, \c shift are not necessary. # # units for each parameters: # # eps, lambda0, Zdih, Zrep -> are in eV # sigma, a, aprime, gamma, flow, fhigh, plow, phigh, c0 -> are in Angstrom # cutoff, skin, shift -> are in Angstrom # BB -> is in Angstrom^4 # beta, Z0, lambdaprime, q, alpha, -> are pure numbers # zbl -> is 0 or 1 # 20.0853862863495 0.827599951322299 0.0490161172713279 1.25714643580808 1.89225338775144 0.169794491000172 3.615 19.86394896867 0.30 1.35419222406125 3.5 1.547 2.270 1.544 1.48 2.000 0.30 0.06 3.2 0 ``` ## References [1] Justo, Jo\~ao F. and Bazant, Martin Z. and Kaxiras, Efthimios and Bulatov, V.V. and Yip, Sidney, "Interatomic potential for silicon defects and disordered phases," Phys. Rev.B, 58(5):2539--2550, Aug 1998. [2] Marks, N.A., "Generalizing the environment-dependent interaction potential for carbon," Phys. Rev.B, 63(3):035401, Dec 2000. [3] Christie, H.J. and Robinson, M. and Roach, D.L. and Ross, D.K. and Suarez-Martinez, I. and Marks, N.A., "Simulating radiation damage cascades in graphite," carbon, 81:105--114, 2015. [4] [https://www.lammps.org](https://www.lammps.org) [5] [https://docs.lammps.org/pair_edip.html](https://docs.lammps.org/pair_edip.html) ## License [LGPLv2.1+](https://www.gnu.org/licenses/old-licenses/lgpl-2.1.html) ## Copyright ```txt LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. ``` ```txt Copyright (c) 2021, Regents of the University of Minnesota. All rights reserved. ``` ## Contributing Contributors:\       Yaser Afshar (UMN) Contributing authors, original code in the LAMMPS MANYBODY (edip):\       Luca Ferraro (CASPUR) Contributing authors, original code in the LAMMPS MANYBODY (edip\multi):\       Chao Jiang Contributing authors, original code (edip\c):\       Nigel Marks (Curtin University)\       Sam McSweeney (Curtin University)