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-2024 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 .

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 MODULE m_fft_mesh
25 
26  use defs_basis
27  use m_errors
28  use m_abicore
29  use m_hide_blas
30  use, intrinsic :: iso_c_binding
31 
32  use defs_fftdata,     only : size_goed_fft
33  use m_fstrings,       only : sjoin, itoa, ltoa
34  use m_numeric_tools,  only : denominator, mincm, iseven, pfactorize
35  use m_symtk,          only : mati3inv
36  use m_geometry,       only : xred2xcart
37  use m_crystal,        only : crystal_t
38 
39  implicit none
40 
41  private
42 
43  public :: setmesh             ! Perform the setup of the FFT mesh for the GW oscillator strengths.
44  public :: check_rot_fft       ! Test whether the mesh is compatible with the rotational part of the space group.
45  public :: fft_check_rotrans   ! Test whether the mesh is compatible with the symmetries of the space group.
46  public :: rotate_fft_mesh     ! Calculate the FFT index of the rotated mesh.
47  public :: denpot_project      ! Compute n(r) + n( $R^{-1}(r-\tau)$ in) / 2
48                                ! Mainly used with R = inversion to select the even/odd part under inversion
49  public :: cigfft              ! Calculate the FFT index of G-G0.
50  public :: ig2gfft             ! Returns the component of a G in the FFT Box from its sequential index.
51  public :: g2ifft              ! Returns the index of the G in the FFT box from its reduced coordinates.
52  public :: get_gftt            ! Calculate the G"s in the FFT box from ngfft
53  public :: calc_ceigr          ! e^{iG.r} on the FFT mesh (complex valued).
54  public :: calc_eigr           ! e^{iG.r} on the FFT mesh (version for real array with RE,IM).
55  public :: calc_ceikr          ! e^{ik.r} on the FFT mesh (complex valued).
56  public :: times_eigr          ! Multiply an array on the real-space mesh by e^{iG0.r}
57  public :: times_eikr          ! Multiply an array on the real-space mesh by e^{ik.r}
58  public :: ctimes_eikr         ! Version for complex array
59  public :: phase               ! Compute ph(ig)=$\exp(\pi\ i \ n/ngfft)$ for n=0,...,ngfft/2,-ngfft/2+1,...,-1
60  public :: mkgrid_fft          ! 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
66 
67  interface calc_ceikr
68    module procedure calc_ceikr_spc
69    module procedure calc_ceikr_dpc
70  end interface calc_ceikr
71 
72 
73  !interface times_eikr
74  !  module procedure times_eikr_dp
75  !  module procedure ctimes_eikr_dpc
76  !end interface times_eikr

ABINIT/mkgrid_fft [ Functions ]

[ Top ] [ Functions ]

NAME

  mkgrid_fft

FUNCTION

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

INPUTS

OUTPUT

SOURCE

1686 subroutine mkgrid_fft(ffti3_local,fftn3_distrib,gridcart,nfft,ngfft,rprimd)
1687 
1688 !Arguments ------------------------------------
1689  integer, intent(in) :: nfft
1690  integer,intent(in) :: ngfft(18)
1691  integer, dimension(*), intent(in) :: ffti3_local,fftn3_distrib
1692  real(dp), dimension(3,nfft), intent(out) :: gridcart
1693  real(dp),intent(in) :: rprimd(3,3)
1694 
1695 !Local variables-------------------------------
1696  integer :: ind,i1,i2,i3,i3loc,me,nproc
1697  integer :: n1,n2,n3
1698  real(dp), dimension(3) :: coord
1699  real(dp), dimension(3,nfft) :: gridred
1700 
1701 ! *************************************************************************
1702 
1703  n1    = ngfft(1)
1704  n2    = ngfft(2)
1705  n3    = ngfft(3)
1706  nproc = ngfft(10)
1707  me    = ngfft(11)
1708 
1709  do i3 = 1, n3, 1
1710    if(fftn3_distrib(i3) == me) then !MPI
1711      i3loc=ffti3_local(i3)
1712      coord(3) = real(i3 - 1, dp) / real(n3, dp)
1713      do i2 = 1, n2, 1
1714        coord(2) = real(i2 - 1, dp) / real(n2, dp)
1715        do i1 = 1, n1, 1
1716          ind=i1+(i2-1)*n1+(i3loc-1)*n1*n2
1717          coord(1) = real(i1 - 1, dp) / real(n1, dp)
1718          gridred(:, ind) = coord(:)
1719        end do
1720      end do
1721    end if
1722  end do
1723  call xred2xcart(nfft, rprimd, gridcart, gridred)
1724 
1725 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.

