! Program to generate LAMMPS format potentials input for ! Mendelev-Underwood-Ackland Temperature-dependent Ti potential, ! Compile the fortran code using, e.g. gfortran LAMMPS_MUA_Tdep.f ! Run the code, e.g. type ./a.out ! it will echo "Input Temperature between 100,999K" ! type in the required temperature. ! Code generates three LAMPS-compatible input ! files Ti1.eam.fs, Ti2.eam.fs, TiT=XXX.eam.fs ! where XXX is the temperature you plan to simulate at (use NVT of NPT ensemble) ! Below 100K use the Ti2 potential, above 1000K, use Ti1 ! Bugs to: gjackland@ed.ac.uk program LAMMPS_MUA_Tdep IMPLICIT REAL*8(A-H,O-Z) character*13 filename real*8:: bier1=0.1818, bier2=0.5099, bier3=0.2802, bier4=0.02817 ! Potential parameters hardcoded later type potparameters real*8::phicor4 real*8:: vcut2,vcut3,vcut4,pcut real*8:: v14,v15,v16,v17,v18 real*8:: v24,v25,v26,v27,v28 real*8:: v34,v35,v36,v37,v38 real*8:: f1,f2,f3,f4,f5 real*8::emcut1,emcut2,emcut3,emcut4,emcut5,emcut6 real*8::splice1,splice2,splice3,splice4 real*8:: ex0,ex1,ex2,ex3,ex4 end type potparameters type(potparameters) :: Ti1, Ti2 ! Set parameters for Ti1 and Ti2 call set_Ti1() call set_Ti2() ! Titanium specific things na1 = 22 amass=47.867 clatt=2.94653381137745 cutoff=Ti1%vcut4 ! Temperature dependent stuff write(*,*)" Input Temperature between 100,999K" read(*,*) ITEMP T0=600 Tw=100 ! Ti1 is the high temperature version, Ti2 the Low T one GT = tanh((ITEMP-T0)/Tw) GT1 = (1d0+GT)/2d0 GT2 = (1d0-GT)/2d0 Nrho=10000 Drho=0.03d0 Nr=10000 Dr=cutoff/Nr ! LAMMPS converterlines 1,2,3 = comments (ignored) !line 4: Nelements Element1 Element2 ... ElementN !line 5: Nrho, drho, Nr, dr, cutoff ! line 1 = atomic number, mass, lattice constant, lattice type ! embedding function F(rho) (Nrho values) ! density function rho(r) for element 1 at element 1 (Nr values) ! Nr values of pair potential phi(r)Â write(filename,27)"TiT",ITEMP,".eam.fs" 27 format(A3,I3,A7) write(*,*)"LAMMPS potential file:",filename open(unit=25, file='Ti1.eam.fs') open(unit=26, file='Ti2.eam.fs') open(unit=27, file=filename) ! write(25,*) "Potentials by Mendelev, Underwood and Ackland", & " Journal of Chemical Physics 2016" write(25,*) "Contact information: mikhail.mendelev@gmail.com" write(25,*) "Titanium potential Ti1, Temperature = High" write(25,*) "1 Ti" write(25,*) Nrho,Drho,Nr,Dr,cutoff write(25,*) na1, amass, clatt, "hcp" write(26,*) "Potentials by Mendelev, Underwood and Ackland", & " Journal of Chemical Physics 2016" write(26,*) "Contact information: mikhail.mendelev@gmail.com" write(26,*) "Titanium potential Ti2, Temperature = Low" write(26,*) "1 Ti" write(26,*) Nrho,Drho,Nr,Dr,cutoff write(26,*) na1, amass, clatt, "hcp" write(27,*) "Ti Sommerfeld potential", & " by Mendelev, Underwood and Ackland,", & " Journal of Chemical Physics 2016" write(27,*) "Contact information: gjackland@ed.ac.uk" write(27,*) "Use in NVT/NPT with Temperature=",ITEMP,"K" write(27,*) "1 Ti" write(27,*) Nrho,Drho,Nr,Dr,cutoff write(27,*) na1, amass, clatt, "hcp" rad=0.0d0 do I = 1,Nrho,5 v1=embed(rad,Ti1) v1x=embed(rad,Ti2) rad=rad+Drho v2=embed(rad,Ti1) v2x=embed(rad,Ti2) rad=rad+Drho v3=embed(rad,Ti1) v3x=embed(rad,Ti2) rad=rad+Drho v4=embed(rad,Ti1) v4x=embed(rad,Ti2) rad=rad+Drho v5=embed(rad,Ti1) v5x=embed(rad,Ti2) rad=rad+Drho write(25,25)v1,v2,v3,v4,v5 write(26,25)v1x,v2x,v3x,v4x,v5x write(27,25)(v1*gt1+v1x*gt2),(v2*gt1+v2x*gt2),(v3*gt1+v3x*gt2), & (v4*gt1+v4x*gt2),(v5*gt1+v5x*gt2) 25 format(5e16.8) enddo rad=0.0d0 do I = 1,Nr,5 v1=phi_src(rad,Ti1) rad=rad+Dr v2=phi_src(rad,Ti1) rad=rad+Dr v3=phi_src(rad,Ti1) rad=rad+Dr v4=phi_src(rad,Ti1) rad=rad+Dr v5=phi_src(rad,Ti1) rad=rad+Dr write(25,25)v1,v2,v3,v4,v5 write(26,25)v1,v2,v3,v4,v5 write(27,25)v1,v2,v3,v4,v5 enddo rad=0.0d0 do I = 1,Nr,5 v1=vee_src(rad,Ti1)*rad v1x=vee_src(rad,Ti2)*rad rad=rad+Dr v2=vee_src(rad,Ti1)*rad v2x=vee_src(rad,Ti2)*rad rad=rad+Dr v3=vee_src(rad,Ti1)*rad v3x=vee_src(rad,Ti2)*rad rad=rad+Dr v4=vee_src(rad,Ti1)*rad v4x=vee_src(rad,Ti2)*rad rad=rad+Dr v5=vee_src(rad,Ti1)*rad v5x=vee_src(rad,Ti2)*rad rad=rad+Dr write(25,25)v1,v2,v3,v4,v5 write(26,25)v1x,v2x,v3x,v4x,v5x write(27,25)(v1*gt1+v1x*gt2),(v2*gt1+v2x*gt2), & (v3*gt1+v3x*gt2),(v4*gt1+v4x*gt2),(v5*gt1+v5x*gt2) enddo rad=0.0 stop contains function phi_src(r,pot) real *8 :: phi_src ! result (eV) real *8, intent(in) :: r ! separation (angstrom) type(potparameters), intent(in) :: pot phi_src=0.0 if(r.lt.pot%pcut) & phi_src=pot%phicor4*exp(0.1d0*r)*(r-pot%pcut)**4 return end function phi_src !---------------------------------------------------------------------------- ! ! emb_src ! !---------------------------------------------------------------------------- function embed(rhoz,pot) real *8 :: embed ! result (eV) real *8, intent(in) :: rhoz type(potparameters), intent(in) :: pot embed = -sqrt(rhoz) IF(rhoz.lt.pot%emcut1)return embed = embed+pot%f1*(rhoz-pot%emcut1)**4 IF(rhoz.lt.pot%emcut2)return embed = embed+pot%f2*(rhoz-pot%emcut2)**4 IF(rhoz.lt.pot%emcut3)return embed = embed+pot%f3*(rhoz-pot%emcut3)**4 IF(rhoz.lt.pot%emcut4)return embed = embed+pot%f4*(rhoz-pot%emcut4)**4 IF(rhoz.lt.pot%emcut5)return embed = embed+pot%f5*(rhoz-pot%emcut5)**4 c Correction to Ti1 at extreme densities IF(RHOZ.gt.pot%emcut6) then embed=pot%ex0+pot%ex1*(rhoz-pot%emcut6) & +pot%ex2*(rhoz-pot%emcut6)**2 & +pot%ex3*(rhoz-pot%emcut6)**3 & +pot%ex4*(rhoz-pot%emcut6)**4 c embed=1.3533741622097+5.9660537422044E-01*(rhoz-pot%emcut6) c & -1.4665711508917E-02*(rhoz-pot%emcut6)**2 c & +6.9056785356062E-04*(rhoz-pot%emcut6)**3 c & -3.5388769339836E-05*(rhoz-pot%emcut6)**4 endif return end function embed !---------------------------------------------------------------------------- ! ! demb_src ! !---------------------------------------------------------------------------- pure function demb_src(rhoz) real *8 :: demb_src ! result (eV) real *8, intent(in) :: rhoz ! Compute the embedding function derivative if(rhoz.gt.1.d-10) then demb_src = -0.5d0/sqrt(rhoz) & -f1*(rhoz-emcut1)**3*4*xh0(rhoz-emcut1) & -f2*(rhoz-emcut2)**3*4*xh0(rhoz-emcut2) & -f3*(rhoz-emcut3)**3*4*xh0(rhoz-emcut3) & -f4*(rhoz-emcut4)**3*4*xh0(rhoz-emcut4) & -f5*(rhoz-emcut5)**3*4*xh0(rhoz-emcut5) else demb_src=0.0d0 endif return end function demb_src function vee_src(r,pot) real *8 :: vee_src, x ! result (eV) real *8, intent(in) :: r ! separation (angstrom) type(potparameters), intent(in) :: pot VEE_SRC = 0.0d0 IF(R.le.0.001d0) RETURN ! Biersack-Zeigler functions for titanium IF(R.LT.1d0) THEN VEE_SRC = (6.9694830017999E+03/r)*( & bier1*exp(-2.7066765729847E+01*r) & +bier2*exp(-7.9703166710110E+00*r) & +bier3*exp(-3.4078749726736E+00*r) & +bier4*exp(-1.7052062409804E+00*r)) ELSEIF(R.LT.2.3d0) THEN VEE_SRC = +exp(pot%splice1+pot%splice2*r & +pot%splice3*r**2+pot%splice4*r**3) ELSEIF(r.lt.pot%vcut4)THEN VEE_SRC = (pot%v14*(pot%vcut2-r)**4 & +pot%v15*(pot%vcut2-r)**5 & +pot%v16*(pot%vcut2-r)**6 & +pot%v17*(pot%vcut2-r)**7 & +pot%v18*(pot%vcut2-r)**8) & *xH0(pot%vcut2-r) + & (pot%v24*(pot%vcut3-r)**4 & +pot%v25*(pot%vcut3-r)**5 & +pot%v26*(pot%vcut3-r)**6 & +pot%v27*(pot%vcut3-r)**7 & +pot%v28*(pot%vcut3-r)**8) & *xH0(pot%vcut3-r) + & (pot%v34*(pot%vcut4-r)**4 & +pot%v35*(pot%vcut4-r)**5 & +pot%v36*(pot%vcut4-r)**6 & +pot%v37*(pot%vcut4-r)**7 & +pot%v38*(pot%vcut4-r)**8) & *xH0(pot%vcut4-r) ENDIF return end function vee_src !---------------------------------------------------------------------------- ! ! Utility Function - Heaviside Step. ! !---------------------------------------------------------------------------- pure function xH0(xin) real *8, intent(in) :: xin ! separation (angstrom) real*8 :: xh0 if(xin.lt.0.0d0)then xH0=0.d0 else xH0=1.d0 endif return end function xh0 subroutine set_Ti1() Ti1%phicor4 =7.9483D-02 Ti1%vcut2=3.2d0 Ti1%vcut3=5.4d0 Ti1%vcut4=6.9d0 Ti1%pcut=5.1d0 Ti1%v14= -19.000009669827 Ti1%v15 = +6.9503328169194 Ti1%v16 = +1.0141463034640E+01 Ti1% v17= -2.3622475018118E+01 Ti1%v18= +8.4715024861295E+00 Ti1%v24= +4.2964815851604E-02 Ti1%v25= -3.1925338689760E+00 Ti1%v26= -8.6997411310845E-01 Ti1%v27= -1.3418235164880E+00 Ti1%v28= -9.3059971805003E-02 Ti1%v34= +7.8873370228139E-01 Ti1%v35= -2.0873451358068 Ti1%v36= +2.0860619467995 Ti1%v37= -9.2402710034568E-01 Ti1%v38= +1.5279931337437E-01 Ti1%f1= 6.4458935128222E-06 Ti1%f2= -4.6861183499214E-05 Ti1%f3= 1.5881962673345E-04 Ti1%f4= -3.3153247278010E-05 Ti1%f5= -1.0727500853945E-04 Ti1%emcut1=18.d0 Ti1%emcut2=22.d0 Ti1%emcut3=26.d0 Ti1%emcut4=28.d0 Ti1%emcut5=30.d0 Ti1%emcut6=50.d0 Ti1%splice1=10.378404275169 Ti1%splice2=-8.2358378101067 Ti1%splice3= +3.0636189315873 Ti1%splice4=-0.58543836194351 Ti1%ex0=-1.3533741622097 Ti1%ex1=+5.9660537422044E-01 Ti1%ex2= 1.4665711508917E-02 Ti1%ex3= -6.9056785356062E-04 Ti1%ex4= +3.5388769339836E-05 end subroutine set_Ti2() Ti2%phicor4 =7.9483D-02 Ti2%vcut2=3.2d0 Ti2%vcut3=5.4d0 Ti2%vcut4=6.9d0 Ti2%pcut=5.1d0 Ti2%v14= -6.4522391868171 Ti2%v15 =-7.4406432070254 Ti2%v16 = -1.1338596723822 Ti2%v17= +11.331538165644 Ti2%v18= -8.4070436781176 Ti2%v24= +.63252414003875 Ti2%v25= -1.3807804217433 Ti2%v26= +.65134886927909 Ti2%v27= -.50882595849649 Ti2%v28= 2.2019845067084E-02 Ti2%v34= +5.5203429076280E-02 Ti2%v35= -2.8338491037854E-01 Ti2%v36= +3.3760474001764E-01 Ti2%v37= -1.5713566606326E-01 Ti2%v38= +2.5880555267166E-02 Ti2%f1= 3.7224420797114E-05 Ti2%f2= -1.4581966752354E-04 Ti2%f3= 2.5307706881546E-04 Ti2%f4= +1.3022021793393E-04 Ti2%f5= -3.4570352330604E-04 Ti2%emcut1=18.d0 Ti2%emcut2=22.d0 Ti2%emcut3=26.d0 Ti2%emcut4=28.d0 Ti2%emcut5=30.d0 Ti2%emcut6=45.d0 Ti2%splice1=9.2294088533131E+00 Ti2%splice2=-5.4768063605710E+00 Ti2%splice3= 9.9254229808463E-01 Ti2%splice4=-1.2439775612028E-01 Ti2%ex0=-1.3758034593242E+00 Ti2%ex1=5.9497853813975E-01 Ti2%ex2=7.6687260515563E-03 Ti2%ex3=-2.0512035653105E-03 Ti2%ex4=3.4861183325532E-04 end end program LAMMPS_MUA_Tdep