#
# 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) 2012, Regents of the University of Minnesota.  All rights reserved.
#
# Contributors:
#    Amit Singh
#

This directory contains the Mistriotis-Flytzanis-Farantos (MFF) model driver written in C.
This can handle a system containing up to two different species.

/*******************************************************************************

* For four-body potential any potential energy function for N-particle system can be written in the
* following form of two-body, three-body and four-body terms:
     F(1, 2, ...., N)   = sum_{i; 1 <= i \neq j <= N} 0.5*phi_two(r; r = r_ij)
                          + sum_{i; 1 <= i \neq j < k <= N} phi_three(r_ij, r_ik, r_jk);
                          + sum_{i; 1 <= i \neq j < k < l <= N} phi_four(r_ij, r_ik, r_il, r_jk, r_jl, r_kl);

* For modified Stillinger-Weber proposed by "Mistriotis et al, Physical Review B 39, 1212 (1989)",
* the two-body term phi_two can be written as:

     phi_two(r) = epsilon * f_2(r_cap; r_cap = r/sigma); where f_2 is

     f_2(r_cap) = A * ( B*r_cap^(-p) - r_cap^-q ) * exp(1/(r_cap - a)); when  r_cap <  a, where a = cutoff/sigma,
                = 0   when r_cap >= a

* And three-body term phi_three can be written as:

     phi_three(r_ij, r_ik, r_jk) = epsilon * f_3(r1_cap, r2_cap, r3_cap); where

     r1_cap = r_ij/sigma, r2_cap = r_ik/sigma, r3_cap = r_jk/sigma,      and

     f_3(r1_cap, r2_cap, r3_cap) = lambda * (exp(gamma((1/(r1_cap - a)) + (1/(r2_cap - a)))))
                                 * {1 - exp[-Q((costheta_jik - costheta_0)^2)]};  when r1_cap < a && r2_cap < a
                                 = 0;                                             otherwise

     costheta_jik = (r_ij^2 + r_ik^2 - r_jk^2) / (2*r_ij*r_ik).

* and four-body term phi_four can be written as:

     phi_four(r_ij, r_ik, r_il, r_jk, r_jl, r_kl) = epsilon * f_4(r1_cap, r2_cap, r3_cap, r4_cap, r5_cap, r6_cap); where

     r1_cap = r_ij/sigma, r2_cap = r_ik/sigma, r3_cap = r_il/sigma, r4_cap = r_jk/sigma, r5_cap = r_jl/sigma, r6_cap = r_kl/sigma    and

     f_4(r1_cap, r2_cap, r3_cap, r4_cap, r5_cap, r6_cap) = g(r1_cap, r2_cap, r3_cap, costheta_jik, costheta_jil, costheta_kil); where

     g(r1_cap, r2_cap, r3_cap, costheta_jik, costheta_jil, costheta_kil) = lambda_2 * (exp(gamma((1/(r1_cap - a)) + (1/(r2_cap - a)) + (1/(r3_cap - a)))))
                                                                         * {1 - exp[-Q((costheta_jik - costheta_0)^2 + (costheta_jil - costheta_0)^2 + (costheta_kil - costheta_0)^2)]}; when r1_cap < a && r2_cap < a && r3_cap < a
                                                                         = 0;                                                                                                       otherwise

     costheta_jik = (r_ij^2 + r_ik^2 - r_jk^2) / (2*r_ij*r_ik);
     costheta_jil = (r_ij^2 + r_il^2 - r_jl^2) / (2*r_ij*r_il);
     costheta_kil = (r_ik^2 + r_il^2 - r_kl^2) / (2*r_ik*r_il).

*******************************************************************************/

For Silicon [Ref 1 and Ref 3]:
-----------------------------------------------------------------------------------------
| A = 7.041365010799137, B = 0.60107948999545, p = 4, q = 0, a = 1.80                   |
| lambda = 1.727861771058315, lambda_2 = 20.302375809935207, gamma = 1.145530046298506, |
| sigma = 2.0951, epsilon = 2.315 ev (for E_coh = 4.63 ev), Q = 5.0, costheta_0 = -1/3  |
-----------------------------------------------------------------------------------------

References:

1. Mistriotis et al, Physical Review B 39, 1212 (1989)
2. Ellad B. Tadmor and Ronald E. Miller, Modeling Materials:
   Continuum, Atomistic and Multiscale Techniques, Cambridge University Press, 2011
3. epsilon = 2.315 eV, lambda = 4.0/epsilon, lambda_2 = 47.0/epsilon.