SOURCE

1254 subroutine calc_ceigr_dpc(gg, nfft, nspinor, ngfft, ceigr)
1255 
1256 !Arguments ------------------------------------
1257 !scalars
1258  integer,intent(in) :: nfft,nspinor
1259 !arrays
1260  integer,intent(in) :: gg(3)
1261  integer,intent(in) :: ngfft(18)
1262  complex(dpc),intent(out) :: ceigr(nfft*nspinor)
1263 
1264 !Local variables-------------------------------
1265 !scalars
1266  integer :: ix,iy,iz,fft_idx,base,isp
1267  real(dp) :: gdotr
1268 
1269 ! *************************************************************************
1270 
1271  if (ALL(gg==0)) then
1272    ceigr=cone; RETURN
1273  end if
1274 
1275  fft_idx=0
1276  do iz=0,ngfft(3)-1
1277    do iy=0,ngfft(2)-1
1278      do ix=0,ngfft(1)-1
1279        gdotr= two_pi*( gg(1)*(ix/DBLE(ngfft(1))) &
1280                       +gg(2)*(iy/DBLE(ngfft(2))) &
1281                       +gg(3)*(iz/DBLE(ngfft(3))) )
1282        fft_idx = fft_idx+1
1283        ceigr(fft_idx)=DCMPLX(DCOS(gdotr),DSIN(gdotr))
1284      end do
1285    end do
1286  end do
1287 
1288  if (nspinor > 1) then
1289    do isp=2,nspinor
1290      base = 1 + (isp-1)*nfft
1291      call xcopy(nfft,ceigr,1,ceigr(base:),1)
1292    end do
1293  end if
1294 
1295 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.

SOURCE

1189 subroutine calc_ceigr_spc(gg, nfft, nspinor, ngfft, ceigr)
1190 
1191 !Arguments ------------------------------------
1192 !scalars
1193  integer,intent(in) :: nfft,nspinor
1194 !arrays
1195  integer,intent(in) :: gg(3)
1196  integer,intent(in) :: ngfft(18)
1197  complex(spc),intent(out) :: ceigr(nfft*nspinor)
1198 
1199 !Local variables-------------------------------
1200 !scalars
1201  integer :: ix,iy,iz,fft_idx,base,isp
1202  real(dp) :: gdotr
1203 
1204 ! *************************************************************************
1205 
1206  if (ALL(gg==0)) then
1207    ceigr=(1._sp,0._sp)
1208    RETURN
1209  end if
1210 
1211  fft_idx=0
1212  do iz=0,ngfft(3)-1
1213    do iy=0,ngfft(2)-1
1214      do ix=0,ngfft(1)-1
1215        gdotr= two_pi*( gg(1)*(ix/DBLE(ngfft(1))) &
1216                       +gg(2)*(iy/DBLE(ngfft(2))) &
1217                       +gg(3)*(iz/DBLE(ngfft(3))) )
1218        fft_idx = fft_idx+1
1219        ceigr(fft_idx)=CMPLX(DCOS(gdotr),DSIN(gdotr), KIND=spc)
1220      end do
1221    end do
1222  end do
1223 
1224  if (nspinor > 1) then
1225    do isp=2,nspinor
1226      base = 1 + (isp-1)*nfft
1227      call xcopy(nfft,ceigr,1,ceigr(base:),1)
1228    end do
1229  end if
1230 
1231 end subroutine calc_ceigr_spc

