fft_interpolation.f90 6.5 KB
Newer Older
Marco Govoni's avatar
Marco Govoni committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
!
! Copyright (C) 2015-2016 M. Govoni 
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! This file is part of WEST.
!
! Contributors to this file: 
! He Ma, Marco Govoni
!
MODULE fourier_interpolation
 ! -----------------------------------------------------------------
 !
 IMPLICIT NONE
 !
 CONTAINS
 !
Marco Govoni's avatar
Marco Govoni committed
20
 SUBROUTINE get_G2R_mapping (n1, n2, n3, ng, ngx, ndim, nl)
Marco Govoni's avatar
Marco Govoni committed
21
22
    ! 
    ! INPUT  : n1, n2, n3 = dimension of arbitrary R grid 
Marco Govoni's avatar
Marco Govoni committed
23
24
    !          ng, ngx    = actual number of PW
    !          ndim       = dimensions of nl, 1: (n), 2:(n,-n)
Marco Govoni's avatar
Marco Govoni committed
25
26
27
28
29
30
31
32
33
34
    ! OUTPUT : nl         = computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] )
    !
    USE kinds,         ONLY : DP
    USE gvect,         ONLY : g
    USE cell_base,     ONLY : at, bg, tpiba
    !
    IMPLICIT NONE
    ! 
    ! I/O
    !
Marco Govoni's avatar
Marco Govoni committed
35
36
    INTEGER, INTENT(IN)  :: n1, n2, n3, ng, ngx, ndim
    INTEGER, INTENT(OUT) :: nl(ngx,ndim) ! mapping from [1, npw] to [1, n1*n2*n3]  
Marco Govoni's avatar
Marco Govoni committed
37
38
39
    !
    ! Workspace
    !
Marco Govoni's avatar
Marco Govoni committed
40
    INTEGER  :: ig, i_dim
Marco Govoni's avatar
Marco Govoni committed
41
42
43
44
45
46
47
48
49
50
    INTEGER  :: h, k, l, hmax, kmax, lmax, hidx, kidx, lidx
    !
    nl = 0
    !
    ! Maximum Miller indexes associated to the R grid 
    !
    hmax = (n1 - 1) / 2
    kmax = (n2 - 1) / 2
    lmax = (n3 - 1) / 2
    !
Marco Govoni's avatar
Marco Govoni committed
51
    DO i_dim = 1, ndim 
Marco Govoni's avatar
Marco Govoni committed
52
       DO ig = 1, ng
Marco Govoni's avatar
Marco Govoni committed
53
54
55
          !
          ! Get the current Miller index
          !
Marco Govoni's avatar
Marco Govoni committed
56
57
58
59
60
61
62
63
64
65
66
67
          SELECT CASE (i_dim)
          CASE(1)  
             h = NINT(SUM(g(:,ig) * at(:,1)))
             k = NINT(SUM(g(:,ig) * at(:,2)))
             l = NINT(SUM(g(:,ig) * at(:,3)))
          CASE(2)
             h = NINT(SUM(-g(:,ig) * at(:,1)))
             k = NINT(SUM(-g(:,ig) * at(:,2)))
             l = NINT(SUM(-g(:,ig) * at(:,3)))
          CASE DEFAULT
             CALL errore("","Bad mapping given: only 1 .OR. 2",1) 
          END SELECT
Marco Govoni's avatar
Marco Govoni committed
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
          !
          IF( ABS(h) <= hmax .AND. ABS(k) <= kmax .AND. ABS(l) <= lmax ) THEN 
             !
             ! derive position of G vector in (n1, n2, n3) grid
             !
             IF ( h >= 0 ) THEN
                hidx = h + 1
             ELSE
                hidx = h + n1 + 1
             ENDIF
             !
             IF ( k >= 0 ) THEN
                kidx = k + 1
             ELSE
                kidx = k + n2 + 1
             ENDIF
             !
             IF ( l >= 0 ) THEN
                lidx = l + 1
             ELSE
                lidx = l + n3 + 1
             ENDIF
             !
