jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2016-04-15 16:16:52 +00:00
parent e6ca2d5e08
commit a4e8eaaf4d
6 changed files with 46 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
lib/meam/Makefile.gfortran
View File

@ -22,12 +22,13 @@ FILES = $(SRC) Makefile
# ------ DEFINITIONS ------
LIB = libmeam.a
OBJ = $(SRC:.F=.o)
OBJ = $(SRC:.F=.o) fm_exp.o
# ------ SETTINGS ------
F90 = gfortran
F90FLAGS = -O -fPIC -fno-second-underscore
CC = gcc
F90FLAGS = -O3 -fPIC -ffast-math -ftree-vectorize -fexpensive-optimizations -fno-second-underscore -g
#F90FLAGS = -O
ARCHIVE = ar
ARCHFLAG = -rc

1
lib/meam/meam_data.F
View File

@ -2,6 +2,7 @@
module meam_data
integer, parameter :: maxelt = 5
real*8 , external :: fm_exp
c cutforce = force cutoff
c cutforcesq = force cutoff squared

12
lib/meam/meam_dens_final.F
View File

@ -202,6 +202,7 @@ c 2 => not implemented
c 3 => G = 2/(1+exp(-Gamma))
c 4 => G = sqrt(1+Gamma)
c -5 => G = +-sqrt(abs(1+Gamma))
use meam_data , only: fm_exp
implicit none
real*8 Gamma,G
real*8 gsmooth_factor, gsmooth_switchpoint
@ -219,7 +220,7 @@ c G = 0.01*(-0.99/Gamma)**99
G = sqrt(1.d0+Gamma)
endif
else if (ibar.eq.1) then
G = exp(Gamma/2.d0)
G = fm_exp(Gamma/2.d0)
else if (ibar.eq.3) then
G = 2.d0/(1.d0+exp(-Gamma))
else if (ibar.eq.-5) then
@ -239,11 +240,12 @@ 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 1 => G = fm_exp(Gamma/2)
c 2 => not implemented
c 3 => G = 2/(1+exp(-Gamma))
c 3 => G = 2/(1+fm_exp(-Gamma))
c 4 => G = sqrt(1+Gamma)
c -5 => G = +-sqrt(abs(1+Gamma))
use meam_data , only: fm_exp
real*8 Gamma,G,dG
real*8 gsmooth_factor, gsmooth_switchpoint
integer ibar
@ -262,10 +264,10 @@ c G = 0.01*(-0.99/Gamma)**99
dG = 1.d0/(2.d0*G)
endif
else if (ibar.eq.1) then
G = exp(Gamma/2.d0)
G = fm_exp(Gamma/2.d0)
dG = G/2.d0
else if (ibar.eq.3) then
G = 2.d0/(1.d0+exp(-Gamma))
G = 2.d0/(1.d0+fm_exp(-Gamma))
dG = G*(2.d0-G)/2
else if (ibar.eq.-5) then
if ((1.d0+Gamma).ge.0) then

16
lib/meam/meam_dens_init.F
View File

@ -245,14 +245,14 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
rhoa0j = ro0j*fm_exp(-beta0_meam(eltj)*aj)*sij
rhoa1j = ro0j*fm_exp(-beta1_meam(eltj)*aj)*sij
rhoa2j = ro0j*fm_exp(-beta2_meam(eltj)*aj)*sij
rhoa3j = ro0j*fm_exp(-beta3_meam(eltj)*aj)*sij
rhoa0i = ro0i*fm_exp(-beta0_meam(elti)*ai)*sij
rhoa1i = ro0i*fm_exp(-beta1_meam(elti)*ai)*sij
rhoa2i = ro0i*fm_exp(-beta2_meam(elti)*ai)*sij
rhoa3i = ro0i*fm_exp(-beta3_meam(elti)*ai)*sij
if (ialloy.eq.1) then
rhoa1j = rhoa1j * t1_meam(eltj)
rhoa2j = rhoa2j * t2_meam(eltj)

16
lib/meam/meam_force.F
View File

@ -144,26 +144,26 @@ 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)
rhoa0i = ro0i*fm_exp(-beta0_meam(elti)*ai)
drhoa0i = -beta0_meam(elti)*invrei*rhoa0i
rhoa1i = ro0i*exp(-beta1_meam(elti)*ai)
rhoa1i = ro0i*fm_exp(-beta1_meam(elti)*ai)
drhoa1i = -beta1_meam(elti)*invrei*rhoa1i
rhoa2i = ro0i*exp(-beta2_meam(elti)*ai)
rhoa2i = ro0i*fm_exp(-beta2_meam(elti)*ai)
drhoa2i = -beta2_meam(elti)*invrei*rhoa2i
rhoa3i = ro0i*exp(-beta3_meam(elti)*ai)
rhoa3i = ro0i*fm_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)
rhoa0j = ro0j*fm_exp(-beta0_meam(eltj)*aj)
drhoa0j = -beta0_meam(eltj)*invrej*rhoa0j
rhoa1j = ro0j*exp(-beta1_meam(eltj)*aj)
rhoa1j = ro0j*fm_exp(-beta1_meam(eltj)*aj)
drhoa1j = -beta1_meam(eltj)*invrej*rhoa1j
rhoa2j = ro0j*exp(-beta2_meam(eltj)*aj)
rhoa2j = ro0j*fm_exp(-beta2_meam(eltj)*aj)
drhoa2j = -beta2_meam(eltj)*invrej*rhoa2j
rhoa3j = ro0j*exp(-beta3_meam(eltj)*aj)
rhoa3j = ro0j*fm_exp(-beta3_meam(eltj)*aj)
drhoa3j = -beta3_meam(eltj)*invrej*rhoa3j
else
rhoa0j = rhoa0i