m_fft_mesh/calc_ceikr_dpc [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

m_fft_mesh/calc_ceikr_spc [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

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.

SOURCE

1317 pure subroutine calc_eigr(gg, nfft, ngfft, eigr)
1318 
1319 !Arguments ------------------------------------
1320 !scalars
1321  integer,intent(in) :: nfft
1322 !arrays
1323  integer,intent(in) :: gg(3)
1324  integer,intent(in) :: ngfft(18)
1325  real(dp),intent(out) :: eigr(2*nfft)
1326 
1327 !Local variables-------------------------------
1328 !scalars
1329  integer :: ix,iy,iz,fft_idx
1330  real(dp) :: gdotr
1331 
1332 ! *************************************************************************
1333 
1334  if (ALL(gg==0)) then
1335    eigr(1:2*nfft:2)=one
1336    eigr(2:2*nfft:2)=zero
1337    RETURN
1338  end if
1339 
1340  fft_idx=1
1341  do iz=0,ngfft(3)-1
1342    do iy=0,ngfft(2)-1
1343      do ix=0,ngfft(1)-1
1344        gdotr= two_pi*( gg(1)*(ix/DBLE(ngfft(1))) &
1345                       +gg(2)*(iy/DBLE(ngfft(2))) &
1346                       +gg(3)*(iz/DBLE(ngfft(3))) )
1347        eigr(fft_idx  )=DCOS(gdotr)
1348        eigr(fft_idx+1)=DSIN(gdotr)
1349        fft_idx = fft_idx + 2
1350      end do
1351    end do
1352  end do
1353 
1354 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

635 pure function check_rot_fft(nsym,symrel,nr1,nr2,nr3)
636 
637 !Arguments
638 !Scalar
639  integer,intent(in) :: nr1,nr2,nr3,nsym
640  logical :: check_rot_fft
641 !Arrays
642  integer,intent(in) :: symrel(3,3,nsym)
643 
644 !local variables
645  integer :: is
646 
647 !************************************************************************
648 
649  ! The grid is compatible with the symmetries (only rotational part) if
650  ! for each symmetry, each n_i and n_j ==> $n_i*R_{ij}/n_j$ is an integer
651  check_rot_fft=.TRUE.
652  do is=1,nsym
653    if ( MOD(symrel(2,1,is)*nr2, nr1) /=0 .or. &
654 &       MOD(symrel(3,1,is)*nr3, nr1) /=0 .or. &
655 &       MOD(symrel(1,2,is)*nr1, nr2) /=0 .or. &
656 &       MOD(symrel(3,2,is)*nr3, nr2) /=0 .or. &
657 &       MOD(symrel(1,3,is)*nr1, nr3) /=0 .or. &
658 &       MOD(symrel(2,3,is)*nr2, nr3) /=0      &
659 &     ) then
660      check_rot_fft=.FALSE.; EXIT
661    end if
662  end do
663 
664 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.

SOURCE

 948 subroutine cigfft(mG0,npwvec,ngfft,gvec,igfft,ierr)
 949 
 950 !Arguments ------------------------------------
 951 !scalars
 952  integer,intent(in) :: npwvec
 953  integer,intent(out) :: ierr
 954 !arrays
 955  integer,intent(in) :: gvec(3,npwvec)
 956  integer,intent(in) :: mg0(3),ngfft(18)
 957  integer,intent(out) :: igfft(npwvec,2*mg0(1)+1,2*mg0(2)+1,2*mg0(3)+1)
 958 
 959 !Local variables ------------------------------
 960 !scalars
 961  integer :: gmg01,gmg02,gmg03,ig,ig01,ig02,ig03,n1,n2,n3
 962  character(len=500) :: msg
 963 !arrays
 964  integer :: gmg0(3)
 965 !************************************************************************
 966 
 967  DBG_ENTER("COLL")
 968 
 969  if (ANY(mg0<0)) then
 970    ABI_BUG(sjoin('Found negative value of mg0:', trim(ltoa(mg0))))
 971  end if
 972 
 973  n1=ngfft(1)
 974  n2=ngfft(2)
 975  n3=ngfft(3)
 976  ierr=0
 977 
 978  do ig=1,npwvec
 979    do ig01=-mg0(1),mG0(1)
 980      gmg0(1) = gvec(1,ig)-ig01
 981      do ig02=-mg0(2),mg0(2)
 982        gmg0(2) = gvec(2,ig)-ig02
 983        do ig03=-mg0(3),mg0(3)
 984          gmg0(3) = gvec(3,ig)-ig03
 985          ! === Calculate FFT index of G-G0 ===
 986          ! * Consider possible wrap around errors.
 987          gmg01=MODULO(gmg0(1),n1)
 988          gmg02=MODULO(gmg0(2),n2)
 989          gmg03=MODULO(gmg0(3),n3)
 990          igfft(ig,ig01+mg0(1)+1,ig02+mg0(2)+1,ig03+mg0(3)+1) = 1+gmg01+gmg02*n1+gmg03*n1*n2
 991          if ( ANY(gmg0>ngfft(1:3)/2) .or. ANY(gmg0<-(ngfft(1:3)-1)/2) ) then
 992            igfft(ig,ig01+mg0(1)+1,ig02+mg0(2)+1,ig03+mg0(3)+1) = 0
 993            ierr=ierr+1
 994          end if
 995        end do
 996      end do
 997    end do
 998  end do !ig
 999 
1000  if (ierr/=0) then
1001    write(msg,'(a,i0,3a)')&
1002     'Found ',ierr,' G-G0 vectors falling outside the FFT box. ',ch10,&
1003     'igfft will be set to zero for these particular G-G0 '
1004    ABI_WARNING(msg)
1005  end if
1006 
1007  DBG_EXIT("COLL")
1008 
1009 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.

SOURCE

689 function fft_check_rotrans(nsym,symrel,tnons,ngfft,err) result(isok)
690 
691 !Arguments ------------------------------------
692 !scalars
693  integer,intent(in) :: nsym
694  logical :: isok
695 !arrays
696  integer,intent(in) :: symrel(3,3,nsym)
697  integer,intent(in) :: ngfft(18)
698  real(dp),intent(in) :: tnons(3,nsym)
699  real(dp),intent(out) :: err(3,nsym)
700 
701 !Local variables-------------------------------
702 !scalars
703  integer :: isym,ix,iy,iz,jx,jy,jz,ngfft1,ngfft2,ngfft3
704  !character(len=500) :: msg
705 !arrays
706  integer :: Rm1(3,3,nsym),r1_FFT(3),red2fft(3,3)
707  real(dp) :: Rm1_FFT(3,3,nsym),fft2red(3,3),r2_FFT(3),tnons_FFT(3,nsym)
708 
709 ! *************************************************************************
710 
711  ! === Precalculate R^-1 and fractional translations in FFT coordinates ===
712  ngfft1=ngfft(1)
713  ngfft2=ngfft(2)
714  ngfft3=ngfft(3)
715 
716  red2fft=RESHAPE((/ngfft1,0,0,0,ngfft2,0,0,0,ngfft3/),(/3,3/))
717  fft2red=RESHAPE((/(one/ngfft1),zero,zero,zero,(one/ngfft2),zero,zero,zero,(one/ngfft3)/),(/3,3/))
718  !
719  ! === For a fully compatible mesh, each Rm1_FFT should be integer ===
720  do isym=1,nsym
721    call mati3inv(symrel(:,:,isym),Rm1(:,:,isym))
722    Rm1(:,:,isym)=TRANSPOSE(Rm1(:,:,isym))
723    Rm1_FFT(:,:,isym)=MATMUL(MATMUL(red2fft,Rm1(:,:,isym)),fft2red)
724    tnons_FFT(:,isym)=MATMUL(red2fft,tnons(:,isym))
725  end do
726 
727  err(:,:)=smallest_real
728  do iz=0,ngfft3-1
729    R1_FFT(3)=DBLE(iz)
730    do iy=0,ngfft2-1
731      R1_FFT(2)=DBLE(iy)
732      do ix=0,ngfft1-1
733        R1_FFT(1)=DBLE(ix)
734        do isym=1,nsym  ! Form R^-1 (r-\tau) in the FFT basis ===
735          R2_FFT(:)=MATMUL(Rm1_FFT(:,:,isym),R1_FFT(:)-tnons_FFT(:,isym))
736          jx=NINT(R2_FFT(1)); err(1,isym)=MAX(err(1,isym),ABS(R2_FFT(1)-jx)/ngfft1)
737          jy=NINT(R2_FFT(2)); err(2,isym)=MAX(err(2,isym),ABS(R2_FFT(2)-jy)/ngfft2)
738          jz=NINT(R2_FFT(3)); err(3,isym)=MAX(err(3,isym),ABS(R2_FFT(3)-jz)/ngfft3)
739        end do
740      end do
741    end do
742  end do
743 
744  isok=.TRUE.
745  do isym=1,nsym
746    if (ANY(err(:,isym)>tol6)) then
747      isok=.FALSE.
748      !write(msg,'(a,i3,a,3es14.6)')' symmetry ',isym,') not compatible with FFT grid, error ',err(:,isym)
749      !ABI_WARNING(msg)
750    end if
751  end do
752 
753 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

1075 pure integer function g2ifft(gg,ngfft) result (gidx)
1076 
1077 !Arguments ------------------------------------
1078  integer,intent(in) :: gg(3),ngfft(3)
1079 
1080 !Local variables-------------------------------
1081 !scalars
1082  integer :: n1,n2,n3,ig1,ig2,ig3
1083 
1084 !************************************************************************
1085 
1086  ! Use the following indexing (N means ngfft of the adequate direction)
1087  ! 0 1 2 3 ... N/2    -(N-1)/2 ... -1    <= gg
1088  ! 1 2 3 4 ....N/2+1  N/2+2    ...  N    <= index
1089  !
1090  if ( ANY(gg>ngfft(1:3)/2) .or. ANY(gg<-(ngfft(1:3)-1)/2) ) then ! out of the box.
1091    gidx=0
1092  else
1093    n1=ngfft(1)
1094    n2=ngfft(2)
1095    n3=ngfft(3)
1096 
1097    ig1=MODULO(gg(1),n1)
1098    ig2=MODULO(gg(2),n2)
1099    ig3=MODULO(gg(3),n3)
1100 
1101    gidx = 1 + ig1 + n1*(ig2+ig3*n2)
1102  end if
1103 
1104 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))

