diff --git a/lib/gpu/lal_coul_long.cpp b/lib/gpu/lal_coul_long.cpp index ab0c8ce665..d6e16a9668 100644 --- a/lib/gpu/lal_coul_long.cpp +++ b/lib/gpu/lal_coul_long.cpp @@ -71,9 +71,6 @@ int CoulLongT::init(const int ntypes, double **host_scale, for (int i=0; iucl_device),UCL_READ_ONLY); - lj3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY); - scale.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY); this->atom->type_pack1(ntypes,lj_types,scale,host_write,host_scale); @@ -88,8 +85,7 @@ int CoulLongT::init(const int ntypes, double **host_scale, _g_ewald=g_ewald; _allocated=true; - this->_max_bytes=lj1.row_bytes()+lj3.row_bytes()+scale.row_bytes()+ - sp_cl.row_bytes(); + this->_max_bytes=scale.row_bytes()+sp_cl.row_bytes(); return 0; } @@ -106,8 +102,6 @@ void CoulLongT::clear() { return; _allocated=false; - lj1.clear(); - lj3.clear(); scale.clear(); sp_cl.clear(); this->clear_atomic(); diff --git a/lib/gpu/lal_coul_long.cu b/lib/gpu/lal_coul_long.cu index 48a9681766..12bbbee7d2 100644 --- a/lib/gpu/lal_coul_long.cu +++ b/lib/gpu/lal_coul_long.cu @@ -124,8 +124,7 @@ texture q_tex; #endif __kernel void k_coul_long(const __global numtyp4 *restrict x_, - const __global numtyp4 *restrict lj1, - const __global numtyp4 *restrict lj3, + const __global numtyp *restrict scale, const int lj_types, const __global numtyp *restrict sp_cl_in, const __global int *dev_nbor, @@ -161,6 +160,7 @@ __kernel void k_coul_long(const __global numtyp4 *restrict x_, n_stride,nbor_end,nbor); numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i]; + int itype=ix.w; numtyp qtmp; fetch(qtmp,i,q_tex); for ( ; nbor0) { - e_coul += prefactor*(_erfc-factor_coul); + e_coul += prefactor*(_erfc-factor_coul); } if (vflag>0) { virial[0] += delx*delx*force; @@ -215,8 +217,7 @@ __kernel void k_coul_long(const __global numtyp4 *restrict x_, } __kernel void k_coul_long_fast(const __global numtyp4 *restrict x_, - const __global numtyp4 *restrict lj1_in, - const __global numtyp4 *restrict lj3_in, + const __global numtyp *restrict scale_in, const __global numtyp *restrict sp_cl_in, const __global int *dev_nbor, const __global int *dev_packed, @@ -230,9 +231,12 @@ __kernel void k_coul_long_fast(const __global numtyp4 *restrict x_, int tid, ii, offset; atom_info(t_per_atom,ii,tid,offset); + __local numtyp scale[MAX_SHARED_TYPES*MAX_SHARED_TYPES]; __local numtyp sp_cl[4]; if (tid<4) sp_cl[tid]=sp_cl_in[tid]; + if (tid0) { - e_coul += prefactor*(_erfc-factor_coul); + e_coul += prefactor*(_erfc-factor_coul); } if (vflag>0) { virial[0] += delx*delx*force; diff --git a/lib/gpu/lal_coul_long.h b/lib/gpu/lal_coul_long.h index 7e3f7f36a7..52ed60111b 100644 --- a/lib/gpu/lal_coul_long.h +++ b/lib/gpu/lal_coul_long.h @@ -59,10 +59,6 @@ class CoulLong : public BaseCharge { // --------------------------- TYPE DATA -------------------------- - /// lj1 dummy - UCL_D_Vec lj1; - /// lj3 dummy - UCL_D_Vec lj3; /// scale UCL_D_Vec scale; /// Special Coul values [0-3] diff --git a/lib/meam/Makefile.mingw32-cross b/lib/meam/Makefile.mingw32-cross index 968ac703f3..d4d2dad093 100644 --- a/lib/meam/Makefile.mingw32-cross +++ b/lib/meam/Makefile.mingw32-cross @@ -23,7 +23,7 @@ FILES = $(SRC) Makefile DIR = Obj_mingw32/ LIB = $(DIR)libmeam.a -OBJ = $(SRC:%.F=$(DIR)%.o) +OBJ = $(SRC:%.F=$(DIR)%.o) $(DIR)fm_exp.o # ------ SETTINGS ------ diff --git a/lib/meam/Makefile.mingw64-cross b/lib/meam/Makefile.mingw64-cross index 2c67837650..1a8e97febe 100644 --- a/lib/meam/Makefile.mingw64-cross +++ b/lib/meam/Makefile.mingw64-cross @@ -23,7 +23,7 @@ FILES = $(SRC) Makefile DIR = Obj_mingw64/ LIB = $(DIR)libmeam.a -OBJ = $(SRC:%.F=$(DIR)%.o) +OBJ = $(SRC:%.F=$(DIR)%.o) $(DIR)fm_exp.o # ------ SETTINGS ------ diff --git a/lib/meam/meam_data.F b/lib/meam/meam_data.F index d21c187dd4..719963bd59 100755 --- a/lib/meam/meam_data.F +++ b/lib/meam/meam_data.F @@ -19,7 +19,7 @@ c A_meam = adjustable parameter c alpha_meam = sqrt(9*Omega*B/Ec) c rho0_meam = density scaling parameter c delta_meam = heat of formation for alloys -c beta[0-3]_meam = electron density constants +c beta[0-3]_meam = electron density constants c t[0-3]_meam = coefficients on densities in Gamma computation c rho_ref_meam = background density for reference structure c ibar_meam(i) = selection parameter for Gamma function for elt i, @@ -71,7 +71,7 @@ c nrar,rdrar = spline coeff array parameters $ phirar3(:,:),phirar4(:,:),phirar5(:,:),phirar6(:,:) real*8 attrac_meam(maxelt,maxelt),repuls_meam(maxelt,maxelt) - + real*8 Cmin_meam(maxelt,maxelt,maxelt) real*8 Cmax_meam(maxelt,maxelt,maxelt) real*8 rc_meam,delr_meam,ebound_meam(maxelt,maxelt) diff --git a/lib/meam/meam_dens_final.F b/lib/meam/meam_dens_final.F index c90ca4b5d3..92195dcaf4 100755 --- a/lib/meam/meam_dens_final.F +++ b/lib/meam/meam_dens_final.F @@ -3,7 +3,7 @@ c c void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, c int *, int *, int *, c double *, double *, double *, double *, double *, double *, -c double *, double *, double *, double *, double *, double *, +c double *, double *, double *, double *, double *, double *, c double *, double *, double *, double *, double *, int *); c c Call from pair_meam.cpp has the form: @@ -69,7 +69,7 @@ c Complete the calculation of density do m = 1,10 rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i) enddo - + if( rho0(i) .gt. 0.0 ) then if (ialloy.eq.1) then t_ave(1,i) = t_ave(1,i)/tsq_ave(1,i) @@ -85,10 +85,10 @@ c Complete the calculation of density t_ave(3,i) = t_ave(3,i)/rho0(i) endif endif - - Gamma(i) = t_ave(1,i)*rho1(i) + + Gamma(i) = t_ave(1,i)*rho1(i) $ + t_ave(2,i)*rho2(i) + t_ave(3,i)*rho3(i) - + if( rho0(i) .gt. 0.0 ) then Gamma(i) = Gamma(i)/(rho0(i)*rho0(i)) end if @@ -137,9 +137,9 @@ c Complete the calculation of density denom = 1.d0/rho_bkgd call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG) - + dGamma1(i) = (G - 2*dG*Gamma(i))*denom - + if( rho0(i) .ne. 0.d0 ) then dGamma2(i) = (dG/rho0(i))*denom else @@ -156,7 +156,7 @@ c included for backward compatibility endif B = A_meam(elti)*Ec_meam(elti,elti) - + if( rhob .ne. 0.d0 ) then if (emb_lin_neg.eq.1 .and. rhob.le.0) then fp(i) = -B @@ -188,7 +188,7 @@ c included for backward compatibility endif endif enddo - + return end @@ -198,7 +198,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Compute G(Gamma) based on selection flag ibar: c 0 => G = sqrt(1+Gamma) c 1 => G = exp(Gamma/2) -c 2 => not implemented +c 2 => not implemented c 3 => G = 2/(1+exp(-Gamma)) c 4 => G = sqrt(1+Gamma) c -5 => G = +-sqrt(abs(1+Gamma)) @@ -241,7 +241,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Compute G(Gamma) and dG(gamma) based on selection flag ibar: c 0 => G = sqrt(1+Gamma) c 1 => G = fm_exp(Gamma/2) -c 2 => not implemented +c 2 => not implemented c 3 => G = 2/(1+fm_exp(-Gamma)) c 4 => G = sqrt(1+Gamma) c -5 => G = +-sqrt(abs(1+Gamma)) diff --git a/lib/meam/meam_dens_init.F b/lib/meam/meam_dens_init.F index 227d3208f3..2ca2558135 100755 --- a/lib/meam/meam_dens_init.F +++ b/lib/meam/meam_dens_init.F @@ -1,23 +1,23 @@ c Extern "C" declaration has the form: -c +c c void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *, c int *, int *, int *, int *, c double *, double *, double *, double *, double *, double *, c double *, double *, double *, double *, double *, int *); -c -c +c +c c Call from pair_meam.cpp has the form: -c +c c meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0], c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], c rho0,&arho1[0][0],&arho2[0][0],arho2b, c &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag); -c +c subroutine meam_dens_init(i, nmax, $ ntype, type, fmap, x, - $ numneigh, firstneigh, + $ numneigh, firstneigh, $ numneigh_full, firstneigh_full, $ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b, $ arho3, arho3b, t_ave, tsq_ave, errorflag) @@ -51,7 +51,7 @@ c Compute screening function and derivatives $ ntype, type, fmap) c Calculate intermediate density terms to be communicated - call calc_rho1(i, nmax, ntype, type, fmap, x, + call calc_rho1(i, nmax, ntype, type, fmap, x, $ numneigh, firstneigh, $ scrfcn, fcpair, rho0, arho1, arho2, arho2b, $ arho3, arho3b, t_ave, tsq_ave) @@ -61,7 +61,7 @@ c Calculate intermediate density terms to be communicated cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - subroutine getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x, + subroutine getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x, $ numneigh, firstneigh, $ numneigh_full, firstneigh_full, $ ntype, type, fmap) @@ -99,13 +99,13 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc xitmp = x(1,i) yitmp = x(2,i) zitmp = x(3,i) - + do jn = 1,numneigh j = firstneigh(jn) eltj = fmap(type(j)) if (eltj.gt.0) then - + c First compute screening function itself, sij xjtmp = x(1,j) yjtmp = x(2,j) @@ -127,7 +127,7 @@ c First compute screening function itself, sij fcij = fc dfcij = dfc*drinv endif - + c Now compute derivatives dscrfcn(jn) = 0.d0 sfcij = sij*fcij @@ -180,23 +180,23 @@ c Note that we never have 0j, set equal to the iebound, -c atom k definitely lies outside the screening function ellipse (so +c atom k definitely lies outside the screening function ellipse (so c there is no need to calculate its effects). Here, compute it for all c triplets (i,j,k) so that ebound(i,j) is the maximized over k do i = 1,neltypes @@ -237,7 +237,7 @@ c phi_aa $ (-Z2*scrn/Z1)**n * phi_meam(rarat*arat**n,a,a) enddo endif - + c phi_bb phibb = phi_meam(rarat,b,b) call get_Zij(Z1,lattce_meam(b,b)) @@ -250,7 +250,7 @@ c phi_bb $ (-Z2*scrn/Z1)**n * phi_meam(rarat*arat**n,b,b) enddo endif - + if (lattce_meam(a,b).eq.'b1'. $ or.lattce_meam(a,b).eq.'b2') then c Add contributions to the B1 or B2 potential @@ -276,11 +276,11 @@ c atoms. call get_sijk(C,b,b,a,s221) S11 = s111 * s111 * s112 * s112 S22 = s221**4 - phir(j,nv2) = phir(j,nv2) - + phir(j,nv2) = phir(j,nv2) - $ 0.75*S11*phiaa - 0.25*S22*phibb endif - + else nmax = 10 do n = 1,nmax @@ -307,7 +307,7 @@ c endif phir(j,nv2) = frac*phir(j,nv2) + (1-frac)*phizbl endif endif - + enddo c call interpolation @@ -320,14 +320,14 @@ c call interpolation end -c----------------------------------------------------------------------c +c----------------------------------------------------------------------c c Compute MEAM pair potential for distance r, element types a and b c real*8 recursive function phi_meam(r,a,b)result(phi_m) use meam_data implicit none - + integer a,b real*8 r real*8 a1,a2,a12 @@ -375,15 +375,15 @@ c calculate average weighting factors for the reference structure t31av = t3_meam(a) t32av = t3_meam(b) else - scalfac = 1.0/(rho01+rho02) + scalfac = 1.0/(rho01+rho02) t11av = scalfac*(t1_meam(a)*rho01 + t1_meam(b)*rho02) - t12av = t11av + t12av = t11av t21av = scalfac*(t2_meam(a)*rho01 + t2_meam(b)*rho02) - t22av = t21av + t22av = t21av t31av = scalfac*(t3_meam(a)*rho01 + t3_meam(b)*rho02) - t32av = t31av + t32av = t31av endif - else + else c average weighting factors for the reference structure, eqn. I.8 call get_tavref(t11av,t21av,t31av,t12av,t22av,t32av, $ t1_meam(a),t2_meam(a),t3_meam(a), @@ -411,10 +411,10 @@ c for c11b structure, calculate background electron densities rhobar1 = (Z12*rho02)**2 + 2.0/3.0*t22av*rho22**2 rhobar1 = sqrt(rhobar1) endif - else + else c for other structures, use formalism developed in Huang's paper c -c composition-dependent scaling, equation I.7 +c composition-dependent scaling, equation I.7 c If using mixing rule for t, apply to reference structure; else c use precomputed values if (mix_ref_t.eq.1) then @@ -508,8 +508,8 @@ c account for second neighbor a-a potential here... enddo endif phi_m = Eu/3. - F1/4. - F2/12. - phiaa - else -c + else +c c potential is computed from Rose function and embedding energy phi_m = (2*Eu - F1 - F2)/Z12 c @@ -519,7 +519,7 @@ c if r = 0, just return 0 if (r.eq.0.d0) then phi_m = 0.d0 endif - + return end @@ -546,7 +546,7 @@ c loop over element types call G_gam(gam,ibar_meam(a),gsmooth_factor, $ Gbar,errorflag) endif - + c The zeroth order density in the reference structure, with c equilibrium spacing, is just the number of first neighbors times c the rho0_meam coefficient... @@ -569,7 +569,7 @@ c screening) return end -c----------------------------------------------------------------------c +c----------------------------------------------------------------------c c Shape factors for various configurations subroutine get_shpfcn(s,latt) implicit none @@ -599,7 +599,7 @@ c call error('Lattice not defined in get_shpfcn.') endif return end -c------------------------------------------------------------------------------c +c------------------------------------------------------------------------------c c Average weighting factors for the reference structure subroutine get_tavref(t11av,t21av,t31av,t12av,t22av,t32av, $ t11,t21,t31,t12,t22,t32, @@ -651,7 +651,7 @@ c call error('Lattice not defined in get_tavref.') endif return end -c------------------------------------------------------------------------------c +c------------------------------------------------------------------------------c c Number of neighbors for the reference structure subroutine get_Zij(Zij,latt) implicit none @@ -681,8 +681,8 @@ c call error('Lattice not defined in get_Zij.') return end -c------------------------------------------------------------------------------c -c Zij2 = number of second neighbors, a = distance ratio R1/R2, and S = second +c------------------------------------------------------------------------------c +c Zij2 = number of second neighbors, a = distance ratio R1/R2, and S = second c neighbor screening function for lattice type "latt" subroutine get_Zij2(Zij2,a,S,latt,cmin,cmax) @@ -734,7 +734,7 @@ c this really shouldn't be allowed; make sure screening is zero c call error('Lattice not defined in get_Zij2.') endif -c Compute screening for each first neighbor +c Compute screening for each first neighbor C = 4.d0/(a*a) - 1.d0 x = (C-cmin)/(cmax-cmin) call fcut(x,sijk) @@ -744,8 +744,8 @@ c There are numscr first neighbors screening the second neighbors return end - -c------------------------------------------------------------------------------c + +c------------------------------------------------------------------------------c subroutine get_sijk(C,i,j,k,sijk) use meam_data implicit none @@ -757,7 +757,7 @@ c------------------------------------------------------------------------------c return end -c------------------------------------------------------------------------------c +c------------------------------------------------------------------------------c c Calculate density functions, assuming reference configuration subroutine get_densref(r,a,b,rho01,rho11,rho21,rho31, $ rho02,rho12,rho22,rho32) @@ -787,7 +787,7 @@ c Calculate density functions, assuming reference configuration rhoa32 = rho0_meam(b)*fm_exp(-beta3_meam(b)*a2) lat = lattce_meam(a,b) - + rho11 = 0.d0 rho21 = 0.d0 rho31 = 0.d0 @@ -807,13 +807,13 @@ c Calculate density functions, assuming reference configuration rho01 = 6*rhoa02 rho02 = 6*rhoa01 else if (lat.eq.'dia') then - rho01 = 4*rhoa02 - rho02 = 4*rhoa01 + rho01 = 4*rhoa02 + rho02 = 4*rhoa01 rho31 = 32.d0/9.d0*rhoa32*rhoa32 rho32 = 32.d0/9.d0*rhoa31*rhoa31 else if (lat.eq.'hcp') then - rho01 = 12*rhoa02 - rho02 = 12*rhoa01 + rho01 = 12*rhoa02 + rho02 = 12*rhoa01 rho31 = 1.d0/3.d0*rhoa32*rhoa32 rho32 = 1.d0/3.d0*rhoa31*rhoa31 else if (lat.eq.'dim') then @@ -830,7 +830,7 @@ c Calculate density functions, assuming reference configuration rho01 = rhoa01 rho02 = rhoa02 rho11 = rhoa11 - rho12 = rhoa12 + rho12 = rhoa12 rho21 = rhoa21 rho22 = rhoa22 rho31 = rhoa31 @@ -858,10 +858,10 @@ c call error('Lattice not defined in get_densref.') call get_Zij2(Zij2nn,arat,scrn,lat, $ Cmin_meam(a,a,b),Cmax_meam(a,a,b)) - + a1 = arat*r/re_meam(a,a) - 1.d0 a2 = arat*r/re_meam(b,b) - 1.d0 - + rhoa01nn = rho0_meam(a)*fm_exp(-beta0_meam(a)*a1) rhoa02nn = rho0_meam(b)*fm_exp(-beta0_meam(b)*a2) @@ -895,41 +895,41 @@ c Assume Zij2nn and arat don't depend on order, but scrn might return end -c--------------------------------------------------------------------- -c Compute ZBL potential -c - real*8 function zbl(r,z1,z2) - use meam_data , only : fm_exp - implicit none - integer i,z1,z2 - real*8 r,c,d,a,azero,cc,x - dimension c(4),d(4) - data c /0.028171,0.28022,0.50986,0.18175/ - data d /0.20162,0.40290,0.94229,3.1998/ - data azero /0.4685/ - data cc /14.3997/ -c azero = (9pi^2/128)^1/3 (0.529) Angstroms - a = azero/(z1**0.23+z2**0.23) - zbl = 0.0 - x = r/a - do i=1,4 - zbl = zbl + c(i)*fm_exp(-d(i)*x) +c--------------------------------------------------------------------- +c Compute ZBL potential +c + real*8 function zbl(r,z1,z2) + use meam_data , only : fm_exp + implicit none + integer i,z1,z2 + real*8 r,c,d,a,azero,cc,x + dimension c(4),d(4) + data c /0.028171,0.28022,0.50986,0.18175/ + data d /0.20162,0.40290,0.94229,3.1998/ + data azero /0.4685/ + data cc /14.3997/ +c azero = (9pi^2/128)^1/3 (0.529) Angstroms + a = azero/(z1**0.23+z2**0.23) + zbl = 0.0 + x = r/a + do i=1,4 + zbl = zbl + c(i)*fm_exp(-d(i)*x) enddo if (r.gt.0.d0) zbl = zbl*z1*z2/r*cc - return - end - + return + end + c--------------------------------------------------------------------- c Compute Rose energy function, I.16 c real*8 function erose(r,re,alpha,Ec,repuls,attrac,form) - use meam_data , only : fm_exp + use meam_data , only : fm_exp implicit none real*8 r,re,alpha,Ec,repuls,attrac,astar,a3 integer form erose = 0.d0 - + if (r.gt.0.d0) then astar = alpha * (r/re - 1.d0) a3 = 0.d0 @@ -950,7 +950,7 @@ c return end - + c ----------------------------------------------------------------------- subroutine interpolate_meam(ind) @@ -980,7 +980,7 @@ c phir interp phirar1(j,ind) = ((phirar(j-2,ind)-phirar(j+2,ind)) + $ 8.0D0*(phirar(j+1,ind)-phirar(j-1,ind)))/12. enddo - + do j = 1,nrar-1 phirar2(j,ind) = 3.0D0*(phirar(j+1,ind)-phirar(j,ind)) - $ 2.0D0*phirar1(j,ind) - phirar1(j+1,ind) @@ -989,7 +989,7 @@ c phir interp enddo phirar2(nrar,ind) = 0.0D0 phirar3(nrar,ind) = 0.0D0 - + do j = 1,nrar phirar4(j,ind) = phirar1(j,ind)/drar phirar5(j,ind) = 2.0D0*phirar2(j,ind)/drar @@ -1014,7 +1014,7 @@ c kk = min(kk,nrar-1) pp = pp - kk pp = min(pp,1.0D0) - compute_phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp + compute_phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp $ + phirar1(kk,ind))*pp + phirar(kk,ind) return diff --git a/lib/meam/meam_setup_global.F b/lib/meam/meam_setup_global.F index 97c4a8427e..d11dec5a4a 100755 --- a/lib/meam/meam_setup_global.F +++ b/lib/meam/meam_setup_global.F @@ -14,7 +14,7 @@ c c subroutine meam_setup_global(nelt, lat, z, ielement, atwt, alpha, - $ b0, b1, b2, b3, alat, esub, asub, + $ b0, b1, b2, b3, alat, esub, asub, $ t0, t1, t2, t3, rozero, ibar) use meam_data diff --git a/lib/meam/meam_setup_param.F b/lib/meam/meam_setup_param.F index e0a755cd6f..cfe7430285 100755 --- a/lib/meam/meam_setup_param.F +++ b/lib/meam/meam_setup_param.F @@ -1,5 +1,5 @@ c -c do a sanity check on index parameters +c do a sanity check on index parameters subroutine meam_checkindex(num,lim,nidx,idx,ierr) implicit none integer i,num,lim,nidx,idx(3),ierr @@ -18,20 +18,20 @@ c do a sanity check on index parameters enddo end -c +c c Declaration in pair_meam.h: -c +c c void meam_setup_param(int *, double *, int *, int *, int *); -c +c c Call in pair_meam.cpp -c +c c meam_setup_param(&which,&value,&nindex,index,&errorflag); -c -c -c +c +c +c c The "which" argument corresponds to the index of the "keyword" array c in pair_meam.cpp: -c +c c 0 = Ec_meam c 1 = alpha_meam c 2 = rho0_meam @@ -54,7 +54,7 @@ c 18 = zbl_meam c 19 = emb_lin_neg c 20 = bkgd_dyn - subroutine meam_setup_param(which, value, nindex, + subroutine meam_setup_param(which, value, nindex, $ index, errorflag) use meam_data @@ -77,7 +77,7 @@ c 1 = alpha_meam call meam_checkindex(2,maxelt,nindex,index,errorflag) if (errorflag.ne.0) return alpha_meam(index(1),index(2)) = value - + c 2 = rho0_meam else if (which.eq.2) then call meam_checkindex(1,maxelt,nindex,index,errorflag)