104 lines
2.5 KiB
Fortran
104 lines
2.5 KiB
Fortran
program madelung
|
|
|
|
implicit none
|
|
|
|
external :: dist
|
|
integer :: nx, ny, nz, ix, iy, iz
|
|
integer :: nat, iat, jat
|
|
real*8 :: r0, r1, r2, rexp, fact
|
|
real*8 :: dx, dy, dz
|
|
real*8 :: qtot, etot, e0, e1, e2, e3
|
|
real*8, allocatable :: q(:), x(:), y(:), z(:)
|
|
real*8 :: xx, yy, zz
|
|
real*8 :: dist, dd
|
|
real*8 :: start, finish
|
|
r0=0.0_8
|
|
r1=1.0_8
|
|
r2=2.0_8
|
|
!
|
|
!Input (free format)
|
|
!
|
|
! Title
|
|
! rexp (minus) potential exponent (1.0 for Coulomb)
|
|
! dx, dy, dz unit-cell dimensions (in bohr)
|
|
! nx, ny, nz unit-cell number
|
|
! nat number of centers in the unit cell
|
|
! loop over nat
|
|
! q x y z atom charge (in a.u.), atom position (in bohr)
|
|
! end loop
|
|
! fact a scale factor to divide the final energy
|
|
! END
|
|
!
|
|
read (*,*)
|
|
read (*,*) rexp
|
|
read (*,*) dx, dy, dz
|
|
read (*,*) nx, ny, nz
|
|
read (*,*) nat
|
|
allocate (q(1:nat),x(1:nat),y(1:nat),z(1:nat))
|
|
do iat=1,nat
|
|
read (5,*) q(iat), x(iat), y(iat), z(iat)
|
|
enddo
|
|
read (*,*) fact
|
|
!
|
|
call cpu_time(start)
|
|
!
|
|
qtot=r0
|
|
etot=r0
|
|
do iat=1,nat ! test cell
|
|
qtot=qtot+q(iat)
|
|
e0=r0
|
|
do ix=0,nx-1
|
|
e1=r0
|
|
do iy=0,ny-1
|
|
e2=r0
|
|
do iz=0,nz-1
|
|
e3=r0
|
|
do jat=1,nat ! general cell
|
|
xx=real(ix,8)+x(iat)-x(jat)
|
|
yy=real(iy,8)+y(iat)-y(jat)
|
|
zz=real(iz,8)+z(iat)-z(jat)
|
|
dd=dist(xx,yy,zz,nx,ny,nz)
|
|
! write (6,*) iat, jat, ix, iy, iz, dd
|
|
! if (dd < 1.q-10) cycle
|
|
if (ix==0.and.iy==0.and.iz==0.and.iat==jat) cycle
|
|
e3=e3+q(iat)*q(jat)/dd**rexp
|
|
enddo
|
|
e2=e2+e3
|
|
enddo
|
|
e1=e1+e2
|
|
enddo
|
|
e0=e0+e1
|
|
enddo
|
|
etot=etot+e0
|
|
enddo
|
|
!
|
|
call cpu_time(finish)
|
|
!
|
|
write (6,*)
|
|
write (6,101) qtot ! total charge - is a check: for neutral systems qtot should be zero
|
|
write (6,102) etot ! unit-cell energy / nat
|
|
write (6,103) etot/real(nat,8) ! unit-cell energy / (nat*fact)
|
|
write (6,104) etot/real(nat,8)/fact
|
|
write (6,*)
|
|
print '("Time = ",f6.3," seconds.")',finish-start
|
|
stop
|
|
101 format (f20.12,' total charge')
|
|
102 format (f20.12,' unit-cell energy')
|
|
103 format (f20.12,' per-atom cell energy')
|
|
104 format (f20.12,' per-atom cell energy / fact')
|
|
end program madelung
|
|
|
|
function dist(x,y,z,nx,ny,nz)
|
|
implicit none
|
|
integer :: nx, ny, nz
|
|
real*8 :: x, y, z
|
|
real*8 :: r1, pi
|
|
real*8 :: dist
|
|
r1=1.0_8
|
|
pi=acos(-r1)
|
|
dist=((real(nx,8)/pi*sin(pi*x/real(nx,8)))**2+(real(ny,8)/pi*sin(pi*y/real(ny,8)))**2 &
|
|
+(real(nz,8)/pi*sin(pi*z/real(nz,8)))**2)
|
|
dist=sqrt(dist)
|
|
return
|
|
end function dist
|