mirror of https://gitlab.com/QEF/q-e.git
420 lines
12 KiB
Fortran
420 lines
12 KiB
Fortran
!
|
|
! Copyright (C) 2002-2008 Quantum ESPRESSO group
|
|
! This file is distributed under the terms of the
|
|
! GNU General Public License. See the file `License'
|
|
! in the root directory of the present distribution,
|
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
|
!
|
|
! Written by Giovanni Bussi
|
|
! Adapted to QE by Andrea Ferretti & Layla Martin Samos
|
|
!
|
|
!----------------------------------
|
|
MODULE coulomb_vcut_module
|
|
!----------------------------------
|
|
!
|
|
IMPLICIT NONE
|
|
PRIVATE
|
|
|
|
!
|
|
! general purpose parameters
|
|
!
|
|
INTEGER, PARAMETER :: DP=KIND(1.0d0)
|
|
REAL(DP), PARAMETER :: PI = 3.14159265358979323846_DP
|
|
REAL(DP), PARAMETER :: TPI = 2.0_DP * pi
|
|
REAL(DP), PARAMETER :: FPI = 4.0_DP * pi
|
|
REAL(DP), PARAMETER :: e2 = 2.0_DP
|
|
REAL(DP), PARAMETER :: eps6 = 1.0E-6_DP
|
|
|
|
!
|
|
! definitions
|
|
!
|
|
TYPE vcut_type
|
|
REAL(DP) :: a(3,3)
|
|
REAL(DP) :: b(3,3)
|
|
REAL(DP) :: a_omega
|
|
REAL(DP) :: b_omega
|
|
REAL(DP), POINTER :: corrected(:,:,:)
|
|
REAL(DP) :: cutoff
|
|
LOGICAL :: orthorombic
|
|
END TYPE vcut_type
|
|
|
|
!
|
|
PUBLIC :: vcut_type
|
|
PUBLIC :: vcut_init
|
|
PUBLIC :: vcut_get
|
|
PUBLIC :: vcut_spheric_get
|
|
PUBLIC :: vcut_destroy
|
|
PUBLIC :: vcut_info
|
|
|
|
CONTAINS
|
|
|
|
!------------------------------------------
|
|
SUBROUTINE vcut_init(vcut,a,cutoff)
|
|
!------------------------------------------
|
|
!
|
|
TYPE(vcut_type), INTENT(OUT) :: vcut
|
|
REAL(DP), INTENT(IN) :: a(3,3)
|
|
REAL(DP), INTENT(IN) :: cutoff
|
|
|
|
INTEGER :: n1,n2,n3
|
|
INTEGER :: i1,i2,i3
|
|
INTEGER :: ierr
|
|
REAL(DP) :: q(3)
|
|
CHARACTER(9) :: subname='vcut_init'
|
|
REAL(DP) :: mod2a(3)
|
|
|
|
vcut%cutoff=cutoff
|
|
|
|
vcut%a=a
|
|
vcut%b= TPI * transpose(num_inverse(vcut%a))
|
|
vcut%b_omega=num_determinant(vcut%b)
|
|
vcut%a_omega=num_determinant(vcut%a)
|
|
|
|
! automatically finds whether the cell is orthorombic or not
|
|
vcut%orthorombic=.false.
|
|
!
|
|
mod2a=sqrt(sum(vcut%a**2,1))
|
|
if(abs(sum(vcut%a(:,1)*vcut%a(:,2)))/(mod2a(1)*mod2a(2))<eps6 .and. &
|
|
abs(sum(vcut%a(:,2)*vcut%a(:,3)))/(mod2a(2)*mod2a(3))<eps6 .and. &
|
|
abs(sum(vcut%a(:,3)*vcut%a(:,1)))/(mod2a(3)*mod2a(1))<eps6) vcut%orthorombic=.true.
|
|
!
|
|
if (.not.vcut%orthorombic) call errore(subname,"'vcut' Coulomb cutoff with non-orthogonal axis untested",1)
|
|
|
|
n1=ceiling(vcut%cutoff*sqrt(sum(vcut%a(1,:)**2))/(2.0*pi))
|
|
n2=ceiling(vcut%cutoff*sqrt(sum(vcut%a(2,:)**2))/(2.0*pi))
|
|
n3=ceiling(vcut%cutoff*sqrt(sum(vcut%a(3,:)**2))/(2.0*pi))
|
|
|
|
ALLOCATE(vcut%corrected(-n1:n1,-n2:n2,-n3:n3), STAT=ierr)
|
|
IF ( ierr/=0 ) CALL errore(subname,'allocating cvut%corrected',ABS(ierr))
|
|
!
|
|
vcut%corrected=0.0
|
|
|
|
!
|
|
! define the Fourier component of the modified Coulomb potential
|
|
!
|
|
DO i1=-n1,n1
|
|
DO i2=-n2,n2
|
|
DO i3=-n3,n3
|
|
!
|
|
q = MATMUL(vcut%b,(/i1,i2,i3/))
|
|
!
|
|
IF( SUM(q**2) > vcut%cutoff**2 ) CYCLE
|
|
!
|
|
vcut%corrected(i1,i2,i3) = &
|
|
vcut_formula(q,vcut%a,vcut%b,vcut%a_omega,vcut%orthorombic)
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
END SUBROUTINE vcut_init
|
|
|
|
!------------------------------------------
|
|
SUBROUTINE vcut_info(iun, vcut)
|
|
!------------------------------------------
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
INTEGER, INTENT(IN) :: iun
|
|
TYPE(vcut_type), INTENT(IN) :: vcut
|
|
!
|
|
INTEGER :: i, n(3)
|
|
!
|
|
IF ( ASSOCIATED( vcut%corrected ) ) THEN
|
|
!
|
|
DO i = 1, 3
|
|
n(i) = ( SIZE( vcut%corrected, i) -1 ) / 2
|
|
ENDDO
|
|
!
|
|
WRITE(iun, "( 2x,'Cutoff: ',f6.2,4x,' n grid: ',3i4,/)") vcut%cutoff, n(:)
|
|
!
|
|
ENDIF
|
|
!
|
|
END SUBROUTINE vcut_info
|
|
|
|
!------------------------------------------
|
|
SUBROUTINE vcut_destroy(vcut)
|
|
!------------------------------------------
|
|
!
|
|
TYPE(vcut_type), INTENT(INOUT) :: vcut
|
|
INTEGER :: ierr
|
|
!
|
|
DEALLOCATE(vcut%corrected, STAT=ierr)
|
|
IF ( ierr/=0 ) CALL errore('vcut_destroy','deallocating vcut',ABS(ierr))
|
|
!
|
|
END SUBROUTINE vcut_destroy
|
|
|
|
!------------------------------------------
|
|
FUNCTION vcut_get(vcut,q) RESULT(res)
|
|
!------------------------------------------
|
|
!
|
|
TYPE(vcut_type), INTENT(IN) :: vcut
|
|
REAL(DP), INTENT(IN) :: q(3)
|
|
REAL(DP) :: res
|
|
!
|
|
REAL(DP) :: i_real(3)
|
|
INTEGER :: i(3)
|
|
CHARACTER(8) :: subname='vcut_get'
|
|
!
|
|
i_real=(MATMUL(TRANSPOSE(vcut%a),q))/ TPI
|
|
i=NINT(i_real)
|
|
!
|
|
! internal check
|
|
IF( SUM( (i-i_real)**2 ) > eps6 ) &
|
|
CALL errore(subname,'q vector out of the grid',10)
|
|
!
|
|
IF( SUM(q**2) > vcut%cutoff**2 ) THEN
|
|
!
|
|
! usual form of Coulomb potential
|
|
res = FPI * e2 / SUM(q**2)
|
|
!
|
|
ELSE
|
|
!
|
|
IF( i(1)>ubound(vcut%corrected,1) .OR. i(1)<lbound(vcut%corrected,1) .OR. &
|
|
i(2)>ubound(vcut%corrected,2) .OR. i(2)<lbound(vcut%corrected,2) .OR. &
|
|
i(3)>ubound(vcut%corrected,3) .OR. i(3)<lbound(vcut%corrected,3)) THEN
|
|
CALL errore(subname,'index out of bound', 10)
|
|
ENDIF
|
|
!
|
|
res=vcut%corrected(i(1),i(2),i(3))
|
|
!
|
|
ENDIF
|
|
!
|
|
END FUNCTION vcut_get
|
|
|
|
!------------------------------------------
|
|
FUNCTION vcut_spheric_get(vcut,q) RESULT(res)
|
|
!------------------------------------------
|
|
!
|
|
TYPE(vcut_type), INTENT(IN) :: vcut
|
|
REAL(DP), INTENT(IN) :: q(3)
|
|
REAL(DP) :: res
|
|
!
|
|
REAl(DP) :: a(3,3), Rcut, kg2
|
|
LOGICAL :: limit
|
|
!
|
|
!
|
|
a = vcut%a
|
|
!
|
|
Rcut=0.5*minval(sqrt(sum(a**2,1)))
|
|
Rcut=Rcut-Rcut/50.0
|
|
limit=.false.
|
|
kg2=sum(q**2)
|
|
if(kg2<eps6) then
|
|
limit=.true.
|
|
endif
|
|
if(.not.limit) then
|
|
res=FPI*e2/kg2*(1.0-cos(Rcut*sqrt(kg2)))
|
|
else
|
|
res=FPI*e2*Rcut**2/2.0
|
|
endif
|
|
!
|
|
END FUNCTION vcut_spheric_get
|
|
|
|
!---------------------------------------------------------
|
|
FUNCTION vcut_formula(q,a,b,a_omega,orthorombic) result(res)
|
|
!---------------------------------------------------------
|
|
!
|
|
! Define the FT of the Coulomb potential according to the
|
|
! current lattice.
|
|
!
|
|
REAL(DP), INTENT(IN) :: q(3)
|
|
REAL(DP), INTENT(IN) :: a(3,3)
|
|
REAL(DP), INTENT(IN) :: b(3,3)
|
|
REAL(DP), INTENT(IN) :: a_omega
|
|
LOGICAL, INTENT(IN) :: orthorombic
|
|
REAL(DP) :: res
|
|
!
|
|
real(dp) :: rwigner
|
|
real(dp) :: sigma
|
|
|
|
rwigner=0.5*sqrt(1.0/maxval(sum(b**2,1)))*2*pi
|
|
|
|
!
|
|
! 3.0 is set to give a short range contribution inside the WS cell
|
|
!
|
|
sigma=3.0/rwigner
|
|
|
|
! compute longrange and shortrange contributions
|
|
res=vcut_formula_longrange(q,a,b,a_omega,sigma,6.0D0,orthorombic) &
|
|
+vcut_formula_shortrange(q,sigma)
|
|
|
|
END FUNCTION vcut_formula
|
|
|
|
!---------------------------------------------------------
|
|
FUNCTION vcut_formula_longrange(q,a,b,a_omega,sigma,security,orthorombic) result(res)
|
|
!---------------------------------------------------------
|
|
! compute the longrange contribution
|
|
real(dp), intent(in) :: q(3)
|
|
real(dp), intent(in) :: a(3,3)
|
|
real(dp), intent(in) :: b(3,3)
|
|
real(dp), intent(in) :: a_omega
|
|
real(dp), intent(in) :: sigma
|
|
real(dp), intent(in) :: security ! it determines the grid for the real-space sum; a reasonable value is 4.0
|
|
logical, intent(in) :: orthorombic
|
|
real(dp) :: res
|
|
integer :: n1,n2,n3
|
|
integer :: i1,i2,i3
|
|
real(dp) :: d1,d2,d3,weight,factor
|
|
real(dp) :: r(3),r2,modr
|
|
logical :: n1_is_even,n1_is_odd
|
|
real(dp) :: tmp, rtmp(3)
|
|
logical, parameter :: shifted=.false.
|
|
integer :: n1max
|
|
real(dp) :: i1_real,i2_real,i3_real
|
|
|
|
n1=security*sqrt(sum(a(:,1)**2))*sigma
|
|
n2=security*sqrt(sum(a(:,2)**2))*sigma
|
|
n3=security*sqrt(sum(a(:,3)**2))*sigma
|
|
|
|
n1_is_even=(n1/2)*2==n1
|
|
n1_is_odd=.not.n1_is_even
|
|
|
|
d1=1.0/real(n1,dp)
|
|
d2=1.0/real(n2,dp)
|
|
d3=1.0/real(n3,dp)
|
|
res=0.0
|
|
weight=a_omega*d1*d2*d3
|
|
! the only symmetry which can be used for any value of q is inversion
|
|
! NON-SHIFTED:
|
|
! if n1 is even: loop between 0 and n1/2, with weight=2.0 for all points except 0 and n1max
|
|
! if n2 is odd: loop between 0 and (n1+1)/2, with weight=2.0 for all points except 0
|
|
! SHIFTED:
|
|
! if n1 is even: loop between 0 and n1/2-1, with weight=2.0 for all points
|
|
! if n2 is odd: loop between 0 and (n1+1)/2, with weight=2.0 for all points except n1max
|
|
|
|
if(shifted)then
|
|
if(n1_is_even) n1max=n1/2-1
|
|
if(n1_is_odd) n1max=(n1+1)/2-1
|
|
else
|
|
if(n1_is_even) n1max=n1/2
|
|
if(n1_is_odd) n1max=(n1+1)/2
|
|
end if
|
|
do i1=0,n1max
|
|
factor=2.0
|
|
if(shifted) then
|
|
if(n1_is_odd .and. i1==n1max) factor=1.0
|
|
else
|
|
if(i1==0) factor=1.0
|
|
if(n1_is_even .and. i1==n1max) factor=1.0
|
|
end if
|
|
i1_real=i1
|
|
if(shifted) i1_real=i1_real+0.5
|
|
do i2=0,n2-1
|
|
i2_real=i2
|
|
if(shifted) i2_real=i2_real+0.5
|
|
do i3=0,n3-1
|
|
i3_real=i3
|
|
if(shifted) i3_real=i3_real+0.5
|
|
rtmp=matmul(a,(/i1_real*d1,i2_real*d2,i3_real*d3/))
|
|
r=vcut_minimal_image(a,b,rtmp,orthorombic)
|
|
r2=sum(r**2)
|
|
modr=sqrt(r2)
|
|
if(modr*sigma<eps6) then
|
|
tmp=e2*sqrt(2.0/pi)*sigma
|
|
else
|
|
tmp=e2*ERF(sigma*sqrt(0.5)*modr)/modr
|
|
end if
|
|
res=res+weight*factor*tmp*cos(sum(r*q))
|
|
end do
|
|
end do
|
|
end do
|
|
END FUNCTION vcut_formula_longrange
|
|
|
|
!---------------------------------------------------------
|
|
FUNCTION vcut_formula_shortrange(q,sigma) result(res)
|
|
!---------------------------------------------------------
|
|
real(dp), intent(in) :: q(3)
|
|
real(dp), intent(in) :: sigma
|
|
real(dp) :: res
|
|
if(sum(q**2/(sigma*sigma))<eps6) then
|
|
! analytic limit for small q
|
|
res=e2*tpi/(sigma*sigma)
|
|
else
|
|
res=e2*fpi/sum(q**2)*(1-exp(-0.5*sum(q**2)/(sigma*sigma)))
|
|
end if
|
|
END FUNCTION vcut_formula_shortrange
|
|
|
|
!---------------------------------------------------------
|
|
FUNCTION vcut_minimal_image(a,b,r,orthorombic) result(res)
|
|
!---------------------------------------------------------
|
|
real(dp), intent(in) :: a(3,3)
|
|
real(dp), intent(in) :: b(3,3)
|
|
real(dp), intent(in) :: r(3)
|
|
logical, intent(in) :: orthorombic
|
|
real(dp) :: res(3)
|
|
real(dp) :: r_minimal(3)
|
|
real(dp) :: r2_minimal
|
|
real(dp) :: r_try(3)
|
|
real(dp) :: r2_try
|
|
real(dp) :: r_components(3)
|
|
integer :: i1,i2,i3
|
|
integer, parameter :: max_displacement=1
|
|
if(orthorombic) then
|
|
! NINT ALGORITHM FOR ORTHOROMBIC CELL
|
|
r_components=(matmul(transpose(b),r))/(2.0*pi)
|
|
r_components=r_components-nint(r_components)
|
|
r_minimal=matmul(a,r_components)
|
|
else
|
|
! POOR MAN ALGORITHM FOR GENERIC CELL
|
|
r_minimal=r
|
|
r2_minimal=sum(r_minimal**2)
|
|
! loop over the possible neighbours
|
|
do i1=-max_displacement,max_displacement
|
|
do i2=-max_displacement,max_displacement
|
|
do i3=-max_displacement,max_displacement
|
|
if(i1==0 .and. i2==0 .and. i3==0) cycle
|
|
r_try=r+matmul(a,(/i1,i2,i3/))
|
|
r2_try=sum(r_try**2)
|
|
if(r2_try<r2_minimal) then
|
|
r2_minimal=r2_try
|
|
r_minimal=r_try
|
|
endif
|
|
end do
|
|
end do
|
|
end do
|
|
end if
|
|
res=r_minimal
|
|
END FUNCTION vcut_minimal_image
|
|
|
|
|
|
|
|
!************************************
|
|
!** tools from sax
|
|
|
|
function num_inverse(a) result(inv)
|
|
real(dp) :: inv(0:2,0:2)
|
|
real(dp), intent(in) :: a(0:2,0:2)
|
|
real(dp) :: tmp(0:2,0:2)
|
|
real(dp) :: det
|
|
real(dp),parameter :: eye3(3,3) = reshape((/ 1,0,0,0,1,0,0,0,1/),(/3,3/))
|
|
integer i,j
|
|
do i=0,2
|
|
do j=0,2
|
|
tmp(i,j) = a(modulo(i+1,3),modulo(j+1,3)) * a(modulo(i+2,3),modulo(j+2,3)) &
|
|
& - a(modulo(i+1,3),modulo(j+2,3)) * a(modulo(i+2,3),modulo(j+1,3))
|
|
end do
|
|
end do
|
|
det = num_determinant(a)
|
|
inv = transpose(tmp) / det
|
|
if(sum((matmul(inv,a))**2-eye3) > 1d-5) then
|
|
write(0,*) "AHIA",sum((matmul(inv,a)-eye3)**2)
|
|
write(0,*) "A",a
|
|
write(0,*) "inv",inv
|
|
write(0,*)">>", matmul(inv,a)
|
|
stop
|
|
end if
|
|
end function num_inverse
|
|
|
|
function num_determinant(a) result(det)
|
|
real(dp), intent(in) :: a(3,3)
|
|
real(dp) :: det
|
|
det = a(1,1)*a(2,2)*a(3,3) + a(1,2)*a(2,3)*a(3,1) + a(1,3)*a(2,1)*a(3,2) &
|
|
- a(1,1)*a(2,3)*a(3,2) - a(1,2)*a(2,1)*a(3,3) - a(1,3)*a(2,2)*a(3,1)
|
|
end function num_determinant
|
|
|
|
!** end tools from sax
|
|
!************************************
|
|
END MODULE coulomb_vcut_module
|
|
|