1! 2! Copyright (C) 2001-2018 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 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 ! If lda_plus_u=.true. compute also the SCF part 19 ! of the response Hubbard potential. 20 ! c) applies P_c^+ (orthogonalization to valence states) 21 ! d) calls cgsolve_all to solve the linear system 22 ! e) computes Delta rho, Delta V_{SCF} and symmetrizes them 23 ! f) If lda_plus_u=.true. compute also the response occupation 24 ! matrices dnsscf 25 ! 26 USE kinds, ONLY : DP 27 USE ions_base, ONLY : nat 28 USE io_global, ONLY : stdout, ionode 29 USE io_files, ONLY : prefix, diropn 30 USE cell_base, ONLY : tpiba2 31 USE klist, ONLY : ltetra, lgauss, xk, wk, ngk, igk_k 32 USE gvect, ONLY : g 33 USE gvecs, ONLY : doublegrid 34 USE fft_base, ONLY : dfftp, dffts 35 USE lsda_mod, ONLY : lsda, nspin, current_spin, isk 36 USE spin_orb, ONLY : domag 37 USE wvfct, ONLY : nbnd, npwx, g2kin, et 38 USE check_stop, ONLY : check_stop_now 39 USE buffers, ONLY : get_buffer, save_buffer 40 USE wavefunctions, ONLY : evc 41 USE uspp, ONLY : okvan, vkb 42 USE uspp_param, ONLY : upf, nhm 43 USE noncollin_module, ONLY : noncolin, npol, nspin_mag 44 USE scf, ONLY : rho 45 USE paw_variables, ONLY : okpaw 46 USE paw_onecenter, ONLY : paw_dpotential 47 USE paw_symmetry, ONLY : paw_desymmetrize 48 49 USE units_ph, ONLY : lrdwf, iudwf, lrdrho, iudrho 50 USE units_lr, ONLY : iuwfc, lrwfc 51 USE output, ONLY : fildrho 52 USE control_ph, ONLY : ext_recover, rec_code, & 53 lnoloc, convt, tr2_ph, nmix_ph, & 54 alpha_mix, lgamma_gamma, niter_ph, & 55 flmixdpot, rec_code_read 56 USE recover_mod, ONLY : read_rec, write_rec 57 USE mp_pools, ONLY : inter_pool_comm 58 USE mp_bands, ONLY : intra_bgrp_comm, ntask_groups 59 USE mp, ONLY : mp_sum 60 USE lrus, ONLY : int3_paw 61 USE qpoint, ONLY : nksq, ikks, ikqs 62 USE eqv, ONLY : dpsi, dvpsi 63 USE control_lr, ONLY : nbnd_occ, lgamma 64 USE dv_of_drho_lr 65 USE fft_helper_subroutines 66 USE fft_interfaces, ONLY : fft_interpolate 67 USE ldaU, ONLY : lda_plus_u 68 69 implicit none 70 71 real(DP) :: thresh, anorm, averlt, dr2 72 ! thresh: convergence threshold 73 ! anorm : the norm of the error 74 ! averlt: average number of iterations 75 ! dr2 : self-consistency error 76 real(DP), allocatable :: h_diag (:,:) 77 ! h_diag: diagonal part of the Hamiltonian 78 79 complex(DP) , allocatable, target :: & 80 dvscfin (:,:,:) ! change of the scf potential (input) 81 complex(DP) , pointer :: & 82 dvscfins (:,:,:) ! change of the scf potential (smooth) 83 complex(DP) , allocatable :: & 84 dvscfout (:,:,:), & ! change of the scf potential (output) 85 dbecsum(:,:,:,:), & ! the becsum with dpsi 86 dbecsum_nc(:,:,:,:,:), & ! the becsum with dpsi 87 mixin(:), mixout(:), & ! auxiliary for paw mixing 88 aux1 (:,:), ps (:,:), & 89 tg_dv(:,:), & 90 tg_psic(:,:), aux2(:,:) 91 92 logical :: conv_root, exst 93 ! conv_root: true if linear system is converged 94 95 integer :: npw, npwq 96 integer :: kter, iter0, ipol, ibnd, iter, lter, ik, ig, is, nrec, ndim, ios 97 ! counters 98 integer :: ltaver, lintercall, incr, jpol, v_siz 99 integer :: ikk, ikq 100 101 real(DP) :: tcpu, get_clock 102 ! timing variables 103 104 external ch_psi_all, cg_psi 105 106 call start_clock ('solve_e') 107 ! 108 ! This routine is task group aware 109 ! 110 allocate (dvscfin( dfftp%nnr, nspin_mag, 3)) 111 if (doublegrid) then 112 allocate (dvscfins(dffts%nnr, nspin_mag, 3)) 113 else 114 dvscfins => dvscfin 115 endif 116 allocate (dvscfout(dfftp%nnr, nspin_mag, 3)) 117 IF (okpaw) THEN 118 ALLOCATE (mixin(dfftp%nnr*nspin_mag*3+(nhm*(nhm+1)*nat*nspin_mag*3)/2) ) 119 ALLOCATE (mixout(dfftp%nnr*nspin_mag*3+(nhm*(nhm+1)*nat*nspin_mag*3)/2) ) 120 ENDIF 121 allocate (dbecsum( nhm*(nhm+1)/2, nat, nspin_mag, 3)) 122 IF (noncolin) allocate (dbecsum_nc (nhm, nhm, nat, nspin, 3)) 123 allocate (aux1(dffts%nnr,npol)) 124 allocate (h_diag(npwx*npol, nbnd)) 125 allocate (aux2(npwx*npol, nbnd)) 126 IF (okpaw) mixin=(0.0_DP,0.0_DP) 127 128 if (rec_code_read == -20.AND.ext_recover) then 129 ! restarting in Electric field calculation 130 IF (okpaw) THEN 131 CALL read_rec(dr2, iter0, 3, dvscfin, dvscfins, dvscfout, dbecsum) 132 CALL setmixout(3*dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag*3)/2, & 133 mixin, dvscfin, dbecsum, ndim, -1 ) 134 ELSE 135 CALL read_rec(dr2, iter0, 3, dvscfin, dvscfins) 136 ENDIF 137 else if (rec_code_read > -20 .AND. rec_code_read <= -10) then 138 ! restarting in Raman: proceed 139 convt = .true. 140 else 141 convt = .false. 142 iter0 = 0 143 endif 144 incr=1 145 IF ( dffts%has_task_groups ) THEN 146 ! 147 v_siz = dffts%nnr_tg 148 ALLOCATE( tg_dv ( v_siz, nspin_mag ) ) 149 ALLOCATE( tg_psic( v_siz, npol ) ) 150 incr = fftx_ntgrp(dffts) 151 ! 152 ENDIF 153 ! 154 IF ( ionode .AND. fildrho /= ' ') THEN 155 INQUIRE (UNIT = iudrho, OPENED = exst) 156 IF (exst) CLOSE (UNIT = iudrho, STATUS='keep') 157 CALL diropn (iudrho, TRIM(fildrho)//'.E', lrdrho, exst) 158 end if 159 IF (rec_code_read > -20) convt=.TRUE. 160 ! 161 if (convt) go to 155 162 ! 163 ! if q=0 for a metal: allocate and compute local DOS at Ef 164 ! 165 if ( (lgauss .or. ltetra) .or..not.lgamma) call errore ('solve_e', & 166 'called in the wrong case', 1) 167 ! 168 ! The outside loop is over the iterations 169 ! 170 do kter = 1, niter_ph 171 ! 172 FLUSH( stdout ) 173 iter = kter + iter0 174 ltaver = 0 175 lintercall = 0 176 ! 177 dvscfout(:,:,:)=(0.d0,0.d0) 178 dbecsum(:,:,:,:)=(0.d0,0.d0) 179 IF (noncolin) dbecsum_nc=(0.d0,0.d0) 180 ! 181 ! DFPT+U: at each iteration calculate dnsscf, 182 ! i.e. the scf variation of the occupation matrix ns. 183 ! 184 IF (lda_plus_u .AND. (iter.NE.1)) & 185 CALL dnsq_scf (3, .false., 0, 1, .false.) 186 ! 187 do ik = 1, nksq 188 ikk=ikks(ik) 189 ikq=ikqs(ik) 190 ! 191 npw = ngk(ikk) 192 npwq= npw ! q=0 always in this routine 193 ! 194 if (lsda) current_spin = isk (ikk) 195 ! 196 ! reads unperturbed wavefunctions psi_k in G_space, for all bands 197 ! 198 if (nksq.ge.1) THEN 199 call get_buffer (evc, lrwfc, iuwfc, ikk) 200 end if 201 ! 202 ! compute beta functions and kinetic energy for k-point ik 203 ! needed by h_psi, called by ch_psi_all, called by cgsolve_all 204 ! 205 CALL init_us_2 (npw, igk_k(1,ikk), xk (1, ikk), vkb) 206 CALL g2_kin(ikk) 207 ! 208 ! compute preconditioning matrix h_diag used by cgsolve_all 209 ! 210 CALL h_prec (ik, evc, h_diag) 211 ! 212 do ipol = 1, 3 213 ! 214 ! computes/reads P_c^+ x psi_kpoint into dvpsi array 215 ! 216 call dvpsi_e (ik, ipol) 217 ! 218 if (iter > 1) then 219 ! 220 ! calculates dvscf_q*psi_k in G_space, for all bands, k=kpoint 221 ! dvscf_q from previous iteration (mix_potential) 222 ! 223 IF( dffts%has_task_groups ) THEN 224 IF (noncolin) THEN 225 CALL tg_cgather( dffts, dvscfins(:,1,ipol), tg_dv(:,1)) 226 IF (domag) THEN 227 DO jpol=2,4 228 CALL tg_cgather( dffts, dvscfins(:,jpol,ipol), tg_dv(:,jpol)) 229 ENDDO 230 ENDIF 231 ELSE 232 CALL tg_cgather( dffts, dvscfins(:,current_spin,ipol), tg_dv(:,1)) 233 ENDIF 234 ENDIF 235 aux2=(0.0_DP,0.0_DP) 236 do ibnd = 1, nbnd_occ (ikk), incr 237 IF ( dffts%has_task_groups ) THEN 238 call cft_wave_tg (ik, evc, tg_psic, 1, v_siz, ibnd, nbnd_occ (ikk) ) 239 call apply_dpot(v_siz, tg_psic, tg_dv, 1) 240 call cft_wave_tg (ik, aux2, tg_psic, -1, v_siz, ibnd, nbnd_occ (ikk)) 241 ELSE 242 call cft_wave (ik, evc (1, ibnd), aux1, +1) 243 call apply_dpot(dffts%nnr, aux1, dvscfins(1,1,ipol), current_spin) 244 call cft_wave (ik, aux2 (1, ibnd), aux1, -1) 245 ENDIF 246 enddo 247 dvpsi=dvpsi+aux2 248 ! 249 ! In the case of US pseudopotentials there is an additional 250 ! selfconsist term which comes from the dependence of D on 251 ! V_{eff} on the bare change of the potential 252 ! 253 call adddvscf(ipol,ik) 254 ! 255 ! DFPT+U: add to dvpsi the scf part of the response 256 ! Hubbard potential dV_hub 257 ! 258 if (lda_plus_u) call adddvhubscf (ipol, ik) 259 ! 260 endif 261 ! 262 ! Orthogonalize dvpsi to valence states: ps = <evc|dvpsi> 263 ! 264 CALL orthogonalize(dvpsi, evc, ikk, ikk, dpsi, npwq, .false.) 265 ! 266 if (iter == 1) then 267 ! 268 ! At the first iteration dpsi and dvscfin are set to zero, 269 ! 270 dpsi(:,:)=(0.d0,0.d0) 271 dvscfin(:,:,:)=(0.d0,0.d0) 272 ! 273 ! starting threshold for the iterative solution of the linear 274 ! system 275 ! 276 thresh = 1.d-2 277 if (lnoloc) thresh = 1.d-5 278 else 279 ! starting value for delta_psi is read from iudwf 280 ! 281 nrec = (ipol - 1) * nksq + ik 282 call get_buffer (dpsi, lrdwf, iudwf, nrec) 283 ! 284 ! threshold for iterative solution of the linear system 285 ! 286 thresh = min (0.1d0 * sqrt (dr2), 1.0d-2) 287 endif 288 ! 289 ! iterative solution of the linear system (H-e)*dpsi=dvpsi 290 ! dvpsi=-P_c+ (dvbare+dvscf)*psi , dvscf fixed. 291 ! 292 293 conv_root = .true. 294 295 call cgsolve_all (ch_psi_all,cg_psi,et(1,ikk),dvpsi,dpsi, & 296 h_diag,npwx,npw,thresh,ik,lter,conv_root,anorm,nbnd_occ(ikk),npol) 297 298 ltaver = ltaver + lter 299 lintercall = lintercall + 1 300 if (.not.conv_root) WRITE( stdout, "(5x,'kpoint',i4,' ibnd',i4, & 301 & ' solve_e: root not converged ',es10.3)") ik & 302 &, ibnd, anorm 303 ! 304 ! writes delta_psi on iunit iudwf, k=kpoint, 305 ! 306 nrec = (ipol - 1) * nksq + ik 307 call save_buffer(dpsi, lrdwf, iudwf, nrec) 308 ! 309 ! calculates dvscf, sum over k => dvscf_q_ipert 310 ! 311 IF (noncolin) THEN 312 call incdrhoscf_nc(dvscfout(1,1,ipol),wk(ikk),ik, & 313 dbecsum_nc(1,1,1,1,ipol), dpsi, 1.0d0) 314 ELSE 315 call incdrhoscf (dvscfout(1,current_spin,ipol), wk(ikk), & 316 ik, dbecsum(1,1,current_spin,ipol), dpsi) 317 ENDIF 318 enddo ! on polarizations 319 enddo ! on k points 320 ! 321 ! The calculation of dbecsum is distributed across processors 322 ! (see addusdbec) - we sum over processors the contributions 323 ! coming from each slice of bands 324 ! 325 IF (noncolin) THEN 326 call mp_sum ( dbecsum_nc, intra_bgrp_comm ) 327 ELSE 328 call mp_sum ( dbecsum, intra_bgrp_comm ) 329 END IF 330 331 if (doublegrid) then 332 do is=1,nspin_mag 333 do ipol=1,3 334 call fft_interpolate (dffts, dvscfout(:,is,ipol), dfftp, dvscfout(:,is,ipol)) 335 enddo 336 enddo 337 endif 338 ! 339 IF (noncolin.and.okvan) CALL set_dbecsum_nc(dbecsum_nc, dbecsum, 3) 340 ! 341 call addusddense (dvscfout, dbecsum) 342 ! 343 ! dvscfout contains the (unsymmetrized) linear charge response 344 ! for the three polarizations - symmetrize it 345 ! 346 call mp_sum ( dvscfout, inter_pool_comm ) 347 IF (okpaw) call mp_sum ( dbecsum, inter_pool_comm ) 348 if (.not.lgamma_gamma) then 349 call psyme (dvscfout) 350 IF ( noncolin.and.domag ) CALL psym_dmage(dvscfout) 351 endif 352 ! 353 ! save the symmetrized linear charge response to file 354 ! calculate the corresponding linear potential response 355 ! 356 do ipol=1,3 357 if (fildrho.ne.' ') call davcio_drho(dvscfout(1,1,ipol),lrdrho, & 358 iudrho,ipol,+1) 359 IF (lnoloc) then 360 dvscfout(:,:,ipol)=(0.d0,0.d0) 361 ELSE 362 call dv_of_drho (dvscfout (1, 1, ipol), .false.) 363 ENDIF 364 enddo 365 ! 366 ! mix the new potential with the old 367 ! 368 IF (okpaw) THEN 369 ! 370 ! In this case we mix also dbecsum 371 ! 372 call setmixout(3*dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag*3)/2, & 373 mixout, dvscfout, dbecsum, ndim, -1 ) 374 call mix_potential (2*3*dfftp%nnr*nspin_mag+2*ndim, mixout, mixin, & 375 alpha_mix(kter), dr2, 3*tr2_ph/npol, iter, & 376 nmix_ph, flmixdpot, convt) 377 call setmixout(3*dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag*3)/2, & 378 mixin, dvscfin, dbecsum, ndim, 1 ) 379 ELSE 380 call mix_potential (2*3*dfftp%nnr*nspin_mag, dvscfout, dvscfin, alpha_mix ( & 381 kter), dr2, 3 * tr2_ph / npol, iter, nmix_ph, flmixdpot, convt) 382 ENDIF 383 if (doublegrid) then 384 do is=1,nspin_mag 385 do ipol = 1, 3 386 call fft_interpolate (dfftp, dvscfin(:,is,ipol), dffts, dvscfins(:,is,ipol)) 387 enddo 388 enddo 389 endif 390 391 IF (okpaw) THEN 392 IF (noncolin) THEN 393! call PAW_dpotential(dbecsum_nc,becsum_nc,int3_paw,3) 394 ELSE 395! 396! The presence of c.c. in the formula gives a factor 2.0 397! 398 dbecsum=2.0_DP * dbecsum 399 IF (.NOT. lgamma_gamma) CALL PAW_desymmetrize(dbecsum) 400 call PAW_dpotential(dbecsum,rho%bec,int3_paw,3) 401 ENDIF 402 ENDIF 403 404 call newdq(dvscfin,3) 405 406 averlt = DBLE (ltaver) / DBLE (lintercall) 407 408 tcpu = get_clock ('PHONON') 409 WRITE( stdout, '(/,5x," iter # ",i3," total cpu time :",f8.1, & 410 & " secs av.it.: ",f5.1)') iter, tcpu, averlt 411 dr2 = dr2 / 3 412 WRITE( stdout, "(5x,' thresh=',es10.3, ' alpha_mix = ',f6.3, & 413 & ' |ddv_scf|^2 = ',es10.3 )") thresh, alpha_mix (kter), dr2 414 ! 415 FLUSH( stdout ) 416 ! 417 ! rec_code: state of the calculation 418 ! rec_code=-20 Electric Field 419 ! 420 rec_code=-20 421 IF (okpaw) THEN 422 CALL write_rec('solve_e...', 0, dr2, iter, convt, 3, dvscfin, & 423 dvscfout, dbecsum) 424 ELSE 425 CALL write_rec('solve_e...', 0, dr2, iter, convt, 3, dvscfin) 426 ENDIF 427 428 if (check_stop_now()) call stop_smoothly_ph (.false.) 429 430 if (convt) goto 155 431 432 enddo 433155 continue 434 ! 435 deallocate (h_diag) 436 deallocate (aux1) 437 deallocate (dbecsum) 438 deallocate (dvscfout) 439 IF (okpaw) THEN 440 DEALLOCATE(mixin) 441 DEALLOCATE(mixout) 442 ENDIF 443 if (doublegrid) deallocate (dvscfins) 444 deallocate (dvscfin) 445 if (noncolin) deallocate(dbecsum_nc) 446 deallocate(aux2) 447 IF ( dffts%has_task_groups ) THEN 448 ! 449 DEALLOCATE( tg_dv ) 450 DEALLOCATE( tg_psic) 451 ! 452 ENDIF 453 454 call stop_clock ('solve_e') 455 return 456end subroutine solve_e 457