SOURCE

1127 pure subroutine get_gftt(ngfft, kpt, gmet, gsq_max, gfft)
1128 
1129 !Arguments ------------------------------------
1130 !scalars
1131  real(dp),intent(out) :: gsq_max
1132 !arrays
1133  integer,intent(in) :: ngfft(18)
1134  integer,intent(out) :: gfft(3,ngfft(1)*ngfft(2)*ngfft(3))
1135  real(dp),intent(in) :: kpt(3),gmet(3,3)
1136 
1137 !Local variables-------------------------------
1138 !scalars
1139  integer :: ifft,g1,g2,g3,i1,i2,i3
1140  real(dp) :: dsq
1141 
1142 !************************************************************************
1143 
1144  ifft=0; gsq_max=smallest_real
1145  do i3=1,ngfft(3)
1146    g3 = ig2gfft(i3,ngfft(3))
1147    do i2=1,ngfft(2)
1148      g2 = ig2gfft(i2,ngfft(2))
1149      do i1=1,ngfft(1)
1150        g1 = ig2gfft(i1,ngfft(1))
1151        ifft = ifft+1
1152        gfft(1,ifft) = g1
1153        gfft(2,ifft) = g2
1154        gfft(3,ifft) = g3
1155        dsq=gmet(1,1)*(kpt(1)+dble(i1))**2 &
1156         +gmet(2,2)*(kpt(2)+dble(i2))**2 &
1157         +gmet(3,3)*(kpt(3)+dble(i3))**2 &
1158         +2._dp*(gmet(1,2)*(kpt(1)+dble(i1))*(kpt(2)+dble(i2)) &
1159         +gmet(2,3)*(kpt(2)+dble(i2))*(kpt(3)+dble(i3)) &
1160         +gmet(3,1)*(kpt(3)+dble(i3))*(kpt(1)+dble(i1)))
1161        gsq_max = MAX(dsq,gsq_max)
1162      end do
1163    end do
1164  end do
1165 
1166 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

