#
# 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
#
#
#
# Author: Daniel S. Karls
#
# Copyright (c) 2012-2021, Regents of the University of Minnesota.  All rights
# reserved.

This C-based implementation of the EDIP bond-order potential is based on
version 1.1 of the C code published by Prof. Martin Z. Bazant at
http://web.mit.edu/bazant/www/EDIP/index.html.  All functional forms are
evaluated analytically, i.e. no splines are used for evaluation.  It computes
energy, forces, and virial and supports flexible unit handling.  It also
publishes all model parameters in the KIM API so that they may be adjusted by a
simulator.

Notes:

(1) The 'delta' parameter appearing in the original version of this model (for
Si this has a value of delta = 78.7590539) has been eliminated, and instead,
'eta' is used as a free parameter.  'delta' and 'eta' are related through eta =
delta/Q0 (where Q0 is another free parameter), but because the latest EDIP
publication lists a value for eta rather than delta, it seemed more appropriate
to use 'eta' as a free parameter in its own right instead of delta.  For
Silicon, eta = 0.2523244 is the value cited in the EDIP papers.

(2) Since forces are required to compute the virial, the forces will be
computed even if the Test asks only for the virial.

-----------------------------------------

A note on how this model is implemented:

This model was adapted from Prof. Martin Z. Bazant's code (thanks also to Prof.
Joao F. Justo for an additional fortran reference code).

Because of the coordination-dependency of the 2-body and 3-body terms, we're
forced to do multiple loops over all of the atoms. The first inner loop (the
"Pre-pass loop") goes over all of the other atoms and determines which ones are
inside the cutoff radius of the atom being considered by the outer loop.  For
those within the cutoff radius, it computes part of the 2-body interaction (the
part that doesn't require the coordination to be known) and the radial parts of
the 3-body interaction.  Finally, the coordination of each of the atoms within
the cutoff is added onto the total coordination, Z, for that atom.  Having Z,
we now have to loop over all other atoms once more in order to calculate the
coordination-dependent portion of the 2-body energy and get the final 2-body
energy for the atom being considered in the outermost loop, V2.  Note also that
we are unable to perform a "half-summation" for the two-body terms (with an
outer loop over i and an inner loop over j>i) because of the asymmetry the
coordination introduces into the 2-body energy (V2(R_ij,Z_i) != V2(R_ji,Z_j)).
We also need to use Z to perform another set of nested loops in order to
calculate the coordination-dependent portion of the 3-body energy, h(l,Z).
This will give us the total 3-body energy for the atom being considered in the
outermost loop, V3.  Since we are forced to loop over the same atoms multiple
times (first to get the coordination and then to compute the energies), it
makes sense to reduce the computational expense of the subsequent looping by
only looping over the atoms which have already been determined to fall within
the cutoff radius.  To this end, this code generates its own *internal*
neighbor list.

The original EDIP publications can be found at:

     1.  M. Z. Bazant and E. Kaxiras, Phys. Rev. Lett. 77, 4370 (1996).
     2.  M. Z. Bazant, E. Kaxiras, J. F. Justo, Phys. Rev. B 56, 8542 (1997).
     3.  J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov, and S. Yip,
	       Phys. Rev. B 58, 2539 (1998).

If you are interested in an explanation of the looping structure, etc., in this
code and why it was chosen, see chapter 6 of Prof.  Bazant's 1997 Ph.D. thesis,
available at http://web.mit.edu/bazant/www/thesis/md.html.  Page 156 of this
thesis includes a pseudocode of what has been implemented here.

Written by Daniel S. Karls (University of Minnesota). Contact: karl0100 |AT|
umn.edu.