1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
|
C author: X. W. Zhou, xzhou@sandia.gov
C updates by: Lucas Hale lucas.hale@nist.gov
c open(unit=5,file='a.i')
call inter
c close(5)
call writeset
stop
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c main subroutine. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine inter
implicit real*8 (a-h,o-z)
implicit integer (i-m)
character*80 atomtype,atommatch,outfile,outelem
namelist /funccard/ atomtype
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
* beta(16),beta1(16),A(16),B(16),cai(16),ramda(16),
* ramda1(16),Fi0(16),Fi1(16),Fi2(16),Fi3(16),
* Fm0(16),Fm1(16),Fm2(16),Fm3(16),Fm4(16),
* fnn(16),Fn(16),rhoin(16),rhoout(16),rhol(16),
* rhoh(16),rhos(16)
common /pass2/ amass(16),Fr(5000,16),rhor(5000,16),
* z2r(5000,16,16),blat(16),drho,dr,rc,outfile,outelem
common /pass3/ ielement(16),ntypes,nrho,nr
ntypes=0
10 continue
atomtype='none'
read(5,funccard)
if (atomtype .eq. 'none') goto 1200
open(unit=10,file='EAM_code',form='FORMATTED',status='OLD')
11 read(10,9501,end=1210)atommatch
9501 format(a80)
if (atomtype .eq. atommatch) then
ntypes=ntypes+1
length=len_trim(outfile)
if (length .eq. len(outfile)) then
outfile = atomtype
else
outfile = outfile(1:length)//atomtype
endif
length=len_trim(outelem)
if (length .eq. len(outelem)) then
outelem = atomtype
else
outelem = outelem(1:length)//' '//atomtype
endif
read(10,*) re(ntypes)
read(10,*) fe(ntypes)
read(10,*) rhoe(ntypes)
read(10,*) rhos(ntypes)
read(10,*) alpha(ntypes)
read(10,*) beta(ntypes)
read(10,*) A(ntypes)
read(10,*) B(ntypes)
read(10,*) cai(ntypes)
read(10,*) ramda(ntypes)
read(10,*) Fi0(ntypes)
read(10,*) Fi1(ntypes)
read(10,*) Fi2(ntypes)
read(10,*) Fi3(ntypes)
read(10,*) Fm0(ntypes)
read(10,*) Fm1(ntypes)
read(10,*) Fm2(ntypes)
read(10,*) Fm3(ntypes)
read(10,*) fnn(ntypes)
read(10,*) Fn(ntypes)
read(10,*) ielement(ntypes)
read(10,*) amass(ntypes)
read(10,*) Fm4(ntypes)
read(10,*) beta1(ntypes)
read(10,*) ramda1(ntypes)
read(10,*) rhol(ntypes)
read(10,*) rhoh(ntypes)
blat(ntypes)=sqrt(2.0_8)*re(ntypes)
rhoin(ntypes)=rhol(ntypes)*rhoe(ntypes)
rhoout(ntypes)=rhoh(ntypes)*rhoe(ntypes)
else
do i=1,27
read(10,*)vtmp
end do
goto 11
endif
close(10)
goto 10
1210 write(6,*)'error: atom type ',atomtype,' not found'
stop
1200 continue
nr=2000
nrho=2000
alatmax=blat(1)
rhoemax=rhoe(1)
do 2 i=2,ntypes
if (alatmax .lt. blat(i)) alatmax=blat(i)
if (rhoemax .lt. rhoe(i)) rhoemax=rhoe(i)
2 continue
rc=sqrt(10.0_8)/2.0_8*alatmax
rst=0.5_8
dr=rc/(nr-1.0_8)
fmax=-1.0_8
do i1=1,ntypes
do i2=1,i1
if ( i1 .eq. i2) then
do i=1,nr
r=(i-1.0_8)*dr
if (r .lt. rst) r=rst
call prof(i1,r,fvalue)
if (fmax .lt. fvalue) fmax=fvalue
rhor(i,i1)=fvalue
call pair(i1,i2,r,psi)
z2r(i,i1,i2)=r*psi
end do
else
do i=1,nr
r=(i-1.0_8)*dr
if (r .lt. rst) r=rst
call pair(i1,i2,r,psi)
z2r(i,i1,i2)=r*psi
z2r(i,i2,i1)=z2r(i,i1,i2)
end do
endif
end do
end do
rhom=fmax
if (rhom .lt. 2.0_8*rhoemax) rhom=2.0_8*rhoemax
if (rhom .lt. 100.0_8) rhom=100.0_8
drho=rhom/(nrho-1.0_8)
do 6 it=1,ntypes
do 7 i=1,nrho
rhoF=(i-1)*drho
if (i .eq. 1) rhoF=0.0_8
call embed(it,rhoF,emb)
Fr(i,it)=emb
7 continue
6 continue
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c This subroutine calculates the electron density. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine prof(it,r,f)
implicit real*8 (a-h,o-z)
implicit integer (i-m)
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
* beta(16),beta1(16),A(16),B(16),cai(16),ramda(16),
* ramda1(16),Fi0(16),Fi1(16),Fi2(16),Fi3(16),
* Fm0(16),Fm1(16),Fm2(16),Fm3(16),Fm4(16),
* fnn(16),Fn(16),rhoin(16),rhoout(16),rhol(16),
* rhoh(16),rhos(16)
f=fe(it)*exp(-beta1(it)*(r/re(it)-1.0_8))
f=f/(1.0_8+(r/re(it)-ramda1(it))**20)
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c This subroutine calculates the pair potential. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine pair(it1,it2,r,psi)
implicit real*8 (a-h,o-z)
implicit integer (i-m)
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
* beta(16),beta1(16),A(16),B(16),cai(16),ramda(16),
* ramda1(16),Fi0(16),Fi1(16),Fi2(16),Fi3(16),
* Fm0(16),Fm1(16),Fm2(16),Fm3(16),Fm4(16),
* fnn(16),Fn(16),rhoin(16),rhoout(16),rhol(16),
* rhoh(16),rhos(16)
if (it1 .eq. it2) then
psi1=A(it1)*exp(-alpha(it1)*(r/re(it1)-1.0_8))
psi1=psi1/(1.0_8+(r/re(it1)-cai(it1))**20)
psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0_8))
psi2=psi2/(1.0_8+(r/re(it1)-ramda(it1))**20)
psi=psi1-psi2
else
psi1=A(it1)*exp(-alpha(it1)*(r/re(it1)-1.0_8))
psi1=psi1/(1.0_8+(r/re(it1)-cai(it1))**20)
psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0_8))
psi2=psi2/(1.0_8+(r/re(it1)-ramda(it1))**20)
psia=psi1-psi2
psi1=A(it2)*exp(-alpha(it2)*(r/re(it2)-1.0_8))
psi1=psi1/(1.0_8+(r/re(it2)-cai(it2))**20)
psi2=B(it2)*exp(-beta(it2)*(r/re(it2)-1.0_8))
psi2=psi2/(1.0_8+(r/re(it2)-ramda(it2))**20)
psib=psi1-psi2
call prof(it1,r,f1)
call prof(it2,r,f2)
psi=0.5_8*(f2/f1*psia+f1/f2*psib)
endif
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c This subroutine calculates the embedding energy. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine embed(it,rho,emb)
implicit real*8 (a-h,o-z)
implicit integer (i-m)
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
* beta(16),beta1(16),A(16),B(16),cai(16),ramda(16),
* ramda1(16),Fi0(16),Fi1(16),Fi2(16),Fi3(16),
* Fm0(16),Fm1(16),Fm2(16),Fm3(16),Fm4(16),
* fnn(16),Fn(16),rhoin(16),rhoout(16),rhol(16),
* rhoh(16),rhos(16)
if (rho .lt. rhoe(it)) then
Fm33=Fm3(it)
else
Fm33=Fm4(it)
endif
if (rho .eq. 0.0_8) then
emb = 0.0_8
else if (rho .lt. rhoin(it)) then
emb=Fi0(it)+
* Fi1(it)*(rho/rhoin(it)-1.0_8)+
* Fi2(it)*(rho/rhoin(it)-1.0_8)**2+
* Fi3(it)*(rho/rhoin(it)-1.0_8)**3
else if (rho .lt. rhoout(it)) then
emb=Fm0(it)+
* Fm1(it)*(rho/rhoe(it)-1.0_8)+
* Fm2(it)*(rho/rhoe(it)-1.0_8)**2+
* Fm33*(rho/rhoe(it)-1.0_8)**3
else
emb=Fn(it)*(1.0_8-fnn(it)*log(rho/rhos(it)))*
* (rho/rhos(it))**fnn(it)
endif
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c write out set file. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine writeset
implicit real*8 (a-h,o-z)
implicit integer (i-m)
character*80 outfile,outelem
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
* beta(16),beta1(16),A(16),B(16),cai(16),ramda(16),
* ramda1(16),Fi0(16),Fi1(16),Fi2(16),Fi3(16),
* Fm0(16),Fm1(16),Fm2(16),Fm3(16),Fm4(16),
* fnn(16),Fn(16),rhoin(16),rhoout(16),rhol(16),
* rhoh(16),rhos(16)
common /pass2/ amass(16),Fr(5000,16),rhor(5000,16),
* z2r(5000,16,16),blat(16),drho,dr,rc,outfile,outelem
common /pass3/ ielement(16),ntypes,nrho,nr
character*80 struc
character(8) date
struc='fcc'
outfile = outfile(1:index(outfile,' ')-1)//'.eam.alloy'
open(unit=1,file=outfile)
call date_and_time(DATE=date)
write(1,*) 'DATE: ',date(1:4),'-',date(5:6),'-',date(7:8),
* ' UNITS: metal CONTRIBUTOR: Xiaowang Zhou xzhou@sandia.gov and ',
* 'Lucas Hale lucas.hale@nist.gov '
write(1,*) 'CITATION: X. W. Zhou, R. A. Johnson, ',
* 'H. N. G. Wadley, Phys. Rev. B, 69, 144113(2004)'
write(1,*) 'Generated by create.f'
write(1,8)ntypes,outelem
8 format(i5,' ',a24)
write(1,9)nrho,drho,nr,dr,rc
9 format(i5,e24.16,i5,2e24.16)
do 10 i=1,ntypes
write(1,11)ielement(i),amass(i),blat(i),struc
write(1,12)(Fr(j,i),j=1,nrho)
write(1,12)(rhor(j,i),j=1,nr)
10 continue
11 format(i5,2g15.5,a8)
12 format(5e24.16)
do i1=1,ntypes
do i2=1,i1
write(1,12)(z2r(i,i1,i2),i=1,nr)
end do
end do
close(1)
return
end
|
