TABLE OF CONTENTS
ABINIT/termcutoff [ Functions ]
NAME
termcutoff
FUNCTION
Apply a cut-off term to the 1/G**2-like terms that appears throughout the code at the ground-state level as follows: Ewald, NC-PSP, Hartree.
INPUTS
gsqcut = cutoff on (k+G)^2 (bohr^-2) (sphere for density and potential) (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2)) icutcoul = Information about the cut-off ngfft(18) = Information on the (fine) FFT grid used for the density. nkpt = Number of k-points in the Brillouin zone rprimd(3,3)=dimensional primitive translations in real space (bohr) vcutgeo(3)= Info on the orientation and extension of the cutoff region.
OUTPUT
gcutoff = Cut-off term applied to 1/G**2 terms
NOTES
1. In order to incur minimal changes in some portions of the code where a cut-off is needed to be applied, one can work only with the cut-off part of the Coulomb potential, unlike what is done in barevcoul module. 2. Fock term has its own legacy cut-off for the moment.
SOURCE
119 subroutine termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo) 120 121 !Arguments ------------------------------------ 122 !scalars 123 integer,intent(in) :: icutcoul, nkpt 124 real(dp),intent(in) :: gsqcut,rcut 125 126 !arrays 127 integer,intent(in) :: ngfft(18) 128 real(dp),intent(in) :: rprimd(3,3),vcutgeo(3) 129 130 !Local variables------------------------------- 131 !scalars 132 integer,parameter :: N0=1000 133 integer,save :: enough 134 integer :: i1,i2,i23,i3,ierr,id(3),ii,ig,ing 135 integer :: c1,c2,opt_cylinder 136 integer :: n1,n2,n3,nfft 137 integer :: test,opt_slab !opt_cylinder 138 real(dp) :: alpha_fac, ap1sqrt, log_alpha 139 real(dp) :: cutoff,rcut_loc,rcut2,check,rmet(3,3) 140 real(dp) :: gvecg2p3,gvecgm12,gvecgm13,gvecgm23,gs2,gs3 141 real(dp) :: gcart_para,gcart_perp,gcart_x,gcart_y,gcart_z 142 real(dp) :: j0,j1,k0,k1 143 real(dp) :: quad,ucvol 144 real(dp) :: hcyl,hcyl2 145 real(dp),parameter :: tolfix=1.0000001_dp,tol999=999.0 146 character(len=50) :: mode 147 character(len=500) :: msg 148 ! type(gcut_t) :: gcut ! 149 150 !arrays 151 integer :: periodic_dir(3) 152 real(dp) :: a1(3),a2(3),a3(3),b1(3),b2(3),b3(3) 153 real(dp) :: gcart(3),gmet(3,3),gprimd(3,3) 154 real(dp) :: alpha(3) 155 real(dp),allocatable :: gvec(:,:),gpq(:),gpq2(:) 156 real(dp),allocatable,intent(inout) :: gcutoff(:) 157 158 ! === Save dimension and other useful quantities in vcut% === 159 ! gcut%nfft = PRODUCT(ngfft(1:3)) ! Number of points in the FFT mesh. 160 ! gcut%ucvol = ucvol ! Unit cell volume. 161 ! gcut%rprimd = rprimd(:,:) ! Dimensional direct lattice. 162 ! gcut%vcutgeo = vcutgeo(:) ! Info on the orientation and extension of the cutoff region. 163 ! 164 !Initialize a few quantities 165 cutoff=gsqcut*tolfix 166 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 167 nfft=n1*n2*n3 168 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 169 170 ! Initialize container 171 ABI_MALLOC(gvec,(3,MAX(n1,n2,n3))) 172 ABI_MALLOC(gpq,(nfft)) 173 ABI_MALLOC(gpq2,(nfft)) 174 ABI_MALLOC(gcutoff,(nfft)) 175 gcart(:) = zero ; gpq = zero ; gpq2 = zero ; gcutoff = zero 176 177 !In order to speed the routine, precompute the components of gvectors 178 !Also check if the booked space was large enough... 179 do ii=1,3 180 id(ii)=ngfft(ii)/2+2 181 do ing=1,ngfft(ii) 182 gvec(ii,ing)=ing-(ing/id(ii))*ngfft(ii)-1 183 end do 184 end do 185 186 ! Get the cut-off method info from the input file 187 ! Assign method to one of the available cases 188 mode='NONE' 189 190 if (icutcoul==0) mode='SPHERE' 191 if (icutcoul==1) mode='CYLINDER' 192 if (icutcoul==2) mode='SLAB' 193 if (icutcoul==3) mode='CRYSTAL' 194 if (icutcoul==4) mode='ERF' 195 if (icutcoul==5) mode='ERFC' 196 197 !Print in log info about the cut-off method at every call: 198 enough = enough + 1 199 if (enough < 5) then 200 write(msg,'(3a)')ch10,' 1/G**2 cut-off applied in the following step : cutoff-mode = ',TRIM(mode) 201 call wrtout(std_out,msg) 202 end if 203 !!! 204 205 do i3=1,n3 206 ! Precompute some products that do not depend on i2 and i1 207 gs3=gvec(3,i3)*gvec(3,i3)*gmet(3,3) 208 gvecgm23=gvec(3,i3)*gmet(2,3)*2 209 gvecgm13=gvec(3,i3)*gmet(1,3)*2 210 211 do i2=1,n2 212 i23=n1*(i2-1 + n2*(i3-1)) 213 gs2=gs3+ gvec(2,i2)*(gvec(2,i2)*gmet(2,2)+gvecgm23) 214 gvecgm12=gvec(2,i2)*gmet(1,2)*2 215 gvecg2p3=gvecgm13+gvecgm12 216 do i1=1,n1 217 ii=i1+i23 218 gpq(ii)=gs2+gvec(1,i1)*(gvec(1,i1)*gmet(1,1)+gvecg2p3) 219 if(gpq(ii)>=tol4) then 220 gpq2(ii) = piinv/gpq(ii) 221 end if 222 end do 223 end do 224 end do 225 226 SELECT CASE (TRIM(mode)) 227 228 CASE('SPHERE') ! Spencer-Alavi method 229 230 ! Calculate rcut for each method 231 if(rcut>tol4) then 232 rcut_loc = rcut 233 else 234 rcut_loc = (three*nkpt*ucvol/four_pi)**(one/three) 235 endif 236 237 do ig=1,nfft 238 if(abs(gpq(ig))<tol4) then 239 gcutoff(ig)=0.0 240 else 241 gcutoff(ig)=one-cos(rcut_loc*sqrt(four_pi/gpq2(ig))) 242 !gcutoff(ig)=one-cos(rcut_loc/sqrt(gpq2(ii))) 243 end if 244 end do 245 246 CASE('CYLINDER') 247 248 test=COUNT(vcutgeo/=zero) 249 ABI_CHECK(test==1,'Wrong cutgeo for cylinder') 250 251 ! === From reduced to Cartesian coordinates === 252 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 253 254 a1=rprimd(:,1); b1=two_pi*gprimd(:,1) 255 a2=rprimd(:,2); b2=two_pi*gprimd(:,2) 256 a3=rprimd(:,3); b3=two_pi*gprimd(:,3) 257 258 !ntasks=nfft 259 !call xmpi_split_work(ntasks,comm,my_start,my_stop) 260 261 !Calculate rcut for each method ! 262 ! 263 264 ! * Check if Bravais lattice is orthorombic and parallel to the Cartesian versors. 265 ! In this case the intersection of the W-S cell with the x-y plane is a rectangle with -ha_<=x<=ha_ and -hb_<=y<=hb_ 266 if ( (ANY(ABS(rprimd(2:3, 1))>tol6)).or.& 267 & (ANY(ABS(rprimd(1:3:2,2))>tol6)).or.& 268 & (ANY(ABS(rprimd(1:2, 3))>tol6)) & 269 & ) then 270 msg = ' Bravais lattice should be orthorombic and parallel to the cartesian versors ' 271 ABI_ERROR(msg) 272 end if 273 274 ! === Beigi method is the default one, i.e infinite cylinder of radius rcut === 275 ! * Negative values to use Rozzi method with finite cylinder of extent hcyl. 276 opt_cylinder=1; hcyl=zero; periodic_dir(:)=0 277 do ii=1,3 278 check=vcutgeo(ii) 279 if (ABS(check)>tol6) then 280 periodic_dir(ii)=1 281 if (check<zero) then ! use Rozzi's method. 282 hcyl=ABS(check)*SQRT(SUM(rprimd(:,ii)**2)) 283 opt_cylinder=2 284 !Check to enter the infinite Rozzi treatment 285 if(vcutgeo(3).le.-tol999) then 286 hcyl=tol12 287 end if 288 end if 289 end if 290 end do 291 292 ! Calculate rcut for each method 293 if(rcut>tol4) then 294 rcut_loc = rcut 295 else 296 rcut_loc = half*SQRT(DOT_PRODUCT(a1,a1)) 297 endif 298 299 if (opt_cylinder==1) then 300 ABI_CHECK(ALL(periodic_dir == (/0,0,1/)),"The cylinder must be along the z-axis") 301 end if 302 303 rcut_= rcut_loc 304 305 ! =================================================== 306 ! === Setup for the quadrature of matrix elements === 307 ! =================================================== 308 qopt_ =6 ! Quadrature method, see quadrature routine. 309 ntrial_ =30 ! Max number of attempts. 310 accuracy_=0.001 ! Fractional accuracy required. 311 npts_ =6 ! Initial number of point (only for Gauss-Legendre method). 312 hcyl_ =hcyl ! Lenght of cylinder along z, only if method==2 313 314 write(msg,'(3a,2(a,i5,a),a,f8.5)')ch10,& 315 ' cutoff_cylinder: Info on the quadrature method : ',ch10,& 316 ' Quadrature scheme = ',qopt_,ch10,& 317 ' Max number of attempts = ',ntrial_,ch10,& 318 ' Fractional accuracy = ',accuracy_ 319 call wrtout(std_out,msg) 320 321 SELECT CASE (opt_cylinder) 322 323 CASE(1) 324 325 ! === Infinite cylinder, interaction is zeroed outside the Wigner-Seitz cell === 326 ! * Beigi"s expression holds only if the BZ is sampled only along z. 327 write(msg,'(2(a,f8.4))')' cutoff_cylinder: Using Beigi''s Infinite cylinder ' 328 call wrtout(std_out,msg) 329 ! * Check if Bravais lattice is orthorombic and parallel to the Cartesian versors. 330 ! In this case the intersection of the W-S cell with the x-y plane is a rectangle with -ha_<=x<=ha_ and -hb_<=y<=hb_ 331 if ( (ANY(ABS(rprimd(2:3, 1))>tol6)).or.& 332 & (ANY(ABS(rprimd(1:3:2,2))>tol6)).or.& 333 & (ANY(ABS(rprimd(1:2, 3))>tol6)) & 334 & ) then 335 msg = ' Bravais lattice should be orthorhombic and parallel to the cartesian verctors ' 336 ABI_ERROR(msg) 337 end if 338 339 ha_=half*SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1))) 340 hb_=half*SQRT(DOT_PRODUCT(rprimd(:,2),rprimd(:,2))) 341 r0_=MIN(ha_,hb_)/N0 342 343 do i3=1,n3 344 do i2=1,n2 345 i23=n1*(i2-1 + n2*(i3-1)) 346 do i1=1,n1 347 ii=i1+i23 348 349 gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3) 350 gcartx_=gcart(1) ; gcarty_=gcart(2) ; gcart_para_=ABS(gcart(3)) 351 gpq(ii)=DOT_PRODUCT(gcart,gcart) 352 353 ! Avoid singularity in K_0{gcart_para_\rho) by using a small g along the periodic dimension. 354 if (gcart_para_<tol8) then 355 gcart_para_ = tol8 356 write(std_out,*)"setting gcart_para to=",gcart_para_ 357 end if 358 ! 359 ! * Calculate $ 2\int_{WS} dxdy K_0{gcart_para_\rho) cos(x.gcartx + y.gcarty) $ 360 ! where WS is the Wigner-Seitz cell. 361 !tmp=zero 362 ! === More stable method: midpoint integration with Romberg extrapolation === 363 call quadrature(K0cos_dy,zero,ha_,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 364 !write(std_out,'(i8,a,es14.6)')ii,' 3 ',quad 365 if (ierr/=0) then 366 ABI_ERROR("Accuracy not reached") 367 end if 368 ! === Store final result === 369 ! * Factor two comes from the replacement WS -> (1,4) quadrant thanks to symmetries of the integrad. 370 !tmp=tmp+quad 371 gcutoff(ii)=quad*gpq(ii)/pi 372 373 end do !i1 374 end do !i2 375 end do !i3 376 377 CASE(2) 378 379 ! === Finite cylinder of length hcyl, from Rozzi et al === 380 ! TODO add check on hcyl value that should be smaller that 1/deltaq 381 if (hcyl_<zero) then 382 write(msg,'(a,f8.4)')' Negative value for cylinder length hcyl=',hcyl_ 383 ABI_BUG(msg) 384 end if 385 386 if (ABS(hcyl_)>tol12) then 387 388 write(msg,'(2(a,f8.4))')' cutoff_cylinder: using finite cylinder of length= ',hcyl,' rcut= ',rcut_loc 389 call wrtout(std_out,msg) 390 hcyl_=hcyl 391 hcyl2=hcyl**2.0_dp 392 rcut2=rcut_loc**2.0_dp 393 394 do i3=1,n3 395 do i2=1,n2 396 i23=n1*(i2-1 + n2*(i3-1)) 397 do i1=1,n1 398 ii=i1+i23 399 400 gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3) 401 gcart_para_=ABS(gcart(3)) ; gcart_perp_=SQRT(gcart(1)**2+gcart(2)**2) 402 gpq(ii)=DOT_PRODUCT(gcart,gcart) 403 404 if (gcart_perp_/=zero.and.gcart_para_/=zero) then 405 call quadrature(F2,zero,rcut_loc,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 406 if (ierr/=0) then 407 ABI_ERROR("Accuracy not reached") 408 end if 409 410 gcutoff(ii)=quad*gpq(ii) 411 412 else if (gcart_perp_==zero.and.gcart_para_/=zero) then 413 414 ! $ \int_0^h sin(qpg_para_.z)/\sqrt(rcut^2+z^2)dz $ 415 call quadrature(F3,zero,hcyl,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 416 if (ierr/=0) then 417 ABI_ERROR("Accuracy not reached") 418 end if 419 420 c1=one/gcart_para_**2-COS(gcart_para_*hcyl_)/gcart_para_**2-hcyl_*SIN(gcart_para_*hcyl_)/gcart_para_ 421 c2=SIN(gcart_para_*hcyl_)*SQRT(hcyl2+rcut2) 422 gcutoff(ii)=(c1+(c2-quad)/gcart_para_)*gpq(ii) 423 424 else if (gcart_perp_/=zero.and.gcart_para_==zero) then 425 ! $ 4pi\int_0^rcut d\rho \rho J_o(qpg_perp_.\rho) ln((h+\sqrt(h^2+\rho^2))/\rho) $ 426 call quadrature(F4,zero,rcut_loc,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 427 if (ierr/=0) then 428 ABI_ERROR("Accuracy not reached") 429 end if 430 431 gcutoff(ii)=quad*gpq(ii) 432 433 else if (gcart_perp_==zero.and.gcart_para_==zero) then 434 ! Use lim q+G --> 0 435 gcutoff(ii)=zero 436 else 437 ABI_BUG('You should not be here!') 438 end if 439 440 end do !i1 441 end do !i2 442 end do !i3 443 444 else 445 446 call wrtout(std_out,'Using Rozzi infinite cut-off cylinder method.') 447 448 do i3=1,n3 449 do i2=1,n2 450 i23=n1*(i2-1 + n2*(i3-1)) 451 do i1=1,n1 452 ii=i1+i23 453 gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3) 454 gcart_x=gcart(1) ; gcart_y=gcart(2) ; gcart_z=ABS(gcart(3)) 455 gcart_perp_ = SQRT(gcart_x**2.0_dp+gcart_y**2.0_dp) ; 456 gpq(ii)=DOT_PRODUCT(gcart,gcart) 457 458 if (gcart_z>tol4) then 459 ! === Analytic expression === 460 call CALJY1(gcart_perp_*rcut_loc,j1,0) 461 call CALCK0(gcart_z*rcut_loc,k0,1) 462 call CALJY0(gcart_perp_*rcut_loc,j0,0) 463 call CALCK1(gcart_z*rcut_loc,k1,1) 464 gcutoff(ii)=one+rcut_loc*gcart_perp_*j1*k0-rcut_loc*gcart_z*j0*k1 465 else 466 if (gcart_perp_>tol4) then 467 ! === Integrate r*Jo(G_xy r)log(r) from 0 up to rcut_ === 468 call quadrature(F5,zero,rcut_loc,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 469 if (ierr/=0) then 470 ABI_ERROR("Accuracy not reached") 471 end if 472 gcutoff(ii)= -quad*gpq(ii) 473 else 474 gcutoff(ii)= zero !-pi*rcut_loc**2*(two*LOG(rcut_loc)-one) 475 end if 476 end if 477 478 end do !i1 479 end do !i2 480 end do !i3 481 end if ! case 2 - selecting Rozzi 482 483 CASE DEFAULT 484 ABI_BUG(sjoin('Wrong value for cylinder method:',itoa(opt_cylinder))) 485 END SELECT 486 487 CASE('SLAB') 488 489 test=COUNT(vcutgeo/=zero) 490 ABI_CHECK(test==2,"Wrong vcutgeo") 491 492 ! === From reduced to cartesian coordinates === 493 a1=rprimd(:,1); b1=two_pi*gprimd(:,1) 494 a2=rprimd(:,2); b2=two_pi*gprimd(:,2) 495 a3=rprimd(:,3); b3=two_pi*gprimd(:,3) 496 497 !SLAB Default - Beigi 498 opt_slab=1; alpha(:)=zero 499 ! Otherwise use Rozzi's method 500 if (ANY(vcutgeo<zero) .or. rcut>tol8) opt_slab=2 501 periodic_dir(:)=0 502 do ii=1,3 503 check=vcutgeo(ii) 504 if (ABS(check)>zero) then 505 periodic_dir(ii)=1 506 !For Rozzi"s method 507 if (check<zero) alpha(ii)=normv(check*rprimd(:,ii),rmet,'R') 508 end if 509 end do 510 511 SELECT CASE (opt_slab) 512 513 !CASE SLAB 1 - Beigi 514 CASE(1) 515 516 ! Calculate rcut for each method ! 517 rcut_loc = half*SQRT(DOT_PRODUCT(a3,a3)) 518 519 do i3=1,n3 520 do i2=1,n2 521 i23=n1*(i2-1 + n2*(i3-1)) 522 do i1=1,n1 523 ii=i1+i23 524 gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3) 525 gcart_para=SQRT(gcart(1)**2+gcart(2)**2) ; gcart_perp = gcart(3) 526 if(gcart_para<tol4.and.ABS(gcart_perp)<tol4) then 527 !if(gcart_para<tol12.and.ABS(gcart_perp)<tol12) then 528 gcutoff(ii)=zero 529 else 530 gcutoff(ii)=one-EXP(-gcart_para*rcut_loc)*COS(gcart_perp*rcut_loc) 531 end if 532 end do !i1 533 end do !i2 534 end do !i3 535 536 !CASE SLAB 2 - Rozzi 537 CASE(2) 538 539 !Set the cut-off radius 540 if(rcut>tol4) then 541 rcut_loc = rcut 542 else 543 rcut_loc = half*SQRT(DOT_PRODUCT(a3,a3)) 544 endif 545 546 !In the case of finite, Rozzi's method provide another parameter 547 !for the cut-off: alpha 548 !!! ATT: alpha = L_x/L_y --> in-plane geometry dependence 549 alpha_fac=SQRT(DOT_PRODUCT(a1,a1))/SQRT(DOT_PRODUCT(a2,a2)) 550 ap1sqrt=SQRT(one+alpha_fac**2) 551 log_alpha=LOG((alpha_fac+ap1sqrt)*(one+ap1sqrt)/alpha_fac) 552 553 do i3=1,n3 554 do i2=1,n2 555 i23=n1*(i2-1 + n2*(i3-1)) 556 do i1=1,n1 557 ii=i1+i23 558 gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3) 559 gcart_para=SQRT(gcart(1)**2+gcart(2)**2) ; gcart_perp = gcart(3) 560 if(gcart_para>tol4) then 561 gcutoff(ii)=one+EXP(-gcart_para*rcut_loc)*(gcart_perp/gcart_para*& 562 & SIN(gcart_perp*rcut_loc)-COS(gcart_perp*rcut_loc)) 563 else 564 if (ABS(gcart_perp)>tol4) then 565 gcutoff(ii)=one-COS(-gcart_perp*rcut_loc)-gcart_perp*rcut_loc*SIN(gcart_perp*rcut_loc) 566 ! gcutoff(ii)=one-COS(-gcart_perp*rcut_loc)-SIN(gcart_perp*rcut_loc) - Altered Rozzi's 567 else 568 gcutoff(ii)=zero 569 endif 570 endif 571 end do !i1 572 end do !i2 573 end do !i3 574 575 CASE DEFAULT 576 write(msg,'(a,i3)')' Wrong value of slab method: ',opt_slab 577 ABI_BUG(msg) 578 END SELECT 579 580 CASE('ERF') 581 582 ! Calculate rcut for each method ! Same as SPHERE 583 if(rcut>tol4) then 584 rcut_loc = rcut 585 else 586 rcut_loc= (three*nkpt*ucvol/four_pi)**(one/three) 587 endif 588 589 do ig=1,nfft 590 if(abs(gpq(ig))<tol4) then 591 gcutoff(ig)=zero ! @Gamma: initialize quantity in each requiered routine 592 else !if(gpq(ig)<=cutoff) then 593 gcutoff(ig)=exp(-pi/(gpq2(ig)*rcut_loc**2)) 594 end if 595 end do !ig 596 597 CASE('ERFC') 598 599 ! Calculate rcut for each method ! Same as SPHERE 600 if(rcut>tol4) then 601 rcut_loc = rcut 602 else 603 rcut_loc= (three*nkpt*ucvol/four_pi)**(one/three) 604 endif 605 606 do ig=1,nfft 607 if(abs(gpq(ig))<tol4) then 608 gcutoff(ig)=zero ! @Gamma: initialize quantity in each requiered routine 609 else 610 gcutoff(ig)=one-exp(-pi/(gpq2(ig)*rcut_loc**2)) 611 end if 612 end do !ig 613 614 CASE('CRYSTAL') 615 gcutoff(:)=one ! Neutral cut-off 616 !write(msg,'(a)')'CRYSTAL method: no cut-off applied to G**2 while CRYSTAL method is implied!' 617 !ABI_WARNING(msg) 618 CASE DEFAULT 619 gcutoff=one ! Neutral cut-off 620 !write(msg,'(a)')'No cut-off applied to G**2!' 621 !ABI_WARNING(msg) 622 END SELECT 623 624 ABI_FREE(gvec) 625 ABI_FREE(gpq) 626 ABI_FREE(gpq2) 627 ! ABI_FREE(gcutoff) 628 629 end subroutine termcutoff
m_gtermcutoff/gtermcut_t [ Types ]
NAME
gtermcut_t
FUNCTION
SOURCE
49 !!! type,public :: gtermcut_t 50 51 !!! integer :: nfft 52 !!! ! Number of points in FFT grid 53 54 !!! integer :: ng 55 !!! ! Number of G-vectors 56 57 !!! real(dp) :: ucvol 58 !!! ! Volume of the unit cell 59 60 !!! ! integer :: periodic_dir(3) 61 !!! ! 1 if the system is periodic along this direction 62 63 !!! ! real(dp) :: boxcenter(3) 64 !!! ! 1 if the point in inside the cutoff region 0 otherwise 65 !!! ! Reduced coordinates of the center of the box (input variable) 66 67 !!! real(dp) :: rprimd(3,3) 68 !!! ! Lattice vectors in real space. 69 70 !!! real(dp),allocatable :: gtermcuoff(:) 71 !!! ! gtermcuoff(nfft) 72 !!! ! G cut-off array on the FFT grid 73 74 !!! end type gtermcut_t 75 76 public :: termcutoff