ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14867 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2016-04-22 17:51:49 +00:00
parent 129796adc2
commit 184d5dc0f0
12 changed files with 168 additions and 171 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
lib/gpu/lal_coul_long.cpp
View File

@ -71,9 +71,6 @@ int CoulLongT::init(const int ntypes, double **host_scale,
for (int i=0; i<lj_types*lj_types; i++) for (int i=0; i<lj_types*lj_types; i++)
host_write[i]=0.0; host_write[i]=0.0;
lj1.alloc(lj_types*lj_types,*(this->ucl_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); scale.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
this->atom->type_pack1(ntypes,lj_types,scale,host_write,host_scale); 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; _g_ewald=g_ewald;
_allocated=true; _allocated=true;
this->_max_bytes=lj1.row_bytes()+lj3.row_bytes()+scale.row_bytes()+ this->_max_bytes=scale.row_bytes()+sp_cl.row_bytes();
sp_cl.row_bytes();
return 0; return 0;
} }
@ -106,8 +102,6 @@ void CoulLongT::clear() {
return; return;
_allocated=false; _allocated=false;
lj1.clear();
lj3.clear();
scale.clear(); scale.clear();
sp_cl.clear(); sp_cl.clear();
this->clear_atomic(); this->clear_atomic();

25
lib/gpu/lal_coul_long.cu
View File

@ -124,8 +124,7 @@ texture<int2> q_tex;
#endif #endif
__kernel void k_coul_long(const __global numtyp4 *restrict x_, __kernel void k_coul_long(const __global numtyp4 *restrict x_,
const __global numtyp4 *restrict lj1, const __global numtyp *restrict scale,
const __global numtyp4 *restrict lj3,
const int lj_types, const int lj_types,
const __global numtyp *restrict sp_cl_in, const __global numtyp *restrict sp_cl_in,
const __global int *dev_nbor, const __global int *dev_nbor,
@ -161,6 +160,7 @@ __kernel void k_coul_long(const __global numtyp4 *restrict x_,
n_stride,nbor_end,nbor); n_stride,nbor_end,nbor);
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i]; numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
int itype=ix.w;
numtyp qtmp; fetch(qtmp,i,q_tex); numtyp qtmp; fetch(qtmp,i,q_tex);
for ( ; nbor<nbor_end; nbor+=n_stride) { for ( ; nbor<nbor_end; nbor+=n_stride) {
@ -171,6 +171,7 @@ __kernel void k_coul_long(const __global numtyp4 *restrict x_,
j &= NEIGHMASK; j &= NEIGHMASK;
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j]; numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
int jtype=jx.w;
// Compute r12 // Compute r12
numtyp delx = ix.x-jx.x; numtyp delx = ix.x-jx.x;
@ -178,6 +179,7 @@ __kernel void k_coul_long(const __global numtyp4 *restrict x_,
numtyp delz = ix.z-jx.z; numtyp delz = ix.z-jx.z;
numtyp rsq = delx*delx+dely*dely+delz*delz; numtyp rsq = delx*delx+dely*dely+delz*delz;
int mtype=itype*lj_types+jtype;
if (rsq < cut_coulsq) { if (rsq < cut_coulsq) {
numtyp r2inv=ucl_recip(rsq); numtyp r2inv=ucl_recip(rsq);
numtyp force, prefactor, _erfc; numtyp force, prefactor, _erfc;
@ -188,7 +190,7 @@ __kernel void k_coul_long(const __global numtyp4 *restrict x_,
numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij); numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
_erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; _erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
fetch(prefactor,j,q_tex); fetch(prefactor,j,q_tex);
prefactor *= qqrd2e * qtmp/r; prefactor *= qqrd2e * scale[mtype] * qtmp/r;
force = prefactor * (_erfc + EWALD_F*grij*expm2-factor_coul) * r2inv; force = prefactor * (_erfc + EWALD_F*grij*expm2-factor_coul) * r2inv;
f.x+=delx*force; f.x+=delx*force;
@ -196,7 +198,7 @@ __kernel void k_coul_long(const __global numtyp4 *restrict x_,
f.z+=delz*force; f.z+=delz*force;
if (eflag>0) { if (eflag>0) {
e_coul += prefactor*(_erfc-factor_coul); e_coul += prefactor*(_erfc-factor_coul);
} }
if (vflag>0) { if (vflag>0) {
virial[0] += delx*delx*force; 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_, __kernel void k_coul_long_fast(const __global numtyp4 *restrict x_,
const __global numtyp4 *restrict lj1_in, const __global numtyp *restrict scale_in,
const __global numtyp4 *restrict lj3_in,
const __global numtyp *restrict sp_cl_in, const __global numtyp *restrict sp_cl_in,
const __global int *dev_nbor, const __global int *dev_nbor,
const __global int *dev_packed, const __global int *dev_packed,
@ -230,9 +231,12 @@ __kernel void k_coul_long_fast(const __global numtyp4 *restrict x_,
int tid, ii, offset; int tid, ii, offset;
atom_info(t_per_atom,ii,tid,offset); atom_info(t_per_atom,ii,tid,offset);
__local numtyp scale[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
__local numtyp sp_cl[4]; __local numtyp sp_cl[4];
if (tid<4) if (tid<4)
sp_cl[tid]=sp_cl_in[tid]; sp_cl[tid]=sp_cl_in[tid];
if (tid<MAX_SHARED_TYPES*MAX_SHARED_TYPES)
scale[tid]=scale_in[tid];
acctyp e_coul=(acctyp)0; acctyp e_coul=(acctyp)0;
acctyp4 f; acctyp4 f;
@ -252,6 +256,8 @@ __kernel void k_coul_long_fast(const __global numtyp4 *restrict x_,
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i]; numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
numtyp qtmp; fetch(qtmp,i,q_tex); numtyp qtmp; fetch(qtmp,i,q_tex);
int iw=ix.w;
int itype=fast_mul((int)MAX_SHARED_TYPES,iw);
for ( ; nbor<nbor_end; nbor+=n_stride) { for ( ; nbor<nbor_end; nbor+=n_stride) {
int j=dev_packed[nbor]; int j=dev_packed[nbor];
@ -261,6 +267,7 @@ __kernel void k_coul_long_fast(const __global numtyp4 *restrict x_,
j &= NEIGHMASK; j &= NEIGHMASK;
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j]; numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
int mtype=itype+jx.w;
// Compute r12 // Compute r12
numtyp delx = ix.x-jx.x; numtyp delx = ix.x-jx.x;
@ -278,15 +285,15 @@ __kernel void k_coul_long_fast(const __global numtyp4 *restrict x_,
numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij); numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
_erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; _erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
fetch(prefactor,j,q_tex); fetch(prefactor,j,q_tex);
prefactor *= qqrd2e * qtmp/r; prefactor *= qqrd2e * scale[mtype] * qtmp/r;
force = prefactor * (_erfc + EWALD_F*grij*expm2-factor_coul) * r2inv; force = prefactor*(_erfc + EWALD_F*grij*expm2-factor_coul) * r2inv;
f.x+=delx*force; f.x+=delx*force;
f.y+=dely*force; f.y+=dely*force;
f.z+=delz*force; f.z+=delz*force;
if (eflag>0) { if (eflag>0) {
e_coul += prefactor*(_erfc-factor_coul); e_coul += prefactor*(_erfc-factor_coul);
} }
if (vflag>0) { if (vflag>0) {
virial[0] += delx*delx*force; virial[0] += delx*delx*force;

4
lib/gpu/lal_coul_long.h
View File

@ -59,10 +59,6 @@ class CoulLong : public BaseCharge<numtyp, acctyp> {
// --------------------------- TYPE DATA -------------------------- // --------------------------- TYPE DATA --------------------------
/// lj1 dummy
UCL_D_Vec<numtyp4> lj1;
/// lj3 dummy
UCL_D_Vec<numtyp4> lj3;
/// scale /// scale
UCL_D_Vec<numtyp> scale; UCL_D_Vec<numtyp> scale;
/// Special Coul values [0-3] /// Special Coul values [0-3]

2
lib/meam/Makefile.mingw32-cross
View File

@ -23,7 +23,7 @@ FILES = $(SRC) Makefile
DIR = Obj_mingw32/ DIR = Obj_mingw32/
LIB = $(DIR)libmeam.a LIB = $(DIR)libmeam.a
OBJ = $(SRC:%.F=$(DIR)%.o) OBJ = $(SRC:%.F=$(DIR)%.o) $(DIR)fm_exp.o
# ------ SETTINGS ------ # ------ SETTINGS ------

2
lib/meam/Makefile.mingw64-cross
View File

@ -23,7 +23,7 @@ FILES = $(SRC) Makefile
DIR = Obj_mingw64/ DIR = Obj_mingw64/
LIB = $(DIR)libmeam.a LIB = $(DIR)libmeam.a
OBJ = $(SRC:%.F=$(DIR)%.o) OBJ = $(SRC:%.F=$(DIR)%.o) $(DIR)fm_exp.o
# ------ SETTINGS ------ # ------ SETTINGS ------

4
lib/meam/meam_data.F
View File

@ -19,7 +19,7 @@ c A_meam = adjustable parameter
c alpha_meam = sqrt(9*Omega*B/Ec) c alpha_meam = sqrt(9*Omega*B/Ec)
c rho0_meam = density scaling parameter c rho0_meam = density scaling parameter
c delta_meam = heat of formation for alloys 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 t[0-3]_meam = coefficients on densities in Gamma computation
c rho_ref_meam = background density for reference structure c rho_ref_meam = background density for reference structure
c ibar_meam(i) = selection parameter for Gamma function for elt i, 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(:,:) $ phirar3(:,:),phirar4(:,:),phirar5(:,:),phirar6(:,:)
real*8 attrac_meam(maxelt,maxelt),repuls_meam(maxelt,maxelt) real*8 attrac_meam(maxelt,maxelt),repuls_meam(maxelt,maxelt)
real*8 Cmin_meam(maxelt,maxelt,maxelt) real*8 Cmin_meam(maxelt,maxelt,maxelt)
real*8 Cmax_meam(maxelt,maxelt,maxelt) real*8 Cmax_meam(maxelt,maxelt,maxelt)
real*8 rc_meam,delr_meam,ebound_meam(maxelt,maxelt) real*8 rc_meam,delr_meam,ebound_meam(maxelt,maxelt)

22
lib/meam/meam_dens_final.F
View File

@ -3,7 +3,7 @@ c
c void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, c void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *,
c int *, int *, int *, 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 *, double *,
c double *, double *, double *, double *, double *, int *); c double *, double *, double *, double *, double *, int *);
c c
c Call from pair_meam.cpp has the form: c Call from pair_meam.cpp has the form:
@ -69,7 +69,7 @@ c Complete the calculation of density
do m = 1,10 do m = 1,10
rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i) rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i)
enddo enddo
if( rho0(i) .gt. 0.0 ) then if( rho0(i) .gt. 0.0 ) then
if (ialloy.eq.1) then if (ialloy.eq.1) then
t_ave(1,i) = t_ave(1,i)/tsq_ave(1,i) 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) t_ave(3,i) = t_ave(3,i)/rho0(i)
endif endif
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) $ + t_ave(2,i)*rho2(i) + t_ave(3,i)*rho3(i)
if( rho0(i) .gt. 0.0 ) then if( rho0(i) .gt. 0.0 ) then
Gamma(i) = Gamma(i)/(rho0(i)*rho0(i)) Gamma(i) = Gamma(i)/(rho0(i)*rho0(i))
end if end if
@ -137,9 +137,9 @@ c Complete the calculation of density
denom = 1.d0/rho_bkgd denom = 1.d0/rho_bkgd
call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG) call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG)
dGamma1(i) = (G - 2*dG*Gamma(i))*denom dGamma1(i) = (G - 2*dG*Gamma(i))*denom
if( rho0(i) .ne. 0.d0 ) then if( rho0(i) .ne. 0.d0 ) then
dGamma2(i) = (dG/rho0(i))*denom dGamma2(i) = (dG/rho0(i))*denom
else else
@ -156,7 +156,7 @@ c included for backward compatibility
endif endif
B = A_meam(elti)*Ec_meam(elti,elti) B = A_meam(elti)*Ec_meam(elti,elti)
if( rhob .ne. 0.d0 ) then if( rhob .ne. 0.d0 ) then
if (emb_lin_neg.eq.1 .and. rhob.le.0) then if (emb_lin_neg.eq.1 .and. rhob.le.0) then
fp(i) = -B fp(i) = -B
@ -188,7 +188,7 @@ c included for backward compatibility
endif endif
endif endif
enddo enddo
return return
end end
@ -198,7 +198,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Compute G(Gamma) based on selection flag ibar: c Compute G(Gamma) based on selection flag ibar:
c 0 => G = sqrt(1+Gamma) c 0 => G = sqrt(1+Gamma)
c 1 => G = exp(Gamma/2) c 1 => G = exp(Gamma/2)
c 2 => not implemented c 2 => not implemented
c 3 => G = 2/(1+exp(-Gamma)) c 3 => G = 2/(1+exp(-Gamma))
c 4 => G = sqrt(1+Gamma) c 4 => G = sqrt(1+Gamma)
c -5 => G = +-sqrt(abs(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 Compute G(Gamma) and dG(gamma) based on selection flag ibar:
c 0 => G = sqrt(1+Gamma) c 0 => G = sqrt(1+Gamma)
c 1 => G = fm_exp(Gamma/2) c 1 => G = fm_exp(Gamma/2)
c 2 => not implemented c 2 => not implemented
c 3 => G = 2/(1+fm_exp(-Gamma)) c 3 => G = 2/(1+fm_exp(-Gamma))
c 4 => G = sqrt(1+Gamma) c 4 => G = sqrt(1+Gamma)
c -5 => G = +-sqrt(abs(1+Gamma)) c -5 => G = +-sqrt(abs(1+Gamma))

44
lib/meam/meam_dens_init.F
View File

@ -1,23 +1,23 @@
c Extern "C" declaration has the form: c Extern "C" declaration has the form:
c c
c void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *, c void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *,
c int *, int *, int *, int *, c int *, int *, int *, int *,
c double *, double *, double *, double *, double *, double *, c double *, double *, double *, double *, double *, double *,
c double *, double *, double *, double *, double *, int *); c double *, double *, double *, double *, double *, int *);
c c
c c
c Call from pair_meam.cpp has the form: c Call from pair_meam.cpp has the form:
c c
c meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0], c meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0],
c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],
c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
c rho0,&arho1[0][0],&arho2[0][0],arho2b, 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 &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag);
c c
subroutine meam_dens_init(i, nmax, subroutine meam_dens_init(i, nmax,
$ ntype, type, fmap, x, $ ntype, type, fmap, x,
$ numneigh, firstneigh, $ numneigh, firstneigh,
$ numneigh_full, firstneigh_full, $ numneigh_full, firstneigh_full,
$ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b, $ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b,
$ arho3, arho3b, t_ave, tsq_ave, errorflag) $ arho3, arho3b, t_ave, tsq_ave, errorflag)
@ -51,7 +51,7 @@ c Compute screening function and derivatives
$ ntype, type, fmap) $ ntype, type, fmap)
c Calculate intermediate density terms to be communicated 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, $ numneigh, firstneigh,
$ scrfcn, fcpair, rho0, arho1, arho2, arho2b, $ scrfcn, fcpair, rho0, arho1, arho2, arho2b,
$ arho3, arho3b, t_ave, tsq_ave) $ arho3, arho3b, t_ave, tsq_ave)
@ -61,7 +61,7 @@ c Calculate intermediate density terms to be communicated
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x, subroutine getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x,
$ numneigh, firstneigh, $ numneigh, firstneigh,
$ numneigh_full, firstneigh_full, $ numneigh_full, firstneigh_full,
$ ntype, type, fmap) $ ntype, type, fmap)
@ -99,13 +99,13 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
xitmp = x(1,i) xitmp = x(1,i)
yitmp = x(2,i) yitmp = x(2,i)
zitmp = x(3,i) zitmp = x(3,i)
do jn = 1,numneigh do jn = 1,numneigh
j = firstneigh(jn) j = firstneigh(jn)
eltj = fmap(type(j)) eltj = fmap(type(j))
if (eltj.gt.0) then if (eltj.gt.0) then
c First compute screening function itself, sij c First compute screening function itself, sij
xjtmp = x(1,j) xjtmp = x(1,j)
yjtmp = x(2,j) yjtmp = x(2,j)
@ -127,7 +127,7 @@ c First compute screening function itself, sij
fcij = fc fcij = fc
dfcij = dfc*drinv dfcij = dfc*drinv
endif endif
c Now compute derivatives c Now compute derivatives
dscrfcn(jn) = 0.d0 dscrfcn(jn) = 0.d0
sfcij = sij*fcij sfcij = sij*fcij
@ -180,23 +180,23 @@ c Note that we never have 0<cikj<Cmin here, else sij=0 (rejected above)
coef2 = sij*dfcij/rij coef2 = sij*dfcij/rij
dscrfcn(jn) = dscrfcn(jn)*coef1 - coef2 dscrfcn(jn) = dscrfcn(jn)*coef1 - coef2
100 continue 100 continue
scrfcn(jn) = sij scrfcn(jn) = sij
fcpair(jn) = fcij fcpair(jn) = fcij
endif endif
enddo enddo
endif endif
return return
end end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine calc_rho1(i, nmax, ntype, type, fmap, x, subroutine calc_rho1(i, nmax, ntype, type, fmap, x,
$ numneigh, firstneigh, $ numneigh, firstneigh,
$ scrfcn, fcpair, rho0, arho1, arho2, arho2b, $ scrfcn, fcpair, rho0, arho1, arho2, arho2b,
$ arho3, arho3b, t_ave, tsq_ave) $ arho3, arho3b, t_ave, tsq_ave)
@ -236,7 +236,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
delij(1) = x(1,j) - xtmp delij(1) = x(1,j) - xtmp
delij(2) = x(2,j) - ytmp delij(2) = x(2,j) - ytmp
delij(3) = x(3,j) - ztmp delij(3) = x(3,j) - ztmp
rij2 = delij(1)*delij(1) + delij(2)*delij(2) rij2 = delij(1)*delij(1) + delij(2)*delij(2)
$ + delij(3)*delij(3) $ + delij(3)*delij(3)
if (rij2.lt.cutforcesq) then if (rij2.lt.cutforcesq) then
eltj = fmap(type(j)) eltj = fmap(type(j))
@ -307,7 +307,7 @@ c For ialloy = 2, use single-element value (not average)
Arho2(nv2,j) = Arho2(nv2,j) + A2i*delij(m)*delij(n) Arho2(nv2,j) = Arho2(nv2,j) + A2i*delij(m)*delij(n)
nv2 = nv2+1 nv2 = nv2+1
do p = n,3 do p = n,3
Arho3(nv3,i) = Arho3(nv3,i) Arho3(nv3,i) = Arho3(nv3,i)
$ + A3j*delij(m)*delij(n)*delij(p) $ + A3j*delij(m)*delij(n)*delij(p)
Arho3(nv3,j) = Arho3(nv3,j) Arho3(nv3,j) = Arho3(nv3,j)
$ - A3i*delij(m)*delij(n)*delij(p) $ - A3i*delij(m)*delij(n)*delij(p)
@ -330,7 +330,7 @@ c Screening function
c Inputs: i = atom 1 id (integer) c Inputs: i = atom 1 id (integer)
c j = atom 2 id (integer) c j = atom 2 id (integer)
c rijsq = squared distance between i and j c rijsq = squared distance between i and j
c Outputs: sij = screening function c Outputs: sij = screening function
use meam_data use meam_data
implicit none implicit none
@ -485,7 +485,7 @@ c cutoff function
a = a*a a = a*a
a = a*a a = a*a
a = 1.d0-a a = 1.d0-a
fc = a*a fc = a*a
c fc = xi c fc = xi
endif endif
return return
@ -558,7 +558,7 @@ c dCikj2 = derivative of Cikj w.r.t. rjk
end end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

