1! 2! Copyright (C) 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!----------------------------------------------------------------------- 10SUBROUTINE orthogonalize_omega(dvpsi, evq, ikk, ikq, dpsi, npwq, w) 11!------------------------------------------------------------------------ 12 ! 13 ! This routine ortogonalizes dvpsi to the valence states: ps = <evq|dvpsi> 14 ! It should be quite general. It works for metals and insulators, with 15 ! NC as well as with US PP, both SR or FR. 16 ! Note that on output it changes sign. So it applies -P^+_c. 17 ! 18 ! NB: IN/OUT is dvpsi ; dpsi is used as work_space 19 ! 20USE kinds, ONLY : DP 21USE klist, ONLY : lgauss, degauss, ngauss 22USE noncollin_module, ONLY : noncolin, npol 23USE wvfct, ONLY : npwx, nbnd, et 24USE ener, ONLY : ef 25USE control_lr, ONLY : alpha_pv, nbnd_occ 26USE becmod, ONLY : bec_type, becp, calbec 27USE uspp, ONLY : vkb, okvan 28USE mp_bands, ONLY : intra_bgrp_comm 29USE mp, ONLY : mp_sum 30USE control_flags, ONLY : gamma_only 31USE gvect, ONLY : gstart 32! 33IMPLICIT NONE 34INTEGER, INTENT(IN) :: ikk, ikq ! the index of the k and k+q points 35INTEGER, INTENT(IN) :: npwq ! the number of plane waves for q 36COMPLEX(DP), INTENT(IN) :: evq(npwx*npol,nbnd) 37COMPLEX(DP), INTENT(INOUT) :: dvpsi(npwx*npol,nbnd) 38COMPLEX(DP), INTENT(INOUT) :: dpsi(npwx*npol,nbnd) ! work space allocated by 39 ! the calling routine 40COMPLEX(DP), INTENT(IN) :: w 41 42COMPLEX(DP), ALLOCATABLE :: ps(:,:) 43REAL(DP), ALLOCATABLE :: ps_r(:,:) 44 45INTEGER :: ibnd, jbnd, nbnd_eff 46REAL(DP) :: wg1, w0g, wgp, deltae, theta 47REAL(DP), EXTERNAL :: w0gauss, wgauss 48COMPLEX(DP) :: wwg 49! functions computing the delta and theta function 50 51CALL start_clock ('ortho') 52IF (gamma_only) THEN 53 ALLOCATE(ps_r(nbnd,nbnd)) 54 ps_r = 0.0_DP 55ENDIF 56 57ALLOCATE(ps(nbnd,nbnd)) 58ps = (0.0_DP, 0.0_DP) 59 60! 61if (lgauss) then 62 ! 63 IF (gamma_only) CALL errore ('orthogonalize', "degauss with gamma & 64 & point algorithms",1) 65 ! 66 ! metallic case 67 ! 68 IF (noncolin) THEN 69 CALL zgemm( 'C', 'N', nbnd, nbnd_occ (ikk), npwx*npol, (1.d0,0.d0), & 70 evq, npwx*npol, dvpsi, npwx*npol, (0.d0,0.d0), ps, nbnd ) 71 ELSE 72 CALL zgemm( 'C', 'N', nbnd, nbnd_occ (ikk), npwq, (1.d0,0.d0), & 73 evq, npwx, dvpsi, npwx, (0.d0,0.d0), ps, nbnd ) 74 END IF 75 ! 76 77 DO ibnd = 1, nbnd_occ (ikk) 78 wg1 = wgauss ((ef-et(ibnd,ikk)) / degauss, ngauss) 79 w0g = w0gauss((ef-et(ibnd,ikk)) / degauss, ngauss) / degauss 80 DO jbnd = 1, nbnd 81 wgp = wgauss ( (ef - et (jbnd, ikq) ) / degauss, ngauss) 82 deltae = et (jbnd, ikq) - et (ibnd, ikk) 83 theta = wgauss (deltae / degauss, 0) 84 wwg = CMPLX(wg1 * (1.d0 - theta) + wgp * theta,0.0_DP) 85 IF (jbnd <= nbnd_occ (ikq) ) THEN 86 wwg = wwg + alpha_pv * theta * (wgp - wg1) / (deltae - w) 87 ENDIF 88 ! 89 ps(jbnd,ibnd) = wwg * ps(jbnd,ibnd) 90 ! 91 ENDDO 92 IF (noncolin) THEN 93 CALL dscal (2*npwx*npol, wg1, dvpsi(1,ibnd), 1) 94 ELSE 95 call dscal (2*npwq, wg1, dvpsi(1,ibnd), 1) 96 END IF 97 END DO 98 nbnd_eff=nbnd 99ELSE 100 ! 101 ! insulators 102 ! 103 IF (noncolin) THEN 104 CALL zgemm( 'C', 'N',nbnd_occ(ikq), nbnd_occ(ikk), npwx*npol, & 105 (1.d0,0.d0), evq, npwx*npol, dvpsi, npwx*npol, & 106 (0.d0,0.d0), ps, nbnd ) 107 ELSEIF (gamma_only) THEN 108 CALL dgemm( 'C', 'N', nbnd_occ(ikq), nbnd_occ (ikk), 2*npwq, & 109 2.0_DP, evq, 2*npwx, dvpsi, 2*npwx, & 110 0.0_DP, ps_r, nbnd ) 111 IF (gstart == 2 ) THEN 112 CALL DGER( nbnd_occ(ikq), nbnd_occ (ikk), -1.0_DP, evq, & 113 & 2*npwq, dvpsi, 2*npwx, ps_r, nbnd ) 114 ENDIF 115 ELSE 116 CALL zgemm( 'C', 'N', nbnd_occ(ikq), nbnd_occ (ikk), npwq, & 117 (1.d0,0.d0), evq, npwx, dvpsi, npwx, & 118 (0.d0,0.d0), ps, nbnd ) 119 END IF 120 nbnd_eff=nbnd_occ(ikk) 121END IF 122IF (gamma_only) THEN 123 call mp_sum(ps_r(:,:),intra_bgrp_comm) 124ELSE 125 call mp_sum(ps(:,1:nbnd_eff),intra_bgrp_comm) 126ENDIF 127! 128! dpsi is used as work space to store S|evc> 129! 130IF (okvan) CALL calbec ( npwq, vkb, evq, becp, nbnd_eff) 131CALL s_psi (npwx, npwq, nbnd_eff, evq, dpsi) 132! 133! |dvspi> = -(|dvpsi> - S|evq><evq|dvpsi>) 134! 135if (lgauss) then 136 ! 137 ! metallic case 138 ! 139 IF (noncolin) THEN 140 CALL zgemm( 'N', 'N', npwx*npol, nbnd_occ(ikk), nbnd, & 141 (1.d0,0.d0), dpsi, npwx*npol, ps, nbnd, (-1.0d0,0.d0), & 142 dvpsi, npwx*npol ) 143 ELSE 144 CALL zgemm( 'N', 'N', npwq, nbnd_occ(ikk), nbnd, & 145 (1.d0,0.d0), dpsi, npwx, ps, nbnd, (-1.0d0,0.d0), & 146 dvpsi, npwx ) 147 END IF 148ELSE 149 ! 150 ! Insulators: note that nbnd_occ(ikk)=nbnd_occ(ikq) in an insulator 151 ! 152 IF (noncolin) THEN 153 CALL zgemm( 'N', 'N', npwx*npol, nbnd_occ(ikk), nbnd_occ(ikk), & 154 (1.d0,0.d0),dpsi,npwx*npol,ps,nbnd,(-1.0d0,0.d0), & 155 dvpsi, npwx*npol ) 156 ELSEIF (gamma_only) THEN 157 ps = CMPLX (ps_r,0.0_DP, KIND=DP) 158 CALL ZGEMM( 'N', 'N', npwq, nbnd_occ(ikk), nbnd_occ(ikk), & 159 (1.d0,0.d0), dpsi, npwx, ps, nbnd, (-1.0d0,0.d0), & 160 dvpsi, npwx ) 161 ELSE 162 CALL zgemm( 'N', 'N', npwq, nbnd_occ(ikk), nbnd_occ(ikk), & 163 (1.d0,0.d0), dpsi, npwx, ps, nbnd, (-1.0d0,0.d0), & 164 dvpsi, npwx ) 165 END IF 166ENDIF 167 168DEALLOCATE(ps) 169CALL stop_clock ('ortho') 170RETURN 171END SUBROUTINE orthogonalize_omega 172