1! 2! Copyright (C) 2001-2020 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!----------------------------------------------------------------------- 10SUBROUTINE vhpsi( ldap, np, mps, psip, hpsi ) 11 !----------------------------------------------------------------------- 12 !! This routine computes the Hubbard potential applied to the electronic 13 !! structure of the current k-point. The result is added to hpsi. 14 ! 15 USE kinds, ONLY : DP 16 USE becmod, ONLY : bec_type, calbec, allocate_bec_type, & 17 deallocate_bec_type 18 USE ldaU, ONLY : Hubbard_lmax, Hubbard_l, is_Hubbard, & 19 nwfcU, wfcU, offsetU, lda_plus_u_kind, & 20 is_hubbard_back, Hubbard_l_back, offsetU_back, & 21 backall, offsetU_back1 22 USE lsda_mod, ONLY : current_spin 23 USE scf, ONLY : v 24 USE ions_base, ONLY : nat, ntyp => nsp, ityp 25 USE control_flags, ONLY : gamma_only 26 USE mp, ONLY : mp_sum 27 ! 28 IMPLICIT NONE 29 ! 30 INTEGER, INTENT(IN) :: ldap 31 !! leading dimension of arrays psip, hpsi 32 INTEGER, INTENT(IN) :: np 33 !! true dimension of psip, hpsi 34 INTEGER, INTENT(IN) :: mps 35 !! number of states psip 36 COMPLEX(DP), INTENT(IN) :: psip(ldap,mps) 37 !! the wavefunction 38 COMPLEX(DP), INTENT(INOUT) :: hpsi(ldap,mps) 39 !! Hamiltonian dot psi 40 ! 41 ! ... local variables 42 ! 43 REAL(DP), ALLOCATABLE :: rtemp(:,:) 44 COMPLEX(DP), ALLOCATABLE :: ctemp(:,:), vaux(:,:) 45 TYPE(bec_type) :: proj 46 ! 47 CALL start_clock('vhpsi') 48 ! 49 ! Offset of atomic wavefunctions initialized in setup and stored in offsetU 50 ! 51 ! Allocate the array proj 52 CALL allocate_bec_type ( nwfcU,mps, proj ) 53 ! 54 ! proj = <wfcU|psip> 55 CALL calbec (np, wfcU, psip, proj) 56 ! 57 IF ( lda_plus_u_kind.EQ.0 .OR. lda_plus_u_kind.EQ.1 ) THEN 58 CALL vhpsi_U () ! DFT+U 59 ELSEIF ( lda_plus_u_kind.EQ.2 ) THEN 60 CALL vhpsi_UV () ! DFT+U+V 61 ENDIF 62 ! 63 CALL deallocate_bec_type (proj) 64 ! 65 CALL stop_clock('vhpsi') 66 ! 67 RETURN 68 ! 69CONTAINS 70 ! 71SUBROUTINE vhpsi_U () 72 ! 73 ! This routine applies the Hubbard potential with U_I 74 ! to the KS wave functions. 75 ! 76 USE ldaU, ONLY : ldim_back, ldmx_b, Hubbard_l1_back 77 ! 78 IMPLICIT NONE 79 INTEGER :: na, nt, ldim, ldim0 80 ! 81 DO nt = 1, ntyp 82 ! 83 ! Compute the action of the Hubbard potential on the KS wave functions: 84 ! V_Hub |psip > = \sum v%ns |wfcU> <wfcU|psip> 85 ! where v%ns = U ( delta/2 - rho%ns ) is computed in v_of_rho 86 ! 87 IF ( is_hubbard(nt) ) THEN 88 ! 89 ldim = 2*Hubbard_l(nt) + 1 90 ! 91 IF (gamma_only) THEN 92 ALLOCATE ( rtemp(ldim,mps) ) 93 ELSE 94 ALLOCATE ( ctemp(ldim,mps) ) 95 ENDIF 96 ! 97 DO na = 1, nat 98 IF ( nt == ityp(na) ) THEN 99 IF (gamma_only) THEN 100 CALL DGEMM ('n','n', ldim,mps,ldim, 1.0_dp, & 101 v%ns(1,1,current_spin,na),2*Hubbard_lmax+1, & 102 proj%r(offsetU(na)+1,1),nwfcU, 0.0_dp, rtemp, ldim) 103 CALL DGEMM ('n','n', 2*np, mps, ldim, 1.0_dp, & 104 wfcU(1,offsetU(na)+1), 2*ldap, rtemp, ldim, & 105 1.0_dp, hpsi, 2*ldap) 106 ELSE 107 ! 108 ALLOCATE(vaux(ldim,ldim)) 109 ! 110 vaux = (0.0_dp,0.0_dp) 111 vaux(:,:) = v%ns(:,:,current_spin,na) 112 ! 113 CALL ZGEMM ('n','n', ldim, mps, ldim, (1.0_dp,0.0_dp), & 114 vaux, ldim, proj%k(offsetU(na)+1,1), nwfcU, & 115 (0.0_dp,0.0_dp), ctemp, ldim) 116 ! 117 DEALLOCATE(vaux) 118 ! 119 CALL ZGEMM ('n','n', np, mps, ldim, (1.0_dp,0.0_dp), & 120 wfcU(1,offsetU(na)+1), ldap, ctemp, ldim, & 121 (1.0_dp,0.0_dp), hpsi, ldap) 122 ! 123 ENDIF 124 ENDIF 125 ENDDO 126 ! 127 IF (gamma_only) THEN 128 DEALLOCATE ( rtemp ) 129 ELSE 130 DEALLOCATE ( ctemp ) 131 ENDIF 132 ! 133 ENDIF 134 ! 135 ! If the background is used then compute extra 136 ! contribution to the Hubbard potential 137 ! 138 IF ( is_hubbard_back(nt) ) THEN 139 ! 140 ldim = ldim_back(nt) 141 ! 142 IF (gamma_only) THEN 143 ALLOCATE ( rtemp(ldim,mps) ) 144 ELSE 145 ALLOCATE ( ctemp(ldim,mps) ) 146 ENDIF 147 ! 148 DO na = 1, nat 149 IF ( nt == ityp(na) ) THEN 150 ! 151 IF (gamma_only) THEN 152 ! 153 ldim = 2*Hubbard_l_back(nt)+1 154 ! 155 CALL DGEMM ('n','n', ldim,mps,ldim, 1.0_dp, & 156 v%nsb(1,1,current_spin,na),ldmx_b, & 157 proj%r(offsetU_back(na)+1,1), & 158 nwfcU, 0.0_dp, rtemp, ldim_back(nt)) 159 ! 160 CALL DGEMM ('n','n', 2*np, mps, ldim, 1.0_dp, & 161 wfcU(1,offsetU_back(na)+1), 2*ldap, rtemp, & 162 ldim_back(nt), 1.0_dp, hpsi, 2*ldap) 163 ! 164 IF (backall(nt)) THEN 165 ! 166 ldim = 2*Hubbard_l1_back(nt)+1 167 ldim0 = 2*Hubbard_l_back(nt)+1 168 ! 169 CALL DGEMM ('n','n', ldim,mps,ldim, 1.0_dp, & 170 v%nsb(ldim0+1,ldim0+1,current_spin,na), & 171 ldim_back(nt), proj%r(offsetU_back1(na)+1,1), & 172 nwfcU, 0.0_dp, rtemp, ldim_back(nt)) 173 ! 174 CALL DGEMM ('n','n', 2*np, mps, ldim, 1.0_dp, & 175 wfcU(1,offsetU_back1(na)+1), 2*ldap, rtemp, & 176 ldim_back(nt), 1.0_dp, hpsi, 2*ldap) 177 ! 178 ENDIF 179 ! 180 ELSE 181 ! 182 ldim = ldim_back(nt) 183 ! 184 ALLOCATE(vaux(ldim,ldim)) 185 ! 186 vaux = (0.0_dp, 0.0_dp) 187 vaux(:,:) = v%nsb(:,:,current_spin,na) 188 ! 189 ldim = 2*Hubbard_l_back(nt)+1 190 ! 191 CALL ZGEMM ('n','n', ldim,mps,ldim, (1.0_dp,0.0_dp), & 192 vaux,ldim_back(nt), proj%k(offsetU_back(na)+1,1), & 193 nwfcU, (0.0_dp,0.0_dp), ctemp, ldim_back(nt)) 194 ! 195 CALL ZGEMM ('n','n', np, mps, ldim, (1.0_dp,0.0_dp), & 196 wfcU(1,offsetU_back(na)+1), ldap, ctemp, & 197 ldim_back(nt), (1.0_dp,0.0_dp), hpsi, ldap) 198 ! 199 IF (backall(nt)) THEN 200 ! 201 ldim = 2*Hubbard_l1_back(nt)+1 202 ldim0 = 2*Hubbard_l_back(nt)+1 203 ! 204 CALL ZGEMM ('n','n', ldim,mps,ldim,(1.0_dp,0.0_dp), & 205 vaux(ldim0+1,ldim0+1),ldim_back(nt), & 206 proj%k(offsetU_back1(na)+1,1), nwfcU, & 207 (0.0_dp,0.0_dp), ctemp, ldim_back(nt)) 208 ! 209 CALL ZGEMM ('n','n', np, mps, ldim, (1.0_dp,0.0_dp),& 210 wfcU(1,offsetU_back1(na)+1), ldap, ctemp, & 211 ldim_back(nt), (1.0_dp,0.0_dp), hpsi, ldap) 212 ! 213 ENDIF 214 ! 215 DEALLOCATE(vaux) 216 ! 217 ENDIF 218 ENDIF 219 ENDDO 220 ! 221 IF (gamma_only) THEN 222 DEALLOCATE ( rtemp ) 223 ELSE 224 DEALLOCATE ( ctemp ) 225 ENDIF 226 ! 227 ENDIF 228 ! 229 ENDDO 230 ! 231 RETURN 232 ! 233END SUBROUTINE vhpsi_U 234!-------------------------------------------------------------------------- 235 236!-------------------------------------------------------------------------- 237SUBROUTINE vhpsi_UV () 238 ! 239 ! This routine applies the Hubbard potential with U_I (=V_II) and V_IJ 240 ! to the KS wave functions. 241 ! By taking a derivative of the Hubbard energy we obtain the potential 242 ! (multiplied by the KS wave function psi_nk): 243 ! - \sum_IJ (J\=I) \sum_{m1,m2} V_IJ/2 244 ! * [ n^IJ_{m1,m2} * |phi^I_m1><phi^J_m2|psi_nk> + 245 ! n^JI_{m2,m1} * |phi^J_m2><phi^I_m1|psi_nk> ] 246 ! Using the symmetry with respect to I/J and m1/m2 we can simplify 247 ! the expression for the Hubbard potential above and write it as: 248 ! - \sum_IJ (J\=I) \sum_{m1,m2} V_IJ 249 ! * n^IJ_{m1,m2} * |phi^I_m1><phi^J_m2|psi_nk> 250 ! Since in practice the indices I and J are not really equivalent 251 ! (due to the way how the generalized Hubbard potential is implemeted), 252 ! instead of the second expression here the first expression is implemented. 253 ! 254 ! Note: This routine assumes that the phase factor phase_fac at a given k point 255 ! has been already computed elsewhere. 256 ! 257 USE ldaU, ONLY : ldim_u, neighood, at_sc, phase_fac, Hubbard_V, v_nsg 258 ! 259 IMPLICIT NONE 260 COMPLEX(DP) :: phase 261 INTEGER :: ldim2, ldimx, ldim1, m1, m2, equiv_na2, & 262 off1, off2, ig, viz, na1, na2, nt1, nt2 263 REAL(DP), ALLOCATABLE :: projauxr(:,:), rvaux(:,:) 264 COMPLEX(DP), ALLOCATABLE :: projauxc(:,:), wfcUaux(:,:) 265 ! 266 ! Find the maximum number of magnetic quantum numbers [i.e. MAX(2l+1)] 267 ! 268 ldimx = 0 269 DO nt1 = 1, ntyp 270 IF ( is_hubbard(nt1) .OR. is_hubbard_back(nt1) ) THEN 271 ldim1 = ldim_u(nt1) 272 ldimx = MAX(ldimx,ldim1) 273 ENDIF 274 ENDDO 275 ! 276 IF (gamma_only) THEN 277 ALLOCATE (rtemp(ldimx,mps)) 278 ALLOCATE (projauxr(ldimx,mps)) 279 ALLOCATE (rvaux(ldimx,ldimx)) 280 ELSE 281 ALLOCATE (ctemp(ldimx,mps)) 282 ALLOCATE (projauxc(ldimx,mps)) 283 ALLOCATE (vaux(ldimx,ldimx)) 284 ENDIF 285 ! 286 ALLOCATE (wfcUaux(np,ldimx)) 287 ! 288 DO nt1 = 1, ntyp 289 ! 290 ldim1 = ldim_u(nt1) 291 ! 292 IF ( is_hubbard(nt1) .OR. is_hubbard_back(nt1) ) THEN 293 ! 294 DO na1 = 1, nat 295 ! 296 IF (ityp(na1).EQ.nt1) THEN 297 ! 298 DO viz = 1, neighood(na1)%num_neigh 299 ! 300 na2 = neighood(na1)%neigh(viz) 301 equiv_na2 = at_sc(na2)%at 302 nt2 = ityp(equiv_na2) 303 phase = phase_fac(na2) 304 ldim2 = ldim_u(nt2) 305 ! 306 ! Note: Below there is a condition on v_nsg 307 ! because it may be that the user specifies 308 ! Hubbard_alpha for some type for which 309 ! Hubbard_V was not specified. 310 ! 311 IF ( (is_hubbard(nt2).OR.is_hubbard_back(nt2)) .AND. & 312 (Hubbard_V(na1,na2,1).NE.0.d0 .OR. & 313 Hubbard_V(na1,na2,2).NE.0.d0 .OR. & 314 Hubbard_V(na1,na2,3).NE.0.d0 .OR. & 315 Hubbard_V(na1,na2,4).NE.0.d0 .OR. & 316 ANY(v_nsg(:,:,viz,na1,current_spin).NE.0.0d0)) ) THEN 317 ! 318 ! Compute the first part of the Hubbard potential, namely: 319 ! - \sum_IJ (J\=I) \sum_{m1,m2} V_IJ/2 320 ! * n^IJ_{m1,m2} * |phi^I_m1><phi^J_m2|psi_nk> 321 ! where 322 ! - V_IJ/2 * n^IJ_{m1,m2} = CONJG(v_nsg) 323 ! <phi^J_m2|Psi_nk> = proj%r (or proj%k) 324 ! |phi^I_m1> = wfcU 325 ! 326 wfcUaux(:,:) = (0.0_dp, 0.0_dp) 327 ! 328 off1 = offsetU(na1) 329 ! 330 DO m1 = 1, ldim_u(nt1) 331 ! 332 IF (m1.GT.2*Hubbard_l(nt1)+1) & 333 off1 = offsetU_back(na1) - 2*Hubbard_l(nt1) - 1 334 ! 335 IF (backall(nt1) .AND. & 336 m1.GT.(2*Hubbard_l(nt1)+1+2*Hubbard_l_back(nt1)+1)) & 337 off1 = offsetU_back1(na1) & 338 - 2*Hubbard_l(nt1) - 2 - 2*Hubbard_l_back(nt1) 339 ! 340 DO ig = 1, np 341 wfcUaux(ig,m1) = wfcU(ig,off1+m1) 342 ENDDO 343 ! 344 ENDDO 345 ! 346 off2 = offsetU(equiv_na2) 347 ! 348 IF (gamma_only) THEN 349 ! 350 rvaux(:,:) = 0.0_dp 351 ! 352 projauxr(:,:) = 0.0_dp 353 ! 354 DO m1 = 1, ldim1 355 ! 356 DO m2 = 1, ldim2 357 ! 358 rvaux(m2,m1) = DBLE( (v_nsg(m2, m1, viz, na1, current_spin))) * 0.5d0 359 ! 360 ENDDO 361 ! 362 ENDDO 363 ! 364 DO m2 = 1, ldim2 365 ! 366 IF (m2.GT.2*Hubbard_l(nt2)+1) & 367 off2 = offsetU_back(equiv_na2) - 2*Hubbard_l(nt2) - 1 368 ! 369 IF (backall(nt2) .AND. & 370 m2.GT.(2*Hubbard_l(nt2)+1+2*Hubbard_l_back(nt2)+1)) & 371 off2 = offsetU_back1(equiv_na2) & 372 - 2*Hubbard_l(nt2) - 2 - 2*Hubbard_l_back(nt2) 373 ! 374 projauxr(m2,:) = DBLE(proj%r(off2+m2,:)) 375 ! 376 ENDDO 377 ! 378 rtemp(:,:) = 0.0_dp 379 ! 380 CALL DGEMM ('t','n', ldim1,mps,ldim2, 1.0_dp, & 381 rvaux,ldimx, projauxr,ldimx, 0.0_dp, rtemp, ldimx) 382 ! 383 CALL DGEMM ('n','n', 2*np, mps, ldim1, 1.0_dp, & 384 wfcUaux, 2*np, rtemp, ldimx, & 385 1.0_dp, hpsi, 2*ldap) 386 ! 387 ELSE 388 ! 389 vaux(:,:) = (0.0_dp, 0.0_dp) 390 ! 391 projauxc(:,:) = (0.0_dp, 0.0_dp) 392 ! 393 DO m1 = 1, ldim1 394 ! 395 DO m2 = 1, ldim2 396 ! 397 vaux(m2,m1) = CONJG( (v_nsg(m2, m1, viz, na1, current_spin))) * 0.5d0 398 ! 399 ENDDO 400 ! 401 ENDDO 402 ! 403 DO m2 = 1, ldim2 404 ! 405 IF (m2.GT.2*Hubbard_l(nt2)+1) & 406 off2 = offsetU_back(equiv_na2) - 2*Hubbard_l(nt2) - 1 407 ! 408 IF (backall(nt2) .AND. & 409 m2.GT.(2*Hubbard_l(nt2)+1+2*Hubbard_l_back(nt2)+1)) & 410 off2 = offsetU_back1(equiv_na2) & 411 - 2*Hubbard_l(nt2) - 2 - 2*Hubbard_l_back(nt2) 412 ! 413 projauxc(m2,:) = proj%k(off2+m2,:) 414 ! 415 ENDDO 416 ! 417 ctemp(:,:) = (0.0_dp,0.0_dp) 418 ! 419 CALL ZGEMM ('t','n', ldim1,mps,ldim2, (1.0_dp,0.0_dp), & 420 vaux,ldimx, projauxc,ldimx, (0.0_dp,0.0_dp), ctemp, ldimx) 421 ! 422 CALL ZGEMM ('n','n', np, mps, ldim1, phase, & 423 wfcUaux, np, ctemp, ldimx, (1.0_dp,0.0_dp), hpsi, ldap) 424 ! 425 ENDIF 426 ! 427 ! 428 ! Compute the second part of the Hubbard potential, namely: 429 ! - \sum_IJ (J\=I) \sum_m1m2 V_IJ/2 * n^JI_m2m1 * |phi^J_m2><phi^I_m1|Psi_nk> 430 ! where 431 ! - V_IJ/2 * n^JI_m2m1 = v_nsg 432 ! <phi^I_m1|Psi_nk> = proj%r (or proj%k) 433 ! |phi^J_m2> = wfcU 434 ! 435 wfcUaux(:,:) = (0.0_dp, 0.0_dp) 436 ! 437 off2 = offsetU(equiv_na2) 438 ! 439 DO m2 = 1, ldim_u(nt2) 440 ! 441 IF (m2.GT.2*Hubbard_l(nt2)+1) & 442 off2 = offsetU_back(equiv_na2) - 2*Hubbard_l(nt2) - 1 443 ! 444 IF (backall(nt2) .AND. & 445 m2.GT.(2*Hubbard_l(nt2)+1+2*Hubbard_l_back(nt2)+1)) & 446 off2 = offsetU_back1(equiv_na2) & 447 - 2*Hubbard_l(nt2) - 2 - 2*Hubbard_l_back(nt2) 448 ! 449 DO ig = 1, np 450 wfcUaux(ig,m2) = wfcU(ig,off2+m2) 451 ENDDO 452 ! 453 ENDDO 454 ! 455 off1 = offsetU(na1) 456 ! 457 IF (gamma_only) THEN 458 ! 459 projauxr(:,:) = 0.0_dp 460 ! 461 DO m1 = 1, ldim1 462 ! 463 IF (m1.GT.2*Hubbard_l(nt1)+1) & 464 off1 = offsetU_back(na1) - 2*Hubbard_l(nt1) - 1 465 ! 466 IF (backall(nt1) .AND. & 467 m1.GT.(2*Hubbard_l(nt1)+1+2*Hubbard_l_back(nt1)+1)) & 468 off1 = offsetU_back1(na1) & 469 - 2*Hubbard_l(nt1) - 2 - 2*Hubbard_l_back(nt1) 470 ! 471 projauxr(m1,:) = DBLE(proj%r(off1+m1,:)) 472 ! 473 ENDDO 474 ! 475 rvaux(:,:) = 0.0_dp 476 ! 477 DO m1 = 1, ldim1 478 ! 479 DO m2 = 1, ldim2 480 ! 481 rvaux(m2,m1) = DBLE(v_nsg(m2, m1, viz, na1, current_spin)) * 0.5d0 482 ! 483 ENDDO 484 ! 485 ENDDO 486 ! 487 rtemp(:,:) = 0.0_dp 488 ! 489 CALL DGEMM ('n','n', ldim2,mps,ldim1, 1.0_dp, & 490 rvaux,ldimx, projauxr,ldimx, 0.0_dp, rtemp, ldimx) 491 ! 492 CALL DGEMM ('n','n', 2*np, mps, ldim2, 1.0_dp, & 493 wfcUaux, 2*np, rtemp, ldimx, & 494 1.0_dp, hpsi, 2*ldap) 495 ! 496 ELSE 497 ! 498 projauxc(:,:) = (0.0_dp,0.0_dp) 499 ! 500 do m1 = 1,ldim1 501 ! 502 IF (m1.GT.2*Hubbard_l(nt1)+1) & 503 off1 = offsetU_back(na1) - 2*Hubbard_l(nt1) - 1 504 ! 505 IF (backall(nt1) .AND. & 506 m1.GT.(2*Hubbard_l(nt1)+1+2*Hubbard_l_back(nt1)+1)) & 507 off1 = offsetU_back1(na1) & 508 - 2*Hubbard_l(nt1) - 2 - 2*Hubbard_l_back(nt1) 509 ! 510 projauxc(m1,:) = proj%k(off1+m1,:) 511 ! 512 end do 513 ! 514 vaux(:,:) = (0.0_dp,0.0_dp) 515 ! 516 DO m1 = 1,ldim1 517 ! 518 DO m2 = 1,ldim2 519 ! 520 vaux(m2,m1) = v_nsg(m2, m1, viz, na1, current_spin) * 0.5d0 521 ! 522 ENDDO 523 ! 524 ENDDO 525 ! 526 ctemp(:,:) = (0.0_dp,0.0_dp) 527 ! 528 CALL ZGEMM ('n','n', ldim2,mps,ldim1, (1.0_dp,0.0_dp), & 529 vaux,ldimx, projauxc,ldimx, (0.0_dp,0.0_dp), ctemp, ldimx) 530 ! 531 CALL ZGEMM ('n','n', np, mps, ldim2, CONJG(phase), & 532 wfcUaux, np, ctemp, ldimx, (1.0_dp,0.0_dp), hpsi, ldap) 533 ! 534 ENDIF 535 ! 536 ENDIF 537 ! 538 ENDDO ! viz 539 ! 540 ENDIF !nt1 = ityp(na1) 541 ! 542 ENDDO !na1 543 ! 544 ENDIF 545 ! 546 ENDDO 547 ! 548 IF (gamma_only) THEN 549 DEALLOCATE (rtemp) 550 DEALLOCATE (projauxr) 551 DEALLOCATE (rvaux) 552 ELSE 553 DEALLOCATE (ctemp) 554 DEALLOCATE (projauxc) 555 DEALLOCATE (vaux) 556 ENDIF 557 ! 558 DEALLOCATE (wfcUaux) 559 ! 560 RETURN 561 ! 562END SUBROUTINE vhpsi_UV 563!------------------------------------------------------------------------- 564END SUBROUTINE vhpsi 565!------------------------------------------------------------------------- 566 567!------------------------------------------------------------------------- 568SUBROUTINE vhpsi_nc( ldap, np, mps, psip, hpsi ) 569 !----------------------------------------------------------------------- 570 !! Noncollinear version of \(\texttt{vhpsi} routine (A. Smogunov). 571 ! 572 USE kinds, ONLY: DP 573 USE ldaU, ONLY: Hubbard_lmax, Hubbard_l, is_Hubbard, nwfcU, & 574 wfcU, offsetU 575 USE scf, ONLY: v 576 USE ions_base, ONLY: nat, ntyp => nsp, ityp 577 USE noncollin_module, ONLY: npol 578 USE mp_bands, ONLY: intra_bgrp_comm 579 USE mp, ONLY: mp_sum 580 ! 581 IMPLICIT NONE 582 ! 583 INTEGER, INTENT(IN) :: ldap 584 !! leading dimension of arrays psip, hpsi 585 INTEGER, INTENT(IN) :: np 586 !! true dimension of psip, hpsi 587 INTEGER, INTENT(IN) :: mps 588 !! number of states psip 589 COMPLEX(DP), INTENT(IN) :: psip(ldap*npol,mps) 590 !! the wavefunction 591 COMPLEX(DP), INTENT(INOUT) :: hpsi(ldap*npol,mps) 592 !! Hamiltonian dot psi 593 ! 594 ! ... local variables 595 ! 596 INTEGER :: ibnd, na, nwfc, is1, is2, nt, m1, m2 597 COMPLEX(DP) :: temp 598 COMPLEX(DP), ALLOCATABLE :: proj(:,:) 599 ! 600 CALL start_clock('vhpsi') 601 ! 602 ALLOCATE( proj(nwfcU, mps) ) 603 ! 604!-- FIXME: to be replaced with ZGEMM 605! calculate <psi_at | phi_k> 606 DO ibnd = 1, mps 607 DO na = 1, nwfcU 608 proj(na, ibnd) = dot_product( wfcU(1:ldap*npol, na), psip(1:ldap*npol, ibnd)) 609 ENDDO 610 ENDDO 611#if defined(__MPI) 612 CALL mp_sum ( proj, intra_bgrp_comm ) 613#endif 614!-- 615 616 do ibnd = 1, mps 617 do na = 1, nat 618 nt = ityp (na) 619 if ( is_hubbard(nt) ) then 620 nwfc = 2 * Hubbard_l(nt) + 1 621 622 do is1 = 1, npol 623 do m1 = 1, nwfc 624 temp = 0.d0 625 do is2 = 1, npol 626 do m2 = 1, nwfc 627 temp = temp + v%ns_nc( m1, m2, npol*(is1-1)+is2, na) * & 628 proj(offsetU(na)+(is2-1)*nwfc+m2, ibnd) 629 enddo 630 enddo 631 call zaxpy (ldap*npol, temp, wfcU(1,offsetU(na)+(is1-1)*nwfc+m1),& 632 1, hpsi(1,ibnd),1) 633 enddo 634 enddo 635 636 endif 637 enddo 638 enddo 639 640 deallocate (proj) 641 CALL stop_clock('vhpsi') 642 643 return 644 ! 645END SUBROUTINE vhpsi_nc 646 647