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