1030 elemental integer function ig2gfft(ig,ng) result (gc)
1031 
1032 !Arguments ------------------------------------
1033 !scalars
1034  integer,intent(in) :: ig,ng
1035 
1036 !************************************************************************
1037 
1038  ! Use the following indexing (N means ngfft of the adequate direction)
1039  ! 0 1 2 3 ... N/2    -(N-1)/2 ... -1    <= gc
1040  ! 1 2 3 4 ....N/2+1  N/2+2    ...  N    <= index ig
1041  !
1042  if (ig<=0 .or. ig > ng) then
1043    ! Wrong ig, returns huge. Parent code will likely crash with SIGSEV.
1044    gc = huge(1)
1045    return
1046  end if
1047 
1048  if ( ig  > ng/2 + 1) then
1049    gc = ig - ng -1
1050  else
1051    gc = ig -1
1052  end if
1053 
1054 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.

SOURCE

1643 subroutine phase(ngfft, ph)
1644 
1645 !Arguments ------------------------------------
1646 !scalars
1647  integer,intent(in) :: ngfft
1648 !arrays
1649  real(dp),intent(out) :: ph(2*ngfft)
1650 
1651 !Local variables-------------------------------
1652 !scalars
1653  integer :: id,ig,nn
1654  real(dp) :: arg,fac
1655 
1656 ! *************************************************************************
1657 
1658  id=ngfft/2+2
1659  fac=pi/dble(ngfft)
1660  do ig=1,ngfft
1661    nn=ig-1-(ig/id)*ngfft
1662    arg=fac*dble(nn)
1663    ph(2*ig-1)=cos(arg)
1664    ph(2*ig)  =sin(arg)
1665 
1666  end do
1667 !XG 990504 Here zero the corresponding sine
1668  if((ngfft/2)*2==ngfft) ph(2*(id-1))=zero
1669 
1670 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)=Indices 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.

SOURCE