74
lib/meam/meam_force.F
View File

@ -1,13 +1,13 @@
c Extern "C" declaration has the form: c Extern "C" declaration has the form:
c c
c void meam_force_(int *, int *, int *, double *, int *, int *, int *, double *, c void meam_force_(int *, int *, int *, double *, int *, int *, int *, double *,
c int *, int *, int *, int *, double *, double *, c int *, int *, int *, int *, double *, double *,
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 *, double *,
c double *, double *, double *, double *, double *, double *, int *); c double *, double *, double *, double *, double *, double *, int *);
c c
c Call from pair_meam.cpp has the form: c Call from pair_meam.cpp has the form:
c c
c meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom, c meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom,
c &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], c &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],
@ -15,15 +15,15 @@ c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
c dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, c dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
c &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag); c &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag);
c c
subroutine meam_force(i, nmax, subroutine meam_force(i, nmax,
$ eflag_either, eflag_global, eflag_atom, vflag_atom, $ eflag_either, eflag_global, eflag_atom, vflag_atom,
$ eng_vdwl, eatom, ntype, type, fmap, x, $ eng_vdwl, eatom, ntype, type, fmap, x,
$ numneigh, firstneigh, numneigh_full, firstneigh_full, $ numneigh, firstneigh, numneigh_full, firstneigh_full,
$ scrfcn, dscrfcn, fcpair, $ scrfcn, dscrfcn, fcpair,
$ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp, $ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp,
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, f, $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, f,
$ vatom, errorflag) $ vatom, errorflag)
use meam_data use meam_data
@ -92,20 +92,20 @@ c
c Compute forces atom i c Compute forces atom i
elti = fmap(type(i)) elti = fmap(type(i))
if (elti.gt.0) then if (elti.gt.0) then
xitmp = x(1,i) xitmp = x(1,i)
yitmp = x(2,i) yitmp = x(2,i)
zitmp = x(3,i) zitmp = x(3,i)
c Treat each pair c Treat each pair
do jn = 1,numneigh do jn = 1,numneigh
j = firstneigh(jn) j = firstneigh(jn)
eltj = fmap(type(j)) eltj = fmap(type(j))
if (scrfcn(jn).ne.0.d0.and.eltj.gt.0) then if (scrfcn(jn).ne.0.d0.and.eltj.gt.0) then
sij = scrfcn(jn)*fcpair(jn) sij = scrfcn(jn)*fcpair(jn)
delij(1) = x(1,j) - xitmp delij(1) = x(1,j) - xitmp
delij(2) = x(2,j) - yitmp delij(2) = x(2,j) - yitmp
@ -115,7 +115,7 @@ c Treat each pair
if (rij2.lt.cutforcesq) then if (rij2.lt.cutforcesq) then
rij = sqrt(rij2) rij = sqrt(rij2)
r = rij r = rij
c Compute phi and phip c Compute phi and phip
ind = eltind(elti,eltj) ind = eltind(elti,eltj)
pp = rij*rdrar + 1.0D0 pp = rij*rdrar + 1.0D0
@ -123,12 +123,12 @@ c Compute phi and phip
kk = min(kk,nrar-1) kk = min(kk,nrar-1)
pp = pp - kk pp = pp - kk
pp = min(pp,1.0D0) pp = min(pp,1.0D0)
phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp
$ + phirar1(kk,ind))*pp + phirar(kk,ind) $ + phirar1(kk,ind))*pp + phirar(kk,ind)
phip = (phirar6(kk,ind)*pp + phirar5(kk,ind))*pp phip = (phirar6(kk,ind)*pp + phirar5(kk,ind))*pp
$ + phirar4(kk,ind) $ + phirar4(kk,ind)
recip = 1.0d0/r recip = 1.0d0/r
if (eflag_either.ne.0) then if (eflag_either.ne.0) then
if (eflag_global.ne.0) eng_vdwl = eng_vdwl + phi*sij if (eflag_global.ne.0) eng_vdwl = eng_vdwl + phi*sij
if (eflag_atom.ne.0) then if (eflag_atom.ne.0) then
@ -136,10 +136,10 @@ c Compute phi and phip
eatom(j) = eatom(j) + 0.5*phi*sij eatom(j) = eatom(j) + 0.5*phi*sij
endif endif
endif endif
c write(1,*) "force_meamf: phi: ",phi c write(1,*) "force_meamf: phi: ",phi
c write(1,*) "force_meamf: phip: ",phip c write(1,*) "force_meamf: phip: ",phip
c Compute pair densities and derivatives c Compute pair densities and derivatives
invrei = 1.d0/re_meam(elti,elti) invrei = 1.d0/re_meam(elti,elti)
ai = rij*invrei - 1.d0 ai = rij*invrei - 1.d0
@ -152,7 +152,7 @@ c Compute pair densities and derivatives
drhoa2i = -beta2_meam(elti)*invrei*rhoa2i drhoa2i = -beta2_meam(elti)*invrei*rhoa2i
rhoa3i = ro0i*fm_exp(-beta3_meam(elti)*ai) rhoa3i = ro0i*fm_exp(-beta3_meam(elti)*ai)
drhoa3i = -beta3_meam(elti)*invrei*rhoa3i drhoa3i = -beta3_meam(elti)*invrei*rhoa3i
if (elti.ne.eltj) then if (elti.ne.eltj) then
invrej = 1.d0/re_meam(eltj,eltj) invrej = 1.d0/re_meam(eltj,eltj)
aj = rij*invrej - 1.d0 aj = rij*invrej - 1.d0
@ -175,7 +175,7 @@ c Compute pair densities and derivatives
rhoa3j = rhoa3i rhoa3j = rhoa3i
drhoa3j = drhoa3i drhoa3j = drhoa3i
endif endif
if (ialloy.eq.1) then if (ialloy.eq.1) then
rhoa1j = rhoa1j * t1_meam(eltj) rhoa1j = rhoa1j * t1_meam(eltj)
rhoa2j = rhoa2j * t2_meam(eltj) rhoa2j = rhoa2j * t2_meam(eltj)
@ -217,15 +217,15 @@ c Compute pair densities and derivatives
arg1i1 = arg1i1 + Arho1(n,i)*delij(n) arg1i1 = arg1i1 + Arho1(n,i)*delij(n)
arg1j1 = arg1j1 - Arho1(n,j)*delij(n) arg1j1 = arg1j1 - Arho1(n,j)*delij(n)
arg3i3 = arg3i3 + Arho3b(n,i)*delij(n) arg3i3 = arg3i3 + Arho3b(n,i)*delij(n)
arg3j3 = arg3j3 - Arho3b(n,j)*delij(n) arg3j3 = arg3j3 - Arho3b(n,j)*delij(n)
enddo enddo
c rho0 terms c rho0 terms
drho0dr1 = drhoa0j * sij drho0dr1 = drhoa0j * sij
drho0dr2 = drhoa0i * sij drho0dr2 = drhoa0i * sij
c rho1 terms c rho1 terms
a1 = 2*sij/rij a1 = 2*sij/rij
drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1 drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1
drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1 drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1
a1 = 2.d0*sij/rij a1 = 2.d0*sij/rij
@ -233,10 +233,10 @@ c rho1 terms
drho1drm1(m) = a1*rhoa1j*Arho1(m,i) drho1drm1(m) = a1*rhoa1j*Arho1(m,i)
drho1drm2(m) = -a1*rhoa1i*Arho1(m,j) drho1drm2(m) = -a1*rhoa1i*Arho1(m,j)
enddo enddo
c rho2 terms c rho2 terms
a2 = 2*sij/rij2 a2 = 2*sij/rij2
drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2 drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2
$ - 2.d0/3.d0*Arho2b(i)*drhoa2j*sij $ - 2.d0/3.d0*Arho2b(i)*drhoa2j*sij
drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2 drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2
$ - 2.d0/3.d0*Arho2b(j)*drhoa2i*sij $ - 2.d0/3.d0*Arho2b(j)*drhoa2i*sij
@ -245,7 +245,7 @@ c rho2 terms
drho2drm1(m) = 0.d0 drho2drm1(m) = 0.d0
drho2drm2(m) = 0.d0 drho2drm2(m) = 0.d0
do n = 1,3 do n = 1,3
drho2drm1(m) = drho2drm1(m) drho2drm1(m) = drho2drm1(m)
$ + Arho2(vind2D(m,n),i)*delij(n) $ + Arho2(vind2D(m,n),i)*delij(n)
drho2drm2(m) = drho2drm2(m) drho2drm2(m) = drho2drm2(m)
$ - Arho2(vind2D(m,n),j)*delij(n) $ - Arho2(vind2D(m,n),j)*delij(n)
@ -253,14 +253,14 @@ c rho2 terms
drho2drm1(m) = a2*rhoa2j*drho2drm1(m) drho2drm1(m) = a2*rhoa2j*drho2drm1(m)
drho2drm2(m) = -a2*rhoa2i*drho2drm2(m) drho2drm2(m) = -a2*rhoa2i*drho2drm2(m)
enddo enddo
c rho3 terms c rho3 terms
rij3 = rij*rij2 rij3 = rij*rij2
a3 = 2*sij/rij3 a3 = 2*sij/rij3
a3a = 6.d0/5.d0*sij/rij a3a = 6.d0/5.d0*sij/rij
drho3dr1 = a3*(drhoa3j - 3*rhoa3j/rij)*arg1i3 drho3dr1 = a3*(drhoa3j - 3*rhoa3j/rij)*arg1i3
$ - a3a*(drhoa3j - rhoa3j/rij)*arg3i3 $ - a3a*(drhoa3j - rhoa3j/rij)*arg3i3
drho3dr2 = a3*(drhoa3i - 3*rhoa3i/rij)*arg1j3 drho3dr2 = a3*(drhoa3i - 3*rhoa3i/rij)*arg1j3
$ - a3a*(drhoa3i - rhoa3i/rij)*arg3j3 $ - a3a*(drhoa3i - rhoa3i/rij)*arg3j3
a3 = 6*sij/rij3 a3 = 6*sij/rij3
a3a = 6*sij/(5*rij) a3a = 6*sij/(5*rij)
@ -327,7 +327,7 @@ c Compute derivatives of weighting functions t wrt rij
dt3dr2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2) dt3dr2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2)
else if (ialloy.eq.2) then else if (ialloy.eq.2) then
dt1dr1 = 0.d0 dt1dr1 = 0.d0
dt1dr2 = 0.d0 dt1dr2 = 0.d0
dt2dr1 = 0.d0 dt2dr1 = 0.d0
@ -352,7 +352,7 @@ c Compute derivatives of weighting functions t wrt rij
dt2dr2 = aj*(t2_meam(elti)-t2j) dt2dr2 = aj*(t2_meam(elti)-t2j)
dt3dr1 = ai*(t3_meam(eltj)-t3i) dt3dr1 = ai*(t3_meam(eltj)-t3i)
dt3dr2 = aj*(t3_meam(elti)-t3j) dt3dr2 = aj*(t3_meam(elti)-t3j)
endif endif
c Compute derivatives of total density wrt rij, sij and rij(3) c Compute derivatives of total density wrt rij, sij and rij(3)
@ -393,9 +393,9 @@ c Compute derivatives wrt sij, but only if necessary
drho1ds1 = a1*rhoa1j*arg1i1 drho1ds1 = a1*rhoa1j*arg1i1
drho1ds2 = a1*rhoa1i*arg1j1 drho1ds2 = a1*rhoa1i*arg1j1
a2 = 2.d0/rij2 a2 = 2.d0/rij2
drho2ds1 = a2*rhoa2j*arg1i2 drho2ds1 = a2*rhoa2j*arg1i2
$ - 2.d0/3.d0*Arho2b(i)*rhoa2j $ - 2.d0/3.d0*Arho2b(i)*rhoa2j
drho2ds2 = a2*rhoa2i*arg1j2 drho2ds2 = a2*rhoa2i*arg1j2
$ - 2.d0/3.d0*Arho2b(j)*rhoa2i $ - 2.d0/3.d0*Arho2b(j)*rhoa2i
a3 = 2.d0/rij3 a3 = 2.d0/rij3
a3a = 6.d0/(5.d0*rij) a3a = 6.d0/(5.d0*rij)
@ -403,7 +403,7 @@ c Compute derivatives wrt sij, but only if necessary
drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3 drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3
if (ialloy.eq.1) then if (ialloy.eq.1) then
a1i = 0.d0 a1i = 0.d0
a1j = 0.d0 a1j = 0.d0
a2i = 0.d0 a2i = 0.d0
@ -455,7 +455,7 @@ c Compute derivatives wrt sij, but only if necessary
if( rho0(j) .ne. 0.d0 ) then if( rho0(j) .ne. 0.d0 ) then
aj = rhoa0i/rho0(j) aj = rhoa0i/rho0(j)
end if end if
dt1ds1 = ai*(t1_meam(eltj)-t1i) dt1ds1 = ai*(t1_meam(eltj)-t1i)
dt1ds2 = aj*(t1_meam(elti)-t1j) dt1ds2 = aj*(t1_meam(elti)-t1j)
dt2ds1 = ai*(t2_meam(eltj)-t2i) dt2ds1 = ai*(t2_meam(eltj)-t2i)
@ -492,7 +492,7 @@ c Compute derivatives of energy wrt rij, sij and rij(3)
do m = 1,3 do m = 1,3
dUdrijm(m) = fp(i)*drhodrm1(m) + fp(j)*drhodrm2(m) dUdrijm(m) = fp(i)*drhodrm1(m) + fp(j)*drhodrm2(m)
enddo enddo
c Add the part of the force due to dUdrij and dUdsij c Add the part of the force due to dUdrij and dUdsij
force = dUdrij*recip + dUdsij*dscrfcn(jn) force = dUdrij*recip + dUdsij*dscrfcn(jn)
@ -603,6 +603,6 @@ c end of j loop
c else if elti=0, this is not a meam atom c else if elti=0, this is not a meam atom
endif endif
return return
end end