Marco Govoni's avatar
Marco Govoni committed
91
             nl(ig,i_dim) = hidx + (kidx - 1) * n1 + (lidx - 1) * n1 * n2
Marco Govoni's avatar
Marco Govoni committed
92
93
94
95
             !
          ENDIF
          !
       ENDDO
Marco Govoni's avatar
Marco Govoni committed
96
    ENDDO
Marco Govoni's avatar
Marco Govoni committed
97
98
99
100
    !
 END SUBROUTINE
 !
 !
Marco Govoni's avatar
Marco Govoni committed
101
 SUBROUTINE single_fwfft_fromArbitraryRGrid (fr, n1, n2, n3, ng, ngx, ndim, nl, fg, igk)
Marco Govoni's avatar
Marco Govoni committed
102
103
104
105
106
107
   !
   ! Note that the interpolation needs to be called by all processors within one band group 
   !
   ! FWFFT : R ---> G
   !
   ! INPUT  : n1, n2, n3 = dimension of arbitrary R grid 
Marco Govoni's avatar
Marco Govoni committed
108
109
   !          ng    = actual number of PW
   !          ngx   = leading dimendion for fg
Marco Govoni's avatar
Marco Govoni committed
110
111
112
113
114
115
116
117
118
119
120
121
   !          ndim  = 1, 2
   !          fr    = ONE COMPLEX array containing ONE function in R space (note that the array is not distributed, i.e. dimension = n1*n2*n3 )
   !          nl    = pre-computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] )
   ! OUTPUT : fg    = ONE COMPLEX array containing ONE functions in G space (note that the array is distributed ) 
   !
   USE kinds,         ONLY : DP
   USE fft_scalar,    ONLY : cfft3d
   USE mp_bands,      ONLY : intra_bgrp_comm, me_bgrp
   USE mp,            ONLY : mp_bcast
   !
   ! I/O
   !
Marco Govoni's avatar
Marco Govoni committed
122
123
   INTEGER,     INTENT(IN)       :: n1, n2, n3, ng, ngx, ndim
   INTEGER,     INTENT(IN)       :: nl(ngx,ndim)
Marco Govoni's avatar
Marco Govoni committed
124
   COMPLEX(DP), INTENT(IN)       :: fr(n1*n2*n3)
Marco Govoni's avatar
Marco Govoni committed
125
126
   COMPLEX(DP), INTENT(OUT)      :: fg(ngx)
   INTEGER, INTENT(IN), OPTIONAL :: igk(ng)
Marco Govoni's avatar
Marco Govoni committed
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
   !
   ! Workspace
   !
   INTEGER :: ig, idx
   !
   ! FFT is serial...
   !
   IF ( me_bgrp == 0 ) THEN
      !
      CALL cfft3d( fr, n1, n2, n3, n1, n2, n3, 1, -1)
      !
   ENDIF
   !
   ! ... result is broadcast
   !
   CALL mp_bcast( fr, 0, intra_bgrp_comm )
   !
   IF( PRESENT(igk) ) THEN
Marco Govoni's avatar
Marco Govoni committed
145
      DO ig=1, ng
Marco Govoni's avatar
Marco Govoni committed
146
147
148
149
150
151
         !
         idx = nl(igk(ig),1)
         IF (idx > 0) fg(ig) = fr(nl(igk(ig),1))
         !
      ENDDO
   ELSE
Marco Govoni's avatar
Marco Govoni committed
152
      DO ig=1, ng
Marco Govoni's avatar
Marco Govoni committed
153
154
155
156
157
158
159
         !
         idx = nl(ig,1)
         IF (idx > 0) fg(ig) = fr(nl(ig,1))
         !
      ENDDO
   ENDIF
   !
Marco Govoni's avatar
Marco Govoni committed
160
   DO ig = (ng+1), ngx
