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