TABLE OF CONTENTS
- ABINIT/m_fft_mesh
- ABINIT/mkgrid_fft
- m_fft_mesh/calc_ceigr_dpc
- m_fft_mesh/calc_ceigr_spc
- m_fft_mesh/calc_ceikr_dpc
- m_fft_mesh/calc_ceikr_spc
- m_fft_mesh/calc_eigr
- m_fft_mesh/check_rot_fft
- m_fft_mesh/cigfft
- m_fft_mesh/fft_check_rotrans
- m_fft_mesh/g2ifft
- m_fft_mesh/get_gftt
- m_fft_mesh/ig2gfft
- m_fft_mesh/phase
- m_fft_mesh/rotate_FFT_mesh
- m_fft_mesh/setmesh
- m_fft_mesh/times_eigr
- m_fft_mesh/times_eikr
- m_fft_mesh/zpad_free
- m_fft_mesh/zpad_init
- m_fft_mesh/zpad_t
- m_numeric_tools/denpot_project
ABINIT/m_fft_mesh [ 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 ]
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