This model driver publishes its parameters and supports optional computation
of `energy', 'forces', 'particleEnergy', 'process_dEdr', and 'process_dE2dr2'.

To create a KIM Model from this Model Driver, a parameter file is required.
This file must have the follwing format:

#  Line 1: Any comments
   Line 2: nSpecies, where nSpecies is the number of species and can be either 1 or 2
   Line 3: chemical symbols of the species (separated by a space if nSpecies = 2)
#  Line 4: Any comments, e.g. parameters for nInteractions below starting with 1st, or blank space. More blank lines may follow.
   Line 5: A
   Line 6: B
   Line 7: p
   Line 8: q
   Line 9: a
   Line 10: lambda
   Line 11: lambda_2
   Line 12: gamma
   Line 13: sigma in Angstrom
   Line 14: epsilon in ev
   Line 15: Q
   Line 16: costheta_0

# If nSpecies = 1, then ignore additional lines otherwise
# continue listing the parameters for 2nd interaction
   Line 19: A
   Line 20: B
   Line 21: p
   Line 22: q
   Line 23: a
   Line 24: lambda
   Line 25: lambda_2
   Line 26: gamma
   Line 27: sigma in Angstrom
   Line 28: epsilon in ev
   Line 29: Q
   Line 30: costheta_0

# If nSpecies = 2, then continue listing the parameters for all
# interactions from 3rd to 16th

Any additional lines will be silently ignored.

The interaction indices are provided as per the following conditions:

if (iSpecies == SPEC1 && jSpecies == SPEC1 && kSpecies == SPEC1 && lSpecies == SPEC1) interaction_index = 0;
else if (iSpecies == SPEC1 && jSpecies == SPEC1 && kSpecies == SPEC1 && lSpecies == SPEC2) interaction_index = 1;
else if (iSpecies == SPEC1 && jSpecies == SPEC1 && kSpecies == SPEC2 && lSpecies == SPEC1) interaction_index = 2;
else if (iSpecies == SPEC1 && jSpecies == SPEC1 && kSpecies == SPEC2 && lSpecies == SPEC2) interaction_index = 3;
else if (iSpecies == SPEC1 && jSpecies == SPEC2 && kSpecies == SPEC1 && lSpecies == SPEC1) interaction_index = 4;
else if (iSpecies == SPEC1 && jSpecies == SPEC2 && kSpecies == SPEC1 && lSpecies == SPEC2) interaction_index = 5;
else if (iSpecies == SPEC1 && jSpecies == SPEC2 && kSpecies == SPEC2 && lSpecies == SPEC1) interaction_index = 6;
else if (iSpecies == SPEC1 && jSpecies == SPEC2 && kSpecies == SPEC2 && lSpecies == SPEC2) interaction_index = 7;
else if (iSpecies == SPEC2 && jSpecies == SPEC1 && kSpecies == SPEC1 && lSpecies == SPEC1) interaction_index = 8;
else if (iSpecies == SPEC2 && jSpecies == SPEC1 && kSpecies == SPEC1 && lSpecies == SPEC2) interaction_index = 9;
else if (iSpecies == SPEC2 && jSpecies == SPEC1 && kSpecies == SPEC2 && lSpecies == SPEC1) interaction_index = 10;
else if (iSpecies == SPEC2 && jSpecies == SPEC1 && kSpecies == SPEC2 && lSpecies == SPEC2) interaction_index = 11;
else if (iSpecies == SPEC2 && jSpecies == SPEC2 && kSpecies == SPEC1 && lSpecies == SPEC1) interaction_index = 12;
else if (iSpecies == SPEC2 && jSpecies == SPEC2 && kSpecies == SPEC1 && lSpecies == SPEC2) interaction_index = 13;
else if (iSpecies == SPEC2 && jSpecies == SPEC2 && kSpecies == SPEC2 && lSpecies == SPEC1) interaction_index = 14;
else interaction_index = 15;

For three-body terms, these sets must have the same phi_three parameters, as only ijk iteraction is involved:
(1111,1112),(1121,1122),(1211,1212),(1221,1222)
(2111,2112),(2121,2122),(2211,2212),(2221,2222)

For two-body terms, these sets have the same A, B, p, q parameters, as only ij interaction is involved:
(1111,1112,1121,1122),(1211,1212,1221,1222,2111,2112,2121,2122)
(2211,2212,2221,2222)

For two-body terms, these sets at least have the same a, sigma parameters:
(1111,1112),(1221,1222,2112,2111),(2221,2222)

There might be cases when cutoffs for ij, ijk, ijkl are different.
In that case the code will need some changes.