/* */ /* 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: */ /* Dipanjan Ghosh /* */ #include /* */ /* Auxiliary files for the Biswas-Hamann model driver */ /* */ /* Functions in this file are only used by cluster.inc. They are not required */ /* by ThreeBodyCluster.c. */ /* */ /* */ /* Define functions used in two-body calculations */ /* */ static double fc(double const * const params, double const r) { double const rc = params[PARAM_rc]; double const mu = params[PARAM_mu]; return 1.0 / (1.0 + exp((r-rc)/mu)); } static double dfc(double const * const params, double const r) { double const rc = params[PARAM_rc]; double const mu = params[PARAM_mu]; double const expterm = exp((r-rc)/mu); double const fc = 1.0 / (1.0 + expterm); return (-1.0/mu) * fc * fc * expterm; } static void f2_df2(double const * const params, double const r, double * const f2, double * const df2_dr) { /* Unpack parameters */ double const A1 = params[PARAM_A1]; double const A2 = params[PARAM_A2]; double const lambda1 = params[PARAM_lambda1]; double const lambda2 = params[PARAM_lambda2]; double const term1 = exp(-lambda1*r*r); double const term2 = exp(-lambda2*r*r); double const cutoff_term = fc(params, r); double const dcutoff = dfc(params, r); double const f2_without_co = A1*term1 + A2*term2; *f2 = f2_without_co*cutoff_term; if (df2_dr != NULL) { *df2_dr = (f2_without_co * dcutoff) + (cutoff_term*((-2*A1*lambda1*r*term1)+(-2*A2*lambda2*r*term2))); } } /* */ /* Define functions used in three-body calculations */ /* */ static void f3_df3(double const * const params, double const Rij, double const Rik, double const Rjk, double * const f3, double * const df3_dRij, double * const df3_dRik, double * const df3_dRjk) { /* Unpack parameters */ double const B1 = params[PARAM_B1]; double const B2 = params[PARAM_B2]; double const alpha1 = params[PARAM_alpha1]; double const alpha2 = params[PARAM_alpha2]; double costheta_jik; double B1_exp1; double B1_exp2; double B2_exp1; double B2_exp2; double costerm; double costermsq; double costermcu; double f2c_Rij; double df2c_dRij; double f2c_Rik; double df2c_dRik; double dcostheta_dRij; double dcostheta_dRik; double dcostheta_dRjk; double common_term_11; double df3_ij_term1; double df3_ij_term2; double df3_ij_term3; double common_term_12; double df3_ij_term4; double df3_ij_term5; double df3_ij_term6; double common_term_21; double df3_ik_term1; double df3_ik_term2; double df3_ik_term3; double common_term_22; double df3_ik_term4; double df3_ik_term5; double df3_ik_term6; double common_term_3; double df3_jk_term1; double df3_jk_term2; /* Law of cosines to get angle from distances */ costheta_jik = (Rij * Rij + Rik * Rik - Rjk * Rjk) / (2 * Rij * Rik); B1_exp1 = exp(-alpha1*Rij*Rij); B1_exp2 = exp(-alpha1*Rik*Rik); B2_exp1 = exp(-alpha2*Rij*Rij); B2_exp2 = exp(-alpha2*Rik*Rik); costerm = (costheta_jik + (1.0/3.0)); costermsq = costerm * costerm; costermcu = costerm * costermsq; f2c_Rij = fc(params, Rij); f2c_Rik = fc(params, Rik); df2c_dRij = dfc(params, Rij); df2c_dRik = dfc(params, Rik); *f3 = (B1*B1_exp1*B1_exp2*costermsq + B2*B2_exp1*B2_exp2*costermcu)*f2c_Rij*f2c_Rik; if (df3_dRij != NULL) { /* Derivatives of costheta_jik wrt Rij, Rjk and Rik */ dcostheta_dRij = (Rij * Rij - Rik * Rik + Rjk * Rjk) / (2 * Rij * Rij * Rik); dcostheta_dRik = (-Rij * Rij + Rik * Rik + Rjk * Rjk) / (2 * Rij * Rik * Rik); dcostheta_dRjk = -Rjk / (Rij * Rik); /* Calculation of df3_dRij*/ common_term_11 = B1*B1_exp1*B1_exp2*f2c_Rik; df3_ij_term1 = costermsq*df2c_dRij; df3_ij_term2 = -2*alpha1*Rij*costermsq*f2c_Rij; df3_ij_term3 = 2*costerm*dcostheta_dRij*f2c_Rij; common_term_12 = B2*B2_exp1*B2_exp2*f2c_Rik; df3_ij_term4 = costermcu*df2c_dRij; df3_ij_term5 = -2*alpha2*Rij*costermcu*f2c_Rij; df3_ij_term6 = 3*costermsq*dcostheta_dRij*f2c_Rij; *df3_dRij = (common_term_11*(df3_ij_term1 + df3_ij_term2 + df3_ij_term3)) + (common_term_12*(df3_ij_term4 + df3_ij_term5 + df3_ij_term6)); /* Calculation of df3_dRik*/ common_term_21 = B1*B1_exp1*B1_exp2*f2c_Rij; df3_ik_term1 = costermsq*df2c_dRik; df3_ik_term2 = -2*alpha1*Rik*costermsq*f2c_Rik; df3_ik_term3 = 2*costerm*dcostheta_dRik*f2c_Rik; common_term_22 = B2*B2_exp1*B2_exp2*f2c_Rij; df3_ik_term4 = costermcu*df2c_dRik; df3_ik_term5 = -2*alpha2*Rik*costermcu*f2c_Rik; df3_ik_term6 = 3*costermsq*dcostheta_dRik*f2c_Rik; *df3_dRik = (common_term_21*(df3_ik_term1 + df3_ik_term2 + df3_ik_term3)) + (common_term_22*(df3_ik_term4 + df3_ik_term5 + df3_ik_term6)); /* Calculation of df3_dRjk*/ common_term_3 = f2c_Rij * f2c_Rik * dcostheta_dRjk; df3_jk_term1 = B1*B1_exp1*B1_exp2*2.*costerm; df3_jk_term2 = B2*B2_exp1*B2_exp2*3.*costermsq; *df3_dRjk = common_term_3 * (df3_jk_term1 + df3_jk_term2); } return ; }