130
lib/meam/meam_setup_done.F
View File

@ -121,7 +121,7 @@ c If i<j and term is unset, use default values (e.g. mean of i-i and j-j)
$ - delta_meam(i,j) $ - delta_meam(i,j)
endif endif
else else
Ec_meam(i,j) = (Ec_meam(i,i)+Ec_meam(j,j))/2.d0 Ec_meam(i,j) = (Ec_meam(i,i)+Ec_meam(j,j))/2.d0
$ - delta_meam(i,j) $ - delta_meam(i,j)
endif endif
endif endif
@ -147,7 +147,7 @@ c where i>j, set equal to the i<j element. Likewise for Cmax.
enddo enddo
c ebound gives the squared distance such that, for rik2 or rjk2>ebound, c ebound gives the squared distance such that, for rik2 or rjk2>ebound,
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 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 c triplets (i,j,k) so that ebound(i,j) is the maximized over k
do i = 1,neltypes do i = 1,neltypes
@ -237,7 +237,7 @@ c phi_aa
$ (-Z2*scrn/Z1)**n * phi_meam(rarat*arat**n,a,a) $ (-Z2*scrn/Z1)**n * phi_meam(rarat*arat**n,a,a)
enddo enddo
endif endif
c phi_bb c phi_bb
phibb = phi_meam(rarat,b,b) phibb = phi_meam(rarat,b,b)
call get_Zij(Z1,lattce_meam(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) $ (-Z2*scrn/Z1)**n * phi_meam(rarat*arat**n,b,b)
enddo enddo
endif endif
if (lattce_meam(a,b).eq.'b1'. if (lattce_meam(a,b).eq.'b1'.
$ or.lattce_meam(a,b).eq.'b2') then $ or.lattce_meam(a,b).eq.'b2') then
c Add contributions to the B1 or B2 potential c Add contributions to the B1 or B2 potential
@ -276,11 +276,11 @@ c atoms.
call get_sijk(C,b,b,a,s221) call get_sijk(C,b,b,a,s221)
S11 = s111 * s111 * s112 * s112 S11 = s111 * s111 * s112 * s112
S22 = s221**4 S22 = s221**4
phir(j,nv2) = phir(j,nv2) - phir(j,nv2) = phir(j,nv2) -
$ 0.75*S11*phiaa - 0.25*S22*phibb $ 0.75*S11*phiaa - 0.25*S22*phibb
endif endif
else else
nmax = 10 nmax = 10
do n = 1,nmax do n = 1,nmax
@ -307,7 +307,7 @@ c endif
phir(j,nv2) = frac*phir(j,nv2) + (1-frac)*phizbl phir(j,nv2) = frac*phir(j,nv2) + (1-frac)*phizbl
endif endif
endif endif
enddo enddo
c call interpolation c call interpolation
@ -320,14 +320,14 @@ c call interpolation
end end
c----------------------------------------------------------------------c c----------------------------------------------------------------------c
c Compute MEAM pair potential for distance r, element types a and b c Compute MEAM pair potential for distance r, element types a and b
c c
real*8 recursive function phi_meam(r,a,b)result(phi_m) real*8 recursive function phi_meam(r,a,b)result(phi_m)
use meam_data use meam_data
implicit none implicit none
integer a,b integer a,b
real*8 r real*8 r
real*8 a1,a2,a12 real*8 a1,a2,a12
@ -375,15 +375,15 @@ c calculate average weighting factors for the reference structure
t31av = t3_meam(a) t31av = t3_meam(a)
t32av = t3_meam(b) t32av = t3_meam(b)
else else
scalfac = 1.0/(rho01+rho02) scalfac = 1.0/(rho01+rho02)
t11av = scalfac*(t1_meam(a)*rho01 + t1_meam(b)*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) t21av = scalfac*(t2_meam(a)*rho01 + t2_meam(b)*rho02)
t22av = t21av t22av = t21av
t31av = scalfac*(t3_meam(a)*rho01 + t3_meam(b)*rho02) t31av = scalfac*(t3_meam(a)*rho01 + t3_meam(b)*rho02)
t32av = t31av t32av = t31av
endif endif
else else
c average weighting factors for the reference structure, eqn. I.8 c average weighting factors for the reference structure, eqn. I.8
call get_tavref(t11av,t21av,t31av,t12av,t22av,t32av, call get_tavref(t11av,t21av,t31av,t12av,t22av,t32av,
$ t1_meam(a),t2_meam(a),t3_meam(a), $ 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 = (Z12*rho02)**2 + 2.0/3.0*t22av*rho22**2
rhobar1 = sqrt(rhobar1) rhobar1 = sqrt(rhobar1)
endif endif
else else
c for other structures, use formalism developed in Huang's paper c for other structures, use formalism developed in Huang's paper
c 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 If using mixing rule for t, apply to reference structure; else
c use precomputed values c use precomputed values
if (mix_ref_t.eq.1) then if (mix_ref_t.eq.1) then
@ -508,8 +508,8 @@ c account for second neighbor a-a potential here...
enddo enddo
endif endif
phi_m = Eu/3. - F1/4. - F2/12. - phiaa phi_m = Eu/3. - F1/4. - F2/12. - phiaa
else else
c c
c potential is computed from Rose function and embedding energy c potential is computed from Rose function and embedding energy
phi_m = (2*Eu - F1 - F2)/Z12 phi_m = (2*Eu - F1 - F2)/Z12
c c
@ -519,7 +519,7 @@ c if r = 0, just return 0
if (r.eq.0.d0) then if (r.eq.0.d0) then
phi_m = 0.d0 phi_m = 0.d0
endif endif
return return
end end
@ -546,7 +546,7 @@ c loop over element types
call G_gam(gam,ibar_meam(a),gsmooth_factor, call G_gam(gam,ibar_meam(a),gsmooth_factor,
$ Gbar,errorflag) $ Gbar,errorflag)
endif endif
c The zeroth order density in the reference structure, with c The zeroth order density in the reference structure, with
c equilibrium spacing, is just the number of first neighbors times c equilibrium spacing, is just the number of first neighbors times
c the rho0_meam coefficient... c the rho0_meam coefficient...
@ -569,7 +569,7 @@ c screening)
return return
end end
c----------------------------------------------------------------------c c----------------------------------------------------------------------c
c Shape factors for various configurations c Shape factors for various configurations
subroutine get_shpfcn(s,latt) subroutine get_shpfcn(s,latt)
implicit none implicit none
@ -599,7 +599,7 @@ c call error('Lattice not defined in get_shpfcn.')
endif endif
return return
end end
c------------------------------------------------------------------------------c c------------------------------------------------------------------------------c
c Average weighting factors for the reference structure c Average weighting factors for the reference structure
subroutine get_tavref(t11av,t21av,t31av,t12av,t22av,t32av, subroutine get_tavref(t11av,t21av,t31av,t12av,t22av,t32av,
$ t11,t21,t31,t12,t22,t32, $ t11,t21,t31,t12,t22,t32,
@ -651,7 +651,7 @@ c call error('Lattice not defined in get_tavref.')
endif endif
return return
end end
c------------------------------------------------------------------------------c c------------------------------------------------------------------------------c
c Number of neighbors for the reference structure c Number of neighbors for the reference structure
subroutine get_Zij(Zij,latt) subroutine get_Zij(Zij,latt)
implicit none implicit none
@ -681,8 +681,8 @@ c call error('Lattice not defined in get_Zij.')
return return
end end
c------------------------------------------------------------------------------c c------------------------------------------------------------------------------c
c Zij2 = number of second neighbors, a = distance ratio R1/R2, and S = second c Zij2 = number of second neighbors, a = distance ratio R1/R2, and S = second
c neighbor screening function for lattice type "latt" c neighbor screening function for lattice type "latt"
subroutine get_Zij2(Zij2,a,S,latt,cmin,cmax) 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.') c call error('Lattice not defined in get_Zij2.')
endif endif
c Compute screening for each first neighbor c Compute screening for each first neighbor
C = 4.d0/(a*a) - 1.d0 C = 4.d0/(a*a) - 1.d0
x = (C-cmin)/(cmax-cmin) x = (C-cmin)/(cmax-cmin)
call fcut(x,sijk) call fcut(x,sijk)
@ -744,8 +744,8 @@ c There are numscr first neighbors screening the second neighbors
return return
end end
c------------------------------------------------------------------------------c c------------------------------------------------------------------------------c
subroutine get_sijk(C,i,j,k,sijk) subroutine get_sijk(C,i,j,k,sijk)
use meam_data use meam_data
implicit none implicit none
@ -757,7 +757,7 @@ c------------------------------------------------------------------------------c
return return
end end
c------------------------------------------------------------------------------c c------------------------------------------------------------------------------c
c Calculate density functions, assuming reference configuration c Calculate density functions, assuming reference configuration
subroutine get_densref(r,a,b,rho01,rho11,rho21,rho31, subroutine get_densref(r,a,b,rho01,rho11,rho21,rho31,
$ rho02,rho12,rho22,rho32) $ 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) rhoa32 = rho0_meam(b)*fm_exp(-beta3_meam(b)*a2)
lat = lattce_meam(a,b) lat = lattce_meam(a,b)
rho11 = 0.d0 rho11 = 0.d0
rho21 = 0.d0 rho21 = 0.d0
rho31 = 0.d0 rho31 = 0.d0
@ -807,13 +807,13 @@ c Calculate density functions, assuming reference configuration
rho01 = 6*rhoa02 rho01 = 6*rhoa02
rho02 = 6*rhoa01 rho02 = 6*rhoa01
else if (lat.eq.'dia') then else if (lat.eq.'dia') then
rho01 = 4*rhoa02 rho01 = 4*rhoa02
rho02 = 4*rhoa01 rho02 = 4*rhoa01
rho31 = 32.d0/9.d0*rhoa32*rhoa32 rho31 = 32.d0/9.d0*rhoa32*rhoa32
rho32 = 32.d0/9.d0*rhoa31*rhoa31 rho32 = 32.d0/9.d0*rhoa31*rhoa31
else if (lat.eq.'hcp') then else if (lat.eq.'hcp') then
rho01 = 12*rhoa02 rho01 = 12*rhoa02
rho02 = 12*rhoa01 rho02 = 12*rhoa01
rho31 = 1.d0/3.d0*rhoa32*rhoa32 rho31 = 1.d0/3.d0*rhoa32*rhoa32
rho32 = 1.d0/3.d0*rhoa31*rhoa31 rho32 = 1.d0/3.d0*rhoa31*rhoa31
else if (lat.eq.'dim') then else if (lat.eq.'dim') then
@ -830,7 +830,7 @@ c Calculate density functions, assuming reference configuration
rho01 = rhoa01 rho01 = rhoa01
rho02 = rhoa02 rho02 = rhoa02
rho11 = rhoa11 rho11 = rhoa11
rho12 = rhoa12 rho12 = rhoa12
rho21 = rhoa21 rho21 = rhoa21
rho22 = rhoa22 rho22 = rhoa22
rho31 = rhoa31 rho31 = rhoa31
@ -858,10 +858,10 @@ c call error('Lattice not defined in get_densref.')
call get_Zij2(Zij2nn,arat,scrn,lat, call get_Zij2(Zij2nn,arat,scrn,lat,
$ Cmin_meam(a,a,b),Cmax_meam(a,a,b)) $ Cmin_meam(a,a,b),Cmax_meam(a,a,b))
a1 = arat*r/re_meam(a,a) - 1.d0 a1 = arat*r/re_meam(a,a) - 1.d0
a2 = arat*r/re_meam(b,b) - 1.d0 a2 = arat*r/re_meam(b,b) - 1.d0
rhoa01nn = rho0_meam(a)*fm_exp(-beta0_meam(a)*a1) rhoa01nn = rho0_meam(a)*fm_exp(-beta0_meam(a)*a1)
rhoa02nn = rho0_meam(b)*fm_exp(-beta0_meam(b)*a2) 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 return
end end
c--------------------------------------------------------------------- c---------------------------------------------------------------------
c Compute ZBL potential c Compute ZBL potential
c c
real*8 function zbl(r,z1,z2) real*8 function zbl(r,z1,z2)
use meam_data , only : fm_exp use meam_data , only : fm_exp
implicit none implicit none
integer i,z1,z2 integer i,z1,z2
real*8 r,c,d,a,azero,cc,x real*8 r,c,d,a,azero,cc,x
dimension c(4),d(4) dimension c(4),d(4)
data c /0.028171,0.28022,0.50986,0.18175/ data c /0.028171,0.28022,0.50986,0.18175/
data d /0.20162,0.40290,0.94229,3.1998/ data d /0.20162,0.40290,0.94229,3.1998/
data azero /0.4685/ data azero /0.4685/
data cc /14.3997/ data cc /14.3997/
c azero = (9pi^2/128)^1/3 (0.529) Angstroms c azero = (9pi^2/128)^1/3 (0.529) Angstroms
a = azero/(z1**0.23+z2**0.23) a = azero/(z1**0.23+z2**0.23)
zbl = 0.0 zbl = 0.0
x = r/a x = r/a
do i=1,4 do i=1,4
zbl = zbl + c(i)*fm_exp(-d(i)*x) zbl = zbl + c(i)*fm_exp(-d(i)*x)
enddo enddo
if (r.gt.0.d0) zbl = zbl*z1*z2/r*cc if (r.gt.0.d0) zbl = zbl*z1*z2/r*cc
return return
end end
c--------------------------------------------------------------------- c---------------------------------------------------------------------
c Compute Rose energy function, I.16 c Compute Rose energy function, I.16
c c
real*8 function erose(r,re,alpha,Ec,repuls,attrac,form) 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 implicit none
real*8 r,re,alpha,Ec,repuls,attrac,astar,a3 real*8 r,re,alpha,Ec,repuls,attrac,astar,a3
integer form integer form
erose = 0.d0 erose = 0.d0
if (r.gt.0.d0) then if (r.gt.0.d0) then
astar = alpha * (r/re - 1.d0) astar = alpha * (r/re - 1.d0)
a3 = 0.d0 a3 = 0.d0
@ -950,7 +950,7 @@ c
return return
end end
c ----------------------------------------------------------------------- c -----------------------------------------------------------------------
subroutine interpolate_meam(ind) subroutine interpolate_meam(ind)
@ -980,7 +980,7 @@ c phir interp
phirar1(j,ind) = ((phirar(j-2,ind)-phirar(j+2,ind)) + phirar1(j,ind) = ((phirar(j-2,ind)-phirar(j+2,ind)) +
$ 8.0D0*(phirar(j+1,ind)-phirar(j-1,ind)))/12. $ 8.0D0*(phirar(j+1,ind)-phirar(j-1,ind)))/12.
enddo enddo
do j = 1,nrar-1 do j = 1,nrar-1
phirar2(j,ind) = 3.0D0*(phirar(j+1,ind)-phirar(j,ind)) - phirar2(j,ind) = 3.0D0*(phirar(j+1,ind)-phirar(j,ind)) -
$ 2.0D0*phirar1(j,ind) - phirar1(j+1,ind) $ 2.0D0*phirar1(j,ind) - phirar1(j+1,ind)
@ -989,7 +989,7 @@ c phir interp
enddo enddo
phirar2(nrar,ind) = 0.0D0 phirar2(nrar,ind) = 0.0D0
phirar3(nrar,ind) = 0.0D0 phirar3(nrar,ind) = 0.0D0
do j = 1,nrar do j = 1,nrar
phirar4(j,ind) = phirar1(j,ind)/drar phirar4(j,ind) = phirar1(j,ind)/drar
phirar5(j,ind) = 2.0D0*phirar2(j,ind)/drar phirar5(j,ind) = 2.0D0*phirar2(j,ind)/drar
@ -1014,7 +1014,7 @@ c
kk = min(kk,nrar-1) kk = min(kk,nrar-1)
pp = pp - kk pp = pp - kk
pp = min(pp,1.0D0) 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) $ + phirar1(kk,ind))*pp + phirar(kk,ind)
return return

2
lib/meam/meam_setup_global.F
View File

@ -14,7 +14,7 @@ c
c c
subroutine meam_setup_global(nelt, lat, z, ielement, atwt, alpha, 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) $ t0, t1, t2, t3, rozero, ibar)
use meam_data use meam_data

22
lib/meam/meam_setup_param.F
View File

@ -1,5 +1,5 @@
c 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) subroutine meam_checkindex(num,lim,nidx,idx,ierr)
implicit none implicit none
integer i,num,lim,nidx,idx(3),ierr integer i,num,lim,nidx,idx(3),ierr
@ -18,20 +18,20 @@ c do a sanity check on index parameters
enddo enddo
end end
c c
c Declaration in pair_meam.h: c Declaration in pair_meam.h:
c c
c void meam_setup_param(int *, double *, int *, int *, int *); c void meam_setup_param(int *, double *, int *, int *, int *);
c c
c Call in pair_meam.cpp c Call in pair_meam.cpp
c c
c meam_setup_param(&which,&value,&nindex,index,&errorflag); 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 The "which" argument corresponds to the index of the "keyword" array
c in pair_meam.cpp: c in pair_meam.cpp:
c c
c 0 = Ec_meam c 0 = Ec_meam
c 1 = alpha_meam c 1 = alpha_meam
c 2 = rho0_meam c 2 = rho0_meam
@ -54,7 +54,7 @@ c 18 = zbl_meam
c 19 = emb_lin_neg c 19 = emb_lin_neg
c 20 = bkgd_dyn c 20 = bkgd_dyn
subroutine meam_setup_param(which, value, nindex, subroutine meam_setup_param(which, value, nindex,
$ index, errorflag) $ index, errorflag)
use meam_data use meam_data
@ -77,7 +77,7 @@ c 1 = alpha_meam
call meam_checkindex(2,maxelt,nindex,index,errorflag) call meam_checkindex(2,maxelt,nindex,index,errorflag)
if (errorflag.ne.0) return if (errorflag.ne.0) return
alpha_meam(index(1),index(2)) = value alpha_meam(index(1),index(2)) = value
c 2 = rho0_meam c 2 = rho0_meam
else if (which.eq.2) then else if (which.eq.2) then
call meam_checkindex(1,maxelt,nindex,index,errorflag) call meam_checkindex(1,maxelt,nindex,index,errorflag)