#
# 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 (Four_Body_Mistriotis_Flytzanis_Farantos) contains  Mistriotis-Flytzanis-Farantos potential Model driver written in C.
This can handle maximum 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. 

It complies with the KIM API interface standard, using all neihgbor handling methods
with full neighbor lists (NEIGH_PURE_F, NEIGH_RVEC_F, MI_OPBC_F and CLUSTER) and
iterator or locator mode. It publishes its parameters and supports optional computation
of `energy', 'forces', 'particleEnergy', 'process_dEdr', and 'process_dE2dr2'.
Its initializaton routine takes a KIM API object pointer and a pointer to a string that
contains all parameters which define a KIM Model based on this KIM Model Driver.

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: Any comments, e.g. parameters for nInteractions below starting with 1st, or blank space. More blank lines may follow.
   Line 4: A 
   Line 5: B 
   Line 6: p
   Line 7: q
   Line 8: a 
   Line 9: lambda 
   Line 10: lambda_2 
   Line 11: gamma
   Line 12: sigma in Angstrom
   Line 13: epsilon in ev 
   Line 14: Q
   Line 15: costheta_0  

# If nSpecies = 1, then ignore additional lines otherwise
# continue listing the parameters for 2nd interaction 
   Line 18: A 
   Line 19: B 
   Line 20: p
   Line 21: q
   Line 22: a 
   Line 23: lambda 
   Line 24: lambda_2 
   Line 25: gamma
   Line 26: sigma in Angstrom
   Line 27: epsilon in ev 
   Line 28: Q
   Line 29: 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 term these sets must have 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 term these sets have 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 term at least these sets have 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. 

The following files are in this directory:

LICENSE.CDDL
     The Common Development and Distribution License (CDDL) Version 1.0 file

Makefile
    makefile to compile and build executables with GNU or Intel compilers

README
    This file

Four_Body_Mistriotis_Flytzanis_Farantos.c
    C file that contains the implementation of the 4-body Mistriotis-Flytzanis-Farantos potential 
    Model Driver which can habdle two different species

Four_Body_Mistriotis_Flytzanis_Farantos.kim.tpl
    Template KIM descriptor file containing all variables and methods that are
    exposed by the driver through the KIM API interface.  This file will be
    used (by the make system) to generate a KIM descriptor file for each KIM
    Model based on this Model Driver.