TABLE OF CONTENTS


ABINIT/cutoff_cylinder [ Functions ]

[ Top ] [ Functions ]

NAME

 cutoff_cylinder

FUNCTION

  Calculate the Fourier components of an effective Coulomb interaction
  zeroed outside a finite cylindrical region. Two methods are implemented:

   method==1: The interaction in the (say) x-y plane is truncated outside the Wigner-Seitz
              cell centered on the wire in the x-y plane. The interaction has infinite
              extent along the z axis and the Fourier transform is singular only at the Gamma point.
              Only orthorombic Bravais lattices are supported.
   method==2: The interaction is truncated outside a cylinder of radius rcut. The cylinder has finite
              extent along z. No singularity occurs.

INPUTS

  boxcenter(3)= center of the wire in the x-y axis
  qpt(3)= q-point
  ng=number of G vectors
  gvec(3,ng)=G vectors in reduced coordinates
  rprimd(3,3)=dimensional real space primitive translations (bohr)
  method=1 for Beigi approach (infinite cylinder with interaction truncated outside the W-S cell)
         2 for Rozzi method (finite cylinder)
  comm=MPI communicator.

OUTPUT

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

SOURCE

 83 subroutine cutoff_cylinder(qpt, ng, gvec, rcut, hcyl, pdir, boxcenter, rprimd, vc_cut, method, comm)
 84 
 85 !Arguments ------------------------------------
 86 !scalars
 87  integer,intent(in) :: ng,method,comm
 88  real(dp),intent(in) :: rcut,hcyl
 89 !arrays
 90  integer,intent(in) :: gvec(3,ng),pdir(3)
 91  real(dp),intent(in) :: boxcenter(3),qpt(3),rprimd(3,3)
 92  real(dp),intent(out) :: vc_cut(ng)
 93 
 94 !Local variables-------------------------------
 95 !scalars
 96  integer,parameter :: N0=1000
 97  integer :: ig,igs,ierr, my_rank, nproc
 98  real(dp) :: j0,j1,k0,k1,qpg2,qpg_xy,tmp
 99  real(dp) :: qpg_z,quad,rcut2,hcyl2,c1,c2,ucvol,SMALL