786 subroutine rotate_fft_mesh(nsym, symrel, tnons, ngfft, irottb, preserve)
787 
788 !Arguments ------------------------------------
789 !scalars
790  integer,intent(in) :: nsym
791  logical,intent(out) :: preserve
792 !arrays
793  integer,intent(in) :: symrel(3,3,nsym)
794  integer,intent(in) :: ngfft(18)
795  integer,intent(out) :: irottb(ngfft(1)*ngfft(2)*ngfft(3),nsym)
796  real(dp),intent(in) :: tnons(3,nsym)
797 
798 !Local variables-------------------------------
799 !scalars
800  integer :: ir1,isym,ix,iy,iz,jx,jy,jz,ngfft1,ngfft2,ngfft3
801  !character(len=500) :: msg
802 !arrays
803  integer :: Rm1(3,3,nsym),r1_FFT(3),red2fft(3,3)
804  real(dp) :: Rm1_FFT(3,3,nsym),err(3,nsym),fft2red(3,3),r2_FFT(3)
805  real(dp) :: tnons_FFT(3,nsym)
806 
807 ! *************************************************************************
808 
809  ! === Precalculate R^-1 and fractional translations in FFT coordinates ===
810  ngfft1=ngfft(1)
811  ngfft2=ngfft(2)
812  ngfft3=ngfft(3)
813 
814  red2fft=RESHAPE((/ngfft1,0,0,0,ngfft2,0,0,0,ngfft3/),(/3,3/))
815  fft2red=RESHAPE((/(one/ngfft1),zero,zero,zero,(one/ngfft2),zero,zero,zero,(one/ngfft3)/),(/3,3/))
816  !
817  ! === For a fully compatible mesh, each Rm1_FFT should be integer ===
818  do isym=1,nsym
819    call mati3inv(symrel(:,:,isym),Rm1(:,:,isym))
820    Rm1(:,:,isym)=TRANSPOSE(Rm1(:,:,isym))
821    Rm1_FFT(:,:,isym)=MATMUL(MATMUL(red2fft,Rm1(:,:,isym)),fft2red)
822    tnons_FFT(:,isym)=MATMUL(red2fft,tnons(:,isym))
823  end do
824 
825  err(:,:)=zero
826 
827 !$OMP PARALLEL DO PRIVATE(R1_FFT,ir1,R2_FFT,jx,jy,jz) reduction(MAX:err)
828  do iz=0,ngfft3-1
829    R1_FFT(3)=DBLE(iz)
830    do iy=0,ngfft2-1
831      R1_FFT(2)=DBLE(iy)
832      do ix=0,ngfft1-1
833        R1_FFT(1)=DBLE(ix)
834        ir1=1+ix+iy*ngfft1+iz*ngfft1*ngfft2
835        do isym=1,nsym
836          ! === Form R^-1 (r-\tau) in the FFT basis ===
837          R2_FFT(:)=MATMUL(Rm1_FFT(:,:,isym),R1_FFT(:)-tnons_FFT(:,isym))
838          jx=NINT(R2_FFT(1)); err(1,isym)=MAX(err(1,isym),ABS(R2_FFT(1)-jx)/ngfft1)
839          jy=NINT(R2_FFT(2)); err(2,isym)=MAX(err(2,isym),ABS(R2_FFT(2)-jy)/ngfft2)
840          jz=NINT(R2_FFT(3)); err(3,isym)=MAX(err(3,isym),ABS(R2_FFT(3)-jz)/ngfft3)
841          jx=MODULO(jx,ngfft1)
842          jy=MODULO(jy,ngfft2)
843          jz=MODULO(jz,ngfft3)
844          irottb(ir1,isym)=1+jx+jy*ngfft1+jz*ngfft1*ngfft2
845        end do
846      end do
847    end do
848  end do
849 
850  preserve=.TRUE.
851  do isym=1,nsym
852    if (ANY(err(:,isym)>tol6)) then
853      preserve=.FALSE.
854      !write(msg,'(a,i3,a,3es14.6)')' symmetry ',isym,') not compatible with FFT grid, error ',err(:,isym)
855      !ABI_WARNING(msg)
856    end if
857  end do
858 
859 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.

SOURCE

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

SOURCE

