1!
2! Copyright (C) 2015 GWW 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  SUBROUTINE adduspos_gamma_r(iw,jw,r_ij,ik,becp_iw,becp_jw)
10  !----------------------------------------------------------------------
11  !
12  !  This routine adds the US term < Psi_iw|r><r|Psi_jw>
13  !  to the array r_ij
14  !  this is a GAMMA only routine (i.e. r_ij is real)
15  !
16  USE kinds,                ONLY : DP
17  USE realus,               ONLY : tabp
18  USE ions_base,            ONLY : nat, ntyp => nsp, ityp
19  USE fft_base,             ONLY : dfftp
20  USE gvect,                ONLY : ngm, gg, g
21  USE lsda_mod,             ONLY : lsda, nspin, current_spin, isk
22  USE uspp,                 ONLY : okvan, becsum, nkb, ijtoh, indv_ijkb0
23  USE uspp_param,           ONLY : upf, lmaxq, nh
24  USE wvfct,                ONLY : wg
25  USE control_flags,        ONLY : gamma_only
26  USE wavefunctions, ONLY : psic
27  USE io_global,            ONLY : stdout
28  USE cell_base,            ONLY : omega
29  !
30  USE mp_bands,         ONLY : intra_bgrp_comm
31  USE mp,               ONLY : mp_sum
32  !
33  IMPLICIT NONE
34  !
35  !
36  INTEGER, INTENT(in) :: iw,jw!the states indices
37  REAL(kind=DP), INTENT(inout) :: r_ij(dfftp%nnr)!where to add the us term
38  INTEGER, INTENT(in) :: ik!spin index for spin polarized calculations NOT IMPLEMENTED YET
39  REAL(kind=DP), INTENT(in) ::  becp_iw( nkb)!overlap of wfcs with us  projectors
40  REAL(kind=DP), INTENT(in) ::  becp_jw( nkb)!overlap of wfcs with us  projectors
41
42  !     here the local variables
43  !
44
45  INTEGER :: na, nt, nhnt, ir, ih, jh, is , ia, mbia, irb, iqs, sizeqsave
46  INTEGER :: ikb, jkb, np
47  ! counters
48
49  ! work space for rho(G,nspin)
50  ! Fourier transform of q
51
52  IF (.not.okvan) RETURN
53  IF( .not.gamma_only) THEN
54     WRITE(stdout,*) ' adduspos_gamma_r is a gamma ONLY routine'
55     STOP
56  ENDIF
57
58  DO is=1,nspin
59     !
60     DO np = 1, ntyp
61        !
62        iqs = 0
63        !
64        IF ( upf(np)%tvanp ) THEN
65           !
66           DO ia = 1, nat
67              !
68              mbia = tabp(ia)%maxbox
69              nt = ityp(ia)
70              nhnt = nh(nt)
71              !
72              IF ( ityp(ia) /= np ) iqs=iqs+(nhnt+1)*nhnt*mbia/2
73              IF ( ityp(ia) /= np ) CYCLE
74              !
75              DO ih = 1, nhnt
76                 !
77                 ikb = indv_ijkb0(ia) + ih
78                 !
79                 DO jh = ih, nhnt
80                    !
81                    jkb = indv_ijkb0(ia) + jh
82                    !
83                    DO ir = 1, mbia
84                       !
85                       irb = tabp(ia)%box(ir)
86                       iqs = iqs + 1
87                       !
88                       r_ij(irb) = r_ij(irb) + tabp(ia)%qr(ir,ijtoh(ih,jh,np))&
89                                              *becp_iw(ikb)*becp_jw(jkb)*omega
90                       !
91                       IF ( ih /= jh ) THEN
92                          r_ij(irb) = r_ij(irb) + tabp(ia)%qr(ir,ijtoh(ih,jh,np))&
93                                                 *becp_iw(jkb)*becp_jw(ikb)*omega
94                       ENDIF
95                    ENDDO
96                 ENDDO
97              ENDDO
98              !
99           ENDDO
100           !
101        ENDIF
102     ENDDO
103  ENDDO
104  !
105  RETURN
106  !
107  END SUBROUTINE adduspos_gamma_r
108  !
109  SUBROUTINE adduspos_r(r_ij,becp_iw,becp_jw)
110  !----------------------------------------------------------------------
111  !
112  !  This routine adds the US term < Psi_iw|r><r|Psi_jw>
113  !  to the array r_ij
114  USE kinds,                ONLY : DP
115  USE realus,               ONLY : tabp
116  USE ions_base,            ONLY : nat, ntyp => nsp, ityp
117  USE fft_base,             ONLY : dfftp
118  USE gvect,                ONLY : ngm, gg, g
119  USE lsda_mod,             ONLY : lsda, nspin, current_spin, isk
120  USE uspp,                 ONLY : okvan, becsum, nkb, ijtoh, indv_ijkb0
121  USE uspp_param,           ONLY : upf, lmaxq, nh
122  USE wvfct,                ONLY : wg
123  USE control_flags,        ONLY : gamma_only
124  USE wavefunctions, ONLY : psic
125  USE cell_base,            ONLY : omega
126  !
127  IMPLICIT NONE
128  !
129  COMPLEX(kind=DP), INTENT(inout) :: r_ij(dfftp%nnr)!where to add the us term
130  COMPLEX(kind=DP), INTENT(in) ::  becp_iw( nkb)!overlap of wfcs with us  projectors
131  COMPLEX(kind=DP), INTENT(in) ::  becp_jw( nkb)!overlap of wfcs with us  projectors
132  !     here the local variables
133  !
134  INTEGER :: na, ia, nt, nhnt, ir, ih, jh, is, mbia, irb, iqs
135  INTEGER :: ikb, jkb, np
136  ! counters
137
138  ! work space for rho(G,nspin)
139  ! Fourier transform of q
140
141  IF (.not.okvan) RETURN
142
143  DO is=1,nspin
144     !
145     DO np = 1, ntyp
146        !
147        iqs = 0
148        !
149        IF ( upf(np)%tvanp ) THEN
150           !
151           DO ia = 1, nat
152              !
153              mbia = tabp(ia)%maxbox
154              nt = ityp(ia)
155              nhnt = nh(nt)
156              !
157              IF ( ityp(ia) /= np ) iqs=iqs+(nhnt+1)*nhnt*mbia/2
158              IF ( ityp(ia) /= np ) CYCLE
159              !
160              DO ih = 1, nhnt
161                 !
162                 ikb = indv_ijkb0(ia) + ih
163                 DO jh = ih, nhnt
164                    !
165                    jkb = indv_ijkb0(ia) + jh
166                    !
167                    DO ir = 1, mbia
168                       !
169                       irb = tabp(ia)%box(ir)
170                       iqs = iqs + 1
171                       !
172                       r_ij(irb) = r_ij(irb) + tabp(ia)%qr(ir,ijtoh(ih,jh,np))&
173                                              *conjg(becp_iw(ikb))*becp_jw(jkb)*omega
174                       !
175                       IF ( ih /= jh ) THEN
176                          r_ij(irb) = r_ij(irb) + tabp(ia)%qr(ir,ijtoh(ih,jh,np))&
177                                                 *conjg(becp_iw(jkb))*becp_jw(ikb)*omega
178                       ENDIF
179                    ENDDO
180                 ENDDO
181              ENDDO
182              !
183           ENDDO
184           !
185        ENDIF
186     ENDDO
187  ENDDO
188  !
189  RETURN
190  END SUBROUTINE adduspos_r
191  !
192  SUBROUTINE adduspos_real(sca,qq_op,becp_iw,becp_jw)
193  !----------------------------------------------------------------------
194  !
195  !  This routine adds the US term < Psi_iw|r><r|Psi_jw>
196  !  to the array r_ij
197  USE kinds,                ONLY : DP
198  USE realus,               ONLY : tabp
199  USE ions_base,            ONLY : nat, ntyp => nsp, ityp
200  USE gvect,                ONLY : ngm, gg, g
201  USE lsda_mod,             ONLY : lsda, nspin, current_spin, isk
202  USE uspp,                 ONLY : okvan, becsum, nkb, indv_ijkb0
203  USE uspp_param,           ONLY : upf, lmaxq, nh, nhm
204  USE wvfct,                ONLY : wg
205  USE control_flags,        ONLY : gamma_only
206  USE wavefunctions, ONLY : psic
207  USE cell_base,            ONLY : omega
208  !
209  USE mp_bands,             ONLY : intra_bgrp_comm
210  USE mp,                   ONLY : mp_sum
211  !
212  IMPLICIT NONE
213  !
214  REAL(kind=DP), INTENT(inout) :: sca!where to add the us term
215  REAL(kind=DP), INTENT(in) ::  becp_iw( nkb)!overlap of wfcs with us  projectors
216  REAL(kind=DP), INTENT(in) ::  becp_jw( nkb)!overlap of wfcs with us  projectors
217  REAL(kind=DP), INTENT(in) ::  qq_op(nhm, nhm,nat)!US augmentation charges
218
219  !     here the local variables
220  !
221
222  INTEGER :: na, ia, nhnt, nt, ih, jh, is, mbia
223  INTEGER :: ikb, jkb, np
224  ! counters
225
226  ! work space for rho(G,nspin)
227  ! Fourier transform of q
228
229  IF (.not.okvan) RETURN
230
231  DO is=1,nspin
232     !
233     DO np = 1, ntyp
234        !
235        IF ( upf(np)%tvanp ) THEN
236           !
237           DO ia = 1, nat
238              !
239              IF ( ityp(ia) /= np ) CYCLE
240              !
241              mbia = tabp(ia)%maxbox
242              nt = ityp(ia)
243              nhnt = nh(nt)
244              !
245              DO ih = 1, nhnt
246                 !
247                 ikb = indv_ijkb0(ia) + ih
248                 DO jh = ih, nhnt
249                    !
250                    jkb = indv_ijkb0(ia) + jh
251                    !
252                    sca = sca + qq_op(ih,jh,ia) * becp_iw(ikb)*becp_jw(jkb)
253                    !
254                    IF ( ih /= jh ) THEN
255                       sca = sca + qq_op(jh,ih,ia) * becp_iw(ikb)*becp_jw(jkb)
256                    ENDIF
257                    !
258                 ENDDO
259              ENDDO
260              !
261           ENDDO
262           !
263        ENDIF
264     ENDDO
265  ENDDO
266  !
267  RETURN
268  !
269  END SUBROUTINE adduspos_real
270  !
271