TABLE OF CONTENTS
ABINIT/barevcoul [ Functions ]
NAME
barevcoul
FUNCTION
Compute bare coulomb term in G-space on the FFT mesh i.e. 4pi/(G+q)**2
INPUTS
qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). gsqcut=cutoff value on G**2 for sphere inside fft box. (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2)) icutcoul=Option for the Coulomb potential cutoff technique divgq0= value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq. Used if q = Gamma gmet(3,3)=metrix tensor in G space in Bohr**-2. izero=if 1, unbalanced components of V(q,g) are set to zero # Used by the PAW library nfft=Total number of FFT grid points. ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft comm=MPI communicator.
OUTPUT
barev(nfft)=4pi/(G+q)**2, G=0 component is set to divgq0/pi if q = Gamma.
NOTES
This routine operates on the full FFT mesh. DO NOT PASS MPI_TYPE One can easily implemente MPI-FFT by just calling this routine and then extracting the G-vectors treated by the node.
SOURCE
143 subroutine barevcoul(rcut,qphon,gsqcut,gmet,nfft,nkpt_bz,ngfft,ucvol,barev,shortrange) 144 145 !Arguments ------------------------------------ 146 !scalars 147 integer,intent(in) :: nfft,nkpt_bz 148 real(dp),intent(in) :: rcut,gsqcut,ucvol 149 logical,intent(in),optional:: shortrange 150 !arrays 151 integer,intent(in) :: ngfft(18) 152 integer :: ng!!!! 153 real(dp),intent(in) :: qphon(3) 154 real(dp),intent(inout) :: gmet(3,3) 155 real(dp),intent(inout) :: barev(nfft) 156 !real(dp) :: a1(3),a2(3),a3(3) 157 real(dp) :: b1(3),b2(3),b3(3),rmet(3,3) !,gprimd(3,3), 158 type(dataset_type) :: dtset 159 type(MPI_type) :: mpi_enreg !!!! 160 type(crystal_t) :: Cryst !!!! 161 !type(gsphere_t) :: Gsph 162 type(vcut_t) :: vcut !!!! 163 !Local variables------------------------------- 164 !scalars 165 integer,parameter :: empty(3,3)=zero 166 integer :: comm 167 integer :: i1,i2,i23,i3,id1,id2,id3,icutcoul_local 168 integer :: ig,ig1min,ig1max,ig2min,ig2max,ig3min,ig3max 169 integer :: ii,ing,n1,n2,n3,npar,npt 170 integer :: opt_cylinder,opt_slab,test 171 real(dp),parameter :: tolfix=1.000000001e0_dp ! Same value as the one used in hartre 172 real(dp) :: check,step 173 real(dp) :: cutoff,gqg2p3,gqgm12,gqgm13,gqgm23,gs2,gs3,divgq0,rcut0 174 real(dp) :: bz_plane,dx,integ,q0_vol,q0_volsph 175 character(len=500) :: msg 176 !arrays 177 integer :: id(3), gamma_pt(3,1) 178 real(dp),allocatable :: gq(:,:),gpq(:),gpq2(:) 179 real(dp),allocatable :: vcfit(:,:),xx(:),yy(:) 180 real(dp),allocatable :: cov(:,:),par(:),qfit(:,:),sigma(:),var(:),qcart(:,:) 181 ! 182 comm=mpi_enreg%comm_world 183 ! 184 ! === Save dimension and other useful quantities in vcut% === 185 vcut%nfft = PRODUCT(ngfft(1:3)) ! Number of points in the FFT mesh. 186 ! ng and gvec are not used yet we don't want to use them without being defined 187 ng = -1 188 vcut%ng = ng ! Number of G-vectors in the Coulomb matrix elements. 189 vcut%rcut = rcut ! Cutoff radius for cylinder. 190 vcut%hcyl = zero ! Length of finite cylinder (Rozzi"s method, default is Beigi). 191 vcut%ucvol = ucvol ! Unit cell volume. 192 193 vcut%rprimd = Cryst%rprimd(:,:) ! Dimensional direct lattice. 194 vcut%boxcenter = dtset%boxcenter ! boxcenter at the moment is supposed to be at the origin. 195 vcut%vcutgeo = dtset%vcutgeo(:) ! Info on the orientation and extension of the cutoff region. 196 ! 197 ! === Define geometry and cutoff radius (if used) === 198 vcut%mode='NONE' 199 !icutcoul_local=dtset%icutcoul 200 201 ! BG: Temporary to circumvent the tests 202 if(shortrange) then 203 icutcoul_local=5 204 else 205 icutcoul_local=0 206 end if 207 ! ------------------------------------- 208 209 if (icutcoul_local==0) vcut%mode='SPHERE' 210 if (icutcoul_local==1) vcut%mode='CYLINDER' 211 if (icutcoul_local==2) vcut%mode='SLAB' 212 if (icutcoul_local==4) vcut%mode='ERF' 213 if (icutcoul_local==5) vcut%mode='ERFC' 214 ! 215 ! Treatment of the divergence at q+g=zero 216 rcut0= (three*nkpt_bz*ucvol/four_pi)**(one/three) 217 divgq0= two_pi*rcut0**two 218 219 !Initialize a few quantities 220 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 221 cutoff=gsqcut*tolfix 222 barev=zero 223 224 !In order to speed the routine, precompute the components of g+q 225 !Also check if the booked space was large enough... 226 227 ABI_MALLOC(gq,(3,max(n1,n2,n3))) 228 ABI_MALLOC(gpq,(nfft)) 229 ABI_MALLOC(gpq2,(nfft)) 230 231 do ii=1,3 232 id(ii)=ngfft(ii)/2+2 233 do ing=1,ngfft(ii) 234 ig=ing-(ing/id(ii))*ngfft(ii)-1 235 gq(ii,ing)=ig+qphon(ii) 236 end do 237 end do 238 ig1max=-1;ig2max=-1;ig3max=-1 239 ig1min=n1;ig2min=n2;ig3min=n3 240 241 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 242 243 do i3=1,n3 244 ! Precompute some products that do not depend on i2 and i1 245 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 246 gqgm23=gq(3,i3)*gmet(2,3)*2 247 gqgm13=gq(3,i3)*gmet(1,3)*2 248 do i2=1,n2 249 i23=n1*(i2-1 +(n2)*(i3-1)) 250 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 251 gqgm12=gq(2,i2)*gmet(1,2)*2 252 gqg2p3=gqgm13+gqgm12 253 do i1=1,n1 254 ii=i1+i23 255 gpq(ii)=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 256 if(gpq(ii)>=tol4) then 257 gpq2(ii) = piinv/gpq(ii) 258 end if 259 end do 260 end do 261 end do 262 263 ! Old version of the code extracted from m_Fock 264 ! do ig=1,nfft 265 ! if(abs(gpq(ig))<tol4) then 266 ! barev(ig)=barev(ig)+divgq0 267 ! else if(gpq(ig)<=cutoff) then 268 ! if(shortrange) then 269 ! barev(ig)=barev(ig)+gpq2(ig)*(one-exp(-pi/(gpq2(ig)*rcut**2))) 270 ! else 271 ! barev(ig)=barev(ig)+gpq2(ig)*(one-cos(rcut*sqrt(four_pi/gpq2(ig)))) 272 ! end if 273 ! end if 274 ! end do 275 276 barev(:)=zero 277 278 ! MG: This triggers SIGFPE as cryst is not initialized 279 !a1=Cryst%rprimd(:,1); b1=two_pi*gprimd(:,1) 280 !a2=Cryst%rprimd(:,2); b2=two_pi*gprimd(:,2) 281 !a3=Cryst%rprimd(:,3); b3=two_pi*gprimd(:,3) 282 283 SELECT CASE (TRIM(vcut%mode)) 284 CASE('SPHERE') ! Spencer-Alavi method 285 286 do ig=1,nfft 287 if(abs(gpq(ig))<tol4) then 288 barev(ig)=barev(ig)+divgq0 289 else if(gpq(ig)<=cutoff) then 290 barev(ig)=barev(ig)+gpq2(ig)*(one-cos(rcut*sqrt(four_pi/gpq2(ig)))) 291 end if 292 end do 293 294 CASE('CYLINDER') 295 296 test=COUNT(ABS(vcut%vcutgeo)>tol6) 297 ABI_CHECK(test==1,'Wrong cutgeo for cylinder') 298 299 ! === Beigi method is the default one, i.e infinite cylinder of radius rcut === 300 ! * Negative values to use Rozzi method with finite cylinder of extent hcyl. 301 opt_cylinder=1; vcut%hcyl=zero; vcut%pdir(:)=0 302 do ii=1,3 303 check=vcut%vcutgeo(ii) 304 if (ABS(check)>tol6) then 305 vcut%pdir(ii)=1 306 if (check<zero) then ! use Rozzi's method. 307 vcut%hcyl=ABS(check)*SQRT(SUM(Cryst%rprimd(:,ii)**2)) 308 opt_cylinder=2 309 end if 310 end if 311 end do 312 313 test=COUNT(vcut%pdir==1) 314 ABI_CHECK((test==1),'Wrong pdir for cylinder') 315 if (vcut%pdir(3)/=1) then 316 ABI_ERROR("The cylinder must be along the z-axis") 317 end if 318 319 ABI_BUG("cutoff cylinder API has changed!") 320 321 ! call cutoff_cylinder(nfft,gq,ng,Gsph%gvec,vcut%rcut,vcut%hcyl,vcut%pdir,& 322 !& vcut%boxcenter,Cryst%rprimd,barev,opt_cylinder,comm) 323 324 ! === If Beigi, treat the limit q--> 0 === 325 if (opt_cylinder==1) then 326 npar=8; npt=100 ; gamma_pt=RESHAPE((/0,0,0/),(/3,1/)) 327 ABI_MALLOC(qfit,(3,npt)) 328 ABI_MALLOC(vcfit,(1,npt)) 329 if (nfft==1) then 330 ABI_ERROR("nfft == 1 not supported when Beigi's method is used") 331 endif 332 qfit(:,:)=zero 333 step=half/(npt*(nfft-1)) ; qfit(3,:)=arth(tol6,step,npt) 334 335 !call cutoff_cylinder(npt,qfit,1,gamma_pt,vcut%rcut,vcut%hcyl,vcut%pdir,& 336 ! vcut%boxcenter,Cryst%rprimd,vcfit,opt_cylinder,comm) 337 338 ABI_MALLOC(xx,(npt)) 339 ABI_MALLOC(yy,(npt)) 340 ABI_MALLOC(sigma,(npt)) 341 ABI_MALLOC(par,(npar)) 342 ABI_MALLOC(var,(npar)) 343 ABI_MALLOC(cov,(npar,npar)) 344 do ii=1,npt 345 xx(ii)=normv(qfit(:,ii),gmet,'G') 346 end do 347 ABI_FREE(qfit) 348 sigma=one ; yy(:)=vcfit(1,:) 349 ABI_FREE(vcfit) 350 351 bz_plane=l2norm(b1.x.b2) 352 dx=(xx(2)-xx(1)) 353 integ=yy(2)*dx*3.0/2.0 354 355 do ii=3,npt-2 356 integ=integ+yy(ii)*dx 357 end do 358 integ=integ+yy(npt-1)*dx*3.0/2.0 359 write(std_out,*)' simple integral',integ 360 q0_volsph=(two_pi)**3/(nkpt_bz*ucvol) 361 q0_vol=bz_plane*two*xx(npt) 362 write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol 363 vcut%i_sz=bz_plane*two*integ/q0_vol 364 write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds) 365 write(std_out,*)' Cylindrical cutoff value ',vcut%i_sz 366 367 ABI_FREE(xx) 368 ABI_FREE(yy) 369 ABI_FREE(sigma) 370 ABI_FREE(par) 371 ABI_FREE(var) 372 ABI_FREE(cov) 373 374 else 375 ! In Rozzi"s method the lim q+G --> 0 is finite. 376 vcut%i_sz=barev(1) 377 end if 378 379 CASE('SLAB') 380 381 test=COUNT(vcut%vcutgeo/=zero) 382 ABI_CHECK(test==2,"Wrong vcutgeo") 383 ! 384 ! Two methods available 385 ! 386 ! === Default is Beigi"s method === 387 opt_slab=1; vcut%alpha(:)=zero 388 if (ANY(vcut%vcutgeo<zero)) opt_slab=2 389 vcut%pdir(:)=zero 390 do ii=1,3 391 check=vcut%vcutgeo(ii) 392 if (ABS(check)>zero) then ! Use Rozzi"s method with a finite slab along x-y 393 vcut%pdir(ii)=1 394 if (check<zero) vcut%alpha(ii)=normv(check*Cryst%rprimd(:,ii),rmet,'R') 395 end if 396 end do 397 398 ! Beigi"s method: the slab must be along x-y and R must be L_Z/2. 399 if (opt_slab==1) then 400 ABI_CHECK(ALL(vcut%pdir == (/1,1,0/)),"Surface must be in the x-y plane") 401 !vcut%rcut = half*SQRT(DOT_PRODUCT(a3,a3)) 402 end if 403 404 ABI_BUG("cutoff surface API has changed!") 405 !call cutoff_slab(nfft,gq,ng,Gsph%gvec,gprimd,vcut%rcut,& 406 ! vcut%boxcenter,vcut%pdir,vcut%alpha,barev,opt_slab) 407 408 ! 409 ! === If Beigi, treat the limit q--> 0 === 410 if (opt_slab==1) then 411 ! Integrate numerically in the plane close to 0 412 npt=100 ! Number of points in 1D 413 gamma_pt=RESHAPE((/0,0,0/),(/3,1/)) ! Gamma point 414 ABI_MALLOC(qfit,(3,npt)) 415 ABI_MALLOC(qcart,(3,npt)) 416 ABI_MALLOC(vcfit,(1,npt)) 417 if (nfft==1) then 418 ABI_ERROR("nfft == 1 not supported when Beigi's method is used") 419 endif 420 qfit(:,:)=zero 421 qcart(:,:)=zero 422 ! Size of the third vector 423 bz_plane=l2norm(b3) 424 q0_volsph=(two_pi)**3/(nkpt_bz*ucvol) 425 ! radius that gives the same volume as q0_volsph 426 ! Let's assume that c is perpendicular to the plane 427 ! We also assume isotropic BZ around gamma 428 step=sqrt((q0_volsph/bz_plane)/pi)/npt 429 430 !step=half/(npt*(Qmesh%nibz-1)) 431 ! Let's take qpoints along 1 line, the vcut does depend only on the norm 432 qcart(1,:)=arth(tol6,step,npt) 433 434 do ii = 1,npt 435 qfit(:,ii) = MATMUL(TRANSPOSE(Cryst%rprimd),qcart(:,ii))/(2*pi) 436 end do 437 438 ABI_BUG("cutoff surface API has changed!") 439 440 ! call cutoff_slab(npt,qfit,1,gamma_pt,gprimd,vcut%rcut,& 441 ! vcut%boxcenter,vcut%pdir,vcut%alpha,vcfit,opt_slab) 442 443 ABI_MALLOC(xx,(npt)) 444 ABI_MALLOC(yy,(npt)) 445 ABI_MALLOC(sigma,(npt)) 446 do ii=1,npt 447 !xx(ii)=qfit(1,:) 448 xx(ii)=normv(qfit(:,ii),gmet,'G') 449 end do 450 ABI_FREE(qfit) 451 sigma=one 452 yy(:)=vcfit(1,:) 453 !yy(:)=one 454 ABI_FREE(vcfit) 455 dx=(xx(2)-xx(1)) 456 ! integ = \int dr r f(r) 457 integ=xx(2)*yy(2)*dx*3.0/2.0 458 do ii=3,npt-2 459 integ=integ+xx(ii)*yy(ii)*dx 460 end do 461 integ=integ+xx(npt-1)*yy(npt-1)*dx*3.0/2.0 462 write(std_out,*)' simple integral',integ 463 q0_vol=bz_plane*pi*xx(npt)**2 464 write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol 465 vcut%i_sz=bz_plane*2*pi*integ/q0_vol 466 write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds) 467 write(std_out,*)' Cylindrical cutoff value ',vcut%i_sz 468 469 ABI_FREE(xx) 470 ABI_FREE(yy) 471 else 472 ! In Rozzi"s method the lim q+G --> 0 is finite. 473 vcut%i_sz=barev(1) 474 end if 475 476 CASE('ERF') 477 478 do ig=1,nfft 479 if(abs(gpq(ig))<tol4) then 480 barev(ig)=barev(ig)+divgq0 481 else if(gpq(ig)<=cutoff) then 482 if(shortrange) then 483 barev(ig)=barev(ig)+gpq2(ig)*exp(-pi/(gpq2(ig)*rcut**2)) 484 end if 485 end if 486 end do 487 488 CASE('ERFC') 489 490 do ig=1,nfft 491 if(abs(gpq(ig))<tol4) then 492 barev(ig)=barev(ig)+divgq0 493 else if(gpq(ig)<=cutoff) then 494 if(shortrange) then 495 barev(ig)=barev(ig)+gpq2(ig)*(one-exp(-pi/(gpq2(ig)*rcut**2))) 496 end if 497 end if 498 end do 499 500 CASE DEFAULT 501 write(msg,'(3a)')'No cut-off applied to the Coulomb Potential.', ch10, & 502 'Either icutcoul value not allowed or not defined.' 503 ABI_WARNING(msg) 504 END SELECT 505 506 ABI_FREE(gq) 507 ABI_FREE(gpq) 508 ABI_FREE(gpq2) 509 510 end subroutine barevcoul
ABINIT/m_barevcoul [ Modules ]
NAME
m_barevcoul
FUNCTION
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group () This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_barevcoul 23 24 use defs_basis 25 use m_abicore 26 use m_dtset 27 use m_errors 28 use m_xmpi 29 30 !use m_fstrings, only : sjoin, itoa 31 use defs_abitypes, only : MPI_type 32 use m_numeric_tools, only : arth, l2norm, OPERATOR(.x.) 33 use m_geometry, only : normv 34 use m_crystal, only : crystal_t 35 !use m_gsphere, only : gsphere_t 36 37 ! Cut-off methods modules 38 !use m_cutoff_sphere, only : cutoff_sphere 39 !use m_cutoff_slab, only : cutoff_slab 40 !use m_cutoff_cylinder, only : cutoff_cylinder 41 42 implicit none 43 44 private
m_barevcoul/vcut_t [ Types ]
[ Top ] [ m_barevcoul ] [ Types ]
NAME
vcoul_t
FUNCTION
SOURCE
55 type,public :: vcut_t 56 57 integer :: nfft 58 ! Number of points in FFT grid 59 60 integer :: ng 61 ! Number of G-vectors 62 63 real(dp) :: alpha(3) 64 ! Lenght of the finite slab 65 66 real(dp) :: rcut 67 ! Cutoff radius 68 69 real(dp) :: i_sz 70 ! Value of the integration of the Coulomb singularity 4\pi/V_BZ \int_BZ d^3q 1/q^2 71 72 real(dp) :: hcyl 73 ! Length of the finite cylinder along the periodic dimension 74 75 real(dp) :: ucvol 76 ! Volume of the unit cell 77 78 character(len=50) :: mode 79 ! String defining the cutoff mode, possible values are: sphere,cylinder,slab,crystal 80 81 integer :: pdir(3) 82 ! 1 if the system is periodic along this direction 83 84 real(dp) :: boxcenter(3) 85 ! 1 if the point in inside the cutoff region 0 otherwise 86 ! Reduced coordinates of the center of the box (input variable) 87 88 real(dp) :: vcutgeo(3) 89 ! For each reduced direction gives the length of the finite system 90 ! 0 if the system is infinite along that particular direction 91 ! negative value to indicate that a finite size has to be used 92 93 real(dp) :: rprimd(3,3) 94 ! Lattice vectors in real space. 95 96 real(dp),allocatable :: qibz(:,:) 97 ! qibz(3,nqibz) 98 ! q-points in the IBZ. 99 100 real(dp),allocatable :: barev(:) 101 ! barev(nfft) 102 ! Bare Coulomb potential on the FFT grid 103 ! A cut might be applied. 104 105 end type vcut_t 106 107 public :: barevcoul