TABLE OF CONTENTS
ABINIT/cutoff_cylinder [ 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 ]
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