git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@265 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
51
lib/meam/Makefile
Executable file
51
lib/meam/Makefile
Executable file
@ -0,0 +1,51 @@
|
||||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = gfortran
|
||||
F90FLAGS = -g -fno-second-underscore
|
||||
#F90FLAGS = -O
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
LINK = g++
|
||||
LINKFLAGS = -O
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
48
lib/meam/Makefile.g95
Normal file
48
lib/meam/Makefile.g95
Normal file
@ -0,0 +1,48 @@
|
||||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = g95
|
||||
F90FLAGS = -O
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
51
lib/meam/Makefile.gfortran
Normal file
51
lib/meam/Makefile.gfortran
Normal file
@ -0,0 +1,51 @@
|
||||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = gfortran
|
||||
F90FLAGS = -g -fno-second-underscore
|
||||
#F90FLAGS = -O
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
LINK = g++
|
||||
LINKFLAGS = -O
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
48
lib/meam/Makefile.ifort
Normal file
48
lib/meam/Makefile.ifort
Normal file
@ -0,0 +1,48 @@
|
||||
# *
|
||||
# *_________________________________________________________________________*
|
||||
# * MEAM: MODEFIED EMBEDDED ATOM METHOD *
|
||||
# * DESCRIPTION: SEE READ-ME *
|
||||
# * FILE NAME: Makefile *
|
||||
# * AUTHORS: Greg Wagner, Sandia National Laboratories *
|
||||
# * CONTACT: gjwagne@sandia.gov *
|
||||
# *_________________________________________________________________________*/
|
||||
|
||||
SHELL = /bin/sh
|
||||
|
||||
# ------ FILES ------
|
||||
|
||||
SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F
|
||||
|
||||
FILES = $(SRC) Makefile
|
||||
|
||||
# ------ DEFINITIONS ------
|
||||
|
||||
LIB = libmeam.a
|
||||
OBJ = $(SRC:.F=.o)
|
||||
|
||||
# ------ SETTINGS ------
|
||||
|
||||
F90 = ifort
|
||||
F90FLAGS = -O
|
||||
ARCHIVE = ar
|
||||
ARCHFLAG = -rc
|
||||
USRLIB =
|
||||
SYSLIB =
|
||||
|
||||
# ------ MAKE PROCEDURE ------
|
||||
|
||||
lib: $(OBJ)
|
||||
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
|
||||
|
||||
# ------ COMPILE RULES ------
|
||||
|
||||
%.o:%.F
|
||||
$(F90) $(F90FLAGS) -c $<
|
||||
|
||||
# ------ CLEAN ------
|
||||
|
||||
clean:
|
||||
-rm *.o *.mod *~ $(LIB)
|
||||
|
||||
tar:
|
||||
-tar -cvf ../MEAM.tar $(FILES)
|
||||
20
lib/meam/README
Normal file
20
lib/meam/README
Normal file
@ -0,0 +1,20 @@
|
||||
MEAM (modified embedded atom method) library
|
||||
|
||||
Greg Wagner, Sandia National Labs
|
||||
gjwagne at sandia.gov
|
||||
Jan 2007
|
||||
|
||||
--------------
|
||||
|
||||
This library is in implementation of the MEAM potential, specifically
|
||||
designed to work with LAMMPS.
|
||||
|
||||
This library must be built with a F90 compiler, before LAMMPS is
|
||||
built, so LAMMPS can link against it.
|
||||
|
||||
Build the library using one of the provided Makefiles or creating your
|
||||
own, specific to your compiler and system. For example:
|
||||
|
||||
make -f Makefile.g95
|
||||
|
||||
If the build is successful, you should end up with a libmeam.a file.
|
||||
76
lib/meam/meam_data.F
Executable file
76
lib/meam/meam_data.F
Executable file
@ -0,0 +1,76 @@
|
||||
|
||||
module meam_data
|
||||
|
||||
integer, parameter :: maxelt = 5
|
||||
|
||||
c cutforce = force cutoff
|
||||
c cutforcesq = force cutoff squared
|
||||
|
||||
real*8 cutforce,cutforcesq
|
||||
|
||||
c Ec_meam = cohesive energy
|
||||
c re_meam = nearest-neighbor distance
|
||||
c Omega_meam = atomic volume
|
||||
c B_meam = bulk modulus
|
||||
c Z_meam = number of first neighbors for reference structure
|
||||
c ielt_meam = atomic number of element
|
||||
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 t[0-3]_meam = coefficients on densities in Gamma computation
|
||||
c ibar_meam(i) = selection parameter for Gamma function for elt i,
|
||||
c lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j)
|
||||
c neltypes = maximum number of element type defined
|
||||
c eltind = index number of pair (similar to Voigt notation; ij = ji)
|
||||
c phir = pair potential function array
|
||||
c phirar[1-6] = spline coeffs
|
||||
c attrac_meam = attraction parameter in Rose energy
|
||||
c repuls_meam = repulsion parameter in Rose energy
|
||||
c nn2_meam = 1 if second nearest neighbors are to be computed, else 0
|
||||
c Cmin_meam, Cmax_meam = min and max values in screening cutoff
|
||||
c rc_meam = cutoff distance for meam
|
||||
c delr_meam = cutoff region for meam
|
||||
c ebound_meam = factor giving maximum boundary of sceen fcn ellipse
|
||||
c augt1 = flag for whether t1 coefficient should be augmented
|
||||
c gsmooth_factor = factor determining length of G smoothing region
|
||||
c vind[23]D = Voight notation index maps for 2 and 3D
|
||||
c v2D,v3D = array of factors to apply for Voight notation
|
||||
|
||||
c nr,dr = pair function discretization parameters
|
||||
c nrar,rdrar = spline coeff array parameters
|
||||
|
||||
real*8 Ec_meam(maxelt,maxelt),re_meam(maxelt,maxelt)
|
||||
real*8 Omega_meam(maxelt),Z_meam(maxelt)
|
||||
real*8 A_meam(maxelt),alpha_meam(maxelt,maxelt),rho0_meam(maxelt)
|
||||
real*8 delta_meam(maxelt,maxelt)
|
||||
real*8 beta0_meam(maxelt),beta1_meam(maxelt)
|
||||
real*8 beta2_meam(maxelt),beta3_meam(maxelt)
|
||||
real*8 t0_meam(maxelt),t1_meam(maxelt)
|
||||
real*8 t2_meam(maxelt),t3_meam(maxelt)
|
||||
integer ibar_meam(maxelt),ielt_meam(maxelt)
|
||||
character*3 lattce_meam(maxelt,maxelt)
|
||||
integer nn2_meam(maxelt,maxelt)
|
||||
integer eltind(maxelt,maxelt)
|
||||
integer neltypes
|
||||
|
||||
real*8, allocatable :: phir(:,:)
|
||||
|
||||
real*8, allocatable :: phirar(:,:),phirar1(:,:),phirar2(:,:),
|
||||
$ 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)
|
||||
integer augt1
|
||||
real*8 gsmooth_factor
|
||||
integer vind2D(3,3),vind3D(3,3,3)
|
||||
integer v2D(6),v3D(10)
|
||||
|
||||
integer nr,nrar
|
||||
real*8 dr,rdrar
|
||||
|
||||
end module
|
||||
201
lib/meam/meam_dens_final.F
Executable file
201
lib/meam/meam_dens_final.F
Executable file
@ -0,0 +1,201 @@
|
||||
c Extern "C" declaration has the form:
|
||||
c
|
||||
c void meam_dens_final_(int *, int *, int *, double *, int *, int *, int *,
|
||||
c double *, double *, double *, double *, double *, double *,
|
||||
c double *, double *, double *, double *, double *, double *,
|
||||
c double *, double *, double *, double *, int *);
|
||||
c
|
||||
c Call from pair_meam.cpp has the form:
|
||||
c
|
||||
c meam_dens_final_(&nlocal,&nmax,&eflag,&eng_vdwl,ntype,type,fmap,
|
||||
c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
|
||||
c &arho3b[0][0],&t_ave[0][0],gamma,dgamma1,
|
||||
c dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
|
||||
c
|
||||
|
||||
subroutine meam_dens_final(nlocal, nmax, eflag, eng_vdwl,
|
||||
$ ntype, type, fmap,
|
||||
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave,
|
||||
$ Gamma, dGamma1, dGamma2, dGamma3,
|
||||
$ rho, rho0, rho1, rho2, rho3, fp, errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer nlocal, nmax, eflag, ntype, type, fmap
|
||||
real*8 eng_vdwl, Arho1, Arho2
|
||||
real*8 Arho2b, Arho3, Arho3b
|
||||
real*8 t_ave
|
||||
real*8 Gamma, dGamma1, dGamma2, dGamma3
|
||||
real*8 rho, rho0, rho1, rho2, rho3
|
||||
real*8 fp
|
||||
integer errorflag
|
||||
|
||||
dimension type(nmax), fmap(ntype)
|
||||
dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax)
|
||||
dimension Arho3(10,nmax), Arho3b(3,nmax), t_ave(3,nmax)
|
||||
dimension Gamma(nmax), dGamma1(nmax), dGamma2(nmax)
|
||||
dimension dGamma3(nmax), rho(nmax), rho0(nmax)
|
||||
dimension rho1(nmax), rho2(nmax), rho3(nmax)
|
||||
dimension fp(nmax)
|
||||
|
||||
integer i, elti
|
||||
integer m
|
||||
real*8 rhob, G, dG, Gbar, dGbar, gam, shp(3), shpi(3), Z
|
||||
real*8 B, denom
|
||||
|
||||
c Complete the calculation of density
|
||||
|
||||
do i = 1,nlocal
|
||||
|
||||
elti = fmap(type(i))
|
||||
rho1(i) = 0.d0
|
||||
rho2(i) = -1.d0/3.d0*Arho2b(i)*Arho2b(i)
|
||||
rho3(i) = 0.d0
|
||||
do m = 1,3
|
||||
rho1(i) = rho1(i) + Arho1(m,i)*Arho1(m,i)
|
||||
rho3(i) = rho3(i) - 3.d0/5.d0*Arho3b(m,i)*Arho3b(m,i)
|
||||
enddo
|
||||
do m = 1,6
|
||||
rho2(i) = rho2(i) + v2D(m)*Arho2(m,i)*Arho2(m,i)
|
||||
enddo
|
||||
do m = 1,10
|
||||
rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i)
|
||||
enddo
|
||||
|
||||
if( rho0(i) .gt. 0.0 ) then
|
||||
t_ave(1,i) = t_ave(1,i)/rho0(i)
|
||||
t_ave(2,i) = t_ave(2,i)/rho0(i)
|
||||
t_ave(3,i) = t_ave(3,i)/rho0(i)
|
||||
endif
|
||||
|
||||
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
|
||||
|
||||
call G_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,errorflag)
|
||||
if (errorflag.ne.0) return
|
||||
if (ibar_meam(elti).le.0) then
|
||||
Gbar = 1.d0
|
||||
else
|
||||
call get_shpfcn(shp,lattce_meam(elti,elti))
|
||||
gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2)
|
||||
$ +t_ave(3,i)*shp(3))/(Z_meam(elti)**2)
|
||||
call G_gam(gam,ibar_meam(elti),gsmooth_factor,
|
||||
$ Gbar,errorflag)
|
||||
endif
|
||||
rho(i) = rho0(i) * G
|
||||
rhob = rho(i)/(rho0_meam(elti)*Z_meam(elti)*Gbar)
|
||||
|
||||
Z = Z_meam(elti)
|
||||
call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG)
|
||||
if (ibar_meam(elti).le.0) then
|
||||
Gbar = 1.d0
|
||||
dGbar = 0.d0
|
||||
else
|
||||
call get_shpfcn(shpi,lattce_meam(elti,elti))
|
||||
gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2)
|
||||
$ +t_ave(3,i)*shpi(3))/(Z*Z)
|
||||
call dG_gam(gam,ibar_meam(elti),gsmooth_factor,
|
||||
$ Gbar,dGbar)
|
||||
endif
|
||||
denom = 1.d0/(rho0_meam(elti)*Z*Gbar)
|
||||
dGamma1(i) = (G - 2*dG*Gamma(i))*denom
|
||||
|
||||
if( rho0(i) .ne. 0.d0 ) then
|
||||
dGamma2(i) = (dG/rho0(i))*denom
|
||||
else
|
||||
dGamma2(i) = 0.d0
|
||||
end if
|
||||
|
||||
dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom
|
||||
|
||||
B = A_meam(elti)*Ec_meam(elti,elti)
|
||||
|
||||
if( rhob .ne. 0.d0 ) then
|
||||
fp(i) = B*(log(rhob)+1.d0)
|
||||
if (eflag.eq.1) eng_vdwl = eng_vdwl + B*rhob*log(rhob)
|
||||
else
|
||||
fp(i) = B
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine G_gam(Gamma,ibar,gsmooth_factor,G,errorflag)
|
||||
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 3 => G = 2/(1+exp(-Gamma))
|
||||
c 4 => G = sqrt(1+Gamma)
|
||||
implicit none
|
||||
real*8 Gamma,G
|
||||
real*8 gsmooth_factor, gsmooth_switchpoint
|
||||
integer ibar, errorflag
|
||||
if (ibar.eq.0.or.ibar.eq.4) then
|
||||
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1)
|
||||
if (Gamma.lt.gsmooth_switchpoint) then
|
||||
c e.g. gsmooth_factor is 99, then:
|
||||
c gsmooth_switchpoint = -0.99
|
||||
c G = 0.01*(-0.99/Gamma)**99
|
||||
G = 1/(gsmooth_factor+1)
|
||||
$ *(gsmooth_switchpoint/Gamma)**gsmooth_factor
|
||||
G = sqrt(G)
|
||||
else
|
||||
G = sqrt(1.d0+Gamma)
|
||||
endif
|
||||
else if (ibar.eq.1) then
|
||||
G = exp(Gamma/2.d0)
|
||||
else if (ibar.eq.3) then
|
||||
G = 2.d0/(1.d0+exp(-Gamma))
|
||||
else
|
||||
errorflag = 1
|
||||
endif
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dG_gam(Gamma,ibar,gsmooth_factor,G,dG)
|
||||
c Compute G(Gamma) and dG(gamma) based on selection flag ibar:
|
||||
c 0 => G = sqrt(1+Gamma)
|
||||
c 1 => G = exp(Gamma/2)
|
||||
c 2 => not implemented
|
||||
c 3 => G = 2/(1+exp(-Gamma))
|
||||
c 4 => G = sqrt(1+Gamma)
|
||||
real*8 Gamma,G,dG
|
||||
real*8 gsmooth_factor, gsmooth_switchpoint
|
||||
integer ibar
|
||||
if (ibar.eq.0.or.ibar.eq.4) then
|
||||
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1)
|
||||
if (Gamma.lt.gsmooth_switchpoint) then
|
||||
c e.g. gsmooth_factor is 99, then:
|
||||
c gsmooth_switchpoint = -0.99
|
||||
c G = 0.01*(-0.99/Gamma)**99
|
||||
G = 1/(gsmooth_factor+1)
|
||||
$ *(gsmooth_switchpoint/Gamma)**gsmooth_factor
|
||||
G = sqrt(G)
|
||||
dG = -gsmooth_factor*G/(2.0*Gamma)
|
||||
else
|
||||
G = sqrt(1.d0+Gamma)
|
||||
dG = 1.d0/(2.d0*G)
|
||||
endif
|
||||
else if (ibar.eq.1) then
|
||||
G = exp(Gamma/2.d0)
|
||||
dG = G/2.d0
|
||||
else if (ibar.eq.3) then
|
||||
G = 2.d0/(1.d0+exp(-Gamma))
|
||||
dG = G*(2.d0-G)/2
|
||||
endif
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
521
lib/meam/meam_dens_init.F
Executable file
521
lib/meam/meam_dens_init.F
Executable file
@ -0,0 +1,521 @@
|
||||
c Extern "C" declaration has the form:
|
||||
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 *, int *);
|
||||
c
|
||||
c
|
||||
c Call from pair_meam.cpp has the form:
|
||||
c
|
||||
c meam_dens_init_(&i,&nmax,&eflag,&eng_vdwl,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],&errorflag);
|
||||
c
|
||||
|
||||
subroutine meam_dens_init(i, nmax, eflag, eng_vdwl,
|
||||
$ ntype, type, fmap, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ numneigh_full, firstneigh_full,
|
||||
$ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b,
|
||||
$ arho3, arho3b, t_ave, errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer i, nmax, eflag, ntype, type, fmap
|
||||
real*8 eng_vdwl
|
||||
real*8 x
|
||||
integer numneigh, firstneigh, numneigh_full, firstneigh_full
|
||||
real*8 scrfcn, dscrfcn, fcpair
|
||||
real*8 rho0, arho1, arho2
|
||||
real*8 arho2b, arho3, arho3b, t_ave
|
||||
integer errorflag
|
||||
integer j,jn
|
||||
|
||||
dimension x(3,nmax)
|
||||
dimension type(nmax), fmap(ntype)
|
||||
dimension firstneigh(numneigh), firstneigh_full(numneigh_full)
|
||||
dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh)
|
||||
dimension rho0(nmax), arho1(3,nmax), arho2(6,nmax)
|
||||
dimension arho2b(nmax), arho3(10,nmax), arho3b(3,nmax)
|
||||
dimension t_ave(3,nmax)
|
||||
|
||||
errorflag = 0
|
||||
|
||||
c Compute screening function and derivatives
|
||||
call getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ numneigh_full, firstneigh_full,
|
||||
$ ntype, type, fmap)
|
||||
|
||||
c Calculate intermediate density terms to be communicated
|
||||
call calc_rho1(i, nmax, ntype, type, fmap, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ scrfcn, fcpair, rho0, arho1, arho2, arho2b,
|
||||
$ arho3, arho3b, t_ave)
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ numneigh_full, firstneigh_full,
|
||||
$ ntype, type, fmap)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer i, nmax
|
||||
real*8 scrfcn, dscrfcn, fcpair, x
|
||||
integer numneigh, firstneigh, numneigh_full, firstneigh_full
|
||||
integer ntype, type, fmap
|
||||
|
||||
dimension scrfcn(numneigh), dscrfcn(numneigh)
|
||||
dimension fcpair(numneigh), x(3,nmax)
|
||||
dimension firstneigh(numneigh), firstneigh_full(numneigh_full)
|
||||
dimension type(nmax), fmap(ntype)
|
||||
|
||||
integer jn,j,kn,k
|
||||
integer elti,eltj,eltk
|
||||
real*8 xitmp,yitmp,zitmp,delxij,delyij,delzij,rij2,rij
|
||||
real*8 xjtmp,yjtmp,zjtmp,delxik,delyik,delzik,rik2,rik
|
||||
real*8 xktmp,yktmp,zktmp,delxjk,delyjk,delzjk,rjk2,rjk
|
||||
real*8 xik,xjk,sij,fcij,sfcij,dfcij,sikj,dfikj,cikj
|
||||
real*8 Cmin,Cmax,delc,ebound,rbound,a,coef1,coef2
|
||||
real*8 coef1a,coef1b,coef2a,coef2b
|
||||
real*8 dcikj
|
||||
real*8 dC1a,dC1b,dC2a,dC2b
|
||||
real*8 rnorm,fc,dfc,drinv
|
||||
|
||||
drinv = 1.d0/delr_meam
|
||||
elti = fmap(type(i))
|
||||
|
||||
xitmp = x(1,i)
|
||||
yitmp = x(2,i)
|
||||
zitmp = x(3,i)
|
||||
|
||||
do jn = 1,numneigh
|
||||
j = firstneigh(jn)
|
||||
|
||||
c First compute screening function itself, sij
|
||||
xjtmp = x(1,j)
|
||||
yjtmp = x(2,j)
|
||||
zjtmp = x(3,j)
|
||||
delxij = xjtmp - xitmp
|
||||
delyij = yjtmp - yitmp
|
||||
delzij = zjtmp - zitmp
|
||||
rij2 = delxij*delxij + delyij*delyij + delzij*delzij
|
||||
rij = sqrt(rij2)
|
||||
if (rij.gt.rc_meam) then
|
||||
fcij = 0.0
|
||||
dfcij = 0.d0
|
||||
sij = 0.d0
|
||||
else
|
||||
rnorm = (rc_meam-rij)*drinv
|
||||
call screen(i, j, nmax, x, rij2, sij,
|
||||
$ numneigh_full, firstneigh_full, ntype, type, fmap)
|
||||
call dfcut(rnorm,fc,dfc)
|
||||
fcij = fc
|
||||
dfcij = dfc*drinv
|
||||
endif
|
||||
|
||||
c Now compute derivatives
|
||||
dscrfcn(jn) = 0.d0
|
||||
sfcij = sij*fcij
|
||||
if (sfcij.eq.0.d0.or.sfcij.eq.1.d0) goto 100
|
||||
eltj = fmap(type(j))
|
||||
rbound = ebound_meam(elti,eltj) * rij2
|
||||
do kn = 1,numneigh_full
|
||||
k = firstneigh_full(kn)
|
||||
if (k.eq.j) goto 10
|
||||
eltk = fmap(type(k))
|
||||
xktmp = x(1,k)
|
||||
yktmp = x(2,k)
|
||||
zktmp = x(3,k)
|
||||
delxjk = xktmp - xjtmp
|
||||
delyjk = yktmp - yjtmp
|
||||
delzjk = zktmp - zjtmp
|
||||
rjk2 = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk
|
||||
if (rjk2.gt.rbound) goto 10
|
||||
delxik = xktmp - xitmp
|
||||
delyik = yktmp - yitmp
|
||||
delzik = zktmp - zitmp
|
||||
rik2 = delxik*delxik + delyik*delyik + delzik*delzik
|
||||
if (rik2.gt.rbound) goto 10
|
||||
xik = rik2/rij2
|
||||
xjk = rjk2/rij2
|
||||
a = 1 - (xik-xjk)*(xik-xjk)
|
||||
if (a.eq.0.d0) goto 10
|
||||
cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
|
||||
Cmax = Cmax_meam(elti,eltj,eltk)
|
||||
Cmin = Cmin_meam(elti,eltj,eltk)
|
||||
if (cikj.ge.Cmax.or.cikj.lt.0.d0) then
|
||||
goto 10
|
||||
c Note that we never have 0<cikj<Cmin here, else sij=0 (rejected above)
|
||||
else
|
||||
delc = Cmax - Cmin
|
||||
cikj = (cikj-Cmin)/delc
|
||||
call dfcut(cikj,sikj,dfikj)
|
||||
coef1 = dfikj/(delc*sikj)
|
||||
call dCfunc(rij2,rik2,rjk2,dCikj)
|
||||
dscrfcn(jn) = dscrfcn(jn) + coef1*dCikj
|
||||
endif
|
||||
10 continue
|
||||
enddo
|
||||
coef1 = sfcij
|
||||
coef2 = sij*dfcij/rij
|
||||
dscrfcn(jn) = dscrfcn(jn)*coef1 - coef2
|
||||
100 continue
|
||||
|
||||
scrfcn(jn) = sij
|
||||
fcpair(jn) = fcij
|
||||
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine calc_rho1(i, nmax, ntype, type, fmap, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ scrfcn, fcpair, rho0, arho1, arho2, arho2b,
|
||||
$ arho3, arho3b, t_ave)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer i, nmax, ntype, type, fmap
|
||||
real*8 x
|
||||
integer numneigh, firstneigh
|
||||
real*8 scrfcn, fcpair, rho0, arho1, arho2
|
||||
real*8 arho2b, arho3, arho3b, t_ave
|
||||
|
||||
dimension type(nmax), fmap(ntype), x(3,nmax)
|
||||
dimension firstneigh(numneigh)
|
||||
dimension scrfcn(numneigh), fcpair(numneigh)
|
||||
dimension rho0(nmax), arho1(3,nmax), arho2(6,nmax)
|
||||
dimension arho2b(nmax), arho3(10,nmax), arho3b(3,nmax)
|
||||
dimension t_ave(3,nmax)
|
||||
|
||||
integer jn,j,m,n,p,elti,eltj
|
||||
integer nv2,nv3
|
||||
real*8 xtmp,ytmp,ztmp,delij(3),rij2,rij,sij
|
||||
real*8 ai,aj,rhoa0j,rhoa1j,rhoa2j,rhoa3j,A1j,A2j,A3j
|
||||
real*8 G,Gbar,gam,shp(3)
|
||||
real*8 ro0i,ro0j
|
||||
real*8 rhoa0i,rhoa1i,rhoa2i,rhoa3i,A1i,A2i,A3i
|
||||
|
||||
elti = fmap(type(i))
|
||||
xtmp = x(1,i)
|
||||
ytmp = x(2,i)
|
||||
ztmp = x(3,i)
|
||||
do jn = 1,numneigh
|
||||
if (scrfcn(jn).ne.0.d0) then
|
||||
j = firstneigh(jn)
|
||||
sij = scrfcn(jn)*fcpair(jn)
|
||||
delij(1) = x(1,j) - xtmp
|
||||
delij(2) = x(2,j) - ytmp
|
||||
delij(3) = x(3,j) - ztmp
|
||||
rij2 = delij(1)*delij(1) + delij(2)*delij(2)
|
||||
$ + delij(3)*delij(3)
|
||||
if (rij2.lt.cutforcesq) then
|
||||
eltj = fmap(type(j))
|
||||
rij = sqrt(rij2)
|
||||
ai = rij/re_meam(elti,elti) - 1.d0
|
||||
aj = rij/re_meam(eltj,eltj) - 1.d0
|
||||
ro0i = rho0_meam(elti)
|
||||
ro0j = rho0_meam(eltj)
|
||||
rhoa0j = ro0j*exp(-beta0_meam(eltj)*aj)*sij
|
||||
rhoa1j = ro0j*exp(-beta1_meam(eltj)*aj)*sij
|
||||
rhoa2j = ro0j*exp(-beta2_meam(eltj)*aj)*sij
|
||||
rhoa3j = ro0j*exp(-beta3_meam(eltj)*aj)*sij
|
||||
rhoa0i = ro0i*exp(-beta0_meam(elti)*ai)*sij
|
||||
rhoa1i = ro0i*exp(-beta1_meam(elti)*ai)*sij
|
||||
rhoa2i = ro0i*exp(-beta2_meam(elti)*ai)*sij
|
||||
rhoa3i = ro0i*exp(-beta3_meam(elti)*ai)*sij
|
||||
rho0(i) = rho0(i) + rhoa0j
|
||||
rho0(j) = rho0(j) + rhoa0i
|
||||
t_ave(1,i) = t_ave(1,i) + t1_meam(eltj)*rhoa0j
|
||||
t_ave(2,i) = t_ave(2,i) + t2_meam(eltj)*rhoa0j
|
||||
t_ave(3,i) = t_ave(3,i) + t3_meam(eltj)*rhoa0j
|
||||
t_ave(1,j) = t_ave(1,j) + t1_meam(elti)*rhoa0i
|
||||
t_ave(2,j) = t_ave(2,j) + t2_meam(elti)*rhoa0i
|
||||
t_ave(3,j) = t_ave(3,j) + t3_meam(elti)*rhoa0i
|
||||
Arho2b(i) = Arho2b(i) + rhoa2j
|
||||
Arho2b(j) = Arho2b(j) + rhoa2i
|
||||
|
||||
A1j = rhoa1j/rij
|
||||
A2j = rhoa2j/rij2
|
||||
A3j = rhoa3j/(rij2*rij)
|
||||
A1i = rhoa1i/rij
|
||||
A2i = rhoa2i/rij2
|
||||
A3i = rhoa3i/(rij2*rij)
|
||||
nv2 = 1
|
||||
nv3 = 1
|
||||
do m = 1,3
|
||||
Arho1(m,i) = Arho1(m,i) + A1j*delij(m)
|
||||
Arho1(m,j) = Arho1(m,j) - A1i*delij(m)
|
||||
Arho3b(m,i) = Arho3b(m,i) + rhoa3j*delij(m)/rij
|
||||
Arho3b(m,j) = Arho3b(m,j) - rhoa3i*delij(m)/rij
|
||||
do n = m,3
|
||||
Arho2(nv2,i) = Arho2(nv2,i) + A2j*delij(m)*delij(n)
|
||||
Arho2(nv2,j) = Arho2(nv2,j) + A2i*delij(m)*delij(n)
|
||||
nv2 = nv2+1
|
||||
do p = n,3
|
||||
Arho3(nv3,i) = Arho3(nv3,i)
|
||||
$ + A3j*delij(m)*delij(n)*delij(p)
|
||||
Arho3(nv3,j) = Arho3(nv3,j)
|
||||
$ - A3i*delij(m)*delij(n)*delij(p)
|
||||
nv3 = nv3+1
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine screen(i, j, nmax, x, rijsq, sij,
|
||||
$ numneigh_full, firstneigh_full, ntype, type, fmap)
|
||||
c Screening function
|
||||
c Inputs: i = atom 1 id (integer)
|
||||
c j = atom 2 id (integer)
|
||||
c rijsq = squared distance between i and j
|
||||
c Outputs: sij = screening function
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer i,j,nmax,k,nk,m
|
||||
real*8 x,rijsq,sij
|
||||
integer numneigh_full, firstneigh_full
|
||||
integer ntype, type, fmap
|
||||
|
||||
dimension x(3,nmax), firstneigh_full(numneigh_full)
|
||||
dimension type(nmax), fmap(ntype)
|
||||
|
||||
integer elti,eltj,eltk
|
||||
real*8 delxik,delyik,delzik
|
||||
real*8 delxjk,delyjk,delzjk
|
||||
real*8 riksq,rjksq,xik,xjk,cikj,a,delc,sikj,fcij,rij
|
||||
real*8 Cmax,Cmin,rbound
|
||||
|
||||
sij = 1.d0
|
||||
elti = fmap(type(i))
|
||||
eltj = fmap(type(j))
|
||||
|
||||
c if rjksq > ebound*rijsq, atom k is definitely outside the ellipse
|
||||
rbound = ebound_meam(elti,eltj)*rijsq
|
||||
|
||||
do nk = 1,numneigh_full
|
||||
k = firstneigh_full(nk)
|
||||
eltk = fmap(type(k))
|
||||
if (k.eq.j) goto 10
|
||||
delxjk = x(1,k) - x(1,j)
|
||||
delyjk = x(2,k) - x(2,j)
|
||||
delzjk = x(3,k) - x(3,j)
|
||||
rjksq = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk
|
||||
if (rjksq.gt.rbound) goto 10
|
||||
delxik = x(1,k) - x(1,i)
|
||||
delyik = x(2,k) - x(2,i)
|
||||
delzik = x(3,k) - x(3,i)
|
||||
riksq = delxik*delxik + delyik*delyik + delzik*delzik
|
||||
if (riksq.gt.rbound) goto 10
|
||||
xik = riksq/rijsq
|
||||
xjk = rjksq/rijsq
|
||||
a = 1 - (xik-xjk)*(xik-xjk)
|
||||
if (a.eq.0.d0) goto 10
|
||||
cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
|
||||
Cmax = Cmax_meam(elti,eltj,eltk)
|
||||
Cmin = Cmin_meam(elti,eltj,eltk)
|
||||
if (cikj.ge.Cmax.or.cikj.lt.0.d0) then
|
||||
goto 10
|
||||
else if (cikj.le.Cmin) then
|
||||
sij = 0.d0
|
||||
goto 20
|
||||
else
|
||||
delc = Cmax - Cmin
|
||||
cikj = (cikj-Cmin)/delc
|
||||
call fcut(cikj,sikj)
|
||||
endif
|
||||
sij = sij * sikj
|
||||
10 continue
|
||||
enddo
|
||||
|
||||
20 continue
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2,
|
||||
$ ntype,type,fmap,x,scrfcn,fcpair)
|
||||
c Inputs: i,j,k = id's of 3 atom triplet
|
||||
c jn = id of i-j pair
|
||||
c rij2 = squared distance between i and j
|
||||
c Outputs: dsij1 = deriv. of sij w.r.t. rik
|
||||
c dsij2 = deriv. of sij w.r.t. rjk
|
||||
use meam_data
|
||||
implicit none
|
||||
integer i,j,k,jn,nmax,numneigh
|
||||
integer elti,eltj,eltk
|
||||
real*8 rij2,rik2,rjk2,dsij1,dsij2
|
||||
integer ntype, type, fmap
|
||||
real*8 x, scrfcn, fcpair
|
||||
|
||||
dimension type(nmax), fmap(ntype)
|
||||
dimension x(3,nmax), scrfcn(numneigh), fcpair(numneigh)
|
||||
|
||||
real*8 dxik,dyik,dzik
|
||||
real*8 dxjk,dyjk,dzjk
|
||||
real*8 rbound,delc,sij,xik,xjk,cikj,sikj,dfc,a
|
||||
real*8 Cmax,Cmin,dCikj1,dCikj2
|
||||
|
||||
sij = scrfcn(jn)*fcpair(jn)
|
||||
elti = fmap(type(i))
|
||||
eltj = fmap(type(j))
|
||||
eltk = fmap(type(k))
|
||||
Cmax = Cmax_meam(elti,eltj,eltk)
|
||||
Cmin = Cmin_meam(elti,eltj,eltk)
|
||||
|
||||
dsij1 = 0.d0
|
||||
dsij2 = 0.d0
|
||||
if ((sij.ne.0.d0).and.(sij.ne.1.d0)) then
|
||||
rbound = rij2*ebound_meam(elti,eltj)
|
||||
delc = Cmax-Cmin
|
||||
dxjk = x(1,k) - x(1,j)
|
||||
dyjk = x(2,k) - x(2,j)
|
||||
dzjk = x(3,k) - x(3,j)
|
||||
rjk2 = dxjk*dxjk + dyjk*dyjk + dzjk*dzjk
|
||||
if (rjk2.le.rbound) then
|
||||
dxik = x(1,k) - x(1,i)
|
||||
dyik = x(2,k) - x(2,i)
|
||||
dzik = x(3,k) - x(3,i)
|
||||
rik2 = dxik*dxik + dyik*dyik + dzik*dzik
|
||||
if (rik2.le.rbound) then
|
||||
xik = rik2/rij2
|
||||
xjk = rjk2/rij2
|
||||
a = 1 - (xik-xjk)*(xik-xjk)
|
||||
if (a.ne.0.d0) then
|
||||
cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
|
||||
if (cikj.ge.Cmin.and.cikj.le.Cmax) then
|
||||
cikj = (cikj-Cmin)/delc
|
||||
call dfcut(cikj,sikj,dfc)
|
||||
call dCfunc2(rij2,rik2,rjk2,dCikj1,dCikj2)
|
||||
a = sij/delc*dfc/sikj
|
||||
dsij1 = a*dCikj1
|
||||
dsij2 = a*dCikj2
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine fcut(xi,fc)
|
||||
c cutoff function
|
||||
implicit none
|
||||
real*8 xi,fc
|
||||
real*8 a
|
||||
if (xi.ge.1.d0) then
|
||||
fc = 1.d0
|
||||
else if (xi.le.0.d0) then
|
||||
fc = 0.d0
|
||||
else
|
||||
a = 1.d0-xi
|
||||
a = a*a
|
||||
a = a*a
|
||||
a = 1.d0-a
|
||||
fc = a*a
|
||||
c fc = xi
|
||||
endif
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dfcut(xi,fc,dfc)
|
||||
c cutoff function and its derivative
|
||||
implicit none
|
||||
real*8 xi,fc,dfc,a,a3,a4
|
||||
if (xi.ge.1.d0) then
|
||||
fc = 1.d0
|
||||
dfc = 0.d0
|
||||
else if (xi.le.0.d0) then
|
||||
fc = 0.d0
|
||||
dfc = 0.d0
|
||||
else
|
||||
a = 1.d0-xi
|
||||
a3 = a*a*a
|
||||
a4 = a*a3
|
||||
fc = (1.d0-a4)**2
|
||||
dfc = 8*(1.d0-a4)*a3
|
||||
c fc = xi
|
||||
c dfc = 1.d0
|
||||
endif
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dCfunc(rij2,rik2,rjk2,dCikj)
|
||||
c Inputs: rij,rij2,rik2,rjk2
|
||||
c Outputs: dCikj = derivative of Cikj w.r.t. rij
|
||||
implicit none
|
||||
real*8 rij2,rik2,rjk2,dCikj
|
||||
real*8 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
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
subroutine dCfunc2(rij2,rik2,rjk2,dCikj1,dCikj2)
|
||||
c Inputs: rij,rij2,rik2,rjk2
|
||||
c Outputs: dCikj1 = derivative of Cikj w.r.t. rik
|
||||
c dCikj2 = derivative of Cikj w.r.t. rjk
|
||||
implicit none
|
||||
real*8 rij2,rik2,rjk2,dCikj1,dCikj2
|
||||
real*8 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
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
|
||||
|
||||
|
||||
|
||||
438
lib/meam/meam_force.F
Executable file
438
lib/meam/meam_force.F
Executable file
@ -0,0 +1,438 @@
|
||||
c Extern "C" declaration has the form:
|
||||
c
|
||||
c void meam_force_(int *, int *, int *, double *, int *, int *, int *, 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 *, int *);
|
||||
c
|
||||
c Call from pair_meam.cpp has the form:
|
||||
c
|
||||
c meam_force_(&i,&nmax,&eflag,&eng_vdwl,&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 dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
|
||||
c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
|
||||
c &t_ave[0][0],&f[0][0],&strssa[0][0][0],&errorflag);
|
||||
c
|
||||
|
||||
subroutine meam_force(i, nmax, eflag, eng_vdwl,
|
||||
$ ntype, type, fmap, x,
|
||||
$ numneigh, firstneigh, numneigh_full, firstneigh_full,
|
||||
$ scrfcn, dscrfcn, fcpair,
|
||||
$ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp,
|
||||
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, f, strssa,
|
||||
$ errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer eflag, nmax, ntype, type, fmap
|
||||
real*8 eng_vdwl, x
|
||||
integer numneigh, firstneigh, numneigh_full, firstneigh_full
|
||||
real*8 scrfcn, dscrfcn, fcpair
|
||||
real*8 dGamma1, dGamma2, dGamma3
|
||||
real*8 rho0, rho1, rho2, rho3, fp
|
||||
real*8 Arho1, Arho2, Arho2b
|
||||
real*8 Arho3, Arho3b
|
||||
real*8 t_ave, f, strssa
|
||||
integer errorflag
|
||||
|
||||
dimension type(nmax), fmap(ntype)
|
||||
dimension x(3,nmax)
|
||||
dimension firstneigh(numneigh), firstneigh_full(numneigh_full)
|
||||
dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh)
|
||||
dimension dGamma1(nmax), dGamma2(nmax), dGamma3(nmax)
|
||||
dimension rho0(nmax), rho1(nmax), rho2(nmax), rho3(nmax), fp(nmax)
|
||||
dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax)
|
||||
dimension Arho3(10,nmax), Arho3b(3,nmax)
|
||||
dimension t_ave(3,nmax), f(3,nmax), strssa(3,3,nmax)
|
||||
|
||||
integer i,j,jn,k,kn,kk,m,n,p,q,strscomp
|
||||
integer nv2,nv3,elti,eltj,ind
|
||||
real*8 xitmp,yitmp,zitmp,delij(3),delref(3),rij2,rij,rij3
|
||||
real*8 delik(3),deljk(3)
|
||||
real*8 Eu,astar,astarp
|
||||
real*8 pp,phiforce,dUdrij,dUdsij,dUdrijm(3),force,forcem
|
||||
real*8 B,r,recip,phi,phip,rhop,a
|
||||
real*8 sij,fcij,dfcij,ds(3)
|
||||
real*8 a0,a1,a1i,a1j,a2,a2i,a2j
|
||||
real*8 a3i,a3j,a3i1,a3i2,a3j1,a3j2
|
||||
real*8 G,dG,Gbar,dGbar,gam,shpi(3),shpj(3),Z,denom
|
||||
real*8 ai,aj,ro0i,ro0j,invrei,invrej
|
||||
real*8 b0,rhoa0j,drhoa0j,rhoa0i,drhoa0i
|
||||
real*8 b1,rhoa1j,drhoa1j,rhoa1i,drhoa1i
|
||||
real*8 b2,rhoa2j,drhoa2j,rhoa2i,drhoa2i
|
||||
real*8 a3,a3a,b3,rhoa3j,drhoa3j,rhoa3i,drhoa3i
|
||||
real*8 drho0dr1,drho0dr2,drho0ds1,drho0ds2
|
||||
real*8 drho1dr1,drho1dr2,drho1ds1,drho1ds2
|
||||
real*8 drho1drm1(3),drho1drm2(3)
|
||||
real*8 drho2dr1,drho2dr2,drho2ds1,drho2ds2
|
||||
real*8 drho2drm1(3),drho2drm2(3)
|
||||
real*8 drho3dr1,drho3dr2,drho3ds1,drho3ds2
|
||||
real*8 drho3drm1(3),drho3drm2(3)
|
||||
real*8 dt1dr1,dt1dr2,dt1ds1,dt1ds2
|
||||
real*8 dt2dr1,dt2dr2,dt2ds1,dt2ds2
|
||||
real*8 dt3dr1,dt3dr2,dt3ds1,dt3ds2
|
||||
real*8 drhodr1,drhodr2,drhods1,drhods2,drhodrm1(3),drhodrm2(3)
|
||||
real*8 arg,arg1,arg2
|
||||
real*8 arg1i1,arg1j1,arg1i2,arg1j2,arg2i2,arg2j2
|
||||
real*8 arg1i3,arg1j3,arg2i3,arg2j3,arg3i3,arg3j3
|
||||
real*8 dsij1,dsij2,force1,force2
|
||||
real*8 t1i,t2i,t3i,t1j,t2j,t3j
|
||||
|
||||
errorflag = 0
|
||||
|
||||
strscomp = 0
|
||||
|
||||
c Compute forces atom i
|
||||
|
||||
elti = fmap(type(i))
|
||||
xitmp = x(1,i)
|
||||
yitmp = x(2,i)
|
||||
zitmp = x(3,i)
|
||||
|
||||
c Treat each pair
|
||||
do jn = 1,numneigh
|
||||
|
||||
if (scrfcn(jn).ne.0.d0) then
|
||||
j = firstneigh(jn)
|
||||
|
||||
sij = scrfcn(jn)*fcpair(jn)
|
||||
delij(1) = x(1,j) - xitmp
|
||||
delij(2) = x(2,j) - yitmp
|
||||
delij(3) = x(3,j) - zitmp
|
||||
rij2 = delij(1)*delij(1) + delij(2)*delij(2)
|
||||
$ + delij(3)*delij(3)
|
||||
if (rij2.lt.cutforcesq) then
|
||||
rij = sqrt(rij2)
|
||||
r = rij
|
||||
eltj = fmap(type(j))
|
||||
|
||||
c Compute phi and phip
|
||||
ind = eltind(elti,eltj)
|
||||
pp = rij*rdrar + 1.0D0
|
||||
kk = pp
|
||||
kk = min(kk,nrar-1)
|
||||
pp = pp - kk
|
||||
pp = min(pp,1.0D0)
|
||||
phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp
|
||||
$ + phirar1(kk,ind))*pp + phirar(kk,ind)
|
||||
phip = (phirar6(kk,ind)*pp + phirar5(kk,ind))*pp
|
||||
$ + phirar4(kk,ind)
|
||||
recip = 1.0d0/r
|
||||
|
||||
if (eflag.eq.1) eng_vdwl = eng_vdwl + phi*sij
|
||||
|
||||
c write(1,*) "force_meamf: phi: ",phi
|
||||
c write(1,*) "force_meamf: phip: ",phip
|
||||
|
||||
c Compute pair densities and derivatives
|
||||
invrei = 1.d0/re_meam(elti,elti)
|
||||
ai = rij*invrei - 1.d0
|
||||
ro0i = rho0_meam(elti)
|
||||
rhoa0i = ro0i*exp(-beta0_meam(elti)*ai)
|
||||
drhoa0i = -beta0_meam(elti)*invrei*rhoa0i
|
||||
rhoa1i = ro0i*exp(-beta1_meam(elti)*ai)
|
||||
drhoa1i = -beta1_meam(elti)*invrei*rhoa1i
|
||||
rhoa2i = ro0i*exp(-beta2_meam(elti)*ai)
|
||||
drhoa2i = -beta2_meam(elti)*invrei*rhoa2i
|
||||
rhoa3i = ro0i*exp(-beta3_meam(elti)*ai)
|
||||
drhoa3i = -beta3_meam(elti)*invrei*rhoa3i
|
||||
|
||||
if (elti.ne.eltj) then
|
||||
invrej = 1.d0/re_meam(eltj,eltj)
|
||||
aj = rij*invrej - 1.d0
|
||||
ro0j = rho0_meam(eltj)
|
||||
rhoa0j = ro0j*exp(-beta0_meam(eltj)*aj)
|
||||
drhoa0j = -beta0_meam(eltj)*invrej*rhoa0j
|
||||
rhoa1j = ro0j*exp(-beta1_meam(eltj)*aj)
|
||||
drhoa1j = -beta1_meam(eltj)*invrej*rhoa1j
|
||||
rhoa2j = ro0j*exp(-beta2_meam(eltj)*aj)
|
||||
drhoa2j = -beta2_meam(eltj)*invrej*rhoa2j
|
||||
rhoa3j = ro0j*exp(-beta3_meam(eltj)*aj)
|
||||
drhoa3j = -beta3_meam(eltj)*invrej*rhoa3j
|
||||
else
|
||||
rhoa0j = rhoa0i
|
||||
drhoa0j = drhoa0i
|
||||
rhoa1j = rhoa1i
|
||||
drhoa1j = drhoa1i
|
||||
rhoa2j = rhoa2i
|
||||
drhoa2j = drhoa2i
|
||||
rhoa3j = rhoa3i
|
||||
drhoa3j = drhoa3i
|
||||
endif
|
||||
|
||||
nv2 = 1
|
||||
nv3 = 1
|
||||
arg1i1 = 0.d0
|
||||
arg1j1 = 0.d0
|
||||
arg1i2 = 0.d0
|
||||
arg1j2 = 0.d0
|
||||
arg1i3 = 0.d0
|
||||
arg1j3 = 0.d0
|
||||
arg3i3 = 0.d0
|
||||
arg3j3 = 0.d0
|
||||
do n = 1,3
|
||||
do p = n,3
|
||||
do q = p,3
|
||||
arg = delij(n)*delij(p)*delij(q)*v3D(nv3)
|
||||
arg1i3 = arg1i3 + Arho3(nv3,i)*arg
|
||||
arg1j3 = arg1j3 - Arho3(nv3,j)*arg
|
||||
nv3 = nv3+1
|
||||
enddo
|
||||
arg = delij(n)*delij(p)*v2D(nv2)
|
||||
arg1i2 = arg1i2 + Arho2(nv2,i)*arg
|
||||
arg1j2 = arg1j2 + Arho2(nv2,j)*arg
|
||||
nv2 = nv2+1
|
||||
enddo
|
||||
arg1i1 = arg1i1 + Arho1(n,i)*delij(n)
|
||||
arg1j1 = arg1j1 - Arho1(n,j)*delij(n)
|
||||
arg3i3 = arg3i3 + Arho3b(n,i)*delij(n)
|
||||
arg3j3 = arg3j3 - Arho3b(n,j)*delij(n)
|
||||
enddo
|
||||
|
||||
c rho0 terms
|
||||
drho0dr1 = drhoa0j * sij
|
||||
drho0dr2 = drhoa0i * sij
|
||||
|
||||
c rho1 terms
|
||||
a1 = 2*sij/rij
|
||||
drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1
|
||||
drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1
|
||||
a1 = 2.d0*sij/rij
|
||||
do m = 1,3
|
||||
drho1drm1(m) = a1*rhoa1j*Arho1(m,i)
|
||||
drho1drm2(m) = -a1*rhoa1i*Arho1(m,j)
|
||||
enddo
|
||||
|
||||
c rho2 terms
|
||||
a2 = 2*sij/rij2
|
||||
drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2
|
||||
$ - 2.d0/3.d0*Arho2b(i)*drhoa2j*sij
|
||||
drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2
|
||||
$ - 2.d0/3.d0*Arho2b(j)*drhoa2i*sij
|
||||
a2 = 4*sij/rij2
|
||||
do m = 1,3
|
||||
drho2drm1(m) = 0.d0
|
||||
drho2drm2(m) = 0.d0
|
||||
do n = 1,3
|
||||
drho2drm1(m) = drho2drm1(m)
|
||||
$ + Arho2(vind2D(m,n),i)*delij(n)
|
||||
drho2drm2(m) = drho2drm2(m)
|
||||
$ - Arho2(vind2D(m,n),j)*delij(n)
|
||||
enddo
|
||||
drho2drm1(m) = a2*rhoa2j*drho2drm1(m)
|
||||
drho2drm2(m) = -a2*rhoa2i*drho2drm2(m)
|
||||
enddo
|
||||
|
||||
c rho3 terms
|
||||
rij3 = rij*rij2
|
||||
a3 = 2*sij/rij3
|
||||
a3a = 6.d0/5.d0*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)
|
||||
do m = 1,3
|
||||
drho3drm1(m) = 0.d0
|
||||
drho3drm2(m) = 0.d0
|
||||
nv2 = 1
|
||||
do n = 1,3
|
||||
do p = n,3
|
||||
arg = delij(n)*delij(p)*v2D(nv2)
|
||||
drho3drm1(m) = drho3drm1(m)
|
||||
$ + Arho3(vind3D(m,n,p),i)*arg
|
||||
drho3drm2(m) = drho3drm2(m)
|
||||
$ + Arho3(vind3D(m,n,p),j)*arg
|
||||
nv2 = nv2 + 1
|
||||
enddo
|
||||
enddo
|
||||
drho3drm1(m) = (a3*drho3drm1(m) - a3a*Arho3b(m,i))
|
||||
$ *rhoa3j
|
||||
drho3drm2(m) = (-a3*drho3drm2(m) + a3a*Arho3b(m,j))
|
||||
$ *rhoa3i
|
||||
enddo
|
||||
|
||||
c Compute derivatives of weighting functions t wrt rij
|
||||
t1i = t_ave(1,i)
|
||||
t2i = t_ave(2,i)
|
||||
t3i = t_ave(3,i)
|
||||
t1j = t_ave(1,j)
|
||||
t2j = t_ave(2,j)
|
||||
t3j = t_ave(3,j)
|
||||
|
||||
ai = 0.d0
|
||||
if( rho0(i) .ne. 0.d0 ) then
|
||||
ai = drhoa0j*sij/rho0(i)
|
||||
end if
|
||||
aj = 0.d0
|
||||
if( rho0(j) .ne. 0.d0 ) then
|
||||
aj = drhoa0i*sij/rho0(j)
|
||||
end if
|
||||
|
||||
dt1dr1 = ai*(t1_meam(eltj)-t1i)
|
||||
dt1dr2 = aj*(t1_meam(elti)-t1j)
|
||||
dt2dr1 = ai*(t2_meam(eltj)-t2i)
|
||||
dt2dr2 = aj*(t2_meam(elti)-t2j)
|
||||
dt3dr1 = ai*(t3_meam(eltj)-t3i)
|
||||
dt3dr2 = aj*(t3_meam(elti)-t3j)
|
||||
|
||||
c Compute derivatives of total density wrt rij, sij and rij(3)
|
||||
call get_shpfcn(shpi,lattce_meam(elti,elti))
|
||||
call get_shpfcn(shpj,lattce_meam(eltj,eltj))
|
||||
drhodr1 = dGamma1(i)*drho0dr1
|
||||
$ + dGamma2(i)*
|
||||
$ (dt1dr1*rho1(i)+t1i*drho1dr1
|
||||
$ + dt2dr1*rho2(i)+t2i*drho2dr1
|
||||
$ + dt3dr1*rho3(i)+t3i*drho3dr1)
|
||||
$ - dGamma3(i)*
|
||||
$ (shpi(1)*dt1dr1+shpi(2)*dt2dr1+shpi(3)*dt3dr1)
|
||||
drhodr2 = dGamma1(j)*drho0dr2
|
||||
$ + dGamma2(j)*
|
||||
$ (dt1dr2*rho1(j)+t1j*drho1dr2
|
||||
$ + dt2dr2*rho2(j)+t2j*drho2dr2
|
||||
$ + dt3dr2*rho3(j)+t3j*drho3dr2)
|
||||
$ - dGamma3(j)*
|
||||
$ (shpj(1)*dt1dr2+shpj(2)*dt2dr2+shpj(3)*dt3dr2)
|
||||
do m = 1,3
|
||||
drhodrm1(m) = 0.d0
|
||||
drhodrm2(m) = 0.d0
|
||||
drhodrm1(m) = dGamma2(i)*
|
||||
$ (t1i*drho1drm1(m)
|
||||
$ + t2i*drho2drm1(m)
|
||||
$ + t3i*drho3drm1(m))
|
||||
drhodrm2(m) = dGamma2(j)*
|
||||
$ (t1j*drho1drm2(m)
|
||||
$ + t2j*drho2drm2(m)
|
||||
$ + t3j*drho3drm2(m))
|
||||
enddo
|
||||
|
||||
c Compute derivatives wrt sij, but only if necessary
|
||||
if (dscrfcn(jn).ne.0.d0) then
|
||||
drho0ds1 = rhoa0j
|
||||
drho0ds2 = rhoa0i
|
||||
a1 = 2.d0/rij
|
||||
drho1ds1 = a1*rhoa1j*arg1i1
|
||||
drho1ds2 = a1*rhoa1i*arg1j1
|
||||
a2 = 2.d0/rij2
|
||||
drho2ds1 = a2*rhoa2j*arg1i2
|
||||
$ - 2.d0/3.d0*Arho2b(i)*rhoa2j
|
||||
drho2ds2 = a2*rhoa2i*arg1j2
|
||||
$ - 2.d0/3.d0*Arho2b(j)*rhoa2i
|
||||
a3 = 2.d0/rij3
|
||||
a3a = 6.d0/(5.d0*rij)
|
||||
drho3ds1 = a3*rhoa3j*arg1i3 - a3a*rhoa3j*arg3i3
|
||||
drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3
|
||||
|
||||
ai = 0.d0
|
||||
if( rho0(i) .ne. 0.d0 ) then
|
||||
ai = rhoa0j/rho0(i)
|
||||
end if
|
||||
aj = 0.d0
|
||||
if( rho0(j) .ne. 0.d0 ) then
|
||||
aj = rhoa0i/rho0(j)
|
||||
end if
|
||||
|
||||
dt1ds1 = ai*(t1_meam(eltj)-t1i)
|
||||
dt1ds2 = aj*(t1_meam(elti)-t1j)
|
||||
dt2ds1 = ai*(t2_meam(eltj)-t2i)
|
||||
dt2ds2 = aj*(t2_meam(elti)-t2j)
|
||||
dt3ds1 = ai*(t3_meam(eltj)-t3i)
|
||||
dt3ds2 = aj*(t3_meam(elti)-t3j)
|
||||
drhods1 = dGamma1(i)*drho0ds1
|
||||
$ + dGamma2(i)*
|
||||
$ (dt1ds1*rho1(i)+t1i*drho1ds1
|
||||
$ + dt2ds1*rho2(i)+t2i*drho2ds1
|
||||
$ + dt3ds1*rho3(i)+t3i*drho3ds1)
|
||||
$ - dGamma3(i)*
|
||||
$ (shpi(1)*dt1ds1+shpi(2)*dt2ds1+shpi(3)*dt3ds1)
|
||||
drhods2 = dGamma1(j)*drho0ds2
|
||||
$ + dGamma2(j)*
|
||||
$ (dt1ds2*rho1(j)+t1j*drho1ds2
|
||||
$ + dt2ds2*rho2(j)+t2j*drho2ds2
|
||||
$ + dt3ds2*rho3(j)+t3j*drho3ds2)
|
||||
$ - dGamma3(j)*
|
||||
$ (shpj(1)*dt1ds2+shpj(2)*dt2ds2+shpj(3)*dt3ds2)
|
||||
endif
|
||||
|
||||
c Compute derivatives of energy wrt rij, sij and rij(3)
|
||||
dUdrij = phip*sij
|
||||
$ + fp(i)*drhodr1 + fp(j)*drhodr2
|
||||
dUdsij = 0.d0
|
||||
if (dscrfcn(jn).ne.0.d0) then
|
||||
dUdsij = phi
|
||||
$ + fp(i)*drhods1 + fp(j)*drhods2
|
||||
endif
|
||||
do m = 1,3
|
||||
dUdrijm(m) = fp(i)*drhodrm1(m) + fp(j)*drhodrm2(m)
|
||||
enddo
|
||||
|
||||
c Add the part of the force due to dUdrij and dUdsij
|
||||
force = dUdrij*recip + dUdsij*dscrfcn(jn)
|
||||
do m = 1,3
|
||||
forcem = delij(m)*force + dUdrijm(m)
|
||||
f(m,i) = f(m,i) + forcem
|
||||
f(m,j) = f(m,j) - forcem
|
||||
if (strscomp>0) then
|
||||
do n = 1,3
|
||||
arg = delij(n)*forcem
|
||||
strssa(n,m,i) = strssa(n,m,i) + arg
|
||||
strssa(n,m,j) = strssa(n,m,j) + arg
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
c write(1,*)"forcei1: ",f(1,i)," ",f(2,i)," ",f(3,i)
|
||||
c write(1,*)"forcej1: ",f(1,j)," ",f(2,j)," ",f(3,j)
|
||||
|
||||
c Now compute forces on other atoms k due to change in sij
|
||||
if (sij.eq.0.d0.or.sij.eq.1.d0) goto 100
|
||||
do kn = 1,numneigh_full
|
||||
k = firstneigh_full(kn)
|
||||
if (k.ne.j) then
|
||||
call dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2,
|
||||
$ ntype,type,fmap,x,scrfcn,fcpair)
|
||||
if (dsij1.ne.0.d0.or.dsij2.ne.0.d0) then
|
||||
force1 = dUdsij*dsij1
|
||||
force2 = dUdsij*dsij2
|
||||
do m = 1,3
|
||||
delik(m) = x(m,k) - x(m,i)
|
||||
deljk(m) = x(m,k) - x(m,j)
|
||||
enddo
|
||||
do m = 1,3
|
||||
f(m,i) = f(m,i) + force1*delik(m)
|
||||
f(m,j) = f(m,j) + force2*deljk(m)
|
||||
f(m,k) = f(m,k) - force1*delik(m)
|
||||
$ - force2*deljk(m)
|
||||
if (strscomp>0) then
|
||||
do n = m,3
|
||||
arg1 = force1*delik(m)*delik(n)
|
||||
arg2 = force2*deljk(m)*deljk(n)
|
||||
strssa(m,n,i) = strssa(m,n,i) + arg1
|
||||
strssa(m,n,j) = strssa(m,n,j) + arg2
|
||||
strssa(m,n,k) = strssa(m,n,k) +
|
||||
$ arg1 + arg2
|
||||
if (n.ne.m) then
|
||||
strssa(n,m,i) = strssa(n,m,i) + arg1
|
||||
strssa(n,m,j) = strssa(n,m,j) + arg2
|
||||
strssa(n,m,k) = strssa(n,m,k)
|
||||
$ + arg1 + arg2
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
|
||||
endif
|
||||
endif
|
||||
c end of k loop
|
||||
enddo
|
||||
endif
|
||||
100 continue
|
||||
endif
|
||||
c end of j loop
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
733
lib/meam/meam_setup_done.F
Executable file
733
lib/meam/meam_setup_done.F
Executable file
@ -0,0 +1,733 @@
|
||||
c Declaration in pair_meam.h:
|
||||
c
|
||||
c void meam_setup_done(double *)
|
||||
c
|
||||
c Call from pair_meam.cpp:
|
||||
c
|
||||
c meam_setup_done(&cutmax)
|
||||
c
|
||||
|
||||
subroutine meam_setup_done(cutmax)
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
real*8 cutmax
|
||||
|
||||
integer nv2, nv3, m, n, p
|
||||
|
||||
c Force cutoff
|
||||
cutforce = rc_meam
|
||||
cutforcesq = cutforce*cutforce
|
||||
|
||||
c Pass cutoff back to calling program
|
||||
cutmax = cutforce
|
||||
|
||||
c Augment t1 term
|
||||
t1_meam(:) = t1_meam(:) + augt1 * 3.d0/5.d0 * t3_meam(:)
|
||||
|
||||
c Compute off-diagonal alloy parameters
|
||||
call alloyparams
|
||||
|
||||
c indices and factors for Voight notation
|
||||
nv2 = 1
|
||||
nv3 = 1
|
||||
do m = 1,3
|
||||
do n = m,3
|
||||
vind2D(m,n) = nv2
|
||||
vind2D(n,m) = nv2
|
||||
nv2 = nv2+1
|
||||
do p = n,3
|
||||
vind3D(m,n,p) = nv3
|
||||
vind3D(m,p,n) = nv3
|
||||
vind3D(n,m,p) = nv3
|
||||
vind3D(n,p,m) = nv3
|
||||
vind3D(p,m,n) = nv3
|
||||
vind3D(p,n,m) = nv3
|
||||
nv3 = nv3+1
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
v2D(1) = 1
|
||||
v2D(2) = 2
|
||||
v2D(3) = 2
|
||||
v2D(4) = 1
|
||||
v2D(5) = 2
|
||||
v2D(6) = 1
|
||||
|
||||
v3D(1) = 1
|
||||
v3D(2) = 3
|
||||
v3D(3) = 3
|
||||
v3D(4) = 3
|
||||
v3D(5) = 6
|
||||
v3D(6) = 3
|
||||
v3D(7) = 1
|
||||
v3D(8) = 3
|
||||
v3D(9) = 3
|
||||
v3D(10) = 1
|
||||
|
||||
nv2 = 1
|
||||
do m = 1,neltypes
|
||||
do n = m,neltypes
|
||||
eltind(m,n) = nv2
|
||||
eltind(n,m) = nv2
|
||||
nv2 = nv2+1
|
||||
enddo
|
||||
enddo
|
||||
|
||||
c Compute pair potentials and setup arrays for interpolation
|
||||
nr = 1000
|
||||
dr = 1.1*rc_meam/nr
|
||||
call compute_pair_meam
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||
c Fill off-diagonal alloy parameters
|
||||
subroutine alloyparams
|
||||
use meam_data
|
||||
implicit none
|
||||
integer i,j,k
|
||||
real*8 eb
|
||||
|
||||
c Loop over pairs
|
||||
do i = 1,neltypes
|
||||
do j = 1,neltypes
|
||||
c Treat off-diagonal pairs
|
||||
c If i>j, set all equal to i<j case (which has aready been set,
|
||||
c here or in the input file)
|
||||
if (i.gt.j) then
|
||||
re_meam(i,j) = re_meam(j,i)
|
||||
Ec_meam(i,j) = Ec_meam(j,i)
|
||||
alpha_meam(i,j) = alpha_meam(j,i)
|
||||
lattce_meam(i,j) = lattce_meam(j,i)
|
||||
c If i<j and term is unset, use default values (e.g. mean of i-i and j-j)
|
||||
else if (j.gt.i) then
|
||||
if (Ec_meam(i,j).eq.0.d0) then
|
||||
Ec_meam(i,j) = (Ec_meam(i,i)+Ec_meam(j,j))/2.d0
|
||||
$ - delta_meam(i,j)
|
||||
endif
|
||||
if (alpha_meam(i,j).eq.0.d0) then
|
||||
alpha_meam(i,j) = (alpha_meam(i,i)+alpha_meam(j,j))/2.d0
|
||||
endif
|
||||
if (re_meam(i,j).eq.0.d0) then
|
||||
re_meam(i,j) = (re_meam(i,i)+re_meam(j,j))/2.d0
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
c Cmin(i,k,j) is symmetric in i-j, but not k. For all triplets
|
||||
c where i>j, set equal to the i<j element. Likewise for Cmax.
|
||||
do i = 2,neltypes
|
||||
do j = 1,i-1
|
||||
do k = 1,neltypes
|
||||
Cmin_meam(i,j,k) = Cmin_meam(j,i,k)
|
||||
Cmax_meam(i,j,k) = Cmax_meam(j,i,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
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 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
|
||||
do j = 1,neltypes
|
||||
do k = 1,neltypes
|
||||
eb = (Cmax_meam(i,j,k)*Cmax_meam(i,j,k))
|
||||
$ /(4.d0*(Cmax_meam(i,j,k)-1.d0))
|
||||
ebound_meam(i,j) = max(ebound_meam(i,j),eb)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
c-----------------------------------------------------------------------
|
||||
c compute MEAM pair potential for each pair of element types
|
||||
c
|
||||
|
||||
subroutine compute_pair_meam
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
real*8 r, temp
|
||||
integer j,a,b,nv2
|
||||
real*8 astar,frac,phizbl
|
||||
integer n,nmax,Z1,Z2
|
||||
real*8 arat,scrn
|
||||
real*8 phitmp
|
||||
|
||||
real*8, external :: phi_meam
|
||||
real*8, external :: zbl
|
||||
|
||||
|
||||
c allocate memory for array that defines the potential
|
||||
allocate(phir(nr,(neltypes*(neltypes+1))/2))
|
||||
|
||||
c allocate coeff memory
|
||||
|
||||
allocate(phirar(nr,(neltypes*(neltypes+1))/2))
|
||||
allocate(phirar1(nr,(neltypes*(neltypes+1))/2))
|
||||
allocate(phirar2(nr,(neltypes*(neltypes+1))/2))
|
||||
allocate(phirar3(nr,(neltypes*(neltypes+1))/2))
|
||||
allocate(phirar4(nr,(neltypes*(neltypes+1))/2))
|
||||
allocate(phirar5(nr,(neltypes*(neltypes+1))/2))
|
||||
allocate(phirar6(nr,(neltypes*(neltypes+1))/2))
|
||||
|
||||
c loop over pairs of element types
|
||||
nv2 = 0
|
||||
do a = 1,neltypes
|
||||
do b = a,neltypes
|
||||
nv2 = nv2 + 1
|
||||
|
||||
c loop over r values and compute
|
||||
do j = 1,nr
|
||||
|
||||
r = (j-1)*dr
|
||||
|
||||
phir(j,nv2) = phi_meam(r,a,b)
|
||||
|
||||
c if using second-nearest neighbor, solve recursive problem
|
||||
c (see Lee and Baskes, PRB 62(13):8564 eqn.(21))
|
||||
if (nn2_meam(a,b).eq.1) then
|
||||
if (a.ne.b) then
|
||||
write(6,*) 'Second-neighbor currently only valid '
|
||||
write(6,*) 'if element types are the same.'
|
||||
c call error(' ')
|
||||
endif
|
||||
call get_Zij(Z1,lattce_meam(a,b))
|
||||
call get_Zij2(Z2,arat,scrn,lattce_meam(a,b),
|
||||
$ Cmin_meam(a,a,b),Cmax_meam(a,a,b))
|
||||
|
||||
nmax = 10
|
||||
do n = 1,nmax
|
||||
phir(j,nv2) = phir(j,nv2) +
|
||||
$ (-Z2*scrn/Z1)**n * phi_meam(r*arat**n,a,b)
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
c if astar <= -3
|
||||
c potential is zbl potential
|
||||
c else if -3 < astar < -1
|
||||
c potential is linear combination with zbl potential
|
||||
c endif
|
||||
astar = alpha_meam(a,b) * (r/re_meam(a,b) - 1.d0)
|
||||
if (astar.le.-3.d0) then
|
||||
phir(j,nv2) = zbl(r,ielt_meam(a),ielt_meam(b))
|
||||
else if (astar.gt.-3.d0.and.astar.lt.-1.d0) then
|
||||
call fcut(1-(astar+1.d0)/(-3.d0+1.d0),frac)
|
||||
phizbl = zbl(r,ielt_meam(a),ielt_meam(b))
|
||||
phir(j,nv2) = frac*phir(j,nv2) + (1-frac)*phizbl
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
c call interpolation
|
||||
call interpolate_meam(nv2)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
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
|
||||
real*8 t11av,t21av,t31av,t12av,t22av,t32av
|
||||
real*8 G1,G2,s1(3),s2(3),s12(3),rho0_1,rho0_2
|
||||
real*8 Gam1,Gam2,Z1,Z2
|
||||
real*8 rhobar1,rhobar2,F1,F2
|
||||
real*8 rhoa01,rhoa11,rhoa21,rhoa31
|
||||
real*8 rhoa02,rhoa12,rhoa22,rhoa32
|
||||
real*8 rho01,rho11,rho21,rho31
|
||||
real*8 rho02,rho12,rho22,rho32
|
||||
real*8 scalfac,phiaa,phibb
|
||||
real*8 Eu
|
||||
integer Z12
|
||||
character*3 latta,lattb
|
||||
|
||||
real*8, external :: erose
|
||||
|
||||
c Equation numbers below refer to:
|
||||
c I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615
|
||||
|
||||
c get number of neighbors in the reference structure
|
||||
c Nref(i,j) = # of i's neighbors of type j
|
||||
call get_Zij(Z12,lattce_meam(a,b))
|
||||
|
||||
call get_densref(r,a,b,rho01,rho11,rho21,rho31,
|
||||
$ rho02,rho12,rho22,rho32)
|
||||
|
||||
c calculate average weighting factors for the reference structure
|
||||
if (lattce_meam(a,b).eq.'c11') then
|
||||
scalfac = 1.0/(rho01+rho02)
|
||||
t11av = scalfac*(t1_meam(a)*rho01 + t1_meam(b)*rho02)
|
||||
t12av = t11av
|
||||
t21av = scalfac*(t2_meam(a)*rho01 + t2_meam(b)*rho02)
|
||||
t22av = t21av
|
||||
t31av = scalfac*(t3_meam(a)*rho01 + t3_meam(b)*rho02)
|
||||
t32av = t31av
|
||||
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),
|
||||
$ t1_meam(b),t2_meam(b),t3_meam(b),lattce_meam(a,b))
|
||||
endif
|
||||
|
||||
c for c11b structure, calculate background electron densities
|
||||
if (lattce_meam(a,b).eq.'c11') then
|
||||
latta = lattce_meam(a,a)
|
||||
if (latta.eq.'dia') then
|
||||
rhobar1 = ((Z12/2)*(rho02+rho01))**2 +
|
||||
$ t11av*(rho12-rho11)**2 +
|
||||
$ t21av/6.0*(rho22+rho21)**2 +
|
||||
$ 121.0/40.*t31av*(rho32-rho31)**2
|
||||
rhobar1 = sqrt(rhobar1)
|
||||
rhobar2 = (Z12*rho01)**2 + 2.0/3.0*t21av*rho21**2
|
||||
rhobar2 = sqrt(rhobar2)
|
||||
else
|
||||
rhobar2 = ((Z12/2)*(rho01+rho02))**2 +
|
||||
$ t12av*(rho11-rho12)**2 +
|
||||
$ t22av/6.0*(rho21+rho22)**2 +
|
||||
$ 121.0/40.*t32av*(rho31-rho32)**2
|
||||
rhobar2 = sqrt(rhobar2)
|
||||
rhobar1 = (Z12*rho02)**2 + 2.0/3.0*t22av*rho22**2
|
||||
rhobar1 = sqrt(rhobar1)
|
||||
endif
|
||||
else
|
||||
c for other structures, use formalism developed in Huang's paper
|
||||
c
|
||||
c composition-dependent scaling, equation I.7
|
||||
c -- note: The shape factors for the individual elements are used, since
|
||||
c this function should cancel the F(rho) terms when the atoms are
|
||||
c in the reference structure.
|
||||
Z1 = Z_meam(a)
|
||||
Z2 = Z_meam(b)
|
||||
if (ibar_meam(a).le.0) then
|
||||
G1 = 1.d0
|
||||
else
|
||||
call get_shpfcn(s1,lattce_meam(a,a))
|
||||
Gam1 = (s1(1)*t11av+s1(2)*t21av+s1(3)*t31av)/(Z1*Z1)
|
||||
call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1)
|
||||
endif
|
||||
if (ibar_meam(b).le.0) then
|
||||
G2 = 1.d0
|
||||
else
|
||||
call get_shpfcn(s2,lattce_meam(b,b))
|
||||
Gam2 = (s2(1)*t12av+s2(2)*t22av+s2(3)*t32av)/(Z2*Z2)
|
||||
call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2)
|
||||
endif
|
||||
rho0_1 = rho0_meam(a)*Z1*G1
|
||||
rho0_2 = rho0_meam(b)*Z2*G2
|
||||
|
||||
c get number of neighbors in the reference structure
|
||||
c Nref(i,j) = # of i's neighbors of type j
|
||||
c call get_Zij(Z12,lattce_meam(a,b))
|
||||
c
|
||||
c call get_densref(r,a,b,rho01,rho11,rho21,rho31,
|
||||
c $ rho02,rho12,rho22,rho32)
|
||||
|
||||
c compute total background density, eqn I.6
|
||||
Gam1 = (t11av*rho11+t21av*rho21+t31av*rho31)
|
||||
Gam1 = Gam1/(rho01*rho01)
|
||||
Gam2 = (t12av*rho12+t22av*rho22+t32av*rho32)
|
||||
Gam2 = Gam2/(rho02*rho02)
|
||||
call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1)
|
||||
call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2)
|
||||
rhobar1 = rho01/rho0_1*G1
|
||||
rhobar2 = rho02/rho0_2*G2
|
||||
endif
|
||||
c compute embedding functions, eqn I.5
|
||||
if (rhobar1.eq.0.d0) then
|
||||
F1 = 0.d0
|
||||
else
|
||||
F1 = A_meam(a)*Ec_meam(a,a)*rhobar1*log(rhobar1)
|
||||
endif
|
||||
if (rhobar2.eq.0.d0) then
|
||||
F2 = 0.d0
|
||||
else
|
||||
F2 = A_meam(b)*Ec_meam(b,b)*rhobar2*log(rhobar2)
|
||||
endif
|
||||
|
||||
c compute Rose function, I.16
|
||||
Eu = erose(r,re_meam(a,b),alpha_meam(a,b),
|
||||
$ Ec_meam(a,b),repuls_meam(a,b),attrac_meam(a,b))
|
||||
|
||||
c calculate the pair energy
|
||||
if (lattce_meam(a,b).eq.'c11') then
|
||||
latta = lattce_meam(a,a)
|
||||
if (latta.eq.'dia') then
|
||||
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
|
||||
endif
|
||||
c phi_m = 0.d0
|
||||
else
|
||||
c
|
||||
c potential is computed from Rose function and embedding energy
|
||||
phi_m = (2*Eu - F1 - F2)/Z12
|
||||
c
|
||||
endif
|
||||
|
||||
c if r = 0, just return 0
|
||||
if (r.eq.0.d0) then
|
||||
phi_m = 0.d0
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
c----------------------------------------------------------------------c
|
||||
c Shape factors for various configurations
|
||||
subroutine get_shpfcn(s,latt)
|
||||
implicit none
|
||||
real*8 s(3)
|
||||
character*3 latt
|
||||
if (latt.eq.'fcc'.or.latt.eq.'bcc'.or.latt.eq.'b1') then
|
||||
s(1) = 0.d0
|
||||
s(2) = 0.d0
|
||||
s(3) = 0.d0
|
||||
else if (latt.eq.'hcp') then
|
||||
s(1) = 0.d0
|
||||
s(2) = 0.d0
|
||||
s(3) = 1.d0/3.d0
|
||||
else if (latt.eq.'dia') then
|
||||
s(1) = 0.d0
|
||||
s(2) = 0.d0
|
||||
s(3) = 32.d0/9.d0
|
||||
else if (latt.eq.'dim') then
|
||||
s(1) = 1.d0
|
||||
s(2) = 2.d0/3.d0
|
||||
c s(3) = 1.d0
|
||||
s(3) = 0.4d0
|
||||
else
|
||||
s(1) = 0.0
|
||||
c call error('Lattice not defined in get_shpfcn.')
|
||||
endif
|
||||
return
|
||||
end
|
||||
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,latt)
|
||||
implicit none
|
||||
real*8 t11av,t21av,t31av,t12av,t22av,t32av
|
||||
real*8 t11,t21,t31,t12,t22,t32
|
||||
character*3 latt
|
||||
if (latt.eq.'fcc'.or.latt.eq.'bcc'.or.latt.eq.'dia'
|
||||
$ .or.latt.eq.'hcp'.or.latt.eq.'b1'
|
||||
$ .or.latt.eq.'dim') then
|
||||
c all neighbors are of the opposite type
|
||||
t11av = t12
|
||||
t21av = t22
|
||||
t31av = t32
|
||||
t12av = t11
|
||||
t22av = t21
|
||||
t32av = t31
|
||||
else
|
||||
c call error('Lattice not defined in get_tavref.')
|
||||
endif
|
||||
return
|
||||
end
|
||||
c------------------------------------------------------------------------------c
|
||||
c Number of neighbors for the reference structure
|
||||
subroutine get_Zij(Zij,latt)
|
||||
implicit none
|
||||
integer Zij
|
||||
character*3 latt
|
||||
if (latt.eq.'fcc') then
|
||||
Zij = 12
|
||||
else if (latt.eq.'bcc') then
|
||||
Zij = 8
|
||||
else if (latt.eq.'hcp') then
|
||||
Zij = 12
|
||||
else if (latt.eq.'b1') then
|
||||
Zij = 6
|
||||
else if (latt.eq.'dia') then
|
||||
Zij = 4
|
||||
else if (latt.eq.'dim') then
|
||||
Zij = 1
|
||||
else if (latt.eq.'c11') then
|
||||
Zij = 10
|
||||
else
|
||||
c call error('Lattice not defined in get_Zij.')
|
||||
endif
|
||||
return
|
||||
end
|
||||
|
||||
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)
|
||||
implicit none
|
||||
integer Zij2
|
||||
real*8 a,S,cmin,cmax
|
||||
character*3 latt
|
||||
real*8 rratio,C,x,sijk
|
||||
integer numscr
|
||||
|
||||
if (latt.eq.'bcc') then
|
||||
Zij2 = 6
|
||||
a = 2.d0/sqrt(3.d0)
|
||||
numscr = 4
|
||||
else if (latt.eq.'fcc') then
|
||||
Zij2 = 6
|
||||
a = sqrt(2.d0)
|
||||
numscr = 4
|
||||
else if (latt.eq.'dia') then
|
||||
Zij2 = 0
|
||||
a = sqrt(8.d0/3.d0)
|
||||
numscr = 4
|
||||
if (cmin.lt.0.500001) then
|
||||
c call error('can not do 2NN MEAM for dia')
|
||||
endif
|
||||
else if (latt.eq.'hcp') then
|
||||
Zij2 = 6
|
||||
a = sqrt(2.d0)
|
||||
numscr = 4
|
||||
else if (latt.eq.'b1') then
|
||||
Zij2 = 6
|
||||
a = sqrt(2.d0)
|
||||
numscr = 4
|
||||
else
|
||||
c call error('Lattice not defined in get_Zij2.')
|
||||
endif
|
||||
|
||||
c Compute screening for each first neighbor
|
||||
C = 4.d0/(a*a) - 1.d0
|
||||
x = (C-cmin)/(cmax-cmin)
|
||||
call fcut(x,sijk)
|
||||
c There are numscr first neighbors screening the second neighbors
|
||||
S = sijk**numscr
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
c------------------------------------------------------------------------------c
|
||||
c Calculate density functions, assuming reference configuration
|
||||
subroutine get_densref(r,a,b,rho01,rho11,rho21,rho31,
|
||||
$ rho02,rho12,rho22,rho32)
|
||||
use meam_data
|
||||
implicit none
|
||||
real*8 r,rho01,rho11,rho21,rho31,rho02,rho12,rho22,rho32
|
||||
real*8 a1,a2
|
||||
real*8 rhoa01,rhoa11,rhoa21,rhoa31,rhoa02,rhoa12,rhoa22,rhoa32
|
||||
real*8 C,Cmin,Cmax,fc,s(3)
|
||||
character*3 lat
|
||||
integer a,b
|
||||
integer Zij1nn,Zij2nn
|
||||
real*8 rhoa01nn,rhoa02nn
|
||||
real*8 arat,scrn
|
||||
|
||||
a1 = r/re_meam(a,a) - 1.d0
|
||||
a2 = r/re_meam(b,b) - 1.d0
|
||||
|
||||
rhoa01 = rho0_meam(a)*exp(-beta0_meam(a)*a1)
|
||||
rhoa11 = rho0_meam(a)*exp(-beta1_meam(a)*a1)
|
||||
rhoa21 = rho0_meam(a)*exp(-beta2_meam(a)*a1)
|
||||
rhoa31 = rho0_meam(a)*exp(-beta3_meam(a)*a1)
|
||||
rhoa02 = rho0_meam(b)*exp(-beta0_meam(b)*a2)
|
||||
rhoa12 = rho0_meam(b)*exp(-beta1_meam(b)*a2)
|
||||
rhoa22 = rho0_meam(b)*exp(-beta2_meam(b)*a2)
|
||||
rhoa32 = rho0_meam(b)*exp(-beta3_meam(b)*a2)
|
||||
|
||||
lat = lattce_meam(a,b)
|
||||
|
||||
rho11 = 0.d0
|
||||
rho21 = 0.d0
|
||||
rho31 = 0.d0
|
||||
rho12 = 0.d0
|
||||
rho22 = 0.d0
|
||||
rho32 = 0.d0
|
||||
|
||||
call get_Zij(Zij1nn,lat)
|
||||
|
||||
if (lat.eq.'fcc') then
|
||||
rho01 = 12.d0*rhoa02
|
||||
rho02 = 12.d0*rhoa01
|
||||
else if (lat.eq.'bcc') then
|
||||
rho01 = 8.d0*rhoa02
|
||||
rho02 = 8.d0*rhoa01
|
||||
else if (lat.eq.'b1') then
|
||||
c Warning: this assumes that atoms of element 1 are not completely
|
||||
c screened from each other, but that atoms of element 2 are. So
|
||||
c Cmin_meam(1,1,2) can be less than 1.0, but Cmin_meam(2,2,1) cannot.
|
||||
c This should be made more general in the future.
|
||||
c Cmin = Cmin_meam(1,1,2)
|
||||
c Cmax = Cmax_meam(1,1,2)
|
||||
c C = (1.0-Cmin)/(Cmax-Cmin)
|
||||
c call fcut(C,fc)
|
||||
c a1 = sqrt(2.0)*r/re_meam(1,1) - 1.d0
|
||||
c rho01 = 6*rhoa02 + 12*fc*fc*rho0_meam(1)*exp(-beta0_meam(1)*a1)
|
||||
c rho02 = 6*rhoa01
|
||||
rho01 = 6*rhoa02
|
||||
rho02 = 6*rhoa01
|
||||
else if (lat.eq.'dia') then
|
||||
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
|
||||
rho31 = 1.d0/3.d0*rhoa32*rhoa32
|
||||
rho32 = 1.d0/3.d0*rhoa31*rhoa31
|
||||
else if (lat.eq.'dim') then
|
||||
call 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.eq.'c11') then
|
||||
rho01 = rhoa01
|
||||
rho02 = rhoa02
|
||||
rho11 = rhoa11
|
||||
rho12 = rhoa12
|
||||
rho21 = rhoa21
|
||||
rho22 = rhoa22
|
||||
rho31 = rhoa31
|
||||
rho32 = rhoa32
|
||||
else
|
||||
c call error('Lattice not defined in get_densref.')
|
||||
endif
|
||||
|
||||
if (nn2_meam(a,b).eq.1) then
|
||||
c Assume that second neighbor is of same type, first neighbor
|
||||
c may be of different type
|
||||
|
||||
if (lat.ne.'bcc') then
|
||||
c call error('Second-neighbor not defined for this lattice')
|
||||
endif
|
||||
|
||||
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)*exp(-beta0_meam(a)*a1)
|
||||
rhoa02nn = rho0_meam(b)*exp(-beta0_meam(b)*a2)
|
||||
|
||||
rho01 = rho01 + Zij2nn*scrn*rhoa02nn
|
||||
rho02 = rho02 + Zij2nn*scrn*rhoa01nn
|
||||
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
c---------------------------------------------------------------------
|
||||
c Compute ZBL potential
|
||||
c
|
||||
real*8 function zbl(r,z1,z2)
|
||||
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)*exp(-d(i)*x)
|
||||
enddo
|
||||
if (r.gt.0.d0) zbl = zbl*z1*z2/r*cc
|
||||
return
|
||||
end
|
||||
|
||||
c---------------------------------------------------------------------
|
||||
c Compute Rose energy function, I.16
|
||||
c
|
||||
real*8 function erose(r,re,alpha,Ec,repuls,attrac)
|
||||
implicit none
|
||||
real*8 r,re,alpha,Ec,repuls,attrac,astar,a3
|
||||
|
||||
erose = 0.d0
|
||||
|
||||
if (r.gt.0.d0) then
|
||||
astar = alpha * (r/re - 1.d0)
|
||||
a3 = 0.d0
|
||||
if (astar.ge.0) then
|
||||
a3 = attrac
|
||||
else if (astar.lt.0) then
|
||||
a3 = repuls
|
||||
endif
|
||||
c erose = -Ec * (1 + astar + a3*(astar**3)/(r/re)) * exp(-astar)
|
||||
erose = -Ec*(1+astar+(-attrac+repuls/r)*(astar**3))*exp(-astar)
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
c -----------------------------------------------------------------------
|
||||
|
||||
subroutine interpolate_meam(ind)
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer j,ind
|
||||
real*8 drar
|
||||
|
||||
c map to coefficient space
|
||||
|
||||
nrar = nr
|
||||
drar = dr
|
||||
rdrar = 1.0D0/drar
|
||||
|
||||
c phir interp
|
||||
do j = 1,nrar
|
||||
phirar(j,ind) = phir(j,ind)
|
||||
enddo
|
||||
|
||||
phirar1(1,ind) = phirar(2,ind)-phirar(1,ind)
|
||||
phirar1(2,ind) = 0.5D0*(phirar(3,ind)-phirar(1,ind))
|
||||
phirar1(nrar-1,ind) = 0.5D0*(phirar(nrar,ind)
|
||||
$ -phirar(nrar-2,ind))
|
||||
phirar1(nrar,ind) = 0.0D0
|
||||
do j = 3,nrar-2
|
||||
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)
|
||||
phirar3(j,ind) = phirar1(j,ind) + phirar1(j+1,ind) -
|
||||
$ 2.0D0*(phirar(j+1,ind)-phirar(j,ind))
|
||||
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
|
||||
phirar6(j,ind) = 3.0D0*phirar3(j,ind)/drar
|
||||
enddo
|
||||
|
||||
end
|
||||
105
lib/meam/meam_setup_global.F
Executable file
105
lib/meam/meam_setup_global.F
Executable file
@ -0,0 +1,105 @@
|
||||
c
|
||||
c declaration in pair_meam.h:
|
||||
c
|
||||
c void meam_setup_global(int *, int *, double *, int *, double *, double *,
|
||||
c double *, double *, double *, double *, double *,
|
||||
c double *, double *, double *, double *, double *,
|
||||
c double *, double *, int *);
|
||||
c
|
||||
c call in pair_meam.cpp:
|
||||
c
|
||||
c meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
|
||||
c alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
|
||||
c
|
||||
c
|
||||
|
||||
subroutine meam_setup_global(nelt, lat, z, ielement, atwt, alpha,
|
||||
$ b0, b1, b2, b3, alat, esub, asub,
|
||||
$ t0, t1, t2, t3, rozero, ibar)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer nelt, lat, ielement, ibar
|
||||
real*8 z, atwt, alpha, b0, b1, b2, b3
|
||||
real*8 alat, esub, asub, t0, t1, t2, t3
|
||||
real*8 rozero
|
||||
|
||||
dimension lat(nelt), ielement(nelt), ibar(nelt)
|
||||
dimension z(nelt), atwt(nelt), alpha(nelt)
|
||||
dimension b0(nelt), b1(nelt), b2(nelt), b3(nelt)
|
||||
dimension alat(nelt), esub(nelt), asub(nelt)
|
||||
dimension t0(nelt), t1(nelt), t2(nelt), t3(nelt), rozero(nelt)
|
||||
|
||||
integer i
|
||||
real*8 tmplat(maxelt)
|
||||
|
||||
neltypes = nelt
|
||||
|
||||
do i = 1,nelt
|
||||
|
||||
if (lat(i).eq.0) then
|
||||
lattce_meam(i,i) = 'fcc'
|
||||
else if (lat(i).eq.1) then
|
||||
lattce_meam(i,i) = 'bcc'
|
||||
else if (lat(i).eq.2) then
|
||||
lattce_meam(i,i) = 'hcp'
|
||||
else if (lat(i).eq.3) then
|
||||
lattce_meam(i,i) = 'dim'
|
||||
else if (lat(i).eq.4) then
|
||||
lattce_meam(i,i) = 'dia'
|
||||
else
|
||||
c unknown
|
||||
endif
|
||||
|
||||
Z_meam(i) = z(i)
|
||||
ielt_meam(i) = ielement(i)
|
||||
alpha_meam(i,i) = alpha(i)
|
||||
beta0_meam(i) = b0(i)
|
||||
beta1_meam(i) = b1(i)
|
||||
beta2_meam(i) = b2(i)
|
||||
beta3_meam(i) = b3(i)
|
||||
tmplat(i) = alat(i)
|
||||
Ec_meam(i,i) = esub(i)
|
||||
A_meam(i) = asub(i)
|
||||
t0_meam(i) = t0(i)
|
||||
t1_meam(i) = t1(i)
|
||||
t2_meam(i) = t2(i)
|
||||
t3_meam(i) = t3(i)
|
||||
rho0_meam(i) = rozero(i)
|
||||
ibar_meam(i) = ibar(i)
|
||||
|
||||
if (lattce_meam(i,i).eq.'fcc') then
|
||||
re_meam(i,i) = tmplat(i)/sqrt(2.d0)
|
||||
elseif (lattce_meam(i,i).eq.'bcc') then
|
||||
re_meam(i,i) = tmplat(i)*sqrt(3.d0)/2.d0
|
||||
elseif (lattce_meam(i,i).eq.'hcp') then
|
||||
re_meam(i,i) = tmplat(i)
|
||||
elseif (lattce_meam(i,i).eq.'dim') then
|
||||
re_meam(i,i) = tmplat(i)
|
||||
elseif (lattce_meam(i,i).eq.'dia') then
|
||||
re_meam(i,i) = tmplat(i)*sqrt(3.d0)/4.d0
|
||||
else
|
||||
c error
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
c Set some defaults
|
||||
rc_meam = 4.0
|
||||
delr_meam = 0.1
|
||||
attrac_meam(:,:) = 0.0
|
||||
repuls_meam(:,:) = 0.0
|
||||
Cmax_meam(:,:,:) = 2.8
|
||||
Cmin_meam(:,:,:) = 2.0
|
||||
ebound_meam(:,:) = (2.8d0**2)/(4.d0*(2.8d0-1.d0))
|
||||
delta_meam(:,:) = 0.0
|
||||
nn2_meam(:,:) = 0
|
||||
gsmooth_factor = 99.0
|
||||
augt1 = 1
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
121
lib/meam/meam_setup_param.F
Executable file
121
lib/meam/meam_setup_param.F
Executable file
@ -0,0 +1,121 @@
|
||||
c
|
||||
c Declaration in pair_meam.h:
|
||||
c
|
||||
c void meam_setup_param(int *, double *, int *, int *, int *);
|
||||
c
|
||||
c Call in pair_meam.cpp
|
||||
c
|
||||
c meam_setup_param(&which,&value,&nindex,index,&errorflag);
|
||||
c
|
||||
c
|
||||
c
|
||||
c The "which" argument corresponds to the index of the "keyword" array
|
||||
c in pair_meam.cpp:
|
||||
c
|
||||
c 0 = Ec_meam
|
||||
c 1 = alpha_meam
|
||||
c 2 = rho0_meam
|
||||
c 3 = delta_meam
|
||||
c 4 = lattce_meam
|
||||
c 5 = attrac_meam
|
||||
c 6 = repuls_meam
|
||||
c 7 = nn2_meam
|
||||
c 8 = Cmin_meam
|
||||
c 9 = Cmax_meam
|
||||
c 10 = rc_meam
|
||||
c 11 = delr_meam
|
||||
c 12 = augt1
|
||||
c 13 = gsmooth_factor
|
||||
c 14 = re_meam
|
||||
|
||||
subroutine meam_setup_param(which, value, nindex,
|
||||
$ index, errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
||||
integer which, nindex, index(3), errorflag
|
||||
real*8 value
|
||||
|
||||
errorflag = 0
|
||||
|
||||
c 0 = Ec_meam
|
||||
if (which.eq.0) then
|
||||
Ec_meam(index(1),index(2)) = value
|
||||
|
||||
c 1 = alpha_meam
|
||||
else if (which.eq.1) then
|
||||
alpha_meam(index(1),index(2)) = value
|
||||
|
||||
c 2 = rho0_meam
|
||||
else if (which.eq.2) then
|
||||
rho0_meam(index(1)) = value
|
||||
|
||||
c 3 = delta_meam
|
||||
else if (which.eq.3) then
|
||||
delta_meam(index(1),index(2)) = value
|
||||
|
||||
c 4 = lattce_meam
|
||||
else if (which.eq.4) then
|
||||
if (value.eq.0) then
|
||||
lattce_meam(index(1),index(2)) = "fcc"
|
||||
else if (value.eq.1) then
|
||||
lattce_meam(index(1),index(2)) = "bcc"
|
||||
else if (value.eq.2) then
|
||||
lattce_meam(index(1),index(2)) = "hcp"
|
||||
else if (value.eq.3) then
|
||||
lattce_meam(index(1),index(2)) = "dim"
|
||||
else if (value.eq.4) then
|
||||
lattce_meam(index(1),index(2)) = "dia"
|
||||
else if (value.eq.5) then
|
||||
lattce_meam(index(1),index(2)) = 'b1'
|
||||
else if (value.eq.6) then
|
||||
lattce_meam(index(1),index(2)) = 'c11'
|
||||
endif
|
||||
|
||||
c 5 = attrac_meam
|
||||
else if (which.eq.5) then
|
||||
attrac_meam(index(1),index(2)) = value
|
||||
|
||||
c 6 = repuls_meam
|
||||
else if (which.eq.6) then
|
||||
repuls_meam(index(1),index(2)) = value
|
||||
|
||||
c 7 = nn2_meam
|
||||
else if (which.eq.7) then
|
||||
nn2_meam(index(1),index(2)) = value
|
||||
|
||||
c 8 = Cmin_meam
|
||||
else if (which.eq.8) then
|
||||
Cmin_meam(index(1),index(2),index(3)) = value
|
||||
|
||||
c 9 = Cmax_meam
|
||||
else if (which.eq.9) then
|
||||
Cmax_meam(index(1),index(2),index(3)) = value
|
||||
|
||||
c 10 = rc_meam
|
||||
else if (which.eq.10) then
|
||||
rc_meam = value
|
||||
|
||||
c 11 = delr_meam
|
||||
else if (which.eq.11) then
|
||||
delr_meam = value
|
||||
|
||||
c 12 = augt1
|
||||
else if (which.eq.12) then
|
||||
augt1 = value
|
||||
|
||||
c 13 = gsmooth
|
||||
else if (which.eq.13) then
|
||||
gsmooth_factor = value
|
||||
|
||||
c 14 = re_meam
|
||||
else if (which.eq.14) then
|
||||
re_meam(index(1),index(2)) = value
|
||||
|
||||
else
|
||||
errorflag = 1
|
||||
endif
|
||||
|
||||
return
|
||||
end
|
||||
Reference in New Issue
Block a user