TABLE OF CONTENTS
ABINIT/m_fourier_interpol [ Modules ]
NAME
m_fourier_interpol
FUNCTION
This module contains routines used to perform a Fourier interpolation. Mainly used in PAW to interpol data from/to the coarse FFT grid from/to the fine FFT grid.
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (FJ, MT, MG) 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
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_fourier_interpol 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 29 use defs_abitypes, only : MPI_type 30 use m_fft, only : zerosym, indirect_parallel_Fourier, fourdp 31 use m_pawfgr, only : pawfgr_type,pawfgr_destroy,indgrid 32 33 implicit none 34 35 private 36 37 !public procedures. 38 public :: transgrid ! Convert a density/potential from the coarse to the fine grid and vice versa 39 public :: fourier_interpol ! Perform a Fourier interpolation. Just a wrapper for transgrid. 40 41 CONTAINS !========================================================================================
m_fourier_interpol/fourier_interpol [ Functions ]
[ Top ] [ m_fourier_interpol ] [ Functions ]
NAME
fourier_interpol
FUNCTION
Perform a Fourier interpolation. Just a wrapper for transgrid, the table giving the correspondence between the coarse and the mesh FFT grids are constructed inside the routine. This allows to specify an arbitrary FFT mesh to be used for the interpolation. Besides the routine works also in case of NC calculations since it does not require Pawfgr.
INPUTS
cplex=1 if rhor[f] is real, 2 if rhor[f] is complex MPI_enreg<MPI_type>=Information about MPI parallelization nspden=number of spin-density components nfft_in =number of points in the input FFT box (WARNING no FFT parallelism) nfft_out=number of points in the output FFT box ngfft_in(18)=all needed information about 3D FFT, for the input grid ngfft_out(18) =all needed information about 3D FFT, for the output grid optin= 0: input density/potential is taken from rhor_in(:,nspden) 1: input density/potential is taken from rhog_in(:) (ispden=1) and rhor_in(:,2:4) (ispden=2,3,4) optout= 0: output density/potential is given in r space in rhor_out(:,nspden) 1: output density/potential is given in r space in rhor_out(:,nspden) and in g space in rhog_out(:) ngfft_asked(18)=All info on the required FFT mesh.
OUTPUT
rhor_out(cplex*nfft_out,nspden)=output density/potential in r space on the required FFT mesh. if optout=1: rhog_out(2,nfftc)=Fourier transform of output density/potential on the coarse grid
SOURCE
414 subroutine fourier_interpol(cplex,nspden,optin,optout,nfft_in,ngfft_in,nfft_out,ngfft_out,& 415 MPI_enreg,rhor_in,rhor_out,rhog_in,rhog_out) 416 417 !Arguments ------------------------------------ 418 !scalars 419 integer,intent(in) :: cplex,nspden,optin,optout 420 integer,intent(in) :: nfft_in,nfft_out 421 type(MPI_type),intent(in) :: MPI_enreg 422 !arrays 423 integer,intent(in) :: ngfft_in(18),ngfft_out(18) 424 real(dp),intent(inout) :: rhor_in(cplex*nfft_in,nspden) 425 real(dp),intent(inout) :: rhog_in(2,nfft_in) 426 real(dp),intent(out) :: rhor_out(cplex*nfft_out,nspden) 427 real(dp),intent(out) :: rhog_out(2,nfft_out) 428 429 !Local variables --------------------------------------- 430 !scalars 431 integer :: nfftf,nfftc,nfftc_tot,nfftf_tot,optgrid 432 logical :: ltest 433 type(Pawfgr_type) :: Pawfgr 434 !arrays 435 integer :: ngfftc(18),ngfftf(18) 436 437 ! ************************************************************************* 438 439 ltest= ALL(ngfft_in(7:) == ngfft_out(7:) ) 440 ABI_CHECK(ltest,'ngfftf_in(7:18)/=ngfftf_out(7:18)') 441 442 !================================ 443 !=== Which one is the coarse? === 444 !================================ 445 if (nfft_out>=nfft_in) then 446 ! From coarse to fine grid. If meshes are equivalent, call transgrid anyway because of optout, optin. 447 nfftf =nfft_out 448 ngfftf(:)=ngfft_out(:) 449 nfftf_tot =PRODUCT(ngfft_out(1:3)) 450 451 nfftc =nfft_in 452 ngfftc(:)=ngfft_in(:) 453 nfftc_tot =PRODUCT(ngfft_in (1:3)) 454 455 Pawfgr%usefinegrid=1 456 if (ALL(ngfft_in(1:3)==ngfft_out(1:3))) Pawfgr%usefinegrid=0 457 optgrid=1 458 459 else 460 ! From fine towards coarse. 461 nfftf =nfft_in 462 ngfftf(:)=ngfft_in(:) 463 nfftf_tot =PRODUCT(ngfft_in(1:3)) 464 465 nfftc =nfft_out 466 ngfftc(:)=ngfft_out(:) 467 nfftc_tot =PRODUCT(ngfft_out (1:3)) 468 Pawfgr%usefinegrid=1 469 optgrid=-1 470 end if 471 472 ABI_MALLOC(Pawfgr%coatofin,(nfftc_tot)) 473 ABI_MALLOC(Pawfgr%fintocoa,(nfftf_tot)) 474 475 call indgrid(Pawfgr%coatofin,Pawfgr%fintocoa,nfftc_tot,nfftf_tot,ngfftc,ngfftf) 476 477 Pawfgr%mgfft =MAXVAL (ngfftf(1:3)) 478 Pawfgr%nfft =PRODUCT(ngfftf(1:3)) ! no FFT parallelism! 479 Pawfgr%ngfft =ngfftf 480 481 Pawfgr%mgfftc =MAXVAL (ngfftc(1:3)) 482 Pawfgr%nfftc =PRODUCT(ngfftc(1:3)) ! no FFT parallelism! 483 Pawfgr%ngfftc =ngfftc 484 485 if (optgrid==1) then 486 call transgrid(cplex,MPI_enreg,nspden,optgrid,optin,optout,MPI_enreg%paral_kgb,Pawfgr,rhog_in ,rhog_out,rhor_in ,rhor_out) 487 else 488 call transgrid(cplex,MPI_enreg,nspden,optgrid,optin,optout,MPI_enreg%paral_kgb,Pawfgr,rhog_out,rhog_in ,rhor_out,rhor_in) 489 end if 490 491 call pawfgr_destroy(Pawfgr) 492 493 end subroutine fourier_interpol
m_fourier_interpol/transgrid [ Functions ]
[ Top ] [ m_fourier_interpol ] [ Functions ]
NAME
transgrid
FUNCTION
Convert a given density (or potential) from the coarse to the fine rectangular grid and vice versa Used in PAW calculations
INPUTS
cplex=1 if rhor[f] is real, 2 if rhor[f] is complex mpi_enreg=information about MPI parallelization nspden=number of spin-density components optgrid=+1 to go from the coarse grid towards the fine grid -1 to go from the fine grid towards the coarse grid optin= 0: input density/potential is taken from rhor(:,nspden) 1: input density/potential is taken from rhog(:) (ispden=1) and rhor(:,2:4) (ispden=2,3,4) optout= 0: output density/potential is given in r space in rhor(:,nspden) 1: output density/potential is given in r space in rhor(:,nspden) and in g space in rhog(:) pawfgr <type(paw_fgr_type)>=fine rectangular grid parameters %nfftc=number of points in the coarse FFT box %nfft =number of points in the fine FFT box %ngfftc(18)=all needed information about 3D FFT, for the coarse grid %ngfft(18) =all needed information about 3D FFT, for the fine grid %coatofin(nfftc)=Index of the points of the coarse grid on the fine grid %fintocoa(nfft) =Index of the points of the fine grid on the coarse grid %usefinegrid= 1 if a fine FFT grid is used (0 otherwise) if optgrid=+1 and optin=1: rhog(2,nfftc)=Fourier transform of input density/potential on the coarse grid if optgrid=-1 and optin=1: rhogf(2,nfftf)=Fourier transform of input density/potential on the fine grid if optgrid=+1 rhor(cplex*nfftc,nspden)=input density/potential in r space on the coarse grid if optgrid=-1: rhorf(cplex*nfftf,nspden)=input density/potential in r space on the fine grid
OUTPUT
if optgrid=-1 and optout=1: rhog(2,nfftc)=Fourier transform of output density/potential on the coarse grid if optgrid=+1 and optout=1: rhogf(2,nfftf)=Fourier transform of output density/potential on the fine grid if optgrid=-1 rhor(cplex*nfftc,nspden)=output density/potential in r space on the coarse grid if optgrid=+1: rhorf(cplex*nfftf,nspden)=output density/potential in r space on the fine grid
SOURCE
95 subroutine transgrid(cplex,mpi_enreg,nspden,optgrid,optin,optout,paral_kgb,pawfgr,rhog,rhogf,rhor,rhorf) 96 97 !Arguments --------------------------------------------- 98 !scalars 99 integer,intent(in) :: cplex,nspden,optgrid,optin,optout,paral_kgb 100 type(MPI_type),intent(in) :: mpi_enreg 101 type(pawfgr_type),intent(in) :: pawfgr 102 !arrays 103 real(dp),intent(inout) :: rhog(2,pawfgr%nfftc),rhogf(2,pawfgr%nfft) 104 real(dp),intent(inout) :: rhor(cplex*pawfgr%nfftc,nspden),rhorf(cplex*pawfgr%nfft,nspden) 105 106 !Local variables --------------------------------------- 107 !scalars 108 integer :: i1,ispden,nfftc,nfftctot,nfftf,nfftftot 109 character(len=500) :: msg 110 !arrays 111 integer :: ngfftc(18),ngfftf(18) 112 real(dp),allocatable :: vectg(:,:),work(:,:),workfft(:) 113 114 ! ************************************************************************* 115 116 DBG_ENTER("COLL") 117 118 !Tests 119 if(pawfgr%nfft<pawfgr%nfftc) then 120 write(msg,'(a,2(i0,1x))')' nfft (fine grid) must be >= nfft (coarse grid) while: ',pawfgr%nfft, pawfgr%nfftc 121 ABI_ERROR(msg) 122 end if 123 124 !Store FFT dimensions 125 nfftc=pawfgr%nfftc;ngfftc(:)=pawfgr%ngfftc(:);nfftctot=ngfftc(1)*ngfftc(2)*ngfftc(3) 126 nfftf=pawfgr%nfft ;ngfftf(:)=pawfgr%ngfft (:);nfftftot=ngfftf(1)*ngfftf(2)*ngfftf(3) 127 128 !If no fine FFT grid is used, this is only a simple transfer 129 if (pawfgr%usefinegrid==0) then 130 if (optgrid==1) then 131 rhorf=rhor 132 if (optout==1.and.optin==1) rhogf=rhog 133 if (optout==1.and.optin/=1) then 134 ABI_MALLOC(workfft,(cplex*nfftc)) 135 workfft(:)=rhor(:,1) 136 call fourdp(cplex,rhogf,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0) 137 ABI_FREE(workfft) 138 end if 139 end if 140 if (optgrid==-1) then 141 rhor=rhorf 142 if (optout==1.and.optin==1) rhog=rhogf 143 if (optout==1.and.optin/=1) then 144 ABI_MALLOC(workfft,(cplex*nfftc)) 145 workfft(:)=rhorf(:,1) 146 call fourdp(cplex,rhog,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0) 147 ABI_FREE(workfft) 148 end if 149 end if 150 return 151 end if 152 153 !====== FROM THE COARSE GRID TOWARDS THE FINE GRID ============= 154 !=============================================================== 155 !Calculate the FT of rhor to have it in the g space on the coarse grid 156 !Transfer the FT of rhor on the coarse grid towards the fine grid 157 !Then calculate the FT back to get rhorf on the fine grid 158 if (optgrid==1) then 159 160 ABI_MALLOC(work,(2,nfftc)) 161 162 ! First spin component 163 ! -------------------------------------------------------------- 164 if (optout==0) then 165 ! if optout=0, rhog on the fine grid is temporary (in vectg) 166 ABI_MALLOC(vectg,(2,nfftf)) 167 vectg(:,:)=zero 168 if (optin==1) then 169 call zerosym(rhog,2,ngfftc(1),ngfftc(2),ngfftc(3),& 170 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 171 if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then 172 call indirect_parallel_Fourier& 173 & (pawfgr%coatofin,vectg,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,rhog,nfftctot) 174 else 175 do i1=1,nfftc 176 vectg(:,pawfgr%coatofin(i1))=rhog(:,i1) 177 end do 178 end if 179 else 180 ABI_MALLOC(workfft,(cplex*nfftc)) 181 workfft(:)=rhor(:,1) 182 call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0) 183 ABI_FREE(workfft) 184 call zerosym(work,2,ngfftc(1),ngfftc(2),ngfftc(3),& 185 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 186 if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then 187 call indirect_parallel_Fourier& 188 & (pawfgr%coatofin,vectg,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,work,nfftctot) 189 else 190 do i1=1,nfftc 191 vectg(:,pawfgr%coatofin(i1))=work(:,i1) 192 end do 193 end if 194 end if 195 ! call zerosym(vectg,2,ngfftf(1),ngfftf(2),ngfftf(3),& 196 ! & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 197 ABI_MALLOC(workfft,(cplex*nfftf)) 198 call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftf,1,ngfftf,0) 199 rhorf(:,1)=workfft(:) 200 ABI_FREE(workfft) 201 ABI_FREE(vectg) 202 else 203 ! if optout=1, rhog on the fine grid is saved 204 call zerosym(rhog,2,ngfftc(1),ngfftc(2),ngfftc(3),& 205 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 206 rhogf(:,:)=zero 207 if (optin==1) then 208 if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then 209 call indirect_parallel_Fourier& 210 & (pawfgr%coatofin,rhogf,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,rhog,nfftctot) 211 else 212 do i1=1,nfftc 213 rhogf(:,pawfgr%coatofin(i1))=rhog(:,i1) 214 end do 215 end if 216 else 217 ABI_MALLOC(workfft,(cplex*nfftc)) 218 workfft(:)=rhor(:,1) 219 call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0) 220 ABI_FREE(workfft) 221 call zerosym(work,2,ngfftc(1),ngfftc(2),ngfftc(3),& 222 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 223 if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then 224 call indirect_parallel_Fourier& 225 & (pawfgr%coatofin,rhogf,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,work,nfftctot) 226 else 227 do i1=1,nfftc 228 rhogf(:,pawfgr%coatofin(i1))=work(:,i1) 229 end do 230 end if 231 end if 232 ! call zerosym(rhogf,2,ngfftf(1),ngfftf(2),ngfftf(3),& 233 ! & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 234 ABI_MALLOC(workfft,(cplex*nfftf)) 235 call fourdp(cplex,rhogf,workfft,1,mpi_enreg,nfftf,1,ngfftf,0) 236 rhorf(:,1)=workfft(:) 237 ABI_FREE(workfft) 238 end if 239 240 ! Additional spin components 241 ! ---------------------------------------------------- 242 if (nspden>=2) then 243 ABI_MALLOC(vectg,(2,nfftf)) 244 do ispden=2,nspden 245 vectg(:,:)=zero 246 ABI_MALLOC(workfft,(cplex*nfftc)) 247 workfft(:)=rhor(:,ispden) 248 call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0) 249 ABI_FREE(workfft) 250 call zerosym(work,2,ngfftc(1),ngfftc(2),ngfftc(3),& 251 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 252 if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then 253 call indirect_parallel_Fourier& 254 & (pawfgr%coatofin,vectg,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,work,nfftctot) 255 else 256 do i1=1,nfftc 257 vectg(:,pawfgr%coatofin(i1))=work(:,i1) 258 end do 259 end if 260 ! call zerosym(vectg,2,ngfftf(1),ngfftf(2),ngfftf(3),& 261 ! & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 262 ABI_MALLOC(workfft,(cplex*nfftf)) 263 call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftf,1,ngfftf,0) 264 rhorf(:,ispden)=workfft(:) 265 ABI_FREE(workfft) 266 end do 267 ABI_FREE(vectg) 268 end if 269 270 ABI_FREE(work) 271 272 273 ! ====== FROM THE FINE GRID TOWARDS THE COARSE GRID ============= 274 ! ============================================================== 275 ! Calculate the FT of rhorf to have it in the g space on the fine grid 276 ! Transfer the FT of rhorf on the fine grid towards the coarse grid 277 ! Then calculate the FT back to get rhor on the coarse grid 278 else if (optgrid==-1) then 279 280 ABI_MALLOC(work,(2,nfftf)) 281 282 ! First spin component 283 ! -------------------------------------------------------------- 284 if (optout==0) then 285 ! if optout=0, rhog on the fine grid is temporary (in vectg) 286 ABI_MALLOC(vectg,(2,nfftc)) 287 vectg(:,:)=zero 288 if (optin==1) then 289 do i1=1,nfftf 290 if (pawfgr%fintocoa(i1)/=0) vectg(:,pawfgr%fintocoa(i1))=rhogf(:,i1) 291 end do 292 else 293 ABI_MALLOC(workfft,(cplex*nfftf)) 294 workfft(:)=rhorf(:,1) 295 call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftf,1,ngfftf,0) 296 ABI_FREE(workfft) 297 if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then 298 call indirect_parallel_Fourier& 299 & (pawfgr%fintocoa,vectg,mpi_enreg,ngfftc,ngfftf,nfftc,nfftf,paral_kgb,work,nfftftot) 300 else 301 do i1=1,nfftf 302 if (pawfgr%fintocoa(i1)/=0) vectg(:,pawfgr%fintocoa(i1))=work(:,i1) 303 end do 304 end if 305 end if 306 call zerosym(vectg,2,ngfftc(1),ngfftc(2),ngfftc(3),& 307 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 308 ABI_MALLOC(workfft,(cplex*nfftc)) 309 call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftc,1,ngfftc,0) 310 rhor(:,1)=workfft(:) 311 ABI_FREE(workfft) 312 ABI_FREE(vectg) 313 else 314 ! if optout=1, rhog on the fine grid is saved 315 rhog(:,:)=zero 316 if (optin==1) then 317 do i1=1,nfftf 318 if (pawfgr%fintocoa(i1)/=0) rhog(:,pawfgr%fintocoa(i1))=rhogf(:,i1) 319 end do 320 else 321 ABI_MALLOC(workfft,(cplex*nfftf)) 322 workfft(:)=rhorf(:,1) 323 call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftf,1,ngfftf,0) 324 ABI_FREE(workfft) 325 if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then 326 call indirect_parallel_Fourier& 327 & (pawfgr%fintocoa,rhog,mpi_enreg,ngfftc,ngfftf,nfftc,nfftf,paral_kgb,work,nfftftot) 328 else 329 do i1=1,nfftf 330 if (pawfgr%fintocoa(i1)/=0) rhog(:,pawfgr%fintocoa(i1))=work(:,i1) 331 end do 332 end if 333 end if 334 call zerosym(rhog,2,ngfftc(1),ngfftc(2),ngfftc(3),& 335 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 336 ABI_MALLOC(workfft,(cplex*nfftc)) 337 call fourdp(cplex,rhog,workfft,1,mpi_enreg,nfftc,1,ngfftc,0) 338 rhor(:,1)=workfft(:) 339 ABI_FREE(workfft) 340 end if 341 342 ! Additional spin components 343 ! ---------------------------------------------------- 344 if (nspden>=2) then 345 ABI_MALLOC(vectg,(2,nfftc)) 346 do ispden=2,nspden 347 vectg(:,:)=zero 348 ABI_MALLOC(workfft,(cplex*nfftf)) 349 workfft(:)=rhorf(:,ispden) 350 call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftf,1,ngfftf,0) 351 ABI_FREE(workfft) 352 if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then 353 call indirect_parallel_Fourier& 354 & (pawfgr%fintocoa,vectg,mpi_enreg,ngfftc,ngfftf,nfftc,nfftf,paral_kgb,work,nfftftot) 355 else 356 do i1=1,nfftf 357 if (pawfgr%fintocoa(i1)/=0) vectg(:,pawfgr%fintocoa(i1))=work(:,i1) 358 end do 359 end if 360 call zerosym(vectg,2,ngfftc(1),ngfftc(2),ngfftc(3),& 361 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 362 ABI_MALLOC(workfft,(cplex*nfftc)) 363 call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftc,1,ngfftc,0) 364 rhor(:,ispden)=workfft(:) 365 ABI_FREE(workfft) 366 end do 367 ABI_FREE(vectg) 368 end if 369 370 ABI_FREE(work) 371 372 end if 373 374 DBG_EXIT("COLL") 375 376 end subroutine transgrid