TABLE OF CONTENTS


ABINIT/m_vcoul/m_cutoff_slab [ Modules ]

[ Top ] [ Modules ]

NAME

  m_cutoff_slab

FUNCTION

SOURCE

 9 #if defined HAVE_CONFIG_H
10 #include "config.h"
11 #endif
12 
13 #include "abi_common.h"
14 
15 module m_cutoff_slab
16 
17  use defs_basis
18  use m_abicore
19  use m_errors
20 
21  use m_fstrings, only : sjoin, itoa
22 
23  implicit none
24 
25  private

m_vcoul/cutoff_slab [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 cutoff_slab

FUNCTION

  Calculate the Fourier components of an effective Coulomb interaction
  within a slab of thickness 2*rcut which is symmetric with respect to the xy plane.
  In this implementation rcut=L_z/2 where L_z is the periodicity along z

INPUTS

  qpt(3)=q-point
  ng=Number of G vectors.
  gvec(3,ng)=G vectors in reduced coordinates.
  gprimd(3,3)=Dimensional primitive translations in reciprocal space ($\textrm{bohr}^{-1}$).
  gmet(3,3)=Metric in reciprocal space.

OUTPUT

  vc_cut(ng)=Fourier components of the effective Coulomb interaction.

NOTES

  The Fourier expression for an interaction truncated along the z-direction (i.e non-zero only if |z|<R) is:

  vc(q.G) = 4pi/|q+G|^2 * [ 1 + e^{-((q+G)_xy)*R} * ( (q_z+G_z)/(q+G)_xy * sin((q_z+G_z)R) -
   - cos((q_z+G_Z)R)) ]  (1)

  Equation (1) diverges when q_xy+G_xy --> 0 for any non zero q_z+G_z
  However if we choose R=L/2, where L defines the periodicity along z,
  and we limit ourselves to consider q-points with q_z==0, then sin((q_z+G_z)R)=sin(G_Z 2pi/L)=0 for every G.
  Under these assumptions we obtain

  v(q,G) = 4pi/|q+G|^2 [1-e^{-(q+G)_xy*L/2}\cos((q_z+G_z)R)]

  which is always finite when G_z /=0 while it diverges as 4piR/(q+G)_xy as (q+G)_xy -->0
  but only in the x-y plane.

SOURCE

 75 subroutine cutoff_slab(qpt, ng, gvec, gprimd, rcut, boxcenter, pdir, alpha, vc_cut, method)
 76 
 77 !Arguments ------------------------------------
 78 !scalars
 79  integer,intent(in) :: method,ng
 80  real(dp),intent(in) :: rcut
 81 !arrays
 82  integer,intent(in) :: gvec(3,ng),pdir(3)
 83  real(dp),intent(in) :: alpha(3),boxcenter(3),gprimd(3,3),qpt(3)
 84  real(dp),intent(out) :: vc_cut(ng)
 85 
 86 !Local variables-------------------------------
 87 !scalars
 88  integer :: ig,igs
 89  real(dp),parameter :: SMALL=tol4  !@WC: was tol6
 90  real(dp) :: qpg2,qpg_para,qpg_perp
 91  character(len=500) :: msg
 92 !arrays
 93  real(dp) :: b1(3),b2(3),b3(3),gcart(3),qc(3),qpg(3)
 94 
 95 ! *************************************************************************
 96 
 97  ABI_UNUSED(pdir)
 98  ABI_UNUSED(boxcenter)
 99 
100  ! From reduced to cartesian coordinates.
101  b1(:)=two_pi*gprimd(:,1)
102  b2(:)=two_pi*gprimd(:,2)
103  b3(:)=two_pi*gprimd(:,3)
104 
105  qc = b1*qpt(1) + b2*qpt(2) + b3*qpt(3)
106 
107  ! Different approaches according to method
108  vc_cut = zero
109 
110  select case (method)
111 
112  case (1)
113    ! Beigi's expression.
114    ! q-points with non-zero component along the z-axis are not allowed if
115    ! the simplified Eq.1 for the Coulomb interaction is used.
116    if (ANY(ABS(qc) > SMALL)) then
117      write(std_out,*)qc
118      write(msg,'(5a)')&
119       'Found q-points with non-zero component along non-periodic direction ',ch10,&
120       'This is not allowed, see Notes in cutoff_slab.F90 ',ch10,&
121       'ACTION: Modify the q-point sampling '
122      ABI_ERROR(msg)
123    end if
124 
125    ! Calculate truncated Coulomb interaction for a infinite surface
126    ! supposing input q-points are different from zero.
127    igs=1; if (SQRT(DOT_PRODUCT(qc,qc))<tol16) igs=2 ! avoid (q=0, G=0)
128    do ig=igs,ng
129      gcart(:) = b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
130      qpg(:) = qc(:) + gcart(:)
131      qpg2  = DOT_PRODUCT(qpg(:),qpg(:))
132      qpg_para = SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp=qpg(3)
133      vc_cut(ig) = four_pi/qpg2*(one-EXP(-qpg_para*rcut)*COS(qpg_perp*rcut))
134    end do
135 
136  case (2)
137    ! Rozzi's method
138    ABI_ERROR("Work in progress")
139    ABI_UNUSED(alpha) ! just to keep alpha as an argument
140    !alpha=?? ; ap1sqrt=SQRT(one+alpha**2)
141    do ig=1,ng
142      gcart(:) = b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
143      qpg(:) = qc(:) + gcart(:)
144      qpg2  =DOT_PRODUCT(qpg(:),qpg(:))
145      qpg_para=SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp =qpg(3)
146      if (qpg_para>SMALL) then
147       vc_cut(ig)=four_pi/qpg2*(one+EXP(-qpg_para*rcut)*(qpg_perp/qpg_para*SIN(qpg_perp*rcut)-COS(qpg_perp*rcut)))
148      else
149        if (ABS(qpg_perp)>SMALL) then
150          vc_cut(ig)=four_pi/qpg_perp**2*(one-COS(qpg_perp*rcut)-qpg_perp*rcut*SIN(qpg_perp*rcut)) ! &
151          ! contribution due to finite slab
152          ! + 8*rcut*SIN(qpg_perp*rcut)/qpg_perp*LOG((alpha+ap1sqrt)*(one+ap1sqrt)/alpha)
153        else
154          vc_cut(ig)=-two_pi*rcut**2
155        end if
156      end if
157    end do !ig
158 
159  case default
160    ABI_BUG(sjoin('Wrong value for method:', itoa(method)))
161  end select
162 
163 end subroutine cutoff_slab