1495 pure subroutine times_eigr(gg, ngfft, nfft, ndat, ur)
1496 
1497 !Arguments ------------------------------------
1498 !scalars
1499  integer,intent(in) :: nfft,ndat
1500 !arrays
1501  integer,intent(in) :: gg(3)
1502  integer,intent(in) :: ngfft(18)
1503  real(dp),intent(inout) :: ur(2,nfft,ndat)
1504 
1505 !Local variables-------------------------------
1506 !scalars
1507  integer :: ix,iy,iz,ifft,idat
1508  real(dp) :: gr
1509 !arrays
1510  real(dp) :: ph(2),val(2)
1511 
1512 ! *************************************************************************
1513 
1514  if (all(gg==0)) return
1515 
1516  do idat=1,ndat
1517    ifft = 0
1518    do iz=0,ngfft(3)-1
1519      do iy=0,ngfft(2)-1
1520        do ix=0,ngfft(1)-1
1521          ifft = ifft + 1
1522          gr = two_pi*(gg(1)*(ix/dble(ngfft(1))) &
1523                      +gg(2)*(iy/dble(ngfft(2))) &
1524                      +gg(3)*(iz/dble(ngfft(3))) )
1525          ph(1) = cos(gr); ph(2) = sin(gr)
1526          val(1) = ur(1,ifft,idat); val(2) = ur(2,ifft,idat)
1527 
1528          ur(1,ifft,idat) = ph(1) * val(1) - ph(2) * val(2)
1529          ur(2,ifft,idat) = ph(1) * val(2) + ph(2) * val(1)
1530        end do
1531      end do
1532    end do
1533  end do ! idat
1534 
1535 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,ndat)= contains u(r) in input. output: u(r) e^{ig.r} on the real-space FFT mesh.

SOURCE

1559 subroutine times_eikr(kk, ngfft, nfft, ndat, ur)
1560 
1561 !Arguments ------------------------------------
1562 !scalars
1563  integer,intent(in) :: nfft,ndat
1564 !arrays
1565  real(dp),intent(in) :: kk(3)
1566  integer,intent(in) :: ngfft(18)
1567  real(dp),intent(inout) :: ur(2,nfft,ndat)
1568 
1569 !Local variables-------------------------------
1570  integer :: ix,iy,iz,ifft,idat
1571  real(dp) :: kr, ph(2),val(2)
1572 
1573 ! *************************************************************************
1574 
1575  if (all(abs(kk) < tol12)) return
1576 
1577  !$OMP PARALLEL DO IF (ndat > 1) PRIVATE(ifft, kr, ph, val)
1578  do idat=1,ndat
1579    ifft = 0
1580    do iz=0,ngfft(3)-1
1581      do iy=0,ngfft(2)-1
1582        do ix=0,ngfft(1)-1
1583          ifft = ifft + 1
1584          kr = two_pi*(kk(1)*(ix/dble(ngfft(1))) &
1585                      +kk(2)*(iy/dble(ngfft(2))) &
1586                      +kk(3)*(iz/dble(ngfft(3))) )
1587          ph(1) = cos(kr); ph(2) = sin(kr)
1588          val(1) = ur(1,ifft,idat); val(2) = ur(2,ifft,idat)
1589 
1590          ur(1,ifft,idat) = ph(1) * val(1) - ph(2) * val(2)
1591          ur(2,ifft,idat) = ph(1) * val(2) + ph(2) * val(1)
1592        end do
1593      end do
1594    end do
1595  end do
1596 
1597 end subroutine times_eikr

m_fft_mesh/zpad_free [ Functions ]

[ Top ] [ m_fft_mesh ] [ Functions ]

NAME

  zpad_free

FUNCTION

INPUTS

OUTPUT

SOURCE

214 subroutine zpad_free(zpad)
215 
216 !Arguments ------------------------------------
217 !scalars
218  type(zpad_t),intent(inout) :: zpad
219 
220 ! *************************************************************************
221 
222  ABI_SFREE(zpad%zplane)
223  ABI_SFREE(zpad%linex2ifft_yz)
224 
225 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)>

SOURCE

