/* */ /* 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. */ /* Portions Copyright (c) Year, Organization */ /* */ /* Contributors: */ /* Ramanish Singh */ /* */ #include /* */ /* Model driver for the three-body bond order Wang-Rockett potential */ /* */ /* Source: */ /* J. Wang and A. Rockett (1991), */ /* Simulating diffusion on Si (001) 2×1 surfaces using a modified interatomic */ /* potential, */ /* Phys. Rev. B 43, 12 571 (1991). */ /* */ /* 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_A, PARAM_B, PARAM_lambda1, PARAM_lambda2, PARAM_lambda3, PARAM_alpha, PARAM_beta, PARAM_n, PARAM_c, PARAM_d, PARAM_h, PARAM_R, PARAM_D, PARAM_Ra, PARAM_Da, NUM_PARAMS }; #include "bondorder_aux.inc" /* 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] = {{"A", "Amplitude in the repulsive pairwise interaction function f_R"}, {"B", "Amplitude in the attractive pairwise interaction function f_A"}, {"lambda1", "Exponent in the repulsive pairwise interaction function f_R"}, {"lambda2", "Exponent in the attractive pairwise interaction function f_A"}, {"lambda3", "Exponent in the zeta_ij term in the bond-order function b_ij"}, {"alpha", "Parameter in the prefactor a_ij to the repulsive pairwise interaction"}, {"beta", "Parameter in the prefactor b_ij to the attractive pairwise interaction"}, {"n", "Exponent in the prefactor b_ij to the attractive pairwise interaction"}, {"c", "Parameter in the angular function g(theta)"}, {"d", "Parameter in the angular function g(theta)"}, {"h", "Parameter in the angular function g(theta)"}, {"R", "Internal cutoff radius below which cutoff function f_C equals one"}, {"D", "Width of which cutoff function f_C goes from one to zero"}, {"Ra", "Internal cutoff radius below which Aij goes from one to Bij"}, {"Da", "Width of which function Aij ges from one to Bij"}}; /* 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}, /* A */ {0.0, 1.0}, /* B */ {-1.0, 0.0}, /* lambda1 */ {-1.0, 0.0}, /* lambda2 */ {-1.0, 0.0}, /* lambda3 */ {0.0, 0.0}, /* alpha */ {0.0, 0.0}, /* beta */ {0.0, 0.0}, /* n */ {0.0, 0.0}, /* c */ {0.0, 0.0}, /* d */ {0.0, 0.0}, /* h */ {1.0, 0.0}, /* R */ {1.0, 0.0}, /* D */ {1.0, 0.0}, /* Ra */ {1.0, 0.0}}; /* Da */ static double get_influence_distance(double const * const params) { double const R = params[PARAM_R]; double const D = params[PARAM_D]; return R + D; } /* Calculate bond order bij and its derivative w.r.t. zeta_ij */ static void calc_bij_dbij(double const * const params, double const zeta_ij, double * const bij, double * const dbij_dzeta_ij) { /* Unpack parameters */ double const beta = params[PARAM_beta]; double const n = params[PARAM_n]; double const betazeta_pow_n = pow(beta * zeta_ij, n); *bij = pow(1.0 + betazeta_pow_n, -0.5 / n); if (dbij_dzeta_ij != NULL) { *dbij_dzeta_ij = -(0.5 / n) * pow(1 + betazeta_pow_n, (-0.5 / n) - 1.0) * n * betazeta_pow_n / zeta_ij; } return; } /* Calculate bond order aij and its derivative w.r.t. zeta_ij */ static void calc_aij_daij(double const * const params, double const Rij, double const bij, double * const aij, double * const daij_dRij, double * const daij_dbij) { /* Unpack parameters */ double const A = params[PARAM_A]; double const Ra = params[PARAM_Ra]; double const Da = params[PARAM_Da]; if (Rij <= Ra - Da) { *aij = 1.0; if (daij_dRij != NULL) { *daij_dRij = 0.0; *daij_dbij = 0.0; } } else if (Ra - Da < Rij && Rij < Ra+ Da) { *aij = 1 * (bij + (1-bij) * (0.5 - 0.5*sin(0.5*PI*(Rij-Ra)/Da))); if (daij_dRij != NULL) { *daij_dRij = -0.25*(PI/Da) * (1 - bij)*cos(0.5*PI*(Rij-Ra)/Da); *daij_dbij = 0.5 + 0.5 * sin(0.5*PI*(Rij-Ra)/Da); /* *daij_dRij *= 1.0; */ /* *daij_dbij *= 1.0; */ } } else { *aij = 1 * bij; if (daij_dRij != NULL) { *daij_dRij = 0.0; *daij_dbij = 1.0; } } return; } /* 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 aij, double const bij, double * const phi2, double * const dphi2_dRij, double * const dphi2_daij, double * const dphi2_dbij) { *phi2 = fc(params, Rij) * (aij * fR(params, Rij) + bij * fA(params, Rij)); if (dphi2_dRij != NULL) { *dphi2_dRij = fc(params, Rij) * ( aij * dfR_dRij(params, Rij) + bij * dfA_dRij(params, Rij)) + dfc_dR(params, Rij) * (aij * fR(params, Rij) + bij * fA(params, Rij)); *dphi2_daij = fc(params, Rij) * fR(params, Rij); *dphi2_dbij = fc(params, Rij) * fA(params, Rij); } 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 lambda3 = params[PARAM_lambda3]; double costheta_jik; double gval; double arg_exp; double exponential; double dcosthetajik_dRij; double dcosthetajik_dRik; double dcosthetajik_dRjk; double dexponential_dRij; double dexponential_dRik; double dg_dcosthetaval; costheta_jik = (Rij * Rij + Rik * Rik - Rjk * Rjk) / (2 * Rij * Rik); gval = g(params, costheta_jik); arg_exp = lambda3 * (Rij - Rik); arg_exp = arg_exp * arg_exp * arg_exp; exponential = exp(arg_exp); *phi3 = fc(params, Rik) * gval * exponential; if (dphi3_dRij != NULL) { dg_dcosthetaval = dg_dcostheta(params, costheta_jik); dcosthetajik_dRij = (Rij * Rij - Rik * Rik + Rjk * Rjk) / (2 * Rij * Rij * Rik); dcosthetajik_dRik = (-Rij * Rij + Rik * Rik + Rjk * Rjk) / (2 * Rij * Rik * Rik); dcosthetajik_dRjk = -Rjk / (Rij * Rik); dexponential_dRij = (lambda3 * lambda3 * lambda3) * 3 * (Rij - Rik) * (Rij - Rik) * exponential; dexponential_dRik = -dexponential_dRij; *dphi3_dRij = fc(params, Rik) * (gval * dexponential_dRij + exponential * dg_dcosthetaval * dcosthetajik_dRij); *dphi3_dRik = fc(params, Rik) * (gval * dexponential_dRik + exponential * dg_dcosthetaval * dcosthetajik_dRik) + dfc_dR(params, Rik) * gval * exponential; *dphi3_dRjk = fc(params, Rik) * dg_dcosthetaval * dcosthetajik_dRjk * exponential; } return; }