TABLE OF CONTENTS


ABINIT/m_fft_mesh [ Modules ]

[ Top ] [ Modules ]

NAME

  m_fft_mesh

FUNCTION

  This module contains routines and helper functions to perform the setup of the FFT mesh
  It also provides a set of tools to test the grid, rotate the mesh according to the symmetry
  operations of the space group etc.

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (MG, XG, GMR, VO, LR, RWG, YMN, RS, TR, DC)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

23 #if defined HAVE_CONFIG_H
24 #include "config.h"
25 #endif
26 
27 #include "abi_common.h"
28 
29 MODULE m_fft_mesh
30 
31  use defs_basis
32  use m_errors
33  use m_abicore
34  use m_hide_blas
35 
36  use defs_fftdata,     only : size_goed_fft
37  use m_numeric_tools,  only : denominator, mincm, iseven, pfactorize
38  use m_symtk,          only : mati3inv
39  use m_geometry,       only : xred2xcart
40  use m_crystal,        only : crystal_t
41 
42  implicit none
43 
44  private
45 
46  public :: setmesh             ! Perform the setup of the FFT mesh for the GW oscillator strengths.
47  public :: check_rot_fft       ! Test whether the mesh is compatible with the rotational part of the space group.
48  public :: fft_check_rotrans   ! Test whether the mesh is compatible with the symmetries of the space group.
49  public :: rotate_fft_mesh     ! Calculate the FFT index of the rotated mesh.
50  public :: cigfft              ! Calculate the FFT index of G-G0.
51  public :: ig2gfft             ! Returns the component of a G in the FFT Box from its sequential index.
52  public :: g2ifft              ! Returns the index of the G in the FFT box from its reduced coordinates.
53  public :: get_gftt            ! Calculate the G"s in the FFT box from ngfft
54  public :: calc_ceigr          ! e^{iG.r} on the FFT mesh (complex valued).
55  public :: calc_eigr           ! e^{iG.r} on the FFT mesh (version for real array with RE,IM).
56  public :: calc_ceikr          ! e^{ik.r} on the FFT mesh (complex valued).
57  public :: times_eigr          ! Multiply an array on the real-space mesh by e^{iG0.r}
58  public :: times_eikr          ! Multiply an array on the real-space mesh by e^{ik.r}
59  public :: phase               ! Compute ph(ig)=$\exp(\pi\ i \ n/ngfft)$ for n=0,...,ngfft/2,-ngfft/2+1,...,-1
60  public :: mkgrid_fft          !  It sets the grid of fft (or real space) points to be treated.
61 
62  interface calc_ceigr
63    module procedure calc_ceigr_spc
64    module procedure calc_ceigr_dpc
65  end interface calc_ceigr

ABINIT/mkgrid_fft [ Functions ]

[ Top ] [ Functions ]

NAME

  mkgrid_fft

FUNCTION

  It sets the grid of fft (or real space) points to be treated.

INPUTS

OUTPUT

PARENTS

      mkcore_paw,mklocl_realspace

CHILDREN

      xred2xcart

SOURCE

1774 subroutine mkgrid_fft(ffti3_local,fftn3_distrib,gridcart,nfft,ngfft,rprimd)
1775 
1776 
1777 !This section has been created automatically by the script Abilint (TD).
1778 !Do not modify the following lines by hand.
1779 #undef ABI_FUNC
1780 #define ABI_FUNC 'mkgrid_fft'
1781 !End of the abilint section
1782 
1783  implicit none
1784 
1785 !Arguments ------------------------------------
1786  integer, intent(in) :: nfft
1787  integer,intent(in) :: ngfft(18)
1788  integer, dimension(*), intent(in) :: ffti3_local,fftn3_distrib
1789  real(dp), dimension(3,nfft), intent(out) :: gridcart
1790  real(dp),intent(in) :: rprimd(3,3)
1791 
1792 !Local variables-------------------------------
1793  integer :: ind,i1,i2,i3,i3loc,me,nproc
1794  integer :: n1,n2,n3
1795  real(dp), dimension(3) :: coord
1796  real(dp), dimension(3,nfft) :: gridred
1797 
1798 ! *************************************************************************
1799 
1800  n1    = ngfft(1)
1801  n2    = ngfft(2)
1802  n3    = ngfft(3)
1803  nproc = ngfft(10)
1804  me    = ngfft(11)
1805 
1806  do i3 = 1, n3, 1
1807    if(fftn3_distrib(i3) == me) then !MPI
1808      i3loc=ffti3_local(i3)
1809      coord(3) = real(i3 - 1, dp) / real(n3, dp)
1810      do i2 = 1, n2, 1
1811        coord(2) = real(i2 - 1, dp) / real(n2, dp)
1812        do i1 = 1, n1, 1
1813          ind=i1+(i2-1)*n1+(i3loc-1)*n1*n2
1814          coord(1) = real(i1 - 1, dp) / real(n1, dp)
1815          gridred(:, ind) = coord(:)
1816        end do
1817      end do
1818    end if
1819  end do
1820  call xred2xcart(nfft, rprimd, gridcart, gridred)
1821 
1822 end subroutine mkgrid_fft

