TABLE OF CONTENTS


ABINIT/m_fourier_interpol [ Modules ]

[ Top ] [ 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-2022 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