set_npwq.f90 10.1 KB
Newer Older
1
!
Marco Govoni's avatar
Marco Govoni committed
2
! Copyright (C) 2015-2021 M. Govoni 
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
! 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: 
! Marco Govoni
!
!-----------------------------------------------------------------------
SUBROUTINE set_npwq()
  !-----------------------------------------------------------------------
  !
  USE kinds,           ONLY : DP
Marco Govoni's avatar
Marco Govoni committed
18
  USE westcom,         ONLY : npwq,npwq_g,npwqx,ngq,ngq_g,igq_q,l_use_ecutrho,fftdriver,dfft_io
19
20
21
22
23
24
25
  USE mp,              ONLY : mp_max, mp_sum
  USE mp_global,       ONLY : intra_bgrp_comm, inter_bgrp_comm, nbgrp, inter_pool_comm, intra_pool_comm
  USE gvect,           ONLY : ig_l2g,ngm,ngmx,g
  USE gvecw,           ONLY : gcutw
  USE pwcom,           ONLY : npw,npwx
  USE control_flags,   ONLY : gamma_only
  USE types_bz_grid,   ONLY : q_grid
Marco Govoni's avatar
Marco Govoni committed
26
  USE fft_base,        ONLY : dffts
27
28
29
30
31
32
33
  !
  IMPLICIT NONE
  !
  ! Workspace
  !
  INTEGER, EXTERNAL :: n_plane_waves
  REAL(DP), ALLOCATABLE :: gq2kin(:)
34
  INTEGER :: iq, ig
35
36
!  INTEGER :: npwqx_g
!  INTEGER, ALLOCATABLE :: igq_l2g(:)
37
38
39
40
41
42
43
  !
  IF ( gamma_only ) THEN
     !
     IF( l_use_ecutrho ) THEN
        npwq      = ngm
        npwqx     = ngmx
        fftdriver = 'Dense'
Marco Govoni's avatar
Marco Govoni committed
44
        dfft_io   = dffts
45
46
47
48
     ELSE
        npwq      = npw
        npwqx     = npwx
        fftdriver = 'Wave'
Marco Govoni's avatar
Marco Govoni committed
49
        CALL set_dfft_ecutwf( dfft_io )
50
51
52
53
54
55
56
57
58
     ENDIF
     !ALLOCATE(q0ig_l2g(npwq0))
     !q0ig_l2g(1:npwq0) = ig_l2g(1:npwq0)
     !npwq0_g=MAXVAL(q0ig_l2g(1:npwq0))
     npwq_g=MAXVAL(ig_l2g(1:npwq))
     CALL mp_max(npwq_g,intra_bgrp_comm)
     !
  ELSE
     !
59
     IF( l_use_ecutrho ) CALL errore("set_npwq", "Dense grid not implemented with q-points",1)
60
61
     fftdriver = 'Wave'
     !
Marco Govoni's avatar
Marco Govoni committed
62
63
     CALL set_dfft_ecutwf( dfft_io )
     !
64
     npwqx = n_plane_waves( gcutw, q_grid%np, q_grid%p_cart, g, ngm )
65
66
     !
     ALLOCATE( gq2kin(npwqx) )
67
68
     ALLOCATE( ngq(q_grid%np) )
     ALLOCATE( igq_q(npwqx,q_grid%np) )
69
     !ALLOCATE( igq_l2g(npwqx,q_grid%nps) )
70
     igq_q(:,:) = 0
71
     !igq_l2g(:,:) = 0
72
73
     DO iq = 1, q_grid%np
        CALL gq_sort( q_grid%p_cart(:,iq), ngm, g, gcutw, ngq(iq), igq_q(:,iq), gq2kin )
74
        !CALL gq_l2gmap( ngm, ig_l2g(1), ngq(iq), igq_q(1,iq), igq_l2g(1,iq) )
75
76
     ENDDO
     !
77
78
     DEALLOCATE(gq2kin)
     !
79
80
     ! ... compute the global number of q+G vectors for each q-point
     !
81
     ALLOCATE( ngq_g(q_grid%np) )
82
83
84
     !
     ngq_g = 0
     ngq_g(:) = ngq(:)
85
86
     !CALL mp_sum( ngq_g, inter_pool_comm )
     !CALL mp_sum( ngq_g, intra_pool_comm )
87
     CALL mp_sum( ngq_g, intra_bgrp_comm )
88
     !ngq_g = ngq_g / nbgrp
89
90
91
     !
     ! ... compute the maximum G vector index among all q+G in processors
     !
92
     npwq_g = 0 
93
     DO iq = 1, q_grid%np
94
95
96
97
98
        DO ig = 1, ngq(iq)
           npwq_g = MAX( npwq_g, ig_l2g(igq_q(ig,iq)) )
        ENDDO
     ENDDO
     !npwq_g = MAXVAL( igq_l2g(:,:) )
99
     !
100
     !CALL mp_max( npwq_g, intra_pool_comm )
101
102
103
104
     CALL mp_max( npwq_g, intra_bgrp_comm )
     !
     ! ... compute the maximum number of G vectors among all q-points
     !
