1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
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
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
|
!
! Copyright (C) 2011 Quantum ESPRESSO group
! 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 .
!
!
!=----------------------------------------------------------------------=
MODULE recvec_subs
!=----------------------------------------------------------------------=
! ... subroutines generating G-vectors and variables nl* needed to map
! ... G-vector components onto the FFT grid(s) in reciprocal space
!
USE kinds, ONLY : dp
USE fft_types, ONLY: fft_stick_index, fft_type_descriptor
USE fft_ggen, ONLY : fft_set_nl
!
PRIVATE
SAVE
PUBLIC :: ggen, ggens
!=----------------------------------------------------------------------=
CONTAINS
!=----------------------------------------------------------------------=
!
!-----------------------------------------------------------------------
SUBROUTINE ggen ( dfftp, gamma_only, at, bg, gcutm, ngm_g, ngm, &
g, gg, mill, ig_l2g, gstart, no_global_sort )
!----------------------------------------------------------------------
!
! This routine generates all the reciprocal lattice vectors
! contained in the sphere of radius gcutm. Furthermore it
! computes the indices nl which give the correspondence
! between the fft mesh points and the array of g vectors.
!
USE mp, ONLY: mp_rank, mp_size, mp_sum
USE constants, ONLY : eps8
!
IMPLICIT NONE
!
TYPE(fft_type_descriptor),INTENT(INOUT) :: dfftp
LOGICAL, INTENT(IN) :: gamma_only
REAL(DP), INTENT(IN) :: at(3,3), bg(3,3), gcutm
INTEGER, INTENT(IN) :: ngm_g
INTEGER, INTENT(INOUT) :: ngm
REAL(DP), INTENT(OUT) :: g(:,:), gg(:)
INTEGER, INTENT(OUT) :: mill(:,:), ig_l2g(:), gstart
! if no_global_sort is present (and it is true) G vectors are sorted only
! locally and not globally. In this case no global array needs to be
! allocated and sorted: saves memory and a lot of time for large systems.
!
LOGICAL, OPTIONAL, INTENT(IN) :: no_global_sort
!
! here a few local variables
!
REAL(DP) :: tx(3), ty(3), t(3)
REAL(DP), ALLOCATABLE :: tt(:)
INTEGER :: ngm_save, n1, n2, n3, ngm_offset, ngm_max, ngm_local
!
REAL(DP), ALLOCATABLE :: g2sort_g(:)
! array containing only g vectors for the current processor
INTEGER, ALLOCATABLE :: mill_unsorted(:,:)
! array containing all g vectors generators, on all processors
! (replicated data). When no_global_sort is present and .true.,
! only g-vectors for the current processor are stored
INTEGER, ALLOCATABLE :: igsrt(:), g2l(:)
!
INTEGER :: ni, nj, nk, i, j, k, ipol, ng, igl, indsw
INTEGER :: istart, jstart, kstart
INTEGER :: mype, npe
LOGICAL :: global_sort, is_local
INTEGER, ALLOCATABLE :: ngmpe(:)
!
global_sort = .TRUE.
IF( PRESENT( no_global_sort ) ) THEN
global_sort = .NOT. no_global_sort
END IF
!
IF( .NOT. global_sort ) THEN
ngm_max = ngm
ELSE
ngm_max = ngm_g
END IF
!
! save current value of ngm
!
ngm_save = ngm
!
ngm = 0
ngm_local = 0
!
! set the total number of fft mesh points and and initial value of gg
! The choice of gcutm is due to the fact that we have to order the
! vectors after computing them.
!
gg(:) = gcutm + 1.d0
!
! and computes all the g vectors inside a sphere
!
ALLOCATE( mill_unsorted( 3, ngm_save ) )
ALLOCATE( igsrt( ngm_max ) )
ALLOCATE( g2l( ngm_max ) )
ALLOCATE( g2sort_g( ngm_max ) )
!
g2sort_g(:) = 1.0d20
!
! allocate temporal array
!
ALLOCATE( tt( dfftp%nr3 ) )
!
! max miller indices (same convention as in module stick_set)
!
ni = (dfftp%nr1-1)/2
nj = (dfftp%nr2-1)/2
nk = (dfftp%nr3-1)/2
!
! gamma-only: exclude space with x < 0
!
IF ( gamma_only ) THEN
istart = 0
ELSE
istart = -ni
ENDIF
!
iloop: DO i = istart, ni
!
! gamma-only: exclude plane with x = 0, y < 0
!
IF ( gamma_only .and. i == 0 ) THEN
jstart = 0
ELSE
jstart = -nj
ENDIF
!
tx(1:3) = i * bg(1:3,1)
!
jloop: DO j = jstart, nj
!
IF ( .NOT. global_sort ) THEN
IF ( fft_stick_index( dfftp, i, j ) == 0 ) CYCLE jloop
is_local = .TRUE.
ELSE
IF ( dfftp%lpara .AND. fft_stick_index( dfftp, i, j ) == 0) THEN
is_local = .FALSE.
ELSE
is_local = .TRUE.
END IF
END IF
!
! gamma-only: exclude line with x = 0, y = 0, z < 0
!
IF ( gamma_only .and. i == 0 .and. j == 0 ) THEN
kstart = 0
ELSE
kstart = -nk
ENDIF
!
ty(1:3) = tx(1:3) + j * bg(1:3,2)
!
! compute all the norm square
!
DO k = kstart, nk
!
t(1) = ty(1) + k * bg(1,3)
t(2) = ty(2) + k * bg(2,3)
t(3) = ty(3) + k * bg(3,3)
tt(k-kstart+1) = t(1)**2 + t(2)**2 + t(3)**2
ENDDO
!
! save all the norm square within cutoff
!
DO k = kstart, nk
IF (tt(k-kstart+1) <= gcutm) THEN
ngm = ngm + 1
IF (ngm > ngm_max) CALL errore ('ggen 1', 'too many g-vectors', ngm)
IF ( tt(k-kstart+1) > eps8 ) THEN
g2sort_g(ngm) = tt(k-kstart+1)
ELSE
g2sort_g(ngm) = 0.d0
ENDIF
IF (is_local) THEN
ngm_local = ngm_local + 1
mill_unsorted( :, ngm_local ) = (/ i,j,k /)
g2l(ngm) = ngm_local
ELSE
g2l(ngm) = 0
ENDIF
ENDIF
ENDDO
ENDDO jloop
ENDDO iloop
IF (ngm /= ngm_max) &
CALL errore ('ggen', 'g-vectors missing !', abs(ngm - ngm_max))
!
igsrt(1) = 0
IF( .NOT. global_sort ) THEN
CALL hpsort_eps( ngm, g2sort_g, igsrt, eps8 )
ELSE
CALL hpsort_eps( ngm_g, g2sort_g, igsrt, eps8 )
END IF
DEALLOCATE( g2sort_g, tt )
IF( .NOT. global_sort ) THEN
!
! compute adeguate offsets in order to avoid overlap between
! g vectors once they are gathered on a single (global) array
!
mype = mp_rank( dfftp%comm )
npe = mp_size( dfftp%comm )
ALLOCATE( ngmpe( npe ) )
ngmpe = 0
ngmpe( mype + 1 ) = ngm
CALL mp_sum( ngmpe, dfftp%comm )
ngm_offset = 0
DO ng = 1, mype
ngm_offset = ngm_offset + ngmpe( ng )
END DO
DEALLOCATE( ngmpe )
!
END IF
ngm = 0
!
ngloop: DO ng = 1, ngm_max
!
IF (g2l(igsrt(ng))>0) THEN
! fetch the indices
i = mill_unsorted(1, g2l(igsrt(ng)))
j = mill_unsorted(2, g2l(igsrt(ng)))
k = mill_unsorted(3, g2l(igsrt(ng)))
!
ngm = ngm + 1
!
! Here map local and global g index !!! N.B: :
! the global G vectors arrangement depends on the number of processors
!
IF( .NOT. global_sort ) THEN
ig_l2g( ngm ) = ng + ngm_offset
ELSE
ig_l2g( ngm ) = ng
END IF
g(1:3, ngm) = i * bg (:, 1) + j * bg (:, 2) + k * bg (:, 3)
gg(ngm) = sum(g(1:3, ngm)**2)
ENDIF
ENDDO ngloop
DEALLOCATE( igsrt, g2l )
IF (ngm /= ngm_save) &
CALL errore ('ggen', 'g-vectors (ngm) missing !', abs(ngm - ngm_save))
!
! determine first nonzero g vector
!
IF (gg(1).le.eps8) THEN
gstart=2
ELSE
gstart=1
ENDIF
!
! Now set nl and nls with the correct fft correspondence
!
CALL fft_set_nl( dfftp, at, g, mill )
!
END SUBROUTINE ggen
!
!-----------------------------------------------------------------------
SUBROUTINE ggens( dffts, gamma_only, at, g, gg, mill, gcutms, ngms, &
gs, ggs )
!-----------------------------------------------------------------------
!
! Initialize number and indices of g-vectors for a subgrid,
! for exactly the same ordering as for the original FFT grid
!
!--------------------------------------------------------------------
!
IMPLICIT NONE
!
LOGICAL, INTENT(IN) :: gamma_only
TYPE (fft_type_descriptor), INTENT(INOUT) :: dffts
! primitive lattice vectors
REAL(dp), INTENT(IN) :: at(3,3)
! G-vectors in FFT grid
REAL(dp), INTENT(IN) :: g(:,:), gg(:)
! Miller indices for G-vectors of FFT grid
INTEGER, INTENT(IN) :: mill(:,:)
! cutoff for subgrid
REAL(DP), INTENT(IN):: gcutms
! Local number of G-vectors in subgrid
INTEGER, INTENT(OUT):: ngms
! Optionally: G-vectors and modules
REAL(DP), INTENT(INOUT), POINTER, OPTIONAL:: gs(:,:), ggs(:)
!
INTEGER :: i, ng, ngm
!
ngm = SIZE(gg)
ngms = dffts%ngm
IF ( ngms > ngm ) CALL errore ('ggens','wrong number of G-vectors',1)
!
IF ( PRESENT(gs) ) ALLOCATE ( gs(3,ngms) )
IF ( PRESENT(ggs)) ALLOCATE ( ggs(ngms) )
ng = 0
DO i = 1, ngm
IF ( gg(i) > gcutms ) exit
IF ( PRESENT(gs) ) gs (:,i) = g(:,i)
IF ( PRESENT(ggs)) ggs(i) = gg(i)
ng = i
END DO
IF ( ng /= ngms ) CALL errore ('ggens','mismatch in number of G-vectors',2)
!
CALL fft_set_nl ( dffts, at, g )
!
END SUBROUTINE ggens
!
!=----------------------------------------------------------------------=
END MODULE recvec_subs
!=----------------------------------------------------------------------=
|
