diff --git a/src/USER-MEAMC/meam_cleanup.cpp b/src/USER-MEAMC/meam_cleanup.cpp index 1beedc9fa0..3b1edc8cd9 100644 --- a/src/USER-MEAMC/meam_cleanup.cpp +++ b/src/USER-MEAMC/meam_cleanup.cpp @@ -1,16 +1,17 @@ extern "C" { #include "meam.h" -void meam_cleanup_(void) { - - deallocate(meam_data.phirar6); - deallocate(meam_data.phirar5); - deallocate(meam_data.phirar4); - deallocate(meam_data.phirar3); - deallocate(meam_data.phirar2); - deallocate(meam_data.phirar1); - deallocate(meam_data.phirar); - deallocate(meam_data.phir); +void +meam_cleanup_(void) +{ + deallocate(meam_data.phirar6); + deallocate(meam_data.phirar5); + deallocate(meam_data.phirar4); + deallocate(meam_data.phirar3); + deallocate(meam_data.phirar2); + deallocate(meam_data.phirar1); + deallocate(meam_data.phirar); + deallocate(meam_data.phir); } } \ No newline at end of file diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp index bdcff4a7f6..587d3bbb5f 100644 --- a/src/USER-MEAMC/meam_dens_final.cpp +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -1,12 +1,14 @@ -extern "C"{ +extern "C" { #include "meam.h" #include - void G_gam(double Gamma, int ibar, double gsmooth_factor, double *G, int *errorflag); - void dG_gam(double Gamma, int ibar, double gsmooth_factor,double *G, double *dG); +void G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, + int* errorflag); +void dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, + double* dG); // in meam_setup_done - void get_shpfcn(double *s /* s(3) */, lattice_t latt); +void get_shpfcn(double* s /* s(3) */, lattice_t latt); // Extern "C" declaration has the form: // @@ -25,241 +27,266 @@ extern "C"{ // dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); // - void meam_dens_final_(int *nlocal, int *nmax, - int *eflag_either, int *eflag_global, int *eflag_atom, double *eng_vdwl, double *eatom, - int *ntype, int *type, int *fmap, - double *Arho1, double *Arho2, double *Arho2b, double *Arho3, double *Arho3b, double *t_ave, double *tsq_ave, - double *Gamma, double *dGamma1, double *dGamma2, double *dGamma3, - double *rho, double *rho0, double *rho1, double *rho2, double *rho3, double *fp, int *errorflag) - { - - arrdim2v(Arho1,3,*nmax) - arrdim2v(Arho2,6,*nmax) - arrdim2v(Arho3,10,*nmax) - arrdim2v(Arho3b,3,*nmax) - arrdim2v(t_ave,3,*nmax) - arrdim2v(tsq_ave,3,*nmax) +void +meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, + int* eflag_atom, double* eng_vdwl, double* eatom, int* ntype, + int* type, int* fmap, double* Arho1, double* Arho2, + double* Arho2b, double* Arho3, double* Arho3b, double* t_ave, + double* tsq_ave, double* Gamma, double* dGamma1, + double* dGamma2, double* dGamma3, double* rho, double* rho0, + double* rho1, double* rho2, double* rho3, double* fp, + int* errorflag) +{ - int i, elti; - int m; - double rhob, G, dG, Gbar, dGbar, gam, shp[3+1], Z; - double B, denom, rho_bkgd; + arrdim2v(Arho1, 3, *nmax); + arrdim2v(Arho2, 6, *nmax); + arrdim2v(Arho3, 10, *nmax); + arrdim2v(Arho3b, 3, *nmax); + arrdim2v(t_ave, 3, *nmax); + arrdim2v(tsq_ave, 3, *nmax); -// Complete the calculation of density + int i, elti; + int m; + double rhob, G, dG, Gbar, dGbar, gam, shp[3 + 1], Z; + double B, denom, rho_bkgd; - for(i=1; i<=*nlocal; i++) { - elti = arr1v(fmap,arr1v(type,i)); - if (elti > 0) { - arr1v(rho1,i) = 0.0; - arr1v(rho2,i) = -1.0/3.0*arr1v(Arho2b,i)*arr1v(Arho2b,i); - arr1v(rho3,i) = 0.0; - for(m=1; m<=3; m++) { - arr1v(rho1,i) = arr1v(rho1,i) + arr2v(Arho1,m,i)*arr2v(Arho1,m,i); - arr1v(rho3,i) = arr1v(rho3,i) - 3.0/5.0*arr2v(Arho3b,m,i)*arr2v(Arho3b,m,i); - } - for(m=1; m<=6; m++) { - arr1v(rho2,i) = arr1v(rho2,i) + meam_data.v2D[m]*arr2v(Arho2,m,i)*arr2v(Arho2,m,i); - } - for(m=1; m<=10; m++) { - arr1v(rho3,i) = arr1v(rho3,i) + meam_data.v3D[m]*arr2v(Arho3,m,i)*arr2v(Arho3,m,i); - } + // Complete the calculation of density - if( arr1v(rho0,i) > 0.0 ) { - if (meam_data.ialloy==1) { - arr2v(t_ave,1,i) = arr2v(t_ave,1,i)/arr2v(tsq_ave,1,i); - arr2v(t_ave,2,i) = arr2v(t_ave,2,i)/arr2v(tsq_ave,2,i); - arr2v(t_ave,3,i) = arr2v(t_ave,3,i)/arr2v(tsq_ave,3,i); - } else if (meam_data.ialloy==2) { - arr2v(t_ave,1,i) = meam_data.t1_meam[elti]; - arr2v(t_ave,2,i) = meam_data.t2_meam[elti]; - arr2v(t_ave,3,i) = meam_data.t3_meam[elti]; - } else { - arr2v(t_ave,1,i) = arr2v(t_ave,1,i)/arr1v(rho0,i); - arr2v(t_ave,2,i) = arr2v(t_ave,2,i)/arr1v(rho0,i); - arr2v(t_ave,3,i) = arr2v(t_ave,3,i)/arr1v(rho0,i); - } - } + for (i = 1; i <= *nlocal; i++) { + elti = arr1v(fmap, arr1v(type, i)); + if (elti > 0) { + arr1v(rho1, i) = 0.0; + arr1v(rho2, i) = -1.0 / 3.0 * arr1v(Arho2b, i) * arr1v(Arho2b, i); + arr1v(rho3, i) = 0.0; + for (m = 1; m <= 3; m++) { + arr1v(rho1, i) = + arr1v(rho1, i) + arr2v(Arho1, m, i) * arr2v(Arho1, m, i); + arr1v(rho3, i) = arr1v(rho3, i) - + 3.0 / 5.0 * arr2v(Arho3b, m, i) * arr2v(Arho3b, m, i); + } + for (m = 1; m <= 6; m++) { + arr1v(rho2, i) = + arr1v(rho2, i) + + meam_data.v2D[m] * arr2v(Arho2, m, i) * arr2v(Arho2, m, i); + } + for (m = 1; m <= 10; m++) { + arr1v(rho3, i) = + arr1v(rho3, i) + + meam_data.v3D[m] * arr2v(Arho3, m, i) * arr2v(Arho3, m, i); + } - arr1v(Gamma,i) = arr2v(t_ave,1,i)*arr1v(rho1,i) + arr2v(t_ave,2,i)*arr1v(rho2,i) + arr2v(t_ave,3,i)*arr1v(rho3,i); - - if( arr1v(rho0,i) > 0.0 ) { - arr1v(Gamma,i) = arr1v(Gamma,i)/(arr1v(rho0,i)*arr1v(rho0,i)); - } - - Z = meam_data.Z_meam[elti]; - - G_gam(arr1v(Gamma,i),meam_data.ibar_meam[elti], meam_data.gsmooth_factor,&G,errorflag); - if (*errorflag!=0) return; - get_shpfcn(shp,meam_data.lattce_meam[elti][elti]); - if (meam_data.ibar_meam[elti]<=0) { - Gbar = 1.0; - dGbar = 0.0; - } else { - if (meam_data.mix_ref_t==1) { - gam = (arr2v(t_ave,1,i)*shp[1]+arr2v(t_ave,2,i)*shp[2] +arr2v(t_ave,3,i)*shp[3])/(Z*Z); - } else { - gam = (meam_data.t1_meam[elti]*shp[1]+meam_data.t2_meam[elti]*shp[2] +meam_data.t3_meam[elti]*shp[3])/(Z*Z); - } - G_gam(gam,meam_data.ibar_meam[elti],meam_data.gsmooth_factor,&Gbar,errorflag); - } - arr1v(rho,i) = arr1v(rho0,i) * G; - - if (meam_data.mix_ref_t==1) { - if (meam_data.ibar_meam[elti]<=0) { - Gbar = 1.0; - dGbar = 0.0; - } else { - gam = (arr2v(t_ave,1,i)*shp[1]+arr2v(t_ave,2,i)*shp[2] +arr2v(t_ave,3,i)*shp[3])/(Z*Z); - dG_gam(gam,meam_data.ibar_meam[elti],meam_data.gsmooth_factor, &Gbar,&dGbar); - } - rho_bkgd = meam_data.rho0_meam[elti]*Z*Gbar; - } else { - if (meam_data.bkgd_dyn==1) { - rho_bkgd = meam_data.rho0_meam[elti]*Z; - } else { - rho_bkgd = meam_data.rho_ref_meam[elti]; - } - } - rhob = arr1v(rho,i)/rho_bkgd; - denom = 1.0/rho_bkgd; - - dG_gam(arr1v(Gamma,i),meam_data.ibar_meam[elti],meam_data.gsmooth_factor,&G,&dG); - - arr1v(dGamma1,i) = (G - 2*dG*arr1v(Gamma,i))*denom; - - if( !iszero(arr1v(rho0,i)) ) { - arr1v(dGamma2,i) = (dG/arr1v(rho0,i))*denom; - } else { - arr1v(dGamma2,i) = 0.0; - } - -// dGamma3 is nonzero only if we are using the "mixed" rule for -// computing t in the reference system (which is not correct, but -// included for backward compatibility - if (meam_data.mix_ref_t==1) { - arr1v(dGamma3,i) = arr1v(rho0,i)*G*dGbar/(Gbar*Z*Z)*denom; - } else { - arr1v(dGamma3,i) = 0.0; - } - - B = meam_data.A_meam[elti]*meam_data.Ec_meam[elti][elti]; - - if( !iszero(rhob) ) { - if (meam_data.emb_lin_neg==1 && rhob<=0) { - arr1v(fp,i) = -B; - } else { - arr1v(fp,i) = B*(log(rhob)+1.0); - } - if (*eflag_either!=0) { - if (*eflag_global!=0) { - if (meam_data.emb_lin_neg==1 && rhob<=0) { - *eng_vdwl = *eng_vdwl - B*rhob; - } else { - *eng_vdwl = *eng_vdwl + B*rhob*log(rhob); - } - } - if (*eflag_atom!=0) { - if (meam_data.emb_lin_neg==1 && rhob<=0) { - arr1v(eatom,i) = arr1v(eatom,i) - B*rhob; - } else { - arr1v(eatom,i) = arr1v(eatom,i) + B*rhob*log(rhob); - } - } - } - } else { - if (meam_data.emb_lin_neg==1) { - arr1v(fp,i) = -B; - } else { - arr1v(fp,i) = B; - } - } + if (arr1v(rho0, i) > 0.0) { + if (meam_data.ialloy == 1) { + arr2v(t_ave, 1, i) = arr2v(t_ave, 1, i) / arr2v(tsq_ave, 1, i); + arr2v(t_ave, 2, i) = arr2v(t_ave, 2, i) / arr2v(tsq_ave, 2, i); + arr2v(t_ave, 3, i) = arr2v(t_ave, 3, i) / arr2v(tsq_ave, 3, i); + } else if (meam_data.ialloy == 2) { + arr2v(t_ave, 1, i) = meam_data.t1_meam[elti]; + arr2v(t_ave, 2, i) = meam_data.t2_meam[elti]; + arr2v(t_ave, 3, i) = meam_data.t3_meam[elti]; + } else { + arr2v(t_ave, 1, i) = arr2v(t_ave, 1, i) / arr1v(rho0, i); + arr2v(t_ave, 2, i) = arr2v(t_ave, 2, i) / arr1v(rho0, i); + arr2v(t_ave, 3, i) = arr2v(t_ave, 3, i) / arr1v(rho0, i); } } + arr1v(Gamma, i) = arr2v(t_ave, 1, i) * arr1v(rho1, i) + + arr2v(t_ave, 2, i) * arr1v(rho2, i) + + arr2v(t_ave, 3, i) * arr1v(rho3, i); + + if (arr1v(rho0, i) > 0.0) { + arr1v(Gamma, i) = arr1v(Gamma, i) / (arr1v(rho0, i) * arr1v(rho0, i)); } -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + Z = meam_data.Z_meam[elti]; - void G_gam(double Gamma, int ibar, double gsmooth_factor, double *G, int *errorflag) - { -// Compute G(Gamma) based on selection flag ibar: -// 0 => G = sqrt(1+Gamma) -// 1 => G = exp(Gamma/2) -// 2 => not implemented -// 3 => G = 2/(1+exp(-Gamma)) -// 4 => G = sqrt(1+Gamma) -// -5 => G = +-sqrt(abs(1+Gamma)) - - double gsmooth_switchpoint; - if (ibar==0 || ibar==4) { - gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1); - if (Gamma < gsmooth_switchpoint) { -// e.g. gsmooth_factor is 99, {: -// gsmooth_switchpoint = -0.99 -// G = 0.01*(-0.99/Gamma)**99 - *G = 1/(gsmooth_factor+1) * pow((gsmooth_switchpoint/Gamma),gsmooth_factor); - *G = sqrt(*G); + G_gam(arr1v(Gamma, i), meam_data.ibar_meam[elti], + meam_data.gsmooth_factor, &G, errorflag); + if (*errorflag != 0) + return; + get_shpfcn(shp, meam_data.lattce_meam[elti][elti]); + if (meam_data.ibar_meam[elti] <= 0) { + Gbar = 1.0; + dGbar = 0.0; + } else { + if (meam_data.mix_ref_t == 1) { + gam = (arr2v(t_ave, 1, i) * shp[1] + arr2v(t_ave, 2, i) * shp[2] + + arr2v(t_ave, 3, i) * shp[3]) / + (Z * Z); } else { - *G = sqrt(1.0+Gamma); + gam = (meam_data.t1_meam[elti] * shp[1] + + meam_data.t2_meam[elti] * shp[2] + + meam_data.t3_meam[elti] * shp[3]) / + (Z * Z); } - } else if (ibar==1) { - *G = fm_exp(Gamma/2.0); - } else if (ibar==3) { - *G = 2.0/(1.0+exp(-Gamma)); - } else if (ibar==-5) { - if ((1.0+Gamma)>=0) { - *G = sqrt(1.0+Gamma); + G_gam(gam, meam_data.ibar_meam[elti], meam_data.gsmooth_factor, &Gbar, + errorflag); + } + arr1v(rho, i) = arr1v(rho0, i) * G; + + if (meam_data.mix_ref_t == 1) { + if (meam_data.ibar_meam[elti] <= 0) { + Gbar = 1.0; + dGbar = 0.0; } else { - *G = -sqrt(-1.0-Gamma); + gam = (arr2v(t_ave, 1, i) * shp[1] + arr2v(t_ave, 2, i) * shp[2] + + arr2v(t_ave, 3, i) * shp[3]) / + (Z * Z); + dG_gam(gam, meam_data.ibar_meam[elti], meam_data.gsmooth_factor, + &Gbar, &dGbar); + } + rho_bkgd = meam_data.rho0_meam[elti] * Z * Gbar; + } else { + if (meam_data.bkgd_dyn == 1) { + rho_bkgd = meam_data.rho0_meam[elti] * Z; + } else { + rho_bkgd = meam_data.rho_ref_meam[elti]; + } + } + rhob = arr1v(rho, i) / rho_bkgd; + denom = 1.0 / rho_bkgd; + + dG_gam(arr1v(Gamma, i), meam_data.ibar_meam[elti], + meam_data.gsmooth_factor, &G, &dG); + + arr1v(dGamma1, i) = (G - 2 * dG * arr1v(Gamma, i)) * denom; + + if (!iszero(arr1v(rho0, i))) { + arr1v(dGamma2, i) = (dG / arr1v(rho0, i)) * denom; + } else { + arr1v(dGamma2, i) = 0.0; + } + + // dGamma3 is nonzero only if we are using the "mixed" rule for + // computing t in the reference system (which is not correct, but + // included for backward compatibility + if (meam_data.mix_ref_t == 1) { + arr1v(dGamma3, i) = arr1v(rho0, i) * G * dGbar / (Gbar * Z * Z) * denom; + } else { + arr1v(dGamma3, i) = 0.0; + } + + B = meam_data.A_meam[elti] * meam_data.Ec_meam[elti][elti]; + + if (!iszero(rhob)) { + if (meam_data.emb_lin_neg == 1 && rhob <= 0) { + arr1v(fp, i) = -B; + } else { + arr1v(fp, i) = B * (log(rhob) + 1.0); + } + if (*eflag_either != 0) { + if (*eflag_global != 0) { + if (meam_data.emb_lin_neg == 1 && rhob <= 0) { + *eng_vdwl = *eng_vdwl - B * rhob; + } else { + *eng_vdwl = *eng_vdwl + B * rhob * log(rhob); + } + } + if (*eflag_atom != 0) { + if (meam_data.emb_lin_neg == 1 && rhob <= 0) { + arr1v(eatom, i) = arr1v(eatom, i) - B * rhob; + } else { + arr1v(eatom, i) = arr1v(eatom, i) + B * rhob * log(rhob); + } + } } } else { - *errorflag = 1; - } - } - -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - void dG_gam(double Gamma, int ibar, double gsmooth_factor,double *G, double *dG) - { -// Compute G(Gamma) and dG(gamma) based on selection flag ibar: -// 0 => G = sqrt(1+Gamma) -// 1 => G = fm_exp(Gamma/2) -// 2 => not implemented -// 3 => G = 2/(1+fm_exp(-Gamma)) -// 4 => G = sqrt(1+Gamma) -// -5 => G = +-sqrt(abs(1+Gamma)) - double gsmooth_switchpoint; - - if (ibar==0 || ibar==4) { - gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1); - if (Gamma < gsmooth_switchpoint) { -// e.g. gsmooth_factor is 99, {: -// gsmooth_switchpoint = -0.99 -// G = 0.01*(-0.99/Gamma)**99 - *G = 1/(gsmooth_factor+1) * pow((gsmooth_switchpoint/Gamma), gsmooth_factor); - *G = sqrt(*G); - *dG = -gsmooth_factor* *G/(2.0*Gamma); + if (meam_data.emb_lin_neg == 1) { + arr1v(fp, i) = -B; } else { - *G = sqrt(1.0+Gamma); - *dG = 1.0/(2.0* *G); - } - } else if (ibar==1) { - *G = fm_exp(Gamma/2.0); - *dG = *G/2.0; - } else if (ibar==3) { - *G = 2.0/(1.0+fm_exp(-Gamma)); - *dG = *G*(2.0-*G)/2; - } else if (ibar==-5) { - if ((1.0+Gamma)>=0) { - *G = sqrt(1.0+Gamma); - *dG = 1.0/(2.0* *G); - } else { - *G = -sqrt(-1.0-Gamma); - *dG = -1.0/(2.0* *G); + arr1v(fp, i) = B; } } - } + } + } +} -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +void +G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, int* errorflag) +{ + // Compute G(Gamma) based on selection flag ibar: + // 0 => G = sqrt(1+Gamma) + // 1 => G = exp(Gamma/2) + // 2 => not implemented + // 3 => G = 2/(1+exp(-Gamma)) + // 4 => G = sqrt(1+Gamma) + // -5 => G = +-sqrt(abs(1+Gamma)) + + double gsmooth_switchpoint; + if (ibar == 0 || ibar == 4) { + gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1); + if (Gamma < gsmooth_switchpoint) { + // e.g. gsmooth_factor is 99, {: + // gsmooth_switchpoint = -0.99 + // G = 0.01*(-0.99/Gamma)**99 + *G = 1 / (gsmooth_factor + 1) * + pow((gsmooth_switchpoint / Gamma), gsmooth_factor); + *G = sqrt(*G); + } else { + *G = sqrt(1.0 + Gamma); + } + } else if (ibar == 1) { + *G = fm_exp(Gamma / 2.0); + } else if (ibar == 3) { + *G = 2.0 / (1.0 + exp(-Gamma)); + } else if (ibar == -5) { + if ((1.0 + Gamma) >= 0) { + *G = sqrt(1.0 + Gamma); + } else { + *G = -sqrt(-1.0 - Gamma); + } + } else { + *errorflag = 1; + } +} + +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +void +dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, double* dG) +{ + // Compute G(Gamma) and dG(gamma) based on selection flag ibar: + // 0 => G = sqrt(1+Gamma) + // 1 => G = fm_exp(Gamma/2) + // 2 => not implemented + // 3 => G = 2/(1+fm_exp(-Gamma)) + // 4 => G = sqrt(1+Gamma) + // -5 => G = +-sqrt(abs(1+Gamma)) + double gsmooth_switchpoint; + + if (ibar == 0 || ibar == 4) { + gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1); + if (Gamma < gsmooth_switchpoint) { + // e.g. gsmooth_factor is 99, {: + // gsmooth_switchpoint = -0.99 + // G = 0.01*(-0.99/Gamma)**99 + *G = 1 / (gsmooth_factor + 1) * + pow((gsmooth_switchpoint / Gamma), gsmooth_factor); + *G = sqrt(*G); + *dG = -gsmooth_factor * *G / (2.0 * Gamma); + } else { + *G = sqrt(1.0 + Gamma); + *dG = 1.0 / (2.0 * *G); + } + } else if (ibar == 1) { + *G = fm_exp(Gamma / 2.0); + *dG = *G / 2.0; + } else if (ibar == 3) { + *G = 2.0 / (1.0 + fm_exp(-Gamma)); + *dG = *G * (2.0 - *G) / 2; + } else if (ibar == -5) { + if ((1.0 + Gamma) >= 0) { + *G = sqrt(1.0 + Gamma); + *dG = 1.0 / (2.0 * *G); + } else { + *G = -sqrt(-1.0 - Gamma); + *dG = -1.0 / (2.0 * *G); + } + } +} + +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc } \ No newline at end of file diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp index 8a17f4e7e5..4ea56ded9c 100644 --- a/src/USER-MEAMC/meam_dens_init.cpp +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -2,25 +2,30 @@ extern "C" { #include "meam.h" #include - void getscreen(int i, int nmax, double *scrfcn, double *dscrfcn, double *fcpair, double *x, - int numneigh, int *firstneigh, - int numneigh_full, int *firstneigh_full, - int ntype, int *type, int *fmap); - void calc_rho1(int i, int nmax, int ntype, int *type, int *fmap, double *x, - int numneigh, int *firstneigh, - double *scrfcn, double *fcpair, double *rho0, double *arho1, double *arho2, double *arho2b, - double *arho3, double *arho3b, double *t_ave, double *tsq_ave); - - void screen(int i, int j, int nmax, double *x, double rijsq, double *sij, int numneigh_full, int *firstneigh_full, int ntype, int *type, int *fmap); - void dsij(int i,int j,int k,int jn,int nmax,int numneigh,double rij2,double *dsij1,double *dsij2, int ntype, int *type, int *fmap,double *x,double *scrfcn, double*fcpair); - void fcut(double xi, double * fc); - void dfcut(double xi, double *fc, double *dfc); - void dCfunc(double rij2,double rik2,double rjk2,double *dCikj); - void dCfunc2(double rij2,double rik2,double rjk2,double *dCikj1,double *dCikj2); +void getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair, + double* x, int numneigh, int* firstneigh, int numneigh_full, + int* firstneigh_full, int ntype, int* type, int* fmap); +void calc_rho1(int i, int nmax, int ntype, int* type, int* fmap, double* x, + int numneigh, int* firstneigh, double* scrfcn, double* fcpair, + double* rho0, double* arho1, double* arho2, double* arho2b, + double* arho3, double* arho3b, double* t_ave, double* tsq_ave); + +void screen(int i, int j, int nmax, double* x, double rijsq, double* sij, + int numneigh_full, int* firstneigh_full, int ntype, int* type, + int* fmap); +void dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2, + double* dsij1, double* dsij2, int ntype, int* type, int* fmap, + double* x, double* scrfcn, double* fcpair); +void fcut(double xi, double* fc); +void dfcut(double xi, double* fc, double* dfc); +void dCfunc(double rij2, double rik2, double rjk2, double* dCikj); +void dCfunc2(double rij2, double rik2, double rjk2, double* dCikj1, + double* dCikj2); // Extern "C" declaration has the form: // -// void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *, +// void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, +// double *, // int *, int *, int *, int *, // double *, double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, int *); @@ -35,470 +40,505 @@ extern "C" { // &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag); // - void meam_dens_init_(int *i, int *nmax, - int *ntype, int *type, int *fmap, double *x, - int *numneigh, int *firstneigh, - int *numneigh_full, int *firstneigh_full, - double *scrfcn, double *dscrfcn, double *fcpair, double *rho0, double *arho1, double *arho2, double *arho2b, - double *arho3, double *arho3b, double *t_ave, double *tsq_ave, int *errorflag) - { - *errorflag = 0; +void +meam_dens_init_(int* i, int* nmax, int* ntype, int* type, int* fmap, double* x, + int* numneigh, int* firstneigh, int* numneigh_full, + int* firstneigh_full, double* scrfcn, double* dscrfcn, + double* fcpair, double* rho0, double* arho1, double* arho2, + double* arho2b, double* arho3, double* arho3b, double* t_ave, + double* tsq_ave, int* errorflag) +{ + *errorflag = 0; -// Compute screening function and derivatives - getscreen(*i, *nmax, scrfcn, dscrfcn, fcpair, x, - *numneigh, firstneigh, - *numneigh_full, firstneigh_full, - *ntype, type, fmap); + // Compute screening function and derivatives + getscreen(*i, *nmax, scrfcn, dscrfcn, fcpair, x, *numneigh, firstneigh, + *numneigh_full, firstneigh_full, *ntype, type, fmap); -// Calculate intermediate density terms to be communicated - calc_rho1(*i, *nmax, *ntype, type, fmap, x, - *numneigh, firstneigh, - scrfcn, fcpair, rho0, arho1, arho2, arho2b, - arho3, arho3b, t_ave, tsq_ave); + // Calculate intermediate density terms to be communicated + calc_rho1(*i, *nmax, *ntype, type, fmap, x, *numneigh, firstneigh, scrfcn, + fcpair, rho0, arho1, arho2, arho2b, arho3, arho3b, t_ave, tsq_ave); +} - } +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +void +getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair, + double* x, int numneigh, int* firstneigh, int numneigh_full, + int* firstneigh_full, int ntype, int* type, int* fmap) +{ - void getscreen(int i, int nmax, double *scrfcn, double *dscrfcn, double *fcpair, double *x, - int numneigh, int *firstneigh, - int numneigh_full, int *firstneigh_full, - int ntype, int *type, int *fmap) - { - - arrdim2v(x,3,nmax) + arrdim2v(x, 3, nmax); - int jn,j,kn,k; - int elti,eltj,eltk; - double xitmp,yitmp,zitmp,delxij,delyij,delzij,rij2,rij; - double xjtmp,yjtmp,zjtmp,delxik,delyik,delzik,rik2/*,rik*/; - double xktmp,yktmp,zktmp,delxjk,delyjk,delzjk,rjk2/*,rjk*/; - double xik,xjk,sij,fcij,sfcij,dfcij,sikj,dfikj,cikj; - double Cmin,Cmax,delc,/*ebound,*/rbound,a,coef1,coef2; - //double coef1a,coef1b,coef2a,coef2b; - double dCikj; - /*unused:double dC1a,dC1b,dC2a,dC2b;*/ - double rnorm,fc,dfc,drinv; + int jn, j, kn, k; + int elti, eltj, eltk; + double xitmp, yitmp, zitmp, delxij, delyij, delzij, rij2, rij; + double xjtmp, yjtmp, zjtmp, delxik, delyik, delzik, rik2 /*,rik*/; + double xktmp, yktmp, zktmp, delxjk, delyjk, delzjk, rjk2 /*,rjk*/; + double xik, xjk, sij, fcij, sfcij, dfcij, sikj, dfikj, cikj; + double Cmin, Cmax, delc, /*ebound,*/ rbound, a, coef1, coef2; + // double coef1a,coef1b,coef2a,coef2b; + double dCikj; + /*unused:double dC1a,dC1b,dC2a,dC2b;*/ + double rnorm, fc, dfc, drinv; - drinv = 1.0/meam_data.delr_meam; - elti = arr1v(fmap,arr1v(type,i)); + drinv = 1.0 / meam_data.delr_meam; + elti = arr1v(fmap, arr1v(type, i)); - if (elti > 0) { + if (elti > 0) { - xitmp = arr2v(x,1,i); - yitmp = arr2v(x,2,i); - zitmp = arr2v(x,3,i); + xitmp = arr2v(x, 1, i); + yitmp = arr2v(x, 2, i); + zitmp = arr2v(x, 3, i); - for(jn=1; jn<=numneigh; jn++) { - j = arr1v(firstneigh,jn); + for (jn = 1; jn <= numneigh; jn++) { + j = arr1v(firstneigh, jn); - eltj = arr1v(fmap,arr1v(type,j)); - if (eltj > 0) { - -// First compute screening function itself, sij - xjtmp = arr2v(x,1,j); - yjtmp = arr2v(x,2,j); - zjtmp = arr2v(x,3,j); - delxij = xjtmp - xitmp; - delyij = yjtmp - yitmp; - delzij = zjtmp - zitmp; - rij2 = delxij*delxij + delyij*delyij + delzij*delzij; - rij = sqrt(rij2); - if (rij > meam_data.rc_meam) { - fcij = 0.0; - dfcij = 0.0; - sij = 0.0; - } else { - rnorm = (meam_data.rc_meam-rij)*drinv; - screen(i, j, nmax, x, rij2, &sij, numneigh_full, firstneigh_full, ntype, type, fmap); - dfcut(rnorm,&fc,&dfc); - fcij = fc; - dfcij = dfc*drinv; - } - -// Now compute derivatives - arr1v(dscrfcn,jn) = 0.0; - sfcij = sij*fcij; - if (iszero(sfcij) || iszero(sfcij-1.0)) goto LABEL_100; - rbound = meam_data.ebound_meam[elti][eltj] * rij2; - for(kn=1; kn<=numneigh_full; kn++) { - k = arr1v(firstneigh_full,kn); - if (k==j) continue; - eltk = arr1v(fmap,arr1v(type,k)); - if (eltk==0) continue; - xktmp = arr2v(x,1,k); - yktmp = arr2v(x,2,k); - zktmp = arr2v(x,3,k); - delxjk = xktmp - xjtmp; - delyjk = yktmp - yjtmp; - delzjk = zktmp - zjtmp; - rjk2 = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk; - if (rjk2 > rbound) continue; - delxik = xktmp - xitmp; - delyik = yktmp - yitmp; - delzik = zktmp - zitmp; - rik2 = delxik*delxik + delyik*delyik + delzik*delzik; - if (rik2 > rbound) continue; - xik = rik2/rij2; - xjk = rjk2/rij2; - a = 1 - (xik-xjk)*(xik-xjk); -// if a < 0, then ellipse equation doesn't describe this case and -// atom k can't possibly screen i-j - if (a<=0.0) continue; - cikj = (2.0*(xik+xjk) + a - 2.0)/a; - Cmax = meam_data.Cmax_meam[elti][eltj][eltk]; - Cmin = meam_data.Cmin_meam[elti][eltj][eltk]; - if (cikj>=Cmax) { - continue; -// Note that cikj may be slightly negative (within numerical -// tolerance) if atoms are colinear, so don't reject that case here -// (other negative cikj cases were handled by the test on "a" above) -// Note that we never have 0 0) { + // First compute screening function itself, sij + xjtmp = arr2v(x, 1, j); + yjtmp = arr2v(x, 2, j); + zjtmp = arr2v(x, 3, j); + delxij = xjtmp - xitmp; + delyij = yjtmp - yitmp; + delzij = zjtmp - zitmp; + rij2 = delxij * delxij + delyij * delyij + delzij * delzij; + rij = sqrt(rij2); + if (rij > meam_data.rc_meam) { + fcij = 0.0; + dfcij = 0.0; + sij = 0.0; + } else { + rnorm = (meam_data.rc_meam - rij) * drinv; + screen(i, j, nmax, x, rij2, &sij, numneigh_full, firstneigh_full, + ntype, type, fmap); + dfcut(rnorm, &fc, &dfc); + fcij = fc; + dfcij = dfc * drinv; } + // Now compute derivatives + arr1v(dscrfcn, jn) = 0.0; + sfcij = sij * fcij; + if (iszero(sfcij) || iszero(sfcij - 1.0)) + goto LABEL_100; + rbound = meam_data.ebound_meam[elti][eltj] * rij2; + for (kn = 1; kn <= numneigh_full; kn++) { + k = arr1v(firstneigh_full, kn); + if (k == j) + continue; + eltk = arr1v(fmap, arr1v(type, k)); + if (eltk == 0) + continue; + xktmp = arr2v(x, 1, k); + yktmp = arr2v(x, 2, k); + zktmp = arr2v(x, 3, k); + delxjk = xktmp - xjtmp; + delyjk = yktmp - yjtmp; + delzjk = zktmp - zjtmp; + rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk; + if (rjk2 > rbound) + continue; + delxik = xktmp - xitmp; + delyik = yktmp - yitmp; + delzik = zktmp - zitmp; + rik2 = delxik * delxik + delyik * delyik + delzik * delzik; + if (rik2 > rbound) + continue; + xik = rik2 / rij2; + xjk = rjk2 / rij2; + a = 1 - (xik - xjk) * (xik - xjk); + // if a < 0, then ellipse equation doesn't describe this case and + // atom k can't possibly screen i-j + if (a <= 0.0) + continue; + cikj = (2.0 * (xik + xjk) + a - 2.0) / a; + Cmax = meam_data.Cmax_meam[elti][eltj][eltk]; + Cmin = meam_data.Cmin_meam[elti][eltj][eltk]; + if (cikj >= Cmax) { + continue; + // Note that cikj may be slightly negative (within numerical + // tolerance) if atoms are colinear, so don't reject that case + // here + // (other negative cikj cases were handled by the test on "a" + // above) + // Note that we never have 0 ebound*rijsq, atom k is definitely outside the ellipse - rbound = meam_data.ebound_meam[elti][eltj]*rijsq; + *sij = 1.0; + eltj = arr1v(fmap, arr1v(type, j)); + elti = arr1v(fmap, arr1v(type, j)); - for(nk=1; nk<=numneigh_full; nk++) { - k = arr1v(firstneigh_full,nk); - eltk = arr1v(fmap,arr1v(type,k)); - if (k==j) continue; - delxjk = arr2v(x,1,k) - arr2v(x,1,j); - delyjk = arr2v(x,2,k) - arr2v(x,2,j); - delzjk = arr2v(x,3,k) - arr2v(x,3,j); - rjksq = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk; - if (rjksq > rbound) continue; - delxik = arr2v(x,1,k) - arr2v(x,1,i); - delyik = arr2v(x,2,k) - arr2v(x,2,i); - delzik = arr2v(x,3,k) - arr2v(x,3,i); - riksq = delxik*delxik + delyik*delyik + delzik*delzik; - if (riksq > rbound) continue; - xik = riksq/rijsq; - xjk = rjksq/rijsq; - a = 1 - (xik-xjk)*(xik-xjk); -// if a < 0, then ellipse equation doesn't describe this case and -// atom k can't possibly screen i-j - if (a<=0.0) continue; - cikj = (2.0*(xik+xjk) + a - 2.0)/a; - Cmax = meam_data.Cmax_meam[elti][eltj][eltk]; - Cmin = meam_data.Cmin_meam[elti][eltj][eltk]; - if (cikj>=Cmax) - continue; -// note that cikj may be slightly negative (within numerical -// tolerance) if atoms are colinear, so don't reject that case here -// (other negative cikj cases were handled by the test on "a" above) - else if (cikj<=Cmin) { - *sij = 0.0; - break; - } else { - delc = Cmax - Cmin; - cikj = (cikj-Cmin)/delc; - fcut(cikj,&sikj); - } - *sij = *sij * sikj; - } + // if rjksq > ebound*rijsq, atom k is definitely outside the ellipse + rbound = meam_data.ebound_meam[elti][eltj] * rijsq; - } + for (nk = 1; nk <= numneigh_full; nk++) { + k = arr1v(firstneigh_full, nk); + eltk = arr1v(fmap, arr1v(type, k)); + if (k == j) + continue; + delxjk = arr2v(x, 1, k) - arr2v(x, 1, j); + delyjk = arr2v(x, 2, k) - arr2v(x, 2, j); + delzjk = arr2v(x, 3, k) - arr2v(x, 3, j); + rjksq = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk; + if (rjksq > rbound) + continue; + delxik = arr2v(x, 1, k) - arr2v(x, 1, i); + delyik = arr2v(x, 2, k) - arr2v(x, 2, i); + delzik = arr2v(x, 3, k) - arr2v(x, 3, i); + riksq = delxik * delxik + delyik * delyik + delzik * delzik; + if (riksq > rbound) + continue; + xik = riksq / rijsq; + xjk = rjksq / rijsq; + a = 1 - (xik - xjk) * (xik - xjk); + // if a < 0, then ellipse equation doesn't describe this case and + // atom k can't possibly screen i-j + if (a <= 0.0) + continue; + cikj = (2.0 * (xik + xjk) + a - 2.0) / a; + Cmax = meam_data.Cmax_meam[elti][eltj][eltk]; + Cmin = meam_data.Cmin_meam[elti][eltj][eltk]; + if (cikj >= Cmax) + continue; + // note that cikj may be slightly negative (within numerical + // tolerance) if atoms are colinear, so don't reject that case here + // (other negative cikj cases were handled by the test on "a" above) + else if (cikj <= Cmin) { + *sij = 0.0; + break; + } else { + delc = Cmax - Cmin; + cikj = (cikj - Cmin) / delc; + fcut(cikj, &sikj); + } + *sij = *sij * sikj; + } +} -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - void dsij(int i,int j,int k,int jn,int nmax,int numneigh,double rij2,double *dsij1,double *dsij2, int ntype, int *type, int *fmap,double *x,double *scrfcn, double*fcpair) - { -// Inputs: i,j,k = id's of 3 atom triplet -// jn = id of i-j pair -// rij2 = squared distance between i and j -// Outputs: dsij1 = deriv. of sij w.r.t. rik -// dsij2 = deriv. of sij w.r.t. rjk +void +dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2, + double* dsij1, double* dsij2, int ntype, int* type, int* fmap, double* x, + double* scrfcn, double* fcpair) +{ + // Inputs: i,j,k = id's of 3 atom triplet + // jn = id of i-j pair + // rij2 = squared distance between i and j + // Outputs: dsij1 = deriv. of sij w.r.t. rik + // dsij2 = deriv. of sij w.r.t. rjk - arrdim2v(x,3,nmax) - int elti,eltj,eltk; - double rik2,rjk2; + arrdim2v(x, 3, nmax) int elti, eltj, eltk; + double rik2, rjk2; - double dxik,dyik,dzik; - double dxjk,dyjk,dzjk; - double rbound,delc,sij,xik,xjk,cikj,sikj,dfc,a; - double Cmax,Cmin,dCikj1,dCikj2; + double dxik, dyik, dzik; + double dxjk, dyjk, dzjk; + double rbound, delc, sij, xik, xjk, cikj, sikj, dfc, a; + double Cmax, Cmin, dCikj1, dCikj2; - sij = arr1v(scrfcn,jn)*arr1v(fcpair,jn); - elti = arr1v(fmap,arr1v(type,i)); - eltj = arr1v(fmap,arr1v(type,j)); - eltk = arr1v(fmap,arr1v(type,k)); - Cmax = meam_data.Cmax_meam[elti][eltj][eltk]; - Cmin = meam_data.Cmin_meam[elti][eltj][eltk]; + sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn); + elti = arr1v(fmap, arr1v(type, i)); + eltj = arr1v(fmap, arr1v(type, j)); + eltk = arr1v(fmap, arr1v(type, k)); + Cmax = meam_data.Cmax_meam[elti][eltj][eltk]; + Cmin = meam_data.Cmin_meam[elti][eltj][eltk]; - *dsij1 = 0.0; - *dsij2 = 0.0; - if (!iszero(sij) && !iszero(sij-1.0)) { - rbound = rij2*meam_data.ebound_meam[elti][eltj]; - delc = Cmax-Cmin; - dxjk = arr2v(x,1,k) - arr2v(x,1,j); - dyjk = arr2v(x,2,k) - arr2v(x,2,j); - dzjk = arr2v(x,3,k) - arr2v(x,3,j); - rjk2 = dxjk*dxjk + dyjk*dyjk + dzjk*dzjk; - if (rjk2<=rbound) { - dxik = arr2v(x,1,k) - arr2v(x,1,i); - dyik = arr2v(x,2,k) - arr2v(x,2,i); - dzik = arr2v(x,3,k) - arr2v(x,3,i); - rik2 = dxik*dxik + dyik*dyik + dzik*dzik; - if (rik2<=rbound) { - xik = rik2/rij2; - xjk = rjk2/rij2; - a = 1 - (xik-xjk)*(xik-xjk); - if (!iszero(a)) { - cikj = (2.0*(xik+xjk) + a - 2.0)/a; - if (cikj>=Cmin && cikj<=Cmax) { - cikj = (cikj-Cmin)/delc; - dfcut(cikj,&sikj,&dfc); - dCfunc2(rij2,rik2,rjk2,&dCikj1,&dCikj2); - a = sij/delc*dfc/sikj; - *dsij1 = a*dCikj1; - *dsij2 = a*dCikj2; - } - } + *dsij1 = 0.0; + *dsij2 = 0.0; + if (!iszero(sij) && !iszero(sij - 1.0)) { + rbound = rij2 * meam_data.ebound_meam[elti][eltj]; + delc = Cmax - Cmin; + dxjk = arr2v(x, 1, k) - arr2v(x, 1, j); + dyjk = arr2v(x, 2, k) - arr2v(x, 2, j); + dzjk = arr2v(x, 3, k) - arr2v(x, 3, j); + rjk2 = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk; + if (rjk2 <= rbound) { + dxik = arr2v(x, 1, k) - arr2v(x, 1, i); + dyik = arr2v(x, 2, k) - arr2v(x, 2, i); + dzik = arr2v(x, 3, k) - arr2v(x, 3, i); + rik2 = dxik * dxik + dyik * dyik + dzik * dzik; + if (rik2 <= rbound) { + xik = rik2 / rij2; + xjk = rjk2 / rij2; + a = 1 - (xik - xjk) * (xik - xjk); + if (!iszero(a)) { + cikj = (2.0 * (xik + xjk) + a - 2.0) / a; + if (cikj >= Cmin && cikj <= Cmax) { + cikj = (cikj - Cmin) / delc; + dfcut(cikj, &sikj, &dfc); + dCfunc2(rij2, rik2, rjk2, &dCikj1, &dCikj2); + a = sij / delc * dfc / sikj; + *dsij1 = a * dCikj1; + *dsij2 = a * dCikj2; } } } - } + } + } +} +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +void +fcut(double xi, double* fc) +{ + // cutoff function + double a; + if (xi >= 1.0) + *fc = 1.0; + else if (xi <= 0.0) + *fc = 0.0; + else { + a = 1.0 - xi; + a = a * a; + a = a * a; + a = 1.0 - a; + *fc = a * a; + // fc = xi + } +} - void fcut(double xi, double * fc) - { -// cutoff function - double a; - if (xi>=1.0) - *fc = 1.0; - else if (xi<=0.0) - *fc = 0.0; - else { - a = 1.0-xi; - a = a*a; - a = a*a; - a = 1.0-a; - *fc = a*a; -// fc = xi - } - } +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +void +dfcut(double xi, double* fc, double* dfc) +{ + // cutoff function and its derivative + double a, a3, a4; + if (xi >= 1.0) { + *fc = 1.0; + *dfc = 0.0; + } else if (xi <= 0.0) { + *fc = 0.0; + *dfc = 0.0; + } else { + a = 1.0 - xi; + a3 = a * a * a; + a4 = a * a3; + *fc = pow((1.0 - a4), 2); + *dfc = 8 * (1.0 - a4) * a3; + // fc = xi + // dfc = 1.d0 + } +} - void dfcut(double xi, double *fc, double *dfc) - { -// cutoff function and its derivative - double a,a3,a4; - if (xi>=1.0) { - *fc = 1.0; - *dfc = 0.0; - } else if (xi<=0.0) { - *fc = 0.0; - *dfc = 0.0; - } else { - a = 1.0-xi; - a3 = a*a*a; - a4 = a*a3; - *fc = pow((1.0-a4), 2); - *dfc = 8*(1.0-a4)*a3; -// fc = xi -// dfc = 1.d0 - } - } +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +void +dCfunc(double rij2, double rik2, double rjk2, double* dCikj) +{ + // Inputs: rij,rij2,rik2,rjk2 + // Outputs: dCikj = derivative of Cikj w.r.t. rij + double rij4, a, b, denom; - void dCfunc(double rij2,double rik2,double rjk2,double *dCikj) - { -// Inputs: rij,rij2,rik2,rjk2 -// Outputs: dCikj = derivative of Cikj w.r.t. rij - double rij4,a,b,denom; + rij4 = rij2 * rij2; + a = rik2 - rjk2; + b = rik2 + rjk2; + denom = rij4 - a * a; + denom = denom * denom; + *dCikj = -4 * (-2 * rij2 * a * a + rij4 * b + a * a * b) / denom; +} - rij4 = rij2*rij2; - a = rik2-rjk2; - b = rik2+rjk2; - denom = rij4 - a*a; - denom = denom*denom; - *dCikj = -4*(-2*rij2*a*a + rij4*b + a*a*b)/denom; - } +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +void +dCfunc2(double rij2, double rik2, double rjk2, double* dCikj1, double* dCikj2) +{ + // Inputs: rij,rij2,rik2,rjk2 + // Outputs: dCikj1 = derivative of Cikj w.r.t. rik + // dCikj2 = derivative of Cikj w.r.t. rjk + double rij4, rik4, rjk4, a, b, denom; - void dCfunc2(double rij2,double rik2,double rjk2,double *dCikj1,double *dCikj2) - { -// Inputs: rij,rij2,rik2,rjk2 -// Outputs: dCikj1 = derivative of Cikj w.r.t. rik -// dCikj2 = derivative of Cikj w.r.t. rjk - double rij4,rik4,rjk4,a,b,denom; - - rij4 = rij2*rij2; - rik4 = rik2*rik2; - rjk4 = rjk2*rjk2; - a = rik2-rjk2; - b = rik2+rjk2; - denom = rij4 - a*a; - denom = denom*denom; - *dCikj1 = 4*rij2*(rij4 + rik4 + 2*rik2*rjk2 - 3*rjk4 - 2*rij2*a)/ denom; - *dCikj2 = 4*rij2*(rij4 - 3*rik4 + 2*rik2*rjk2 + rjk4 + 2*rij2*a)/ denom; - - (void)(b); - } - -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + rij4 = rij2 * rij2; + rik4 = rik2 * rik2; + rjk4 = rjk2 * rjk2; + a = rik2 - rjk2; + b = rik2 + rjk2; + denom = rij4 - a * a; + denom = denom * denom; + *dCikj1 = 4 * rij2 * + (rij4 + rik4 + 2 * rik2 * rjk2 - 3 * rjk4 - 2 * rij2 * a) / denom; + *dCikj2 = 4 * rij2 * + (rij4 - 3 * rik4 + 2 * rik2 * rjk2 + rjk4 + 2 * rij2 * a) / denom; + (void)(b); +} +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc } \ No newline at end of file diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp index b10a766716..8f2efa2daa 100644 --- a/src/USER-MEAMC/meam_force.cpp +++ b/src/USER-MEAMC/meam_force.cpp @@ -1,21 +1,24 @@ -extern "C"{ +extern "C" { #include "meam.h" #include // in meam_setup_done -void get_shpfcn(double *s /* s(3) */, lattice_t latt); +void get_shpfcn(double* s /* s(3) */, lattice_t latt); // in meam_dens_init -void dsij(int i,int j,int k,int jn,int nmax,int numneigh,double rij2,double *dsij1,double *dsij2, int ntype, int *type, int *fmap,double *x,double *scrfcn, double*fcpair); - +void dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2, + double* dsij1, double* dsij2, int ntype, int* type, int* fmap, + double* x, double* scrfcn, double* fcpair); // Extern "C" declaration has the form: // -// void meam_force_(int *, int *, int *, double *, int *, int *, int *, double *, +// void meam_force_(int *, int *, int *, double *, int *, int *, int *, double +// *, // int *, int *, int *, int *, double *, double *, // double *, double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, double *, -// double *, double *, double *, double *, double *, double *, int *); +// double *, double *, double *, double *, double *, double *, int +//*); // // Call from pair_meam.cpp has the form: // @@ -28,527 +31,575 @@ void dsij(int i,int j,int k,int jn,int nmax,int numneigh,double rij2,double *dsi // &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag); // - void meam_force_(int *iptr, int *nmax, - int *eflag_either, int *eflag_global, int *eflag_atom, int *vflag_atom, - double *eng_vdwl, double *eatom, int *ntype, int *type, int *fmap, double *x, - int *numneigh, int *firstneigh, int *numneigh_full, int *firstneigh_full, - double *scrfcn, double *dscrfcn, double *fcpair, - double *dGamma1, double *dGamma2, double *dGamma3, double *rho0, double *rho1, double *rho2, double *rho3, double *fp, - double *Arho1, double *Arho2, double *Arho2b, double *Arho3, double *Arho3b, double *t_ave, double *tsq_ave, double *f, - double *vatom, int *errorflag) - { - arrdim2v(x,3,*nmax) - arrdim2v(Arho1,3,*nmax) - arrdim2v(Arho2,6,*nmax) - arrdim2v(Arho3,10,*nmax) - arrdim2v(Arho3b,3,*nmax) - arrdim2v(t_ave,3,*nmax) - arrdim2v(tsq_ave,3,*nmax) - arrdim2v(f,3,*nmax) - arrdim2v(vatom,6,*nmax) +void +meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, + int* eflag_atom, int* vflag_atom, double* eng_vdwl, double* eatom, + int* ntype, int* type, int* fmap, double* x, int* numneigh, + int* firstneigh, int* numneigh_full, int* firstneigh_full, + double* scrfcn, double* dscrfcn, double* fcpair, double* dGamma1, + double* dGamma2, double* dGamma3, double* rho0, double* rho1, + double* rho2, double* rho3, double* fp, double* Arho1, + double* Arho2, double* Arho2b, double* Arho3, double* Arho3b, + double* t_ave, double* tsq_ave, double* f, double* vatom, + int* errorflag) +{ + arrdim2v(x, 3, *nmax); + arrdim2v(Arho1, 3, *nmax); + arrdim2v(Arho2, 6, *nmax); + arrdim2v(Arho3, 10, *nmax); + arrdim2v(Arho3b, 3, *nmax); + arrdim2v(t_ave, 3, *nmax); + arrdim2v(tsq_ave, 3, *nmax); + arrdim2v(f, 3, *nmax); + arrdim2v(vatom, 6, *nmax); - int i,j,jn,k,kn,kk,m,n,p,q; - int nv2,nv3,elti,eltj,eltk,ind; - double xitmp,yitmp,zitmp,delij[3+1],rij2,rij,rij3; - double delik[3+1],deljk[3+1],v[6+1],fi[3+1],fj[3+1]; - double third,sixth; - double pp,dUdrij,dUdsij,dUdrijm[3+1],force,forcem; - double r,recip,phi,phip; - double sij; - double a1,a1i,a1j,a2,a2i,a2j; - double a3i,a3j; - double shpi[3+1],shpj[3+1]; - double ai,aj,ro0i,ro0j,invrei,invrej; - double rhoa0j,drhoa0j,rhoa0i,drhoa0i; - double rhoa1j,drhoa1j,rhoa1i,drhoa1i; - double rhoa2j,drhoa2j,rhoa2i,drhoa2i; - double a3,a3a,rhoa3j,drhoa3j,rhoa3i,drhoa3i; - double drho0dr1,drho0dr2,drho0ds1,drho0ds2; - double drho1dr1,drho1dr2,drho1ds1,drho1ds2; - double drho1drm1[3+1],drho1drm2[3+1]; - double drho2dr1,drho2dr2,drho2ds1,drho2ds2; - double drho2drm1[3+1],drho2drm2[3+1]; - double drho3dr1,drho3dr2,drho3ds1,drho3ds2; - double drho3drm1[3+1],drho3drm2[3+1]; - double dt1dr1,dt1dr2,dt1ds1,dt1ds2; - double dt2dr1,dt2dr2,dt2ds1,dt2ds2; - double dt3dr1,dt3dr2,dt3ds1,dt3ds2; - double drhodr1,drhodr2,drhods1,drhods2,drhodrm1[3+1],drhodrm2[3+1]; - double arg; - double arg1i1,arg1j1,arg1i2,arg1j2,arg1i3,arg1j3,arg3i3,arg3j3; - double dsij1,dsij2,force1,force2; - double t1i,t2i,t3i,t1j,t2j,t3j; + int i, j, jn, k, kn, kk, m, n, p, q; + int nv2, nv3, elti, eltj, eltk, ind; + double xitmp, yitmp, zitmp, delij[3 + 1], rij2, rij, rij3; + double delik[3 + 1], deljk[3 + 1], v[6 + 1], fi[3 + 1], fj[3 + 1]; + double third, sixth; + double pp, dUdrij, dUdsij, dUdrijm[3 + 1], force, forcem; + double r, recip, phi, phip; + double sij; + double a1, a1i, a1j, a2, a2i, a2j; + double a3i, a3j; + double shpi[3 + 1], shpj[3 + 1]; + double ai, aj, ro0i, ro0j, invrei, invrej; + double rhoa0j, drhoa0j, rhoa0i, drhoa0i; + double rhoa1j, drhoa1j, rhoa1i, drhoa1i; + double rhoa2j, drhoa2j, rhoa2i, drhoa2i; + double a3, a3a, rhoa3j, drhoa3j, rhoa3i, drhoa3i; + double drho0dr1, drho0dr2, drho0ds1, drho0ds2; + double drho1dr1, drho1dr2, drho1ds1, drho1ds2; + double drho1drm1[3 + 1], drho1drm2[3 + 1]; + double drho2dr1, drho2dr2, drho2ds1, drho2ds2; + double drho2drm1[3 + 1], drho2drm2[3 + 1]; + double drho3dr1, drho3dr2, drho3ds1, drho3ds2; + double drho3drm1[3 + 1], drho3drm2[3 + 1]; + double dt1dr1, dt1dr2, dt1ds1, dt1ds2; + double dt2dr1, dt2dr2, dt2ds1, dt2ds2; + double dt3dr1, dt3dr2, dt3ds1, dt3ds2; + double drhodr1, drhodr2, drhods1, drhods2, drhodrm1[3 + 1], drhodrm2[3 + 1]; + double arg; + double arg1i1, arg1j1, arg1i2, arg1j2, arg1i3, arg1j3, arg3i3, arg3j3; + double dsij1, dsij2, force1, force2; + double t1i, t2i, t3i, t1j, t2j, t3j; - *errorflag = 0; - third = 1.0/3.0; - sixth = 1.0/6.0; - - //: aliased - i = *iptr; + *errorflag = 0; + third = 1.0 / 3.0; + sixth = 1.0 / 6.0; -// Compute forces atom i + //: aliased + i = *iptr; - elti = arr1v(fmap,arr1v(type,i)); + // Compute forces atom i - if (elti > 0) { - xitmp = arr2v(x,1,i); - yitmp = arr2v(x,2,i); - zitmp = arr2v(x,3,i); + elti = arr1v(fmap, arr1v(type, i)); -// Treat each pair - for(jn=1; jn<=*numneigh; jn++) { - j = arr1v(firstneigh,jn); - eltj = arr1v(fmap,arr1v(type,j)); + if (elti > 0) { + xitmp = arr2v(x, 1, i); + yitmp = arr2v(x, 2, i); + zitmp = arr2v(x, 3, i); - if (!iszero(arr1v(scrfcn,jn)) && eltj > 0) { + // Treat each pair + for (jn = 1; jn <= *numneigh; jn++) { + j = arr1v(firstneigh, jn); + eltj = arr1v(fmap, arr1v(type, j)); - sij = arr1v(scrfcn,jn)*arr1v(fcpair,jn); - delij[1] = arr2v(x,1,j) - xitmp; - delij[2] = arr2v(x,2,j) - yitmp; - delij[3] = arr2v(x,3,j) - zitmp; - rij2 = delij[1]*delij[1] + delij[2]*delij[2] + delij[3]*delij[3]; - if (rij2 < meam_data.cutforcesq) { - rij = sqrt(rij2); - r = rij; + if (!iszero(arr1v(scrfcn, jn)) && eltj > 0) { -// Compute phi and phip - ind = meam_data.eltind[elti][eltj]; - pp = rij*meam_data.rdrar + 1.0; - kk = (int)pp; - kk = min(kk,meam_data.nrar-1); - pp = pp - kk; - pp = min(pp,1.0); - phi = ((arr2(meam_data.phirar3,kk,ind)*pp + arr2(meam_data.phirar2,kk,ind))*pp + arr2(meam_data.phirar1,kk,ind))*pp + arr2(meam_data.phirar,kk,ind); - phip = (arr2(meam_data.phirar6,kk,ind)*pp + arr2(meam_data.phirar5,kk,ind))*pp + arr2(meam_data.phirar4,kk,ind); - recip = 1.0/r; + sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn); + delij[1] = arr2v(x, 1, j) - xitmp; + delij[2] = arr2v(x, 2, j) - yitmp; + delij[3] = arr2v(x, 3, j) - zitmp; + rij2 = delij[1] * delij[1] + delij[2] * delij[2] + delij[3] * delij[3]; + if (rij2 < meam_data.cutforcesq) { + rij = sqrt(rij2); + r = rij; - if (*eflag_either!=0) { - if (*eflag_global!=0) *eng_vdwl = *eng_vdwl + phi*sij; - if (*eflag_atom!=0) { - arr1v(eatom,i) = arr1v(eatom,i) + 0.5*phi*sij; - arr1v(eatom,j) = arr1v(eatom,j) + 0.5*phi*sij; - } - } + // Compute phi and phip + ind = meam_data.eltind[elti][eltj]; + pp = rij * meam_data.rdrar + 1.0; + kk = (int)pp; + kk = min(kk, meam_data.nrar - 1); + pp = pp - kk; + pp = min(pp, 1.0); + phi = ((arr2(meam_data.phirar3, kk, ind) * pp + + arr2(meam_data.phirar2, kk, ind)) * + pp + + arr2(meam_data.phirar1, kk, ind)) * + pp + + arr2(meam_data.phirar, kk, ind); + phip = (arr2(meam_data.phirar6, kk, ind) * pp + + arr2(meam_data.phirar5, kk, ind)) * + pp + + arr2(meam_data.phirar4, kk, ind); + recip = 1.0 / r; -// write(1,*) "force_meamf: phi: ",phi -// write(1,*) "force_meamf: phip: ",phip - -// Compute pair densities and derivatives - invrei = 1.0/meam_data.re_meam[elti][elti]; - ai = rij*invrei - 1.0; - ro0i = meam_data.rho0_meam[elti]; - rhoa0i = ro0i*fm_exp(-meam_data.beta0_meam[elti]*ai); - drhoa0i = -meam_data.beta0_meam[elti]*invrei*rhoa0i; - rhoa1i = ro0i*fm_exp(-meam_data.beta1_meam[elti]*ai); - drhoa1i = -meam_data.beta1_meam[elti]*invrei*rhoa1i; - rhoa2i = ro0i*fm_exp(-meam_data.beta2_meam[elti]*ai); - drhoa2i = -meam_data.beta2_meam[elti]*invrei*rhoa2i; - rhoa3i = ro0i*fm_exp(-meam_data.beta3_meam[elti]*ai); - drhoa3i = -meam_data.beta3_meam[elti]*invrei*rhoa3i; - - if (elti!=eltj) { - invrej = 1.0/meam_data.re_meam[eltj][eltj]; - aj = rij*invrej - 1.0; - ro0j = meam_data.rho0_meam[eltj]; - rhoa0j = ro0j*fm_exp(-meam_data.beta0_meam[eltj]*aj); - drhoa0j = -meam_data.beta0_meam[eltj]*invrej*rhoa0j; - rhoa1j = ro0j*fm_exp(-meam_data.beta1_meam[eltj]*aj); - drhoa1j = -meam_data.beta1_meam[eltj]*invrej*rhoa1j; - rhoa2j = ro0j*fm_exp(-meam_data.beta2_meam[eltj]*aj); - drhoa2j = -meam_data.beta2_meam[eltj]*invrej*rhoa2j; - rhoa3j = ro0j*fm_exp(-meam_data.beta3_meam[eltj]*aj); - drhoa3j = -meam_data.beta3_meam[eltj]*invrej*rhoa3j; - } else { - rhoa0j = rhoa0i; - drhoa0j = drhoa0i; - rhoa1j = rhoa1i; - drhoa1j = drhoa1i; - rhoa2j = rhoa2i; - drhoa2j = drhoa2i; - rhoa3j = rhoa3i; - drhoa3j = drhoa3i; - } - - if (meam_data.ialloy==1) { - rhoa1j = rhoa1j * meam_data.t1_meam[eltj]; - rhoa2j = rhoa2j * meam_data.t2_meam[eltj]; - rhoa3j = rhoa3j * meam_data.t3_meam[eltj]; - rhoa1i = rhoa1i * meam_data.t1_meam[elti]; - rhoa2i = rhoa2i * meam_data.t2_meam[elti]; - rhoa3i = rhoa3i * meam_data.t3_meam[elti]; - drhoa1j = drhoa1j * meam_data.t1_meam[eltj]; - drhoa2j = drhoa2j * meam_data.t2_meam[eltj]; - drhoa3j = drhoa3j * meam_data.t3_meam[eltj]; - drhoa1i = drhoa1i * meam_data.t1_meam[elti]; - drhoa2i = drhoa2i * meam_data.t2_meam[elti]; - drhoa3i = drhoa3i * meam_data.t3_meam[elti]; - } - - nv2 = 1; - nv3 = 1; - arg1i1 = 0.0; - arg1j1 = 0.0; - arg1i2 = 0.0; - arg1j2 = 0.0; - arg1i3 = 0.0; - arg1j3 = 0.0; - arg3i3 = 0.0; - arg3j3 = 0.0; - for(n=1; n<=3; n++) { - for(p=n; p<=3; p++) { - for(q=p; q<=3; q++) { - arg = delij[n]*delij[p]*delij[q]*meam_data.v3D[nv3]; - arg1i3 = arg1i3 + arr2v(Arho3,nv3,i)*arg; - arg1j3 = arg1j3 - arr2v(Arho3,nv3,j)*arg; - nv3 = nv3+1; - } - arg = delij[n]*delij[p]*meam_data.v2D[nv2]; - arg1i2 = arg1i2 + arr2v(Arho2,nv2,i)*arg; - arg1j2 = arg1j2 + arr2v(Arho2,nv2,j)*arg; - nv2 = nv2+1; - } - arg1i1 = arg1i1 + arr2v(Arho1,n,i)*delij[n]; - arg1j1 = arg1j1 - arr2v(Arho1,n,j)*delij[n]; - arg3i3 = arg3i3 + arr2v(Arho3b,n,i)*delij[n]; - arg3j3 = arg3j3 - arr2v(Arho3b,n,j)*delij[n]; - } - -// rho0 terms - drho0dr1 = drhoa0j * sij; - drho0dr2 = drhoa0i * sij; - -// rho1 terms - a1 = 2*sij/rij; - drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1; - drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1; - a1 = 2.0*sij/rij; - for(m=1; m<=3; m++) { - drho1drm1[m] = a1*rhoa1j*arr2v(Arho1,m,i); - drho1drm2[m] = -a1*rhoa1i*arr2v(Arho1,m,j); - } - -// rho2 terms - a2 = 2*sij/rij2; - drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2 - 2.0/3.0*arr1v(Arho2b,i)*drhoa2j*sij; - drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2 - 2.0/3.0*arr1v(Arho2b,j)*drhoa2i*sij; - a2 = 4*sij/rij2; - for(m=1; m<=3; m++) { - drho2drm1[m] = 0.0; - drho2drm2[m] = 0.0; - for(n=1; n<=3; n++) { - drho2drm1[m] = drho2drm1[m] + arr2v(Arho2,meam_data.vind2D[m][n],i)*delij[n]; - drho2drm2[m] = drho2drm2[m] - arr2v(Arho2,meam_data.vind2D[m][n],j)*delij[n]; - } - drho2drm1[m] = a2*rhoa2j*drho2drm1[m]; - drho2drm2[m] = -a2*rhoa2i*drho2drm2[m]; - } - -// rho3 terms - rij3 = rij*rij2; - a3 = 2*sij/rij3; - a3a = 6.0/5.0*sij/rij; - drho3dr1 = a3*(drhoa3j - 3*rhoa3j/rij)*arg1i3 - a3a*(drhoa3j - rhoa3j/rij)*arg3i3; - drho3dr2 = a3*(drhoa3i - 3*rhoa3i/rij)*arg1j3 - a3a*(drhoa3i - rhoa3i/rij)*arg3j3; - a3 = 6*sij/rij3; - a3a = 6*sij/(5*rij); - for(m=1; m<=3; m++) { - drho3drm1[m] = 0.0; - drho3drm2[m] = 0.0; - nv2 = 1; - for(n=1; n<=3; n++) { - for(p=n; p<=3; p++) { - arg = delij[n]*delij[p]*meam_data.v2D[nv2]; - drho3drm1[m] = drho3drm1[m] + arr2v(Arho3,meam_data.vind3D[m][n][p],i)*arg; - drho3drm2[m] = drho3drm2[m] + arr2v(Arho3,meam_data.vind3D[m][n][p],j)*arg; - nv2 = nv2 + 1; - } - } - drho3drm1[m] = ( a3*drho3drm1[m] - a3a*arr2v(Arho3b,m,i)) *rhoa3j; - drho3drm2[m] = (-a3*drho3drm2[m] + a3a*arr2v(Arho3b,m,j)) *rhoa3i; - } - -// Compute derivatives of weighting functions t wrt rij - t1i = arr2v(t_ave,1,i); - t2i = arr2v(t_ave,2,i); - t3i = arr2v(t_ave,3,i); - t1j = arr2v(t_ave,1,j); - t2j = arr2v(t_ave,2,j); - t3j = arr2v(t_ave,3,j); - - if (meam_data.ialloy==1) { - - a1i = 0.0; - a1j = 0.0; - a2i = 0.0; - a2j = 0.0; - a3i = 0.0; - a3j = 0.0; - if ( !iszero(arr2v(tsq_ave,1,i)) ) a1i = drhoa0j*sij/arr2v(tsq_ave,1,i); - if ( !iszero(arr2v(tsq_ave,1,j)) ) a1j = drhoa0i*sij/arr2v(tsq_ave,1,j); - if ( !iszero(arr2v(tsq_ave,2,i)) ) a2i = drhoa0j*sij/arr2v(tsq_ave,2,i); - if ( !iszero(arr2v(tsq_ave,2,j)) ) a2j = drhoa0i*sij/arr2v(tsq_ave,2,j); - if ( !iszero(arr2v(tsq_ave,3,i)) ) a3i = drhoa0j*sij/arr2v(tsq_ave,3,i); - if ( !iszero(arr2v(tsq_ave,3,j)) ) a3j = drhoa0i*sij/arr2v(tsq_ave,3,j); - - dt1dr1 = a1i*(meam_data.t1_meam[eltj]-t1i*pow(meam_data.t1_meam[eltj],2)); - dt1dr2 = a1j*(meam_data.t1_meam[elti]-t1j*pow(meam_data.t1_meam[elti],2)); - dt2dr1 = a2i*(meam_data.t2_meam[eltj]-t2i*pow(meam_data.t2_meam[eltj],2)); - dt2dr2 = a2j*(meam_data.t2_meam[elti]-t2j*pow(meam_data.t2_meam[elti],2)); - dt3dr1 = a3i*(meam_data.t3_meam[eltj]-t3i*pow(meam_data.t3_meam[eltj],2)); - dt3dr2 = a3j*(meam_data.t3_meam[elti]-t3j*pow(meam_data.t3_meam[elti],2)); - - } else if (meam_data.ialloy==2) { - - dt1dr1 = 0.0; - dt1dr2 = 0.0; - dt2dr1 = 0.0; - dt2dr2 = 0.0; - dt3dr1 = 0.0; - dt3dr2 = 0.0; - - } else { - - ai = 0.0; - if( !iszero(arr1v(rho0,i)) ) - ai = drhoa0j*sij/arr1v(rho0,i); - aj = 0.0; - if( !iszero(arr1v(rho0,j)) ) - aj = drhoa0i*sij/arr1v(rho0,j); - - dt1dr1 = ai*(meam_data.t1_meam[elti]-t1i); - dt1dr2 = aj*(meam_data.t1_meam[elti]-t1j); - dt2dr1 = ai*(meam_data.t2_meam[elti]-t2i); - dt2dr2 = aj*(meam_data.t2_meam[elti]-t2j); - dt3dr1 = ai*(meam_data.t3_meam[elti]-t3i); - dt3dr2 = aj*(meam_data.t3_meam[elti]-t3j); - - } - -// Compute derivatives of total density wrt rij, sij and rij(3) - get_shpfcn(shpi,meam_data.lattce_meam[elti][elti]); - get_shpfcn(shpj,meam_data.lattce_meam[eltj][eltj]); - drhodr1 = arr1v(dGamma1,i)*drho0dr1 - + arr1v(dGamma2,i)* - (dt1dr1*arr1v(rho1,i)+t1i*drho1dr1 - + dt2dr1*arr1v(rho2,i)+t2i*drho2dr1 - + dt3dr1*arr1v(rho3,i)+t3i*drho3dr1) - - arr1v(dGamma3,i)* - (shpi[1]*dt1dr1+shpi[2]*dt2dr1+shpi[3]*dt3dr1); - drhodr2 = arr1v(dGamma1,j)*drho0dr2 - + arr1v(dGamma2,j)* - (dt1dr2*arr1v(rho1,j)+t1j*drho1dr2 - + dt2dr2*arr1v(rho2,j)+t2j*drho2dr2 - + dt3dr2*arr1v(rho3,j)+t3j*drho3dr2) - - arr1v(dGamma3,j)* - (shpj[1]*dt1dr2+shpj[2]*dt2dr2+shpj[3]*dt3dr2); - for(m=1; m<=3; m++) { - drhodrm1[m] = 0.0; - drhodrm2[m] = 0.0; - drhodrm1[m] = arr1v(dGamma2,i)* - (t1i*drho1drm1[m] - + t2i*drho2drm1[m] - + t3i*drho3drm1[m]); - drhodrm2[m] = arr1v(dGamma2,j)* - (t1j*drho1drm2[m] - + t2j*drho2drm2[m] - + t3j*drho3drm2[m]); - } - -// Compute derivatives wrt sij, but only if necessary - if ( !iszero(arr1v(dscrfcn,jn))) { - drho0ds1 = rhoa0j; - drho0ds2 = rhoa0i; - a1 = 2.0/rij; - drho1ds1 = a1*rhoa1j*arg1i1; - drho1ds2 = a1*rhoa1i*arg1j1; - a2 = 2.0/rij2; - drho2ds1 = a2*rhoa2j*arg1i2 - 2.0/3.0*arr1v(Arho2b,i)*rhoa2j; - drho2ds2 = a2*rhoa2i*arg1j2 - 2.0/3.0*arr1v(Arho2b,j)*rhoa2i; - a3 = 2.0/rij3; - a3a = 6.0/(5.0*rij); - drho3ds1 = a3*rhoa3j*arg1i3 - a3a*rhoa3j*arg3i3; - drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3; - - if (meam_data.ialloy==1) { - - a1i = 0.0; - a1j = 0.0; - a2i = 0.0; - a2j = 0.0; - a3i = 0.0; - a3j = 0.0; - if ( !iszero(arr2v(tsq_ave,1,i)) ) a1i = rhoa0j/arr2v(tsq_ave,1,i); - if ( !iszero(arr2v(tsq_ave,1,j)) ) a1j = rhoa0i/arr2v(tsq_ave,1,j); - if ( !iszero(arr2v(tsq_ave,2,i)) ) a2i = rhoa0j/arr2v(tsq_ave,2,i); - if ( !iszero(arr2v(tsq_ave,2,j)) ) a2j = rhoa0i/arr2v(tsq_ave,2,j); - if ( !iszero(arr2v(tsq_ave,3,i)) ) a3i = rhoa0j/arr2v(tsq_ave,3,i); - if ( !iszero(arr2v(tsq_ave,3,j)) ) a3j = rhoa0i/arr2v(tsq_ave,3,j); - - dt1ds1 = a1i*(meam_data.t1_meam[eltj]-t1i*pow(meam_data.t1_meam[eltj],2)); - dt1ds2 = a1j*(meam_data.t1_meam[elti]-t1j*pow(meam_data.t1_meam[elti],2)); - dt2ds1 = a2i*(meam_data.t2_meam[eltj]-t2i*pow(meam_data.t2_meam[eltj],2)); - dt2ds2 = a2j*(meam_data.t2_meam[elti]-t2j*pow(meam_data.t2_meam[elti],2)); - dt3ds1 = a3i*(meam_data.t3_meam[eltj]-t3i*pow(meam_data.t3_meam[eltj],2)); - dt3ds2 = a3j*(meam_data.t3_meam[elti]-t3j*pow(meam_data.t3_meam[elti],2)); - - } else if (meam_data.ialloy==2) { - - dt1ds1 = 0.0; - dt1ds2 = 0.0; - dt2ds1 = 0.0; - dt2ds2 = 0.0; - dt3ds1 = 0.0; - dt3ds2 = 0.0; - - } else { - - ai = 0.0; - if( !iszero(arr1v(rho0,i)) ) - ai = rhoa0j/arr1v(rho0,i); - aj = 0.0; - if( !iszero(arr1v(rho0,j)) ) - aj = rhoa0i/arr1v(rho0,j); - - dt1ds1 = ai*(meam_data.t1_meam[eltj]-t1i); - dt1ds2 = aj*(meam_data.t1_meam[elti]-t1j); - dt2ds1 = ai*(meam_data.t2_meam[eltj]-t2i); - dt2ds2 = aj*(meam_data.t2_meam[elti]-t2j); - dt3ds1 = ai*(meam_data.t3_meam[eltj]-t3i); - dt3ds2 = aj*(meam_data.t3_meam[elti]-t3j); - - } - - drhods1 = arr1v(dGamma1,i)*drho0ds1 - + arr1v(dGamma2,i)* - (dt1ds1*arr1v(rho1,i)+t1i*drho1ds1 - + dt2ds1*arr1v(rho2,i)+t2i*drho2ds1 - + dt3ds1*arr1v(rho3,i)+t3i*drho3ds1) - - arr1v(dGamma3,i)* - (shpi[1]*dt1ds1+shpi[2]*dt2ds1+shpi[3]*dt3ds1); - drhods2 = arr1v(dGamma1,j)*drho0ds2 - + arr1v(dGamma2,j)* - (dt1ds2*arr1v(rho1,j)+t1j*drho1ds2 - + dt2ds2*arr1v(rho2,j)+t2j*drho2ds2 - + dt3ds2*arr1v(rho3,j)+t3j*drho3ds2) - - arr1v(dGamma3,j)* - (shpj[1]*dt1ds2+shpj[2]*dt2ds2+shpj[3]*dt3ds2); - } - -// Compute derivatives of energy wrt rij, sij and rij[3] - dUdrij = phip*sij + arr1v(fp,i)*drhodr1 + arr1v(fp,j)*drhodr2; - dUdsij = 0.0; - if ( !iszero(arr1v(dscrfcn,jn)) ) { - dUdsij = phi + arr1v(fp,i)*drhods1 + arr1v(fp,j)*drhods2; - } - for(m=1; m<=3; m++) { - dUdrijm[m] = arr1v(fp,i)*drhodrm1[m] + arr1v(fp,j)*drhodrm2[m]; - } - -// Add the part of the force due to dUdrij and dUdsij - - force = dUdrij*recip + dUdsij*arr1v(dscrfcn,jn); - for(m=1; m<=3; m++) { - forcem = delij[m]*force + dUdrijm[m]; - arr2v(f,m,i) = arr2v(f,m,i) + forcem; - arr2v(f,m,j) = arr2v(f,m,j) - forcem; - } - -// Tabulate per-atom virial as symmetrized stress tensor - - if (*vflag_atom!=0) { - fi[1] = delij[1]*force + dUdrijm[1]; - fi[2] = delij[2]*force + dUdrijm[2]; - fi[3] = delij[3]*force + dUdrijm[3]; - v[1] = -0.5 * (delij[1] * fi[1]); - v[2] = -0.5 * (delij[2] * fi[2]); - v[3] = -0.5 * (delij[3] * fi[3]); - v[4] = -0.25 * (delij[1]*fi[2] + delij[2]*fi[1]); - v[5] = -0.25 * (delij[1]*fi[3] + delij[3]*fi[1]); - v[6] = -0.25 * (delij[2]*fi[3] + delij[3]*fi[2]); - - arr2v(vatom,1,i) = arr2v(vatom,1,i) + v[1]; - arr2v(vatom,2,i) = arr2v(vatom,2,i) + v[2]; - arr2v(vatom,3,i) = arr2v(vatom,3,i) + v[3]; - arr2v(vatom,4,i) = arr2v(vatom,4,i) + v[4]; - arr2v(vatom,5,i) = arr2v(vatom,5,i) + v[5]; - arr2v(vatom,6,i) = arr2v(vatom,6,i) + v[6]; - arr2v(vatom,1,j) = arr2v(vatom,1,j) + v[1]; - arr2v(vatom,2,j) = arr2v(vatom,2,j) + v[2]; - arr2v(vatom,3,j) = arr2v(vatom,3,j) + v[3]; - arr2v(vatom,4,j) = arr2v(vatom,4,j) + v[4]; - arr2v(vatom,5,j) = arr2v(vatom,5,j) + v[5]; - arr2v(vatom,6,j) = arr2v(vatom,6,j) + v[6]; - } - -// Now compute forces on other atoms k due to change in sij - - if (iszero(sij) || iszero(sij-1.0)) continue; //: cont jn loop - for(kn=1; kn<=*numneigh_full; kn++) { - k = arr1v(firstneigh_full,kn); - eltk = arr1v(fmap,arr1v(type,k)); - if (k!=j && eltk > 0) { - dsij(i,j,k,jn,*nmax,*numneigh,rij2,&dsij1,&dsij2,*ntype,type,fmap,x,scrfcn,fcpair); - if (!iszero(dsij1) || !iszero(dsij2)) { - force1 = dUdsij*dsij1; - force2 = dUdsij*dsij2; - for(m=1; m<=3; m++) { - delik[m] = arr2v(x,m,k) - arr2v(x,m,i); - deljk[m] = arr2v(x,m,k) - arr2v(x,m,j); - } - for(m=1; m<=3; m++) { - arr2v(f,m,i) = arr2v(f,m,i) + force1*delik[m]; - arr2v(f,m,j) = arr2v(f,m,j) + force2*deljk[m]; - arr2v(f,m,k) = arr2v(f,m,k) - force1*delik[m] - force2*deljk[m]; - } - -// Tabulate per-atom virial as symmetrized stress tensor - - if (*vflag_atom!=0) { - fi[1] = force1*delik[1]; - fi[2] = force1*delik[2]; - fi[3] = force1*delik[3]; - fj[1] = force2*deljk[1]; - fj[2] = force2*deljk[2]; - fj[3] = force2*deljk[3]; - v[1] = -third * (delik[1]*fi[1] + deljk[1]*fj[1]); - v[2] = -third * (delik[2]*fi[2] + deljk[2]*fj[2]); - v[3] = -third * (delik[3]*fi[3] + deljk[3]*fj[3]); - v[4] = -sixth * (delik[1]*fi[2] + deljk[1]*fj[2] + delik[2]*fi[1] + deljk[2]*fj[1]); - v[5] = -sixth * (delik[1]*fi[3] + deljk[1]*fj[3] + delik[3]*fi[1] + deljk[3]*fj[1]); - v[6] = -sixth * (delik[2]*fi[3] + deljk[2]*fj[3] + delik[3]*fi[2] + deljk[3]*fj[2]); - - arr2v(vatom,1,i) = arr2v(vatom,1,i) + v[1]; - arr2v(vatom,2,i) = arr2v(vatom,2,i) + v[2]; - arr2v(vatom,3,i) = arr2v(vatom,3,i) + v[3]; - arr2v(vatom,4,i) = arr2v(vatom,4,i) + v[4]; - arr2v(vatom,5,i) = arr2v(vatom,5,i) + v[5]; - arr2v(vatom,6,i) = arr2v(vatom,6,i) + v[6]; - arr2v(vatom,1,j) = arr2v(vatom,1,j) + v[1]; - arr2v(vatom,2,j) = arr2v(vatom,2,j) + v[2]; - arr2v(vatom,3,j) = arr2v(vatom,3,j) + v[3]; - arr2v(vatom,4,j) = arr2v(vatom,4,j) + v[4]; - arr2v(vatom,5,j) = arr2v(vatom,5,j) + v[5]; - arr2v(vatom,6,j) = arr2v(vatom,6,j) + v[6]; - arr2v(vatom,1,k) = arr2v(vatom,1,k) + v[1]; - arr2v(vatom,2,k) = arr2v(vatom,2,k) + v[2]; - arr2v(vatom,3,k) = arr2v(vatom,3,k) + v[3]; - arr2v(vatom,4,k) = arr2v(vatom,4,k) + v[4]; - arr2v(vatom,5,k) = arr2v(vatom,5,k) + v[5]; - arr2v(vatom,6,k) = arr2v(vatom,6,k) + v[6]; - } - - } - } -// end of k loop - } + if (*eflag_either != 0) { + if (*eflag_global != 0) + *eng_vdwl = *eng_vdwl + phi * sij; + if (*eflag_atom != 0) { + arr1v(eatom, i) = arr1v(eatom, i) + 0.5 * phi * sij; + arr1v(eatom, j) = arr1v(eatom, j) + 0.5 * phi * sij; } } -// end of j loop + + // write(1,*) "force_meamf: phi: ",phi + // write(1,*) "force_meamf: phip: ",phip + + // Compute pair densities and derivatives + invrei = 1.0 / meam_data.re_meam[elti][elti]; + ai = rij * invrei - 1.0; + ro0i = meam_data.rho0_meam[elti]; + rhoa0i = ro0i * fm_exp(-meam_data.beta0_meam[elti] * ai); + drhoa0i = -meam_data.beta0_meam[elti] * invrei * rhoa0i; + rhoa1i = ro0i * fm_exp(-meam_data.beta1_meam[elti] * ai); + drhoa1i = -meam_data.beta1_meam[elti] * invrei * rhoa1i; + rhoa2i = ro0i * fm_exp(-meam_data.beta2_meam[elti] * ai); + drhoa2i = -meam_data.beta2_meam[elti] * invrei * rhoa2i; + rhoa3i = ro0i * fm_exp(-meam_data.beta3_meam[elti] * ai); + drhoa3i = -meam_data.beta3_meam[elti] * invrei * rhoa3i; + + if (elti != eltj) { + invrej = 1.0 / meam_data.re_meam[eltj][eltj]; + aj = rij * invrej - 1.0; + ro0j = meam_data.rho0_meam[eltj]; + rhoa0j = ro0j * fm_exp(-meam_data.beta0_meam[eltj] * aj); + drhoa0j = -meam_data.beta0_meam[eltj] * invrej * rhoa0j; + rhoa1j = ro0j * fm_exp(-meam_data.beta1_meam[eltj] * aj); + drhoa1j = -meam_data.beta1_meam[eltj] * invrej * rhoa1j; + rhoa2j = ro0j * fm_exp(-meam_data.beta2_meam[eltj] * aj); + drhoa2j = -meam_data.beta2_meam[eltj] * invrej * rhoa2j; + rhoa3j = ro0j * fm_exp(-meam_data.beta3_meam[eltj] * aj); + drhoa3j = -meam_data.beta3_meam[eltj] * invrej * rhoa3j; + } else { + rhoa0j = rhoa0i; + drhoa0j = drhoa0i; + rhoa1j = rhoa1i; + drhoa1j = drhoa1i; + rhoa2j = rhoa2i; + drhoa2j = drhoa2i; + rhoa3j = rhoa3i; + drhoa3j = drhoa3i; + } + + if (meam_data.ialloy == 1) { + rhoa1j = rhoa1j * meam_data.t1_meam[eltj]; + rhoa2j = rhoa2j * meam_data.t2_meam[eltj]; + rhoa3j = rhoa3j * meam_data.t3_meam[eltj]; + rhoa1i = rhoa1i * meam_data.t1_meam[elti]; + rhoa2i = rhoa2i * meam_data.t2_meam[elti]; + rhoa3i = rhoa3i * meam_data.t3_meam[elti]; + drhoa1j = drhoa1j * meam_data.t1_meam[eltj]; + drhoa2j = drhoa2j * meam_data.t2_meam[eltj]; + drhoa3j = drhoa3j * meam_data.t3_meam[eltj]; + drhoa1i = drhoa1i * meam_data.t1_meam[elti]; + drhoa2i = drhoa2i * meam_data.t2_meam[elti]; + drhoa3i = drhoa3i * meam_data.t3_meam[elti]; + } + + nv2 = 1; + nv3 = 1; + arg1i1 = 0.0; + arg1j1 = 0.0; + arg1i2 = 0.0; + arg1j2 = 0.0; + arg1i3 = 0.0; + arg1j3 = 0.0; + arg3i3 = 0.0; + arg3j3 = 0.0; + for (n = 1; n <= 3; n++) { + for (p = n; p <= 3; p++) { + for (q = p; q <= 3; q++) { + arg = delij[n] * delij[p] * delij[q] * meam_data.v3D[nv3]; + arg1i3 = arg1i3 + arr2v(Arho3, nv3, i) * arg; + arg1j3 = arg1j3 - arr2v(Arho3, nv3, j) * arg; + nv3 = nv3 + 1; + } + arg = delij[n] * delij[p] * meam_data.v2D[nv2]; + arg1i2 = arg1i2 + arr2v(Arho2, nv2, i) * arg; + arg1j2 = arg1j2 + arr2v(Arho2, nv2, j) * arg; + nv2 = nv2 + 1; + } + arg1i1 = arg1i1 + arr2v(Arho1, n, i) * delij[n]; + arg1j1 = arg1j1 - arr2v(Arho1, n, j) * delij[n]; + arg3i3 = arg3i3 + arr2v(Arho3b, n, i) * delij[n]; + arg3j3 = arg3j3 - arr2v(Arho3b, n, j) * delij[n]; + } + + // rho0 terms + drho0dr1 = drhoa0j * sij; + drho0dr2 = drhoa0i * sij; + + // rho1 terms + a1 = 2 * sij / rij; + drho1dr1 = a1 * (drhoa1j - rhoa1j / rij) * arg1i1; + drho1dr2 = a1 * (drhoa1i - rhoa1i / rij) * arg1j1; + a1 = 2.0 * sij / rij; + for (m = 1; m <= 3; m++) { + drho1drm1[m] = a1 * rhoa1j * arr2v(Arho1, m, i); + drho1drm2[m] = -a1 * rhoa1i * arr2v(Arho1, m, j); + } + + // rho2 terms + a2 = 2 * sij / rij2; + drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 - + 2.0 / 3.0 * arr1v(Arho2b, i) * drhoa2j * sij; + drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 - + 2.0 / 3.0 * arr1v(Arho2b, j) * drhoa2i * sij; + a2 = 4 * sij / rij2; + for (m = 1; m <= 3; m++) { + drho2drm1[m] = 0.0; + drho2drm2[m] = 0.0; + for (n = 1; n <= 3; n++) { + drho2drm1[m] = drho2drm1[m] + + arr2v(Arho2, meam_data.vind2D[m][n], i) * delij[n]; + drho2drm2[m] = drho2drm2[m] - + arr2v(Arho2, meam_data.vind2D[m][n], j) * delij[n]; + } + drho2drm1[m] = a2 * rhoa2j * drho2drm1[m]; + drho2drm2[m] = -a2 * rhoa2i * drho2drm2[m]; + } + + // rho3 terms + rij3 = rij * rij2; + a3 = 2 * sij / rij3; + a3a = 6.0 / 5.0 * sij / rij; + drho3dr1 = a3 * (drhoa3j - 3 * rhoa3j / rij) * arg1i3 - + a3a * (drhoa3j - rhoa3j / rij) * arg3i3; + drho3dr2 = a3 * (drhoa3i - 3 * rhoa3i / rij) * arg1j3 - + a3a * (drhoa3i - rhoa3i / rij) * arg3j3; + a3 = 6 * sij / rij3; + a3a = 6 * sij / (5 * rij); + for (m = 1; m <= 3; m++) { + drho3drm1[m] = 0.0; + drho3drm2[m] = 0.0; + nv2 = 1; + for (n = 1; n <= 3; n++) { + for (p = n; p <= 3; p++) { + arg = delij[n] * delij[p] * meam_data.v2D[nv2]; + drho3drm1[m] = drho3drm1[m] + + arr2v(Arho3, meam_data.vind3D[m][n][p], i) * arg; + drho3drm2[m] = drho3drm2[m] + + arr2v(Arho3, meam_data.vind3D[m][n][p], j) * arg; + nv2 = nv2 + 1; + } + } + drho3drm1[m] = + (a3 * drho3drm1[m] - a3a * arr2v(Arho3b, m, i)) * rhoa3j; + drho3drm2[m] = + (-a3 * drho3drm2[m] + a3a * arr2v(Arho3b, m, j)) * rhoa3i; + } + + // Compute derivatives of weighting functions t wrt rij + t1i = arr2v(t_ave, 1, i); + t2i = arr2v(t_ave, 2, i); + t3i = arr2v(t_ave, 3, i); + t1j = arr2v(t_ave, 1, j); + t2j = arr2v(t_ave, 2, j); + t3j = arr2v(t_ave, 3, j); + + if (meam_data.ialloy == 1) { + + a1i = 0.0; + a1j = 0.0; + a2i = 0.0; + a2j = 0.0; + a3i = 0.0; + a3j = 0.0; + if (!iszero(arr2v(tsq_ave, 1, i))) + a1i = drhoa0j * sij / arr2v(tsq_ave, 1, i); + if (!iszero(arr2v(tsq_ave, 1, j))) + a1j = drhoa0i * sij / arr2v(tsq_ave, 1, j); + if (!iszero(arr2v(tsq_ave, 2, i))) + a2i = drhoa0j * sij / arr2v(tsq_ave, 2, i); + if (!iszero(arr2v(tsq_ave, 2, j))) + a2j = drhoa0i * sij / arr2v(tsq_ave, 2, j); + if (!iszero(arr2v(tsq_ave, 3, i))) + a3i = drhoa0j * sij / arr2v(tsq_ave, 3, i); + if (!iszero(arr2v(tsq_ave, 3, j))) + a3j = drhoa0i * sij / arr2v(tsq_ave, 3, j); + + dt1dr1 = a1i * (meam_data.t1_meam[eltj] - + t1i * pow(meam_data.t1_meam[eltj], 2)); + dt1dr2 = a1j * (meam_data.t1_meam[elti] - + t1j * pow(meam_data.t1_meam[elti], 2)); + dt2dr1 = a2i * (meam_data.t2_meam[eltj] - + t2i * pow(meam_data.t2_meam[eltj], 2)); + dt2dr2 = a2j * (meam_data.t2_meam[elti] - + t2j * pow(meam_data.t2_meam[elti], 2)); + dt3dr1 = a3i * (meam_data.t3_meam[eltj] - + t3i * pow(meam_data.t3_meam[eltj], 2)); + dt3dr2 = a3j * (meam_data.t3_meam[elti] - + t3j * pow(meam_data.t3_meam[elti], 2)); + + } else if (meam_data.ialloy == 2) { + + dt1dr1 = 0.0; + dt1dr2 = 0.0; + dt2dr1 = 0.0; + dt2dr2 = 0.0; + dt3dr1 = 0.0; + dt3dr2 = 0.0; + + } else { + + ai = 0.0; + if (!iszero(arr1v(rho0, i))) + ai = drhoa0j * sij / arr1v(rho0, i); + aj = 0.0; + if (!iszero(arr1v(rho0, j))) + aj = drhoa0i * sij / arr1v(rho0, j); + + dt1dr1 = ai * (meam_data.t1_meam[elti] - t1i); + dt1dr2 = aj * (meam_data.t1_meam[elti] - t1j); + dt2dr1 = ai * (meam_data.t2_meam[elti] - t2i); + dt2dr2 = aj * (meam_data.t2_meam[elti] - t2j); + dt3dr1 = ai * (meam_data.t3_meam[elti] - t3i); + dt3dr2 = aj * (meam_data.t3_meam[elti] - t3j); + } + + // Compute derivatives of total density wrt rij, sij and rij(3) + get_shpfcn(shpi, meam_data.lattce_meam[elti][elti]); + get_shpfcn(shpj, meam_data.lattce_meam[eltj][eltj]); + drhodr1 = + arr1v(dGamma1, i) * drho0dr1 + + arr1v(dGamma2, i) * (dt1dr1 * arr1v(rho1, i) + t1i * drho1dr1 + + dt2dr1 * arr1v(rho2, i) + t2i * drho2dr1 + + dt3dr1 * arr1v(rho3, i) + t3i * drho3dr1) - + arr1v(dGamma3, i) * + (shpi[1] * dt1dr1 + shpi[2] * dt2dr1 + shpi[3] * dt3dr1); + drhodr2 = + arr1v(dGamma1, j) * drho0dr2 + + arr1v(dGamma2, j) * (dt1dr2 * arr1v(rho1, j) + t1j * drho1dr2 + + dt2dr2 * arr1v(rho2, j) + t2j * drho2dr2 + + dt3dr2 * arr1v(rho3, j) + t3j * drho3dr2) - + arr1v(dGamma3, j) * + (shpj[1] * dt1dr2 + shpj[2] * dt2dr2 + shpj[3] * dt3dr2); + for (m = 1; m <= 3; m++) { + drhodrm1[m] = 0.0; + drhodrm2[m] = 0.0; + drhodrm1[m] = + arr1v(dGamma2, i) * + (t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]); + drhodrm2[m] = + arr1v(dGamma2, j) * + (t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]); + } + + // Compute derivatives wrt sij, but only if necessary + if (!iszero(arr1v(dscrfcn, jn))) { + drho0ds1 = rhoa0j; + drho0ds2 = rhoa0i; + a1 = 2.0 / rij; + drho1ds1 = a1 * rhoa1j * arg1i1; + drho1ds2 = a1 * rhoa1i * arg1j1; + a2 = 2.0 / rij2; + drho2ds1 = + a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arr1v(Arho2b, i) * rhoa2j; + drho2ds2 = + a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arr1v(Arho2b, j) * rhoa2i; + a3 = 2.0 / rij3; + a3a = 6.0 / (5.0 * rij); + drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3; + drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3; + + if (meam_data.ialloy == 1) { + + a1i = 0.0; + a1j = 0.0; + a2i = 0.0; + a2j = 0.0; + a3i = 0.0; + a3j = 0.0; + if (!iszero(arr2v(tsq_ave, 1, i))) + a1i = rhoa0j / arr2v(tsq_ave, 1, i); + if (!iszero(arr2v(tsq_ave, 1, j))) + a1j = rhoa0i / arr2v(tsq_ave, 1, j); + if (!iszero(arr2v(tsq_ave, 2, i))) + a2i = rhoa0j / arr2v(tsq_ave, 2, i); + if (!iszero(arr2v(tsq_ave, 2, j))) + a2j = rhoa0i / arr2v(tsq_ave, 2, j); + if (!iszero(arr2v(tsq_ave, 3, i))) + a3i = rhoa0j / arr2v(tsq_ave, 3, i); + if (!iszero(arr2v(tsq_ave, 3, j))) + a3j = rhoa0i / arr2v(tsq_ave, 3, j); + + dt1ds1 = a1i * (meam_data.t1_meam[eltj] - + t1i * pow(meam_data.t1_meam[eltj], 2)); + dt1ds2 = a1j * (meam_data.t1_meam[elti] - + t1j * pow(meam_data.t1_meam[elti], 2)); + dt2ds1 = a2i * (meam_data.t2_meam[eltj] - + t2i * pow(meam_data.t2_meam[eltj], 2)); + dt2ds2 = a2j * (meam_data.t2_meam[elti] - + t2j * pow(meam_data.t2_meam[elti], 2)); + dt3ds1 = a3i * (meam_data.t3_meam[eltj] - + t3i * pow(meam_data.t3_meam[eltj], 2)); + dt3ds2 = a3j * (meam_data.t3_meam[elti] - + t3j * pow(meam_data.t3_meam[elti], 2)); + + } else if (meam_data.ialloy == 2) { + + dt1ds1 = 0.0; + dt1ds2 = 0.0; + dt2ds1 = 0.0; + dt2ds2 = 0.0; + dt3ds1 = 0.0; + dt3ds2 = 0.0; + + } else { + + ai = 0.0; + if (!iszero(arr1v(rho0, i))) + ai = rhoa0j / arr1v(rho0, i); + aj = 0.0; + if (!iszero(arr1v(rho0, j))) + aj = rhoa0i / arr1v(rho0, j); + + dt1ds1 = ai * (meam_data.t1_meam[eltj] - t1i); + dt1ds2 = aj * (meam_data.t1_meam[elti] - t1j); + dt2ds1 = ai * (meam_data.t2_meam[eltj] - t2i); + dt2ds2 = aj * (meam_data.t2_meam[elti] - t2j); + dt3ds1 = ai * (meam_data.t3_meam[eltj] - t3i); + dt3ds2 = aj * (meam_data.t3_meam[elti] - t3j); + } + + drhods1 = + arr1v(dGamma1, i) * drho0ds1 + + arr1v(dGamma2, i) * (dt1ds1 * arr1v(rho1, i) + t1i * drho1ds1 + + dt2ds1 * arr1v(rho2, i) + t2i * drho2ds1 + + dt3ds1 * arr1v(rho3, i) + t3i * drho3ds1) - + arr1v(dGamma3, i) * + (shpi[1] * dt1ds1 + shpi[2] * dt2ds1 + shpi[3] * dt3ds1); + drhods2 = + arr1v(dGamma1, j) * drho0ds2 + + arr1v(dGamma2, j) * (dt1ds2 * arr1v(rho1, j) + t1j * drho1ds2 + + dt2ds2 * arr1v(rho2, j) + t2j * drho2ds2 + + dt3ds2 * arr1v(rho3, j) + t3j * drho3ds2) - + arr1v(dGamma3, j) * + (shpj[1] * dt1ds2 + shpj[2] * dt2ds2 + shpj[3] * dt3ds2); + } + + // Compute derivatives of energy wrt rij, sij and rij[3] + dUdrij = phip * sij + arr1v(fp, i) * drhodr1 + arr1v(fp, j) * drhodr2; + dUdsij = 0.0; + if (!iszero(arr1v(dscrfcn, jn))) { + dUdsij = phi + arr1v(fp, i) * drhods1 + arr1v(fp, j) * drhods2; + } + for (m = 1; m <= 3; m++) { + dUdrijm[m] = + arr1v(fp, i) * drhodrm1[m] + arr1v(fp, j) * drhodrm2[m]; + } + + // Add the part of the force due to dUdrij and dUdsij + + force = dUdrij * recip + dUdsij * arr1v(dscrfcn, jn); + for (m = 1; m <= 3; m++) { + forcem = delij[m] * force + dUdrijm[m]; + arr2v(f, m, i) = arr2v(f, m, i) + forcem; + arr2v(f, m, j) = arr2v(f, m, j) - forcem; + } + + // Tabulate per-atom virial as symmetrized stress tensor + + if (*vflag_atom != 0) { + fi[1] = delij[1] * force + dUdrijm[1]; + fi[2] = delij[2] * force + dUdrijm[2]; + fi[3] = delij[3] * force + dUdrijm[3]; + v[1] = -0.5 * (delij[1] * fi[1]); + v[2] = -0.5 * (delij[2] * fi[2]); + v[3] = -0.5 * (delij[3] * fi[3]); + v[4] = -0.25 * (delij[1] * fi[2] + delij[2] * fi[1]); + v[5] = -0.25 * (delij[1] * fi[3] + delij[3] * fi[1]); + v[6] = -0.25 * (delij[2] * fi[3] + delij[3] * fi[2]); + + arr2v(vatom, 1, i) = arr2v(vatom, 1, i) + v[1]; + arr2v(vatom, 2, i) = arr2v(vatom, 2, i) + v[2]; + arr2v(vatom, 3, i) = arr2v(vatom, 3, i) + v[3]; + arr2v(vatom, 4, i) = arr2v(vatom, 4, i) + v[4]; + arr2v(vatom, 5, i) = arr2v(vatom, 5, i) + v[5]; + arr2v(vatom, 6, i) = arr2v(vatom, 6, i) + v[6]; + arr2v(vatom, 1, j) = arr2v(vatom, 1, j) + v[1]; + arr2v(vatom, 2, j) = arr2v(vatom, 2, j) + v[2]; + arr2v(vatom, 3, j) = arr2v(vatom, 3, j) + v[3]; + arr2v(vatom, 4, j) = arr2v(vatom, 4, j) + v[4]; + arr2v(vatom, 5, j) = arr2v(vatom, 5, j) + v[5]; + arr2v(vatom, 6, j) = arr2v(vatom, 6, j) + v[6]; + } + + // Now compute forces on other atoms k due to change in sij + + if (iszero(sij) || iszero(sij - 1.0)) + continue; //: cont jn loop + for (kn = 1; kn <= *numneigh_full; kn++) { + k = arr1v(firstneigh_full, kn); + eltk = arr1v(fmap, arr1v(type, k)); + if (k != j && eltk > 0) { + dsij(i, j, k, jn, *nmax, *numneigh, rij2, &dsij1, &dsij2, *ntype, + type, fmap, x, scrfcn, fcpair); + if (!iszero(dsij1) || !iszero(dsij2)) { + force1 = dUdsij * dsij1; + force2 = dUdsij * dsij2; + for (m = 1; m <= 3; m++) { + delik[m] = arr2v(x, m, k) - arr2v(x, m, i); + deljk[m] = arr2v(x, m, k) - arr2v(x, m, j); + } + for (m = 1; m <= 3; m++) { + arr2v(f, m, i) = arr2v(f, m, i) + force1 * delik[m]; + arr2v(f, m, j) = arr2v(f, m, j) + force2 * deljk[m]; + arr2v(f, m, k) = + arr2v(f, m, k) - force1 * delik[m] - force2 * deljk[m]; + } + + // Tabulate per-atom virial as symmetrized stress tensor + + if (*vflag_atom != 0) { + fi[1] = force1 * delik[1]; + fi[2] = force1 * delik[2]; + fi[3] = force1 * delik[3]; + fj[1] = force2 * deljk[1]; + fj[2] = force2 * deljk[2]; + fj[3] = force2 * deljk[3]; + v[1] = -third * (delik[1] * fi[1] + deljk[1] * fj[1]); + v[2] = -third * (delik[2] * fi[2] + deljk[2] * fj[2]); + v[3] = -third * (delik[3] * fi[3] + deljk[3] * fj[3]); + v[4] = -sixth * (delik[1] * fi[2] + deljk[1] * fj[2] + + delik[2] * fi[1] + deljk[2] * fj[1]); + v[5] = -sixth * (delik[1] * fi[3] + deljk[1] * fj[3] + + delik[3] * fi[1] + deljk[3] * fj[1]); + v[6] = -sixth * (delik[2] * fi[3] + deljk[2] * fj[3] + + delik[3] * fi[2] + deljk[3] * fj[2]); + + arr2v(vatom, 1, i) = arr2v(vatom, 1, i) + v[1]; + arr2v(vatom, 2, i) = arr2v(vatom, 2, i) + v[2]; + arr2v(vatom, 3, i) = arr2v(vatom, 3, i) + v[3]; + arr2v(vatom, 4, i) = arr2v(vatom, 4, i) + v[4]; + arr2v(vatom, 5, i) = arr2v(vatom, 5, i) + v[5]; + arr2v(vatom, 6, i) = arr2v(vatom, 6, i) + v[6]; + arr2v(vatom, 1, j) = arr2v(vatom, 1, j) + v[1]; + arr2v(vatom, 2, j) = arr2v(vatom, 2, j) + v[2]; + arr2v(vatom, 3, j) = arr2v(vatom, 3, j) + v[3]; + arr2v(vatom, 4, j) = arr2v(vatom, 4, j) + v[4]; + arr2v(vatom, 5, j) = arr2v(vatom, 5, j) + v[5]; + arr2v(vatom, 6, j) = arr2v(vatom, 6, j) + v[6]; + arr2v(vatom, 1, k) = arr2v(vatom, 1, k) + v[1]; + arr2v(vatom, 2, k) = arr2v(vatom, 2, k) + v[2]; + arr2v(vatom, 3, k) = arr2v(vatom, 3, k) + v[3]; + arr2v(vatom, 4, k) = arr2v(vatom, 4, k) + v[4]; + arr2v(vatom, 5, k) = arr2v(vatom, 5, k) + v[5]; + arr2v(vatom, 6, k) = arr2v(vatom, 6, k) + v[6]; + } + } + } + // end of k loop + } } - -// else if elti=0, this is not a meam atom - } - } + // end of j loop + } + // else if elti=0, this is not a meam atom + } +} } \ No newline at end of file diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp index 0d445f84bc..0b20827420 100644 --- a/src/USER-MEAMC/meam_setup_done.cpp +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -4,23 +4,31 @@ extern "C" { void alloyparams(); void compute_pair_meam(); -double phi_meam(double r,int a, int b); +double phi_meam(double r, int a, int b); void compute_reference_density(); -void get_shpfcn(double *s /* s(3) */, lattice_t latt); -void get_tavref(double *t11av,double *t21av,double *t31av,double *t12av,double *t22av,double *t32av, double t11,double t21,double t31,double t12,double t22,double t32, double r,int a,int b,lattice_t latt); -void get_Zij(int *Zij, lattice_t latt); -void get_Zij2(int *Zij2, double *a, double*S, lattice_t latt,double cmin,double cmax); -void get_sijk(double C,int i,int j,int k, double *sijk); -void get_densref(double r,int a,int b,double *rho01,double *rho11,double *rho21,double *rho31, double *rho02,double *rho12,double *rho22,double *rho32); +void get_shpfcn(double* s /* s(3) */, lattice_t latt); +void get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, + double* t22av, double* t32av, double t11, double t21, + double t31, double t12, double t22, double t32, double r, int a, + int b, lattice_t latt); +void get_Zij(int* Zij, lattice_t latt); +void get_Zij2(int* Zij2, double* a, double* S, lattice_t latt, double cmin, + double cmax); +void get_sijk(double C, int i, int j, int k, double* sijk); +void get_densref(double r, int a, int b, double* rho01, double* rho11, + double* rho21, double* rho31, double* rho02, double* rho12, + double* rho22, double* rho32); double zbl(double r, int z1, int z2); -double erose(double r,double re,double alpha,double Ec,double repuls,double attrac,int form); +double erose(double r, double re, double alpha, double Ec, double repuls, + double attrac, int form); void interpolate_meam(int ind); double compute_phi(double rij, int elti, int eltj); // in meam_dens_init -void fcut(double xi, double *fc); +void fcut(double xi, double* fc); // in meam_dens_final -void G_gam(double Gamma,int ibar,double gsmooth_factor, double *G, int *errorflag); +void G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, + int* errorflag); // Declaration in pair_meam.h: // @@ -31,919 +39,1047 @@ void G_gam(double Gamma,int ibar,double gsmooth_factor, double *G, int *errorfla // meam_setup_done(&cutmax) // - void meam_setup_done_(double *cutmax) - { - int nv2, nv3, m, n, p; +void +meam_setup_done_(double* cutmax) +{ + int nv2, nv3, m, n, p; -// Force cutoff - meam_data.cutforce = meam_data.rc_meam; - meam_data.cutforcesq = meam_data.cutforce*meam_data.cutforce; + // Force cutoff + meam_data.cutforce = meam_data.rc_meam; + meam_data.cutforcesq = meam_data.cutforce * meam_data.cutforce; -// Pass cutoff back to calling program - *cutmax = meam_data.cutforce; + // Pass cutoff back to calling program + *cutmax = meam_data.cutforce; -// Augment t1 term - for (int i=1; i<=maxelt; i++) - meam_data.t1_meam[i] = meam_data.t1_meam[i] + meam_data.augt1 * 3.0/5.0 * meam_data.t3_meam[i]; + // Augment t1 term + for (int i = 1; i <= maxelt; i++) + meam_data.t1_meam[i] = + meam_data.t1_meam[i] + meam_data.augt1 * 3.0 / 5.0 * meam_data.t3_meam[i]; -// Compute off-diagonal alloy parameters - alloyparams(); + // Compute off-diagonal alloy parameters + alloyparams(); -// indices and factors for Voight notation - nv2 = 1; - nv3 = 1; - for(m = 1; m<=3; m++) { - for(n = m; n<=3; n++) { - meam_data.vind2D[m][n] = nv2; - meam_data.vind2D[n][m] = nv2; - nv2 = nv2+1; - for (p = n; p<=3; p++) { - meam_data.vind3D[m][n][p] = nv3; - meam_data.vind3D[m][p][n] = nv3; - meam_data.vind3D[n][m][p] = nv3; - meam_data.vind3D[n][p][m] = nv3; - meam_data.vind3D[p][m][n] = nv3; - meam_data.vind3D[p][n][m] = nv3; - nv3 = nv3+1; - } - } + // indices and factors for Voight notation + nv2 = 1; + nv3 = 1; + for (m = 1; m <= 3; m++) { + for (n = m; n <= 3; n++) { + meam_data.vind2D[m][n] = nv2; + meam_data.vind2D[n][m] = nv2; + nv2 = nv2 + 1; + for (p = n; p <= 3; p++) { + meam_data.vind3D[m][n][p] = nv3; + meam_data.vind3D[m][p][n] = nv3; + meam_data.vind3D[n][m][p] = nv3; + meam_data.vind3D[n][p][m] = nv3; + meam_data.vind3D[p][m][n] = nv3; + meam_data.vind3D[p][n][m] = nv3; + nv3 = nv3 + 1; } - - meam_data.v2D[1] = 1; - meam_data.v2D[2] = 2; - meam_data.v2D[3] = 2; - meam_data.v2D[4] = 1; - meam_data.v2D[5] = 2; - meam_data.v2D[6] = 1; - - meam_data.v3D[1] = 1; - meam_data.v3D[2] = 3; - meam_data.v3D[3] = 3; - meam_data.v3D[4] = 3; - meam_data.v3D[5] = 6; - meam_data.v3D[6] = 3; - meam_data.v3D[7] = 1; - meam_data.v3D[8] = 3; - meam_data.v3D[9] = 3; - meam_data.v3D[10] = 1; - - nv2 = 1; - for(m = 1; m<=meam_data.neltypes; m++) { - for(n = m; n<=meam_data.neltypes; n++) { - meam_data.eltind[m][n] = nv2; - meam_data.eltind[n][m] = nv2; - nv2 = nv2+1; - } - } - -// Compute background densities for reference structure - compute_reference_density(); - -// Compute pair potentials and setup arrays for interpolation - meam_data.nr = 1000; - meam_data.dr = 1.1*meam_data.rc_meam/meam_data.nr; - compute_pair_meam(); } + } -//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + meam_data.v2D[1] = 1; + meam_data.v2D[2] = 2; + meam_data.v2D[3] = 2; + meam_data.v2D[4] = 1; + meam_data.v2D[5] = 2; + meam_data.v2D[6] = 1; + + meam_data.v3D[1] = 1; + meam_data.v3D[2] = 3; + meam_data.v3D[3] = 3; + meam_data.v3D[4] = 3; + meam_data.v3D[5] = 6; + meam_data.v3D[6] = 3; + meam_data.v3D[7] = 1; + meam_data.v3D[8] = 3; + meam_data.v3D[9] = 3; + meam_data.v3D[10] = 1; + + nv2 = 1; + for (m = 1; m <= meam_data.neltypes; m++) { + for (n = m; n <= meam_data.neltypes; n++) { + meam_data.eltind[m][n] = nv2; + meam_data.eltind[n][m] = nv2; + nv2 = nv2 + 1; + } + } + + // Compute background densities for reference structure + compute_reference_density(); + + // Compute pair potentials and setup arrays for interpolation + meam_data.nr = 1000; + meam_data.dr = 1.1 * meam_data.rc_meam / meam_data.nr; + compute_pair_meam(); +} + +// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc // Fill off-diagonal alloy parameters - void alloyparams(void) - { - - int i,j,k; - double eb; +void +alloyparams(void) +{ -// Loop over pairs - for(i = 1; i<=meam_data.neltypes; i++) { - for(j = 1;i<=meam_data.neltypes; i++) { -// Treat off-diagonal pairs -// If i>j, set all equal to i j) { - meam_data.re_meam[i][j] = meam_data.re_meam[j][i]; - meam_data.Ec_meam[i][j] = meam_data.Ec_meam[j][i]; - meam_data.alpha_meam[i][j] = meam_data.alpha_meam[j][i]; - meam_data.lattce_meam[i][j] = meam_data.lattce_meam[j][i]; - meam_data.nn2_meam[i][j] = meam_data.nn2_meam[j][i]; -// If i i) { - if (iszero(meam_data.Ec_meam[i][j])) { - if (meam_data.lattce_meam[i][j]==L12) - meam_data.Ec_meam[i][j] = (3*meam_data.Ec_meam[i][i]+meam_data.Ec_meam[j][j])/4.0 - meam_data.delta_meam[i][j]; - else if (meam_data.lattce_meam[i][j]==C11) { - if (meam_data.lattce_meam[i][i]==DIA) - meam_data.Ec_meam[i][j] = (2*meam_data.Ec_meam[i][i]+meam_data.Ec_meam[j][j])/3.0 - meam_data.delta_meam[i][j]; - else - meam_data.Ec_meam[i][j] = (meam_data.Ec_meam[i][i]+2*meam_data.Ec_meam[j][j])/3.0 - meam_data.delta_meam[i][j]; - } else - meam_data.Ec_meam[i][j] = (meam_data.Ec_meam[i][i]+meam_data.Ec_meam[j][j])/2.0 - meam_data.delta_meam[i][j]; - } - if (iszero(meam_data.alpha_meam[i][j])) - meam_data.alpha_meam[i][j] = (meam_data.alpha_meam[i][i]+meam_data.alpha_meam[j][j])/2.0; - if (iszero(meam_data.re_meam[i][j])) - meam_data.re_meam[i][j] = (meam_data.re_meam[i][i]+meam_data.re_meam[j][j])/2.0; - } - } - } + int i, j, k; + double eb; -// Cmin[i][k][j] is symmetric in i-j, but not k. For all triplets -// where i>j, set equal to the iebound, -// atom k definitely lies outside the screening function ellipse (so -// there is no need to calculate its effects). Here, compute it for all -// triplets [i][j][k] so that ebound[i][j] is the maximized over k - for(i = 2;i<=meam_data.neltypes;i++){ - for(j = 1;j<=meam_data.neltypes;j++){ - for(k = 1;k<=meam_data.neltypes;k++){ - eb = (meam_data.Cmax_meam[i][j][k]*meam_data.Cmax_meam[i][j][k]) / (4.0*(meam_data.Cmax_meam[i][j][k]-1.0)); - meam_data.ebound_meam[i][j] = max(meam_data.ebound_meam[i][j],eb); - } + // Loop over pairs + for (i = 1; i <= meam_data.neltypes; i++) { + for (j = 1; i <= meam_data.neltypes; i++) { + // Treat off-diagonal pairs + // If i>j, set all equal to i j) { + meam_data.re_meam[i][j] = meam_data.re_meam[j][i]; + meam_data.Ec_meam[i][j] = meam_data.Ec_meam[j][i]; + meam_data.alpha_meam[i][j] = meam_data.alpha_meam[j][i]; + meam_data.lattce_meam[i][j] = meam_data.lattce_meam[j][i]; + meam_data.nn2_meam[i][j] = meam_data.nn2_meam[j][i]; + // If i i) { + if (iszero(meam_data.Ec_meam[i][j])) { + if (meam_data.lattce_meam[i][j] == L12) + meam_data.Ec_meam[i][j] = + (3 * meam_data.Ec_meam[i][i] + meam_data.Ec_meam[j][j]) / 4.0 - + meam_data.delta_meam[i][j]; + else if (meam_data.lattce_meam[i][j] == C11) { + if (meam_data.lattce_meam[i][i] == DIA) + meam_data.Ec_meam[i][j] = + (2 * meam_data.Ec_meam[i][i] + meam_data.Ec_meam[j][j]) / 3.0 - + meam_data.delta_meam[i][j]; + else + meam_data.Ec_meam[i][j] = + (meam_data.Ec_meam[i][i] + 2 * meam_data.Ec_meam[j][j]) / 3.0 - + meam_data.delta_meam[i][j]; + } else + meam_data.Ec_meam[i][j] = + (meam_data.Ec_meam[i][i] + meam_data.Ec_meam[j][j]) / 2.0 - + meam_data.delta_meam[i][j]; } + if (iszero(meam_data.alpha_meam[i][j])) + meam_data.alpha_meam[i][j] = + (meam_data.alpha_meam[i][i] + meam_data.alpha_meam[j][j]) / 2.0; + if (iszero(meam_data.re_meam[i][j])) + meam_data.re_meam[i][j] = + (meam_data.re_meam[i][i] + meam_data.re_meam[j][j]) / 2.0; } } + } + + // Cmin[i][k][j] is symmetric in i-j, but not k. For all triplets + // where i>j, set equal to the iebound, + // atom k definitely lies outside the screening function ellipse (so + // there is no need to calculate its effects). Here, compute it for all + // triplets [i][j][k] so that ebound[i][j] is the maximized over k + for (i = 2; i <= meam_data.neltypes; i++) { + for (j = 1; j <= meam_data.neltypes; j++) { + for (k = 1; k <= meam_data.neltypes; k++) { + eb = (meam_data.Cmax_meam[i][j][k] * meam_data.Cmax_meam[i][j][k]) / + (4.0 * (meam_data.Cmax_meam[i][j][k] - 1.0)); + meam_data.ebound_meam[i][j] = max(meam_data.ebound_meam[i][j], eb); + } + } + } +} //----------------------------------------------------------------------- // compute MEAM pair potential for each pair of element types // - void compute_pair_meam(void) - { +void +compute_pair_meam(void) +{ - double r/*ununsed:, temp*/; - int j,a,b,nv2; - double astar,frac,phizbl; - int n,nmax,Z1,Z2; - double arat,rarat,scrn,scrn2; - double phiaa,phibb/*unused:,phitmp*/; - double C,s111,s112,s221,S11,S22; + double r /*ununsed:, temp*/; + int j, a, b, nv2; + double astar, frac, phizbl; + int n, nmax, Z1, Z2; + double arat, rarat, scrn, scrn2; + double phiaa, phibb /*unused:,phitmp*/; + double C, s111, s112, s221, S11, S22; -// check for previously allocated arrays and free them - if(allocated(meam_data.phir)) deallocate(meam_data.phir); - if(allocated(meam_data.phirar)) deallocate(meam_data.phirar); - if(allocated(meam_data.phirar1)) deallocate(meam_data.phirar1); - if(allocated(meam_data.phirar2)) deallocate(meam_data.phirar2); - if(allocated(meam_data.phirar3)) deallocate(meam_data.phirar3); - if(allocated(meam_data.phirar4)) deallocate(meam_data.phirar4); - if(allocated(meam_data.phirar5)) deallocate(meam_data.phirar5); - if(allocated(meam_data.phirar6)) deallocate(meam_data.phirar6); + // check for previously allocated arrays and free them + if (allocated(meam_data.phir)) + deallocate(meam_data.phir); + if (allocated(meam_data.phirar)) + deallocate(meam_data.phirar); + if (allocated(meam_data.phirar1)) + deallocate(meam_data.phirar1); + if (allocated(meam_data.phirar2)) + deallocate(meam_data.phirar2); + if (allocated(meam_data.phirar3)) + deallocate(meam_data.phirar3); + if (allocated(meam_data.phirar4)) + deallocate(meam_data.phirar4); + if (allocated(meam_data.phirar5)) + deallocate(meam_data.phirar5); + if (allocated(meam_data.phirar6)) + deallocate(meam_data.phirar6); -// allocate memory for array that defines the potential - allocate_2d(meam_data.phir,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2); + // allocate memory for array that defines the potential + allocate_2d(meam_data.phir, meam_data.nr, + (meam_data.neltypes * (meam_data.neltypes + 1)) / 2); -// allocate coeff memory + // allocate coeff memory - allocate_2d(meam_data.phirar,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2); - allocate_2d(meam_data.phirar1,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2); - allocate_2d(meam_data.phirar2,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2); - allocate_2d(meam_data.phirar3,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2); - allocate_2d(meam_data.phirar4,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2); - allocate_2d(meam_data.phirar5,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2); - allocate_2d(meam_data.phirar6,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2); + allocate_2d(meam_data.phirar, meam_data.nr, + (meam_data.neltypes * (meam_data.neltypes + 1)) / 2); + allocate_2d(meam_data.phirar1, meam_data.nr, + (meam_data.neltypes * (meam_data.neltypes + 1)) / 2); + allocate_2d(meam_data.phirar2, meam_data.nr, + (meam_data.neltypes * (meam_data.neltypes + 1)) / 2); + allocate_2d(meam_data.phirar3, meam_data.nr, + (meam_data.neltypes * (meam_data.neltypes + 1)) / 2); + allocate_2d(meam_data.phirar4, meam_data.nr, + (meam_data.neltypes * (meam_data.neltypes + 1)) / 2); + allocate_2d(meam_data.phirar5, meam_data.nr, + (meam_data.neltypes * (meam_data.neltypes + 1)) / 2); + allocate_2d(meam_data.phirar6, meam_data.nr, + (meam_data.neltypes * (meam_data.neltypes + 1)) / 2); -// loop over pairs of element types - nv2 = 0; - for(a=1; a<=meam_data.neltypes; a++) { - for(b=a; b<=meam_data.neltypes; b++) { - nv2 = nv2 + 1; + // loop over pairs of element types + nv2 = 0; + for (a = 1; a <= meam_data.neltypes; a++) { + for (b = a; b <= meam_data.neltypes; b++) { + nv2 = nv2 + 1; -// loop over r values and compute - for(j=1; j<=meam_data.nr; j++) { - r = (j-1)*meam_data.dr; + // loop over r values and compute + for (j = 1; j <= meam_data.nr; j++) { + r = (j - 1) * meam_data.dr; - arr2(meam_data.phir,j,nv2) = phi_meam(r,a,b); + arr2(meam_data.phir, j, nv2) = phi_meam(r, a, b); -// if using second-nearest neighbor, solve recursive problem -// (see Lee and Baskes, PRB 62(13):8564 eqn.(21)) - if (meam_data.nn2_meam[a][b]==1) { - get_Zij(&Z1,meam_data.lattce_meam[a][b]); - get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[a][b],meam_data.Cmin_meam[a][a][b],meam_data.Cmax_meam[a][a][b]); + // if using second-nearest neighbor, solve recursive problem + // (see Lee and Baskes, PRB 62(13):8564 eqn.(21)) + if (meam_data.nn2_meam[a][b] == 1) { + get_Zij(&Z1, meam_data.lattce_meam[a][b]); + get_Zij2(&Z2, &arat, &scrn, meam_data.lattce_meam[a][b], + meam_data.Cmin_meam[a][a][b], meam_data.Cmax_meam[a][a][b]); -// The B1, B2, and L12 cases with NN2 have a trick to them; we need to -// compute the contributions from second nearest neighbors, like a-a -// pairs, but need to include NN2 contributions to those pairs as -// well. - if (meam_data.lattce_meam[a][b]==B1 || meam_data.lattce_meam[a][b]==B2 || meam_data.lattce_meam[a][b]==L12) { - rarat = r*arat; + // The B1, B2, and L12 cases with NN2 have a trick to them; we + // need to + // compute the contributions from second nearest neighbors, like + // a-a + // pairs, but need to include NN2 contributions to those pairs as + // well. + if (meam_data.lattce_meam[a][b] == B1 || + meam_data.lattce_meam[a][b] == B2 || + meam_data.lattce_meam[a][b] == L12) { + rarat = r * arat; -// phi_aa - phiaa = phi_meam(rarat,a,a); - get_Zij(&Z1,meam_data.lattce_meam[a][a]); - get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[a][a], meam_data.Cmin_meam[a][a][a],meam_data.Cmax_meam[a][a][a]); - nmax = 10; - if (scrn > 0.0) { - for(n=1; n<=nmax; n++) { - phiaa = phiaa + pow((-Z2*scrn/Z1),n) * phi_meam(rarat*pow(arat,n),a,a); - } - } - -// phi_bb - phibb = phi_meam(rarat,b,b); - get_Zij(&Z1,meam_data.lattce_meam[b][b]); - get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[b][b], meam_data.Cmin_meam[b][b][b],meam_data.Cmax_meam[b][b][b]); - nmax = 10; - if (scrn > 0.0) { - for(n=1; n<=nmax; n++) { - phibb = phibb + pow((-Z2*scrn/Z1),n) * phi_meam(rarat*pow(arat,n),b,b); - } - } - - if (meam_data.lattce_meam[a][b]==B1 || meam_data.lattce_meam[a][b]==B2) { -// Add contributions to the B1 or B2 potential - get_Zij(&Z1,meam_data.lattce_meam[a][b]); - get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[a][b], meam_data.Cmin_meam[a][a][b],meam_data.Cmax_meam[a][a][b]); - arr2(meam_data.phir,j,nv2) = arr2(meam_data.phir,j,nv2) - Z2*scrn/(2*Z1) * phiaa; - get_Zij2(&Z2,&arat,&scrn2,meam_data.lattce_meam[a][b], meam_data.Cmin_meam[b][b][a],meam_data.Cmax_meam[b][b][a]); - arr2(meam_data.phir,j,nv2) = arr2(meam_data.phir,j,nv2) - Z2*scrn2/(2*Z1) * phibb; - - } else if (meam_data.lattce_meam[a][b]==L12) { -// The L12 case has one last trick; we have to be careful to compute -// the correct screening between 2nd-neighbor pairs. 1-1 -// second-neighbor pairs are screened by 2 type 1 atoms and two type -// 2 atoms. 2-2 second-neighbor pairs are screened by 4 type 1 -// atoms. - C = 1.0; - get_sijk(C,a,a,a,&s111); - get_sijk(C,a,a,b,&s112); - get_sijk(C,b,b,a,&s221); - S11 = s111 * s111 * s112 * s112; - S22 = pow(s221,4); - arr2(meam_data.phir,j,nv2) = arr2(meam_data.phir,j,nv2) - 0.75*S11*phiaa - 0.25*S22*phibb; - - } - - } else { - nmax = 10; - for(n=1; n<=nmax; n++) { - arr2(meam_data.phir,j,nv2) = arr2(meam_data.phir,j,nv2) + pow((-Z2*scrn/Z1),n) * phi_meam(r*pow(arat,n),a,b); - } - } - - } - -// For Zbl potential: -// if astar <= -3 -// potential is zbl potential -// else if -3 < astar < -1 -// potential is linear combination with zbl potential -// endif - if (meam_data.zbl_meam[a][b]==1) { - astar = meam_data.alpha_meam[a][b] * (r/meam_data.re_meam[a][b] - 1.0); - if (astar <= -3.0) - arr2(meam_data.phir,j,nv2) = zbl(r,meam_data.ielt_meam[a],meam_data.ielt_meam[b]); - else if (astar > -3.0 && astar < -1.0) { - fcut(1-(astar+1.0)/(-3.0+1.0),&frac); - phizbl = zbl(r,meam_data.ielt_meam[a],meam_data.ielt_meam[b]); - arr2(meam_data.phir,j,nv2) = frac*arr2(meam_data.phir,j,nv2) + (1-frac)*phizbl; + // phi_aa + phiaa = phi_meam(rarat, a, a); + get_Zij(&Z1, meam_data.lattce_meam[a][a]); + get_Zij2(&Z2, &arat, &scrn, meam_data.lattce_meam[a][a], + meam_data.Cmin_meam[a][a][a], + meam_data.Cmax_meam[a][a][a]); + nmax = 10; + if (scrn > 0.0) { + for (n = 1; n <= nmax; n++) { + phiaa = phiaa + + pow((-Z2 * scrn / Z1), n) * + phi_meam(rarat * pow(arat, n), a, a); } } + // phi_bb + phibb = phi_meam(rarat, b, b); + get_Zij(&Z1, meam_data.lattce_meam[b][b]); + get_Zij2(&Z2, &arat, &scrn, meam_data.lattce_meam[b][b], + meam_data.Cmin_meam[b][b][b], + meam_data.Cmax_meam[b][b][b]); + nmax = 10; + if (scrn > 0.0) { + for (n = 1; n <= nmax; n++) { + phibb = phibb + + pow((-Z2 * scrn / Z1), n) * + phi_meam(rarat * pow(arat, n), b, b); + } + } + + if (meam_data.lattce_meam[a][b] == B1 || + meam_data.lattce_meam[a][b] == B2) { + // Add contributions to the B1 or B2 potential + get_Zij(&Z1, meam_data.lattce_meam[a][b]); + get_Zij2(&Z2, &arat, &scrn, meam_data.lattce_meam[a][b], + meam_data.Cmin_meam[a][a][b], + meam_data.Cmax_meam[a][a][b]); + arr2(meam_data.phir, j, nv2) = + arr2(meam_data.phir, j, nv2) - Z2 * scrn / (2 * Z1) * phiaa; + get_Zij2(&Z2, &arat, &scrn2, meam_data.lattce_meam[a][b], + meam_data.Cmin_meam[b][b][a], + meam_data.Cmax_meam[b][b][a]); + arr2(meam_data.phir, j, nv2) = + arr2(meam_data.phir, j, nv2) - Z2 * scrn2 / (2 * Z1) * phibb; + + } else if (meam_data.lattce_meam[a][b] == L12) { + // The L12 case has one last trick; we have to be careful to + // compute + // the correct screening between 2nd-neighbor pairs. 1-1 + // second-neighbor pairs are screened by 2 type 1 atoms and + // two type + // 2 atoms. 2-2 second-neighbor pairs are screened by 4 type + // 1 + // atoms. + C = 1.0; + get_sijk(C, a, a, a, &s111); + get_sijk(C, a, a, b, &s112); + get_sijk(C, b, b, a, &s221); + S11 = s111 * s111 * s112 * s112; + S22 = pow(s221, 4); + arr2(meam_data.phir, j, nv2) = arr2(meam_data.phir, j, nv2) - + 0.75 * S11 * phiaa - + 0.25 * S22 * phibb; + } + + } else { + nmax = 10; + for (n = 1; n <= nmax; n++) { + arr2(meam_data.phir, j, nv2) = + arr2(meam_data.phir, j, nv2) + + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b); + } } + } -// call interpolation - interpolate_meam(nv2); - + // For Zbl potential: + // if astar <= -3 + // potential is zbl potential + // else if -3 < astar < -1 + // potential is linear combination with zbl potential + // endif + if (meam_data.zbl_meam[a][b] == 1) { + astar = + meam_data.alpha_meam[a][b] * (r / meam_data.re_meam[a][b] - 1.0); + if (astar <= -3.0) + arr2(meam_data.phir, j, nv2) = + zbl(r, meam_data.ielt_meam[a], meam_data.ielt_meam[b]); + else if (astar > -3.0 && astar < -1.0) { + fcut(1 - (astar + 1.0) / (-3.0 + 1.0), &frac); + phizbl = zbl(r, meam_data.ielt_meam[a], meam_data.ielt_meam[b]); + arr2(meam_data.phir, j, nv2) = + frac * arr2(meam_data.phir, j, nv2) + (1 - frac) * phizbl; + } } } - } - + // call interpolation + interpolate_meam(nv2); + } + } +} //----------------------------------------------------------------------c // Compute MEAM pair potential for distance r, element types a and b // - double phi_meam(double r,int a, int b) - { - /*unused:double a1,a2,a12;*/ - double t11av,t21av,t31av,t12av,t22av,t32av; - double G1,G2,s1[3+1],s2[3+1]/*,s12[3+1]*/,rho0_1,rho0_2; - double Gam1,Gam2,Z1,Z2; - double rhobar1,rhobar2,F1,F2; - double rho01,rho11,rho21,rho31; - double rho02,rho12,rho22,rho32; - double scalfac,phiaa,phibb; - double Eu; - double arat,scrn/*unused:,scrn2*/; - int Z12, errorflag; - int n,nmax,Z1nn,Z2nn; - lattice_t latta/*unused:,lattb*/; - double rho_bkgd1, rho_bkgd2; - - double phi_m = 0.0; +double +phi_meam(double r, int a, int b) +{ + /*unused:double a1,a2,a12;*/ + double t11av, t21av, t31av, t12av, t22av, t32av; + double G1, G2, s1[3 + 1], s2[3 + 1] /*,s12[3+1]*/, rho0_1, rho0_2; + double Gam1, Gam2, Z1, Z2; + double rhobar1, rhobar2, F1, F2; + double rho01, rho11, rho21, rho31; + double rho02, rho12, rho22, rho32; + double scalfac, phiaa, phibb; + double Eu; + double arat, scrn /*unused:,scrn2*/; + int Z12, errorflag; + int n, nmax, Z1nn, Z2nn; + lattice_t latta /*unused:,lattb*/; + double rho_bkgd1, rho_bkgd2; -// Equation numbers below refer to: -// I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615 + double phi_m = 0.0; -// get number of neighbors in the reference structure -// Nref[i][j] = # of i's neighbors of type j - get_Zij(&Z12,meam_data.lattce_meam[a][b]); + // Equation numbers below refer to: + // I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615 - get_densref(r,a,b,&rho01,&rho11,&rho21,&rho31,&rho02,&rho12,&rho22,&rho32); + // get number of neighbors in the reference structure + // Nref[i][j] = # of i's neighbors of type j + get_Zij(&Z12, meam_data.lattce_meam[a][b]); -// if densities are too small, numerical problems may result; just return zero - if (rho01<=1e-14 && rho02<=1e-14) - return 0.0; + get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, + &rho32); -// calculate average weighting factors for the reference structure - if (meam_data.lattce_meam[a][b]==C11) { - if (meam_data.ialloy==2) { - t11av = meam_data.t1_meam[a]; - t12av = meam_data.t1_meam[b]; - t21av = meam_data.t2_meam[a]; - t22av = meam_data.t2_meam[b]; - t31av = meam_data.t3_meam[a]; - t32av = meam_data.t3_meam[b]; - } else { - scalfac = 1.0/(rho01+rho02); - t11av = scalfac*(meam_data.t1_meam[a]*rho01 + meam_data.t1_meam[b]*rho02); - t12av = t11av; - t21av = scalfac*(meam_data.t2_meam[a]*rho01 + meam_data.t2_meam[b]*rho02); - t22av = t21av; - t31av = scalfac*(meam_data.t3_meam[a]*rho01 + meam_data.t3_meam[b]*rho02); - t32av = t31av; - } - } else { -// average weighting factors for the reference structure, eqn. I.8 - get_tavref(&t11av,&t21av,&t31av,&t12av,&t22av,&t32av, meam_data.t1_meam[a],meam_data.t2_meam[a],meam_data.t3_meam[a], meam_data.t1_meam[b],meam_data.t2_meam[b],meam_data.t3_meam[b], r,a,b,meam_data.lattce_meam[a][b]); - } + // if densities are too small, numerical problems may result; just return zero + if (rho01 <= 1e-14 && rho02 <= 1e-14) + return 0.0; -// for c11b structure, calculate background electron densities - if (meam_data.lattce_meam[a][b]==C11) { - latta = meam_data.lattce_meam[a][a]; - if (latta==DIA) { - rhobar1 = pow(((Z12/2)*(rho02+rho01)),2) + t11av*pow((rho12-rho11),2) + t21av/6.0*pow(rho22+rho21,2)+ 121.0/40.0*t31av*pow((rho32-rho31),2); - rhobar1 = sqrt(rhobar1); - rhobar2 = pow(Z12*rho01,2) + 2.0/3.0*t21av*pow(rho21,2); - rhobar2 = sqrt(rhobar2); - } else { - rhobar2 = pow(((Z12/2)*(rho01+rho02)),2) + t12av*pow((rho11-rho12),2) + t22av/6.0*pow(rho21+rho22,2) + 121.0/40.0*t32av*pow((rho31-rho32),2); - rhobar2 = sqrt(rhobar2); - rhobar1 = pow(Z12*rho02,2) + 2.0/3.0*t22av*pow(rho22,2); - rhobar1 = sqrt(rhobar1); - } - } else { -// for other structures, use formalism developed in Huang's paper -// -// composition-dependent scaling, equation I.7 -// If using mixing rule for t, apply to reference structure; else -// use precomputed values - if (meam_data.mix_ref_t==1) { - Z1 = meam_data.Z_meam[a]; - Z2 = meam_data.Z_meam[b]; - if (meam_data.ibar_meam[a]<=0) - G1 = 1.0; - else { - get_shpfcn(s1,meam_data.lattce_meam[a][a]); - Gam1 = (s1[1]*t11av+s1[2]*t21av+s1[3]*t31av)/(Z1*Z1); - G_gam(Gam1,meam_data.ibar_meam[a],meam_data.gsmooth_factor,&G1,&errorflag); - } - if (meam_data.ibar_meam[b]<=0) - G2 = 1.0; - else { - get_shpfcn(s2,meam_data.lattce_meam[b][b]); - Gam2 = (s2[1]*t12av+s2[2]*t22av+s2[3]*t32av)/(Z2*Z2); - G_gam(Gam2,meam_data.ibar_meam[b],meam_data.gsmooth_factor,&G2,&errorflag); - } - rho0_1 = meam_data.rho0_meam[a]*Z1*G1; - rho0_2 = meam_data.rho0_meam[b]*Z2*G2; - } - Gam1 = (t11av*rho11+t21av*rho21+t31av*rho31); - if (rho01 < 1.0e-14) - Gam1 = 0.0; - else - Gam1 = Gam1/(rho01*rho01); - - Gam2 = (t12av*rho12+t22av*rho22+t32av*rho32); - if (rho02 < 1.0e-14) - Gam2 = 0.0; - else - Gam2 = Gam2/(rho02*rho02); - - G_gam(Gam1,meam_data.ibar_meam[a],meam_data.gsmooth_factor,&G1,&errorflag); - G_gam(Gam2,meam_data.ibar_meam[b],meam_data.gsmooth_factor,&G2,&errorflag); - if (meam_data.mix_ref_t==1) { - rho_bkgd1 = rho0_1; - rho_bkgd2 = rho0_2; - } else { - if (meam_data.bkgd_dyn==1) { - rho_bkgd1 = meam_data.rho0_meam[a]*meam_data.Z_meam[a]; - rho_bkgd2 = meam_data.rho0_meam[b]*meam_data.Z_meam[b]; - } else { - rho_bkgd1 = meam_data.rho_ref_meam[a]; - rho_bkgd2 = meam_data.rho_ref_meam[b]; - } - } - rhobar1 = rho01/rho_bkgd1*G1; - rhobar2 = rho02/rho_bkgd2*G2; - - } - -// compute embedding functions, eqn I.5 - if (iszero(rhobar1)) - F1 = 0.0; - else { - if (meam_data.emb_lin_neg==1 && rhobar1<=0) - F1 = -meam_data.A_meam[a]*meam_data.Ec_meam[a][a]*rhobar1; - else - F1 = meam_data.A_meam[a]*meam_data.Ec_meam[a][a]*rhobar1*log(rhobar1); - } - if (iszero(rhobar2)) - F2 = 0.0; - else { - if (meam_data.emb_lin_neg==1 && rhobar2<=0) - F2 = -meam_data.A_meam[b]*meam_data.Ec_meam[b][b]*rhobar2; - else - F2 = meam_data.A_meam[b]*meam_data.Ec_meam[b][b]*rhobar2*log(rhobar2); - } - -// compute Rose function, I.16 - Eu = erose(r,meam_data.re_meam[a][b],meam_data.alpha_meam[a][b], meam_data.Ec_meam[a][b],meam_data.repuls_meam[a][b],meam_data.attrac_meam[a][b],meam_data.erose_form); - -// calculate the pair energy - if (meam_data.lattce_meam[a][b]==C11) { - latta = meam_data.lattce_meam[a][a]; - if (latta==DIA) { - phiaa = phi_meam(r,a,a); - phi_m = (3*Eu - F2 - 2*F1 - 5*phiaa)/Z12; - } else { - phibb = phi_meam(r,b,b); - phi_m = (3*Eu - F1 - 2*F2 - 5*phibb)/Z12; - } - } else if (meam_data.lattce_meam[a][b]==L12) { - phiaa = phi_meam(r,a,a); -// account for second neighbor a-a potential here... - get_Zij(&Z1nn,meam_data.lattce_meam[a][a]); - get_Zij2(&Z2nn,&arat,&scrn,meam_data.lattce_meam[a][a],meam_data.Cmin_meam[a][a][a],meam_data.Cmax_meam[a][a][a]); - nmax = 10; - if (scrn > 0.0) { - for(n=1; n<=nmax; n++) { - phiaa = phiaa + pow((-Z2nn*scrn/Z1nn),n) * phi_meam(r*pow(arat,n),a,a); - } - } - phi_m = Eu/3.0 - F1/4.0 - F2/12.0 - phiaa; - - } else { -// -// potential is computed from Rose function and embedding energy - phi_m = (2*Eu - F1 - F2)/Z12; -// - } - -// if r = 0, just return 0 - if (iszero(r)) { - phi_m = 0.0; - } - - return phi_m; + // calculate average weighting factors for the reference structure + if (meam_data.lattce_meam[a][b] == C11) { + if (meam_data.ialloy == 2) { + t11av = meam_data.t1_meam[a]; + t12av = meam_data.t1_meam[b]; + t21av = meam_data.t2_meam[a]; + t22av = meam_data.t2_meam[b]; + t31av = meam_data.t3_meam[a]; + t32av = meam_data.t3_meam[b]; + } else { + scalfac = 1.0 / (rho01 + rho02); + t11av = + scalfac * (meam_data.t1_meam[a] * rho01 + meam_data.t1_meam[b] * rho02); + t12av = t11av; + t21av = + scalfac * (meam_data.t2_meam[a] * rho01 + meam_data.t2_meam[b] * rho02); + t22av = t21av; + t31av = + scalfac * (meam_data.t3_meam[a] * rho01 + meam_data.t3_meam[b] * rho02); + t32av = t31av; } + } else { + // average weighting factors for the reference structure, eqn. I.8 + get_tavref(&t11av, &t21av, &t31av, &t12av, &t22av, &t32av, + meam_data.t1_meam[a], meam_data.t2_meam[a], meam_data.t3_meam[a], + meam_data.t1_meam[b], meam_data.t2_meam[b], meam_data.t3_meam[b], + r, a, b, meam_data.lattce_meam[a][b]); + } + + // for c11b structure, calculate background electron densities + if (meam_data.lattce_meam[a][b] == C11) { + latta = meam_data.lattce_meam[a][a]; + if (latta == DIA) { + rhobar1 = pow(((Z12 / 2) * (rho02 + rho01)), 2) + + t11av * pow((rho12 - rho11), 2) + + t21av / 6.0 * pow(rho22 + rho21, 2) + + 121.0 / 40.0 * t31av * pow((rho32 - rho31), 2); + rhobar1 = sqrt(rhobar1); + rhobar2 = pow(Z12 * rho01, 2) + 2.0 / 3.0 * t21av * pow(rho21, 2); + rhobar2 = sqrt(rhobar2); + } else { + rhobar2 = pow(((Z12 / 2) * (rho01 + rho02)), 2) + + t12av * pow((rho11 - rho12), 2) + + t22av / 6.0 * pow(rho21 + rho22, 2) + + 121.0 / 40.0 * t32av * pow((rho31 - rho32), 2); + rhobar2 = sqrt(rhobar2); + rhobar1 = pow(Z12 * rho02, 2) + 2.0 / 3.0 * t22av * pow(rho22, 2); + rhobar1 = sqrt(rhobar1); + } + } else { + // for other structures, use formalism developed in Huang's paper + // + // composition-dependent scaling, equation I.7 + // If using mixing rule for t, apply to reference structure; else + // use precomputed values + if (meam_data.mix_ref_t == 1) { + Z1 = meam_data.Z_meam[a]; + Z2 = meam_data.Z_meam[b]; + if (meam_data.ibar_meam[a] <= 0) + G1 = 1.0; + else { + get_shpfcn(s1, meam_data.lattce_meam[a][a]); + Gam1 = (s1[1] * t11av + s1[2] * t21av + s1[3] * t31av) / (Z1 * Z1); + G_gam(Gam1, meam_data.ibar_meam[a], meam_data.gsmooth_factor, &G1, + &errorflag); + } + if (meam_data.ibar_meam[b] <= 0) + G2 = 1.0; + else { + get_shpfcn(s2, meam_data.lattce_meam[b][b]); + Gam2 = (s2[1] * t12av + s2[2] * t22av + s2[3] * t32av) / (Z2 * Z2); + G_gam(Gam2, meam_data.ibar_meam[b], meam_data.gsmooth_factor, &G2, + &errorflag); + } + rho0_1 = meam_data.rho0_meam[a] * Z1 * G1; + rho0_2 = meam_data.rho0_meam[b] * Z2 * G2; + } + Gam1 = (t11av * rho11 + t21av * rho21 + t31av * rho31); + if (rho01 < 1.0e-14) + Gam1 = 0.0; + else + Gam1 = Gam1 / (rho01 * rho01); + + Gam2 = (t12av * rho12 + t22av * rho22 + t32av * rho32); + if (rho02 < 1.0e-14) + Gam2 = 0.0; + else + Gam2 = Gam2 / (rho02 * rho02); + + G_gam(Gam1, meam_data.ibar_meam[a], meam_data.gsmooth_factor, &G1, + &errorflag); + G_gam(Gam2, meam_data.ibar_meam[b], meam_data.gsmooth_factor, &G2, + &errorflag); + if (meam_data.mix_ref_t == 1) { + rho_bkgd1 = rho0_1; + rho_bkgd2 = rho0_2; + } else { + if (meam_data.bkgd_dyn == 1) { + rho_bkgd1 = meam_data.rho0_meam[a] * meam_data.Z_meam[a]; + rho_bkgd2 = meam_data.rho0_meam[b] * meam_data.Z_meam[b]; + } else { + rho_bkgd1 = meam_data.rho_ref_meam[a]; + rho_bkgd2 = meam_data.rho_ref_meam[b]; + } + } + rhobar1 = rho01 / rho_bkgd1 * G1; + rhobar2 = rho02 / rho_bkgd2 * G2; + } + + // compute embedding functions, eqn I.5 + if (iszero(rhobar1)) + F1 = 0.0; + else { + if (meam_data.emb_lin_neg == 1 && rhobar1 <= 0) + F1 = -meam_data.A_meam[a] * meam_data.Ec_meam[a][a] * rhobar1; + else + F1 = + meam_data.A_meam[a] * meam_data.Ec_meam[a][a] * rhobar1 * log(rhobar1); + } + if (iszero(rhobar2)) + F2 = 0.0; + else { + if (meam_data.emb_lin_neg == 1 && rhobar2 <= 0) + F2 = -meam_data.A_meam[b] * meam_data.Ec_meam[b][b] * rhobar2; + else + F2 = + meam_data.A_meam[b] * meam_data.Ec_meam[b][b] * rhobar2 * log(rhobar2); + } + + // compute Rose function, I.16 + Eu = erose(r, meam_data.re_meam[a][b], meam_data.alpha_meam[a][b], + meam_data.Ec_meam[a][b], meam_data.repuls_meam[a][b], + meam_data.attrac_meam[a][b], meam_data.erose_form); + + // calculate the pair energy + if (meam_data.lattce_meam[a][b] == C11) { + latta = meam_data.lattce_meam[a][a]; + if (latta == DIA) { + phiaa = phi_meam(r, a, a); + phi_m = (3 * Eu - F2 - 2 * F1 - 5 * phiaa) / Z12; + } else { + phibb = phi_meam(r, b, b); + phi_m = (3 * Eu - F1 - 2 * F2 - 5 * phibb) / Z12; + } + } else if (meam_data.lattce_meam[a][b] == L12) { + phiaa = phi_meam(r, a, a); + // account for second neighbor a-a potential here... + get_Zij(&Z1nn, meam_data.lattce_meam[a][a]); + get_Zij2(&Z2nn, &arat, &scrn, meam_data.lattce_meam[a][a], + meam_data.Cmin_meam[a][a][a], meam_data.Cmax_meam[a][a][a]); + nmax = 10; + if (scrn > 0.0) { + for (n = 1; n <= nmax; n++) { + phiaa = + phiaa + + pow((-Z2nn * scrn / Z1nn), n) * phi_meam(r * pow(arat, n), a, a); + } + } + phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa; + + } else { + // + // potential is computed from Rose function and embedding energy + phi_m = (2 * Eu - F1 - F2) / Z12; + // + } + + // if r = 0, just return 0 + if (iszero(r)) { + phi_m = 0.0; + } + + return phi_m; +} //----------------------------------------------------------------------c // Compute background density for reference structure of each element - void compute_reference_density(void) - { - int a,Z,Z2,errorflag; - double gam,Gbar,shp[3+1]; - double rho0,rho0_2nn,arat,scrn; +void +compute_reference_density(void) +{ + int a, Z, Z2, errorflag; + double gam, Gbar, shp[3 + 1]; + double rho0, rho0_2nn, arat, scrn; -// loop over element types - for(a=1; a<=meam_data.neltypes; a++) { - Z = (int)meam_data.Z_meam[a]; - if (meam_data.ibar_meam[a]<=0) - Gbar = 1.0; - else { - get_shpfcn(shp,meam_data.lattce_meam[a][a]); - gam = (meam_data.t1_meam[a]*shp[1]+meam_data.t2_meam[a]*shp[2]+meam_data.t3_meam[a]*shp[3])/(Z*Z); - G_gam(gam,meam_data.ibar_meam[a],meam_data.gsmooth_factor,&Gbar,&errorflag); - } + // loop over element types + for (a = 1; a <= meam_data.neltypes; a++) { + Z = (int)meam_data.Z_meam[a]; + if (meam_data.ibar_meam[a] <= 0) + Gbar = 1.0; + else { + get_shpfcn(shp, meam_data.lattce_meam[a][a]); + gam = (meam_data.t1_meam[a] * shp[1] + meam_data.t2_meam[a] * shp[2] + + meam_data.t3_meam[a] * shp[3]) / + (Z * Z); + G_gam(gam, meam_data.ibar_meam[a], meam_data.gsmooth_factor, &Gbar, + &errorflag); + } -// The zeroth order density in the reference structure, with -// equilibrium spacing, is just the number of first neighbors times -// the rho0_meam coefficient... - rho0 = meam_data.rho0_meam[a]*Z; + // The zeroth order density in the reference structure, with + // equilibrium spacing, is just the number of first neighbors times + // the rho0_meam coefficient... + rho0 = meam_data.rho0_meam[a] * Z; -// ...unless we have unscreened second neighbors, in which case we -// add on the contribution from those (accounting for partial -// screening) - if (meam_data.nn2_meam[a][a]==1) { - get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[a][a],meam_data.Cmin_meam[a][a][a],meam_data.Cmax_meam[a][a][a]); - rho0_2nn = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta0_meam[a]*(arat-1)); - rho0 = rho0 + Z2*rho0_2nn*scrn; - } + // ...unless we have unscreened second neighbors, in which case we + // add on the contribution from those (accounting for partial + // screening) + if (meam_data.nn2_meam[a][a] == 1) { + get_Zij2(&Z2, &arat, &scrn, meam_data.lattce_meam[a][a], + meam_data.Cmin_meam[a][a][a], meam_data.Cmax_meam[a][a][a]); + rho0_2nn = + meam_data.rho0_meam[a] * fm_exp(-meam_data.beta0_meam[a] * (arat - 1)); + rho0 = rho0 + Z2 * rho0_2nn * scrn; + } - meam_data.rho_ref_meam[a] = rho0*Gbar; - } - } + meam_data.rho_ref_meam[a] = rho0 * Gbar; + } +} //----------------------------------------------------------------------c // Shape factors for various configurations - void get_shpfcn(double *s /* s(3) */, lattice_t latt) - { - if (latt==FCC || latt==BCC || latt==B1 || latt==B2) { - s[1] = 0.0; - s[2] = 0.0; - s[3] = 0.0; - } else if (latt==HCP) { - s[1] = 0.0; - s[2] = 0.0; - s[3] = 1.0/3.0; - } else if (latt==DIA) { - s[1] = 0.0; - s[2] = 0.0; - s[3] = 32.0/9.0; - } else if (latt==DIM) { - s[1] = 1.0; - s[2] = 2.0/3.0; -// s(3) = 1.d0 - s[3] = 0.40; - } else { - s[1] = 0.0; -// call error('Lattice not defined in get_shpfcn.') - } - } - +void +get_shpfcn(double* s /* s(3) */, lattice_t latt) +{ + if (latt == FCC || latt == BCC || latt == B1 || latt == B2) { + s[1] = 0.0; + s[2] = 0.0; + s[3] = 0.0; + } else if (latt == HCP) { + s[1] = 0.0; + s[2] = 0.0; + s[3] = 1.0 / 3.0; + } else if (latt == DIA) { + s[1] = 0.0; + s[2] = 0.0; + s[3] = 32.0 / 9.0; + } else if (latt == DIM) { + s[1] = 1.0; + s[2] = 2.0 / 3.0; + // s(3) = 1.d0 + s[3] = 0.40; + } else { + s[1] = 0.0; + // call error('Lattice not defined in get_shpfcn.') + } +} + //------------------------------------------------------------------------------c // Average weighting factors for the reference structure - void get_tavref(double *t11av,double *t21av,double *t31av,double *t12av,double *t22av,double *t32av, double t11,double t21,double t31,double t12,double t22,double t32, double r,int a,int b,lattice_t latt) - { - double rhoa01,rhoa02,a1,a2,rho01/*,rho02*/; +void +get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, + double* t22av, double* t32av, double t11, double t21, double t31, + double t12, double t22, double t32, double r, int a, int b, + lattice_t latt) +{ + double rhoa01, rhoa02, a1, a2, rho01 /*,rho02*/; -// For ialloy = 2, no averaging is done - if (meam_data.ialloy==2) { - *t11av = t11; - *t21av = t21; - *t31av = t31; - *t12av = t12; - *t22av = t22; - *t32av = t32; + // For ialloy = 2, no averaging is done + if (meam_data.ialloy == 2) { + *t11av = t11; + *t21av = t21; + *t31av = t31; + *t12av = t12; + *t22av = t22; + *t32av = t32; + } else { + if (latt == FCC || latt == BCC || latt == DIA || latt == HCP || + latt == B1 || latt == DIM || latt == B2) { + // all neighbors are of the opposite type + *t11av = t12; + *t21av = t22; + *t31av = t32; + *t12av = t11; + *t22av = t21; + *t32av = t31; + } else { + a1 = r / meam_data.re_meam[a][a] - 1.0; + a2 = r / meam_data.re_meam[b][b] - 1.0; + rhoa01 = meam_data.rho0_meam[a] * fm_exp(-meam_data.beta0_meam[a] * a1); + rhoa02 = meam_data.rho0_meam[b] * fm_exp(-meam_data.beta0_meam[b] * a2); + if (latt == L12) { + rho01 = 8 * rhoa01 + 4 * rhoa02; + *t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01; + *t12av = t11; + *t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01; + *t22av = t21; + *t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01; + *t32av = t31; } else { - if (latt==FCC || latt==BCC || latt==DIA || latt==HCP || latt==B1 || latt==DIM || latt==B2) { -// all neighbors are of the opposite type - *t11av = t12; - *t21av = t22; - *t31av = t32; - *t12av = t11; - *t22av = t21; - *t32av = t31; - } else { - a1 = r/meam_data.re_meam[a][a] - 1.0; - a2 = r/meam_data.re_meam[b][b] - 1.0; - rhoa01 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta0_meam[a]*a1); - rhoa02 = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta0_meam[b]*a2); - if (latt==L12) { - rho01 = 8*rhoa01 + 4*rhoa02; - *t11av = (8*t11*rhoa01 + 4*t12*rhoa02)/rho01; - *t12av = t11; - *t21av = (8*t21*rhoa01 + 4*t22*rhoa02)/rho01; - *t22av = t21; - *t31av = (8*t31*rhoa01 + 4*t32*rhoa02)/rho01; - *t32av = t31; - } else { -// call error('Lattice not defined in get_tavref.') - } - } + // call error('Lattice not defined in get_tavref.') } - } - + } + } +} + //------------------------------------------------------------------------------c // Number of neighbors for the reference structure - void get_Zij(int *Zij, lattice_t latt) - { - if (latt==FCC) - *Zij = 12; - else if (latt==BCC) - *Zij = 8; - else if (latt==HCP) - *Zij = 12; - else if (latt==B1) - *Zij = 6; - else if (latt==DIA) - *Zij = 4; - else if (latt==DIM) - *Zij = 1; - else if (latt==C11) - *Zij = 10; - else if (latt==L12) - *Zij = 12; - else if (latt==B2) - *Zij = 8; - else { -// call error('Lattice not defined in get_Zij.') - } - } - +void +get_Zij(int* Zij, lattice_t latt) +{ + if (latt == FCC) + *Zij = 12; + else if (latt == BCC) + *Zij = 8; + else if (latt == HCP) + *Zij = 12; + else if (latt == B1) + *Zij = 6; + else if (latt == DIA) + *Zij = 4; + else if (latt == DIM) + *Zij = 1; + else if (latt == C11) + *Zij = 10; + else if (latt == L12) + *Zij = 12; + else if (latt == B2) + *Zij = 8; + else { + // call error('Lattice not defined in get_Zij.') + } +} + //------------------------------------------------------------------------------c // Zij2 = number of second neighbors, a = distance ratio R1/R2, and S = second // neighbor screening function for lattice type "latt" - void get_Zij2(int *Zij2, double *a, double*S, lattice_t latt,double cmin,double cmax) - { - - double /*rratio,*/C,x,sijk; - int numscr = 0; +void +get_Zij2(int* Zij2, double* a, double* S, lattice_t latt, double cmin, + double cmax) +{ - if (latt==BCC) { - *Zij2 = 6; - *a = 2.0/sqrt(3.0); - numscr = 4; - } else if (latt==FCC) { - *Zij2 = 6; - *a = sqrt(2.0); - numscr = 4; - } else if (latt==DIA) { - *Zij2 = 0; - *a = sqrt(8.0/3.0); - numscr = 4; - if (cmin < 0.500001) { -// call error('can not do 2NN MEAM for dia') - } - } else if (latt==HCP) { - *Zij2 = 6; - *a = sqrt(2.0); - numscr = 4; - } else if (latt==B1) { - *Zij2 = 12; - *a = sqrt(2.0); - numscr = 2; - } else if (latt==L12) { - *Zij2 = 6; - *a = sqrt(2.0); - numscr = 4; - } else if (latt==B2) { - *Zij2 = 6; - *a = 2.0/sqrt(3.0); - numscr = 4; - } else if (latt==DIM) { -// this really shouldn't be allowed; make sure screening is zero - *Zij2 = 0; - *a = 1; - *S = 0; - return; - } else { -// call error('Lattice not defined in get_Zij2.') - } + double /*rratio,*/ C, x, sijk; + int numscr = 0; -// Compute screening for each first neighbor - C = 4.0/(*a * *a) - 1.0; - x = (C-cmin)/(cmax-cmin); - fcut(x,&sijk); -// There are numscr first neighbors screening the second neighbors - *S = pow(sijk,numscr); - - } - + if (latt == BCC) { + *Zij2 = 6; + *a = 2.0 / sqrt(3.0); + numscr = 4; + } else if (latt == FCC) { + *Zij2 = 6; + *a = sqrt(2.0); + numscr = 4; + } else if (latt == DIA) { + *Zij2 = 0; + *a = sqrt(8.0 / 3.0); + numscr = 4; + if (cmin < 0.500001) { + // call error('can not do 2NN MEAM for dia') + } + } else if (latt == HCP) { + *Zij2 = 6; + *a = sqrt(2.0); + numscr = 4; + } else if (latt == B1) { + *Zij2 = 12; + *a = sqrt(2.0); + numscr = 2; + } else if (latt == L12) { + *Zij2 = 6; + *a = sqrt(2.0); + numscr = 4; + } else if (latt == B2) { + *Zij2 = 6; + *a = 2.0 / sqrt(3.0); + numscr = 4; + } else if (latt == DIM) { + // this really shouldn't be allowed; make sure screening is zero + *Zij2 = 0; + *a = 1; + *S = 0; + return; + } else { + // call error('Lattice not defined in get_Zij2.') + } + // Compute screening for each first neighbor + C = 4.0 / (*a * *a) - 1.0; + x = (C - cmin) / (cmax - cmin); + fcut(x, &sijk); + // There are numscr first neighbors screening the second neighbors + *S = pow(sijk, numscr); +} //------------------------------------------------------------------------------c - void get_sijk(double C,int i,int j,int k, double *sijk) - { - double x; - x = (C-meam_data.Cmin_meam[i][j][k])/(meam_data.Cmax_meam[i][j][k]-meam_data.Cmin_meam[i][j][k]); - fcut(x,sijk); - } +void +get_sijk(double C, int i, int j, int k, double* sijk) +{ + double x; + x = (C - meam_data.Cmin_meam[i][j][k]) / + (meam_data.Cmax_meam[i][j][k] - meam_data.Cmin_meam[i][j][k]); + fcut(x, sijk); +} //------------------------------------------------------------------------------c // Calculate density functions, assuming reference configuration - void get_densref(double r,int a,int b,double *rho01,double *rho11,double *rho21,double *rho31, double *rho02,double *rho12,double *rho22,double *rho32) - { - double a1,a2; - double s[3+1]; - lattice_t lat; - int Zij1nn,Zij2nn; - double rhoa01nn,rhoa02nn; - double rhoa01,rhoa11,rhoa21,rhoa31; - double rhoa02,rhoa12,rhoa22,rhoa32; - double arat,scrn,denom; - double C,s111,s112,s221,S11,S22; +void +get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, + double* rho31, double* rho02, double* rho12, double* rho22, + double* rho32) +{ + double a1, a2; + double s[3 + 1]; + lattice_t lat; + int Zij1nn, Zij2nn; + double rhoa01nn, rhoa02nn; + double rhoa01, rhoa11, rhoa21, rhoa31; + double rhoa02, rhoa12, rhoa22, rhoa32; + double arat, scrn, denom; + double C, s111, s112, s221, S11, S22; - a1 = r/meam_data.re_meam[a][a] - 1.0; - a2 = r/meam_data.re_meam[b][b] - 1.0; + a1 = r / meam_data.re_meam[a][a] - 1.0; + a2 = r / meam_data.re_meam[b][b] - 1.0; - rhoa01 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta0_meam[a]*a1); - rhoa11 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta1_meam[a]*a1); - rhoa21 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta2_meam[a]*a1); - rhoa31 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta3_meam[a]*a1); - rhoa02 = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta0_meam[b]*a2); - rhoa12 = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta1_meam[b]*a2); - rhoa22 = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta2_meam[b]*a2); - rhoa32 = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta3_meam[b]*a2); + rhoa01 = meam_data.rho0_meam[a] * fm_exp(-meam_data.beta0_meam[a] * a1); + rhoa11 = meam_data.rho0_meam[a] * fm_exp(-meam_data.beta1_meam[a] * a1); + rhoa21 = meam_data.rho0_meam[a] * fm_exp(-meam_data.beta2_meam[a] * a1); + rhoa31 = meam_data.rho0_meam[a] * fm_exp(-meam_data.beta3_meam[a] * a1); + rhoa02 = meam_data.rho0_meam[b] * fm_exp(-meam_data.beta0_meam[b] * a2); + rhoa12 = meam_data.rho0_meam[b] * fm_exp(-meam_data.beta1_meam[b] * a2); + rhoa22 = meam_data.rho0_meam[b] * fm_exp(-meam_data.beta2_meam[b] * a2); + rhoa32 = meam_data.rho0_meam[b] * fm_exp(-meam_data.beta3_meam[b] * a2); - lat = meam_data.lattce_meam[a][b]; + lat = meam_data.lattce_meam[a][b]; - *rho11 = 0.0; - *rho21 = 0.0; - *rho31 = 0.0; - *rho12 = 0.0; - *rho22 = 0.0; - *rho32 = 0.0; + *rho11 = 0.0; + *rho21 = 0.0; + *rho31 = 0.0; + *rho12 = 0.0; + *rho22 = 0.0; + *rho32 = 0.0; - get_Zij(&Zij1nn,lat); + get_Zij(&Zij1nn, lat); - if (lat==FCC) { - *rho01 = 12.0*rhoa02; - *rho02 = 12.0*rhoa01; - } else if (lat==BCC) { - *rho01 = 8.0*rhoa02; - *rho02 = 8.0*rhoa01; - } else if (lat==B1) { - *rho01 = 6.0*rhoa02; - *rho02 = 6.0*rhoa01; - } else if (lat==DIA) { - *rho01 = 4.0*rhoa02; - *rho02 = 4.0*rhoa01; - *rho31 = 32.0/9.0*rhoa32*rhoa32; - *rho32 = 32.0/9.0*rhoa31*rhoa31; - } else if (lat==HCP) { - *rho01 = 12*rhoa02; - *rho02 = 12*rhoa01; - *rho31 = 1.0/3.0*rhoa32*rhoa32; - *rho32 = 1.0/3.0*rhoa31*rhoa31; - } else if (lat==DIM) { - get_shpfcn(s,DIM); - *rho01 = rhoa02; - *rho02 = rhoa01; - *rho11 = s[1]*rhoa12*rhoa12; - *rho12 = s[1]*rhoa11*rhoa11; - *rho21 = s[2]*rhoa22*rhoa22; - *rho22 = s[2]*rhoa21*rhoa21; - *rho31 = s[3]*rhoa32*rhoa32; - *rho32 = s[3]*rhoa31*rhoa31; - } else if (lat==C11) { - *rho01 = rhoa01; - *rho02 = rhoa02; - *rho11 = rhoa11; - *rho12 = rhoa12; - *rho21 = rhoa21; - *rho22 = rhoa22; - *rho31 = rhoa31; - *rho32 = rhoa32; - } else if (lat==L12) { - *rho01 = 8*rhoa01 + 4*rhoa02; - *rho02 = 12*rhoa01; - if (meam_data.ialloy==1) { - *rho21 = 8./3.*pow(rhoa21*meam_data.t2_meam[a]-rhoa22*meam_data.t2_meam[b],2); - denom = 8*rhoa01*pow(meam_data.t2_meam[a],2) + 4*rhoa02*pow(meam_data.t2_meam[b],2); - if (denom > 0.) - *rho21 = *rho21/denom * *rho01; - } else - *rho21 = 8./3.*(rhoa21-rhoa22)*(rhoa21-rhoa22); - } else if (lat==B2) { - *rho01 = 8.0*rhoa02; - *rho02 = 8.0*rhoa01; - } else { -// call error('Lattice not defined in get_densref.') - } + if (lat == FCC) { + *rho01 = 12.0 * rhoa02; + *rho02 = 12.0 * rhoa01; + } else if (lat == BCC) { + *rho01 = 8.0 * rhoa02; + *rho02 = 8.0 * rhoa01; + } else if (lat == B1) { + *rho01 = 6.0 * rhoa02; + *rho02 = 6.0 * rhoa01; + } else if (lat == DIA) { + *rho01 = 4.0 * rhoa02; + *rho02 = 4.0 * rhoa01; + *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32; + *rho32 = 32.0 / 9.0 * rhoa31 * rhoa31; + } else if (lat == HCP) { + *rho01 = 12 * rhoa02; + *rho02 = 12 * rhoa01; + *rho31 = 1.0 / 3.0 * rhoa32 * rhoa32; + *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31; + } else if (lat == DIM) { + get_shpfcn(s, DIM); + *rho01 = rhoa02; + *rho02 = rhoa01; + *rho11 = s[1] * rhoa12 * rhoa12; + *rho12 = s[1] * rhoa11 * rhoa11; + *rho21 = s[2] * rhoa22 * rhoa22; + *rho22 = s[2] * rhoa21 * rhoa21; + *rho31 = s[3] * rhoa32 * rhoa32; + *rho32 = s[3] * rhoa31 * rhoa31; + } else if (lat == C11) { + *rho01 = rhoa01; + *rho02 = rhoa02; + *rho11 = rhoa11; + *rho12 = rhoa12; + *rho21 = rhoa21; + *rho22 = rhoa22; + *rho31 = rhoa31; + *rho32 = rhoa32; + } else if (lat == L12) { + *rho01 = 8 * rhoa01 + 4 * rhoa02; + *rho02 = 12 * rhoa01; + if (meam_data.ialloy == 1) { + *rho21 = + 8. / 3. * + pow(rhoa21 * meam_data.t2_meam[a] - rhoa22 * meam_data.t2_meam[b], 2); + denom = 8 * rhoa01 * pow(meam_data.t2_meam[a], 2) + + 4 * rhoa02 * pow(meam_data.t2_meam[b], 2); + if (denom > 0.) + *rho21 = *rho21 / denom * *rho01; + } else + *rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22); + } else if (lat == B2) { + *rho01 = 8.0 * rhoa02; + *rho02 = 8.0 * rhoa01; + } else { + // call error('Lattice not defined in get_densref.') + } - if (meam_data.nn2_meam[a][b]==1) { + if (meam_data.nn2_meam[a][b] == 1) { - get_Zij2(&Zij2nn,&arat,&scrn,lat,meam_data.Cmin_meam[a][a][b],meam_data.Cmax_meam[a][a][b]); + get_Zij2(&Zij2nn, &arat, &scrn, lat, meam_data.Cmin_meam[a][a][b], + meam_data.Cmax_meam[a][a][b]); - a1 = arat*r/meam_data.re_meam[a][a] - 1.0; - a2 = arat*r/meam_data.re_meam[b][b] - 1.0; + a1 = arat * r / meam_data.re_meam[a][a] - 1.0; + a2 = arat * r / meam_data.re_meam[b][b] - 1.0; - rhoa01nn = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta0_meam[a]*a1); - rhoa02nn = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta0_meam[b]*a2); + rhoa01nn = meam_data.rho0_meam[a] * fm_exp(-meam_data.beta0_meam[a] * a1); + rhoa02nn = meam_data.rho0_meam[b] * fm_exp(-meam_data.beta0_meam[b] * a2); - if (lat==L12) { -// As usual, L12 thinks it's special; we need to be careful computing -// the screening functions - C = 1.0; - get_sijk(C,a,a,a,&s111); - get_sijk(C,a,a,b,&s112); - get_sijk(C,b,b,a,&s221); - S11 = s111 * s111 * s112 * s112; - S22 = pow(s221,4); - *rho01 = *rho01 + 6*S11*rhoa01nn; - *rho02 = *rho02 + 6*S22*rhoa02nn; + if (lat == L12) { + // As usual, L12 thinks it's special; we need to be careful computing + // the screening functions + C = 1.0; + get_sijk(C, a, a, a, &s111); + get_sijk(C, a, a, b, &s112); + get_sijk(C, b, b, a, &s221); + S11 = s111 * s111 * s112 * s112; + S22 = pow(s221, 4); + *rho01 = *rho01 + 6 * S11 * rhoa01nn; + *rho02 = *rho02 + 6 * S22 * rhoa02nn; - } else { -// For other cases, assume that second neighbor is of same type, -// first neighbor may be of different type + } else { + // For other cases, assume that second neighbor is of same type, + // first neighbor may be of different type - *rho01 = *rho01 + Zij2nn*scrn*rhoa01nn; + *rho01 = *rho01 + Zij2nn * scrn * rhoa01nn; -// Assume Zij2nn and arat don't depend on order, but scrn might - get_Zij2(&Zij2nn,&arat,&scrn,lat,meam_data.Cmin_meam[b][b][a],meam_data.Cmax_meam[b][b][a]); - *rho02 = *rho02 + Zij2nn*scrn*rhoa02nn; - } - } - } + // Assume Zij2nn and arat don't depend on order, but scrn might + get_Zij2(&Zij2nn, &arat, &scrn, lat, meam_data.Cmin_meam[b][b][a], + meam_data.Cmax_meam[b][b][a]); + *rho02 = *rho02 + Zij2nn * scrn * rhoa02nn; + } + } +} //--------------------------------------------------------------------- // Compute ZBL potential // - double zbl(double r, int z1, int z2) - { - int i; - const double c[]={0.028171,0.28022,0.50986,0.18175}; - const double d[]={0.20162,0.40290,0.94229,3.1998}; - const double azero=0.4685; - const double cc=14.3997; - double a,x; -// azero = (9pi^2/128)^1/3 (0.529) Angstroms - a = azero/(pow(z1,0.23)+pow(z2,0.23)); - double result = 0.0; - x = r/a; - for(i=0; i<=3; i++) { - result = result + c[i]*fm_exp(-d[i]*x); - } - if (r > 0.0) result = result*z1*z2/r*cc; - return result; - } +double +zbl(double r, int z1, int z2) +{ + int i; + const double c[] = { 0.028171, 0.28022, 0.50986, 0.18175 }; + const double d[] = { 0.20162, 0.40290, 0.94229, 3.1998 }; + const double azero = 0.4685; + const double cc = 14.3997; + double a, x; + // azero = (9pi^2/128)^1/3 (0.529) Angstroms + a = azero / (pow(z1, 0.23) + pow(z2, 0.23)); + double result = 0.0; + x = r / a; + for (i = 0; i <= 3; i++) { + result = result + c[i] * fm_exp(-d[i] * x); + } + if (r > 0.0) + result = result * z1 * z2 / r * cc; + return result; +} //--------------------------------------------------------------------- // Compute Rose energy function, I.16 // - double erose(double r,double re,double alpha,double Ec,double repuls,double attrac,int form) - { - double astar,a3; - double result = 0.0; +double +erose(double r, double re, double alpha, double Ec, double repuls, + double attrac, int form) +{ + double astar, a3; + double result = 0.0; - if (r > 0.0) { - astar = alpha * (r/re - 1.0); - a3 = 0.0; - if (astar>=0) - a3 = attrac; - else if (astar < 0) - a3 = repuls; + if (r > 0.0) { + astar = alpha * (r / re - 1.0); + a3 = 0.0; + if (astar >= 0) + a3 = attrac; + else if (astar < 0) + a3 = repuls; - if (form==1) - result = -Ec*(1+astar+(-attrac+repuls/r)* pow(astar,3))*fm_exp(-astar); - else if (form==2) - result = -Ec * (1 +astar + a3*pow(astar,3))*fm_exp(-astar); - else - result = -Ec * (1+ astar + a3*pow(astar,3)/(r/re))*fm_exp(-astar); - - } - return result; - } + if (form == 1) + result = -Ec * (1 + astar + (-attrac + repuls / r) * pow(astar, 3)) * + fm_exp(-astar); + else if (form == 2) + result = -Ec * (1 + astar + a3 * pow(astar, 3)) * fm_exp(-astar); + else + result = + -Ec * (1 + astar + a3 * pow(astar, 3) / (r / re)) * fm_exp(-astar); + } + return result; +} // ----------------------------------------------------------------------- - void interpolate_meam(int ind) - { - int j; - double drar; +void +interpolate_meam(int ind) +{ + int j; + double drar; -// map to coefficient space + // map to coefficient space - meam_data.nrar = meam_data.nr; - drar = meam_data.dr; - meam_data.rdrar = 1.0/drar; + meam_data.nrar = meam_data.nr; + drar = meam_data.dr; + meam_data.rdrar = 1.0 / drar; -// phir interp - for(j=1; j<=meam_data.nrar; j++) { - arr2(meam_data.phirar,j,ind) = arr2(meam_data.phir,j,ind); - } - arr2(meam_data.phirar1,1,ind) = arr2(meam_data.phirar,2,ind)-arr2(meam_data.phirar,1,ind); - arr2(meam_data.phirar1,2,ind) = 0.5*(arr2(meam_data.phirar,3,ind)-arr2(meam_data.phirar,1,ind)); - arr2(meam_data.phirar1,meam_data.nrar-1,ind) = 0.5*(arr2(meam_data.phirar,meam_data.nrar,ind) -arr2(meam_data.phirar,meam_data.nrar-2,ind)); - arr2(meam_data.phirar1,meam_data.nrar,ind) = 0.0; - for(j=3; j<=meam_data.nrar-2; j++) { - arr2(meam_data.phirar1,j,ind) = ((arr2(meam_data.phirar,j-2,ind)-arr2(meam_data.phirar,j+2,ind)) + 8.0*(arr2(meam_data.phirar,j+1,ind)-arr2(meam_data.phirar,j-1,ind)))/12.; - } + // phir interp + for (j = 1; j <= meam_data.nrar; j++) { + arr2(meam_data.phirar, j, ind) = arr2(meam_data.phir, j, ind); + } + arr2(meam_data.phirar1, 1, ind) = + arr2(meam_data.phirar, 2, ind) - arr2(meam_data.phirar, 1, ind); + arr2(meam_data.phirar1, 2, ind) = + 0.5 * (arr2(meam_data.phirar, 3, ind) - arr2(meam_data.phirar, 1, ind)); + arr2(meam_data.phirar1, meam_data.nrar - 1, ind) = + 0.5 * (arr2(meam_data.phirar, meam_data.nrar, ind) - + arr2(meam_data.phirar, meam_data.nrar - 2, ind)); + arr2(meam_data.phirar1, meam_data.nrar, ind) = 0.0; + for (j = 3; j <= meam_data.nrar - 2; j++) { + arr2(meam_data.phirar1, j, ind) = + ((arr2(meam_data.phirar, j - 2, ind) - + arr2(meam_data.phirar, j + 2, ind)) + + 8.0 * (arr2(meam_data.phirar, j + 1, ind) - + arr2(meam_data.phirar, j - 1, ind))) / + 12.; + } - for(j=1; j<=meam_data.nrar-1; j++) { - arr2(meam_data.phirar2,j,ind) = 3.0*(arr2(meam_data.phirar,j+1,ind)-arr2(meam_data.phirar,j,ind)) - 2.0*arr2(meam_data.phirar1,j,ind) - arr2(meam_data.phirar1,j+1,ind); - arr2(meam_data.phirar3,j,ind) = arr2(meam_data.phirar1,j,ind) + arr2(meam_data.phirar1,j+1,ind) - 2.0*(arr2(meam_data.phirar,j+1,ind)-arr2(meam_data.phirar,j,ind)); - } - arr2(meam_data.phirar2,meam_data.nrar,ind) = 0.0; - arr2(meam_data.phirar3,meam_data.nrar,ind) = 0.0; + for (j = 1; j <= meam_data.nrar - 1; j++) { + arr2(meam_data.phirar2, j, ind) = + 3.0 * + (arr2(meam_data.phirar, j + 1, ind) - arr2(meam_data.phirar, j, ind)) - + 2.0 * arr2(meam_data.phirar1, j, ind) - + arr2(meam_data.phirar1, j + 1, ind); + arr2(meam_data.phirar3, j, ind) = + arr2(meam_data.phirar1, j, ind) + arr2(meam_data.phirar1, j + 1, ind) - + 2.0 * + (arr2(meam_data.phirar, j + 1, ind) - arr2(meam_data.phirar, j, ind)); + } + arr2(meam_data.phirar2, meam_data.nrar, ind) = 0.0; + arr2(meam_data.phirar3, meam_data.nrar, ind) = 0.0; - for(j=1; j<=meam_data.nrar; j++) { - arr2(meam_data.phirar4,j,ind) = arr2(meam_data.phirar1,j,ind)/drar; - arr2(meam_data.phirar5,j,ind) = 2.0*arr2(meam_data.phirar2,j,ind)/drar; - arr2(meam_data.phirar6,j,ind) = 3.0*arr2(meam_data.phirar3,j,ind)/drar; - } - - } + for (j = 1; j <= meam_data.nrar; j++) { + arr2(meam_data.phirar4, j, ind) = arr2(meam_data.phirar1, j, ind) / drar; + arr2(meam_data.phirar5, j, ind) = + 2.0 * arr2(meam_data.phirar2, j, ind) / drar; + arr2(meam_data.phirar6, j, ind) = + 3.0 * arr2(meam_data.phirar3, j, ind) / drar; + } +} //--------------------------------------------------------------------- // Compute Rose energy function, I.16 // - double compute_phi(double rij, int elti, int eltj) - { - double pp; - int ind, kk; - - ind = meam_data.eltind[elti][eltj]; - pp = rij*meam_data.rdrar + 1.0; - kk = (int)pp; - kk = min(kk,meam_data.nrar-1); - pp = pp - kk; - pp = min(pp,1.0); - double result = ((arr2(meam_data.phirar3,kk,ind)*pp + arr2(meam_data.phirar2,kk,ind))*pp + arr2(meam_data.phirar1,kk,ind))*pp + arr2(meam_data.phirar,kk,ind); - - return result; - } +double +compute_phi(double rij, int elti, int eltj) +{ + double pp; + int ind, kk; + ind = meam_data.eltind[elti][eltj]; + pp = rij * meam_data.rdrar + 1.0; + kk = (int)pp; + kk = min(kk, meam_data.nrar - 1); + pp = pp - kk; + pp = min(pp, 1.0); + double result = ((arr2(meam_data.phirar3, kk, ind) * pp + + arr2(meam_data.phirar2, kk, ind)) * + pp + + arr2(meam_data.phirar1, kk, ind)) * + pp + + arr2(meam_data.phirar, kk, ind); + return result; +} } diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp index 5f6dfcb317..eac03277f2 100644 --- a/src/USER-MEAMC/meam_setup_global.cpp +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -17,85 +17,83 @@ extern "C" { // // - void meam_setup_global_(int *nelt, int *lat, double *z, int *ielement, double *atwt, double *alpha, - double *b0, double *b1, double *b2, double *b3, double *alat, double *esub, double *asub, - double *t0, double *t1, double *t2, double *t3, double *rozero, int *ibar) - { - - int i; - double tmplat[maxelt]; +void +meam_setup_global_(int* nelt, int* lat, double* z, int* ielement, double* atwt, + double* alpha, double* b0, double* b1, double* b2, + double* b3, double* alat, double* esub, double* asub, + double* t0, double* t1, double* t2, double* t3, + double* rozero, int* ibar) +{ - meam_data.neltypes = *nelt; + int i; + double tmplat[maxelt]; - for(i=1; i<=*nelt; i++) { - if (arr1v(lat,i)==0) - meam_data.lattce_meam[i][i] = FCC; - else if (arr1v(lat,i)==1) - meam_data.lattce_meam[i][i] = BCC; - else if (arr1v(lat,i)==2) - meam_data.lattce_meam[i][i] = HCP; - else if (arr1v(lat,i)==3) - meam_data.lattce_meam[i][i] = DIM; - else if (arr1v(lat,i)==4) - meam_data.lattce_meam[i][i] = DIA; - else { -// unknown - } + meam_data.neltypes = *nelt; - meam_data.Z_meam[i] = arr1v(z,i); - meam_data.ielt_meam[i] = arr1v(ielement,i); - meam_data.alpha_meam[i][i] = arr1v(alpha,i); - meam_data.beta0_meam[i] = arr1v(b0,i); - meam_data.beta1_meam[i] = arr1v(b1,i); - meam_data.beta2_meam[i] = arr1v(b2,i); - meam_data.beta3_meam[i] = arr1v(b3,i); - tmplat[i] = arr1v(alat,i); - meam_data.Ec_meam[i][i] = arr1v(esub,i); - meam_data.A_meam[i] = arr1v(asub,i); - meam_data.t0_meam[i] = arr1v(t0,i); - meam_data.t1_meam[i] = arr1v(t1,i); - meam_data.t2_meam[i] = arr1v(t2,i); - meam_data.t3_meam[i] = arr1v(t3,i); - meam_data.rho0_meam[i] = arr1v(rozero,i); - meam_data.ibar_meam[i] = arr1v(ibar,i); + for (i = 1; i <= *nelt; i++) { + if (arr1v(lat, i) == 0) + meam_data.lattce_meam[i][i] = FCC; + else if (arr1v(lat, i) == 1) + meam_data.lattce_meam[i][i] = BCC; + else if (arr1v(lat, i) == 2) + meam_data.lattce_meam[i][i] = HCP; + else if (arr1v(lat, i) == 3) + meam_data.lattce_meam[i][i] = DIM; + else if (arr1v(lat, i) == 4) + meam_data.lattce_meam[i][i] = DIA; + else { + // unknown + } - if (meam_data.lattce_meam[i][i]==FCC) - meam_data.re_meam[i][i] = tmplat[i]/sqrt(2.0); - else if (meam_data.lattce_meam[i][i]==BCC) - meam_data.re_meam[i][i] = tmplat[i]*sqrt(3.0)/2.0; - else if (meam_data.lattce_meam[i][i]==HCP) - meam_data.re_meam[i][i] = tmplat[i]; - else if (meam_data.lattce_meam[i][i]==DIM) - meam_data.re_meam[i][i] = tmplat[i]; - else if (meam_data.lattce_meam[i][i]==DIA) - meam_data.re_meam[i][i] = tmplat[i]*sqrt(3.0)/4.0; - else { -// error - } - - } - - -// Set some defaults - meam_data.rc_meam = 4.0; - meam_data.delr_meam = 0.1; - setall2d(meam_data.attrac_meam, 0.0); - setall2d(meam_data.repuls_meam, 0.0); - setall3d(meam_data.Cmax_meam, 2.8); - setall3d(meam_data.Cmin_meam, 2.0); - setall2d(meam_data.ebound_meam, pow(2.8,2)/(4.0*(2.8-1.0))); - setall2d(meam_data.delta_meam, 0.0); - setall2d(meam_data.nn2_meam, 0); - setall2d(meam_data.zbl_meam, 1); - meam_data.gsmooth_factor = 99.0; - meam_data.augt1 = 1; - meam_data.ialloy = 0; - meam_data.mix_ref_t = 0; - meam_data.emb_lin_neg = 0; - meam_data.bkgd_dyn = 0; - meam_data.erose_form = 0; - - } + meam_data.Z_meam[i] = arr1v(z, i); + meam_data.ielt_meam[i] = arr1v(ielement, i); + meam_data.alpha_meam[i][i] = arr1v(alpha, i); + meam_data.beta0_meam[i] = arr1v(b0, i); + meam_data.beta1_meam[i] = arr1v(b1, i); + meam_data.beta2_meam[i] = arr1v(b2, i); + meam_data.beta3_meam[i] = arr1v(b3, i); + tmplat[i] = arr1v(alat, i); + meam_data.Ec_meam[i][i] = arr1v(esub, i); + meam_data.A_meam[i] = arr1v(asub, i); + meam_data.t0_meam[i] = arr1v(t0, i); + meam_data.t1_meam[i] = arr1v(t1, i); + meam_data.t2_meam[i] = arr1v(t2, i); + meam_data.t3_meam[i] = arr1v(t3, i); + meam_data.rho0_meam[i] = arr1v(rozero, i); + meam_data.ibar_meam[i] = arr1v(ibar, i); + if (meam_data.lattce_meam[i][i] == FCC) + meam_data.re_meam[i][i] = tmplat[i] / sqrt(2.0); + else if (meam_data.lattce_meam[i][i] == BCC) + meam_data.re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0; + else if (meam_data.lattce_meam[i][i] == HCP) + meam_data.re_meam[i][i] = tmplat[i]; + else if (meam_data.lattce_meam[i][i] == DIM) + meam_data.re_meam[i][i] = tmplat[i]; + else if (meam_data.lattce_meam[i][i] == DIA) + meam_data.re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0; + else { + // error + } + } + // Set some defaults + meam_data.rc_meam = 4.0; + meam_data.delr_meam = 0.1; + setall2d(meam_data.attrac_meam, 0.0); + setall2d(meam_data.repuls_meam, 0.0); + setall3d(meam_data.Cmax_meam, 2.8); + setall3d(meam_data.Cmin_meam, 2.0); + setall2d(meam_data.ebound_meam, pow(2.8, 2) / (4.0 * (2.8 - 1.0))); + setall2d(meam_data.delta_meam, 0.0); + setall2d(meam_data.nn2_meam, 0); + setall2d(meam_data.zbl_meam, 1); + meam_data.gsmooth_factor = 99.0; + meam_data.augt1 = 1; + meam_data.ialloy = 0; + meam_data.mix_ref_t = 0; + meam_data.emb_lin_neg = 0; + meam_data.bkgd_dyn = 0; + meam_data.erose_form = 0; +} } \ No newline at end of file diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp index ff89eda96a..4be38a0cb1 100644 --- a/src/USER-MEAMC/meam_setup_param.cpp +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -3,22 +3,23 @@ extern "C" { // // do a sanity check on index parameters - void meam_checkindex(int num, int lim, int nidx, int *idx /*idx(3)*/, int *ierr) - { - //: idx[0..2] - *ierr = 0; - if (nidx < num) { - *ierr = 2; - return; - } +void +meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr) +{ + //: idx[0..2] + *ierr = 0; + if (nidx < num) { + *ierr = 2; + return; + } - for (int i=0; i lim)) { - *ierr = 3; - return; - } - } - } + for (int i = 0; i < num; i++) { + if ((idx[i] < 1) || (idx[i] > lim)) { + *ierr = 3; + return; + } + } +} // // Declaration in pair_meam.h: @@ -56,171 +57,184 @@ extern "C" { // 19 = emb_lin_neg // 20 = bkgd_dyn - void meam_setup_param_(int *which_p, double *value_p, int *nindex_p, int *index /*index(3)*/, int *errorflag) - { - //: index[0..2] - int i1, i2, val; - *errorflag = 0; - int which = *which_p; - double value = *value_p; - int nindex = *nindex_p; +void +meam_setup_param_(int* which_p, double* value_p, int* nindex_p, + int* index /*index(3)*/, int* errorflag) +{ + //: index[0..2] + int i1, i2, val; + *errorflag = 0; + int which = *which_p; + double value = *value_p; + int nindex = *nindex_p; - switch(which) { - // 0 = Ec_meam - case 0: - meam_checkindex(2,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - meam_data.Ec_meam[index[0]][index[1]] = value; - break; + switch (which) { + // 0 = Ec_meam + case 0: + meam_checkindex(2, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + meam_data.Ec_meam[index[0]][index[1]] = value; + break; - // 1 = alpha_meam - case 1: - meam_checkindex(2,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - meam_data.alpha_meam[index[0]][index[1]] = value; - break; + // 1 = alpha_meam + case 1: + meam_checkindex(2, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + meam_data.alpha_meam[index[0]][index[1]] = value; + break; - // 2 = rho0_meam - case 2: - meam_checkindex(1,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - meam_data.rho0_meam[index[0]] = value; - break; + // 2 = rho0_meam + case 2: + meam_checkindex(1, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + meam_data.rho0_meam[index[0]] = value; + break; - // 3 = delta_meam - case 3: - meam_checkindex(2,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - meam_data.delta_meam[index[0]][index[1]] = value; - break; + // 3 = delta_meam + case 3: + meam_checkindex(2, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + meam_data.delta_meam[index[0]][index[1]] = value; + break; - // 4 = lattce_meam - case 4: - meam_checkindex(2,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - val = (int)value; + // 4 = lattce_meam + case 4: + meam_checkindex(2, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + val = (int)value; - if (val==0) - meam_data.lattce_meam[index[0]][index[1]] = FCC; - else if (val==1) - meam_data.lattce_meam[index[0]][index[1]] = BCC; - else if (val==2) - meam_data.lattce_meam[index[0]][index[1]] = HCP; - else if (val==3) - meam_data.lattce_meam[index[0]][index[1]] = DIM; - else if (val==4) - meam_data.lattce_meam[index[0]][index[1]] = DIA; - else if (val==5) - meam_data.lattce_meam[index[0]][index[1]] = B1; - else if (val==6) - meam_data.lattce_meam[index[0]][index[1]] = C11; - else if (val==7) - meam_data.lattce_meam[index[0]][index[1]] = L12; - else if (val==8) - meam_data.lattce_meam[index[0]][index[1]] = B2; - break; + if (val == 0) + meam_data.lattce_meam[index[0]][index[1]] = FCC; + else if (val == 1) + meam_data.lattce_meam[index[0]][index[1]] = BCC; + else if (val == 2) + meam_data.lattce_meam[index[0]][index[1]] = HCP; + else if (val == 3) + meam_data.lattce_meam[index[0]][index[1]] = DIM; + else if (val == 4) + meam_data.lattce_meam[index[0]][index[1]] = DIA; + else if (val == 5) + meam_data.lattce_meam[index[0]][index[1]] = B1; + else if (val == 6) + meam_data.lattce_meam[index[0]][index[1]] = C11; + else if (val == 7) + meam_data.lattce_meam[index[0]][index[1]] = L12; + else if (val == 8) + meam_data.lattce_meam[index[0]][index[1]] = B2; + break; - // 5 = attrac_meam - case 5: - meam_checkindex(2,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - meam_data.attrac_meam[index[0]][index[1]] = value; - break; + // 5 = attrac_meam + case 5: + meam_checkindex(2, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + meam_data.attrac_meam[index[0]][index[1]] = value; + break; - // 6 = repuls_meam - case 6: - meam_checkindex(2,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - meam_data.repuls_meam[index[0]][index[1]] = value; - break; + // 6 = repuls_meam + case 6: + meam_checkindex(2, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + meam_data.repuls_meam[index[0]][index[1]] = value; + break; - // 7 = nn2_meam - case 7: - meam_checkindex(2,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - i1 = min(index[0],index[1]); - i2 = max(index[0],index[1]); - meam_data.nn2_meam[i1][i2] = (int)value; - break; + // 7 = nn2_meam + case 7: + meam_checkindex(2, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + i1 = min(index[0], index[1]); + i2 = max(index[0], index[1]); + meam_data.nn2_meam[i1][i2] = (int)value; + break; - // 8 = Cmin_meam - case 8: - meam_checkindex(3,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - meam_data.Cmin_meam[index[0]][index[1]][index[2]] = value; - break; + // 8 = Cmin_meam + case 8: + meam_checkindex(3, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + meam_data.Cmin_meam[index[0]][index[1]][index[2]] = value; + break; - // 9 = Cmax_meam - case 9: - meam_checkindex(3,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - meam_data.Cmax_meam[index[0]][index[1]][index[2]] = value; - break; + // 9 = Cmax_meam + case 9: + meam_checkindex(3, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + meam_data.Cmax_meam[index[0]][index[1]][index[2]] = value; + break; - // 10 = rc_meam - case 10: - meam_data.rc_meam = value; - break; + // 10 = rc_meam + case 10: + meam_data.rc_meam = value; + break; - // 11 = delr_meam - case 11: - meam_data.delr_meam = value; - break; + // 11 = delr_meam + case 11: + meam_data.delr_meam = value; + break; - // 12 = augt1 - case 12: - meam_data.augt1 = (int)value; - break; + // 12 = augt1 + case 12: + meam_data.augt1 = (int)value; + break; - // 13 = gsmooth - case 13: - meam_data.gsmooth_factor = value; - break; + // 13 = gsmooth + case 13: + meam_data.gsmooth_factor = value; + break; - // 14 = re_meam - case 14: - meam_checkindex(2,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - meam_data.re_meam[index[0]][index[1]] = value; - break; + // 14 = re_meam + case 14: + meam_checkindex(2, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + meam_data.re_meam[index[0]][index[1]] = value; + break; - // 15 = ialloy - case 15: - meam_data.ialloy = (int)value; - break; + // 15 = ialloy + case 15: + meam_data.ialloy = (int)value; + break; - // 16 = mixture_ref_t - case 16: - meam_data.mix_ref_t = (int)value; - break; + // 16 = mixture_ref_t + case 16: + meam_data.mix_ref_t = (int)value; + break; - // 17 = erose_form - case 17: - meam_data.erose_form = (int)value; - break; + // 17 = erose_form + case 17: + meam_data.erose_form = (int)value; + break; - // 18 = zbl_meam - case 18: - meam_checkindex(2,maxelt,nindex,index,errorflag); - if (*errorflag!=0) return; - i1 = min(index[0],index[1]); - i2 = max(index[0],index[1]); - meam_data.zbl_meam[i1][i2] = (int)value; - break; + // 18 = zbl_meam + case 18: + meam_checkindex(2, maxelt, nindex, index, errorflag); + if (*errorflag != 0) + return; + i1 = min(index[0], index[1]); + i2 = max(index[0], index[1]); + meam_data.zbl_meam[i1][i2] = (int)value; + break; - // 19 = emb_lin_neg - case 19: - meam_data.emb_lin_neg = (int)value; - break; + // 19 = emb_lin_neg + case 19: + meam_data.emb_lin_neg = (int)value; + break; - // 20 = bkgd_dyn - case 20: - meam_data.bkgd_dyn = (int)value; - break; - - default: - *errorflag = 1; - } - } + // 20 = bkgd_dyn + case 20: + meam_data.bkgd_dyn = (int)value; + break; + default: + *errorflag = 1; + } +} } \ No newline at end of file