105
!     npwqx_g = MAXVAL( ngq_g(1:q_grid%nps) )
106
107
108
     !
     ! ... define a further l2g map
     !
109
110
111
112
113
114
115
116
117
118
119
!     ALLOCATE( igq_l2g_kdip(npwqx_g, q_grid%nps) )
!     igq_l2g_kdip(:,:) = 0
!     !
!     DO iq = 1, q_grid%nps
!        ALLOCATE( igq_l2g(ngq(iq)) )
!        DO ig = 1, ngq(iq)
!           igq_l2g(ig) = ig_l2g( igq_q(ig,iq) ) 
!        ENDDO 
!        CALL gq_l2gmap_kdip( npwq_g, ngq_g(iq), ngq(iq), igq_l2g, igq_l2g_kdip(1,iq) )
!        DEALLOCATE( igq_l2g ) 
!     ENDDO
120
121
122
123
124
     !
  ENDIF
  !
END SUBROUTINE
!
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
!----------------------------------------------------------------------------
SUBROUTINE gq_sort( q, ngm, g, ecut, ngq, igq, gq )
   !----------------------------------------------------------------------------
   !
   ! ... sorts q+g in order of increasing magnitude, up to ecut
   ! ... NB: this version should yield the same ordering for different ecut
   ! ...     and the same ordering in all machines
   !
   USE kinds,     ONLY : DP
   USE constants, ONLY : eps8
   USE westcom,   ONLY : npwqx
   !
   IMPLICIT NONE
   !
   REAL(DP), INTENT(in) :: q(3)       ! the q point
   INTEGER, INTENT(in) :: ngm         ! the number of g vectors
   REAL(DP), INTENT(in) :: g(3,ngm)   ! the coordinates of G vectors
   REAL(DP), INTENT(in) :: ecut       ! the cut-off energy
   INTEGER, INTENT(out) :: ngq        ! the number of q+G vectors inside the "ecut sphere"
   INTEGER, INTENT(out) :: igq(npwqx) ! the correspondence q+G <-> G
   REAL(DP), INTENT(out) :: gq(npwqx) ! the moduli of q+G
   !
   INTEGER :: ng   ! counter on   G vectors
   INTEGER :: nk   ! counter on k+G vectors
   REAL(DP) :: qq  ! |k+G|^2
   REAL(DP) :: q2x ! upper bound for |G|
   !
   ! ... first we count the number of q+G vectors inside the cut-off sphere
   !
   q2x = ( SQRT( SUM(q(:)**2) ) + SQRT( ecut ) )**2
   !
   ngq = 0
   igq(:) = 0
   gq (:) = 0.0_dp
   !
   DO ng = 1, ngm
      qq = SUM( ( q(:) + g(:,ng) )**2 )
      IF ( qq <= eps8 ) qq = 0.d0
      !
      ! ... here if |q+G|^2 <= Ecut
      !
      IF ( qq <= ecut ) THEN
         ngq = ngq + 1
         IF ( ngq > npwqx ) &
            CALL errore( 'gq_sort', 'array gq out-of-bounds', 1 )
         !
         gq(ngq) = qq
         !
         ! set the initial value of index array
         igq(ngq) = ng
      ELSE
         ! if |G| > |q| + SQRT( Ecut )  stop search and order vectors
         IF ( SUM( g(:,ng)**2 ) > ( q2x + eps8 ) ) EXIT
      ENDIF
   ENDDO
   !
   IF ( ng > ngm ) &
      CALL infomsg( 'gq_sort', 'unexpected exit from do-loop')
   !
   ! ... order vector gq keeping initial position in index
   !
   CALL hpsort_eps( ngq, gq, igq, eps8 )
   !
   ! ... now order true |q+G|
   !
   DO nk = 1, ngq
      gq(nk) = SUM( (q(:) + g(:,igq(nk)) )**2 )
   ENDDO
   !
