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