From 60d03f34cc3c7ef519857a14da5d50d0698f8850 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 9 Feb 2022 10:54:15 -0500 Subject: [PATCH] update create.f with changes from NIST database also add parameters for Cr and document in README file and change the code to create output files with .eam.alloy extension --- doc/src/Tools.rst | 10 +++- tools/eam_database/EAM_code | 28 ++++++++++ tools/eam_database/README | 21 +++++-- tools/eam_database/create.f | 106 +++++++++++++++++++++--------------- 4 files changed, 114 insertions(+), 51 deletions(-) diff --git a/doc/src/Tools.rst b/doc/src/Tools.rst index 0d8ffa3f45..0788a27623 100644 --- a/doc/src/Tools.rst +++ b/doc/src/Tools.rst @@ -278,16 +278,20 @@ eam database tool ----------------------------- The tools/eam_database directory contains a Fortran program that will -generate EAM alloy setfl potential files for any combination of 16 +generate EAM alloy setfl potential files for any combination of 17 elements: Cu, Ag, Au, Ni, Pd, Pt, Al, Pb, Fe, Mo, Ta, W, Mg, Co, Ti, -Zr. The files can then be used with the :doc:`pair_style eam/alloy ` command. +Zr, Cr. The files can then be used with the :doc:`pair_style eam/alloy ` command. The tool is authored by Xiaowang Zhou (Sandia), xzhou at sandia.gov, -and is based on his paper: +with updates from Lucas Hale (NIST) lucas.hale at nist.gov and is based on his paper: X. W. Zhou, R. A. Johnson, and H. N. G. Wadley, Phys. Rev. B, 69, 144113 (2004). +The parameters for Cr were taken from: + +Lin Z B, Johnson R A and Zhigilei L V, Phys. Rev. B 77 214108 (2008). + ---------- .. _eamgn: diff --git a/tools/eam_database/EAM_code b/tools/eam_database/EAM_code index 6968a03aff..d339c87ab4 100644 --- a/tools/eam_database/EAM_code +++ b/tools/eam_database/EAM_code @@ -446,3 +446,31 @@ Zr 1.0 0.85 1.15 +Cr +2.493879 +1.793835 +17.641302 +19.60545 +8.604593 +7.170494 +1.551848 +1.827556 +0.18533 +0.277995 +-2.022754 +0.039608 +-0.183611 +-2.245972 +-2.02 +0.00 +-0.056517 +0.439144 +0.456 +-2.020038 +24 +51.9961 +0.439144 +7.170494 +0.277995 +0.85 +1.15 diff --git a/tools/eam_database/README b/tools/eam_database/README index 8d42eb7a12..fa868b583b 100644 --- a/tools/eam_database/README +++ b/tools/eam_database/README @@ -1,22 +1,35 @@ EAM database tool Xiaowang Zhou (Sandia), xzhou at sandia.gov +Revised version including fixes from https://www.ctcms.nist.gov/potentials/entry/2004--Zhou-X-W-Johnson-R-A-Wadley-H-N-G--Al/ + based on this paper: X. W. Zhou, R. A. Johnson, and H. N. G. Wadley, Phys. Rev. B, 69, 144113 (2004). +Parameters for Cr were taken from: +Lin Z B, Johnson R A and Zhigilei L V, Phys. Rev. B 77 214108 (2008) + This tool can be used to create an DYNAMO-formatted EAM -setfl file for alloy systems, using any combination +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 +Cu, Ag, Au, Ni, Pd, Pt, Al, Pb, Fe, Mo, Ta, W, Mg, Co, Ti, Zr, Cr + +WARNING: Please note that the parameter sets used here are all optimized +for the pure metals of the individual elements and that mixing rules will +be applied for creating the inter-element interactions. Those are inferior +to models where the mixed terms were specifically optimized for particular +alloys. Thus any potential files created with this tool should be used +with care and test calculations (e.g. on multiple binary mixtures) performed +to gauge the error. 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 +3) a.out < EAM.input will create an *.eam.alloy potential file +4) in DYNAMO or LAMMPS context this file is referred to as 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 index 528b4251db..ce6f03a7ac 100644 --- a/tools/eam_database/create.f +++ b/tools/eam_database/create.f @@ -1,4 +1,5 @@ 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) @@ -9,6 +10,8 @@ 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), @@ -17,9 +20,9 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc * 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 + 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' @@ -69,7 +72,7 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc read(10,*) ramda1(ntypes) read(10,*) rhol(ntypes) read(10,*) rhoh(ntypes) - blat(ntypes)=sqrt(2.0)*re(ntypes) + blat(ntypes)=sqrt(2.0_8)*re(ntypes) rhoin(ntypes)=rhol(ntypes)*rhoe(ntypes) rhoout(ntypes)=rhoh(ntypes)*rhoe(ntypes) else @@ -91,15 +94,15 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 + 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)*dr + r=(i-1.0_8)*dr if (r .lt. rst) r=rst call prof(i1,r,fvalue) if (fmax .lt. fvalue) fmax=fvalue @@ -109,7 +112,7 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc end do else do i=1,nr - r=(i-1.0)*dr + r=(i-1.0_8)*dr if (r .lt. rst) r=rst call pair(i1,i2,r,psi) z2r(i,i1,i2)=r*psi @@ -119,12 +122,13 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc end do end do 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) + 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.0)*drho + rhoF=(i-1)*drho + if (i .eq. 1) rhoF=0.0_8 call embed(it,rhoF,emb) Fr(i,it)=emb 7 continue @@ -135,20 +139,24 @@ 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)) - f=f/(1.0+(r/re(it)-ramda1(it))**20) + 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), @@ -156,25 +164,25 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc * 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) + 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)) - 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) + 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)) - 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) + 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*(f2/f1*psia+f1/f2*psib) + psi=0.5_8*(f2/f1*psia+f1/f2*psib) endif return end @@ -182,6 +190,8 @@ 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), @@ -193,18 +203,20 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc else Fm33=Fm4(it) endif - if (rho .lt. rhoin(it)) then + 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)+ - * Fi2(it)*(rho/rhoin(it)-1.0)**2+ - * Fi3(it)*(rho/rhoin(it)-1.0)**3 + * 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)+ - * Fm2(it)*(rho/rhoe(it)-1.0)**2+ - * Fm33*(rho/rhoe(it)-1.0)**3 + * 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-fnn(it)*log(rho/rhos(it)))* + emb=Fn(it)*(1.0_8-fnn(it)*log(rho/rhos(it)))* * (rho/rhos(it))**fnn(it) endif return @@ -213,6 +225,8 @@ 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), @@ -220,16 +234,20 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc * 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 + 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)//'.set' + outfile = outfile(1:index(outfile,' ')-1)//'.eam.alloy' open(unit=1,file=outfile) - write(1,*) - write(1,*) - write(1,*) + call date_and_time(DATE=date) + write(1,*) ' DATE: ',date(1:4),'-',date(5:6),'-',date(7:8),' ', + * 'CONTRIBUTOR: Xiaowang Zhou xzhou@sandia.gov and ', + * 'Lucas Hale lucas.hale@nist.gov ', + * 'CITATION: X. W. Zhou, R. A. Johnson, ', + * 'H. N. G. Wadley, Phys. Rev. B, 69, 144113(2004)' write(1,8)ntypes,outelem 8 format(i5,' ',a24) write(1,9)nrho,drho,nr,dr,rc