1! 2! Copyright (C) 2001-2007 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_e_fpol ( iw ) 11 !----------------------------------------------------------------------- 12 ! 13 ! This routine is a driver for the solution of the linear system which 14 ! defines the change of the wavefunction due to an electric field. 15 ! It performs the following tasks: 16 ! a) computes the bare potential term x | psi > 17 ! b) adds to it the screening term Delta V_{SCF} | psi > 18 ! c) applies P_c^+ (orthogonalization to valence states) 19 ! d) calls cgsolve_all to solve the linear system 20 ! e) computes Delta rho, Delta V_{SCF} and symmetrizes them 21 ! 22 USE kinds, ONLY : DP 23 USE ions_base, ONLY : nat 24 USE io_global, ONLY : stdout, ionode 25 USE io_files, ONLY : prefix, diropn 26 USE buffers, ONLY : get_buffer, save_buffer 27 USE check_stop, ONLY : check_stop_now 28 USE wavefunctions, ONLY : evc 29 USE cell_base, ONLY : tpiba2 30 USE klist, ONLY : ltetra, lgauss, nkstot, wk, xk, ngk, igk_k 31 USE lsda_mod, ONLY : lsda, nspin, current_spin, isk 32 USE fft_base, ONLY : dffts, dfftp 33 USE fft_interfaces, ONLY : fwfft, invfft, fft_interpolate 34 USE gvect, ONLY : g 35 USE gvecs, ONLY : doublegrid 36 USE becmod, ONLY : becp, calbec 37 USE wvfct, ONLY : npwx, nbnd, g2kin, et 38 USE uspp, ONLY : okvan, vkb 39 USE uspp_param, ONLY : nhm 40 USE control_ph, ONLY : nmix_ph, tr2_ph, alpha_mix, convt, & 41 niter_ph, & 42 rec_code, flmixdpot 43 USE output, ONLY : fildrho 44 USE qpoint, ONLY : nksq 45 USE units_ph, ONLY : lrdwf, iudwf, iudrho, lrdrho 46 USE units_lr, ONLY : iuwfc, lrwfc 47 USE mp_pools, ONLY : inter_pool_comm 48 USE mp_bands, ONLY : intra_bgrp_comm 49 USE mp, ONLY : mp_sum 50 51 USE eqv, ONLY : dpsi, dvpsi 52 USE control_lr, ONLY : nbnd_occ, lgamma 53 USE dv_of_drho_lr 54 55 implicit none 56 57 real(DP) :: thresh, anorm, averlt, dr2 58 ! thresh: convergence threshold 59 ! anorm : the norm of the error 60 ! averlt: average number of iterations 61 ! dr2 : self-consistency error 62 63 complex(kind=DP), allocatable :: etc(:,:), h_diag(:,:) 64 ! the eigenvalues plus imaginary frequency 65 ! the diagonal part of the Hamiltonian which becomes complex now 66 67 68 complex(DP) , allocatable, target :: & 69 dvscfin (:,:,:) ! change of the scf potential (input) 70 complex(DP) , pointer :: & 71 dvscfins (:,:,:) ! change of the scf potential (smooth) 72 complex(DP) , allocatable :: & 73 dvscfout (:,:,:), & ! change of the scf potential (output) 74 dbecsum(:,:,:,:), & ! the becsum with dpsi 75 auxg (:), aux1 (:), ps (:,:) 76 77 logical :: conv_root, exst 78 ! conv_root: true if linear system is converged 79 80 integer :: npw, npwq 81 integer :: kter, iter0, ipol, ibnd, jbnd, iter, lter, & 82 ik, ig, irr, ir, is, nrec, ios 83 ! counters 84 integer :: ltaver, lintercall 85 86 real(DP) :: tcpu 87 real(DP) :: eprec1 ! 1.35<ek>, for preconditioning 88 real(DP) :: iw !frequency 89 real(dp), external :: ddot, get_clock 90 91 external cch_psi_all, ccg_psi 92 93 if (lsda) call errore ('solve_e_fpol', ' LSDA not implemented', 1) 94 95 call start_clock ('solve_e') 96 allocate (dvscfin( dfftp%nnr, nspin, 3)) 97 if (doublegrid) then 98 allocate (dvscfins( dffts%nnr, nspin, 3)) 99 else 100 dvscfins => dvscfin 101 endif 102 allocate (dvscfout( dfftp%nnr, nspin, 3)) 103 allocate (dbecsum( nhm*(nhm+1)/2, nat, nspin, 3)) 104 allocate (auxg(npwx)) 105 allocate (aux1(dffts%nnr)) 106 allocate (ps (nbnd,nbnd)) 107 ps (:,:) = (0.d0, 0.d0) 108 allocate (h_diag(npwx, nbnd)) 109 110 allocate (etc(nbnd, nkstot)) 111 etc(:,:) = CMPLX( et(:,:), iw ,kind=DP) 112 113 ! restart NOT IMPLEMENTED 114 115 if (rec_code == -20) then 116 !read (iunrec) iter0, convt, dr2 117 !read (iunrec) dvscfin 118 !if (okvan) read (iunrec) int3 119 !close (unit = iunrec, status = 'keep') 120 !if (doublegrid) then 121 ! do is=1,nspin 122 ! do ipol=1,3 123 ! call fft_interpolate (dfftp, dvscfin(:,is,ipol), dffts, dvscfins(:,is,ipol)) 124 ! enddo 125 ! enddo 126 !endif 127 else if (rec_code > -20 .AND. rec_code <= -10) then 128 ! restarting in Raman: proceed 129 convt = .true. 130 else 131 convt = .false. 132 iter0 = 0 133 endif 134 ! 135 IF (ionode .AND. fildrho /= ' ') THEN 136 INQUIRE (UNIT = iudrho, OPENED = exst) 137 IF (exst) CLOSE (UNIT = iudrho, STATUS='keep') 138 CALL diropn (iudrho, TRIM(fildrho)//'.E', lrdrho, exst) 139 end if 140 ! 141 if (convt) go to 155 142 ! 143 ! if q=0 for a metal: allocate and compute local DOS at Ef 144 ! 145 if ( (lgauss .or. ltetra) .or. .not.lgamma) call errore ('solve_e_fpol', & 146 'called in the wrong case', 1) 147 ! 148 ! The outside loop is over the iterations 149 ! 150 do kter = 1, niter_ph 151 152 iter = kter + iter0 153 ltaver = 0 154 lintercall = 0 155 156 dvscfout(:,:,:)=(0.d0,0.d0) 157 dbecsum(:,:,:,:)=(0.d0,0.d0) 158 159 do ik = 1, nksq 160 if (lsda) current_spin = isk (ik) 161 npw = ngk(ik) 162 npwq = npw ! q=0 in ths routine 163 ! 164 ! read unperturbed wavefunctions psi_k in G_space, for all bands 165 ! 166 if (nksq.gt.1) call get_buffer(evc, lrwfc, iuwfc, ik) 167 ! 168 ! compute beta functions and kinetic energy for k-point ik 169 ! needed by h_psi, called by cch_psi_all, called by gmressolve_all 170 ! 171 CALL init_us_2 (npw, igk_k(1,ik), xk (1, ik), vkb) 172 CALL g2_kin(ik) 173 ! 174 do ipol = 1, 3 175 ! 176 ! computes/reads P_c^+ x psi_kpoint into dvpsi array 177 ! 178 call dvpsi_e (ik, ipol) 179 ! 180 if (iter > 1) then 181 ! 182 ! calculates dvscf_q*psi_k in G_space, for all bands, k=kpoint 183 ! dvscf_q from previous iteration (mix_potential) 184 ! 185 do ibnd = 1, nbnd_occ (ik) 186 aux1(:) = (0.d0, 0.d0) 187 do ig = 1, npw 188 aux1 (dffts%nl(igk_k(ig,ik)))=evc(ig,ibnd) 189 enddo 190 CALL invfft ('Wave', aux1, dffts) 191 do ir = 1, dffts%nnr 192 aux1(ir)=aux1(ir)*dvscfins(ir,current_spin,ipol) 193 enddo 194 CALL fwfft ('Wave', aux1, dffts) 195 do ig = 1, npwq 196 dvpsi(ig,ibnd)=dvpsi(ig,ibnd)+aux1(dffts%nl(igk_k(ig,ik))) 197 enddo 198 enddo 199 ! 200 call adddvscf(ipol,ik) 201 ! 202 endif 203 ! 204 ! Orthogonalize dvpsi to valence states: ps = <evc|dvpsi> 205 ! 206 CALL zgemm( 'C', 'N', nbnd_occ (ik), nbnd_occ (ik), npw, & 207 (1.d0,0.d0), evc(1,1), npwx, dvpsi(1,1), npwx, (0.d0,0.d0), & 208 ps(1,1), nbnd ) 209 210 call mp_sum ( ps( :, 1:nbnd_occ(ik) ), intra_bgrp_comm ) 211 ! dpsi is used as work space to store S|evc> 212 ! 213 CALL calbec (npw, vkb, evc, becp, nbnd_occ(ik) ) 214 CALL s_psi (npwx, npw, nbnd_occ(ik), evc, dpsi) 215 ! 216 ! |dvpsi> = - (|dvpsi> - S|evc><evc|dvpsi>) 217 ! note the change of sign! 218 ! 219 CALL zgemm( 'N', 'N', npw, nbnd_occ(ik), nbnd_occ(ik), & 220 (1.d0,0.d0), dpsi(1,1), npwx, ps(1,1), nbnd, (-1.d0,0.d0), & 221 dvpsi(1,1), npwx ) 222 ! 223 if (iter == 1) then 224 ! 225 ! At the first iteration dpsi and dvscfin are set to zero, 226 ! 227 dpsi(:,:)=(0.d0,0.d0) 228 dvscfin(:,:,:)=(0.d0,0.d0) 229 ! 230 ! starting threshold for the iterative solution of the linear 231 ! system 232 ! 233 thresh = 1.d-2 234 else 235 ! starting value for delta_psi is read from iudwf 236 ! 237 nrec = (ipol - 1) * nksq + ik 238 call get_buffer(dpsi, lrdwf, iudwf, nrec) 239 ! 240 ! threshold for iterative solution of the linear system 241 ! 242 thresh = min (0.1d0 * sqrt (dr2), 1.0d-2) 243 endif 244 ! 245 ! iterative solution of the linear system (H-e)*dpsi=dvpsi 246 ! dvpsi=-P_c+ (dvbare+dvscf)*psi , dvscf fixed. 247 ! 248 do ibnd = 1, nbnd_occ (ik) 249 ! 250 if ( (abs(iw).lt.0.05) .or. (abs(iw).gt.1.d0) ) then 251 ! 252 DO ig = 1, npwq 253 auxg (ig) = g2kin (ig) * evc (ig, ibnd) 254 END DO 255 eprec1 = 1.35_dp*ddot(2*npwq,evc(1,ibnd),1,auxg,1) 256 ! 257 do ig = 1, npw 258 h_diag(ig,ibnd)=CMPLX(1.d0, 0.d0,kind=DP) / & 259 CMPLX( max(1.0d0,g2kin(ig)/eprec1)-et(ibnd,ik),-iw ,kind=DP) 260 end do 261 else 262 do ig = 1, npw 263 h_diag(ig,ibnd)=CMPLX(1.d0, 0.d0,kind=DP) 264 end do 265 endif 266 ! 267 enddo 268 269 conv_root = .true. 270 271 call gmressolve_all (cch_psi_all,ccg_psi,etc(1,ik),dvpsi,dpsi, & 272 h_diag,npwx,npw,thresh,ik,lter,conv_root,anorm,nbnd_occ(ik), 4 ) 273 274 ltaver = ltaver + lter 275 lintercall = lintercall + 1 276 if (.not.conv_root) WRITE( stdout, "(5x,'kpoint',i4,' ibnd',i4, & 277 & ' solve_e: root not converged ',es10.3)") ik & 278 &, ibnd, anorm 279 ! 280 ! writes delta_psi on iunit iudwf, k=kpoint, 281 ! 282 nrec = (ipol - 1) * nksq + ik 283 call save_buffer (dpsi, lrdwf, iudwf, nrec) 284 ! 285 ! calculates dvscf, sum over k => dvscf_q_ipert 286 ! 287 call incdrhoscf (dvscfout(1,current_spin,ipol), wk(ik), & 288 ik, dbecsum(1,1,current_spin,ipol), dpsi) 289 enddo ! on polarizations 290 enddo ! on k points 291 ! 292 ! The calculation of dbecsum is distributed across processors 293 ! (see addusdbec) - we sum over processors the contributions 294 ! coming from each slice of bands 295 ! 296 call mp_sum ( dbecsum, intra_bgrp_comm ) 297 298 if (doublegrid) then 299 do is=1,nspin 300 do ipol=1,3 301 call fft_interpolate (dffts, dvscfout(:,is,ipol), dfftp, dvscfout(:,is,ipol)) 302 enddo 303 enddo 304 endif 305 306 call addusddense (dvscfout, dbecsum) 307 ! 308 ! dvscfout contains the (unsymmetrized) linear charge response 309 ! for the three polarizations - symmetrize it 310 ! 311 call mp_sum ( dvscfout, inter_pool_comm ) 312 call psyme (dvscfout) 313 ! 314 ! save the symmetrized linear charge response to file 315 ! calculate the corresponding linear potential response 316 ! 317 do ipol=1,3 318 if (fildrho.ne.' ') call davcio_drho(dvscfout(1,1,ipol),lrdrho, & 319 iudrho,ipol,+1) 320 call dv_of_drho (dvscfout (1, 1, ipol), .false.) 321 enddo 322 ! 323 ! mix the new potential with the old 324 ! 325 call mix_potential (2 * 3 * dfftp%nnr *nspin, dvscfout, dvscfin, alpha_mix ( & 326 kter), dr2, 3 * tr2_ph, iter, nmix_ph, flmixdpot, convt) 327 if (doublegrid) then 328 do is=1,nspin 329 do ipol = 1, 3 330 call fft_interpolate (dfftp,dvscfin(:,is,ipol),dffts,dvscfins(:,is,ipol)) 331 enddo 332 enddo 333 endif 334 335 call newdq(dvscfin,3) 336 337 averlt = DBLE (ltaver) / DBLE (lintercall) 338 339 tcpu = get_clock ('PHONON') 340 WRITE( stdout, '(/,5x," iter # ",i3," total cpu time :",f8.1, & 341 & " secs av.it.: ",f5.1)') iter, tcpu, averlt 342 dr2 = dr2 / 3 343 WRITE( stdout, "(5x,' thresh=',es10.3, ' alpha_mix = ',f6.3, & 344 & ' |ddv_scf|^2 = ',es10.3 )") thresh, alpha_mix (kter), dr2 345 ! 346 FLUSH( stdout ) 347 ! 348 ! restart NOT IMPLEMENTED 349 ! 350 !call seqopn (iunrec, 'recover', 'unformatted', exst) 351 ! 352 ! irr: state of the calculation 353 ! irr=-20 Electric Field 354 ! 355 !irr = -20 356 ! 357 !write (iunrec) irr 358 ! 359 ! partially calculated results 360 ! 361 !write (iunrec) dyn, dyn00 362 !write (iunrec) epsilon, zstareu, zstarue, zstareu0, zstarue0 363 ! 364 ! info on current iteration (iter=0 if potential mixing not available) 365 ! 366 !if (reduce_io) then 367 ! write (iunrec) 0, convt, dr2 368 !else 369 ! write (iunrec) iter, convt, dr2 370 !end if 371 !write (iunrec) dvscfin 372 !if (okvan) write (iunrec) int3 373 374 !close (unit = iunrec, status = 'keep') 375 if (check_stop_now()) then 376 call stop_smoothly_ph (.false.) 377 goto 155 378 endif 379 if (convt) goto 155 380 enddo 381155 continue 382 deallocate (h_diag) 383 deallocate (ps) 384 deallocate (aux1) 385 deallocate (auxg) 386 deallocate (dbecsum) 387 deallocate (dvscfout) 388 if (doublegrid) deallocate (dvscfins) 389 deallocate (dvscfin) 390 deallocate(etc) 391 392 call stop_clock ('solve_e') 393 return 394end subroutine solve_e_fpol 395