/* */ /* 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 */ /* */ /* */ /* Portions Copyright (c) 2019, Regents of the University of Minnesota. */ /* All rights reserved. */ /* */ /* Contributors: */ /* Daniel S. Karls */ /* Ellad B. Tadmor */ /* Ryan S. Elliott */ /* */ /* Portions Copyright (c) 2019, Regents of the University of Minnesota. */ /* All rights reserved. */ /* */ /* Contributors: Josh Vievering */ /* */ #include /* */ /* Model driver for the Kaxiras-Pandey (KP) potential */ /* */ /* Source: */ /* Kaxiras E, Pandey KC (1988), */ /* New classical potential for accurate simulation of atomic processes in Si. */ /* Physical Review B 38 12736-12739. */ /* doi: 10.1103/PhysRevB.38.12736 */ /* */ /* The following flag indicates whether or not to do an explicit check to */ /* determine if the distance between neighbors j and k of atom i in a */ /* three-body computations exceeds the (single) cutoff distance. The check */ /* is skipped if the flag is set to 1, and performed if the flag is set to 0. */ /* This feature is provided for efficiency reasons. */ /* */ #define SKIP_CHECK_ON_RJK 1 /* The following enumeration provides a named index to each element in your */ /* parameter array. The last entry *must* be NUM_PARAMS (by construction by */ /* being at the end of the list, its value will be the number of parameters). */ /* These names will be used in the functions below to unpack and access */ /* parameters. The order of parameters *is* important. It must match the */ /* ordering in the parameter file read in by the model driver. */ /* */ enum PARAMS { PARAM_A1, PARAM_A2, PARAM_alpha1, PARAM_alpha2, PARAM_B, PARAM_D, PARAM_beta, PARAM_C, NUM_PARAMS }; /* The following array of strings contains names and descriptions of the */ /* parameters that will be registered in the KIM API object. The number of */ /* entries must match the number of PARAM_* lines above, and the order must */ /* correspond to the ordering above. */ /* */ static char const * const param_strings[][2] = {{"A1", "Prefactor of first exponential term in pairwise interaction"}, {"A2", "Prefactor of second exponential term in pairwise interaction"}, {"alpha1", "Prefactor of Rij^2 in first exponential term of pairwise interaction"}, {"alpha2", "Prefactor of Rij^2 in second exponential term of pairwise interaction"}, {"B", "Prefactor of g^2(theta) in three-body term"}, {"D", "Prefactor of g^4(theta) in three-body term"}, {"beta", "Prefactor of (Rij^2 + Rik^2) in three-body term"}, {"C", "Cutoff radius"}}; /* The following two variables define the base unit system required to be */ /* used by the parameter file. */ /* */ static KIM_LengthUnit const * const length_unit = &KIM_LENGTH_UNIT_A; static KIM_EnergyUnit const * const energy_unit = &KIM_ENERGY_UNIT_eV; /* The following array of double precision values define the unit */ /* exponents for each parameter. The first number is the length exponent */ /* and the second number is the energy exponent. The number of entries must */ /* match the number of PARAM_* lines above, and the order must correspond */ /* to the ordering above. Use two zeros for unitless parameters. */ /* */ static double const param_units[][2] = {{0.0, 1.0}, /* A1 */ {0.0, 1.0}, /* A2 */ {-2.0, 0.0}, /* alpha1 */ {-2.0, 0.0}, /* alpha2 */ {0.0, 1.0}, /* B */ {0.0, 1.0}, /* D */ {-2.0, 0.0}, /* beta */ {1.0, 0.0}}; /* C */ static double get_influence_distance(double const * const params) { double const C = params[PARAM_C]; return C; } /* Calculate two-body energy phi2(r) and its derivative w.r.t. Rij */ static void calc_phi2_dphi2(double const * const params, double const Rij, double * const phi2, double * const dphi2_dRij) { /* Unpack parameters */ double const A1 = params[PARAM_A1]; double const A2 = params[PARAM_A2]; double const alpha1 = params[PARAM_alpha1]; double const alpha2 = params[PARAM_alpha2]; double const C = params[PARAM_C]; if (Rij < C) { *phi2 = (A1 * exp(-1*alpha1*pow(Rij,2))) - (A2 * exp(-1*alpha2*pow(Rij,2))); } else { *phi2 = 0.0; } if (Rij < C) { if (dphi2_dRij != NULL) { *dphi2_dRij = (A1 * exp(-1*alpha1*pow(Rij,2)) * -2 * alpha1 * Rij) - (A2 * exp(-1*alpha2*pow(Rij,2)) * -2 * alpha2 * Rij); }} else { *dphi2_dRij = 0.0; } return; } /* Calculate three-body energy phi3(Rij, Rik, Rjk) and its derivatives w.r.t Rij, Rik, and Rjk */ static void calc_phi3_dphi3(double const * const params, double const Rij, double const Rik, double const Rjk, double * const phi3, double * const dphi3_dRij, double * const dphi3_dRik, double * const dphi3_dRjk) { /* Unpack parameters */ double const beta = params[PARAM_beta]; double const B = params[PARAM_B]; double const D = params[PARAM_D]; double const C = params[PARAM_C]; if (Rij < C && Rik < C) { *phi3 = exp(-1 * beta * (pow(Rij,2) + pow(Rik,2))) * ((B * pow(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) +(1.0 / 3.0),2)) + (D * pow(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) + (1.0 / 3.0),4))); } else { *phi3 = 0.0; } if (Rij < C && Rik < C) { if (dphi3_dRij != NULL) { *dphi3_dRij = (exp(-1 * beta * (pow(Rij,2) + pow(Rik,2))) * ((B * pow(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij))+(1.0 / 3.0),2)) + (D * pow(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) + (1.0 / 3.0),4))) * (-2 * beta * Rij)) + (exp(-1 * beta * (pow(Rij,2) + pow(Rik,2))) * ((2*B*(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) + (1.0/3.0))) + (4*D*pow(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) + (1.0/3.0),3)))*((pow(Rij,2) - pow(Rik,2) + pow(Rjk,2))/(2*pow(Rij,2)*Rik))); *dphi3_dRik = (exp(-1 * beta * (pow(Rij,2) + pow(Rik,2))) * ((B * pow(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij))+(1.0 / 3.0),2)) + (D * pow(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) + (1.0 / 3.0),4))) * (-2 * beta * Rik)) + (exp(-1 * beta * (pow(Rij,2) + pow(Rik,2))) * ((2*B*(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) + (1.0/3.0))) + (4*D*pow(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) + (1.0/3.0),3)))*((-pow(Rij,2) + pow(Rik,2) + pow(Rjk,2))/(2*pow(Rik,2)*Rij))); *dphi3_dRjk = (exp(-1 * beta * (pow(Rij,2) + pow(Rik,2))) * ((2*B*(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) + (1.0/3.0))) + (4*D*pow(((pow(Rjk,2) - pow(Rik,2) - pow(Rij,2)) / (-2*Rik*Rij)) + (1.0/3.0),3)))*(-Rjk / (Rij*Rik))); }} else { *dphi3_dRij = 0.0; *dphi3_dRik = 0.0; *dphi3_dRjk = 0.0; } return; }