TABLE OF CONTENTS

ABINIT/m_vcoul/m_cutoff_surface [ Modules ]

[ Top ] [ Modules ]

NAME

  m_cutoff_surface

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_surface
16 
17  use defs_basis
18  use m_abicore
19  use m_errors
20 
21  implicit none
22 
23  private

m_vcoul/cutoff_surface [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 cutoff_surface

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

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

OUTPUT

  vc_cut(ng,nq)=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 such as 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

 73 subroutine cutoff_surface(nq,qpt,ng,gvec,gprimd,rcut,boxcenter,pdir,alpha,vc_cut,method)
 74 
 75 !Arguments ------------------------------------
 76 !scalars
 77  integer,intent(in) :: method,ng,nq
 78  real(dp),intent(in) :: rcut
 79 !arrays
 80  integer,intent(in) :: gvec(3,ng),pdir(3)
 81  real(dp),intent(in) :: alpha(3),boxcenter(3),gprimd(3,3),qpt(3,nq)
 82  real(dp),intent(out) :: vc_cut(ng,nq)
 83 
 84 !Local variables-------------------------------
 85 !scalars
 86  integer :: ig,iq,igs
 87  real(dp),parameter :: SMALL=tol4  !@WC: was tol6
 88  real(dp) :: qpg2,qpg_para,qpg_perp
 89  character(len=500) :: msg
 90 !arrays
 91  real(dp) :: b1(3),b2(3),b3(3),gcart(3),qc(3),qpg(3)
 92  real(dp),allocatable :: qcart(:,:)
 93 
 94 ! *************************************************************************
 95 
 96  ABI_UNUSED(pdir)
 97  ABI_UNUSED(boxcenter)
 98 
 99  ! === From reduced to cartesian coordinates ===
100  b1(:)=two_pi*gprimd(:,1)
101  b2(:)=two_pi*gprimd(:,2)
102  b3(:)=two_pi*gprimd(:,3)
103  ABI_MALLOC(qcart,(3,nq))
104  do iq=1,nq
105    qcart(:,iq) = b1*qpt(1,iq) + b2*qpt(2,iq) + b3*qpt(3,iq)
106  end do
107 
108  ! === Different approaches according to method ===
109  vc_cut(:,:)=zero
110 
111  SELECT CASE (method)
112 
113  CASE (1) ! 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(qcart(3,:))>SMALL)) then
117      write(std_out,*)qcart(:,:)
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_surface.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    ! * Here I suppose that all the input q-points are different from zero
127    do iq=1,nq
128      qc(:)=qcart(:,iq)
129      igs=1; if (SQRT(DOT_PRODUCT(qc,qc))<tol16) igs=2 ! avoid (q=0, G=0)
130      do ig=igs,ng
131        gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
132        qpg(:)=qc(:)+gcart(:)
133        qpg2  =DOT_PRODUCT(qpg(:),qpg(:))
134        qpg_para=SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp=qpg(3)
135        ! if (abs(qpg_perp)<SMALL.and.qpg_para<SMALL) stop 'SMALL in cutoff_surface
136        vc_cut(ig,iq)=four_pi/qpg2*(one-EXP(-qpg_para*rcut)*COS(qpg_perp*rcut))
137      end do
138    end do
139 
140  CASE (2) ! Rozzi"s method
141    ABI_ERROR("Work in progress")
142    ABI_UNUSED(alpha) ! just to keep alpha as an argument
143    !alpha=?? ; ap1sqrt=SQRT(one+alpha**2)
144    do iq=1,nq
145      qc(:)=qcart(:,iq)
146      do ig=1,ng
147        gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
148        qpg(:)=qc(:)+gcart(:)
149        qpg2  =DOT_PRODUCT(qpg(:),qpg(:))
150        qpg_para=SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp =qpg(3)
151        if (qpg_para>SMALL) then
152         vc_cut(ig,iq)=four_pi/qpg2*(one+EXP(-qpg_para*rcut)*(qpg_perp/qpg_para*SIN(qpg_perp*rcut)-COS(qpg_perp*rcut)))
153        else
154          if (ABS(qpg_perp)>SMALL) then
155            vc_cut(ig,iq)=four_pi/qpg_perp**2*(one-COS(qpg_perp*rcut)-qpg_perp*rcut*SIN(qpg_perp*rcut)) !&
156 !&          + 8*rcut*SIN(qpg_perp*rcut)/qpg_perp*LOG((alpha+ap1sqrt)*(one+ap1sqrt)/alpha) ! contribution due to finite surface
157          else
158            vc_cut(ig,iq)=-two_pi*rcut**2
159          end if
160        end if
161      end do !ig
162    end do !iq
163 
164  CASE DEFAULT
165    write(msg,'(a,i3)')' Wrong value of method: ',method
166    ABI_BUG(msg)
167  END SELECT
168 
169  ABI_FREE(qcart)
170 
171 end subroutine cutoff_surface