END SUBROUTINE
!
!----------------------------------------------------------------------------
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
!SUBROUTINE gq_l2gmap( ngm, ig_l2g, ngk, igk, igk_l2g )
!  !----------------------------------------------------------------------------
!  !
!  ! ... This subroutine maps local G+k index to the global G vector index
!  ! ... the mapping is used to collect wavefunctions subsets distributed
!  ! ... across processors.
!  ! ... Written by Carlo Cavazzoni
!  !
!  IMPLICIT NONE
!  !
!  ! ... Here the dummy variables
!  !
!  INTEGER, INTENT(IN)  :: ngm, ngk, igk(ngk), ig_l2g(ngm)
!  INTEGER, INTENT(OUT) :: igk_l2g(ngk)
!  INTEGER              :: ig
!  !
!  ! ... input: mapping between local and global G vector index
!  !
!  DO ig = 1, ngk
!     !
!     igk_l2g(ig) = ig_l2g(igk(ig))
!     !
!  END DO
!  !
!END SUBROUTINE
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
!
!-----------------------------------------------------------------------
SUBROUTINE gq_l2gmap_kdip( npw_g, ngk_g, ngk, igk_l2g, igk_l2g_kdip )
  !-----------------------------------------------------------------------
  !
  ! ... This subroutine maps local G+k index to the global G vector index
  ! ... the mapping is used to collect wavefunctions subsets distributed
  ! ... across processors.
  ! ... This map is used to obtained the G+k grids related to each kpt
  !
  USE mp_bands, ONLY : intra_bgrp_comm
  USE mp,       ONLY : mp_sum
  !
  IMPLICIT NONE
  !
  ! ... Here the dummy variables
  !
  INTEGER, INTENT(IN)  :: npw_g, ngk_g, ngk
  INTEGER, INTENT(IN)  :: igk_l2g(ngk)
  INTEGER, INTENT(OUT) :: igk_l2g_kdip(ngk)
  !
  INTEGER, ALLOCATABLE :: igwk_(:), itmp(:), igwk_lup(:)
  INTEGER              :: ig, ig_, ngg
  !
  !
  ALLOCATE( itmp( npw_g ) )
  ALLOCATE( igwk_( ngk_g ) )
  !
  itmp(:)  = 0
  !
  !
  DO ig = 1, ngk
     !
     itmp(igk_l2g(ig)) = igk_l2g(ig)
     !
  END DO
  !
  CALL mp_sum( itmp, intra_bgrp_comm )
  !
  ngg = 0
  DO ig = 1, npw_g
     !
     IF ( itmp(ig) == ig ) THEN
        !
        ngg = ngg + 1
        !
        igwk_(ngg) = ig
        !
     END IF
     !
  END DO
  !
  IF ( ngg /= ngk_g ) CALL errore( 'gk_l2gmap_kdip', 'unexpected dimension in ngg', 1 )
  !
  ALLOCATE( igwk_lup( npw_g ) )
  !
!$omp parallel private(ig_, ig)
!$omp workshare
  igwk_lup = 0
!$omp end workshare
!$omp do
  DO ig_ = 1, ngk_g
     igwk_lup(igwk_(ig_)) = ig_
  ENDDO
!$omp end do
!$omp do
  DO ig = 1, ngk
     igk_l2g_kdip(ig) = igwk_lup(igk_l2g(ig))
  ENDDO
!$omp end do
!$omp end parallel
  !
  DEALLOCATE( igwk_lup )
  !
  DEALLOCATE( itmp, igwk_ )
  !
  RETURN
  !
END SUBROUTINE
Marco Govoni's avatar
Marco Govoni committed
302
303
304
305
306
307
308
309
310
311
312
313
 !
 SUBROUTINE set_dfft_ecutwf (dfft)
    ! 
    ! OUTPUT  : dfft = FFT descriptor with real space set accorfing to ecutwfc
    !
    USE kinds,                ONLY : DP
    USE klist,                ONLY : xk, nks
    USE cell_base,            ONLY : at,bg
    USE control_flags,        ONLY : gamma_only
    USE fft_types,            ONLY : fft_type_descriptor, fft_type_init
    USE gvecw,                ONLY : gcutw
    USE mp,                   ONLY : mp_max
Victor Yu's avatar
Victor Yu committed
314
    USE mp_bands,             ONLY : ntask_groups
Marco Govoni's avatar
Marco Govoni committed
315
316
317
318
    USE mp_global,            ONLY : intra_bgrp_comm,inter_pool_comm
    USE stick_base,           ONLY : sticks_map
    USE gvecs,                ONLY : gcutms 
    !
Victor Yu's avatar
Victor Yu committed
319
320
321
    IMPLICIT NONE
    !
    !
Marco Govoni's avatar
Marco Govoni committed
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
    ! I/O
    !
    TYPE ( fft_type_descriptor ), INTENT(OUT) :: dfft ! customized fft descriptor
    !
    ! Workspace 
    !
    REAL(DP) :: gkcut
    TYPE( sticks_map ) :: smap
    INTEGER :: ik
    !
#if defined (__MPI)
  LOGICAL :: lpara = .true.
#else
  LOGICAL :: lpara = .false.
#endif
    !
    ! ... calculate gkcut = max |k+G|^2, in (2pi/a)^2 units
    !
    IF (nks == 0) THEN
       !
       ! if k-points are automatically generated (which happens later)
       ! use max(bg)/2 as an estimate of the largest k-point
       !
       gkcut = 0.5_DP * MAX ( &
          SQRT (SUM(bg (1:3, 1)**2) ), &
          SQRT (SUM(bg (1:3, 2)**2) ), &
          SQRT (SUM(bg (1:3, 3)**2) ) )
    ELSE
       gkcut = 0._DP
       DO ik = 1, nks
          gkcut = MAX (gkcut, SQRT ( SUM(xk (1:3, ik)**2) ) )
       ENDDO
    ENDIF
    !
    gkcut = (SQRT (gcutw) + gkcut)**2
    !
    ! ... find maximum value among all the processors
    !
    CALL mp_max (gkcut, inter_pool_comm )
    !
    CALL fft_type_init( dfft, smap, "wave", gamma_only, lpara, intra_bgrp_comm, at, bg, gkcut, MAX(gcutms/gkcut/4.0_DP,1.0_DP) )
    !
END SUBROUTINE