135 subroutine zpad_init(zpad,nx,ny,nz,ldx,ldy,ldz,mgfft,gbound)
136 
137 !Arguments ------------------------------------
138 !scalars
139  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,mgfft
140  type(zpad_t),intent(out) :: zpad
141 !arrays
142  integer,intent(in) :: gbound(2*mgfft+8,2)
143 
144 !Local variables-------------------------------
145 !scalars
146  integer :: jj,g3_max,g3_min,gg3,ifft_g3,igb,g2min,g2max,nlinex
147 
148 ! *************************************************************************
149 
150  g3_min = gbound(3, 2)
151  g3_max = gbound(4, 2)
152 
153  zpad%n_zplanes = g3_max - g3_min + 1
154 
155  ABI_MALLOC(zpad%zplane,      (2, nz))
156  ABI_MALLOC(zpad%linex2ifft_yz, (2, nx*ny*nz))
157 
158  ! Loop over the z-planes intersecting the G-sphere.
159  nlinex = 0
160  do gg3=1,zpad%n_zplanes
161    !
162    if (gg3<=g3_max+1) then
163      ifft_g3 = gg3
164    else
165      ifft_g3 = gg3 + nz - zpad%n_zplanes ! Wrap around for negative gg3.
166    end if
167    !
168    ! Select the set of y for this z-plane.
169    igb=2*gg3+3
170    g2min = gbound(igb  ,2)
171    g2max = gbound(igb+1,2)
172 
173    zpad%zplane(1,gg3) = ifft_g3
174    zpad%zplane(2,gg3) = igb
175 
176    !(1:g2max+1,ifft_g3)     ! Positive g_y.
177    !(g2min+ny+1:ny,ifft_g3) ! Negative g_y.
178 
179    do jj=1,g2max+1
180      nlinex = nlinex + 1
181      zpad%linex2ifft_yz(1,nlinex) = jj
182      zpad%linex2ifft_yz(2,nlinex) = ifft_g3
183    end do
184 
185    do jj=g2min+ny+1,ny
186      nlinex = nlinex + 1
187      zpad%linex2ifft_yz(1,nlinex) = jj
188      zpad%linex2ifft_yz(2,nlinex) = ifft_g3
189    end do
190  end do
191 
192  zpad%nlinex = nlinex
193 
194  RETURN
195  ABI_UNUSED((/ldx,ldy,ldz/))
196 
197 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

 90  type,public :: zpad_t
 91 
 92    integer :: nlinex
 93    ! Total number of 1D transforms.
 94 
 95    integer :: n_zplanes
 96    ! Number of z-planes intersecting the sphere.
 97 
 98    integer,allocatable :: zplane(:,:)
 99    ! zplane(3,n_zplanes)
100    ! zplane(1,zpl) : mapping z-plane index -> FFT index_z
101    ! zplane(2,zpl) : mapping z-plane index -> igb index in array gbound
102 
103    integer,allocatable :: linex2ifft_yz(:,:)
104    ! linex2ifft_yz(2,nlinex)
105    ! mapping 1D-FFT -> (FFT_index_y, FFT index_z)
106 
107  end type zpad_t
108 
109  public :: zpad_init
110  public :: zpad_free

m_numeric_tools/denpot_project [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

FUNCTION

  Compute n(r) + n( $R^{-1}(r-\tau)$ in) / 2
  Mainly used with R = inversion to select the even/odd part under inversion

INPUTS

  cplex=1 for real, 2 for complex data.
  ngfft(3)=Mesh divisions of input array
  nspden=Number of density components.
  in_rhor(cplex * nfftot * nspden)=Input array
  one_symrel(3,3)= R operation
  tau(3)=Fractional translation.

OUTPUT

  out_rhor(cplex * nfftot * nspden)=Output array

SOURCE

884 subroutine denpot_project(cplex,  ngfft, nspden, in_rhor, one_symrel, one_tnons, out_rhor)
885 
886 !Arguments-------------------------------------------------------------
887 !scalars
888  integer,intent(in) :: cplex, nspden
889 !arrays
890  integer,intent(in) :: ngfft(18), one_symrel(3,3)
891  real(dp),intent(in) :: in_rhor(cplex, product(ngfft(1:3)), nspden)
892  real(dp),intent(in) :: one_tnons(3)
893  real(dp),intent(out) :: out_rhor(cplex, product(ngfft(1:3)), nspden)
894 
895 !Local variables--------------------------------------------------------
896 !scalars
897  integer,parameter :: nsym1 = 1, isgn = 1
898  integer :: ispden, ii, ifft, ifft_rot, nfft
899  logical :: preserve
900 !arrays
901  integer,allocatable :: irottb(:)
902 
903 ! *************************************************************************
904 
905  nfft = product(ngfft(1:3))
906  ABI_MALLOC(irottb, (nfft))
907 
908  call rotate_fft_mesh(nsym1, one_symrel, one_tnons, ngfft, irottb, preserve)
909  ABI_CHECK(preserve, "FFT mesh is not compatible with {R, tau}")
910 
911  do ispden=1,nspden
912    do ifft=1,nfft
913      ifft_rot = irottb(ifft)
914      do ii=1,cplex
915        out_rhor(cplex, ifft, ispden) = (in_rhor(cplex, ifft, ispden) + isgn * in_rhor(cplex, ifft_rot, ispden)) * half
916      end do
917    end do
918  end do
919 
920  ABI_FREE(irottb)
921 
922 end subroutine denpot_project