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

@@ -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 <pair_eam>` command.
+Zr, Cr.  The files can then be used with the :doc:`pair_style eam/alloy <pair_eam>` 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:

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

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
+3) a.out < EAM.input will create an *.eam.alloy potential file
-4) in DYNAMO or LAMMPS lingo, this is a setfl 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

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),
+      common /pass2/ amass(16),Fr(5000,16),rhor(5000,16),
-     *   rhor(5000,16),z2r(5000,16,16),ntypes,blat(16),
+     *   z2r(5000,16,16),blat(16),drho,dr,rc,outfile,outelem
-     *   nrho,drho,nr,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
+      rc=sqrt(10.0_8)/2.0_8*alatmax
-      rst=0.5
+      rst=0.5_8
-      dr=rc/(nr-1.0)
+      dr=rc/(nr-1.0_8)
-      fmax=-1.0
+      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. 2.0_8*rhoemax) rhom=2.0_8*rhoemax
-      if (rhom .lt. 100.0) rhom=100.0
+      if (rhom .lt. 100.0_8) rhom=100.0_8
-      drho=rhom/(nrho-1.0)
+      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=fe(it)*exp(-beta1(it)*(r/re(it)-1.0_8))
-      f=f/(1.0+(r/re(it)-ramda1(it))**20)
+      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=A(it1)*exp(-alpha(it1)*(r/re(it1)-1.0_8))
-         psi1=psi1/(1.0+(r/re(it1)-cai(it1))**20)
+         psi1=psi1/(1.0_8+(r/re(it1)-cai(it1))**20)
-         psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0))
+         psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0_8))
-         psi2=psi2/(1.0+(r/re(it1)-ramda(it1))**20)
+         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=A(it1)*exp(-alpha(it1)*(r/re(it1)-1.0_8))
-         psi1=psi1/(1.0+(r/re(it1)-cai(it1))**20)
+         psi1=psi1/(1.0_8+(r/re(it1)-cai(it1))**20)
-         psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0))
+         psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0_8))
-         psi2=psi2/(1.0+(r/re(it1)-ramda(it1))**20)
+         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=A(it2)*exp(-alpha(it2)*(r/re(it2)-1.0_8))
-         psi1=psi1/(1.0+(r/re(it2)-cai(it2))**20)
+         psi1=psi1/(1.0_8+(r/re(it2)-cai(it2))**20)
-         psi2=B(it2)*exp(-beta(it2)*(r/re(it2)-1.0))
+         psi2=B(it2)*exp(-beta(it2)*(r/re(it2)-1.0_8))
-         psi2=psi2/(1.0+(r/re(it2)-ramda(it2))**20)
+         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)+
+     *       Fi1(it)*(rho/rhoin(it)-1.0_8)+
-     *       Fi2(it)*(rho/rhoin(it)-1.0)**2+
+     *       Fi2(it)*(rho/rhoin(it)-1.0_8)**2+
-     *       Fi3(it)*(rho/rhoin(it)-1.0)**3
+     *       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)+
+     *       Fm1(it)*(rho/rhoe(it)-1.0_8)+
-     *       Fm2(it)*(rho/rhoe(it)-1.0)**2+
+     *       Fm2(it)*(rho/rhoe(it)-1.0_8)**2+
-     *       Fm33*(rho/rhoe(it)-1.0)**3
+     *       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),
+      common /pass2/ amass(16),Fr(5000,16),rhor(5000,16),
-     *   rhor(5000,16),z2r(5000,16,16),ntypes,blat(16),
+     *   z2r(5000,16,16),blat(16),drho,dr,rc,outfile,outelem
-     *   nrho,drho,nr,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,*)
+      call date_and_time(DATE=date)
-      write(1,*)
+      write(1,*) ' DATE: ',date(1:4),'-',date(5:6),'-',date(7:8),'  ',
-      write(1,*)
+     *   '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