1! 2! Copyright (C) 2001-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 solve_e2 11 !----------------------------------------------------------------------- 12 ! 13 ! Self consistent cycle to compute the second order derivatives 14 ! of the wavefunctions with respect to electric fields 15 ! 16 USe kinds, ONLY : DP 17 USE io_global, ONLY : stdout 18 USE cell_base, ONLY : tpiba2 19 USE klist, ONLY : ltetra, lgauss, wk, xk, ngk, igk_k 20 USE lsda_mod, ONLY : lsda, nspin 21 USE gvect, ONLY : g 22 USE gvecs, ONLY : doublegrid 23 USE fft_base, ONLY : dfftp, dffts 24 USE fft_interfaces, ONLY : fft_interpolate 25 USE wvfct, ONLY : npwx, nbnd, et 26 USE buffers, ONLY: get_buffer 27 USE ions_base, ONLY: nat 28 USE uspp, ONLY: okvan, nkb, vkb 29 USE uspp_param,ONLY : nhm 30 USE wavefunctions, ONLY: evc 31 USE control_ph, ONLY : convt, nmix_ph, alpha_mix, tr2_ph, & 32 niter_ph, rec_code, flmixdpot, rec_code_read 33 USE units_lr, ONLY : lrwfc, iuwfc 34 USE ramanm, ONLY : lrba2, iuba2, lrd2w, iud2w 35 USE recover_mod, ONLY : read_rec, write_rec 36 37 USE check_stop, ONLY: check_stop_now 38 USE mp_pools, ONLY : inter_pool_comm 39 USE mp_bands, ONLY : intra_bgrp_comm 40 USE mp, ONLY : mp_sum 41 42 USE eqv, ONLY : dpsi, dvpsi 43 USE qpoint, ONLY : nksq, ikks, ikqs 44 USE control_lr, ONLY : nbnd_occ, lgamma 45 USE dv_of_drho_lr 46 47 implicit none 48 49 real(DP) :: thresh, weight, avg_iter, dr2 50 ! convergence threshold for the solution of the 51 ! linear system 52 ! used for summation over k points 53 ! average number of iterations 54 ! convergence limit 55 56 complex(DP) , pointer :: dvscfin (:,:,:), dvscfins (:,:,:) 57 ! change of the scf potential (input) 58 ! change of the scf potential (smooth) 59 60 complex(DP) , allocatable :: dvscfout (:,:,:), dbecsum (:,:), & 61 aux1 (:) 62 ! change of the scf potential (output) 63 ! auxiliary space 64 ! auxiliary space 65 66 logical :: exst 67 ! used to open the recover file 68 69 integer :: npw, npwq, ikk, ikq 70 integer :: kter, iter0, ipol, ibnd, iter, ik, is, ig, iig, irr, ir, nrec, ios 71 ! counter on iterations 72 ! counter on perturbations 73 ! counter on bands 74 ! counter on iterations 75 ! counter on k points 76 ! counter on G vectors 77 ! counter on g vectors 78 ! counter on mesh points 79 ! the record number 80 ! integer variable for I/O control 81 82 if (lsda) call errore ('solve_e2', ' LSDA not implemented', 1) 83 if (okvan) call errore ('solve_e2', ' Ultrasoft PP not implemented', 1) 84 85 call start_clock('solve_e2') 86 allocate (dvscfin( dfftp%nnr, nspin, 6)) 87 if (doublegrid) then 88 allocate (dvscfins(dffts%nnr, nspin, 6)) 89 else 90 dvscfins => dvscfin 91 endif 92 allocate (dvscfout( dfftp%nnr , nspin, 6)) 93 allocate (dbecsum( nhm*(nhm+1)/2, nat)) 94 allocate (aux1(dffts%nnr)) 95 convt=.FALSE. 96 if (rec_code_read == -10) then 97 ! restarting in Raman 98 CALL read_rec(dr2, iter0, 6, dvscfin, dvscfins) 99 else 100 iter0 = 0 101 end if 102 if (convt) goto 155 103 ! 104 if ( (lgauss .or. ltetra) .or..not.lgamma) & 105 call errore ('solve_e2', 'called in the wrong case', 1) 106 ! 107 ! The outside loop is over the iterations 108 ! 109 110 do kter = 1, niter_ph 111 112 iter = kter + iter0 113 avg_iter = 0.d0 114 115 dvscfout (:,:,:) = (0.d0, 0.d0) 116 dbecsum (:,:) = (0.d0, 0.d0) 117 118 do ik = 1, nksq 119 ! in this routine, ikk=ikq=ik always (q=0) 120 ikk = ikks(ik) 121 ikq = ikqs(ik) 122 ! 123 ! reads unperturbed wavefunctions psi_k in G_space, for all bands 124 ! 125 if (nksq.gt.1) call get_buffer(evc, lrwfc, iuwfc, ikk) 126 npw = ngk(ikk) 127 npwq= ngk(ikq) 128 ! 129 ! compute beta functions and kinetic energy for k-point ikk 130 ! needed by h_psi, called by pcgreen 131 ! 132 CALL init_us_2 (npw, igk_k(1,ikk), xk (1, ikk), vkb) 133 CALL g2_kin(ikk) 134 ! 135 ! The counter on the polarizations runs only on the 6 inequivalent 136 ! indexes --see the comment on raman.F-- 137 ! 138 do ipol = 1, 6 139 nrec = (ipol - 1) * nksq + ik 140 141 if (kter.eq.1) then 142 dpsi (:,:) = (0.d0, 0.d0) 143 else 144 call davcio (dpsi, lrd2w, iud2w, nrec, -1) 145 endif 146 147 if (iter.eq.1) then 148 dvscfin (:,:,:) = (0.d0, 0.d0) 149 call davcio (dvpsi, lrba2, iuba2, nrec, -1) 150 thresh = 1.0d-2 151 else 152 call davcio (dvpsi, lrba2, iuba2, nrec, -1) 153 do ibnd = 1, nbnd_occ (ik) 154 call cft_wave (ik, evc (1, ibnd), aux1, +1) 155 do ir = 1, dffts%nnr 156 aux1 (ir) = aux1 (ir) * dvscfins (ir, 1, ipol) 157 enddo 158 call cft_wave (ik, dvpsi (1, ibnd), aux1, -1) 159 enddo 160 thresh = min (0.1d0 * sqrt(dr2), 1.0d-2) 161 endif 162 163 call pcgreen (avg_iter, thresh, ik, et (1, ik) ) 164 call davcio ( dpsi, lrd2w, iud2w, nrec, +1) 165 ! 166 ! calculates dvscf, sum over k => dvscf_q_ipert 167 ! 168 weight = wk (ik) 169 call incdrhoscf (dvscfout (1,1,ipol), weight, ik, & 170 dbecsum (1, 1), dpsi) 171 enddo ! on perturbations 172 enddo ! on k points 173 call mp_sum ( dbecsum, intra_bgrp_comm ) 174 if (doublegrid) then 175 do is = 1, nspin 176 do ipol = 1, 6 177 call fft_interpolate (dffts, dvscfout (:, is, ipol), dfftp, dvscfout (:, is, ipol)) 178 enddo 179 enddo 180 endif 181 ! 182 ! After the loop over the perturbations we have the change of the pote 183 ! for all the modes, and we symmetrize this potential 184 ! 185 call mp_sum ( dvscfout, inter_pool_comm ) 186 do ipol = 1, 6 187 call dv_of_drho (dvscfout (1, 1, ipol), .false.) 188 enddo 189 190 call psyme2(dvscfout) 191 ! 192 ! Mixing with the old potential 193 ! 194 call mix_potential (2 * 6 * dfftp%nnr* nspin, dvscfout, dvscfin, & 195 alpha_mix (kter), dr2, 6 * tr2_ph, iter, & 196 nmix_ph, flmixdpot, convt) 197 198 if (doublegrid) then 199 do is = 1, nspin 200 do ipol = 1, 6 201 call fft_interpolate (dfftp,dvscfin (:, is, ipol), dffts, dvscfins (:, is, ipol)) 202 enddo 203 enddo 204 end if 205 206 write (6, "(//,5x,' iter # ',i3, & 207 & ' av.it.: ',f5.1)") iter, avg_iter / (6.d0 * nksq) 208 dr2 = dr2 / 6 209 write (6, "(5x,' thresh=',es10.3, ' alpha_mix = ',f6.3, & 210 & ' |ddv_scf|^2 = ',es10.3 )") thresh, alpha_mix (kter), dr2 211 ! 212 FLUSH( stdout ) 213 ! 214 ! rec_code: state of the calculation 215 ! rec_code=-10 to -19 Raman 216 rec_code=-10 217 CALL write_rec('solve_e2..', irr, dr2, iter, convt, 6, dvscfin) 218 219 if ( check_stop_now() ) call stop_smoothly_ph (.false.) 220 if ( convt ) goto 155 221 222 enddo 223155 continue 224 deallocate (dvscfin ) 225 if (doublegrid) deallocate (dvscfins ) 226 deallocate (dvscfout ) 227 deallocate (dbecsum ) 228 deallocate (aux1 ) 229 230 call stop_clock('solve_e2') 231 232 return 233end subroutine solve_e2 234