博文
Fortran program writing *.eam.alloy potential for lammps
||||
This blog is about a program written by my self for generating the *.eam.alloy file for lammps, which has been finished 2012.04.18. Potential bugs must exist, therefore, advice and suggestion about this program is highly appreciated.
https://m.sciencenet.cn/blog-646691-560818.html
上一篇：Interesting potential strategy from artificial intelligence
!#############################################################################################################################
!EAM_Table Program Version 1.0
!Generating the setfl file for eam/alloy potential in lammps
!Programmer: Bin Ouyang, McGill University
!Supervisor: Jun Song, McGill uUniversity
!Department of Mining and Materials Engineering, McGill University
!Date: 2012.04.17
!#############################################################################################################################
program eam_table
implicit none
character*50 :: para_file, out_file
real :: re(2), fe(2), rou_e(2), arfa(2), beta(2), A(2), B(2), kai(2), lambda(2), Fn0(2), Fn1(2)
real :: Fn2(2), Fn3(2), F0(2), F1(2), F2(2), F3(2), ita(2), F_e(2)
real :: drho, dr, cutoff, latt_c(2), mass(2)
integer :: Nrho, Nr, atom_num, ind(2)
character* 50 :: comment(3)
character* 7 :: elem(2), latt(2)
real,allocatable :: rho(:, :), F(:, :), phi(:, :), radius(:, :)
integer :: i, j
real :: rou_n(2), rou_0(2)
namelist /eam_para/ re, fe, rou_e, arfa, beta, A, B, kai, lambda, Fn0, Fn1, Fn2&
&, Fn3, F0, F1, F2, F3, ita, F_e
namelist /eam_format/ drho, dr, cutoff, latt_c, Nrho, Nr, atom_num, ind, elem, comment, latt, mass
do while(1)
!display the welcome screen
write(*,*)'EAM_Table Program Version 1.0'
write(*,*)'====================================================================='
write(*,*)'Generating the setfl format file for eam/alloy potential in lammps'
write(*,*)'Programmer: Bin Ouyang, McGill university'
write(*,*)'Supervisor: Jun Song, McGill university'
write(*,*)'Department of Mining and Materials Engineering, McGill University'
write(*,*)'Date: 2012.04.17'
write(*,*)'======================================================================'
write(*,*)' '
write(*,*)'Please enter the name of eam parameter file'
read(*,*) para_file
!reading the parameter file of eam potential
open(1,file = para_file,status = 'old')
read(1,eam_para)
read(1,eam_format)
close(1)
allocate(rho(Nrho, atom_num), F(Nrho, atom_num), radius(Nrho, atom_num))
allocate(phi(Nrho, atom_num + 1))
write(*,*)' '
write(*,*)'Reading parameter file successfully!'
write(*,*)'Please enter the name of output file'
write(*,*)'Format: ***.eam.alloy'
read(*,*)out_file
write(*,*)' '
!Calculating the required table of parameters for EAM
do i = 1, Nrho
radius(i, 1) = (i-1) * dr
radius(i, 2) = radius(i, 1)
enddo !discrete the radius
do i = 1, 2
rho(:, i) = fe(i) * exp(-1 * beta(i) * (radius(:, i) / re(i) - 1)) / (1&
&+ (radius(:, i) / re(i) - lambda(i)) ** 20)
enddo !calculating the density of electonic
rou_n = rou_e * 0.85
rou_0 = rou_e * 1.15
do i = 1, 2
do j = 1, Nr
if(rho(j, i) < rou_n(i))then
F(j, i) = Fn0(i) + Fn1(i) * (rho(j, i) / 0.85 / rou_e(i)) + Fn2(i)&
& * (rho(j, i) / 0.85 / rou_e(i)) ** 2 + Fn3(i) * (rho(j, i) / 0.85 / rou_e(i)) ** 3
elseif(rho(j, i) < rou_0(i))then
F(j, i) = F0(i) + F1(i) * (rho(j, i) / rou_e(i)) + F2(i) * (rho(j, i)&
& / rou_e(i)) ** 2 + F3(i) * (rho(j, i) / rou_e(i)) ** 3
else
F(j, i) = F_e(i) * (1 - ita(i) * log(rho(j, i) / rou_e(i))) * (rho(j, i) / rou_e(i)) ** ita(i)
endif
enddo
enddo
do i = 1, Nr
phi(i, 1) = A(1) * exp(arfa(1) * (1 - radius(i,1) / re(1))) / (1 + (radius(i,1) / re(1) - kai(1)) ** 20)
enddo
do i = 1, Nr
phi(i, 3) = A(2) * exp(arfa(2) * (1 - radius(2,2) / re(2))) / (1 + (radius(i,2) / re(2) - kai(2)) ** 20)
enddo
do i = 1, Nr
phi(i, 2) = 0.5 * (rho(i, 1) / rho(i, 2) * phi(i, 3) + rho(i, 2) / rho(i, 1) * phi(i, 1))
enddo
!End of calculation
!Wrting the eam/alloy file
open(2, file = out_file)
write(2,*) comment(1)
write(2,*) comment(2)
write(2,*) comment(3)
write(2,*) atom_num, ' ', elem(1), elem(2)
write(2,*) Nrho, drho, Nr, dr, cutoff
write(2,*) ind(1), mass(1), latt_c(1), latt(1)
do i = 1, Nrho / 5
write(2,*) F(5 * i - 4, 1), F(5 * i - 3, 1), F(5 * i - 2, 1), F(5 * i - 1, 1), F(5 * i, 1)
enddo
do i = 1, Nr / 5
write(2,*) rho(5 * i - 4, 1), rho(5 * i - 3, 1), rho(5 * i - 2, 1), rho(5 * i - 1, 1), rho(5 * i, 1)
enddo
write(2,*) ind(2), mass(2), latt_c(2), latt(2)
do i = 1, Nrho / 5
write(2,*) F(5 * i - 4, 2), F(5 * i - 3, 2), F(5 * i - 2, 2), F(5 * i - 1, 2), F(5 * i, 2)
enddo
do i = 1, Nr / 5
write(2,*) rho(5 * i - 4, 2), rho(5 * i - 3, 2), rho(5 * i - 2, 2), rho(5 * i - 1, 2), rho(5 * i, 2)
enddo
do i = 1, Nr / 5
write(2,*) phi(5 * i - 4, 1), phi(5 * i - 3, 1), phi(5 * i - 2, 1), phi(5 * i - 1, 1), phi(5 * i, 1)
enddo
do i = 1, Nr / 5
write(2,*) phi(5 * i - 4, 2), phi(5 * i - 3, 2), phi(5 * i - 2, 2), phi(5 * i - 1, 2), phi(5 * i, 2)
enddo
do i = 1, Nr / 5
write(2,*) phi(5 * i - 4, 3), phi(5 * i - 3, 3), phi(5 * i - 2, 3), phi(5 * i - 1, 3), phi(5 * i, 3)
enddo
close(2)
write(*,*) 'The eam potential file has been output succesfully!'
deallocate(rho, F, radius, phi)
enddo
end program
https://m.sciencenet.cn/blog-646691-560818.html
上一篇：Interesting potential strategy from artificial intelligence