m_fft_mesh/calc_ceigr_dpc [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 calc_ceigr_dpc

FUNCTION

  Helper function to calculate e^{iG.r} on the FFT mesh.

INPUTS

  gg(3)=G vector in reduced coordinates.
  nfft=Total number of points in the FFT mesh.
  nspinor=Number of spinors
  ngfft(18)=information about 3D FFT,

OUTPUT

  ceigr(nfft*nspinor)=e^{ik.r} on the FFT mesh.

PARENTS

CHILDREN

      xcopy

SOURCE

1341 subroutine calc_ceigr_dpc(gg,nfft,nspinor,ngfft,ceigr)
1342 
1343 
1344 !This section has been created automatically by the script Abilint (TD).
1345 !Do not modify the following lines by hand.
1346 #undef ABI_FUNC
1347 #define ABI_FUNC 'calc_ceigr_dpc'
1348 !End of the abilint section
1349 
1350  implicit none
1351 
1352 !Arguments ------------------------------------
1353 !scalars
1354  integer,intent(in) :: nfft,nspinor
1355 !arrays
1356  integer,intent(in) :: gg(3)
1357  integer,intent(in) :: ngfft(18)
1358  complex(dpc),intent(out) :: ceigr(nfft*nspinor)
1359 
1360 !Local variables-------------------------------
1361 !scalars
1362  integer :: ix,iy,iz,fft_idx,base,isp
1363  real(dp) :: gdotr
1364 
1365 ! *************************************************************************
1366 
1367  if (ALL(gg==0)) then
1368    ceigr=cone; RETURN
1369  end if
1370 
1371  fft_idx=0
1372  do iz=0,ngfft(3)-1
1373    do iy=0,ngfft(2)-1
1374      do ix=0,ngfft(1)-1
1375        gdotr= two_pi*( gg(1)*(ix/DBLE(ngfft(1))) &
1376 &                     +gg(2)*(iy/DBLE(ngfft(2))) &
1377 &                     +gg(3)*(iz/DBLE(ngfft(3))) )
1378        fft_idx = fft_idx+1
1379        ceigr(fft_idx)=DCMPLX(DCOS(gdotr),DSIN(gdotr))
1380      end do
1381    end do
1382  end do
1383 
1384  if (nspinor > 1) then
1385    do isp=2,nspinor
1386      base = 1 + (isp-1)*nfft
1387      call xcopy(nfft,ceigr,1,ceigr(base:),1)
1388    end do
1389  end if
1390 
1391 end subroutine calc_ceigr_dpc

m_fft_mesh/calc_ceigr_spc [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 calc_ceigr_spc

FUNCTION

  Helper function to calculate e^{iG.r} on the FFT mesh.

INPUTS

  gg(3)=G vector in reduced coordinates.
  nfft=Total number of points in the FFT mesh.
  nspinor=Number of spinors
  ngfft(18)=information about 3D FFT,

OUTPUT

  ceigr(nfft*nspinor)=e^{ik.r} on the FFT mesh.

PARENTS

CHILDREN

      xcopy

SOURCE

1262 subroutine calc_ceigr_spc(gg,nfft,nspinor,ngfft,ceigr)
1263 
1264 
1265 !This section has been created automatically by the script Abilint (TD).
1266 !Do not modify the following lines by hand.
1267 #undef ABI_FUNC
1268 #define ABI_FUNC 'calc_ceigr_spc'
1269 !End of the abilint section
1270 
1271  implicit none
1272 
1273 !Arguments ------------------------------------
1274 !scalars
1275  integer,intent(in) :: nfft,nspinor
1276 !arrays
1277  integer,intent(in) :: gg(3)
1278  integer,intent(in) :: ngfft(18)
1279  complex(spc),intent(out) :: ceigr(nfft*nspinor)
1280 
1281 !Local variables-------------------------------
1282 !scalars
1283  integer :: ix,iy,iz,fft_idx,base,isp
1284  real(dp) :: gdotr
1285 
1286 ! *************************************************************************
1287 
1288  if (ALL(gg==0)) then
1289    ceigr=(1._sp,0._sp)
1290    RETURN
1291  end if
1292 
1293  fft_idx=0
1294  do iz=0,ngfft(3)-1
1295    do iy=0,ngfft(2)-1
1296      do ix=0,ngfft(1)-1
1297        gdotr= two_pi*( gg(1)*(ix/DBLE(ngfft(1))) &
1298 &                     +gg(2)*(iy/DBLE(ngfft(2))) &
1299 &                     +gg(3)*(iz/DBLE(ngfft(3))) )
1300        fft_idx = fft_idx+1
1301        ceigr(fft_idx)=CMPLX(DCOS(gdotr),DSIN(gdotr), KIND=spc)
1302      end do
1303    end do
1304  end do
1305 
1306  if (nspinor > 1) then
1307    do isp=2,nspinor
1308      base = 1 + (isp-1)*nfft
1309      call xcopy(nfft,ceigr,1,ceigr(base:),1)
1310    end do
1311  end if
1312 
1313 end subroutine calc_ceigr_spc

m_fft_mesh/calc_ceikr [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 calc_ceikr

FUNCTION

  Helper function to calculate e^{ik.r} on the FFT mesh.

INPUTS

  kk(3)=k-point in reduced coordinates.
  nfft=Total number of points in the FFT mesh.
  ngfft(18)=information about 3D FFT,

OUTPUT

  ceikr(nfft)=e^{ik.r} on the FFT mesh.

PARENTS

      m_wfd

CHILDREN

SOURCE

1491 pure subroutine calc_ceikr(kk,nfft,ngfft,ceikr)
1492 
1493 
1494 !This section has been created automatically by the script Abilint (TD).
1495 !Do not modify the following lines by hand.
1496 #undef ABI_FUNC
1497 #define ABI_FUNC 'calc_ceikr'
1498 !End of the abilint section
1499 
1500  implicit none
1501 
1502 !Arguments ------------------------------------
1503 !scalars
1504  integer,intent(in) :: nfft
1505 !arrays
1506  real(dp),intent(in) :: kk(3)
1507  integer,intent(in) :: ngfft(18)
1508  complex(dpc),intent(out) :: ceikr(nfft)
1509 
1510 !Local variables-------------------------------
1511 !scalars
1512  integer :: ix,iy,iz,fft_idx
1513  real(dp) :: kdotr
1514 
1515 ! *************************************************************************
1516 
1517  !if (ALL(ABS(kk<tol12)) then
1518  !  ceikr=cone; RETURN
1519  !end if
1520 
1521  fft_idx=0
1522  do iz=0,ngfft(3)-1
1523    do iy=0,ngfft(2)-1
1524      do ix=0,ngfft(1)-1
1525        kdotr= two_pi*( kk(1)*(ix/DBLE(ngfft(1))) &
1526 &                     +kk(2)*(iy/DBLE(ngfft(2))) &
1527 &                     +kk(3)*(iz/DBLE(ngfft(3))) )
1528        fft_idx = fft_idx+1
1529        ceikr(fft_idx)=DCMPLX(DCOS(kdotr),DSIN(kdotr))
1530      end do
1531    end do
1532  end do
1533 
1534 end subroutine calc_ceikr

m_fft_mesh/calc_eigr [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 calc_eigr

FUNCTION

  Helper function to calculate e^{iG.r} on the FFT mesh.

INPUTS

  gg(3)=G vector in reduced coordinates.
  nfft=Total number of points in the FFT mesh.
  ngfft(18)=information about 3D FFT,

OUTPUT

  eigr(2*nfft)=e^{ig.r} on the FFT mesh.

PARENTS

      m_fft_prof

CHILDREN

SOURCE

1418 pure subroutine calc_eigr(gg,nfft,ngfft,eigr)
1419 
1420 
1421 !This section has been created automatically by the script Abilint (TD).
1422 !Do not modify the following lines by hand.
1423 #undef ABI_FUNC
1424 #define ABI_FUNC 'calc_eigr'
1425 !End of the abilint section
1426 
1427  implicit none
1428 
1429 !Arguments ------------------------------------
1430 !scalars
1431  integer,intent(in) :: nfft
1432 !arrays
1433  integer,intent(in) :: gg(3)
1434  integer,intent(in) :: ngfft(18)
1435  real(dp),intent(out) :: eigr(2*nfft)
1436 
1437 !Local variables-------------------------------
1438 !scalars
1439  integer :: ix,iy,iz,fft_idx
1440  real(dp) :: gdotr
1441 
1442 ! *************************************************************************
1443 
1444  if (ALL(gg==0)) then
1445    eigr(1:2*nfft:2)=one
1446    eigr(2:2*nfft:2)=zero
1447    RETURN
1448  end if
1449 
1450  fft_idx=1
1451  do iz=0,ngfft(3)-1
1452    do iy=0,ngfft(2)-1
1453      do ix=0,ngfft(1)-1
1454        gdotr= two_pi*( gg(1)*(ix/DBLE(ngfft(1))) &
1455 &                     +gg(2)*(iy/DBLE(ngfft(2))) &
1456 &                     +gg(3)*(iz/DBLE(ngfft(3))) )
1457        eigr(fft_idx  )=DCOS(gdotr)
1458        eigr(fft_idx+1)=DSIN(gdotr)
1459        fft_idx = fft_idx+2
1460      end do
1461    end do
1462  end do
1463 
1464 end subroutine calc_eigr

m_fft_mesh/check_rot_fft [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

   check_rot_fft

FUNCTION

  Return .TRUE. if the given grid in real space is compatible
  with the rotational part of the space group symmetries.

INPUTS

  nsym=Number of symmetry operations
  symrel(3,3,nsym)=Symmetry operations in real space.
  nr1,nr2,nr3=FFT divisions.

SOURCE

666 pure function check_rot_fft(nsym,symrel,nr1,nr2,nr3)
667 
668 
669 !This section has been created automatically by the script Abilint (TD).
670 !Do not modify the following lines by hand.
671 #undef ABI_FUNC
672 #define ABI_FUNC 'check_rot_fft'
673 !End of the abilint section
674 
675  implicit none
676 
677 !Arguments
678 !Scalar
679  integer,intent(in) :: nr1,nr2,nr3,nsym
680  logical :: check_rot_fft
681 !Arrays
682  integer,intent(in) :: symrel(3,3,nsym)
683 
684 !local variables
685  integer :: is
686 
687 !************************************************************************
688 
689  ! The grid is compatible with the symmetries (only rotational part) if
690  ! for each symmetry, each n_i and n_j ==> $n_i*R_{ij}/n_j$ is an integer
691  check_rot_fft=.TRUE.
692  do is=1,nsym
693    if ( MOD(symrel(2,1,is)*nr2, nr1) /=0 .or. &
694 &       MOD(symrel(3,1,is)*nr3, nr1) /=0 .or. &
695 &       MOD(symrel(1,2,is)*nr1, nr2) /=0 .or. &
696 &       MOD(symrel(3,2,is)*nr3, nr2) /=0 .or. &
697 &       MOD(symrel(1,3,is)*nr1, nr3) /=0 .or. &
698 &       MOD(symrel(2,3,is)*nr2, nr3) /=0      &
699 &     ) then
700      check_rot_fft=.FALSE.; EXIT
701    end if
702  end do
703 
704 end function check_rot_fft

m_fft_mesh/cigfft [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 cigfft

FUNCTION

 For each of the (2*nG0sh+1)**3 vectors G0 around the origin,
 calculate G-G0 and its FFT index number for all the NPWVEC vectors G.

INPUTS

 mG0(3)= For each reduced direction gives the max G0 component to account for umklapp processes.
 npwvec=Number of plane waves
 ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
 gvec(3,npwvec)=Reduced coordinates of G vectors.

OUTPUT

 igfft(npwvec,2*mG0(1)+1,2*mG0(2)+1,2*mG0(3)+1)=For each G, and each G0 vector,
  it gives the FFT grid index of the G-G0 vector.
 ierr=Number of G-G0 vectors falling outside the inout FFT box.

PARENTS

CHILDREN

      xcopy

SOURCE

 960 subroutine cigfft(mG0,npwvec,ngfft,gvec,igfft,ierr)
 961 
 962 
 963 !This section has been created automatically by the script Abilint (TD).
 964 !Do not modify the following lines by hand.
 965 #undef ABI_FUNC
 966 #define ABI_FUNC 'cigfft'
 967 !End of the abilint section
 968 
 969  implicit none
 970 
 971 !Arguments ------------------------------------
 972 !scalars
 973  integer,intent(in) :: npwvec
 974  integer,intent(out) :: ierr
 975 !arrays
 976  integer,intent(in) :: gvec(3,npwvec)
 977  integer,intent(in) :: mg0(3),ngfft(18)
 978  integer,intent(out) :: igfft(npwvec,2*mg0(1)+1,2*mg0(2)+1,2*mg0(3)+1)
 979 
 980 !Local variables ------------------------------
 981 !scalars
 982  integer :: gmg01,gmg02,gmg03,ig,ig01,ig02,ig03,n1,n2,n3
 983  character(len=500) :: msg
 984 !arrays
 985  integer :: gmg0(3)
 986 !************************************************************************
 987 
 988  DBG_ENTER("COLL")
 989 
 990  if (ANY(mg0<0)) then
 991    write(msg,'(a,3i4)')' Found negative value of mg0= ',mg0
 992    MSG_BUG(msg)
 993  end if
 994 
 995  n1=ngfft(1)
 996  n2=ngfft(2)
 997  n3=ngfft(3)
 998  ierr=0
 999 
1000  do ig=1,npwvec
1001    do ig01=-mg0(1),mG0(1)
1002      gmg0(1) = gvec(1,ig)-ig01
1003      do ig02=-mg0(2),mg0(2)
1004        gmg0(2) = gvec(2,ig)-ig02
1005        do ig03=-mg0(3),mg0(3)
1006          gmg0(3) = gvec(3,ig)-ig03
1007          ! === Calculate FFT index of G-G0 ===
1008          ! * Consider possible wrap around errors.
1009          gmg01=MODULO(gmg0(1),n1)
1010          gmg02=MODULO(gmg0(2),n2)
1011          gmg03=MODULO(gmg0(3),n3)
1012          igfft(ig,ig01+mg0(1)+1,ig02+mg0(2)+1,ig03+mg0(3)+1) = 1+gmg01+gmg02*n1+gmg03*n1*n2
1013          if ( ANY(gmg0>ngfft(1:3)/2) .or. ANY(gmg0<-(ngfft(1:3)-1)/2) ) then
1014            igfft(ig,ig01+mg0(1)+1,ig02+mg0(2)+1,ig03+mg0(3)+1) = 0
1015            ierr=ierr+1
1016          end if
1017        end do
1018      end do
1019    end do
1020  end do !ig
1021 
1022  if (ierr/=0) then
1023    write(msg,'(a,i4,3a)')&
1024 &    'Found ',ierr,' G-G0 vectors falling outside the FFT box. ',ch10,&
1025 &    'igfft will be set to zero for these particular G-G0 '
1026    MSG_WARNING(msg)
1027  end if
1028 
1029  DBG_EXIT("COLL")
1030 
1031 end subroutine cigfft

m_fft_mesh/fft_check_rotrans [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 fft_check_rotrans

FUNCTION

  Checks if the real space FFT mesh is compatible both with the rotational
  and the translational part of space group of the crystal.

INPUTS

  nsym=Number of symmetries.
  symrel(3,3,nsym)=Symmetries in real space in reduced coordinates.
  tnons(3,nsym)=Fractional translations.
  ngfft(18)=Information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft

OUTPUT

  err(3,nsym)=The max error for each symmetry. (given in terms of the FFT vectors)
  isok=.FALSE. if the FFT mesh does not fulfil all symmetry properties of the crystal.

PARENTS

      m_shirley

CHILDREN

SOURCE

734 function fft_check_rotrans(nsym,symrel,tnons,ngfft,err) result(isok)
735 
736 
737 !This section has been created automatically by the script Abilint (TD).
738 !Do not modify the following lines by hand.
739 #undef ABI_FUNC
740 #define ABI_FUNC 'fft_check_rotrans'
741 !End of the abilint section
742 
743  implicit none
744 
745 !Arguments ------------------------------------
746 !scalars
747  integer,intent(in) :: nsym
748  logical :: isok
749 !arrays
750  integer,intent(in) :: symrel(3,3,nsym)
751  integer,intent(in) :: ngfft(18)
752  real(dp),intent(in) :: tnons(3,nsym)
753  real(dp),intent(out) :: err(3,nsym)
754 
755 !Local variables-------------------------------
756 !scalars
757  integer :: isym,ix,iy,iz,jx,jy,jz,ngfft1,ngfft2,ngfft3
758  !character(len=500) :: msg
759 !arrays
760  integer :: Rm1(3,3,nsym),r1_FFT(3),red2fft(3,3)
761  real(dp) :: Rm1_FFT(3,3,nsym),fft2red(3,3),r2_FFT(3),tnons_FFT(3,nsym)
762 
763 ! *************************************************************************
764 
765  ! === Precalculate R^-1 and fractional translations in FFT coordinates ===
766  ngfft1=ngfft(1)
767  ngfft2=ngfft(2)
768  ngfft3=ngfft(3)
769 
770  red2fft=RESHAPE((/ngfft1,0,0,0,ngfft2,0,0,0,ngfft3/),(/3,3/))
771  fft2red=RESHAPE((/(one/ngfft1),zero,zero,zero,(one/ngfft2),zero,zero,zero,(one/ngfft3)/),(/3,3/))
772  !
773  ! === For a fully compatible mesh, each Rm1_FFT should be integer ===
774  do isym=1,nsym
775    call mati3inv(symrel(:,:,isym),Rm1(:,:,isym))
776    Rm1(:,:,isym)=TRANSPOSE(Rm1(:,:,isym))
777    Rm1_FFT(:,:,isym)=MATMUL(MATMUL(red2fft,Rm1(:,:,isym)),fft2red)
778    tnons_FFT(:,isym)=MATMUL(red2fft,tnons(:,isym))
779  end do
780 
781  err(:,:)=smallest_real
782  do iz=0,ngfft3-1
783    R1_FFT(3)=DBLE(iz)
784    do iy=0,ngfft2-1
785      R1_FFT(2)=DBLE(iy)
786      do ix=0,ngfft1-1
787        R1_FFT(1)=DBLE(ix)
788        do isym=1,nsym  ! Form R^-1 (r-\tau) in the FFT basis ===
789          R2_FFT(:)=MATMUL(Rm1_FFT(:,:,isym),R1_FFT(:)-tnons_FFT(:,isym))
790          jx=NINT(R2_FFT(1)); err(1,isym)=MAX(err(1,isym),ABS(R2_FFT(1)-jx)/ngfft1)
791          jy=NINT(R2_FFT(2)); err(2,isym)=MAX(err(2,isym),ABS(R2_FFT(2)-jy)/ngfft2)
792          jz=NINT(R2_FFT(3)); err(3,isym)=MAX(err(3,isym),ABS(R2_FFT(3)-jz)/ngfft3)
793        end do
794      end do
795    end do
796  end do
797 
798  isok=.TRUE.
799  do isym=1,nsym
800    if (ANY(err(:,isym)>tol6)) then
801      isok=.FALSE.
802      !write(msg,'(a,i3,a,3es14.6)')' symmetry ',isym,') not compatible with FFT grid, error ',err(:,isym)
803      !MSG_WARNING(msg)
804    end if
805  end do
806 
807 end function fft_check_rotrans

m_fft_mesh/g2ifft [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

  g2ifft

FUNCTION

 Returns the index of G in the FFT box from its reduced coordinates. 0 if not in the BOX.

INPUTS

  gg(3)=Reduced coordinated of the G vector.
  ngfft(18) = Info on the FFT box.

OUTPUT

  gidx=Index in the FFT box. 0 if G is outside the box.

SOURCE

1113 pure integer function g2ifft(gg,ngfft) result (gidx)
1114 
1115 
1116 !This section has been created automatically by the script Abilint (TD).
1117 !Do not modify the following lines by hand.
1118 #undef ABI_FUNC
1119 #define ABI_FUNC 'g2ifft'
1120 !End of the abilint section
1121 
1122  implicit none
1123 
1124 !This section has been created automatically by the script Abilint (TD).
1125 !Do not modify the following lines by hand.
1126 #undef ABI_FUNC
1127 #define ABI_FUNC 'g2ifft'
1128 !End of the abilint section
1129 
1130 !Arguments ------------------------------------
1131 !scalars
1132 !arrays
1133  integer,intent(in) :: gg(3),ngfft(3)
1134 
1135 !Local variables-------------------------------
1136 !scalars
1137  integer :: n1,n2,n3,ig1,ig2,ig3
1138 
1139 !************************************************************************
1140 
1141  ! Use the following indexing (N means ngfft of the adequate direction)
1142  ! 0 1 2 3 ... N/2    -(N-1)/2 ... -1    <= gg
1143  ! 1 2 3 4 ....N/2+1  N/2+2    ...  N    <= index
1144  !
1145  if ( ANY(gg>ngfft(1:3)/2) .or. ANY(gg<-(ngfft(1:3)-1)/2) ) then ! out of the box.
1146    gidx=0
1147  else
1148    n1=ngfft(1)
1149    n2=ngfft(2)
1150    n3=ngfft(3)
1151 
1152    ig1=MODULO(gg(1),n1)
1153    ig2=MODULO(gg(2),n2)
1154    ig3=MODULO(gg(3),n3)
1155 
1156    gidx = 1 + ig1 + n1*(ig2+ig3*n2)
1157  end if
1158 
1159 end function g2ifft

m_fft_mesh/get_gftt [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

  get_gftt

FUNCTION

  Returns the set of G-vectors in the FFT mesh and the maximal kinetic energy of k+G.

INPUTS

 ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
 kpt(3)=input k vector (reduced coordinates --in terms of reciprocal lattice primitive translations)
 gmet(3,3)=reciprocal space metric (bohr^-2)

OUTPUT

  gsq_max=Max value of (k+G)^2 for G in the FFT box
  gfft(3,nfft_tot) = The reduced components of the G in the FFT mesh (nfft_tot=PRODUCT(ngfft(1:3))

PARENTS

      bethe_salpeter,calc_sigc_me,cchi0,cchi0q0,cohsex_me,screening,sigma
      m_dvdb

CHILDREN

      xcopy

SOURCE

1189 pure subroutine get_gftt(ngfft,kpt,gmet,gsq_max,gfft)
1190 
1191 
1192 !This section has been created automatically by the script Abilint (TD).
1193 !Do not modify the following lines by hand.
1194 #undef ABI_FUNC
1195 #define ABI_FUNC 'get_gftt'
1196 !End of the abilint section
1197 
1198  implicit none
1199 
1200 !Arguments ------------------------------------
1201 !scalars
1202  real(dp),intent(out) :: gsq_max
1203 !arrays
1204  integer,intent(in) :: ngfft(18)
1205  integer,intent(out) :: gfft(3,ngfft(1)*ngfft(2)*ngfft(3))
1206  real(dp),intent(in) :: kpt(3),gmet(3,3)
1207 
1208 !Local variables-------------------------------
1209 !scalars
1210  integer :: ifft,g1,g2,g3,i1,i2,i3
1211  real(dp) :: dsq
1212 
1213 !************************************************************************
1214 
1215  ifft=0; gsq_max=smallest_real
1216  do i3=1,ngfft(3)    ; g3 = ig2gfft(i3,ngfft(3))
1217    do i2=1,ngfft(2)  ; g2 = ig2gfft(i2,ngfft(2))
1218      do i1=1,ngfft(1); g1 = ig2gfft(i1,ngfft(1))
1219        ifft=ifft+1
1220        gfft(1,ifft) = g1
1221        gfft(2,ifft) = g2
1222        gfft(3,ifft) = g3
1223        dsq=gmet(1,1)*(kpt(1)+dble(i1))**2&
1224 &        +gmet(2,2)*(kpt(2)+dble(i2))**2&
1225 &        +gmet(3,3)*(kpt(3)+dble(i3))**2&
1226 &        +2._dp*(gmet(1,2)*(kpt(1)+dble(i1))*(kpt(2)+dble(i2))&
1227 &        +gmet(2,3)*(kpt(2)+dble(i2))*(kpt(3)+dble(i3))&
1228 &        +gmet(3,1)*(kpt(3)+dble(i3))*(kpt(1)+dble(i1)))
1229        gsq_max = MAX(dsq,gsq_max)
1230      end do
1231    end do
1232  end do
1233 
1234 end subroutine get_gftt

m_fft_mesh/ig2gfft [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

  ig2gfft

FUNCTION

  Return the reduced component of a G-vector in the FFT mesh starting from is index.

INPUTS

  ig = The index >=1, <=ng
  ng = The number of FFT points along this direction.

OUTPUT

  gc = The reduced component

SOURCE

1052 elemental function ig2gfft(ig,ng) result (gc)
1053 
1054 
1055 !This section has been created automatically by the script Abilint (TD).
1056 !Do not modify the following lines by hand.
1057 #undef ABI_FUNC
1058 #define ABI_FUNC 'ig2gfft'
1059 !End of the abilint section
1060 
1061  implicit none
1062 
1063 !This section has been created automatically by the script Abilint (TD).
1064 !Do not modify the following lines by hand.
1065 #undef ABI_FUNC
1066 #define ABI_FUNC 'ig2gfft'
1067 !End of the abilint section
1068 
1069 !Arguments ------------------------------------
1070 !scalars
1071  integer,intent(in) :: ig,ng
1072  integer :: gc
1073 
1074 !************************************************************************
1075 
1076  ! Use the following indexing (N means ngfft of the adequate direction)
1077  ! 0 1 2 3 ... N/2    -(N-1)/2 ... -1    <= gc
1078  ! 1 2 3 4 ....N/2+1  N/2+2    ...  N    <= index ig
1079  !
1080  if (ig<=0 .or. ig > ng) then
1081    ! Wrong ig, returns huge. Parent code will likely crash with SIGSEV.
1082    gc = huge(1)
1083    return
1084  end if
1085 
1086  if ( ig  > ng/2 + 1) then
1087    gc = ig - ng -1
1088  else
1089    gc = ig -1
1090  end if
1091 
1092 end function ig2gfft

m_fft_mesh/phase [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 phase

FUNCTION

 Compute ph(ig)=$\exp(\pi\ i \ n/ngfft)$ for n=0,...,ngfft/2,-ngfft/2+1,...,-1
 while ig runs from 1 to ngfft.

INPUTS

  ngfft=number of points

OUTPUT

  ph(2*ngfft)=phase array (complex)

NOTES

 XG 990504 : changed the formulation, in order to preserve
 the invariance between n and -n, that was broken for n=ngfft/2 if ngfft even.
 Simply suppresses the corresponding sine.

PARENTS

      xcden,xcpot

CHILDREN

SOURCE

1716 subroutine phase(ngfft,ph)
1717 
1718 
1719 !This section has been created automatically by the script Abilint (TD).
1720 !Do not modify the following lines by hand.
1721 #undef ABI_FUNC
1722 #define ABI_FUNC 'phase'
1723 !End of the abilint section
1724 
1725  implicit none
1726 
1727 !Arguments ------------------------------------
1728 !scalars
1729  integer,intent(in) :: ngfft
1730 !arrays
1731  real(dp),intent(out) :: ph(2*ngfft)
1732 
1733 !Local variables-------------------------------
1734 !scalars
1735  integer :: id,ig,nn
1736  real(dp) :: arg,fac
1737 
1738 ! *************************************************************************
1739 
1740  id=ngfft/2+2
1741  fac=pi/dble(ngfft)
1742  do ig=1,ngfft
1743    nn=ig-1-(ig/id)*ngfft
1744    arg=fac*dble(nn)
1745    ph(2*ig-1)=cos(arg)
1746    ph(2*ig)  =sin(arg)
1747 
1748  end do
1749 !XG 990504 Here zero the corresponding sine
1750  if((ngfft/2)*2==ngfft) ph(2*(id-1))=zero
1751 
1752 end subroutine phase

m_fft_mesh/rotate_FFT_mesh [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 rotate_FFT_mesh

FUNCTION

  Find the FFT index of $ R{-1}(r-\tau) $ for each point in the FFT box.
  $R$ is a symmetry operation in real space, $\tau$ is the associated
  fractional translation.

INPUTS

  nsym=Number of symmetries.
  symrel(3,3,nsym)=Symmetries in real space in reduced coordinates.
  tnons(3,nsym)=Fractional translations.
  ngfft(18)=Information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft

OUTPUT

  irottb(ngfftot,nsym)=Indeces of $R^{-1}(r-\tau)$ in the FFT box.
  preserve=.FALSE. if the FFT mesh does not fulfil all symmetry properties of the crystal.

NOTES

  The evaluation of the rotated point $R^{-1}(r-\tau)$ is done using real arithmetic.
  As a consequence, if the FFT mesh does not fulfil the symmetry properties
  of the crystal, the array irottb will contain the index of the FFT point which
  is the closest one to $R^{-1}(r-\tau)$. This might lead to inaccuracies in the
  final results, in particular in the description of degenerate states.

PARENTS

      bethe_salpeter,calc_sigc_me,calc_sigx_me,cchi0q0_intraband
      classify_bands,cohsex_me,m_dvdb,m_wfd,prep_calc_ucrpa,screening

CHILDREN

      xcopy

SOURCE

847 subroutine rotate_fft_mesh(nsym,symrel,tnons,ngfft,irottb,preserve)
848 
849 
850 !This section has been created automatically by the script Abilint (TD).
851 !Do not modify the following lines by hand.
852 #undef ABI_FUNC
853 #define ABI_FUNC 'rotate_fft_mesh'
854 !End of the abilint section
855 
856  implicit none
857 
858 !Arguments ------------------------------------
859 !scalars
860  integer,intent(in) :: nsym
861  logical,intent(out) :: preserve
862 !arrays
863  integer,intent(in) :: symrel(3,3,nsym)
864  integer,intent(in) :: ngfft(18)
865  integer,intent(out) :: irottb(ngfft(1)*ngfft(2)*ngfft(3),nsym)
866  real(dp),intent(in) :: tnons(3,nsym)
867 
868 !Local variables-------------------------------
869 !scalars
870  integer :: ir1,isym,ix,iy,iz,jx,jy,jz,ngfft1,ngfft2,ngfft3
871  !character(len=500) :: msg
872 !arrays
873  integer :: Rm1(3,3,nsym),r1_FFT(3),red2fft(3,3)
874  real(dp) :: Rm1_FFT(3,3,nsym),err(3,nsym),fft2red(3,3),r2_FFT(3)
875  real(dp) :: tnons_FFT(3,nsym)
876 
877 ! *************************************************************************
878 
879  ! === Precalculate R^-1 and fractional translations in FFT coordinates ===
880  ngfft1=ngfft(1)
881  ngfft2=ngfft(2)
882  ngfft3=ngfft(3)
883 
884  red2fft=RESHAPE((/ngfft1,0,0,0,ngfft2,0,0,0,ngfft3/),(/3,3/))
885  fft2red=RESHAPE((/(one/ngfft1),zero,zero,zero,(one/ngfft2),zero,zero,zero,(one/ngfft3)/),(/3,3/))
886  !
887  ! === For a fully compatible mesh, each Rm1_FFT should be integer ===
888  do isym=1,nsym
889    call mati3inv(symrel(:,:,isym),Rm1(:,:,isym))
890    Rm1(:,:,isym)=TRANSPOSE(Rm1(:,:,isym))
891    Rm1_FFT(:,:,isym)=MATMUL(MATMUL(red2fft,Rm1(:,:,isym)),fft2red)
892    tnons_FFT(:,isym)=MATMUL(red2fft,tnons(:,isym))
893  end do
894 
895  err(:,:)=zero
896 
897 !$OMP PARALLEL DO PRIVATE(R1_FFT,ir1,R2_FFT,jx,jy,jz) reduction(MAX:err)
898  do iz=0,ngfft3-1
899    R1_FFT(3)=DBLE(iz)
900    do iy=0,ngfft2-1
901      R1_FFT(2)=DBLE(iy)
902      do ix=0,ngfft1-1
903        R1_FFT(1)=DBLE(ix)
904        ir1=1+ix+iy*ngfft1+iz*ngfft1*ngfft2
905        do isym=1,nsym
906          ! === Form R^-1 (r-\tau) in the FFT basis ===
907          R2_FFT(:)=MATMUL(Rm1_FFT(:,:,isym),R1_FFT(:)-tnons_FFT(:,isym))
908          jx=NINT(R2_FFT(1)); err(1,isym)=MAX(err(1,isym),ABS(R2_FFT(1)-jx)/ngfft1)
909          jy=NINT(R2_FFT(2)); err(2,isym)=MAX(err(2,isym),ABS(R2_FFT(2)-jy)/ngfft2)
910          jz=NINT(R2_FFT(3)); err(3,isym)=MAX(err(3,isym),ABS(R2_FFT(3)-jz)/ngfft3)
911          jx=MODULO(jx,ngfft1)
912          jy=MODULO(jy,ngfft2)
913          jz=MODULO(jz,ngfft3)
914          irottb(ir1,isym)=1+jx+jy*ngfft1+jz*ngfft1*ngfft2
915        end do
916      end do
917    end do
918  end do
919 
920  preserve=.TRUE.
921  do isym=1,nsym
922    if (ANY(err(:,isym)>tol6)) then
923      preserve=.FALSE.
924      !write(msg,'(a,i3,a,3es14.6)')' symmetry ',isym,') not compatible with FFT grid, error ',err(:,isym)
925      !MSG_WARNING(msg)
926    end if
927  end do
928 
929 end subroutine rotate_fft_mesh

m_fft_mesh/setmesh [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 setmesh

FUNCTION

 Calculate the size of the FFT grid for the GW calculation.

INPUTS

  gmet(3,3)=Reciprocal space metric.
  gvec(3,npwvec)=G-vectors in reduced coordinates.
  npwvec=Number of G vectors in the array gvec max(npwwfn,npwsigx)
  npwsigx=Size of the dielectric or self-energy matrix.
  npwwfn=Number of G-vectors in the wavefunctions.
  method=Integer flag for FFT grid (see below)
  mG0=Number of shells that must be added to take into account umklapp processes.
  Cryst<crystal_t>=Data type gathering information on unit cell and symmetries
    %nsym=Number of symmetry operations in the SG.
    %symrel(3,3,nsym)=Symmetry operations in real space.
    %tnons(3,nsym)=Fractional translations.
  enforce_sym=Flag to enforce a FFT which fulfils all symmetry operations, both the
   rotational part and fractional translations.
  [unit]=Output unit, defaults to std_out

OUTPUT

 ngfft(18)=contain all needed information about 3D FFT,
  see also ~abinit/doc/variables/vargs.htm#ngfft
 nfftot= ngfft(1)*ngfft(2)*ngfft(3)=Total number of points in the FFT grid.

NOTES

 Four methods are implemented for the calculation of the mesh:
  method=0 --> FFT mesh defined by the user, useful for debugging.
  method=1     Roughly takes the FFT box which encloses the larger of the two spheres of radius
               aliasing_factor*rwfn and rsigx, where rwfn and rsigx are the radius of the spheres
               with npwwfn and npwsigx planewaves respectively. The default aliasing_factor is 1.
  method=2 --> Calculates the optimal FFT grid which allows aliasing only outside the sphere of the
               npwsigx planewaves (finer than method=1 with aliasing_factor=1).
  method=3 --> Calculates the FFT grid needed to expand the density.
               (even finer than method=2, roughly corresponds to method=1 with aliasing_factor=2).

  See defs_fftdata for a list of allowed sizes of FFT.

PARENTS

      m_shirley,setup_bse,setup_screening,setup_sigma

CHILDREN

      xcopy

SOURCE

302 subroutine setmesh(gmet,gvec,ngfft,npwvec,npwsigx,npwwfn,nfftot,method,mG0,Cryst,enforce_sym,unit)
303 
304 
305 !This section has been created automatically by the script Abilint (TD).
306 !Do not modify the following lines by hand.
307 #undef ABI_FUNC
308 #define ABI_FUNC 'setmesh'
309 !End of the abilint section
310 
311  implicit none
312 
313 !Arguments ------------------------------------
314 !scalars
315  integer,intent(in) :: enforce_sym,method,npwsigx,npwvec,npwwfn
316  integer,intent(out) :: nfftot
317  integer,optional,intent(in) :: unit
318  type(crystal_t),target,intent(in) :: Cryst
319 !arrays
320  integer,intent(in) :: gvec(3,npwvec),mG0(3)
321  integer,intent(inout) :: ngfft(18)
322  real(dp),intent(in) :: gmet(3,3)
323 
324 !Local variables ------------------------------
325 !scalars
326  integer :: aliasing_factor,fftalg,fftalga,fftalgc,ig,ig1,ig1max,ig2,ig2max,ig3,ig3max,ii,idx,ierr
327  integer :: is,m1,m2,m3,mm1,mm2,mm3,n1,n2,n3,nsym,nt,ount
328  real(dp) :: ecuteff,ecutsigx,ecutwfn,g1,g2,g3,gsq,gsqmax,reff,rsigx,rwfn
329  logical :: fft_ok
330  character(len=500) :: msg
331 !arrays
332  integer :: fftnons(3),fftsym(3),mdum(3)
333  !integer,allocatable :: pfactors(:),powers(:)
334  integer,pointer :: symrel(:,:,:)
335  real(dp),pointer :: tnons(:,:)
336 
337 !************************************************************************
338 
339  DBG_ENTER("COLL")
340 
341  if (ANY(mg0<0)) then
342    write(msg,'(a,3(i0,1x))')' called with wrong value of mG0 = ',mG0
343    MSG_BUG(msg)
344  end if
345 
346  ount = std_out; if (present(unit)) ount = unit
347 
348  nsym   =  Cryst%nsym
349  symrel => Cryst%symrel
350  tnons  => Cryst%tnons
351  !
352  ! Calculate the limits of the sphere of npwwfn G-vectors in each direction.
353  m1=MAXVAL(ABS(gvec(1,1:npwwfn)))
354  m2=MAXVAL(ABS(gvec(2,1:npwwfn)))
355  m3=MAXVAL(ABS(gvec(3,1:npwwfn)))
356  !
357  ! Calculate the limits of the sphere of npsigx G-vectors in each direction.
358  ! Ensure that G+G0 will fit into the FFT grid, where G is any of the npwsigx/npweps vectors
359  ! and G0 is (i,j,k) [-nG0shell<i,j,k<nG0shell]. This is required when npwsigx>npwwfn since
360  ! we have to take into account umklapp G0 vectors to evaluate the oscillator matrix elements
361  ! (see rho_tw_g) or to symmetrize these quantities (see also cigfft).
362  mm1=MAXVAL(ABS(gvec(1,1:npwsigx)))
363  mm2=MAXVAL(ABS(gvec(2,1:npwsigx)))
364  mm3=MAXVAL(ABS(gvec(3,1:npwsigx)))
365 
366  mm1=mm1+mG0(1)
367  mm2=mm2+mG0(2)
368  mm3=mm3+mG0(3)
369 
370  ! To avoid possible wrap-around errors in cigfft, it is safe to start
371  ! with odd divisions so that the FFT box is centered on Gamma
372  ! This holds only if npwsigx > npwwfn.
373  !if (iseven(mm1)) mm1=mm1+1
374  !if (iseven(mm2)) mm2=mm2+1
375  !if (iseven(mm3)) mm3=mm3+1
376 
377  write(msg,'(2(2a,i8,a,3i6),2a,3i3)')ch10,&
378 &  ' setmesh: npwwfn        = ',npwwfn, '; Max (m1,m2,m3)   = ',m1,m2,m3,ch10,&
379 &  '          npweps/npwsigx= ',npwsigx,'; Max (mm1,mm2,mm3)= ',mm1,mm2,mm3,ch10,&
380 &  '          mG0 added     = ',mG0(:)
381  call wrtout(ount,msg,'COLL')
382  !
383  ! === Different FFT grids according to method ==
384  select case (method)
385 
386  case (0)
387    ! * FFT mesh defined by user, useful for testing.
388    n1=ngfft(1)
389    n2=ngfft(2)
390    n3=ngfft(3)
391    write(msg,'(3(a,i3))')' Mesh size enforced by user = ',n1,'x',n2,'x',n3
392    MSG_COMMENT(msg)
393 
394    ngfft(1)=n1
395    ngfft(2)=n2
396    ngfft(3)=n3
397    ngfft(4)=2*(ngfft(1)/2)+1
398    ngfft(5)=2*(ngfft(2)/2)+1
399    ngfft(6)=   ngfft(3)
400    !ngfft(4:6)=ngfft(1:3)
401    nfftot=n1*n2*n3
402    RETURN
403 
404  case (1)
405    aliasing_factor=1
406    write(msg,'(2a,i3)')ch10,' using method 1 with aliasing_factor = ',aliasing_factor
407    call wrtout(ount,msg,'COLL')
408    m1=m1*aliasing_factor
409    m2=m2*aliasing_factor
410    m3=m3*aliasing_factor
411 
412  case (2,3)
413 
414    ecutwfn=-one  ! Calculate the radius of the sphere of npwwfn G-vectors.
415    do ig=1,npwwfn
416      g1=REAL(gvec(1,ig))
417      g2=REAL(gvec(2,ig))
418      g3=REAL(gvec(3,ig))
419      gsq=       gmet(1,1)*g1**2+gmet(2,2)*g2**2+gmet(3,3)*g3**2+ &
420 &          two*(gmet(1,2)*g1*g2+gmet(1,3)*g1*g3+gmet(2,3)*g2*g3)
421      ecutwfn=MAX(ecutwfn,gsq)
422    end do
423    rwfn=SQRT(ecutwfn); ecutwfn=two*ecutwfn*pi**2
424 
425    ! * Calculate the radius of the sphere of (npwsigx|npweps) G-vectors.
426    ecutsigx=-one
427    do ig=1,npwsigx
428      g1=REAL(gvec(1,ig))
429      g2=REAL(gvec(2,ig))
430      g3=REAL(gvec(3,ig))
431      gsq=      gmet(1,1)*g1**2+gmet(2,2)*g2**2+gmet(3,3)*g3**2+ &
432 &         two*(gmet(1,2)*g1*g2+gmet(1,3)*g1*g3+gmet(2,3)*g2*g3)
433      ecutsigx=MAX(ecutsigx,gsq)
434    end do
435    rsigx=SQRT(ecutsigx); ecutsigx=two*ecutsigx*pi**2
436 
437    write(msg,'(a,f7.3,3a,f7.3,a)')&
438 &    ' calculated ecutwfn          = ',ecutwfn, ' [Ha] ',ch10,&
439 &    ' calculated ecutsigx/ecuteps = ',ecutsigx,' [Ha]'
440    call wrtout(ount,msg,'COLL')
441    !
442    ! In the calculation of the GW self-energy or of the RPA dielectric matrix,
443    ! we have products $ \rho_{12}(r)=u_1*(r) u_2(r) $ of wavefunctions whose Fourier
444    ! coefficients lie in the sphere of radius rwfn. Such products will have non
445    ! vanishing Fourier coefficients in the whole sphere of radius 2*rwfn since:
446    !  $ rho_{12}(G) = \sum_T u_1*(T) u_2(T+G) $.
447    ! However, we only need the Fourier coefficients of $rho_{12}$ that lie in the sphere
448    ! of radius rsigx. We can thus allow aliasing outside that sphere, so that the FFT box
449    ! will only enclose a sphere of radius reff given by:
450 
451    reff=rsigx+rwfn
452    if (method==3) reff=two*rwfn ! Yields back the GS FFT grid if full wavefunctions are considered.
453    ecuteff=two*(pi*reff)**2
454    gsqmax=reff**2
455 
456    write(msg,'(a,i2,a,f7.3,a)')' using method = ',method,' with ecuteff = ',ecuteff,' [Ha]'
457    call wrtout(ount,msg,'COLL')
458    !
459    ! === Search the limits of the reff sphere in each direction ===
460    !ig1max=2*m1+1
461    !ig2max=2*m2+1
462    !ig3max=2*m3+1
463    if (method==2) then
464      ig1max=mm1+m1+1
465      ig2max=mm2+m2+1
466      ig3max=mm3+m3+1
467    else if (method==3) then
468      ig1max=MAX(2*m1+1,2*mm1+1,mm1+m1+1)
469      ig2max=MAX(2*m2+1,2*mm2+1,mm2+m2+1)
470      ig3max=MAX(2*m3+1,2*mm3+1,mm3+m3+1)
471    else
472      write(msg,'(a,i0)')"Wrong method : ",method
473      MSG_BUG(msg)
474    end if
475 
476    m1=-1; m2=-1; m3=-1
477    do ig1=0,ig1max
478      do ig2=0,ig2max
479        do ig3=0,ig3max
480          g1=REAL(ig1)
481          g2=REAL(ig2)
482          g3=REAL(ig3)
483          gsq=     gmet(1,1)*g1**2+gmet(2,2)*g2**2+gmet(3,3)*g3**2+ &
484 &            two*(gmet(1,2)*g1*g2+gmet(1,3)*g1*g3+gmet(2,3)*g2*g3)
485          if (gsq>gsqmax+tol6) CYCLE ! tol6 to improve portability
486          m1=MAX(m1,ig1)
487          m2=MAX(m2,ig2)
488          m3=MAX(m3,ig3)
489        end do
490      end do
491    end do
492 
493  case default
494    write(msg,'(a,i0)')' Method > 3 or < 0 not allowed in setmesh while method= ',method
495    MSG_BUG(msg)
496  end select
497  !
498  ! * Warning if low npwwfn.
499  if (m1<mm1 .or. m2<mm2 .or. m3<mm3) then
500    write(msg,'(5a)')&
501 &    'Note that npwwfn is small with respect to npweps or with respect to npwsigx. ',ch10,&
502 &    'Such a small npwwfn is a waste: ',ch10,&
503 &    'You could raise npwwfn without loss in cpu time. '
504    MSG_COMMENT(msg)
505  end if
506  !
507  ! Keep the largest of the m/mm and and find the FFT grid which is compatible
508  ! with the library and, if required, with the symmetry operations.
509  m1=MAX(m1,mm1)
510  m2=MAX(m2,mm2)
511  m3=MAX(m3,mm3)
512 
513  if (enforce_sym==0) then
514    ! === Determine the best size for the FFT grid *without* considering the symm ops ===
515    ! * Ideally n=2*m+1 but this  could not be allowed by the FFT library.
516    call size_goed_fft(m1,n1,ierr)
517    call size_goed_fft(m2,n2,ierr)
518    call size_goed_fft(m3,n3,ierr)
519    ABI_CHECK(ierr==0,"size_goed_fft failed")
520    nfftot=n1*n2*n3
521 
522    ! * Check if the FFT is compatible, write ONLY a warning if it breaks the symmetry
523    fftnons(1)=n1
524    fftnons(2)=n2
525    fftnons(3)=n3
526    fft_ok=.TRUE.
527    rd: do ii=1,3
528      do is=1,nsym
529        nt=denominator(tnons(ii,is), ierr)
530        if (((fftnons(ii)/nt)*nt) /= fftnons(ii)) then
531          fft_ok=.FALSE.; EXIT rd
532        end if
533      end do
534    end do rd
535    !
536    ! * Warn if not compatibile with tnons or rotational part.
537    if (.not.fft_ok) then
538     MSG_WARNING('FFT mesh is not compatible with non-symmorphic translations')
539    end if
540    if (.not.(check_rot_fft(nsym,symrel,n1,n2,n3))) then
541      MSG_WARNING('FFT mesh is not compatible with rotations')
542    end if
543 
544  else
545    ! === Determine the best size for the FFT grid considering symm ops ===
546    ! * Ideally n=2*m+1 but this could not be allowed by the FFT library (at present only Goedecker)
547    call wrtout(ount,' Finding a FFT mesh compatible with all the symmetries','COLL')
548 
549    ! 1) Find a FFT mesh compatible with the non-symmorphic operations
550    fftnons(:)=1
551    do ii=1,3
552      fftnons(ii)=1
553      do is=1,nsym
554        nt=denominator(tnons(ii,is), ierr)
555        if (((fftnons(ii)/nt)*nt)/=fftnons(ii)) fftnons(ii)=mincm(fftnons(ii),nt)
556      end do
557    end do
558    write(msg,'(a,3(i0,1x))')' setmesh: divisor mesh',fftnons(:)
559    call wrtout(ount,msg,'COLL')
560    !
561    ! 2) Check if also rotations preserve the grid.
562    ! * Use previous m values as Initial guess.
563    call size_goed_fft(m1,fftsym(1),ierr)
564    call size_goed_fft(m2,fftsym(2),ierr)
565    call size_goed_fft(m3,fftsym(3),ierr)
566    ABI_CHECK(ierr==0,"size_goed_fft failed")
567    mdum(1)=m1
568    mdum(2)=m2
569    mdum(3)=m3
570 
571    idx=0
572    do ! If a FFT division gets too large the code stops in size_goed_fft.
573      if ( check_rot_fft(nsym,symrel,fftsym(1),fftsym(2),fftsym(3)) .and. &
574 &         (MOD(fftsym(1),fftnons(1))==0) .and.                           &
575 &         (MOD(fftsym(2),fftnons(2))==0) .and.                           &
576 &         (MOD(fftsym(3),fftnons(3))==0)                                 &
577 &       ) EXIT
578      ii=MOD(idx,3)+1
579      mdum(ii)=mdum(ii)+1
580      call size_goed_fft(mdum(ii),fftsym(ii),ierr)
581      ABI_CHECK(ierr==0,"size_goed_fft failed")
582      idx=idx+1
583    end do
584    !
585    ! Got a good FFT grid, Calculate the number of FFT grid points
586    n1=fftsym(1)
587    n2=fftsym(2)
588    n3=fftsym(3); nfftot=n1*n2*n3
589 
590    if (.not.( check_rot_fft(nsym,symrel,n1,n2,n3)) &
591 &       .or.( MOD(fftsym(1),fftnons(1))/=0) .and.  &
592 &           ( MOD(fftsym(2),fftnons(2))/=0) .and.  &
593 &           ( MOD(fftsym(3),fftnons(3))/=0)        &
594 &     ) then
595      MSG_BUG('Not able to generate a symmetric FFT')
596    end if
597  end if ! enforce_sym
598 
599  write(msg,'(3(a,i3),2a,i8,a)')&
600 &  ' setmesh: FFT mesh size selected  = ',n1,'x',n2,'x',n3,ch10,&
601 &  '          total number of points  = ',nfftot,ch10
602  call wrtout(ount,msg,'COLL')
603  if (ount /= dev_null) then
604    call wrtout(ab_out,msg,'COLL')
605  end if
606 
607  ngfft(1)=n1
608  ngfft(2)=n2
609  ngfft(3)=n3
610  ngfft(4)=2*(ngfft(1)/2)+1
611  ngfft(5)=2*(ngfft(2)/2)+1
612  ngfft(6)=   ngfft(3)
613  !ngfft(4:6) = ngfft(1:3)
614  !
615  ! === Check the value of fftalg i.e ngfft(7) ===
616  ! * Presently only Goedecker"s library or FFTW3 are allowed, see size_goed_fft.F90
617  fftalg=ngfft(7); fftalga=fftalg/100; fftalgc=MOD(fftalg,10)
618 
619  if ( ALL(fftalga /= (/FFT_SG,FFT_FFTW3, FFT_DFTI/)) ) then
620    write(msg,'(6a)')ch10,&
621 &    "Only Goedecker's routines with fftalg=1xx or FFTW3/DFTI routines are allowed in GW calculations. ",ch10,&
622 &    "Action : check the value of fftalg in your input file, ",ch10,&
623 &    "or modify setmesh.F90 to make sure the FFT mesh is compatible with the FFT library. "
624    MSG_ERROR(msg)
625  end if
626 
627 ! TODO Had to change setmesh to avoid bad values for FFTW3
628 ! if (fftalga==3) then ! check whether mesh is optimal for FFTW3
629 !   ABI_MALLOC(pfactors,(5))
630 !   ABI_MALLOC(powers,(6))
631 !   pfactors = (/2, 3, 5, 7, 11/)
632 !   do ii=1,3
633 !     call pfactorize(ngfft(ii),5,pfactors,powers)
634 !     if (powers(6)/=1 .or. powers(4)/=0 .or. powers(5)/=0) then
635 !       write(msg,'(a,i0,a)')&
636 !&        "ngfft(ii) ",ngfft(ii)," contains powers of 7-11 or greater; FFTW3 is not optimal "
637 !       MSG_WARNING(msg)
638 !     end if
639 !   end do
640 !   ABI_FREE(pfactors)
641 !   ABI_FREE(powers)
642 ! end if
643 
644  DBG_EXIT("COLL")
645 
646 end subroutine setmesh

m_fft_mesh/times_eigr [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 times_eigr

FUNCTION

  Multiply an array on the real-space mesh by e^{iG0.r} where G0 is a reciprocal lattice vector.

INPUTS

  gg(3)=G vector in reduced coordinates.
  ngfft(18)=information about 3D FFT,
  nfft=Number of points in the FFT mesh.
  ndat=Number of arrays

SIDE EFFECTS

  ur(2,nfft,ndat)= contains u(r) in input. output: u(r) e^{ig.r} on the real-space FFT mesh.

PARENTS

CHILDREN

SOURCE

1561 pure subroutine times_eigr(gg,ngfft,nfft,ndat,ur)
1562 
1563 
1564 !This section has been created automatically by the script Abilint (TD).
1565 !Do not modify the following lines by hand.
1566 #undef ABI_FUNC
1567 #define ABI_FUNC 'times_eigr'
1568 !End of the abilint section
1569 
1570  implicit none
1571 
1572 !Arguments ------------------------------------
1573 !scalars
1574  integer,intent(in) :: nfft,ndat
1575 !arrays
1576  integer,intent(in) :: gg(3)
1577  integer,intent(in) :: ngfft(18)
1578  real(dp),intent(inout) :: ur(2,nfft,ndat)
1579 
1580 !Local variables-------------------------------
1581 !scalars
1582  integer :: ix,iy,iz,ifft,idat
1583  real(dp) :: gr
1584 !arrays
1585  real(dp) :: ph(2),val(2)
1586 
1587 ! *************************************************************************
1588 
1589  if (all(gg==0)) return
1590 
1591  do idat=1,ndat
1592    ifft = 0
1593    do iz=0,ngfft(3)-1
1594      do iy=0,ngfft(2)-1
1595        do ix=0,ngfft(1)-1
1596          ifft = ifft + 1
1597          gr = two_pi*(gg(1)*(ix/dble(ngfft(1))) &
1598                      +gg(2)*(iy/dble(ngfft(2))) &
1599                      +gg(3)*(iz/dble(ngfft(3))) )
1600          ph(1) = cos(gr); ph(2) = sin(gr)
1601          val(1) = ur(1,ifft,idat); val(2) = ur(2,ifft,idat)
1602 
1603          ur(1,ifft,idat) = ph(1) * val(1) - ph(2) * val(2)
1604          ur(2,ifft,idat) = ph(1) * val(2) + ph(2) * val(1)
1605        end do
1606      end do
1607    end do
1608  end do ! idat
1609 
1610 end subroutine times_eigr

m_fft_mesh/times_eikr [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

 times_eikr

FUNCTION

  Multiply an array on the real-space mesh by e^{ik.r} where k is a real(dp) vector in reduced coordinates

INPUTS

  kk(3)=k-vector in reduced coordinates.
  ngfft(18)=information about 3D FFT,
  nfft=Number of points in the FFT mesh.
  ndat=Number of arrays to transform

SIDE EFFECTS

  ur(2,nfft)= contains u(r) in input. output: u(r) e^{ig.r} on the real-space FFT mesh.

PARENTS

  m_dvdb

CHILDREN

SOURCE

1638 pure subroutine times_eikr(kk,ngfft,nfft,ndat,ur)
1639 
1640 
1641 !This section has been created automatically by the script Abilint (TD).
1642 !Do not modify the following lines by hand.
1643 #undef ABI_FUNC
1644 #define ABI_FUNC 'times_eikr'
1645 !End of the abilint section
1646 
1647  implicit none
1648 
1649 !Arguments ------------------------------------
1650 !scalars
1651  integer,intent(in) :: nfft,ndat
1652 !arrays
1653  real(dp),intent(in) :: kk(3)
1654  integer,intent(in) :: ngfft(18)
1655  real(dp),intent(inout) :: ur(2,nfft,ndat)
1656 
1657 !Local variables-------------------------------
1658 !scalars
1659  integer :: ix,iy,iz,ifft,idat
1660  real(dp) :: kr
1661 !arrays
1662  real(dp) :: ph(2),val(2)
1663 
1664 ! *************************************************************************
1665 
1666  if (all(abs(kk) < tol12)) return
1667 
1668  do idat=1,ndat
1669    ifft = 0
1670    do iz=0,ngfft(3)-1
1671      do iy=0,ngfft(2)-1
1672        do ix=0,ngfft(1)-1
1673          ifft = ifft + 1
1674          kr = two_pi*(kk(1)*(ix/dble(ngfft(1))) &
1675                      +kk(2)*(iy/dble(ngfft(2))) &
1676                      +kk(3)*(iz/dble(ngfft(3))) )
1677          ph(1) = cos(kr); ph(2) = sin(kr)
1678          val(1) = ur(1,ifft,idat); val(2) = ur(2,ifft,idat)
1679 
1680          ur(1,ifft,idat) = ph(1) * val(1) - ph(2) * val(2)
1681          ur(2,ifft,idat) = ph(1) * val(2) + ph(2) * val(1)
1682        end do
1683      end do
1684    end do
1685  end do
1686 
1687 end subroutine times_eikr

m_fft_mesh/zpad_free [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

  zpad_free

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

      xcopy

SOURCE

222 subroutine zpad_free(zpad)
223 
224 
225 !This section has been created automatically by the script Abilint (TD).
226 !Do not modify the following lines by hand.
227 #undef ABI_FUNC
228 #define ABI_FUNC 'zpad_free'
229 !End of the abilint section
230 
231  implicit none
232 
233 !Arguments ------------------------------------
234 !scalars
235  type(zpad_t),intent(inout) :: zpad
236 
237 ! *************************************************************************
238 
239  if (allocated(zpad%zplane)) then
240    ABI_FREE(zpad%zplane)
241  end if
242 
243  if (allocated(zpad%linex2ifft_yz)) then
244    ABI_FREE(zpad%linex2ifft_yz)
245  end if
246 
247 end subroutine zpad_free

m_fft_mesh/zpad_init [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

  zpad_init

FUNCTION

  Creation method

INPUTS

   mgfft=MAX(nx,ny,nz), only used to dimension gbound
   gbound(2*mgfft+8,2)= The boundaries of the basis sphere of G vectors at a given k-point.
     See sphereboundary for more info.

OUTPUT

  zpad<type(zpad_t)>

PARENTS

CHILDREN

      xcopy

SOURCE

129 subroutine zpad_init(zpad,nx,ny,nz,ldx,ldy,ldz,mgfft,gbound)
130 
131 
132 !This section has been created automatically by the script Abilint (TD).
133 !Do not modify the following lines by hand.
134 #undef ABI_FUNC
135 #define ABI_FUNC 'zpad_init'
136 !End of the abilint section
137 
138  implicit none
139 
140 !Arguments ------------------------------------
141 !scalars
142  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,mgfft
143  type(zpad_t),intent(out) :: zpad
144 !arrays
145  integer,intent(in) :: gbound(2*mgfft+8,2)
146 
147 !Local variables-------------------------------
148 !scalars
149  integer :: jj,g3_max,g3_min,gg3,ifft_g3,igb,g2min,g2max,nlinex
150 
151 ! *************************************************************************
152 
153  g3_min = gbound(3,2)
154  g3_max = gbound(4,2)
155 
156  zpad%n_zplanes = g3_max - g3_min + 1
157 
158  ABI_MALLOC(zpad%zplane,      (2,nz))
159  ABI_MALLOC(zpad%linex2ifft_yz, (2,nx*ny*nz))
160  !
161  ! Loop over the z-planes intersecting the G-sphere.
162  nlinex = 0
163  do gg3=1,zpad%n_zplanes
164    !
165    if (gg3<=g3_max+1) then
166      ifft_g3 = gg3
167    else
168      ifft_g3 = gg3 + nz - zpad%n_zplanes ! Wrap around for negative gg3.
169    end if
170    !
171    ! Select the set of y for this z-plane.
172    igb=2*gg3+3
173    g2min = gbound(igb  ,2)
174    g2max = gbound(igb+1,2)
175 
176    zpad%zplane(1,gg3) = ifft_g3
177    zpad%zplane(2,gg3) = igb
178 
179    !(1:g2max+1,ifft_g3)     ! Positive g_y.
180    !(g2min+ny+1:ny,ifft_g3) ! Negative g_y.
181 
182    do jj=1,g2max+1
183      nlinex = nlinex + 1
184      zpad%linex2ifft_yz(1,nlinex) = jj
185      zpad%linex2ifft_yz(2,nlinex) = ifft_g3
186    end do
187 
188    do jj=g2min+ny+1,ny
189      nlinex = nlinex + 1
190      zpad%linex2ifft_yz(1,nlinex) = jj
191      zpad%linex2ifft_yz(2,nlinex) = ifft_g3
192    end do
193  end do
194 
195  zpad%nlinex = nlinex
196 
197  RETURN
198  ABI_UNUSED((/ldx,ldy,ldz/))
199 
200 end subroutine zpad_init

m_fft_mesh/zpad_t [ Types ]

[ Top ] [ m_fft_mesh ] [ Types ]

NAME

  zpad_t

FUNCTION

   Store tables used for zero-padded FFTs.

SOURCE

79  type,public :: zpad_t
80 
81    integer :: nlinex
82    ! Total number of 1D transforms.
83 
84    integer :: n_zplanes
85    ! Number of z-planes intersecting the sphere.
86 
87    integer,allocatable :: zplane(:,:)
88    ! zplane(3,n_zplanes)
89    ! zplane(1,zpl) : mapping z-plane index -> FFT index_z
90    ! zplane(2,zpl) : mapping z-plane index -> igb index in array gbound
91 
92    integer,allocatable :: linex2ifft_yz(:,:)
93    ! linex2ifft_yz(2,nlinex)
94    ! mapping 1D-FFT -> (FFT_index_y, FFT index_z)
95 
96  end type zpad_t
97 
98  public :: zpad_init
99  public :: zpad_free