100  logical :: q_is_zero
101  character(len=500) :: msg
102 !arrays
103  real(dp) :: qpg(3),b1(3),b2(3),b3(3),gmet(3,3),rmet(3,3),gprimd(3,3),qc(3),gcart(3)
104 !************************************************************************
105 
106  ABI_UNUSED(pdir)
107  ABI_UNUSED(boxcenter)
108 
109  ! ===================================================
110  ! === Setup for the quadrature of matrix elements ===
111  ! ===================================================
112  qopt_    =6         ! Quadrature method, see quadrature routine.
113  ntrial_  =30        ! Max number of attempts.
114  accuracy_=0.001     ! Fractional accuracy required.
115  npts_    =6         ! Initial number of point (only for Gauss-Legendre method).
116  SMALL    =tol4      ! Below this value (q+G)_i is treated as zero.
117  rcut_    =rcut      ! Radial cutoff, used only if method==2
118  hcyl_    =hcyl      ! Lenght of cylinder along z, only if method==2
119 
120  !write(msg,'(3a,2(a,i5,a),a,f8.5)')ch10,&
121  ! ' cutoff_cylinder: Info on the quadrature method : ',ch10,&
122  ! '  Quadrature scheme      = ',qopt_,ch10,&
123  ! '  Max number of attempts = ',ntrial_,ch10,&
124  ! '  Fractional accuracy    = ',accuracy_
125  !call wrtout(std_out, msg)
126 
127  ! From reduced to Cartesian coordinates.
128  call metric(gmet, gprimd, -1, rmet, rprimd, ucvol)
129  b1(:) =two_pi*gprimd(:,1)
130  b2(:) =two_pi*gprimd(:,2)
131  b3(:) =two_pi*gprimd(:,3)
132 
133  qc = b1(:)*qpt(1) + b2(:)*qpt(2) + b3(:)*qpt(3)
134 
135  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
136 
137  ! ================================================
138  ! === Different approaches according to method ===
139  ! ================================================
140  vc_cut = zero
141 
142  select case (method)
143 
144  case (1)
145    ! Infinite cylinder, interaction is zeroed outside the Wigner-Seitz cell.
146    ! NB: Beigi's expression holds only if the BZ is sampled only along z.
147    !call wrtout(std_out, 'cutoff_cylinder: Using Beigi''s Infinite cylinder')
148 
149    if (ANY(qc(1:2) > SMALL)) then
150      write(msg,'(5a)')&
151       ' found q-points with non zero components in the X-Y plane. ',ch10,&
152       ' This is not allowed, see Notes in cutoff_cylinder.F90. ',ch10,&
153       ' ACTION: Modify the q-point sampling. '
154      ABI_ERROR(msg)
155    end if
156 
157    ! Check if Bravais lattice is orthorombic and parallel to the Cartesian versors.
158    ! In this case the intersection of the WS cell with the x-y plane is a rectangle with -ha_<=x<=ha_ and -hb_<=y<=hb_
159    if ((ANY(ABS(rprimd(2:3,  1)) > tol6)) .or. &
160        (ANY(ABS(rprimd(1:3:2,2)) > tol6)) .or. &
161        (ANY(ABS(rprimd(1:2,  3)) > tol6))) then
162      ABI_ERROR('Bravais lattice should be orthorombic and parallel to the Cartesian versors')
163    end if
164 
165    ha_ = half*SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1)))
166    hb_ = half*SQRT(DOT_PRODUCT(rprimd(:,2),rprimd(:,2)))
167    r0_ = MIN(ha_,hb_)/N0
168 
169    ! For each (q,G) pair evaluate the integral defining the Coulomb cutoff.
170    ! NB: the code assumes that all q-vectors are non zero and q_xy/=0.
171    igs=1
172    ! Skip singularity at Gamma, it will be treated "by hand" in csigme.
173    q_is_zero = (normv(qpt, gmet, 'G') < tol4)
174 
175    do ig=igs,ng
176      if (mod(ig, nproc) /= my_rank) cycle ! MPI parallelism
177 
178      gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
179      qpg(:)=qc(:)+gcart(:)
180      qpgx_=qpg(1); qpgy_=qpg(2); qpg_para_=ABS(qpg(3))
181      !write(std_out,*)"qpgx_=",qpgx_, "qpgy_=",qpgy_, "qpg_para=",qpg_para_
182 
183      ! Avoid singularity in K_0{qpg_para_\rho) by using a small q along the periodic dimension.
184      if (q_is_zero .and. qpg_para_ < tol6) qpg_para_ = tol6
185 
186      ! Calculate $ 2\int_{WS} dxdy K_0{qpg_para_\rho) cos(x.qpg_x + y.qpg_y) $
187      ! where WS is the Wigner-Seitz cell.
188      tmp=zero
189 
190      ! Difficult part, integrate on a small cirle of radius r0 using spherical coordinates
191      !call quadrature(K0cos_dth_r0,zero,r0_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
192      !ABI_CHECK(ierr == 0, "Accuracy not reached")
193      !write(std_out,'(i8,a,es14.6)')ig,' 1 ',quad
194      !tmp=tmp+quad
195      ! Add region with 0<=x<=r0 and y>=+-(SQRT(r0^2-x^2))since WS is rectangular
196      !call quadrature(K0cos_dy_r0,zero,r0_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
197      !ABI_CHECK(ierr == 0, "Accuracy not reached")
198      !write(std_out,'(i8,a,es14.6)')ig,' 2 ',quad
199      !tmp=tmp+quad
200      ! Get the in integral in the rectangle with x>=r0, should be the easiest but sometimes has problems to converge
201      !call quadrature(K0cos_dy,r0_,ha_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
202      !ABI_CHECK(ierr == 0, "Accuracy not reached")
203      !write(std_out,'(i8,a,es14.6)')ig,' 3 ',quad
204      !
205      ! More stable method: midpoint integration with Romberg extrapolation ===
206      call quadrature(K0cos_dy,zero,ha_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
207      !write(std_out,'(i8,a,es14.6)')ig,' 3 ',quad
208      ABI_CHECK(ierr == 0, "Accuracy not reached in quadrature!")
209 
210      ! Store final result
211      ! Factor two comes from the replacement WS -> (1,4) quadrant thanks to symmetries of the integrad.
212      tmp = tmp+quad
213      vc_cut(ig) = two*(tmp*two)
214    end do ! ig
215 
216  case (2)
217    ! Finite cylinder of length hcyl from Rozzi et al.
218    ! TODO add check on hcyl value that should be smaller that 1/deltaq
219    if (hcyl_ < zero) then
220      write(msg,'(a,f8.4)')' Negative value for cylinder length hcyl_=',hcyl_
221      ABI_BUG(msg)
222    end if
223 
224    if (ABS(hcyl_) > tol12) then
225      !write(std_out,'(2(a,f8.4))')' cutoff_cylinder: using finite cylinder of length= ',hcyl_,' rcut= ',rcut_
226      hcyl2=hcyl_**2
227      rcut2=rcut_**2
228 
229      ! No singularity occurs in finite cylinder, thus start from 1.
230      do ig=1,ng
231        if (mod(ig, nproc) /= my_rank) cycle ! MPI parallelism
232 
233        gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
234        qpg(:)=qc(:)+gcart(:)
235        qpg_para_=ABS(qpg(3)) ; qpg_perp_=SQRT(qpg(1)**2+qpg(2)**2)
236 
237        if (qpg_perp_ /= zero .and. qpg_para_ /= zero) then
238          ! $ 4\pi\int_0^{R_c} d\rho\rho j_o(qpg_perp_.\rho)\int_0^hcyl dz\cos(qpg_para_*z)/sqrt(\rho^2+z^2) $
239          call quadrature(F2,zero,rcut_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
240          ABI_CHECK(ierr == 0, "Accuracy not reached")
241          vc_cut(ig) = four_pi*quad
242 
243        else if (qpg_perp_ == zero .and. qpg_para_ /= zero) then
244          ! $ \int_0^h sin(qpg_para_.z)/\sqrt(rcut^2+z^2)dz $
245          call quadrature(F3,zero,hcyl_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
246          ABI_CHECK(ierr == 0, "Accuracy not reached")
247 
248          c1=one/qpg_para_**2-COS(qpg_para_*hcyl_)/qpg_para_**2-hcyl_*SIN(qpg_para_*hcyl_)/qpg_para_
249          c2=SIN(qpg_para_*hcyl_)*SQRT(hcyl2+rcut2)
250          vc_cut(ig) = four_pi*c1+four_pi*(c2-quad)/qpg_para_
251 
252        else if (qpg_perp_ /= zero .and. qpg_para_ == zero) then
253          ! $ 4pi\int_0^rcut d\rho \rho J_o(qpg_perp_.\rho) ln((h+\sqrt(h^2+\rho^2))/\rho) $
254          call quadrature(F4,zero,rcut_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
255          ABI_CHECK(ierr == 0, "Accuracy not reached")
256          vc_cut(ig) = four_pi*quad
257 
258        else if (qpg_perp_ == zero .and. qpg_para_ == zero) then
259          ! Use lim q+G --> 0
260          vc_cut(ig) = two_pi*(-hcyl2+hcyl_*SQRT(hcyl2+rcut2)+rcut2*LOG((hcyl_+SQRT(hcyl_+SQRT(hcyl2+rcut2)))/rcut_))
261 
262        else
263          ABI_BUG('You should not be here!')
264        end if
265 
266      end do !ig
267 
268    else
269      ! Infinite cylinder.
270      !call wrtout(std_out, ' cutoff_cylinder: using Rozzi''s method with infinite cylinder ')
271 
272      do ig=1,ng
273        if (mod(ig, nproc) /= my_rank) cycle ! MPI parallelism
274 
275        gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
276        qpg(:)=qc(:)+gcart(:)
277        qpg2  =DOT_PRODUCT(qpg,qpg)
278        qpg_z =ABS(qpg(3)) ; qpg_xy=SQRT(qpg(1)**2+qpg(2)**2)
279 
280        if (qpg_z > SMALL) then
281          ! Analytic expression.
282          call CALCK0(qpg_z *rcut_, k0, 1)
283          call CALJY1(qpg_xy*rcut_, j1, 0)
284          call CALJY0(qpg_xy*rcut_, j0, 0)
285          call CALCK1(qpg_z *rcut_, k1, 1)
286          vc_cut(ig) = (four_pi/qpg2)*(one+rcut_*qpg_xy*j1*k0-qpg_z*rcut_*j0*k1)
287        else
288          if (qpg_xy > SMALL) then
289            ! Integrate r*Jo(G_xy r)log(r) from 0 up to rcut_
290            call quadrature(F5,zero,rcut_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
291            ABI_CHECK(ierr == 0, "Accuracy not reached")
292            vc_cut(ig)=-four_pi*quad
293          else
294            ! Analytic expression
295            vc_cut(ig)=-pi*rcut_**2*(two*LOG(rcut_)-one)
296          end if
297        end if
298      end do ! ig
299    end if !finite/infinite
300 
301  case default
302    ABI_BUG(sjoin('Wrong value for method:',itoa(method)))
303  end select
304 
305  ! Collect vc_cut on each core
306  call xmpi_sum(vc_cut, comm, ierr)
307 
308 end subroutine cutoff_cylinder

ABINIT/m_cutoff_cylinder [ Modules ]

[ Top ] [ Modules ]

NAME

  m_cutoff_cylinder

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_cylinder
16 
17  use defs_basis
18  use m_abicore
19  use m_errors
20  use m_xmpi
21  use m_splines
22  use m_sort
23 
24  use m_fstrings,        only : sjoin, itoa
25  use m_geometry,        only : normv, metric
26  use m_bessel,          only : CALJY0, CALJY1, CALCK0, CALCK1
27  use m_numeric_tools,   only : OPERATOR(.x.), quadrature
28  use m_paw_numeric,     only : paw_jbessel
29 
30  implicit none
31 
32  private