1! 2! Copyright (C) 2001-2007 PWSCF 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 becmod 11 ! 12 ! ... *bec* contain <beta|psi> - used in h_psi, s_psi, many other places 13 ! ... calbec( npw, beta, psi, betapsi [, nbnd ] ) is an interface calculating 14 ! ... betapsi(i,j) = <beta(i)|psi(j)> (the sum is over npw components) 15 ! ... or betapsi(i,s,j)= <beta(i)|psi(s,j)> (s=polarization index) 16 ! 17 USE kinds, ONLY : DP 18 USE control_flags, ONLY : gamma_only, smallmem 19 USE gvect, ONLY : gstart 20 USE noncollin_module, ONLY : noncolin, npol 21 ! 22 SAVE 23 ! 24 TYPE bec_type 25 REAL(DP), ALLOCATABLE :: r(:,:) ! appropriate for gammaonly 26 COMPLEX(DP),ALLOCATABLE :: k(:,:) ! appropriate for generic k 27 COMPLEX(DP),ALLOCATABLE :: nc(:,:,:) ! appropriate for noncolin 28 INTEGER :: comm 29 INTEGER :: nbnd 30 INTEGER :: nproc 31 INTEGER :: mype 32 INTEGER :: nbnd_loc 33 INTEGER :: ibnd_begin 34 END TYPE bec_type 35 ! 36 TYPE (bec_type) :: becp ! <beta|psi> 37 38 PRIVATE 39 ! 40 INTERFACE calbec 41 ! 42 MODULE PROCEDURE calbec_k, calbec_gamma, calbec_gamma_nocomm, calbec_nc, calbec_bec_type 43 ! 44 END INTERFACE 45 46 INTERFACE becscal 47 ! 48 MODULE PROCEDURE becscal_nck, becscal_gamma 49 ! 50 END INTERFACE 51 ! 52 PUBLIC :: bec_type, becp, allocate_bec_type, deallocate_bec_type, calbec, & 53 beccopy, becscal, is_allocated_bec_type 54 ! 55CONTAINS 56 !----------------------------------------------------------------------- 57 SUBROUTINE calbec_bec_type ( npw, beta, psi, betapsi, nbnd ) 58 !----------------------------------------------------------------------- 59 !_ 60 USE mp_bands, ONLY: intra_bgrp_comm 61 USE mp, ONLY: mp_get_comm_null 62 ! 63 IMPLICIT NONE 64 COMPLEX (DP), INTENT (in) :: beta(:,:), psi(:,:) 65 TYPE (bec_type), INTENT (inout) :: betapsi ! NB: must be INOUT otherwise 66 ! the allocatd array is lost 67 INTEGER, INTENT (in) :: npw 68 INTEGER, OPTIONAL :: nbnd 69 ! 70 INTEGER :: local_nbnd 71 INTEGER, EXTERNAL :: ldim_block, gind_block 72 INTEGER :: m_loc, m_begin, ip 73 REAL(DP), ALLOCATABLE :: dtmp(:,:) 74 ! 75 IF ( present (nbnd) ) THEN 76 local_nbnd = nbnd 77 ELSE 78 local_nbnd = size ( psi, 2) 79 ENDIF 80 81 IF ( gamma_only ) THEN 82 ! 83 IF( betapsi%comm == mp_get_comm_null() ) THEN 84 ! 85 CALL calbec_gamma ( npw, beta, psi, betapsi%r, local_nbnd, intra_bgrp_comm ) 86 ! 87 ELSE 88 ! 89 ALLOCATE( dtmp( SIZE( betapsi%r, 1 ), SIZE( betapsi%r, 2 ) ) ) 90 ! 91 DO ip = 0, betapsi%nproc - 1 92 m_loc = ldim_block( betapsi%nbnd , betapsi%nproc, ip ) 93 m_begin = gind_block( 1, betapsi%nbnd, betapsi%nproc, ip ) 94 IF( ( m_begin + m_loc - 1 ) > local_nbnd ) m_loc = local_nbnd - m_begin + 1 95 IF( m_loc > 0 ) THEN 96 CALL calbec_gamma ( npw, beta, psi(:,m_begin:m_begin+m_loc-1), dtmp, m_loc, betapsi%comm ) 97 IF( ip == betapsi%mype ) THEN 98 betapsi%r(:,1:m_loc) = dtmp(:,1:m_loc) 99 END IF 100 END IF 101 END DO 102 103 DEALLOCATE( dtmp ) 104 ! 105 END IF 106 ! 107 ELSEIF ( noncolin) THEN 108 ! 109 CALL calbec_nc ( npw, beta, psi, betapsi%nc, local_nbnd ) 110 ! 111 ELSE 112 ! 113 CALL calbec_k ( npw, beta, psi, betapsi%k, local_nbnd ) 114 ! 115 ENDIF 116 ! 117 RETURN 118 ! 119 END SUBROUTINE calbec_bec_type 120 !----------------------------------------------------------------------- 121 SUBROUTINE calbec_gamma_nocomm ( npw, beta, psi, betapsi, nbnd ) 122 !----------------------------------------------------------------------- 123 USE mp_bands, ONLY: intra_bgrp_comm 124 IMPLICIT NONE 125 COMPLEX (DP), INTENT (in) :: beta(:,:), psi(:,:) 126 REAL (DP), INTENT (out) :: betapsi(:,:) 127 INTEGER, INTENT (in) :: npw 128 INTEGER, OPTIONAL :: nbnd 129 INTEGER :: m 130 IF ( present (nbnd) ) THEN 131 m = nbnd 132 ELSE 133 m = size ( psi, 2) 134 ENDIF 135 CALL calbec_gamma ( npw, beta, psi, betapsi, m, intra_bgrp_comm ) 136 RETURN 137 ! 138 END SUBROUTINE calbec_gamma_nocomm 139 !----------------------------------------------------------------------- 140 SUBROUTINE calbec_gamma ( npw, beta, psi, betapsi, nbnd, comm ) 141 !----------------------------------------------------------------------- 142 ! 143 ! ... matrix times matrix with summation index (k=1,npw) running on 144 ! ... half of the G-vectors or PWs - assuming k=0 is the G=0 component: 145 ! ... betapsi(i,j) = 2Re(\sum_k beta^*(i,k)psi(k,j)) + beta^*(i,0)psi(0,j) 146 ! 147 USE mp, ONLY : mp_sum 148 149 IMPLICIT NONE 150 COMPLEX (DP), INTENT (in) :: beta(:,:), psi(:,:) 151 REAL (DP), INTENT (out) :: betapsi(:,:) 152 INTEGER, INTENT (in) :: npw 153 INTEGER, INTENT (in) :: nbnd 154 INTEGER, INTENT (in) :: comm 155 ! 156 INTEGER :: nkb, npwx, m 157 ! 158 m = nbnd 159 ! 160 nkb = size (beta, 2) 161 IF ( nkb == 0 ) RETURN 162 ! 163 CALL start_clock( 'calbec' ) 164 IF ( npw == 0 ) betapsi(:,:)=0.0_DP 165 npwx= size (beta, 1) 166 IF ( npwx /= size (psi, 1) ) CALL errore ('calbec', 'size mismatch', 1) 167 IF ( npwx < npw ) CALL errore ('calbec', 'size mismatch', 2) 168#if defined(DEBUG) 169 WRITE (*,*) 'calbec gamma' 170 WRITE (*,*) nkb, size (betapsi,1) , m , size (betapsi, 2) 171#endif 172 IF ( nkb /= size (betapsi,1) .or. m > size (betapsi, 2) ) & 173 CALL errore ('calbec', 'size mismatch', 3) 174 ! 175 IF ( m == 1 ) THEN 176 ! 177 CALL DGEMV( 'C', 2*npw, nkb, 2.0_DP, beta, 2*npwx, psi, 1, 0.0_DP, & 178 betapsi, 1 ) 179 IF ( gstart == 2 ) betapsi(:,1) = betapsi(:,1) - beta(1,:)*psi(1,1) 180 ! 181 ELSE 182 ! 183 CALL DGEMM( 'C', 'N', nkb, m, 2*npw, 2.0_DP, beta, 2*npwx, psi, & 184 2*npwx, 0.0_DP, betapsi, nkb ) 185 IF ( gstart == 2 ) & 186 CALL DGER( nkb, m, -1.0_DP, beta, 2*npwx, psi, 2*npwx, betapsi, nkb ) 187 ! 188 ENDIF 189 ! 190 CALL mp_sum( betapsi( :, 1:m ), comm ) 191 ! 192 CALL stop_clock( 'calbec' ) 193 ! 194 RETURN 195 ! 196 END SUBROUTINE calbec_gamma 197 ! 198 !----------------------------------------------------------------------- 199 SUBROUTINE calbec_k ( npw, beta, psi, betapsi, nbnd ) 200 !----------------------------------------------------------------------- 201 ! 202 ! ... matrix times matrix with summation index (k=1,npw) running on 203 ! ... G-vectors or PWs : betapsi(i,j) = \sum_k beta^*(i,k) psi(k,j) 204 ! 205 USE mp_bands, ONLY : intra_bgrp_comm 206 USE mp, ONLY : mp_sum 207 208 IMPLICIT NONE 209 COMPLEX (DP), INTENT (in) :: beta(:,:), psi(:,:) 210 COMPLEX (DP), INTENT (out) :: betapsi(:,:) 211 INTEGER, INTENT (in) :: npw 212 INTEGER, OPTIONAL :: nbnd 213 ! 214 INTEGER :: nkb, npwx, m 215 ! 216 nkb = size (beta, 2) 217 IF ( nkb == 0 ) RETURN 218 ! 219 CALL start_clock( 'calbec' ) 220 IF ( npw == 0 ) betapsi(:,:)=(0.0_DP,0.0_DP) 221 npwx= size (beta, 1) 222 IF ( npwx /= size (psi, 1) ) CALL errore ('calbec', 'size mismatch', 1) 223 IF ( npwx < npw ) CALL errore ('calbec', 'size mismatch', 2) 224 IF ( present (nbnd) ) THEN 225 m = nbnd 226 ELSE 227 m = size ( psi, 2) 228 ENDIF 229#if defined(DEBUG) 230 WRITE (*,*) 'calbec k' 231 WRITE (*,*) nkb, size (betapsi,1) , m , size (betapsi, 2) 232#endif 233 IF ( nkb /= size (betapsi,1) .or. m > size (betapsi, 2) ) & 234 CALL errore ('calbec', 'size mismatch', 3) 235 ! 236 IF ( m == 1 ) THEN 237 ! 238 CALL ZGEMV( 'C', npw, nkb, (1.0_DP,0.0_DP), beta, npwx, psi, 1, & 239 (0.0_DP, 0.0_DP), betapsi, 1 ) 240 ! 241 ELSE 242 ! 243 CALL ZGEMM( 'C', 'N', nkb, m, npw, (1.0_DP,0.0_DP), & 244 beta, npwx, psi, npwx, (0.0_DP,0.0_DP), betapsi, nkb ) 245 ! 246 ENDIF 247 ! 248 CALL mp_sum( betapsi( :, 1:m ), intra_bgrp_comm ) 249 ! 250 CALL stop_clock( 'calbec' ) 251 ! 252 RETURN 253 ! 254 END SUBROUTINE calbec_k 255 ! 256 !----------------------------------------------------------------------- 257 SUBROUTINE calbec_nc ( npw, beta, psi, betapsi, nbnd ) 258 !----------------------------------------------------------------------- 259 ! 260 ! ... matrix times matrix with summation index (k below) running on 261 ! ... G-vectors or PWs corresponding to two different polarizations: 262 ! ... betapsi(i,1,j) = \sum_k=1,npw beta^*(i,k) psi(k,j) 263 ! ... betapsi(i,2,j) = \sum_k=1,npw beta^*(i,k) psi(k+npwx,j) 264 ! 265 USE mp_bands, ONLY : intra_bgrp_comm 266 USE mp, ONLY : mp_sum 267 268 IMPLICIT NONE 269 COMPLEX (DP), INTENT (in) :: beta(:,:), psi(:,:) 270 COMPLEX (DP), INTENT (out) :: betapsi(:,:,:) 271 INTEGER, INTENT (in) :: npw 272 INTEGER, OPTIONAL :: nbnd 273 ! 274 INTEGER :: nkb, npwx, npol, m 275 ! 276 nkb = size (beta, 2) 277 IF ( nkb == 0 ) RETURN 278 ! 279 CALL start_clock ('calbec') 280 IF ( npw == 0 ) betapsi(:,:,:)=(0.0_DP,0.0_DP) 281 npwx= size (beta, 1) 282 IF ( 2*npwx /= size (psi, 1) ) CALL errore ('calbec', 'size mismatch', 1) 283 IF ( npwx < npw ) CALL errore ('calbec', 'size mismatch', 2) 284 IF ( present (nbnd) ) THEN 285 m = nbnd 286 ELSE 287 m = size ( psi, 2) 288 ENDIF 289 npol= size (betapsi, 2) 290#if defined(DEBUG) 291 WRITE (*,*) 'calbec nc' 292 WRITE (*,*) nkb, size (betapsi,1) , m , size (betapsi, 3) 293#endif 294 IF ( nkb /= size (betapsi,1) .or. m > size (betapsi, 3) ) & 295 CALL errore ('calbec', 'size mismatch', 3) 296 ! 297 CALL ZGEMM ('C', 'N', nkb, m*npol, npw, (1.0_DP, 0.0_DP), beta, & 298 npwx, psi, npwx, (0.0_DP, 0.0_DP), betapsi, nkb) 299 ! 300 CALL mp_sum( betapsi( :, :, 1:m ), intra_bgrp_comm ) 301 ! 302 CALL stop_clock( 'calbec' ) 303 ! 304 RETURN 305 ! 306 END SUBROUTINE calbec_nc 307 ! 308 ! 309 !----------------------------------------------------------------------- 310 FUNCTION is_allocated_bec_type (bec) RESULT (isalloc) 311 !----------------------------------------------------------------------- 312 IMPLICIT NONE 313 TYPE (bec_type) :: bec 314 LOGICAL :: isalloc 315 isalloc = (allocated(bec%r) .or. allocated(bec%nc) .or. allocated(bec%k)) 316 RETURN 317 ! 318 !----------------------------------------------------------------------- 319 END FUNCTION is_allocated_bec_type 320 !----------------------------------------------------------------------- 321 ! 322 !----------------------------------------------------------------------- 323 SUBROUTINE allocate_bec_type ( nkb, nbnd, bec, comm ) 324 !----------------------------------------------------------------------- 325 USE mp, ONLY: mp_size, mp_rank, mp_get_comm_null 326 IMPLICIT NONE 327 TYPE (bec_type) :: bec 328 INTEGER, INTENT (in) :: nkb, nbnd 329 INTEGER, INTENT (in), OPTIONAL :: comm 330 INTEGER :: ierr, nbnd_siz 331 INTEGER, EXTERNAL :: ldim_block, gind_block 332 ! 333 nbnd_siz = nbnd 334 bec%comm = mp_get_comm_null() 335 bec%nbnd = nbnd 336 bec%mype = 0 337 bec%nproc = 1 338 bec%nbnd_loc = nbnd 339 bec%ibnd_begin = 1 340 ! 341 IF( PRESENT( comm ) .AND. gamma_only .AND. smallmem ) THEN 342 bec%comm = comm 343 bec%nproc = mp_size( comm ) 344 IF( bec%nproc > 1 ) THEN 345 nbnd_siz = nbnd / bec%nproc 346 IF( MOD( nbnd, bec%nproc ) /= 0 ) nbnd_siz = nbnd_siz + 1 347 bec%mype = mp_rank( bec%comm ) 348 bec%nbnd_loc = ldim_block( becp%nbnd , bec%nproc, bec%mype ) 349 bec%ibnd_begin = gind_block( 1, becp%nbnd, bec%nproc, bec%mype ) 350 END IF 351 END IF 352 ! 353 IF ( gamma_only ) THEN 354 ! 355 ALLOCATE( bec%r( nkb, nbnd_siz ), STAT=ierr ) 356 IF( ierr /= 0 ) & 357 CALL errore( ' allocate_bec_type ', ' cannot allocate bec%r ', ABS(ierr) ) 358 ! 359 bec%r(:,:)=0.0D0 360 ! 361 ELSEIF ( noncolin) THEN 362 ! 363 ALLOCATE( bec%nc( nkb, npol, nbnd_siz ), STAT=ierr ) 364 IF( ierr /= 0 ) & 365 CALL errore( ' allocate_bec_type ', ' cannot allocate bec%nc ', ABS(ierr) ) 366 ! 367 bec%nc(:,:,:)=(0.0D0,0.0D0) 368 ! 369 ELSE 370 ! 371 ALLOCATE( bec%k( nkb, nbnd_siz ), STAT=ierr ) 372 IF( ierr /= 0 ) & 373 CALL errore( ' allocate_bec_type ', ' cannot allocate bec%k ', ABS(ierr) ) 374 ! 375 bec%k(:,:)=(0.0D0,0.0D0) 376 ! 377 ENDIF 378 ! 379 RETURN 380 ! 381 END SUBROUTINE allocate_bec_type 382 ! 383 !----------------------------------------------------------------------- 384 SUBROUTINE deallocate_bec_type (bec) 385 !----------------------------------------------------------------------- 386 ! 387 USE mp, ONLY: mp_get_comm_null 388 IMPLICIT NONE 389 TYPE (bec_type) :: bec 390 ! 391 bec%comm = mp_get_comm_null() 392 bec%nbnd = 0 393 ! 394 IF (allocated(bec%r)) DEALLOCATE(bec%r) 395 IF (allocated(bec%nc)) DEALLOCATE(bec%nc) 396 IF (allocated(bec%k)) DEALLOCATE(bec%k) 397 ! 398 RETURN 399 ! 400 END SUBROUTINE deallocate_bec_type 401 402 SUBROUTINE beccopy(bec, bec1, nkb, nbnd, comm) 403 USE mp, ONLY: mp_size, mp_sum 404 IMPLICIT NONE 405 TYPE(bec_type), INTENT(in) :: bec 406 TYPE(bec_type) :: bec1 407 INTEGER, INTENT(in) :: nkb, nbnd 408 INTEGER, INTENT (in), OPTIONAL :: comm 409 410 INTEGER :: nbgrp, ib_start, ib_end, this_bgrp_nbnd 411 412 nbgrp = 1; ib_start = 1; ib_end = nbnd ; this_bgrp_nbnd = nbnd 413 IF( PRESENT( comm ) ) THEN 414 nbgrp = mp_size( comm ) 415 call divide( comm, nbnd, ib_start, ib_end) ; this_bgrp_nbnd = ib_end - ib_start + 1 416 END IF 417 418 IF (gamma_only) THEN 419 if(nbgrp>1) bec1%r = 0.d0 420 CALL dcopy(nkb*this_bgrp_nbnd, bec%r, 1, bec1%r(1,ib_start), 1) 421 if (nbgrp > 1) CALL mp_sum( bec1%r, comm ) 422 ELSEIF (noncolin) THEN 423 if(nbgrp>1) bec1%nc = ( 0.d0, 0.d0 ) 424 CALL zcopy(nkb*npol*this_bgrp_nbnd, bec%nc, 1, bec1%nc(1,1,ib_start), 1) 425 if (nbgrp > 1) CALL mp_sum( bec1%nc, comm ) 426 ELSE 427 if(nbgrp>1) bec1%k = ( 0.d0, 0.d0 ) 428 CALL zcopy(nkb*this_bgrp_nbnd, bec%k, 1, bec1%k(1,ib_start), 1) 429 if (nbgrp > 1) CALL mp_sum( bec1%k, comm ) 430 ENDIF 431 432 RETURN 433 END SUBROUTINE beccopy 434 435 SUBROUTINE becscal_nck(alpha, bec, nkb, nbnd) 436 IMPLICIT NONE 437 TYPE(bec_type), INTENT(INOUT) :: bec 438 COMPLEX(DP), INTENT(IN) :: alpha 439 INTEGER, INTENT(IN) :: nkb, nbnd 440 441 IF (gamma_only) THEN 442 CALL errore('becscal_nck','called in the wrong case',1) 443 ELSEIF (noncolin) THEN 444 CALL zscal(nkb*npol*nbnd, alpha, bec%nc, 1) 445 ELSE 446 CALL zscal(nkb*nbnd, alpha, bec%k, 1) 447 ENDIF 448 449 RETURN 450 END SUBROUTINE becscal_nck 451 452 SUBROUTINE becscal_gamma(alpha, bec, nkb, nbnd) 453 IMPLICIT NONE 454 TYPE(bec_type), INTENT(INOUT) :: bec 455 REAL(DP), INTENT(IN) :: alpha 456 INTEGER, INTENT(IN) :: nkb, nbnd 457 458 IF (gamma_only) THEN 459 CALL dscal(nkb*nbnd, alpha, bec%r, 1) 460 ELSE 461 CALL errore('becscal_gamma','called in the wrong case',1) 462 ENDIF 463 464 RETURN 465 END SUBROUTINE becscal_gamma 466 467END MODULE becmod 468