diff --git a/src/USER-MEAMC/Install.sh b/src/USER-MEAMC/Install.sh deleted file mode 100644 index 2b4098a47..000000000 --- a/src/USER-MEAMC/Install.sh +++ /dev/null @@ -1,47 +0,0 @@ -# Install/unInstall package files in LAMMPS -# mode = 0/1/2 for uninstall/install/update - -# this is default Install.sh for all packages -# if package has an auxiliary library or a file with a dependency, -# then package dir has its own customized Install.sh - -mode=$1 - -# enforce using portable C locale -LC_ALL=C -export LC_ALL - -# arg1 = file, arg2 = file it depends on - -action () { - if (test $mode = 0) then - rm -f ../$1 - elif (! cmp -s $1 ../$1) then - if (test -z "$2" || test -e ../$2) then - cp $1 .. - if (test $mode = 2) then - echo " updating src/$1" - fi - fi - elif (test -n "$2") then - if (test ! -e ../$2) then - rm -f ../$1 - fi - fi -} - -# all package files with no dependencies - -for file in *.cpp *.h; do - test -f ${file} && action $file -done - - -# additional files - -for file in meam_*.c; do - test -f ${file} && action ${file} -done - -action fm_exp.c -action meam.h diff --git a/src/USER-MEAMC/fm_exp.c b/src/USER-MEAMC/fm_exp.cpp similarity index 99% rename from src/USER-MEAMC/fm_exp.c rename to src/USER-MEAMC/fm_exp.cpp index 00b957b27..253af868f 100644 --- a/src/USER-MEAMC/fm_exp.c +++ b/src/USER-MEAMC/fm_exp.cpp @@ -1,133 +1,135 @@ +extern "C" { /* Copyright (c) 2012,2013 Axel Kohlmeyer All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. * Neither the name of the nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* faster versions of 2**x, e**x, and 10**x in single and double precision. * * Based on the Cephes math library 2.8 */ #include #include /* internal definitions for the fastermath library */ /* IEEE 754 double precision floating point data manipulation */ typedef union { double f; uint64_t u; struct {int32_t i0,i1;}; } udi_t; #define FM_DOUBLE_BIAS 1023 #define FM_DOUBLE_EMASK 2146435072 #define FM_DOUBLE_MBITS 20 #define FM_DOUBLE_MMASK 1048575 #define FM_DOUBLE_EZERO 1072693248 /* generate 2**num in floating point by bitshifting */ #define FM_DOUBLE_INIT_EXP(var,num) \ var.i0 = 0; \ var.i1 = (((int) num) + FM_DOUBLE_BIAS) << 20 /* double precision constants */ #define FM_DOUBLE_LOG2OFE 1.4426950408889634074 #define FM_DOUBLE_LOGEOF2 6.9314718055994530942e-1 #define FM_DOUBLE_LOG2OF10 3.32192809488736234789 #define FM_DOUBLE_LOG10OF2 3.0102999566398119521e-1 #define FM_DOUBLE_LOG10OFE 4.3429448190325182765e-1 #define FM_DOUBLE_SQRT2 1.41421356237309504880 #define FM_DOUBLE_SQRTH 0.70710678118654752440 /* optimizer friendly implementation of exp2(x). * * strategy: * * split argument into an integer part and a fraction: * ipart = floor(x+0.5); * fpart = x - ipart; * * compute exp2(ipart) from setting the ieee754 exponent * compute exp2(fpart) using a pade' approximation for x in [-0.5;0.5[ * * the result becomes: exp2(x) = exp2(ipart) * exp2(fpart) */ static const double fm_exp2_q[] = { /* 1.00000000000000000000e0, */ 2.33184211722314911771e2, 4.36821166879210612817e3 }; static const double fm_exp2_p[] = { 2.30933477057345225087e-2, 2.02020656693165307700e1, 1.51390680115615096133e3 }; static double fm_exp2(double x) { double ipart, fpart, px, qx; udi_t epart; ipart = floor(x+0.5); fpart = x - ipart; FM_DOUBLE_INIT_EXP(epart,ipart); x = fpart*fpart; px = fm_exp2_p[0]; px = px*x + fm_exp2_p[1]; qx = x + fm_exp2_q[0]; px = px*x + fm_exp2_p[2]; qx = qx*x + fm_exp2_q[1]; px = px * fpart; x = 1.0 + 2.0*(px/(qx-px)); return epart.f*x; } double fm_exp(double x) { #if defined(__BYTE_ORDER__) #if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ return fm_exp2(FM_DOUBLE_LOG2OFE * x); #endif #endif return exp(x); } /* * Local Variables: * mode: c * compile-command: "make -C .." * c-basic-offset: 4 * fill-column: 76 * indent-tabs-mode: nil * End: */ +} diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index ab2c756dc..c8e46903c 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -1,157 +1,157 @@ #include #define maxelt 5 double fm_exp(double); typedef enum {FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2} lattice_t; typedef struct { int dim1, dim2; double *ptr; } allocatable_double_2d; typedef struct { // cutforce = force cutoff // cutforcesq = force cutoff squared double cutforce,cutforcesq; // Ec_meam = cohesive energy // re_meam = nearest-neighbor distance // Omega_meam = atomic volume // B_meam = bulk modulus // Z_meam = number of first neighbors for reference structure // ielt_meam = atomic number of element // A_meam = adjustable parameter // alpha_meam = sqrt(9*Omega*B/Ec) // rho0_meam = density scaling parameter // delta_meam = heat of formation for alloys // beta[0-3]_meam = electron density constants // t[0-3]_meam = coefficients on densities in Gamma computation // rho_ref_meam = background density for reference structure // ibar_meam(i) = selection parameter for Gamma function for elt i, // lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j) // neltypes = maximum number of element type defined // eltind = index number of pair (similar to Voigt notation; ij = ji) // phir = pair potential function array // phirar[1-6] = spline coeffs // attrac_meam = attraction parameter in Rose energy // repuls_meam = repulsion parameter in Rose energy // nn2_meam = 1 if second nearest neighbors are to be computed, else 0 // zbl_meam = 1 if zbl potential for small r to be use, else 0 // emb_lin_neg = 1 if linear embedding function for rhob to be used, else 0 // bkgd_dyn = 1 if reference densities follows Dynamo, else 0 // Cmin_meam, Cmax_meam = min and max values in screening cutoff // rc_meam = cutoff distance for meam // delr_meam = cutoff region for meam // ebound_meam = factor giving maximum boundary of sceen fcn ellipse // augt1 = flag for whether t1 coefficient should be augmented // ialloy = flag for newer alloy formulation (as in dynamo code) // mix_ref_t = flag to recover "old" way of computing t in reference config // erose_form = selection parameter for form of E_rose function // gsmooth_factor = factor determining length of G smoothing region // vind[23]D = Voight notation index maps for 2 and 3D // v2D,v3D = array of factors to apply for Voight notation // nr,dr = pair function discretization parameters // nrar,rdrar = spline coeff array parameters double Ec_meam[maxelt+1][maxelt+1],re_meam[maxelt+1][maxelt+1]; double Omega_meam[maxelt+1],Z_meam[maxelt+1]; double A_meam[maxelt+1],alpha_meam[maxelt+1][maxelt+1],rho0_meam[maxelt+1]; double delta_meam[maxelt+1][maxelt+1]; double beta0_meam[maxelt+1],beta1_meam[maxelt+1]; double beta2_meam[maxelt+1],beta3_meam[maxelt+1]; double t0_meam[maxelt+1],t1_meam[maxelt+1]; double t2_meam[maxelt+1],t3_meam[maxelt+1]; double rho_ref_meam[maxelt+1]; int ibar_meam[maxelt+1],ielt_meam[maxelt+1]; lattice_t lattce_meam[maxelt+1][maxelt+1]; int nn2_meam[maxelt+1][maxelt+1]; int zbl_meam[maxelt+1][maxelt+1]; int eltind[maxelt+1][maxelt+1]; int neltypes; allocatable_double_2d phir; // [:,:] allocatable_double_2d phirar,phirar1,phirar2,phirar3,phirar4,phirar5,phirar6; // [:,:] double attrac_meam[maxelt+1][maxelt+1],repuls_meam[maxelt+1][maxelt+1]; double Cmin_meam[maxelt+1][maxelt+1][maxelt+1]; double Cmax_meam[maxelt+1][maxelt+1][maxelt+1]; double rc_meam,delr_meam,ebound_meam[maxelt+1][maxelt+1]; int augt1, ialloy, mix_ref_t, erose_form; int emb_lin_neg, bkgd_dyn; double gsmooth_factor; int vind2D[3+1][3+1],vind3D[3+1][3+1][3+1]; int v2D[6+1],v3D[10+1]; int nr,nrar; double dr,rdrar; } meam_data_t; -meam_data_t meam_data; +extern meam_data_t meam_data; // Functions we need for compat #ifndef max #define max(a,b) \ ({ __typeof__ (a) _a = (a); \ __typeof__ (b) _b = (b); \ _a > _b ? _a : _b; }) #endif #ifndef min #define min(a,b) \ ({ __typeof__ (a) _a = (a); \ __typeof__ (b) _b = (b); \ _a < _b ? _a : _b; }) #endif #define iszero(f) (fabs(f)<1e-20) #define setall2d(arr, v) {for(int __i=1;__i<=maxelt;__i++) for(int __j=1;__j<=maxelt;__j++) arr[__i][__j] = v;} #define setall3d(arr, v) {for(int __i=1;__i<=maxelt;__i++) for(int __j=1;__j<=maxelt;__j++) for(int __k=1;__k<=maxelt;__k++) arr[__i][__j][__k] = v;} /* Fortran Array Semantics in C. - Stack-Allocated and global arrays are 1-based, declared as foo[N+1] and simply ignoring the first element - Multi-Dimensional MUST be declared in reverse, so that the order when accessing is the same as in Fortran - arrays that are passed externally via the meam_* functions must use the arr*v() functions below (or be used with 0-based indexing) - allocatable arrays (only global phir*) are actually a struct, and must be accessed with arr2() */ // we receive a pointer to the first element, and F dimensions is ptr(a,b,c) // we know c data structure is ptr[c][b][a] #define arrdim2v(ptr,a,b) \ const int DIM1__##ptr = a; const int DIM2__##ptr = b;\ (void)(DIM1__##ptr); (void)(DIM2__##ptr); #define arrdim3v(ptr,a,b,c) \ const int DIM1__##ptr = a; const int DIM2__##ptr = b; const int DIM3__##ptr = c;\ (void)(DIM1__##ptr); (void)(DIM2__##ptr; (void)(DIM3__##ptr); // access data with same index as used in fortran (1-based) #define arr1v(ptr,i) \ ptr[i-1] #define arr2v(ptr,i,j) \ ptr[(DIM1__##ptr)*(j-1) + (i-1)] #define arr3v(ptr,i,j,k) \ ptr[(i-1) + (j-1)*(DIM1__##ptr) + (k-1)*(DIM1__##ptr)*(DIM2__##ptr)] // allocatable arrays via special type #define allocate_2d(arr,cols,rows) \ arr.dim1 = cols; arr.dim2=rows; arr.ptr=(double*)calloc((size_t)(cols)*(size_t)(rows),sizeof(double)); #define allocated(a) \ (a.ptr!=NULL) #define deallocate(a) \ ({free(a.ptr); a.ptr=NULL;}) // access data with same index as used in fortran (1-based) #define arr2(arr,i,j) \ arr.ptr[(arr.dim1)*(j-1) + (i-1)] diff --git a/src/USER-MEAMC/meam_cleanup.c b/src/USER-MEAMC/meam_cleanup.cpp similarity index 95% rename from src/USER-MEAMC/meam_cleanup.c rename to src/USER-MEAMC/meam_cleanup.cpp index 5c81fc3b4..1beedc9fa 100644 --- a/src/USER-MEAMC/meam_cleanup.c +++ b/src/USER-MEAMC/meam_cleanup.cpp @@ -1,14 +1,16 @@ +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); } +} \ No newline at end of file diff --git a/src/USER-MEAMC/meam_data.c b/src/USER-MEAMC/meam_data.c deleted file mode 100644 index 54ddf81b6..000000000 --- a/src/USER-MEAMC/meam_data.c +++ /dev/null @@ -1,3 +0,0 @@ -#include "meam.h" - -meam_data_t meam_data = {}; \ No newline at end of file diff --git a/src/USER-MEAMC/meam_data.cpp b/src/USER-MEAMC/meam_data.cpp new file mode 100644 index 000000000..c5fd15214 --- /dev/null +++ b/src/USER-MEAMC/meam_data.cpp @@ -0,0 +1,5 @@ +extern "C" { +#include "meam.h" + +meam_data_t meam_data = {}; +} \ No newline at end of file diff --git a/src/USER-MEAMC/meam_dens_final.c b/src/USER-MEAMC/meam_dens_final.cpp similarity index 99% rename from src/USER-MEAMC/meam_dens_final.c rename to src/USER-MEAMC/meam_dens_final.cpp index 1e553866e..bdcff4a7f 100644 --- a/src/USER-MEAMC/meam_dens_final.c +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -1,262 +1,265 @@ +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); // in meam_setup_done void get_shpfcn(double *s /* s(3) */, lattice_t latt); // Extern "C" declaration has the form: // // void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, // int *, int *, int *, // double *, double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, int *); // // Call from pair_meam.cpp has the form: // // meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, // &eng_vdwl,eatom,ntype,type,fmap, // &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], // &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1, // 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) int i, elti; int m; double rhob, G, dG, Gbar, dGbar, gam, shp[3+1], Z; double B, denom, rho_bkgd; // Complete the calculation of density 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); } 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)); } 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; } } } } } //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.c b/src/USER-MEAMC/meam_dens_init.cpp similarity index 99% rename from src/USER-MEAMC/meam_dens_init.c rename to src/USER-MEAMC/meam_dens_init.cpp index 8c1d84f6d..8a17f4e7e 100644 --- a/src/USER-MEAMC/meam_dens_init.c +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -1,504 +1,504 @@ +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); // Extern "C" declaration has the form: // // 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 *); // // // Call from pair_meam.cpp has the form: // // meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0], // &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], // &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], // rho0,&arho1[0][0],&arho2[0][0],arho2b, // &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; // 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); } //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) { 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; drinv = 1.0/meam_data.delr_meam; elti = arr1v(fmap,arr1v(type,i)); if (elti > 0) { 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); 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 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 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; 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]; *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 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 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 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; } //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; 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.c b/src/USER-MEAMC/meam_force.cpp similarity index 99% rename from src/USER-MEAMC/meam_force.c rename to src/USER-MEAMC/meam_force.cpp index d27d8ca7e..b10a76671 100644 --- a/src/USER-MEAMC/meam_force.c +++ b/src/USER-MEAMC/meam_force.cpp @@ -1,551 +1,554 @@ +extern "C"{ #include "meam.h" #include // in meam_setup_done 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); // Extern "C" declaration has the form: // // 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 *); // // Call from pair_meam.cpp has the form: // // meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom, // &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], // &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], // &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], // dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, // &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], // &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) 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; // Compute forces atom i elti = arr1v(fmap,arr1v(type,i)); if (elti > 0) { xitmp = arr2v(x,1,i); yitmp = arr2v(x,2,i); zitmp = arr2v(x,3,i); // Treat each pair for(jn=1; jn<=*numneigh; jn++) { j = arr1v(firstneigh,jn); eltj = arr1v(fmap,arr1v(type,j)); if (!iszero(arr1v(scrfcn,jn)) && eltj > 0) { 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; // 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; 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; } } // 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 } } } // 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.c b/src/USER-MEAMC/meam_setup_done.cpp similarity index 99% rename from src/USER-MEAMC/meam_setup_done.c rename to src/USER-MEAMC/meam_setup_done.cpp index e77702506..0d445f84b 100644 --- a/src/USER-MEAMC/meam_setup_done.c +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -1,945 +1,949 @@ +extern "C" { #include "meam.h" #include void alloyparams(); void compute_pair_meam(); 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); double zbl(double r, int z1, int z2); 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); // in meam_dens_final void G_gam(double Gamma,int ibar,double gsmooth_factor, double *G, int *errorflag); // Declaration in pair_meam.h: // // void meam_setup_done(double *) // // Call from pair_meam.cpp: // // meam_setup_done(&cutmax) // 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; // 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]; // 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; } } } 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; // 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) { 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); // 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_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 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); // 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; // 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; } } } // 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; // Equation numbers below refer to: // I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615 // 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]); get_densref(r,a,b,&rho01,&rho11,&rho21,&rho31,&rho02,&rho12,&rho22,&rho32); // if densities are too small, numerical problems may result; just return zero if (rho01<=1e-14 && rho02<=1e-14) return 0.0; // 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; // 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; // ...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; } } //----------------------------------------------------------------------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.') } } //------------------------------------------------------------------------------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*/; // 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 { // 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.') } } //------------------------------------------------------------------------------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; 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); } //------------------------------------------------------------------------------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; 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); 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; 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 (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]); 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); 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 *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; } } } //--------------------------------------------------------------------- // 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; } //--------------------------------------------------------------------- // 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; 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; } // ----------------------------------------------------------------------- void interpolate_meam(int ind) { int j; double drar; // map to coefficient space 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.; } 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; } } //--------------------------------------------------------------------- // 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; } + + +} diff --git a/src/USER-MEAMC/meam_setup_global.c b/src/USER-MEAMC/meam_setup_global.cpp similarity index 99% rename from src/USER-MEAMC/meam_setup_global.c rename to src/USER-MEAMC/meam_setup_global.cpp index c8c4ef846..5f6dfcb31 100644 --- a/src/USER-MEAMC/meam_setup_global.c +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -1,99 +1,101 @@ +extern "C" { #include "meam.h" #include // // declaration in pair_meam.h: // // void meam_setup_global(int *, int *, double *, int *, double *, double *, // double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, // double *, double *, int *); // // call in pair_meam.cpp: // // meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, // alat,esub,asub,t0,t1,t2,t3,rozero,ibar); // // 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]; meam_data.neltypes = *nelt; 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.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.c b/src/USER-MEAMC/meam_setup_param.cpp similarity index 98% rename from src/USER-MEAMC/meam_setup_param.c rename to src/USER-MEAMC/meam_setup_param.cpp index 2c49481a2..ff89eda96 100644 --- a/src/USER-MEAMC/meam_setup_param.c +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -1,223 +1,226 @@ +extern "C" { #include "meam.h" // // 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; } for (int i=0; i lim)) { *ierr = 3; return; } } } // // Declaration in pair_meam.h: // // void meam_setup_param(int *, double *, int *, int *, int *); // // in pair_meam.cpp // // meam_setup_param(&which,&value,&nindex,index,&errorflag); // // // // The "which" argument corresponds to the index of the "keyword" array // in pair_meam.cpp: // // 0 = Ec_meam // 1 = alpha_meam // 2 = rho0_meam // 3 = delta_meam // 4 = lattce_meam // 5 = attrac_meam // 6 = repuls_meam // 7 = nn2_meam // 8 = Cmin_meam // 9 = Cmax_meam // 10 = rc_meam // 11 = delr_meam // 12 = augt1 // 13 = gsmooth_factor // 14 = re_meam // 15 = ialloy // 16 = mixture_ref_t // 17 = erose_form // 18 = zbl_meam // 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; + 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; // 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; // 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; - int val = (int)value; + 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; // 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; // 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; // 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; // 11 = delr_meam case 11: meam_data.delr_meam = value; break; // 12 = augt1 case 12: meam_data.augt1 = (int)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; // 15 = ialloy case 15: meam_data.ialloy = (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; // 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; // 20 = bkgd_dyn case 20: meam_data.bkgd_dyn = (int)value; break; default: *errorflag = 1; } } + +} \ No newline at end of file