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