1! 2! Copyright (C) 2002-2008 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 10 11!=----------------------------------------------------------------------------=! 12 SUBROUTINE interpolate_lambda_x( lambdap, lambda, lambdam ) 13!=----------------------------------------------------------------------------=! 14 USE kinds, ONLY: DP 15 IMPLICIT NONE 16 REAL(DP) :: lambdap(:,:,:), lambda(:,:,:), lambdam(:,:,:) 17 ! 18 ! interpolate new lambda at (t+dt) from lambda(t) and lambda(t-dt): 19 ! 20 lambdap= 2.d0*lambda - lambdam 21 lambdam=lambda 22 lambda =lambdap 23 RETURN 24 END SUBROUTINE interpolate_lambda_x 25 26 27!=----------------------------------------------------------------------------=! 28 SUBROUTINE update_lambda_x( i, lambda, c0, c2, n, noff, tdist ) 29!=----------------------------------------------------------------------------=! 30 USE kinds, ONLY: DP 31 USE electrons_module, ONLY: ib_owner, ib_local 32 USE mp_global, ONLY: me_bgrp, intra_bgrp_comm 33 USE mp, ONLY: mp_sum 34 USE wave_base, ONLY: hpsi 35 USE gvect, ONLY: gstart 36 IMPLICIT NONE 37 INTEGER, INTENT(IN) :: n, noff 38 REAL(DP) :: lambda(:,:) 39 COMPLEX(DP) :: c0(:,:), c2(:) 40 INTEGER, INTENT(IN) :: i 41 LOGICAL, INTENT(IN) :: tdist ! if .true. lambda is distributed 42 ! 43 REAL(DP), ALLOCATABLE :: prod(:) 44 INTEGER :: ibl 45 LOGICAL :: gzero 46 ! 47 gzero = (gstart == 2) 48 ALLOCATE( prod( n ) ) 49 prod = hpsi( gzero, c0, SIZE( c0, 1 ), c2, n, noff ) 50 CALL mp_sum( prod, intra_bgrp_comm ) 51 IF( tdist ) THEN 52 IF( me_bgrp == ib_owner( i ) ) THEN 53 ibl = ib_local( i ) 54 lambda( ibl, : ) = prod( : ) 55 END IF 56 ELSE 57 lambda( i, : ) = prod( : ) 58 END IF 59 DEALLOCATE( prod ) 60 RETURN 61 END SUBROUTINE update_lambda_x 62 63 64 65 66!=----------------------------------------------------------------------------=! 67 subroutine elec_fakekine_x( ekincm, ema0bg, emass, c0, cm, ngw, n, noff, delt ) 68!=----------------------------------------------------------------------------=! 69 ! 70 ! This subroutine computes the CP(fake) wave functions kinetic energy 71 72 USE kinds, only : DP 73 use mp, only : mp_sum 74 use mp_global, only : intra_bgrp_comm, nbgrp, inter_bgrp_comm 75 use gvect, only : gstart 76 use wave_base, only : wave_speed2 77 use mp_world, only: mpime 78 ! 79 IMPLICIT NONE 80 ! 81 integer, intent(in) :: ngw ! number of plane wave coeff. 82 integer, intent(in) :: n ! number of bands 83 integer, intent(in) :: noff ! offset for band index 84 real(DP), intent(out) :: ekincm 85 real(DP), intent(in) :: ema0bg( ngw ), delt, emass 86 complex(DP), intent(in) :: c0( ngw, n ), cm( ngw, n ) 87 ! 88 real(DP), allocatable :: emainv(:) 89 real(DP) :: ftmp 90 integer :: i 91 92 ekincm=0.0d0 93 94 IF( ngw > 0 ) THEN 95 96 ALLOCATE( emainv( ngw ) ) 97 emainv = 1.0d0 / ema0bg 98 ftmp = 1.0d0 99 if( gstart == 2 ) ftmp = 0.5d0 100 101 do i = noff, n + noff - 1 102 ekincm = ekincm + 2.0d0 * wave_speed2( c0(:,i), cm(:,i), emainv, ftmp ) 103 end do 104 ekincm = ekincm * emass / ( delt * delt ) 105 DEALLOCATE( emainv ) 106 107 END IF 108 109 CALL mp_sum( ekincm, intra_bgrp_comm ) 110 IF( nbgrp > 1 ) & 111 CALL mp_sum( ekincm, inter_bgrp_comm ) 112 113 return 114 end subroutine elec_fakekine_x 115 116!=----------------------------------------------------------------------------=! 117 subroutine bandsum( bsum, c0, ngw, tbgrp ) 118!=----------------------------------------------------------------------------=! 119 ! 120 ! This subroutine computes the CP(fake) wave functions kinetic energy 121 122 USE kinds, only : DP 123 use mp, only : mp_sum 124 use mp_global, only : intra_bgrp_comm, nbgrp, inter_bgrp_comm 125 USE electrons_base, ONLY : nbsp, nbsp_bgrp 126 ! 127 IMPLICIT NONE 128 ! 129 integer, intent(in) :: ngw ! number of plane wave coeff. 130 real(DP), intent(out) :: bsum 131 complex(DP), intent(in) :: c0( ngw, * ) 132 logical, intent(in) :: tbgrp 133 ! 134 integer :: i, n 135 136 n = nbsp 137 IF( tbgrp ) n = nbsp_bgrp 138 139 bsum=0.0d0 140 do i = 1, n 141 bsum = bsum + SUM( DBLE( CONJG( c0( :, i ) ) * c0( :, i ) ) ) 142 end do 143 CALL mp_sum( bsum, intra_bgrp_comm ) 144 IF( tbgrp ) & 145 CALL mp_sum( bsum, inter_bgrp_comm ) 146 147 return 148 end subroutine bandsum 149 150 151 152!=----------------------------------------------------------------------------=! 153 SUBROUTINE protate_x ( c0, bec, c0rot, becrot, ngwl, nss, noff, lambda, nrl, & 154 ityp, nat, indv_ijkb0, nh, np_rot, me_rot, comm_rot ) 155!=----------------------------------------------------------------------------=! 156 157 ! this routine rotates the wave functions using the matrix lambda 158 ! it works with a block-like distributed matrix 159 ! of the Lagrange multipliers ( lambda ). 160 ! no replicated data are used, allowing scalability for large problems. 161 ! the layout of lambda is as follows : 162 ! 163 ! (PE 0) (PE 1) .. (PE NPE-1) 164 ! lambda(1 ,1:nx) lambda(2 ,1:nx) .. lambda(NPE ,1:nx) 165 ! lambda(1+ NPE,1:nx) lambda(2+ NPE,1:nx) .. lambda(NPE+ NPE,1:nx) 166 ! lambda(1+2*NPE,1:nx) lambda(2+2*NPE,1:nx) .. lambda(NPE+2*NPE,1:nx) 167 ! 168 ! distributes lambda's rows across processors with a blocking factor 169 ! of 1, ( row 1 to PE 1, row 2 to PE 2, .. row nproc_bgrp+1 to PE 1 and 170 ! so on). 171 ! nrl = local number of rows 172 ! ---------------------------------------------- 173 174 ! ... declare modules 175 176 USE kinds, ONLY: DP 177 USE mp, ONLY: mp_bcast 178 USE mp_global, ONLY: nproc_bgrp, me_bgrp, intra_bgrp_comm 179 180 IMPLICIT NONE 181 182 ! ... declare subroutine arguments 183 184 INTEGER, INTENT(IN) :: ngwl, nss, nrl, noff 185 INTEGER, INTENT(IN) :: ityp(:), nat, indv_ijkb0(:), nh(:) 186 INTEGER, INTENT(IN) :: np_rot, me_rot, comm_rot 187 COMPLEX(DP), INTENT(IN) :: c0(:,:) 188 COMPLEX(DP), INTENT(OUT) :: c0rot(:,:) 189 REAL(DP), INTENT(IN) :: lambda(:,:) 190 REAL(DP), INTENT(IN) :: bec(:,:) 191 REAL(DP), INTENT(OUT) :: becrot(:,:) 192 193 ! ... declare other variables 194 INTEGER :: i, j, k, ip 195 INTEGER :: jl, nrl_ip, is, ia, jv, jnl, nj 196 REAL(DP), ALLOCATABLE :: uu(:,:) 197 198 IF( nss < 1 ) THEN 199 RETURN 200 END IF 201 202 CALL start_clock('protate') 203 204 DO i = 1, nss 205 c0rot( :, i+noff-1 ) = 0.0d0 206 becrot(:,i+noff-1 ) = 0.0d0 207 END DO 208 209 210 211! becrot = 0.0d0 212! c0rot = 0.0d0 213 214 DO ip = 1, np_rot 215 216 nrl_ip = nss / np_rot 217 IF( (ip-1) < mod( nss, np_rot ) ) THEN 218 nrl_ip = nrl_ip + 1 219 END IF 220 221 ALLOCATE( uu( nrl_ip, nss ) ) 222 IF( me_rot .EQ. (ip-1) ) THEN 223 uu = lambda( 1:nrl_ip, 1:nss ) 224 END IF 225 CALL mp_bcast( uu, (ip-1), intra_bgrp_comm) 226 227 j = ip 228 DO jl = 1, nrl_ip 229 DO i = 1, nss 230 CALL daxpy(2*ngwl,uu(jl,i),c0(1,j+noff-1),1,c0rot(1,i+noff-1),1) 231 END DO 232 233 do ia=1,nat 234 is=ityp(ia) 235 do jv=1,nh(is) 236 jnl = indv_ijkb0(ia) + jv 237 do i = 1, nss 238 becrot(jnl,i+noff-1) = becrot(jnl,i+noff-1)+ uu(jl, i) * bec( jnl, j+noff-1 ) 239 end do 240 end do 241 end do 242 243 j = j + np_rot 244 END DO 245 246 DEALLOCATE(uu) 247 248 END DO 249 250 CALL stop_clock('protate') 251 252 RETURN 253 END SUBROUTINE protate_x 254 255 256 257!=----------------------------------------------------------------------------=! 258 SUBROUTINE crot_gamma2 ( c0rot, c0, ngw, n, noffr, noff, lambda, nx, eig ) 259!=----------------------------------------------------------------------------=! 260 261 ! this routine rotates the wave functions to the Kohn-Sham base 262 ! it works with a block-like distributed matrix 263 ! of the Lagrange multipliers ( lambda ). 264 ! 265 ! ... declare modules 266 267 USE kinds, ONLY: DP 268 USE mp, ONLY: mp_bcast 269 USE mp_global, ONLY: nproc_bgrp, me_bgrp, intra_bgrp_comm 270 271 IMPLICIT NONE 272 273 include 'laxlib.fh' 274 275 ! ... declare subroutine arguments 276 277 INTEGER, INTENT(IN) :: ngw, n, nx, noffr, noff 278 COMPLEX(DP), INTENT(INOUT) :: c0rot(:,:) 279 COMPLEX(DP), INTENT(IN) :: c0(:,:) 280 REAL(DP), INTENT(IN) :: lambda(:,:) 281 REAL(DP), INTENT(OUT) :: eig(:) 282 283 ! ... declare other variables 284 ! 285 REAL(DP), ALLOCATABLE :: vv(:,:), ap(:) 286 INTEGER :: i, j, k 287 288 IF( nx < 1 ) THEN 289 RETURN 290 END IF 291 292 ALLOCATE( vv( nx, nx ) ) 293 294 ! NON distributed lambda 295 296 ALLOCATE( ap( nx * ( nx + 1 ) / 2 ) ) 297 298 K = 0 299 DO J = 1, n 300 DO I = J, n 301 K = K + 1 302 ap( k ) = lambda( i, j ) 303 END DO 304 END DO 305 306 CALL dspev_drv( 'V', 'L', n, ap, eig, vv, nx ) 307 308 DEALLOCATE( ap ) 309 310 DO i = 1, n 311 c0rot( :, i+noffr-1 ) = 0.0d0 312 END DO 313 314 DO j = 1, n 315 DO i = 1, n 316 CALL daxpy( 2*ngw, vv(j,i), c0(1,j+noff-1), 1, c0rot(1,i+noffr-1), 1 ) 317 END DO 318 END DO 319 320 DEALLOCATE( vv ) 321 322 RETURN 323 END SUBROUTINE crot_gamma2 324 325 326 327!=----------------------------------------------------------------------------=! 328 SUBROUTINE proj_gamma( a, b, ngw, n, noff, lambda) 329!=----------------------------------------------------------------------------=! 330 331 ! projection A=A-SUM{B}<B|A>B 332 ! no replicated data are used, allowing scalability for large problems. 333 ! The layout of lambda is as follows : 334 ! 335 ! (PE 0) (PE 1) .. (PE NPE-1) 336 ! lambda(1 ,1:nx) lambda(2 ,1:nx) .. lambda(NPE ,1:nx) 337 ! lambda(1+ NPE,1:nx) lambda(2+ NPE,1:nx) .. lambda(NPE+ NPE,1:nx) 338 ! lambda(1+2*NPE,1:nx) lambda(2+2*NPE,1:nx) .. lambda(NPE+2*NPE,1:nx) 339 ! 340 ! distribute lambda's rows across processors with a blocking factor 341 ! of 1, ( row 1 to PE 1, row 2 to PE 2, .. row nproc_bgrp+1 to PE 1 and so on). 342 ! ---------------------------------------------- 343 344! ... declare modules 345 USE kinds, ONLY: DP 346 USE mp_global, ONLY: nproc_bgrp, me_bgrp, intra_bgrp_comm 347 USE wave_base, ONLY: dotp 348 USE gvect, ONLY: gstart 349 350 IMPLICIT NONE 351 352! ... declare subroutine arguments 353 INTEGER, INTENT( IN ) :: ngw, n, noff 354 COMPLEX(DP), INTENT(INOUT) :: a(:,:), b(:,:) 355 REAL(DP), OPTIONAL :: lambda(:,:) 356 357! ... declare other variables 358 REAL(DP), ALLOCATABLE :: ee(:) 359 INTEGER :: i, j, jl 360 COMPLEX(DP) :: alp 361 LOGICAL :: gzero 362 363! ... end of declarations 364! ---------------------------------------------- 365 366 IF( n < 1 ) THEN 367 RETURN 368 END IF 369 gzero = (gstart == 2) 370 ALLOCATE( ee( n ) ) 371 372 DO i = 1, n 373 DO j = 1, n 374 ee(j) = -dotp( gzero, ngw, b(:,j+noff-1), a(:,i+noff-1), intra_bgrp_comm ) 375 END DO 376 IF( PRESENT(lambda) ) THEN 377 IF( MOD( (i-1), nproc_bgrp ) == me_bgrp ) THEN 378 DO j = 1, n 379 lambda( (i-1) / nproc_bgrp + 1, j ) = ee(j) 380 END DO 381 END IF 382 END IF 383 DO j = 1, n 384 alp = CMPLX(ee(j),0.0d0,kind=DP) 385 CALL zaxpy( ngw, alp, b(1,j+noff-1), 1, a(1,i+noff-1), 1 ) 386 END DO 387 END DO 388 DEALLOCATE(ee) 389 390 RETURN 391 END SUBROUTINE proj_gamma 392 393 394 395 396!=----------------------------------------------------------------------------=! 397 SUBROUTINE wave_rand_init_x( cm_bgrp, global ) 398!=----------------------------------------------------------------------------=! 399 400 ! this routine sets the initial wavefunctions at random 401 402! ... declare modules 403 USE kinds, ONLY: DP 404 USE mp, ONLY: mp_sum, mp_max, mp_min 405 USE mp_wave, ONLY: splitwf 406 USE mp_global, ONLY: me_bgrp, nproc_bgrp, root_bgrp, intra_bgrp_comm 407 USE gvect, ONLY: ig_l2g, gstart 408 USE gvect, ONLY: mill, gg 409 USE gvecw, ONLY: ngw, ngw_g 410 USE io_global, ONLY: stdout 411 USE random_numbers, ONLY: randy 412 USE electrons_base, ONLY: nbsp, ibgrp_g2l 413 414 IMPLICIT NONE 415 416 ! ... declare subroutine arguments 417 COMPLEX(DP), INTENT(OUT) :: cm_bgrp(:,:) 418 LOGICAL, OPTIONAL, INTENT(IN) :: global 419 420 ! ... declare other variables 421 INTEGER :: ntest, ig, ib, ibgrp, lb, ub 422 REAL(DP) :: rranf1, rranf2, ampre, ggx, fac, r1, r2, r3 423 COMPLEX(DP), ALLOCATABLE :: pwt( : ) 424 REAL(DP), ALLOCATABLE :: RND( : , : ) 425 INTEGER :: iss, n1, n2, m1, m2, i 426 LOGICAL :: local 427 428 ! ... initialize the wave functions in such a way that the values 429 ! ... of the components are independent on the number of processors 430 ! ... with "local" option the initialization is independend from G sorting too! 431 432 ! ... Check array dimensions 433 434 IF( SIZE( cm_bgrp, 1 ) < ngw ) THEN 435 CALL errore(' wave_rand_init ', ' wrong dimensions ', 3) 436 END IF 437 438 local = .TRUE. 439 IF( PRESENT( global ) ) THEN 440 local = .NOT. global 441 END IF 442 443 ! ... Reset them to zero 444 445 cm_bgrp = 0.0d0 446 447 ampre = 0.01d0 448 449 IF( local ) THEN 450 ggx = MAXVAL( gg( 1:ngw ) ) 451 CALL mp_max( ggx, intra_bgrp_comm ) 452 lb = MINVAL( mill ) 453 CALL mp_min( lb, intra_bgrp_comm ) 454 ub = MAXVAL( mill ) 455 CALL mp_max( ub, intra_bgrp_comm ) 456 ALLOCATE( RND( 3, lb:ub ) ) 457 ELSE 458 ALLOCATE( pwt( ngw_g ) ) 459 ntest = ngw_g / 4 460 IF( ntest < SIZE( cm_bgrp, 2 ) ) THEN 461 ntest = ngw_g 462 END IF 463 END IF 464 ! 465 ! ... assign random values to wave functions 466 ! 467 DO ib = 1, nbsp 468 469 IF( local ) THEN 470 rnd = 0.0d0 471 DO ig = lb, ub 472 rnd( 1, ig ) = 0.5d0 - randy() 473 rnd( 2, ig ) = 0.5d0 - randy() 474 rnd( 3, ig ) = 0.5d0 - randy() 475 END DO 476 ELSE 477 pwt( : ) = 0.0d0 478 DO ig = 3, ntest 479 rranf1 = 0.5d0 - randy() 480 rranf2 = randy() 481 pwt( ig ) = ampre * CMPLX(rranf1, rranf2,kind=DP) 482 END DO 483 END IF 484 ! 485 ibgrp = ibgrp_g2l( ib ) 486 ! 487 IF( ibgrp > 0 ) THEN 488 DO ig = 1, ngw 489 IF( local ) THEN 490 IF( gg(ig) < ggx / 2.519d0 ) THEN ! 2.519 = 4^(2/3), equivalent to keep only (ngw_g/4) values 491 rranf1 = rnd( 1, mill(1,ig) ) * rnd( 2, mill(2,ig) ) * rnd( 3, mill(3,ig) ) 492 rranf2 = 0.0d0 493 cm_bgrp( ig, ibgrp ) = ampre * CMPLX( rranf1, rranf2 ,kind=DP) / ( 1.0d0 + gg(ig) ) 494 END IF 495 ELSE 496 cm_bgrp( ig, ibgrp ) = pwt( ig_l2g( ig ) ) 497 END IF 498 END DO 499 END IF 500 ! 501 END DO 502 IF ( gstart == 2 ) THEN 503 cm_bgrp( 1, : ) = (0.0d0, 0.0d0) 504 END IF 505 506 IF( ALLOCATED( pwt ) ) DEALLOCATE( pwt ) 507 IF( ALLOCATED( rnd ) ) DEALLOCATE( rnd ) 508 509#ifdef PIPPO_DEBUG 510 write(1000+mpime,fmt='(8I5)') nbsp, nbsp_bgrp, nudx, nudx_bgrp, nbsp, nbsp_bgrp, nbspx, nbspx_bgrp 511 DO iss = 1, nspin 512 write(1000+mpime,fmt='(5I5)') nupdwn(iss), iupdwn(iss), nupdwn_bgrp(iss), iupdwn_bgrp(iss), i2gupdwn_bgrp(iss) 513 END DO 514 DO ib = 1, nbsp 515 ! write(1000+mpime,fmt='(2I5)') ib, ibgrp_g2l(ib) 516 END DO 517 518 DO iss = nspin, 1, -1 519 write(1000+mpime,*) 'copy' 520 n1 = iupdwn_bgrp(iss) 521 n2 = n1 + nupdwn_bgrp(iss) - 1 522 m1 = iupdwn(iss)+i2gupdwn_bgrp(iss) - 1 523 m2 = m1 + nupdwn_bgrp(iss) - 1 524 DO i = m2, m1, -1 525 write(1000+mpime,fmt='(2I5)') i, i-m1+n1 526 END DO 527 END DO 528 DO iss = 1, nspin 529 m1 = iupdwn(iss)+i2gupdwn_bgrp(iss) - 1 530 m2 = m1 + nupdwn_bgrp(iss) - 1 531 write(1000+mpime,*) 'zero' 532 DO i = iupdwn(iss), m1-1 533 write(1000+mpime,fmt='(1I5)') i 534 END DO 535 write(1000+mpime,*) 'zero' 536 DO i = m2+1, iupdwn(iss) + nupdwn(iss) - 1 537 write(1000+mpime,fmt='(1I5)') i 538 END DO 539 END DO 540#endif 541 542 543 RETURN 544 END SUBROUTINE wave_rand_init_x 545 546 547 SUBROUTINE c_bgrp_expand_x( c_bgrp ) 548 USE kinds, ONLY: DP 549 USE mp, ONLY: mp_sum 550 USE electrons_base, ONLY: nspin, i2gupdwn_bgrp, nupdwn, iupdwn_bgrp, iupdwn, nupdwn_bgrp 551 USE mp_global, ONLY: nbgrp, inter_bgrp_comm 552 IMPLICIT NONE 553 COMPLEX(DP) :: c_bgrp(:,:) 554 INTEGER :: iss, n1, n2, m1, m2, i 555 IF( nbgrp < 2 ) & 556 RETURN 557 DO iss = nspin, 1, -1 558 n1 = iupdwn_bgrp(iss) 559 n2 = n1 + nupdwn_bgrp(iss) - 1 560 m1 = iupdwn(iss)+i2gupdwn_bgrp(iss) - 1 561 m2 = m1 + nupdwn_bgrp(iss) - 1 562 DO i = m2, m1, -1 563 c_bgrp(:,i) = c_bgrp(:,i-m1+n1) 564 END DO 565 END DO 566 DO iss = 1, nspin 567 m1 = iupdwn(iss)+i2gupdwn_bgrp(iss) - 1 568 m2 = m1 + nupdwn_bgrp(iss) - 1 569 DO i = iupdwn(iss), m1-1 570 c_bgrp(:,i) = 0.0d0 571 END DO 572 DO i = m2+1, iupdwn(iss) + nupdwn(iss) - 1 573 c_bgrp(:,i) = 0.0d0 574 END DO 575 END DO 576 CALL mp_sum( c_bgrp, inter_bgrp_comm ) 577 RETURN 578 END SUBROUTINE c_bgrp_expand_x 579 580 581 SUBROUTINE c_bgrp_pack_x( c_bgrp ) 582 USE kinds, ONLY: DP 583 USE electrons_base, ONLY: nspin, i2gupdwn_bgrp, nupdwn, iupdwn_bgrp, iupdwn, nupdwn_bgrp 584 USE mp_global, ONLY: nbgrp 585 IMPLICIT NONE 586 COMPLEX(DP) :: c_bgrp(:,:) 587 INTEGER :: iss, n1, n2, m1, m2, i 588 IF( nbgrp < 2 ) & 589 RETURN 590 DO iss = 1, nspin 591 n1 = iupdwn_bgrp(iss) 592 n2 = n1 + nupdwn_bgrp(iss) - 1 593 m1 = iupdwn(iss)+i2gupdwn_bgrp(iss) - 1 594 m2 = m1 + nupdwn_bgrp(iss) - 1 595 DO i = n1, n2 596 c_bgrp(:,i) = c_bgrp(:,i-n1+m1) 597 END DO 598 END DO 599 RETURN 600 END SUBROUTINE c_bgrp_pack_x 601 602