Marco Govoni's avatar
Marco Govoni committed
161
162
163
164
165
166
167
168
      !
      fg(ig) = (0.0_DP,0.0_DP)
      !
   ENDDO
   !
 ENDSUBROUTINE
 !
 !
Marco Govoni's avatar
Marco Govoni committed
169
 SUBROUTINE single_invfft_toArbitraryRGrid (fr, n1, n2, n3, ng, ngx, ndim, nl, fg, igk)
Marco Govoni's avatar
Marco Govoni committed
170
171
172
   !
   ! Note that the interpolation needs to be called by all processors within one band group 
   !
Marco Govoni's avatar
Marco Govoni committed
173
   ! INVFFT : R <--- G
Marco Govoni's avatar
Marco Govoni committed
174
175
176
177
178
179
180
181
182
183
184
185
   !
   ! INPUT  : n1, n2, n3 = dimension of arbitrary R grid 
   !          n     = actual number of PW
   !          nx    = leading dimendion for fg
   !          ndim  = 1,2
   !          fr    = ONE COMPLEX array containing ONE function in R space (note that the array is not distributed, i.e. dimension = n1*n2*n3 )
   !          nl    = pre-computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] )
   ! OUTPUT : fg    = ONE COMPLEX array containing ONE functions in G space (note that the array is distributed ) 
   !
   USE kinds,         ONLY : DP
   USE fft_scalar,    ONLY : cfft3d
   USE mp_bands,      ONLY : intra_bgrp_comm, me_bgrp
186
   USE mp,            ONLY : mp_bcast, mp_sum
Marco Govoni's avatar
Marco Govoni committed
187
188
189
   !
   ! I/O
   !
Marco Govoni's avatar
Marco Govoni committed
190
191
   INTEGER,     INTENT(IN) :: n1, n2, n3, ng, ngx, ndim
   INTEGER,     INTENT(IN) :: nl(ng,ndim)
Marco Govoni's avatar
Marco Govoni committed
192
   COMPLEX(DP), INTENT(OUT) :: fr(n1*n2*n3)
Marco Govoni's avatar
Marco Govoni committed
193
194
   COMPLEX(DP), INTENT(IN):: fg(ngx)
   INTEGER, INTENT(IN), OPTIONAL :: igk(ng)
Marco Govoni's avatar
Marco Govoni committed
195
196
197
198
199
200
201
202
   !
   ! Workspace
   !
   INTEGER :: ig, idx
   !
   fr = 0._DP
   !
   IF( PRESENT(igk) ) THEN
Marco Govoni's avatar
Marco Govoni committed
203
      DO ig=1, ng
Marco Govoni's avatar
Marco Govoni committed
204
205
206
207
208
209
210
211
         !
         idx = nl(igk(ig),1)
         IF (idx > 0) THEN 
            fr(nl(igk(ig),1)) = fg(ig)
         ENDIF
         !
      ENDDO
   ELSE 
Marco Govoni's avatar
Marco Govoni committed
212
      DO ig=1, ng
Marco Govoni's avatar
Marco Govoni committed
213
214
215
216
217
218
219
220
         !
         idx = nl(ig,1)
         IF (idx > 0) THEN 
            fr(nl(ig,1)) = fg(ig)
         ENDIF
         !
      ENDDO
      IF ( ndim == 2 ) THEN
Marco Govoni's avatar
Marco Govoni committed
221
         DO ig=1, ng
Marco Govoni's avatar
Marco Govoni committed
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
            !
            idx = nl(ig,2)
            IF (idx > 0) THEN 
               fr(nl(ig,2)) = DCONJG(fg(ig))
            ENDIF
            !
         ENDDO
      ENDIF
   ENDIF
   !
   ! ... result is gathered
   !
   CALL mp_sum( fr, intra_bgrp_comm )
   !
   ! FFT is serial...
   !
   IF ( me_bgrp == 0 ) THEN
      !
      CALL cfft3d( fr, n1, n2, n3, n1, n2, n3, 1, 1)
      !
   ENDIF
   !
 ENDSUBROUTINE
 !
END MODULE