#
# CDDL HEADER START
#
# The contents of this file are subject to the terms of the Common Development
# and Distribution License Version 1.0 (the "License").
#
# You can obtain a copy of the license at
# http://www.opensource.org/licenses/CDDL-1.0.  See the License for the
# specific language governing permissions and limitations under the License.
#
# When distributing Covered Code, include this CDDL HEADER in each file and
# include the License file in a prominent location with the name LICENSE.CDDL.
# If applicable, add the following below this CDDL HEADER, with the fields
# enclosed by brackets "[]" replaced with your own identifying information:
#
# Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
#
# CDDL HEADER END
#

#
# Copyright (c) 2018, Regents of the University of Minnesota.
#
# Contributors:
#    Mingjian Wen
#

This directory contains three-body Stillinger-Weber potential Model
Driver for transition metal dichalcogenide monolayers of the type MX2
(e.g. MoS2 and WSe2).


1. The functional form has been rewritten from the standard SW functional form
   so as to separate the parameters to avoid inconvenience of using this Model
   Driver. The variable 'epsilon' in the standard SW potential is a dummy variable,
   so it is combined with 'A' and 'lambda'. The 'cutoff' of the standard SW potential
   is determined implicitly by 'cutoff' = 'a*sigma', where 'a' and 'sigma' are
   potential parameters. In this Model Driver, the 'cutoff' is defined explicitly.
   The benefit is that if 'sigma' is updated there is no need to update 'cutoff'.
   Above all, compared to the standard SW potential, the following redefinitions
   have been made:

   A      := A*epsilon
   lambda := lambda*epsilon
   gamma  := gamma*sigma
   cutoff := a*sigma

   Finally, it is assumed that the angle 'beta0' is the same for all three-body
   interactions.

   Consequently, the SW Driver for multiple species can be written as:

   E = sum_{i,j>i} phi_2(rij) + sum_{i, j!=i, k>j} phi_3(rij, rik, beta)

   phi_2(rij) = Aij*(Bij(rij/sigma_ij)^(-p) - (rij/sigma_ij)^(-q))
                * exp(sigma_ij/(rij - cutoff_ij))

   phi_3(rij,rik,beta) = lambda_ijk*(cos[beta] - cos[beta0])^2
                          * exp(gamma_ij/(rij - cutoff_ij) + gamma_ik/(rik - cutoff_ik))


2. Three-body bond angle bending is considered only for certain interations. Bond
   angle ijk (i is the vertex atom at which bond ij and bond ik form the angle)
   is considered only when atoms j and k are of the same species, but different from
   the species of atom i.  Explicitly, only bond angle of types spec1-spec2-spec2
   and spec2-spec1-spec1 contribute energy to the system. Besides the usual cutoffs
   cutoff_ij and cutoff_ik employed in the standard SW potential, additional cutoff
   cutoff_jk is applied to bond jk. So, once either of the bond length rij, rik, or
   rjk is larger than its corresponding cutoff, the 3-body term vanishes (equals 0).


3. One parameter file is needed, and it should be in the following format:

   First line: number of species (should be 2)

   Lines 2~4: parameters for 2-body interaction
     species1 species2 A B p q sigma gamma cutoff

   Lines 5~6: parameters for 3-body interaction
     species1 species2 species3 lambda_ijk cos_beta0 cutoff_jk

   species is valid KIM API particle species string
   A and lambda_ijk in [eV]
   sigma, gamma, cutoff, and cutoff_jk in [Angstrom]
   others are unitless


1. F. H. Stillinger and T. A. Weber, "Computer simulation of local order in condensed
   phases of silicon", Phys. Rev. B, vol. 31, 5262-5271, 1985

2. Ellad B. Tadmor and Ronald E. Miller, Modeling Materials: Continuum, Atomistic and
   Multiscale Techniques, Cambridge University Press, 2011

3. Mingjian Wen, Sharmila N. Shirodkar, Petr Plechac, Efthimios Kaxiras, Ryan S. Elliott
   and Ellad B. Tadmor, "Stillinger-Weber potential for MoS2: parameterization and
   sensitivity analysis", J. Appl. Phys., 122, 244301, 2017