1/* */ 2/* CDDL HEADER START */ 3/* */ 4/* The contents of this file are subject to the terms of the Common */ 5/* Development and Distribution License Version 1.0 (the "License"). */ 6/* */ 7/* You can obtain a copy of the license at */ 8/* http://www.opensource.org/licenses/CDDL-1.0. See the License for the */ 9/* specific language governing permissions and limitations under the License. */ 10/* */ 11/* When distributing Covered Code, include this CDDL HEADER in each file and */ 12/* include the License file in a prominent location with the name */ 13/* LICENSE.CDDL. If applicable, add the following below this CDDL HEADER, */ 14/* with the fields enclosed by brackets "[]" replaced with your own */ 15/* identifying information: */ 16/* */ 17/* Portions Copyright (c) [yyyy] [name of copyright owner]. */ 18/* All rights reserved. */ 19/* */ 20/* CDDL HEADER END */ 21/* */ 22 23/* */ 24/* Portions Copyright (c) 2019, Regents of the University of Minnesota. */ 25/* All rights reserved. */ 26/* */ 27/* Contributors: */ 28/* Daniel S. Karls */ 29/* Ellad B. Tadmor */ 30/* Ryan S. Elliott */ 31/* */ 32/* Portions Copyright (c) 2019, Regents of the University of Minnesota. */ 33/* All rights reserved. */ 34/* */ 35/* Contributors: */ 36/* Anshul Chawla */ 37/* */ 38 39#include <math.h> 40 41/* Model Driver for the Khor-Das Sarma three-body bond-order potential */ 42/* */ 43/* Source: */ 44/* K.E.Khor and S. Das Sarma */ 45/* "Proposed universal interatomic potential for elemental tetrahedrally */ 46/* bonded semiconductors", Physical Review B, Vol. 38, 3318-3322, 1988. */ 47 48/* The following flag indicates whether or not to do an explicit check to */ 49/* determine if the distance between neighbors j and k of atom i in a */ 50/* three-body computations exceeds the (single) cutoff distance. The check */ 51/* is skipped if the flag is set to 1, and performed if the flag is set to 0. */ 52/* This feature is provided for efficiency reasons. */ 53/* */ 54#define SKIP_CHECK_ON_RJK 1 55 56/* The following enumeration provides a named index to each element in your */ 57/* parameter array. The last entry *must* be NUM_PARAMS (by construction by */ 58/* being at the end of the list, its value will be the number of parameters). */ 59/* These names will be used in the functions below to unpack and access */ 60/* parameters. The order of parameters *is* important. It must match the */ 61/* ordering in the parameter file read in by the model driver. */ 62/* */ 63enum PARAMS { 64 PARAM_A, 65 PARAM_B, 66 PARAM_theta, 67 PARAM_lambda, 68 PARAM_alpha, 69 PARAM_beta, 70 PARAM_gamma, 71 PARAM_eta, 72 PARAM_R, 73 PARAM_thetanot, 74 NUM_PARAMS 75}; 76 77#include "bondorder_aux.inc" 78 79/* The following array of strings contains names and descriptions of the */ 80/* parameters that will be registered in the KIM API object. The number of */ 81/* entries must match the number of PARAM_* lines above, and the order must */ 82/* correspond to the ordering above. */ 83/* */ 84static char const * const param_strings[][2] 85 = {{"A", "Amplitude of the cutoff function"}, 86 {"B", "Parameter in the prefactor b_ij to the attractive pairwise interaction"}, 87 {"theta", "Exponent in the repulsive pairwise interaction function f_R"}, 88 {"lambda", "Exponent in the attractive pairwise interaction function f_A"}, 89 {"alpha", "Parameter in the prefactor b_ij to the attractive pairwise interaction"}, 90 {"beta", "Parameter in the prefactor a_ij and b_ij to the repulsive pairwise interaction and attractive pairwise interaction respectively"}, 91 {"gamma", "Parameter in the prefactor a_ij and b_ij to the repulsive pairwise interaction and attractive pairwise interaction respectively"}, 92 {"eta", "Parameter in the angular function g(theta)"}, 93 {"R", "minimum interatomic distance"}, 94 {"thetanot","Equilibrium bond angle"}}; 95 96/* The following two variables define the base unit system required to be */ 97/* used by the parameter file. */ 98/* */ 99static KIM_LengthUnit const * const length_unit = &KIM_LENGTH_UNIT_A; 100static KIM_EnergyUnit const * const energy_unit = &KIM_ENERGY_UNIT_eV; 101 102/* The following array of double precision values define the unit */ 103/* exponents for each parameter. The first number is the length exponent */ 104/* and the second number is the energy exponent. The number of entries must */ 105/* match the number of PARAM_* lines above, and the order must correspond */ 106/* to the ordering above. Use two zeros for unitless parameters. */ 107/* */ 108static double const param_units[][2] = {{0.0, 1.0}, /* A length energy */ 109 {0.0, 0.0}, /* B */ 110 {-1.0, 0.0}, /* theta */ 111 {-1.0, 0.0}, /* lambda */ 112 {0.0, 0.0}, /* alpha */ 113 {-1.0, 0.0}, /* beta (gamma) */ 114 {0.0, 0.0}, /* gamma */ 115 {0.0, 0.0}, /* eta */ 116 {1.0, 0.0}, /* R */ 117 {0.0, 0.0}}; /* theta */ 118 119static double get_influence_distance(double const * const params) 120{ 121 double const beta = params[PARAM_beta]; 122 double const gamma = params[PARAM_gamma]; 123 double const R = params[PARAM_R]; 124 125 return (R + pow(-log(1e-16)/beta, 1/gamma)); 126} 127 128/* Calculate coordination Zi and its derivative w.r.t. Rij */ 129static void calc_Zi_dZi(double const * const params, 130 double const Rij, 131 double * const Z_i, 132 double * const dZi_dRij) 133{ 134 /* Unpack parameters */ 135 136 double const beta = params[PARAM_beta]; 137 double const gamma = params[PARAM_gamma]; 138 double const R = params[PARAM_R]; 139 140 if (Rij <= R) 141 { 142 *Z_i = 1.0; 143 144 if (dZi_dRij != NULL) 145 { 146 *dZi_dRij = 0.0; 147 } 148 } 149 else 150 { 151 *Z_i = exp(-beta * pow(Rij - R, gamma)); 152 153 if (dZi_dRij != NULL) 154 { 155 *dZi_dRij = -beta * gamma * pow(Rij - R, gamma - 1.0) * (*Z_i); 156 } 157 } 158 159 return; 160} 161 162/* Calculate bond order bij and its derivative w.r.t. zeta_ij */ 163static void calc_bij_dbij(double const * const params, 164 double const Z_i, 165 double const zeta_ij, 166 double * const bij, 167 double * const dbij_dZi, 168 double * const dbij_dzeta_ij) 169{ 170 /* Unpack parameters */ 171 double const B = params[PARAM_B]; 172 double const alpha = params[PARAM_alpha]; 173 174 double const Zi_pow_alpha = pow(Z_i, alpha); 175 176 *bij = (B/Zi_pow_alpha) *(1+ zeta_ij); 177 178 if (dbij_dZi != NULL) 179 { 180 *dbij_dZi = -(alpha/(Z_i)) * *bij; 181 *dbij_dzeta_ij = (B/Zi_pow_alpha); 182 } 183 184 return; 185} 186 187/* Calculate two-body energy phi2(r) and its derivative w.r.t. Rij */ 188static void calc_phi2_dphi2(double const * const params, 189 double const Rij, 190 double const bij, 191 double * const phi2, 192 double * const dphi2_dRij, 193 double * const dphi2_dbij) 194{ 195 *phi2 = fc(params, Rij) * (fR(params, Rij) + bij * fA(params, Rij)); 196 197 if (dphi2_dRij != NULL) 198 { 199 *dphi2_dRij 200 = fc(params, Rij) 201 * (dfR_dRij(params, Rij) + bij * dfA_dRij(params, Rij)) 202 + dfc_dR(params, Rij) * (fR(params, Rij) + bij * fA(params, Rij)); 203 204 *dphi2_dbij = fc(params, Rij) * fA(params, Rij); 205 } 206 207 return; 208} 209 210/* Calculate three-body energy phi3(Rij, Rik, Rjk) and its derivatives w.r.t 211 * Rij, Rik, and Rjk */ 212static void calc_phi3_dphi3(double const * const params, 213 double const Rij, 214 double const Rik, 215 double const Rjk, 216 double * const phi3, 217 double * const dphi3_dRij, 218 double * const dphi3_dRik, 219 double * const dphi3_dRjk) 220{ 221 /* Unpack parameters */ 222 223 double theta_jik; 224 double gval; 225 double dcosthetajik_dRij; 226 double dcosthetajik_dRik; 227 double dcosthetajik_dRjk; 228 double dg_dcosthetaval; 229 230 theta_jik = acos((Rij * Rij + Rik * Rik - Rjk * Rjk) / (2 * Rij * Rik)); 231 gval = g(params, theta_jik); 232 233 *phi3 = gval; 234 235 if (dphi3_dRij != NULL) 236 { 237 dg_dcosthetaval = dg_dcostheta(params, theta_jik); 238 239 dcosthetajik_dRij 240 = (Rij * Rij - Rik * Rik + Rjk * Rjk) / (2 * Rij * Rij * Rik); 241 dcosthetajik_dRik 242 = (-Rij * Rij + Rik * Rik + Rjk * Rjk) / (2 * Rij * Rik * Rik); 243 dcosthetajik_dRjk = -Rjk / (Rij * Rik); 244 245 *dphi3_dRij = (dg_dcosthetaval * dcosthetajik_dRij); 246 247 *dphi3_dRik = (dg_dcosthetaval * dcosthetajik_dRik); 248 249 *dphi3_dRjk = (dg_dcosthetaval * dcosthetajik_dRjk); 250 251 252 } 253 254 return; 255} 256