Branch data 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 : 55171 : 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 : 55171 : real(dp),allocatable :: vectg(:,:),work(:,:),workfft(:)
115 : :
116 : : ! *************************************************************************
117 : :
118 : : DBG_ENTER("COLL")
119 : :
120 : : !Tests
121 [ - + ]: 55171 : 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 [ + + ]: 1048249 : nfftc=pawfgr%nfftc;ngfftc(:)=pawfgr%ngfftc(:);nfftctot=ngfftc(1)*ngfftc(2)*ngfftc(3)
128 [ + + ]: 1048249 : 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 [ + + ]: 55171 : 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 [ + + ]: 51119 : if (optgrid==1) then
161 : :
162 [ - + + - : 58899 : ABI_MALLOC(work,(2,nfftc))
+ - + - -
+ - + ]
163 : :
164 : : ! First spin component
165 : : ! --------------------------------------------------------------
166 [ + + ]: 19633 : if (optout==0) then
167 : : ! if optout=0, rhog on the fine grid is temporary (in vectg)
168 [ - + + - : 56835 : ABI_MALLOC(vectg,(2,nfftf))
+ - + - -
+ - + ]
169 [ + + + + ]: 2395552236 : vectg(:,:)=zero
170 [ + + ]: 18945 : if (optin==1) then
171 : : call zerosym(rhog,2,ngfftc(1),ngfftc(2),ngfftc(3),&
172 : 18748 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
173 [ + + - + ]: 18748 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
174 : : call indirect_parallel_Fourier&
175 [ + - - - : 1280 : & (pawfgr%coatofin,vectg,mpi_enreg,ngfftf,ngfftc,nfftf,nfftc,paral_kgb,rhog,nfftctot)
- - - + ]
176 : : else
177 [ + + ]: 235836195 : do i1=1,nfftc
178 [ + + ]: 707473649 : 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 [ - + + - : 56835 : ABI_MALLOC(workfft,(cplex*nfftf))
+ - + - -
+ - + ]
200 : 18945 : call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftf,1,ngfftf,0)
201 [ + + ]: 806287034 : rhorf(:,1)=workfft(:)
202 [ - + ]: 18945 : ABI_FREE(workfft)
203 [ - + ]: 18945 : 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 : 688 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
208 [ + + + + ]: 96501448 : rhogf(:,:)=zero
209 [ + - ]: 688 : if (optin==1) then
210 [ + + - + ]: 688 : 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 [ + + ]: 12298156 : do i1=1,nfftc
215 [ + + ]: 36893116 : 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 [ - + + - : 2064 : ABI_MALLOC(workfft,(cplex*nfftf))
+ - + - -
+ - + ]
237 : 688 : call fourdp(cplex,rhogf,workfft,1,mpi_enreg,nfftf,1,ngfftf,0)
238 [ + + ]: 32167608 : rhorf(:,1)=workfft(:)
239 [ - + ]: 688 : ABI_FREE(workfft)
240 : : end if
241 : :
242 : : ! Additional spin components
243 : : ! ----------------------------------------------------
244 [ + + ]: 19633 : if (nspden>=2) then
245 [ - + + - : 11178 : ABI_MALLOC(vectg,(2,nfftf))
+ - + - -
+ - + ]
246 [ + + ]: 9632 : do ispden=2,nspden
247 [ + + + + ]: 412329947 : vectg(:,:)=zero
248 [ - + + - : 17718 : ABI_MALLOC(workfft,(cplex*nfftc))
+ - + - -
+ - + ]
249 [ + + ]: 39602380 : workfft(:)=rhor(:,ispden)
250 : 5906 : call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftc,1,ngfftc,0)
251 [ - + ]: 5906 : ABI_FREE(workfft)
252 : : call zerosym(work,2,ngfftc(1),ngfftc(2),ngfftc(3),&
253 : 5906 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
254 [ + + - + ]: 5906 : 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 [ + + ]: 39038384 : do i1=1,nfftc
259 [ + + ]: 117103764 : 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 [ - + + - : 17718 : ABI_MALLOC(workfft,(cplex*nfftf))
+ - + - -
+ - + ]
265 : 5906 : call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftf,1,ngfftf,0)
266 [ + + ]: 137587221 : rhorf(:,ispden)=workfft(:)
267 [ - + ]: 9632 : ABI_FREE(workfft)
268 : : end do
269 [ - + ]: 3726 : ABI_FREE(vectg)
270 : : end if
271 : :
272 [ - + ]: 19633 : 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 [ + - ]: 31486 : else if (optgrid==-1) then
281 : :
282 [ - + + - : 94458 : ABI_MALLOC(work,(2,nfftf))
+ - + - -
+ - + ]
283 : :
284 : : ! First spin component
285 : : ! --------------------------------------------------------------
286 [ + - ]: 31486 : if (optout==0) then
287 : : ! if optout=0, rhog on the fine grid is temporary (in vectg)
288 [ - + + - : 94458 : ABI_MALLOC(vectg,(2,nfftc))
+ - + - -
+ - + ]
289 [ + + + + ]: 2113129939 : vectg(:,:)=zero
290 [ - + ]: 31486 : 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 [ - + + - : 94458 : ABI_MALLOC(workfft,(cplex*nfftf))
+ - + - -
+ - + ]
296 [ + + ]: 2050491976 : workfft(:)=rhorf(:,1)
297 : 31486 : call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftf,1,ngfftf,0)
298 [ - + ]: 31486 : ABI_FREE(workfft)
299 [ + + - + ]: 31486 : if(mpi_enreg%nproc_fft > 1 .and. mpi_enreg%paral_kgb==1) then
300 : : call indirect_parallel_Fourier&
301 [ + - - - : 1636 : & (pawfgr%fintocoa,vectg,mpi_enreg,ngfftc,ngfftf,nfftc,nfftf,paral_kgb,work,nfftftot)
- - - + ]
302 : : else
303 [ + + ]: 1765026148 : do i1=1,nfftf
304 [ + + + + ]: 3016383826 : 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 : 31486 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
310 [ - + + - : 94458 : ABI_MALLOC(workfft,(cplex*nfftc))
+ - + - -
+ - + ]
311 : 31486 : call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftc,1,ngfftc,0)
312 [ + + ]: 713574054 : rhor(:,1)=workfft(:)
313 [ - + ]: 31486 : ABI_FREE(workfft)
314 [ - + ]: 31486 : 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 [ + + ]: 31486 : if (nspden>=2) then
347 [ - + + - : 24438 : ABI_MALLOC(vectg,(2,nfftc))
+ - + - -
+ - + ]
348 [ + + ]: 18456 : do ispden=2,nspden
349 [ + + + + ]: 211348778 : vectg(:,:)=zero
350 [ - + + - : 30930 : ABI_MALLOC(workfft,(cplex*nfftf))
+ - + - -
+ - + ]
351 [ + + ]: 238847431 : workfft(:)=rhorf(:,ispden)
352 : 10310 : call fourdp(cplex,work,workfft,-1,mpi_enreg,nfftf,1,ngfftf,0)
353 [ - + ]: 10310 : ABI_FREE(workfft)
354 [ + + - + ]: 10310 : 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 [ + + ]: 235099143 : do i1=1,nfftf
359 [ + + + + ]: 374748447 : 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 : 10310 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
364 [ - + + - : 30930 : ABI_MALLOC(workfft,(cplex*nfftc))
+ - + - -
+ - + ]
365 : 10310 : call fourdp(cplex,vectg,workfft,1,mpi_enreg,nfftc,1,ngfftc,0)
366 [ + + ]: 70537466 : rhor(:,ispden)=workfft(:)
367 [ - + ]: 18456 : ABI_FREE(workfft)
368 : : end do
369 [ - + ]: 8146 : ABI_FREE(vectg)
370 : : end if
371 : :
372 [ - + ]: 31486 : 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 : : !!***
|
