Line data Source code
1 : !!****m* ABINIT/m_fourier_interpol
2 : !! NAME
3 : !! m_fourier_interpol
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines used to perform a Fourier interpolation.
7 : !! Mainly used in PAW to interpol data from/to the coarse FFT grid from/to the fine FFT grid.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2018-2026 ABINIT group (FJ, MT, MG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
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 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
33 :
34 : implicit none
35 :
36 : private
37 :
38 : !public procedures.
39 : public :: transgrid ! Convert a density/potential from the coarse to the fine grid and vice versa
40 : public :: fourier_interpol ! Fourier interpolation. Just a wrapper for transgrid.
41 : public :: fourier_interpol_seq ! Fourier interpolation without MPI-FFT.
42 :
43 : CONTAINS !========================================================================================
44 : !!***
45 :
46 : !----------------------------------------------------------------------
47 :
48 : !!****f* m_fourier_interpol/transgrid
49 : !! NAME
50 : !! transgrid
51 : !!
52 : !! FUNCTION
53 : !! Convert a given density (or potential) from the coarse to the fine rectangular grid and vice versa
54 : !! Used in PAW calculations
55 : !!
56 : !! INPUTS
57 : !! cplex=1 if rhor[f] is real, 2 if rhor[f] is complex
58 : !! mpi_enreg=information about MPI parallelization
59 : !! nspden=number of spin-density components
60 : !! optgrid=+1 to go from the coarse grid towards the fine grid
61 : !! -1 to go from the fine grid towards the coarse grid
62 : !! optin= 0: input density/potential is taken from rhor(:,nspden)
63 : !! 1: input density/potential is taken from rhog(:) (ispden=1)
64 : !! and rhor(:,2:4) (ispden=2,3,4)
65 : !! optout= 0: output density/potential is given in r space in rhor(:,nspden)
66 : !! 1: output density/potential is given in r space in rhor(:,nspden)
67 : !! and in g space in rhog(:)
68 : !! pawfgr <type(paw_fgr_type)>=fine rectangular grid parameters
69 : !! %nfftc=number of points in the coarse FFT box
70 : !! %nfft =number of points in the fine FFT box
71 : !! %ngfftc(18)=all needed information about 3D FFT, for the coarse grid
72 : !! %ngfft(18) =all needed information about 3D FFT, for the fine grid
73 : !! %coatofin(nfftc)=Index of the points of the coarse grid on the fine grid
74 : !! %fintocoa(nfft) =Index of the points of the fine grid on the coarse grid
75 : !! %usefinegrid= 1 if a fine FFT grid is used (0 otherwise)
76 : !! if optgrid=+1 and optin=1:
77 : !! rhog(2,nfftc)=Fourier transform of input density/potential on the coarse grid
78 : !! if optgrid=-1 and optin=1:
79 : !! rhogf(2,nfftf)=Fourier transform of input density/potential on the fine grid
80 : !! if optgrid=+1
81 : !! rhor(cplex*nfftc,nspden)=input density/potential in r space on the coarse grid
82 : !! if optgrid=-1:
83 : !! rhorf(cplex*nfftf,nspden)=input density/potential in r space on the fine grid
84 : !!
85 : !! OUTPUT
86 : !! if optgrid=-1 and optout=1:
87 : !! rhog(2,nfftc)=Fourier transform of output density/potential on the coarse grid
88 : !! if optgrid=+1 and optout=1:
89 : !! rhogf(2,nfftf)=Fourier transform of output density/potential on the fine grid
90 : !! if optgrid=-1
91 : !! rhor(cplex*nfftc,nspden)=output density/potential in r space on the coarse grid
92 : !! if optgrid=+1:
93 : !! rhorf(cplex*nfftf,nspden)=output density/potential in r space on the fine grid
94 : !!
95 : !! SOURCE
96 :
97 52554 : subroutine transgrid(cplex,mpi_enreg,nspden,optgrid,optin,optout,paral_kgb,pawfgr,rhog,rhogf,rhor,rhorf)
98 :
99 : !Arguments ---------------------------------------------
100 : !scalars
101 : integer,intent(in) :: cplex,nspden,optgrid,optin,optout,paral_kgb
102 : type(MPI_type),intent(in) :: mpi_enreg
103 : type(pawfgr_type),intent(in) :: pawfgr
104 : !arrays
105 : real(dp),intent(inout) :: rhog(2,pawfgr%nfftc),rhogf(2,pawfgr%nfft)
106 : real(dp),intent(inout) :: rhor(cplex*pawfgr%nfftc,nspden),rhorf(cplex*pawfgr%nfft,nspden)
107 :
108 : !Local variables ---------------------------------------
109 : !scalars
110 : integer :: i1,ispden,nfftc,nfftctot,nfftf,nfftftot
111 : character(len=500) :: msg
112 : !arrays
113 : integer :: ngfftc(18),ngfftf(18)
114 52554 : real(dp),allocatable :: vectg(:,:),work(:,:),workfft(:)
115 :
116 : ! *************************************************************************
117 :
118 : DBG_ENTER("COLL")
119 :
120 : !Tests
121 52554 : if(pawfgr%nfft<pawfgr%nfftc) then
122 0 : write(msg,'(a,2(i0,1x))')' nfft (fine grid) must be >= nfft (coarse grid) while: ',pawfgr%nfft, pawfgr%nfftc
123 0 : ABI_ERROR(msg)
124 : end if
125 :
126 : !Store FFT dimensions
127 998526 : nfftc=pawfgr%nfftc;ngfftc(:)=pawfgr%ngfftc(:);nfftctot=ngfftc(1)*ngfftc(2)*ngfftc(3)
128 998526 : nfftf=pawfgr%nfft ;ngfftf(:)=pawfgr%ngfft (:);nfftftot=ngfftf(1)*ngfftf(2)*ngfftf(3)
129 :
130 : !If no fine FFT grid is used, this is only a simple transfer
131 52554 : if (pawfgr%usefinegrid==0) then
132 4052 : if (optgrid==1) then
133 43212698 : rhorf=rhor
134 1213748 : if (optout==1.and.optin==1) rhogf=rhog
135 4052 : if (optout==1.and.optin/=1) then
136 0 : ABI_MALLOC(workfft,(cplex*nfftc))
137 0 : workfft(:)=rhor(:,1)
138 0 : call fourdp(cplex,rhogf,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0)
139 0 : ABI_FREE(workfft)
140 : end if
141 : end if
142 4052 : if (optgrid==-1) then
143 0 : rhor=rhorf
144 0 : if (optout==1.and.optin==1) rhog=rhogf
145 0 : if (optout==1.and.optin/=1) then
146 0 : ABI_MALLOC(workfft,(cplex*nfftc))
147 0 : workfft(:)=rhorf(:,1)
148 0 : call fourdp(cplex,rhog,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0)
149 0 : ABI_FREE(workfft)
150 : end if
151 : end if
152 : return
153 : end if
154 :
155 : !====== FROM THE COARSE GRID TOWARDS THE FINE GRID =============
156 : !===============================================================
157 : !Calculate the FT of rhor to have it in the g space on the coarse grid
158 : !Transfer the FT of rhor on the coarse grid towards the fine grid
159 : !Then calculate the FT back to get rhorf on the fine grid
160 48502 : if (optgrid==1) then
161 :
162 56232 : ABI_MALLOC(work,(2,nfftc))
163 :
164 : ! First spin component
165 : ! --------------------------------------------------------------
166 18744 : if (optout==0) then
167 : ! if optout=0, rhog on the fine grid is temporary (in vectg)
168 54183 : ABI_MALLOC(vectg,(2,nfftf))
169 2176294246 : vectg(:,:)=zero
170 18061 : if (optin==1) then
171 : call zerosym(rhog,2,ngfftc(1),ngfftc(2),ngfftc(3),&
172 17864 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
173 17864 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
174 : call indirect_parallel_Fourier&
175 1180 : & (pawfgr%coatofin,vectg,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,rhog,nfftctot)
176 : else
177 228995203 : do i1=1,nfftc
178 686952241 : vectg(:,pawfgr%coatofin(i1))=rhog(:,i1)
179 : end do
180 : end if
181 : else
182 591 : ABI_MALLOC(workfft,(cplex*nfftc))
183 5516 : workfft(:)=rhor(:,1)
184 197 : call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0)
185 197 : ABI_FREE(workfft)
186 : call zerosym(work,2,ngfftc(1),ngfftc(2),ngfftc(3),&
187 197 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
188 197 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
189 : call indirect_parallel_Fourier&
190 0 : & (pawfgr%coatofin,vectg,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,work,nfftctot)
191 : else
192 5516 : do i1=1,nfftc
193 16154 : vectg(:,pawfgr%coatofin(i1))=work(:,i1)
194 : end do
195 : end if
196 : end if
197 : ! call zerosym(vectg,2,ngfftf(1),ngfftf(2),ngfftf(3),&
198 : ! & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
199 54183 : ABI_MALLOC(workfft,(cplex*nfftf))
200 18061 : call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftf,1,ngfftf,0)
201 733200448 : rhorf(:,1)=workfft(:)
202 18061 : ABI_FREE(workfft)
203 18061 : ABI_FREE(vectg)
204 : else
205 : ! if optout=1, rhog on the fine grid is saved
206 : call zerosym(rhog,2,ngfftc(1),ngfftc(2),ngfftc(3),&
207 683 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
208 84169667 : rhogf(:,:)=zero
209 683 : if (optin==1) then
210 683 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
211 : call indirect_parallel_Fourier&
212 12 : & (pawfgr%coatofin,rhogf,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,rhog,nfftctot)
213 : else
214 10772391 : do i1=1,nfftc
215 32315831 : rhogf(:,pawfgr%coatofin(i1))=rhog(:,i1)
216 : end do
217 : end if
218 : else
219 0 : ABI_MALLOC(workfft,(cplex*nfftc))
220 0 : workfft(:)=rhor(:,1)
221 0 : call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0)
222 0 : ABI_FREE(workfft)
223 : call zerosym(work,2,ngfftc(1),ngfftc(2),ngfftc(3),&
224 0 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
225 0 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
226 : call indirect_parallel_Fourier&
227 0 : & (pawfgr%coatofin,rhogf,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,work,nfftctot)
228 : else
229 0 : do i1=1,nfftc
230 0 : rhogf(:,pawfgr%coatofin(i1))=work(:,i1)
231 : end do
232 : end if
233 : end if
234 : ! call zerosym(rhogf,2,ngfftf(1),ngfftf(2),ngfftf(3),&
235 : ! & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
236 2049 : ABI_MALLOC(workfft,(cplex*nfftf))
237 683 : call fourdp(cplex,rhogf,workfft,1,mpi_enreg,nfftf,1,ngfftf,0)
238 28057011 : rhorf(:,1)=workfft(:)
239 683 : ABI_FREE(workfft)
240 : end if
241 :
242 : ! Additional spin components
243 : ! ----------------------------------------------------
244 18744 : if (nspden>=2) then
245 10863 : ABI_MALLOC(vectg,(2,nfftf))
246 9422 : do ispden=2,nspden
247 377493362 : vectg(:,:)=zero
248 17403 : ABI_MALLOC(workfft,(cplex*nfftc))
249 36161635 : workfft(:)=rhor(:,ispden)
250 5801 : call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0)
251 5801 : ABI_FREE(workfft)
252 : call zerosym(work,2,ngfftc(1),ngfftc(2),ngfftc(3),&
253 5801 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
254 5801 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
255 : call indirect_parallel_Fourier&
256 212 : & (pawfgr%coatofin,vectg,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,work,nfftctot)
257 : else
258 35597639 : do i1=1,nfftc
259 106781739 : vectg(:,pawfgr%coatofin(i1))=work(:,i1)
260 : end do
261 : end if
262 : ! call zerosym(vectg,2,ngfftf(1),ngfftf(2),ngfftf(3),&
263 : ! & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
264 17403 : ABI_MALLOC(workfft,(cplex*nfftf))
265 5801 : call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftf,1,ngfftf,0)
266 125974956 : rhorf(:,ispden)=workfft(:)
267 9422 : ABI_FREE(workfft)
268 : end do
269 3621 : ABI_FREE(vectg)
270 : end if
271 :
272 18744 : ABI_FREE(work)
273 :
274 :
275 : ! ====== FROM THE FINE GRID TOWARDS THE COARSE GRID =============
276 : ! ==============================================================
277 : ! Calculate the FT of rhorf to have it in the g space on the fine grid
278 : ! Transfer the FT of rhorf on the fine grid towards the coarse grid
279 : ! Then calculate the FT back to get rhor on the coarse grid
280 29758 : else if (optgrid==-1) then
281 :
282 89274 : ABI_MALLOC(work,(2,nfftf))
283 :
284 : ! First spin component
285 : ! --------------------------------------------------------------
286 29758 : if (optout==0) then
287 : ! if optout=0, rhog on the fine grid is temporary (in vectg)
288 89274 : ABI_MALLOC(vectg,(2,nfftc))
289 1394024899 : vectg(:,:)=zero
290 29758 : if (optin==1) then
291 0 : do i1=1,nfftf
292 0 : if (pawfgr%fintocoa(i1)/=0) vectg(:,pawfgr%fintocoa(i1))=rhogf(:,i1)
293 : end do
294 : else
295 89274 : ABI_MALLOC(workfft,(cplex*nfftf))
296 1400442766 : workfft(:)=rhorf(:,1)
297 29758 : call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftf,1,ngfftf,0)
298 29758 : ABI_FREE(workfft)
299 29758 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
300 : call indirect_parallel_Fourier&
301 1236 : & (pawfgr%fintocoa,vectg,mpi_enreg,ngfftc,ngfftf,nfftc,nfftf,paral_kgb,work,nfftftot)
302 : else
303 1314977338 : do i1=1,nfftf
304 2236232008 : if (pawfgr%fintocoa(i1)/=0) vectg(:,pawfgr%fintocoa(i1))=work(:,i1)
305 : end do
306 : end if
307 : end if
308 : call zerosym(vectg,2,ngfftc(1),ngfftc(2),ngfftc(3),&
309 29758 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
310 89274 : ABI_MALLOC(workfft,(cplex*nfftc))
311 29758 : call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftc,1,ngfftc,0)
312 473871222 : rhor(:,1)=workfft(:)
313 29758 : ABI_FREE(workfft)
314 29758 : ABI_FREE(vectg)
315 : else
316 : ! if optout=1, rhog on the fine grid is saved
317 0 : rhog(:,:)=zero
318 0 : if (optin==1) then
319 0 : do i1=1,nfftf
320 0 : if (pawfgr%fintocoa(i1)/=0) rhog(:,pawfgr%fintocoa(i1))=rhogf(:,i1)
321 : end do
322 : else
323 0 : ABI_MALLOC(workfft,(cplex*nfftf))
324 0 : workfft(:)=rhorf(:,1)
325 0 : call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftf,1,ngfftf,0)
326 0 : ABI_FREE(workfft)
327 0 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
328 : call indirect_parallel_Fourier&
329 0 : & (pawfgr%fintocoa,rhog,mpi_enreg,ngfftc,ngfftf,nfftc,nfftf,paral_kgb,work,nfftftot)
330 : else
331 0 : do i1=1,nfftf
332 0 : if (pawfgr%fintocoa(i1)/=0) rhog(:,pawfgr%fintocoa(i1))=work(:,i1)
333 : end do
334 : end if
335 : end if
336 : call zerosym(rhog,2,ngfftc(1),ngfftc(2),ngfftc(3),&
337 0 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
338 0 : ABI_MALLOC(workfft,(cplex*nfftc))
339 0 : call fourdp(cplex,rhog,workfft,1,mpi_enreg,nfftc,1,ngfftc,0)
340 0 : rhor(:,1)=workfft(:)
341 0 : ABI_FREE(workfft)
342 : end if
343 :
344 : ! Additional spin components
345 : ! ----------------------------------------------------
346 29758 : if (nspden>=2) then
347 23718 : ABI_MALLOC(vectg,(2,nfftc))
348 17976 : do ispden=2,nspden
349 187755578 : vectg(:,:)=zero
350 30210 : ABI_MALLOC(workfft,(cplex*nfftf))
351 212305111 : workfft(:)=rhorf(:,ispden)
352 10070 : call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftf,1,ngfftf,0)
353 10070 : ABI_FREE(workfft)
354 10070 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
355 : call indirect_parallel_Fourier&
356 256 : & (pawfgr%fintocoa,vectg,mpi_enreg,ngfftc,ngfftf,nfftc,nfftf,paral_kgb,work,nfftftot)
357 : else
358 208556823 : do i1=1,nfftf
359 332477487 : if (pawfgr%fintocoa(i1)/=0) vectg(:,pawfgr%fintocoa(i1))=work(:,i1)
360 : end do
361 : end if
362 : call zerosym(vectg,2,ngfftc(1),ngfftc(2),ngfftc(3),&
363 10070 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
364 30210 : ABI_MALLOC(workfft,(cplex*nfftc))
365 10070 : call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftc,1,ngfftc,0)
366 62672906 : rhor(:,ispden)=workfft(:)
367 17976 : ABI_FREE(workfft)
368 : end do
369 7906 : ABI_FREE(vectg)
370 : end if
371 :
372 29758 : ABI_FREE(work)
373 :
374 : end if
375 :
376 : DBG_EXIT("COLL")
377 :
378 : end subroutine transgrid
379 : !!***
380 :
381 : !----------------------------------------------------------------------
382 :
383 : !!****f* m_fourier_interpol/fourier_interpol
384 : !! NAME
385 : !! fourier_interpol
386 : !!
387 : !! FUNCTION
388 : !! Perform a Fourier interpolation. Just a wrapper for transgrid, the table giving the correspondence
389 : !! between the coarse and the mesh FFT grids are constructed inside the routine. This allows to specify an
390 : !! arbitrary FFT mesh to be used for the interpolation. Besides the routine works also in
391 : !! case of NC calculations since it does not require Pawfgr.
392 : !!
393 : !! INPUTS
394 : !! cplex=1 if rhor[f] is real, 2 if rhor[f] is complex
395 : !! MPI_enreg<MPI_type>=Information about MPI parallelization
396 : !! nspden=number of spin-density components
397 : !! nfft_in =number of points in the input FFT box (WARNING no FFT parallelism)
398 : !! nfft_out=number of points in the output FFT box
399 : !! ngfft_in(18)=all needed information about 3D FFT, for the input grid
400 : !! ngfft_out(18) =all needed information about 3D FFT, for the output grid
401 : !! optin= 0: input density/potential is taken from rhor_in(:,nspden)
402 : !! 1: input density/potential is taken from rhog_in(:) (ispden=1)
403 : !! and rhor_in(:,2:4) (ispden=2,3,4)
404 : !! optout= 0: output density/potential is given in r space in rhor_out(:,nspden)
405 : !! 1: output density/potential is given in r space in rhor_out(:,nspden)
406 : !! and in g space in rhog_out(:)
407 : !! ngfft_asked(18)=All info on the required FFT mesh.
408 : !!
409 : !! OUTPUT
410 : !! rhor_out(cplex*nfft_out,nspden)=output density/potential in r space on the required FFT mesh.
411 : !! if optout=1:
412 : !! rhog_out(2,nfftc)=Fourier transform of output density/potential on the coarse grid
413 : !!
414 : !! SOURCE
415 :
416 1732 : subroutine fourier_interpol(cplex,nspden,optin,optout,nfft_in,ngfft_in,nfft_out,ngfft_out,&
417 1732 : MPI_enreg,rhor_in,rhor_out,rhog_in,rhog_out)
418 :
419 : !Arguments ------------------------------------
420 : !scalars
421 : integer,intent(in) :: cplex,nspden,optin,optout
422 : integer,intent(in) :: nfft_in,nfft_out
423 : type(MPI_type),intent(in) :: MPI_enreg
424 : !arrays
425 : integer,intent(in) :: ngfft_in(18),ngfft_out(18)
426 : real(dp),intent(inout) :: rhor_in(cplex*nfft_in,nspden)
427 : real(dp),intent(inout) :: rhog_in(2,nfft_in)
428 : real(dp),intent(out) :: rhor_out(cplex*nfft_out,nspden)
429 : real(dp),intent(out) :: rhog_out(2,nfft_out)
430 :
431 : !Local variables ---------------------------------------
432 : !scalars
433 : integer :: nfftf,nfftc,nfftc_tot,nfftf_tot,optgrid
434 : logical :: ltest
435 : type(Pawfgr_type) :: Pawfgr
436 : !arrays
437 : integer :: ngfftc(18),ngfftf(18)
438 :
439 : ! *************************************************************************
440 :
441 22516 : ltest= ALL(ngfft_in(7:) == ngfft_out(7:) )
442 1732 : ABI_CHECK(ltest,'ngfftf_in(7:18)/=ngfftf_out(7:18)')
443 :
444 : !================================
445 : !=== Which one is the coarse? ===
446 : !================================
447 1732 : if (nfft_out>=nfft_in) then
448 : ! From coarse to fine grid. If meshes are equivalent, call transgrid anyway because of optout, optin.
449 0 : nfftf =nfft_out
450 0 : ngfftf(:)=ngfft_out(:)
451 0 : nfftf_tot =PRODUCT(ngfft_out(1:3))
452 :
453 0 : nfftc =nfft_in
454 0 : ngfftc(:)=ngfft_in(:)
455 0 : nfftc_tot =PRODUCT(ngfft_in (1:3))
456 :
457 0 : Pawfgr%usefinegrid=1
458 0 : if (ALL(ngfft_in(1:3)==ngfft_out(1:3))) Pawfgr%usefinegrid=0
459 0 : optgrid=1
460 :
461 : else
462 : ! From fine towards coarse.
463 1732 : nfftf =nfft_in
464 1732 : ngfftf(:)=ngfft_in(:)
465 6928 : nfftf_tot =PRODUCT(ngfft_in(1:3))
466 :
467 1732 : nfftc =nfft_out
468 1732 : ngfftc(:)=ngfft_out(:)
469 6928 : nfftc_tot =PRODUCT(ngfft_out (1:3))
470 1732 : Pawfgr%usefinegrid=1
471 1732 : optgrid=-1
472 : end if
473 :
474 5196 : ABI_MALLOC(Pawfgr%coatofin,(nfftc_tot))
475 5196 : ABI_MALLOC(Pawfgr%fintocoa,(nfftf_tot))
476 :
477 1732 : call indgrid(Pawfgr%coatofin,Pawfgr%fintocoa,nfftc_tot,nfftf_tot,ngfftc,ngfftf)
478 :
479 6928 : Pawfgr%mgfft =MAXVAL (ngfftf(1:3))
480 6928 : Pawfgr%nfft =PRODUCT(ngfftf(1:3)) ! no FFT parallelism!
481 32908 : Pawfgr%ngfft =ngfftf
482 :
483 6928 : Pawfgr%mgfftc =MAXVAL (ngfftc(1:3))
484 6928 : Pawfgr%nfftc =PRODUCT(ngfftc(1:3)) ! no FFT parallelism!
485 32908 : Pawfgr%ngfftc =ngfftc
486 :
487 1732 : if (optgrid==1) then
488 0 : call transgrid(cplex,MPI_enreg,nspden,optgrid,optin,optout,MPI_enreg%paral_kgb,Pawfgr,rhog_in ,rhog_out,rhor_in ,rhor_out)
489 : else
490 1732 : call transgrid(cplex,MPI_enreg,nspden,optgrid,optin,optout,MPI_enreg%paral_kgb,Pawfgr,rhog_out,rhog_in ,rhor_out,rhor_in)
491 : end if
492 :
493 1732 : call pawfgr_destroy(Pawfgr)
494 :
495 1732 : end subroutine fourier_interpol
496 : !!***
497 :
498 : !----------------------------------------------------------------------
499 :
500 : !!****f* m_fourier_interpol/fourier_interpol_seq
501 : !! NAME
502 : !! fourier_interpol_seq
503 : !!
504 : !! FUNCTION
505 : !! Perform a Fourier interpolation in sequential i.e. without MPI-FFT.
506 : !!
507 : !! INPUTS
508 : !! cplex=1 if rhor is real, 2 if rhor is complex
509 : !! nspden=number of spin-density components
510 : !! from_nfft =number of points in the input FFT box (WARNING no FFT parallelism)
511 : !! from_ngfft_(18)=all needed information about 3D FFT, for the input grid
512 : !! to_nfft=number of points in the output FFT box
513 : !! to_ngfft_out=all needed information about 3D FFT, for the output grid
514 : !! from_rhor=Input density
515 : !!
516 : !! OUTPUT
517 : !! to_rhor=Output density in r-space
518 : !! [to_rhog]=Output density in g-space.
519 : !!
520 : !! SOURCE
521 :
522 1732 : subroutine fourier_interpol_seq(cplex, nspden, from_nfft, from_ngfft, to_nfft, to_ngfft, from_rhor, to_rhor, &
523 0 : to_rhog) ! optional
524 :
525 : !Arguments ------------------------------------
526 : !scalars
527 : integer,intent(in) :: cplex, nspden, from_nfft, to_nfft
528 : !arrays
529 : integer,intent(in) :: from_ngfft(18), to_ngfft(18)
530 : real(dp),intent(inout) :: from_rhor(cplex*from_nfft, nspden)
531 : real(dp),intent(out) :: to_rhor(cplex*to_nfft, nspden)
532 : real(dp),optional,intent(out) :: to_rhog(2, to_nfft)
533 :
534 : !Local variables ---------------------------------------
535 : !scalars
536 1732 : type(MPI_type) :: MPI_enreg_seq
537 : !arrays
538 : integer :: ngfft_in(18),ngfft_out(18)
539 1732 : real(dp),allocatable :: from_rhog(:,:), to_rhog__(:,:)
540 : ! *************************************************************************
541 :
542 1732 : ngfft_in = from_ngfft; ngfft_out = to_ngfft
543 13856 : ngfft_in(4:6) = ngfft_in(1:3); ngfft_out(4:6) = ngfft_out(1:3)
544 : ! No MPI-FFT here.
545 38104 : ngfft_in(9:18) = 0; ngfft_out(9:18) = 0
546 1732 : ngfft_in(10) = 1; ngfft_out(10) = 1
547 :
548 1732 : call initmpi_seq(MPI_enreg_seq)
549 :
550 : ! Which one is coarse? Note that this part is not very robust and can fail!
551 1732 : if (ngfft_in(2) * ngfft_in(3) < ngfft_out(2) * ngfft_out(3)) then
552 0 : call MPI_enreg_seq%distribfft%init_seq('c', ngfft_in(2), ngfft_in(3), 'all')
553 0 : call MPI_enreg_seq%distribfft%init_seq('f', ngfft_out(2), ngfft_out(3), 'all')
554 : else
555 1732 : call MPI_enreg_seq%distribfft%init_seq('f', ngfft_in(2), ngfft_in(3), 'all')
556 1732 : call MPI_enreg_seq%distribfft%init_seq('c', ngfft_out(2), ngfft_out(3), 'all')
557 : end if
558 :
559 5196 : ABI_MALLOC(from_rhog, (2, from_nfft))
560 5196 : ABI_MALLOC(to_rhog__, (2, to_nfft))
561 :
562 : call fourier_interpol(cplex, nspden, 0, 0, from_nfft, ngfft_in, to_nfft, ngfft_out, &
563 1732 : MPI_enreg_seq, from_rhor, to_rhor, from_rhog, to_rhog__)
564 :
565 1732 : if (present(to_rhog)) to_rhog = to_rhog__
566 :
567 1732 : ABI_FREE(from_rhog)
568 1732 : ABI_FREE(to_rhog__)
569 1732 : call destroy_mpi_enreg(MPI_enreg_seq)
570 :
571 3464 : end subroutine fourier_interpol_seq
572 : !!***
573 :
574 : !----------------------------------------------------------------------
575 :
576 : END MODULE m_fourier_interpol
577 : !!***
|
