1!
2! Copyright (C) 2001-2016 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!----------------------------------------------------------------------
9subroutine dvpsi_e2
10  !-----------------------------------------------------------------------
11  !
12  ! This routine shold be called before the self-consistent cycle used to
13  ! compute the second derivative of the wavefunctions with respect to
14  ! electric-fields. It computes that part of the potential that remains
15  ! constant during the cycle.
16  !
17  USE kinds,           ONLY : DP
18  USE cell_base,       ONLY : omega
19  USE klist,           ONLY : wk, ngk
20  USE gvecs,           ONLY : doublegrid
21  USE wvfct,           ONLY : npwx, nbnd
22  USE wavefunctions, ONLY: evc
23  USE buffers,         ONLY : get_buffer
24  USE fft_base,        ONLY : dfftp, dffts
25  USE fft_interfaces,  ONLY : fft_interpolate
26  USE scf,             ONLY : rho
27  USE qpoint,          ONLY : nksq
28  USE units_ph,        ONLY : lrdrho, iudrho, lrdwf, iudwf
29  USE units_lr,        ONLY : iuwfc, lrwfc
30  USE control_lr,      ONLY : nbnd_occ
31  USE ramanm,          ONLY : lrba2, iuba2, lrchf, iuchf, a1j, a2j
32  USE mp_pools,        ONLY : my_pool_id, inter_pool_comm
33  USE mp_bands,        ONLY : intra_bgrp_comm
34  USE mp,        ONLY: mp_sum
35  USE dv_of_drho_lr
36
37  implicit none
38
39  INTEGER :: npw, npwq
40  integer :: ik, ipa, ipb, ir, ibnd, jbnd, nrec
41  ! counter on k-points
42  ! counter on polarizations
43  ! counter on points of the real-space mesh
44  ! counter on bands
45  ! the record number
46  real(DP), allocatable  :: raux6 (:,:), d2muxc (:)
47  ! function on the real space smooth-mesh
48  ! second derivative of the XC-potential
49  real(DP) ::  d2mxc, rhotot
50  ! external function
51  ! total charge on a point
52  complex(DP), allocatable :: depsi (:,:,:), auxg (:,:), auxs1 (:), &
53               auxs2 (:), aux3s (:,:), aux3 (:,:), ps (:,:,:,:)
54  ! d |psi> / dE  (E=electric field)
55  ! chi-wavefunction
56  ! function on the real space smooth-mesh
57  ! function on the real space smooth-mesh
58  ! function on the real space smooth-mesh
59  ! function on the real space thick-mesh
60  complex(DP), pointer :: aux6s (:,:), aux6 (:,:)
61  ! function on the real space smooth-mesh
62  ! function on the real space thick-mesh
63  complex(DP) :: tmp, weight
64  ! working space
65  ! weight in k-point summation
66  !
67  call start_clock('dvpsi_e2')
68  !
69  ! First, calculates the second derivative of the charge-density
70  ! -only the part that does not depend on the self-consistent cycle-
71  !
72  allocate (raux6  (dffts%nnr,6))
73  allocate (depsi  (npwx,nbnd,3))
74  allocate (aux3s  (dffts%nnr,3))
75  allocate (ps     (nbnd,nbnd,3,3))
76
77  raux6 (:,:) = 0.d0
78  do ik = 1, nksq
79     npw = ngk(ik)
80     npwq = npw
81     if (nksq.gt.1) call get_buffer (evc, lrwfc, iuwfc, ik)
82     weight = 2.d0 * wk(ik) / omega
83
84     do ipa = 1, 3
85        nrec = (ipa - 1) * nksq + ik
86        call get_buffer (depsi (1, 1, ipa), lrdwf, iudwf, nrec)
87     enddo
88
89     do ibnd = 1, nbnd_occ (ik)
90        do ipa = 1, 3
91           call cft_wave (ik, depsi (1, ibnd, ipa), aux3s (1, ipa), +1)
92        enddo
93        do ipa = 1, 6
94           do ir = 1, dffts%nnr
95              tmp = CONJG(aux3s (ir, a1j (ipa))) * aux3s (ir, a2j (ipa))
96              raux6 (ir, ipa) = raux6 (ir, ipa) + weight *  DBLE (tmp)
97           enddo
98        enddo
99     enddo
100
101     do ipa = 1, 3
102        do ipb = 1, 3
103           CALL zgemm( 'C', 'N', nbnd_occ (ik), nbnd_occ (ik), npwq, &
104                (1.d0,0.d0), depsi(1,1, ipa), npwx, depsi(1,1,ipb), npwx, &
105                (0.d0,0.d0), ps(1,1,ipa,ipb), nbnd )
106        enddo
107     enddo
108     call mp_sum ( ps, intra_bgrp_comm )
109
110     do ibnd = 1, nbnd_occ (ik)
111        call cft_wave (ik, evc (1, ibnd), aux3s (1,1), +1)
112        do jbnd = 1, nbnd_occ (ik)
113           call cft_wave (ik, evc (1, jbnd), aux3s (1,2), +1)
114           do ipa = 1, 6
115              do ir = 1, dffts%nnr
116                 tmp =  aux3s (ir,1) *                           &
117                        ps(ibnd, jbnd, a1j (ipa), a2j (ipa)) *   &
118                        CONJG(aux3s (ir,2))
119                 raux6 (ir, ipa) = raux6 (ir, ipa) - weight *  DBLE (tmp)
120              enddo
121           enddo
122        enddo
123     enddo
124
125  enddo
126
127  deallocate (depsi)
128  deallocate (aux3s)
129  deallocate (ps)
130
131  !
132  ! Multiplies the charge with the potential
133  !
134  if (doublegrid) then
135     allocate (auxs1  (dffts%nnr))
136     allocate (aux6   (dfftp%nnr,6))
137  else
138     allocate (aux6s  (dffts%nnr,6))
139     aux6 => aux6s
140  endif
141
142  do ipa = 1, 6
143     if (doublegrid) then
144        do ir = 1, dffts%nnr
145           auxs1 (ir) = CMPLX(raux6 (ir, ipa), 0.d0,kind=DP)
146        enddo
147        call fft_interpolate (dffts, auxs1, dfftp, aux6 (:, ipa))
148     else
149        do ir = 1, dffts%nnr
150           aux6 (ir, ipa) = CMPLX(raux6 (ir, ipa), 0.d0,kind=DP)
151        enddo
152     endif
153     call dv_of_drho (aux6(:, ipa), .false.)
154  enddo
155
156  if (doublegrid) deallocate (auxs1)
157  deallocate (raux6)
158
159  !
160  ! Calculates the term depending on the third derivative of the
161  !                     Exchange-correlation energy
162  !
163  allocate (d2muxc (dfftp%nnr))
164  allocate (aux3   (dfftp%nnr,3))
165  do ipa = 1, 3
166     call davcio_drho (aux3 (1, ipa), lrdrho, iudrho, ipa, -1)
167  enddo
168
169#if defined(__MPI)
170  if (my_pool_id .ne. 0) goto 100
171#endif
172  d2muxc (:) = 0.d0
173  do ir = 1, dfftp%nnr
174!     rhotot = rho%of_r(ir,1) + rho_core(ir)
175     rhotot = rho%of_r(ir,1)
176     if ( rhotot.gt. 1.d-30 ) d2muxc(ir)= d2mxc( rhotot)
177     if ( rhotot.lt.-1.d-30 ) d2muxc(ir)=-d2mxc(-rhotot)
178  enddo
179
180  do ipa = 1, 6
181     do ir = 1, dfftp%nnr
182        aux6 (ir, ipa) = aux6 (ir, ipa) + d2muxc (ir) * &
183                   aux3 (ir, a1j (ipa)) * aux3 (ir, a2j (ipa))
184     enddo
185  enddo
186
187 100  continue
188  call mp_sum ( aux6, inter_pool_comm )
189  call psyme2 (aux6)
190
191  deallocate (d2muxc)
192  deallocate (aux3)
193
194
195  if (doublegrid) then
196     allocate (aux6s  (dffts%nnr,6))
197     do ipa = 1, 6
198        call fft_interpolate (dfftp, aux6 (:, ipa), dffts, aux6s (:, ipa))
199     enddo
200     deallocate (aux6)
201  endif
202  !
203  ! Multiplies the obtained potential with the wavefunctions and
204  ! writes the results on iuba2; a faster way of proceeding would
205  ! be that of keeping the potential in memory and use it directly in
206  ! solve_e2
207  !
208  allocate (auxg   (npwx,nbnd))
209  allocate (auxs1  (dffts%nnr))
210  allocate (auxs2  (dffts%nnr))
211
212  do ik = 1, nksq
213     npw = ngk(ik)
214     npwq = npw
215     if (nksq.gt.1) call get_buffer (evc, lrwfc, iuwfc, ik)
216     do ipa = 1, 6
217        nrec = (ipa - 1) * nksq + ik
218        call davcio (auxg, lrchf, iuchf, nrec, -1)
219        do ibnd = 1, nbnd_occ (ik)
220           call cft_wave (ik, evc (1, ibnd), auxs1, +1)
221           do ir = 1, dffts%nnr
222              auxs2 (ir) = auxs1 (ir) * aux6s (ir, ipa)
223           enddo
224           call cft_wave (ik, auxg (1, ibnd), auxs2, -1)
225        enddo
226        nrec = (ipa - 1) * nksq + ik
227        call davcio (auxg, lrba2, iuba2, nrec, +1)
228     enddo
229  enddo
230
231  deallocate (auxg)
232  deallocate (auxs1)
233  deallocate (auxs2)
234
235  deallocate (aux6s)
236  call stop_clock('dvpsi_e2')
237
238  return
239end subroutine dvpsi_e2
240