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