36
lib/meam/meam_setup_done.F
View File

@ -558,7 +558,7 @@ c screening)
if (nn2_meam(a,a).eq.1) then
call get_Zij2(Z2,arat,scrn,lattce_meam(a,a),
$ Cmin_meam(a,a,a),Cmax_meam(a,a,a))
rho0_2nn = rho0_meam(a)*exp(-beta0_meam(a)*(arat-1))
rho0_2nn = rho0_meam(a)*fm_exp(-beta0_meam(a)*(arat-1))
rho0 = rho0 + Z2*rho0_2nn*scrn
endif
@ -634,8 +634,8 @@ c all neighbors are of the opposite type
else
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)
rhoa02 = rho0_meam(b)*exp(-beta0_meam(b)*a2)
rhoa01 = rho0_meam(a)*fm_exp(-beta0_meam(a)*a1)
rhoa02 = rho0_meam(b)*fm_exp(-beta0_meam(b)*a2)
if (latt.eq.'l12') then
rho01 = 8*rhoa01 + 4*rhoa02
t11av = (8*t11*rhoa01 + 4*t12*rhoa02)/rho01
@ -777,14 +777,14 @@ c Calculate density functions, assuming reference configuration
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)
rhoa01 = rho0_meam(a)*fm_exp(-beta0_meam(a)*a1)
rhoa11 = rho0_meam(a)*fm_exp(-beta1_meam(a)*a1)
rhoa21 = rho0_meam(a)*fm_exp(-beta2_meam(a)*a1)
rhoa31 = rho0_meam(a)*fm_exp(-beta3_meam(a)*a1)
rhoa02 = rho0_meam(b)*fm_exp(-beta0_meam(b)*a2)
rhoa12 = rho0_meam(b)*fm_exp(-beta1_meam(b)*a2)
rhoa22 = rho0_meam(b)*fm_exp(-beta2_meam(b)*a2)
rhoa32 = rho0_meam(b)*fm_exp(-beta3_meam(b)*a2)
lat = lattce_meam(a,b)
@ -862,8 +862,8 @@ c call error('Lattice not defined in get_densref.')
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)
rhoa01nn = rho0_meam(a)*fm_exp(-beta0_meam(a)*a1)
rhoa02nn = rho0_meam(b)*fm_exp(-beta0_meam(b)*a2)
if (lat.eq.'l12') then
c As usual, L12 thinks it's special; we need to be careful computing
@ -899,6 +899,7 @@ c---------------------------------------------------------------------
c Compute ZBL potential
c
real*8 function zbl(r,z1,z2)
use meam_data , only : fm_exp
implicit none
integer i,z1,z2
real*8 r,c,d,a,azero,cc,x
@ -912,7 +913,7 @@ c azero = (9pi^2/128)^1/3 (0.529) Angstroms
zbl = 0.0
x = r/a
do i=1,4
zbl = zbl + c(i)*exp(-d(i)*x)
zbl = zbl + c(i)*fm_exp(-d(i)*x)
enddo
if (r.gt.0.d0) zbl = zbl*z1*z2/r*cc
return
@ -922,6 +923,7 @@ c---------------------------------------------------------------------
c Compute Rose energy function, I.16
c
real*8 function erose(r,re,alpha,Ec,repuls,attrac,form)
use meam_data , only : fm_exp
implicit none
real*8 r,re,alpha,Ec,repuls,attrac,astar,a3
integer form
@ -938,11 +940,11 @@ c
endif
if (form.eq.1) then
erose = -Ec*(1+astar+(-attrac+repuls/r)*
$ (astar**3))*exp(-astar)
$ (astar**3))*fm_exp(-astar)
else if (form.eq.2) then
erose = -Ec * (1 +astar + a3*(astar**3))*exp(-astar)
erose = -Ec * (1 +astar + a3*(astar**3))*fm_exp(-astar)
else
erose = -Ec * (1 + astar + a3*(astar**3)/(r/re)) * exp(-astar)
erose = -Ec * (1+ astar + a3*(astar**3)/(r/re))*fm_exp(-astar)
endif
endif