diff --git a/tools/eam_database/EAM.input b/tools/eam_database/EAM.input new file mode 100644 index 0000000000..a6529616ba --- /dev/null +++ b/tools/eam_database/EAM.input @@ -0,0 +1,9 @@ + &funccard + atomtype='Ta' + &end + &funccard + atomtype='Cu' + &end + &funccard + &end + diff --git a/tools/eam_database/EAM_code b/tools/eam_database/EAM_code new file mode 100644 index 0000000000..6968a03aff --- /dev/null +++ b/tools/eam_database/EAM_code @@ -0,0 +1,448 @@ +Cu +2.556162 +1.554485 +21.175871 +21.175395 +8.127620 +4.334731 +0.396620 +0.548085 +0.308782 +0.756515 +-2.170269 +-0.263788 +1.088878 +-0.817603 +-2.19 +0.00 +0.561830 +-2.100595 +0.310490 +-2.186568 +29 +63.546 +-2.100595 +4.334731 +0.756515 +0.85 +1.15 +Ag +2.891814 +1.106232 +14.604100 +14.604144 +9.132010 +4.870405 +0.277758 +0.419611 +0.339710 +0.750758 +-1.729364 +-0.255882 +0.912050 +-0.561432 +-1.75 +0.00 +0.744561 +-1.150650 +0.783924 +-1.748423 +47 +107.8682 +-1.150650 +4.870405 +0.750758 +0.85 +1.15 +Au +2.885034 +1.529021 +19.991632 +19.991509 +9.516052 +5.075228 +0.229762 +0.356666 +0.356570 +0.748798 +-2.937772 +-0.500288 +1.601954 +-0.835530 +-2.98 +0.00 +1.706587 +-1.134778 +1.021095 +-2.978815 +79 +196.96654 +-1.134778 +5.075228 +0.748798 +0.85 +1.15 +Ni +2.488746 +2.007018 +27.562015 +27.562031 +8.383453 +4.471175 +0.429046 +0.633531 +0.443599 +0.820658 +-2.693513 +-0.076445 +0.241442 +-2.375626 +-2.70 +0.00 +0.265390 +-0.152856 +0.445470 +-2.7 +28 +58.6934 +-0.152856 +4.471175 +0.820658 +0.85 +1.15 +Pd +2.750897 +1.595417 +21.335246 +21.940073 +8.697397 +4.638612 +0.406763 +0.598880 +0.397263 +0.754799 +-2.321006 +-0.473983 +1.615343 +-0.231681 +-2.36 +0.00 +1.481742 +-1.675615 +1.130000 +-2.352753 +46 +106.42 +-1.675615 +4.638612 +0.754799 +0.85 +1.15 +Pt +2.771916 +2.336509 +33.367564 +35.205357 +7.105782 +3.789750 +0.556398 +0.696037 +0.385255 +0.770510 +-1.455568 +-2.149952 +0.528491 +1.222875 +-4.17 +0.00 +3.010561 +-2.420128 +1.450000 +-4.145597 +78 +195.08 +-2.420128 +3.789750 +0.770510 +0.25 +1.15 +Al +2.863924 +1.403115 +20.418205 +23.195740 +6.613165 +3.527021 +0.314873 +0.365551 +0.379846 +0.759692 +-2.807602 +-0.301435 +1.258562 +-1.247604 +-2.83 +0.00 +0.622245 +-2.488244 +0.785902 +-2.824528 +13 +26.981539 +-2.488244 +3.527021 +0.759692 +0.85 +1.15 +Pb +3.499723 +0.647872 +8.450154 +8.450063 +9.121799 +5.212457 +0.161219 +0.236884 +0.250805 +0.764955 +-1.422370 +-0.210107 +0.682886 +-0.529378 +-1.44 +0.00 +0.702726 +-0.538766 +0.935380 +-1.439436 +82 +207.2 +-0.538766 +5.212457 +0.764955 +0.85 +1.15 +Fe +2.481987 +1.885957 +20.041463 +20.041463 +9.818270 +5.236411 +0.392811 +0.646243 +0.170306 +0.340613 +-2.534992 +-0.059605 +0.193065 +-2.282322 +-2.54 +0.00 +0.200269 +-0.148770 +0.391750 +-2.539945 +26 +55.847 +-0.148770 +5.236411 +0.340613 +0.85 +1.15 +Mo +2.728100 +2.723710 +29.354065 +29.354065 +8.393531 +4.476550 +0.708787 +1.120373 +0.137640 +0.275280 +-3.692913 +-0.178812 +0.380450 +-3.133650 +-3.71 +0.00 +0.875874 +0.776222 +0.790879 +-3.712093 +42 +95.94 +0.776222 +4.476550 +0.275280 +0.85 +1.15 +Ta +2.860082 +3.086341 +33.787168 +33.787168 +8.489528 +4.527748 +0.611679 +1.032101 +0.176977 +0.353954 +-5.103845 +-0.405524 +1.112997 +-3.585325 +-5.14 +0.00 +1.640098 +0.221375 +0.848843 +-5.141526 +73 +180.9479 +0.221375 +4.527748 +0.353954 +0.85 +1.15 +W +2.740840 +3.487340 +37.234847 +37.234847 +8.900114 +4.746728 +0.882435 +1.394592 +0.139209 +0.278417 +-4.946281 +-0.148818 +0.365057 +-4.432406 +-4.96 +0.00 +0.661935 +0.348147 +0.582714 +-4.961306 +74 +183.84 +0.348147 +4.746728 +0.278417 +0.85 +1.15 +Mg +3.196291 +0.544323 +7.132600 +7.132600 +10.228708 +5.455311 +0.137518 +0.225930 +0.5 +1.0 +-0.896473 +-0.044291 +0.162232 +-0.689950 +-0.90 +0.00 +0.122838 +-0.226010 +0.431425 +-0.899702 +12 +24.305 +-0.226010 +5.455311 +1.0 +0.85 +1.15 +Co +2.505979 +1.975299 +27.206789 +27.206789 +8.679625 +4.629134 +0.421378 +0.640107 +0.5 +1.0 +-2.541799 +-0.219415 +0.733381 +-1.589003 +-2.56 +0.00 +0.705845 +-0.687140 +0.694608 +-2.559307 +27 +58.9332 +-0.687140 +4.629134 +1.0 +0.85 +1.15 +Ti +2.933872 +1.863200 +25.565138 +25.565138 +8.775431 +4.680230 +0.373601 +0.570968 +0.5 +1.0 +-3.203773 +-0.198262 +0.683779 +-2.321732 +-3.22 +0.00 +0.608587 +-0.750710 +0.558572 +-3.219176 +22 +47.88 +-0.750710 +4.680230 +1.0 +0.85 +1.15 +Zr +3.199978 +2.230909 +30.879991 +30.879991 +8.559190 +4.564902 +0.424667 +0.640054 +0.5 +1.0 +-4.485793 +-0.293129 +0.990148 +-3.202516 +-4.51 +0.00 +0.928602 +-0.981870 +0.597133 +-4.509025 +40 +91.224 +-0.981870 +4.564902 +1.0 +0.85 +1.15 diff --git a/tools/eam_database/README b/tools/eam_database/README new file mode 100644 index 0000000000..8d42eb7a12 --- /dev/null +++ b/tools/eam_database/README @@ -0,0 +1,22 @@ +EAM database tool +Xiaowang Zhou (Sandia), xzhou at sandia.gov + +based on this paper: + +X. W. Zhou, R. A. Johnson, and H. N. G. Wadley, Phys. Rev. B, 69, +144113 (2004). + +This tool can be used to create an DYNAMO-formatted EAM +setfl file for alloy systems, using any combination +of the elements discussed in the paper and listed in +the EAM_code file, namely: + +Cu, Ag, Au, Ni, Pd, Pt, Al, Pb, Fe, Mo, Ta, W, Mg, Co, Ti, Zr + +Steps: + +1) compile create.f -> a.out (e.g. gfortran create.f) +2) edit the input file EAM.input to list 2 or more desired elements to include +3) a.out < EAM.input will create an EAM *.set file +4) in DYNAMO or LAMMPS lingo, this is a setfl file + that can be used with the LAMMPS pair_style eam/alloy command diff --git a/tools/eam_database/create.f b/tools/eam_database/create.f new file mode 100644 index 0000000000..2bdf49dd0d --- /dev/null +++ b/tools/eam_database/create.f @@ -0,0 +1,248 @@ +C author: X. W. Zhou, xzhou@sandia.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 + 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/ ielement(16),amass(16),Fr(5000,16), + * rhor(5000,16),z2r(5000,16,16),ntypes,blat(16), + * nrho,drho,nr,dr,rc,outfile,outelem + 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)*re(ntypes) + rhoin(ntypes)=rhol(ntypes)*rhoe(ntypes) + rhoout(ntypes)=rhoh(ntypes)*rhoe(ntypes) + else + do 1 i=1,27 +1 read(10,*)vtmp + 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)/2.0*alatmax + rst=0.5 + dr=rc/(nr-1.0) + fmax=-1.0 + do 3 i1=1,ntypes + do 3 i2=1,i1 + if ( i1 .eq. i2) then + do 4 i=1,nr + r=(i-1.0)*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 +4 continue + else + do 5 i=1,nr + r=(i-1.0)*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) +5 continue + endif +3 continue + rhom=fmax + if (rhom .lt. 2.0*rhoemax) rhom=2.0*rhoemax + if (rhom .lt. 100.0) rhom=100.0 + drho=rhom/(nrho-1.0) + do 6 it=1,ntypes + do 7 i=1,nrho + rhoF=(i-1.0)*drho + 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) + 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)) + f=f/(1.0+(r/re(it)-ramda1(it))**20) + return + end +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c This subroutine calculates the pair potential. c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + subroutine pair(it1,it2,r,psi) + 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)) + psi1=psi1/(1.0+(r/re(it1)-cai(it1))**20) + psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0)) + psi2=psi2/(1.0+(r/re(it1)-ramda(it1))**20) + psi=psi1-psi2 + else + psi1=A(it1)*exp(-alpha(it1)*(r/re(it1)-1.0)) + psi1=psi1/(1.0+(r/re(it1)-cai(it1))**20) + psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0)) + psi2=psi2/(1.0+(r/re(it1)-ramda(it1))**20) + psia=psi1-psi2 + psi1=A(it2)*exp(-alpha(it2)*(r/re(it2)-1.0)) + psi1=psi1/(1.0+(r/re(it2)-cai(it2))**20) + psi2=B(it2)*exp(-beta(it2)*(r/re(it2)-1.0)) + psi2=psi2/(1.0+(r/re(it2)-ramda(it2))**20) + psib=psi1-psi2 + call prof(it1,r,f1) + call prof(it2,r,f2) + psi=0.5*(f2/f1*psia+f1/f2*psib) + endif + return + end +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c This subroutine calculates the embedding energy. c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + subroutine embed(it,rho,emb) + 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 .lt. rhoin(it)) then + emb=Fi0(it)+ + * Fi1(it)*(rho/rhoin(it)-1.0)+ + * Fi2(it)*(rho/rhoin(it)-1.0)**2+ + * Fi3(it)*(rho/rhoin(it)-1.0)**3 + else if (rho .lt. rhoout(it)) then + emb=Fm0(it)+ + * Fm1(it)*(rho/rhoe(it)-1.0)+ + * Fm2(it)*(rho/rhoe(it)-1.0)**2+ + * Fm33*(rho/rhoe(it)-1.0)**3 + else + emb=Fn(it)*(1.0-fnn(it)*log(rho/rhos(it)))* + * (rho/rhos(it))**fnn(it) + endif + return + end +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c write out set file. c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + subroutine writeset + 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/ ielement(16),amass(16),Fr(5000,16), + * rhor(5000,16),z2r(5000,16,16),ntypes,blat(16), + * nrho,drho,nr,dr,rc,outfile,outelem + character*80 struc + struc='fcc' + outfile = outfile(1:index(outfile,' ')-1)//'.set' + open(unit=1,file=outfile) + write(1,*) + write(1,*) + write(1,*) + 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 13 i1=1,ntypes + do 13 i2=1,i1 + write(1,12)(z2r(i,i1,i2),i=1,nr) +13 continue + close(1) + return + end