1! 2! Copyright (C) 2001-2020 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! 8MODULE lr_sternheimer 9 ! 10 ! This routine generalizes to finite complex frequencies and 11 ! finite q vectors the routine solve_e of the Quantum ESPRESSO 12 ! distribution. 13 ! 14 ! This routine is a driver for the solution of the linear system which 15 ! defines the change of the wavefunction due to an electric field 16 ! of finite wavevector q and complex frequency omega. 17 ! It performs the following tasks: 18 ! a) computes the bare potential term e^{iqr} | psi > 19 ! b) adds to it the screening term Delta V_{SCF} | psi > 20 ! c) applies P_c^+ (orthogonalization to valence states) 21 ! d) calls cgsolve_all to solve the linear system at zero 22 ! frequency or ccg_many_vectors 23 ! e) computes Delta rho, Delta V_{SCF} and symmetrizes them 24 ! 25CONTAINS 26 27SUBROUTINE one_sternheimer_step(iu, flag) 28 ! 29 USE kinds, ONLY : DP 30 USE constants, ONLY : e2, fpi, rytoev 31 USE ions_base, ONLY : nat 32 USE io_global, ONLY : stdout, ionode 33 USE io_files, ONLY : diropn, nwordwfc, iunwfc 34 USE cell_base, ONLY : tpiba2 35 USE fft_interfaces, ONLY : fwfft 36 USE fft_interfaces, ONLY : fft_interpolate 37 USE klist, ONLY : lgauss, xk, wk 38 USE gvecs, ONLY : doublegrid 39 USE fft_base, ONLY : dfftp, dffts 40 USE lsda_mod, ONLY : lsda, nspin, current_spin, isk 41 USE spin_orb, ONLY : domag 42 USE wvfct, ONLY : nbnd, npwx, g2kin, et 43 USE klist, ONLY : ngk, igk_k 44 USE check_stop, ONLY : check_stop_now 45 USE buffers, ONLY : get_buffer, save_buffer 46 USE wavefunctions, ONLY : evc 47 USE uspp, ONLY : okvan, vkb 48 USE uspp_param, ONLY : nhm 49 USE noncollin_module, ONLY : noncolin, npol, nspin_mag 50 USE scf, ONLY : rho, v_of_0 51 USE gvect, ONLY : gg 52 USE paw_variables, ONLY : okpaw 53 USE paw_onecenter, ONLY : paw_dpotential 54 USE eqv, ONLY : dpsi, dvpsi, evq 55 USE units_lr, ONLY : lrwfc, iuwfc 56 USE control_lr, ONLY : lgamma, alpha_pv, nbnd_occ, & 57 ext_recover, rec_code, & 58 lnoloc, convt, tr2_ph, & 59 alpha_mix, lgamma_gamma, niter_ph, & 60 flmixdpot, rec_code_read 61 62 USE lrus, ONLY : int3_paw, intq, intq_nc 63 USE qpoint, ONLY : xq, nksq, ikks, ikqs 64 USE linear_solvers, ONLY : ccg_many_vectors 65 USE dv_of_drho_lr, ONLY : dv_of_drho 66 USE mp_pools, ONLY : inter_pool_comm 67 USE mp_bands, ONLY : intra_bgrp_comm 68 USE mp_images, ONLY : root_image, my_image_id 69 USE mp, ONLY : mp_sum 70 USE fft_helper_subroutines, ONLY : fftx_ntgrp 71 USE lr_variables, ONLY : fru, fiu, iundvpsi, iudwf, & 72 lrdrho, iudrho, n_ipol, lr_verbosity, & 73 chirr, chirz, chizr, chizz, epsm1, & 74 current_w, lr1dwf, iu1dwf, itermax!, & 75 !intq, intq_nc 76 USE paw_add_symmetry, ONLY : paw_deqsymmetrize 77 USE wavefunctions, ONLY : psic 78 USE lr_sym_mod, ONLY : psymeq 79 ! 80 IMPLICIT NONE 81 ! 82 INTEGER, INTENT(IN) :: iu 83 INTEGER, INTENT(IN) :: flag ! if 1 compute the charge-charge and 84 ! charge magnetization responses 85 ! if 2 and lsda computes the magnetization 86 ! magnetization response 87 REAL(DP) :: thresh, anorm, averlt, dr2 88 ! thresh: convergence threshold 89 ! anorm : the norm of the error 90 ! averlt: average number of iterations 91 ! dr2 : self-consistency error 92 COMPLEX(DP), ALLOCATABLE :: h_diag (:,:) 93 COMPLEX(DP), ALLOCATABLE :: h_diag1 (:,:) 94 REAL(DP), ALLOCATABLE :: h_diagr (:,:) 95 REAL(DP), ALLOCATABLE :: h_dia (:,:), s_dia(:,:) 96 ! h_diag: diagonal part of the Hamiltonian 97 ! 98 COMPLEX(DP) , ALLOCATABLE, TARGET :: & 99 dvscfin (:,:,:) ! change of the scf potential (input) 100 COMPLEX(DP) , POINTER :: & 101 dvscfins (:,:,:) ! change of the scf potential (smooth) 102 COMPLEX(DP) , ALLOCATABLE :: & 103 dpsi1(:,:), & 104 dvscfout (:,:,:), & ! change of the scf potential (output) 105 drhoscfout (:,:), & ! change of the scf charge (output) 106 dbecsum(:,:,:,:), & ! the becsum with dpsi 107 dbecsum_nc(:,:,:,:,:), & ! the becsum with dpsi 108 mixin(:), mixout(:), & ! auxiliary for paw mixing 109 aux1 (:,:), ps (:,:), & 110 tg_dv(:,:), & 111 tg_psic(:,:), aux2(:,:), dvpsi1(:,:) 112 ! 113 COMPLEX(DP), EXTERNAL :: zdotc ! the scalar product function 114 ! 115 LOGICAL :: conv_root, exst, all_done_asyn 116 ! conv_root: true if linear system is converged 117 INTEGER :: kter, iter0, ipol, ibnd, iter, lter, ik, ikk, ikq, & 118 ig, is, nrec, ndim, npw, npwq, ios 119 ! counters 120 INTEGER :: ltaver, lintercall, incr, jpol, v_siz, nwordd0psi 121 REAL(DP) :: xqmod2, alpha_pv0 122 ! 123 REAL(DP) :: tcpu, get_clock 124 ! timing variables 125 ! 126 COMPLEX(DP) :: w !frequency 127 REAL(DP) :: aa, weight 128 LOGICAL :: ldpsi1 129 ! 130 EXTERNAL ch_psi_all, cg_psi 131 EXTERNAL ch_psi_all_complex, ccg_psi 132 COMPLEX(DP) :: scal_prod 133 ! 134 135 CALL start_clock ('stern_step') 136 137 w=CMPLX(fru(iu),fiu(iu)) 138 ldpsi1=ABS(w)>1.D-7 139 alpha_pv0=alpha_pv 140 alpha_pv=alpha_pv0 + REAL(w) 141 ! 142 ALLOCATE (dvscfin( dfftp%nnr, nspin_mag, 1)) 143 IF (doublegrid) THEN 144 ALLOCATE (dvscfins(dffts%nnr, nspin_mag, 1)) 145 ELSE 146 dvscfins => dvscfin 147 ENDIF 148 ALLOCATE (dvscfout(dfftp%nnr, nspin_mag, 1)) 149 ALLOCATE (drhoscfout(dfftp%nnr, nspin_mag)) 150 IF (okpaw) THEN 151 ALLOCATE (mixin(dfftp%nnr*nspin_mag+(nhm*(nhm+1)*nat*nspin_mag)/2) ) 152 ALLOCATE (mixout(dfftp%nnr*nspin_mag+(nhm*(nhm+1)*nat*nspin_mag)/2) ) 153 ENDIF 154 ALLOCATE (dbecsum( nhm*(nhm+1)/2, nat, nspin_mag, 1)) 155 IF (noncolin) ALLOCATE (dbecsum_nc (nhm, nhm, nat, nspin, 1)) 156 IF (ldpsi1) THEN 157 ALLOCATE (dpsi1(npwx*npol,nbnd)) 158 ALLOCATE (dvpsi1(npwx*npol,nbnd)) 159 ALLOCATE (h_diag(npwx*npol, nbnd)) 160 ALLOCATE (h_diag1(npwx*npol, nbnd)) 161 ALLOCATE (h_dia(npwx,npol)) 162 ALLOCATE (s_dia(npwx,npol)) 163 ELSE 164 ALLOCATE (h_diagr(npwx*npol, nbnd)) 165 ENDIF 166 ALLOCATE (aux1(dffts%nnr,npol)) 167 ALLOCATE (aux2(npwx*npol, nbnd)) 168 IF (okpaw) mixin=(0.0_DP,0.0_DP) 169 ! 170 171dvpsi =(0.0d0, 0.0d0) 172 173! IF (rec_code_read == -20.AND.ext_recover) then 174! ! restarting in Electric field calculation 175! IF (okpaw) THEN 176! CALL read_rec(dr2, iter0, 1, dvscfin, dvscfins, dvscfout, dbecsum) 177! CALL setmixout(3*dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag*3)/2, & 178! mixin, dvscfin, dbecsum, ndim, -1 ) 179! ELSE 180! CALL read_rec(dr2, iter0, 1, dvscfin, dvscfins) 181! ENDIF 182! ELSEIF (rec_code_read > -20 .AND. rec_code_read <= -10) then 183! ! restarting in Raman: proceed 184! convt = .true. 185! ELSE 186 convt = .false. 187 iter0 = 0 188! ENDIF 189 ! 190 incr=1 191 IF ( dffts%has_task_groups ) THEN 192 ! 193 v_siz = dffts%nnr_tg 194 ALLOCATE( tg_dv ( v_siz, nspin_mag ) ) 195 ALLOCATE( tg_psic( v_siz, npol ) ) 196 incr = fftx_ntgrp(dffts) 197 ! 198 ENDIF 199 ! 200! IF ( ionode .AND. fildrho /= ' ') THEN 201! INQUIRE (UNIT = iudrho, OPENED = exst) 202! IF (exst) CLOSE (UNIT = iudrho, STATUS='keep') 203! CALL diropn (iudrho, TRIM(fildrho)//'.E', lrdrho, exst) 204! ENDIF 205 IF (rec_code_read > -20) convt=.TRUE. 206 ! 207 IF (convt) go to 155 208 ! 209 IF ((lgauss.and..not.ldpsi1)) & 210 CALL errore ('solve_eq', 'insert a finite frequency', 1) 211 ! 212 IF (lr_verbosity > 5) THEN 213 WRITE(stdout,'("<lr_sternheimer_one_step>")') 214 ENDIF 215 ! 216 IF (.NOT. ALLOCATED(psic)) ALLOCATE(psic(dfftp%nnr)) 217 218! nwordd0psi = 2 * nbnd * npwx * npol * nksq 219! CALL diropn ( iund0psi, 'dvpsi.', nwordd0psi, exst) 220! 221! CALL diropn ( iudwf, 'dwf', nwordd0psi, exst) 222! CALL diropn ( iu1dwf, 'mwf', nwordd0psi, exst) 223 ! 224 ! Loop over the iterations 225 ! 226 DO kter = 1, itermax 227 ! 228 FLUSH( stdout) 229 iter = kter + iter0 230 231 ltaver = 0 232 lintercall = 0 233 ! 234 dvscfout(:,:,:) = (0.0d0, 0.0d0) 235 dbecsum(:,:,:,:) = (0.0d0, 0.0d0) 236 ! 237 IF (noncolin) dbecsum_nc = (0.0d0, 0.0d0) 238 ! 239 DO ik = 1, nksq 240 ! 241 ikk = ikks(ik) 242 ikq = ikqs(ik) 243 npw = ngk(ikk) 244 npwq = ngk(ikq) 245 if (lsda) current_spin = isk (ikk) 246 ! 247 ! Calculate beta-functions vkb at k+q (Kleinman-Bylander projectors) 248 ! The vkb's are needed for the non-local potential in h_psi, 249 ! and for the ultrasoft term. 250 ! 251 CALL init_us_2 (npwq, igk_k(1,ikq), xk(1,ikq), vkb) 252 ! 253 ! Read unperturbed wavefuctions evc (wfct at k) 254 ! and evq (wfct at k+q) 255 ! 256 IF (nksq > 1) THEN 257 CALL get_buffer (evc, nwordwfc, iunwfc, ikk) 258 CALL get_buffer (evq, nwordwfc, iunwfc, ikq) 259 ENDIF 260 ! 261 ! Compute the kinetic energy g2kin: (k+q+G)^2 262 ! 263 CALL g2_kin(ikq) 264 ! 265 ! IF omega non zero 266 ! 267 IF (ldpsi1) THEN 268 h_diag=(0.0_DP,0.0_DP) 269 h_diag1=(0.0_DP,0.0_DP) 270 h_dia=0.0_DP 271 s_dia=0.0_DP 272 CALL usnldiag( npwq, h_dia, s_dia ) 273 ! 274 275 DO ibnd = 1, nbnd_occ (ikk) 276 ! 277 DO ig = 1, npwq 278 aa=g2kin(ig)+v_of_0+h_dia(ig,1)- & 279 (et(ibnd,ikk)+w)*s_dia(ig,1) 280 IF (ABS(aa)<1.0_DP) aa=1.0_DP 281 h_diag(ig,ibnd)=CMPLX(1.0d0, 0.d0,kind=DP) / aa 282 aa=g2kin(ig)+v_of_0+h_dia(ig,1)- & 283 (et(ibnd,ikk)-w)*s_dia(ig,1) 284 IF (ABS(aa)<1.0_DP) aa=1.0_DP 285 h_diag1(ig,ibnd)=CMPLX(1.0d0, 0.d0,kind=DP) / aa 286 ENDDO 287 ! 288 IF (noncolin) THEN 289 DO ig = 1, npwq 290 aa=g2kin(ig)+v_of_0+h_dia(ig,2)- & 291 (et(ibnd,ikk)+w)*s_dia(ig,2) 292 IF (ABS(aa)<1.0_DP) aa=1.0_DP 293 h_diag(ig+npwx,ibnd)=CMPLX(1.d0, 0.d0,kind=DP) / aa 294 aa=g2kin(ig)+v_of_0+h_dia(ig,2)- & 295 (et(ibnd,ikk)-w)*s_dia(ig,2) 296 IF (ABS(aa)<1.0_DP) aa=1.0_DP 297 h_diag1(ig+npwx,ibnd)=CMPLX(1.d0, 0.d0,kind=DP) / aa 298 ENDDO 299 ENDIF 300 ENDDO 301 ELSE 302 CALL h_prec (ik, evc, h_diagr) 303 ! 304 DO ibnd = 1, nbnd_occ (ikk) 305 ! 306 DO ig = 1, npwq 307 aa=1.0_DP / h_diagr(ig,ibnd)-et(ibnd,ikk)-REAL(w,KIND=DP) 308 h_diagr(ig,ibnd)=1.d0 /max(1.0d0,aa) 309 ENDDO 310 IF (noncolin) THEN 311 DO ig = 1, npwq 312 h_diagr(ig+npwx,ibnd)= h_diagr(ig,ibnd) 313 ENDDO 314 ENDIF 315 ENDDO 316 ENDIF 317 ! 318 ! do over polarization 319 ! 320 DO ipol = 1, n_ipol 321 ! 322 nrec = (ipol - 1) * nksq + ik 323 ! 324 IF (iter ==1) THEN 325 ! 326 ! At the first iteration dvbare_q*psi_kpoint is calculated 327 ! and written to file 328 ! 329 CALL dveqpsi_us (ik) 330 ! 331 ! with flag=2 the perturbation is a magnetic field along z 332 ! 333 IF (lsda.AND.current_spin==2.AND.flag==2) dvpsi=-dvpsi 334! call save_buffer (dvpsi, lrbar, iubar, nrec) 335 336 CALL save_buffer (dvpsi, nwordwfc, iundvpsi, nrec) 337 ! 338 ELSE 339 ! 340 ! After the first iteration dvbare_q*psi_kpoint is read from file 341 ! 342! call get_buffer (dvpsi, lrbar, iubar, nrec) 343 CALL get_buffer (dvpsi, nwordwfc, iundvpsi, nrec) 344 ! 345 ! calculates dvscf_q*psi_k in G_space, for all bands, k=kpoint 346 ! dvscf_q from previous iteration (mix_potential) 347 ! 348 IF ( dffts%has_task_groups ) THEN 349 IF (noncolin) THEN 350 CALL tg_cgather( dffts, dvscfins(:,1,ipol), & 351 tg_dv(:,1)) 352 IF (domag) THEN 353 DO jpol=2,4 354 CALL tg_cgather( dffts, dvscfins(:,jpol,ipol), & 355 tg_dv(:,jpol)) 356 ENDDO 357 ENDIF 358 ELSE 359 CALL tg_cgather( dffts, dvscfins(:,current_spin,ipol), & 360 tg_dv(:,1)) 361 ENDIF 362 ENDIF 363 ! 364 aux2=(0.0_DP,0.0_DP) 365 DO ibnd = 1, nbnd_occ (ikk), incr 366 IF ( dffts%has_task_groups ) THEN 367 CALL cft_wave_tg (ik, evc, tg_psic, 1, v_siz, ibnd, & 368 nbnd_occ (ikk) ) 369 CALL apply_dpot(v_siz, tg_psic, tg_dv, 1) 370 CALL cft_wave_tg (ik, aux2, tg_psic, -1, v_siz, ibnd, & 371 nbnd_occ (ikk)) 372 ELSE 373 CALL cft_wave (ik, evc (1, ibnd), aux1, +1) 374 CALL apply_dpot(dffts%nnr, aux1, dvscfins(1,1,ipol), & 375 current_spin) 376 CALL cft_wave (ik, aux2 (1, ibnd), aux1, -1) 377 ENDIF 378 379 ENDDO 380 ! 381 dvpsi=dvpsi+aux2 382 ! 383 CALL adddvscf(ipol,ik) 384 ! 385 ENDIF 386 ! 387 ! Orthogonalize dvpsi to valence states: ps = <evq|dvpsi> 388 ! 389 IF (ldpsi1) THEN 390 dvpsi1=dvpsi 391 CALL orthogonalize_omega(dvpsi1, evq, ikk, ikq, dpsi, npwq, -w) 392 ENDIF 393 CALL orthogonalize_omega(dvpsi, evq, ikk, ikq, dpsi, npwq, w) 394 ! 395 IF (iter == 1) THEN 396 ! 397 ! At the first iteration dpsi and dvscfin are set to zero, 398 ! 399 dpsi(:,:)=(0.d0,0.d0) 400 IF (ldpsi1) dpsi1(:,:)=(0.d0,0.d0) 401 dvscfin(:,:,:)=(0.d0,0.d0) 402 ! 403 ! starting threshold for the iterative solution of the linear 404 ! system 405 ! 406 thresh = 1.d-2 407 IF (lnoloc) thresh = 1.d-5 408 ! 409 ELSE 410 ! 411 ! starting value for delta_psi is read from iudwf 412 ! 413 nrec = (ipol - 1) * nksq + ik 414 CALL get_buffer (dpsi, nwordwfc, iudwf, nrec) 415! call get_buffer (dpsi, lrdwf, iudwf, nrec) 416 IF (ldpsi1) CALL get_buffer (dpsi1, nwordwfc, iu1dwf, nrec) 417 ! 418 ! threshold for iterative solution of the linear system 419 ! 420 thresh = min (0.1d0 * sqrt (dr2), thresh) 421 ! 422 ENDIF 423 ! 424 ! iterative solution of the linear system (H-e)*dpsi=dvpsi 425 ! dvpsi=-P_c+ (dvbare+dvscf)*psi , dvscf fixed. 426 ! 427 conv_root = .true. 428 ! 429 current_w=w 430 IF (ldpsi1) THEN 431 ! 432 ! Complex or imaginary frequency. Use bicojugate gradient. 433 ! 434 435 CALL ccgsolve_all (ch_psi_all_complex,ccg_psi,et(1,ikk),dvpsi,dpsi, & 436 h_diag,npwx,npwq,thresh,ik,lter,conv_root,anorm,& 437 nbnd_occ(ikk),npol,current_w) 438 439 ! 440 ELSE 441 ! 442 ! zero frequency. The standard QE solver 443 ! 444 CALL cgsolve_all (ch_psi_all,cg_psi,et(1,ikk),dvpsi,dpsi, & 445 h_diagr,npwx,npwq,thresh,ik,lter,conv_root,anorm,& 446 nbnd_occ(ikk),npol) 447 ! 448 ENDIF 449 ! 450 ltaver = ltaver + lter 451 lintercall = lintercall + 1 452 IF (.not.conv_root) WRITE( stdout, "(5x,'kpoint',i4,' ibnd',i4, & 453 & ' solve_e: root not converged ',es10.3)") ik & 454 &, ibnd, anorm 455 ! 456 ! writes delta_psi on iunit iudwf, k=kpoint, 457 ! 458 nrec = (ipol - 1) * nksq + ik 459 CALL save_buffer (dpsi, nwordwfc, iudwf, nrec) 460! call save_buffer(dpsi, lrdwf, iudwf, nrec) 461 ! 462 463 ! calculates dvscf, sum over k => dvscf_q_ipert 464 ! 465 weight=wk(ikk) 466 IF (ldpsi1) THEN 467 ! 468 ! complex frequency, two wavefunctions must be computed 469 ! 470 weight=wk(ikk)/2.0_DP 471 ! 472 ! In this case compute also the wavefunction at frequency -w. 473 ! 474 current_w=-w 475 476 CALL ccgsolve_all (ch_psi_all_complex,ccg_psi,et(1,ikk),dvpsi1,dpsi1, & 477 h_diag1,npwx,npwq,thresh,ik,lter,conv_root,anorm,& 478 nbnd_occ(ikk),npol,current_w) 479 480 ltaver = ltaver + lter 481 lintercall = lintercall + 1 482 IF (.not.conv_root) WRITE( stdout, "(5x,'kpoint',i4,' ibnd',i4, & 483 & ' solve_e: root not converged ',es10.3)") ik & 484 &, ibnd, anorm 485 ! 486 ! writes delta_psi on iunit iudwf, k=kpoint, 487 ! 488 nrec = (ipol - 1) * nksq + ik 489 CALL save_buffer(dpsi1, nwordwfc, iu1dwf, nrec) 490 ! 491 ! calculates dvscf, sum over k => dvscf_q_ipert 492 ! 493 CALL DAXPY(npwx*nbnd_occ(ikk)*npol*2, 1.0_DP, dpsi1, 1, dpsi, 1) 494 ! 495 ENDIF 496 IF (noncolin) THEN 497 CALL incdrhoscf_nc(dvscfout(1,1,ipol), weight, ik, & 498 dbecsum_nc(1,1,1,1,ipol), dpsi, 1.0d0) 499 ELSE 500 501 CALL incdrhoscf (dvscfout(1,current_spin,ipol), weight, & 502 ik, dbecsum(1,1,current_spin,ipol), dpsi) 503 ENDIF 504 ! 505 ENDDO ! on polarizations 506 ! 507! IF ( with_asyn_images.AND.my_image_id==root_image.AND.ionode ) & 508! CALL asyn_master(all_done_asyn) 509 ENDDO ! on k points 510 ! 511 512 current_w=w 513 ! 514 ! The calculation of dbecsum is distributed across processors 515 ! (see addusdbec) - we sum over processors the contributions 516 ! coming from each slice of bands 517 ! 518 IF (noncolin) THEN 519 CALL mp_sum ( dbecsum_nc, intra_bgrp_comm ) 520 ELSE 521 CALL mp_sum ( dbecsum, intra_bgrp_comm ) 522 ENDIF 523 ! 524 IF (doublegrid) then 525 DO is=1,nspin_mag 526 CALL fft_interpolate (dffts, dvscfout(:,is,1), dfftp, & 527 dvscfout(:,is,1)) 528 ENDDO 529 ENDIF 530 ! 531 IF (noncolin.and.okvan) CALL set_dbecsum_nc(dbecsum_nc, dbecsum, 1) 532 ! 533 CALL addusddenseq (dvscfout, dbecsum) 534 ! 535 ! dvscfout contains the (unsymmetrized) linear charge response 536 ! for the three polarizations - symmetrize it 537 ! 538 CALL mp_sum ( dvscfout, inter_pool_comm ) 539 ! 540 IF (okpaw) CALL mp_sum ( dbecsum, inter_pool_comm ) 541 ! 542 IF (.not. lgamma_gamma) THEN 543 CALL psymeq (dvscfout) 544 ENDIF 545 ! 546 drhoscfout(:,:)=dvscfout(:,:,1) 547 ! 548 ! save the symmetrized linear charge response to file 549 ! calculate the corresponding linear potential response 550 ! 551 IF (lnoloc) THEN 552! CALL dv_of_drho_nlf (dvscfout (1, 1, 1)) 553 ELSE 554 CALL dv_of_drho (dvscfout (1, 1, 1), .FALSE.) 555 ENDIF 556 ! 557 ! mix the new potential with the old 558 ! 559 IF (okpaw) THEN 560 ! 561 ! In this case we mix also dbecsum 562 ! 563 CALL setmixout(dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag)/2, & 564 mixout, dvscfout, dbecsum, ndim, -1 ) 565 CALL mix_potential_eels (2*dfftp%nnr*nspin_mag+2*ndim, mixout, mixin, & 566 alpha_mix(kter), dr2, tr2_ph/npol, iter, flmixdpot, & 567 convt) 568 CALL setmixout(dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag)/2, & 569 mixin, dvscfin, dbecsum, ndim, 1 ) 570 ELSE 571 572 CALL mix_potential_eels(2*dfftp%nnr*nspin_mag, dvscfout, dvscfin, & 573 alpha_mix (kter), dr2, tr2_ph/npol, iter, flmixdpot, convt) 574 575 ENDIF 576 ! 577 IF (doublegrid) then 578 DO is=1,nspin_mag 579 CALL fft_interpolate (dfftp, dvscfin(:,is,1), dffts, & 580 dvscfins(:,is,1)) 581 ENDDO 582 ENDIF 583 ! 584 IF (okpaw) THEN 585 IF (noncolin.AND.domag) THEN 586! CALL PAW_dpotential(dbecsum_nc,becsum_nc,int3_paw,3) 587 ELSE 588 ! 589 ! The presence of c.c. in the formula gives a factor 2.0 590 ! 591 dbecsum=2.0_DP * dbecsum 592 IF (.NOT. lgamma_gamma) CALL paw_deqsymmetrize(dbecsum) 593 CALL PAW_dpotential(dbecsum,rho%bec,int3_paw,1) 594 ENDIF 595 ENDIF 596 ! 597 CALL newdq(dvscfin,1) 598 ! 599 CALL mp_sum(ltaver,inter_pool_comm) 600 CALL mp_sum(lintercall,inter_pool_comm) 601 ! 602 averlt = DBLE (ltaver) / DBLE (lintercall) 603 ! 604! tcpu = get_clock ('PHONON') 605 tcpu = get_clock ('ccgsolve') 606 WRITE( stdout, '(/,5x," iter # ",i3," total cpu time :",f8.1, & 607 & " secs av.it.: ",f5.1)') iter, tcpu, averlt 608 WRITE( stdout, "(5x,' thresh=',es10.3, ' alpha_mix = ',f6.3, & 609 & ' |ddv_scf|^2 = ',es10.3 )") thresh, alpha_mix (kter), dr2 610 ! 611 FLUSH( stdout ) 612 ! 613 ! rec_code: state of the calculation 614 ! rec_code=-20 Electric Field 615 ! 616 rec_code=-20 617! IF (okpaw) THEN 618! CALL write_rec('solve_e...', 0, dr2, iter, convt, 1, dvscfin, & 619! dvscfout, dbecsum) 620! ELSE 621! CALL write_rec('solve_e...', 0, dr2, iter, convt, 1, dvscfin) 622! ENDIF 623 ! 624! IF (check_stop_now()) CALL stop_smoothly_ph (.false.) 625 ! 626 IF (convt) goto 155 627 ! 628 ENDDO ! do over iter_ph 629 ! 630155 CONTINUE 631 ! 632 ! compute here the susceptibility and the inverse of the dielectric 633 ! constant 634 ! 635 ! CALL compute_susceptibility(drhoscfout) 636 ! 637! CLOSE( UNIT = iund0psi) 638 639 DO is=1,nspin_mag 640 CALL fwfft ('Rho', drhoscfout(:,is), dfftp) 641 ENDDO 642 ! 643 IF (flag==1) THEN 644 chirr(iu)=(0.0_DP,0.0_DP) 645 chizr(iu)=(0.0_DP,0.0_DP) 646 epsm1(iu)=(0.0_DP,0.0_DP) 647 ELSE 648 chirz(iu)=(0.0_DP,0.0_DP) 649 chizz(iu)=(0.0_DP,0.0_DP) 650 ENDIF 651 ! 652 xqmod2=(xq(1)**2+xq(2)**2+xq(3)**2)*tpiba2 653 ! 654 IF (ABS(gg(1))<1.d-8) THEN 655 IF (flag==1) THEN 656 chirr(iu) = drhoscfout(dfftp%nl(1),1) 657 IF (lsda) chirr(iu) = chirr(iu) + drhoscfout(dfftp%nl(1),2) 658 epsm1(iu) = CMPLX(1.0_DP,0.0_DP)+ chirr(iu)*fpi*e2/xqmod2 659 IF (lsda) chizr(iu) = drhoscfout(dfftp%nl(1),1) - & 660 drhoscfout(dfftp%nl(1),2) 661 ELSEIF (lsda) THEN 662 chizz(iu)=drhoscfout(dfftp%nl(1),1)-drhoscfout(dfftp%nl(1),2) 663 chirz(iu)=drhoscfout(dfftp%nl(1),1)+drhoscfout(dfftp%nl(1),2) 664 ENDIF 665 ENDIF 666 ! 667 IF (flag==1) THEN 668 CALL mp_sum(epsm1(iu),intra_bgrp_comm) 669 CALL mp_sum(chirr(iu),intra_bgrp_comm) 670 CALL mp_sum(chizr(iu),intra_bgrp_comm) 671 ELSE 672 CALL mp_sum(chizz(iu),intra_bgrp_comm) 673 CALL mp_sum(chirz(iu),intra_bgrp_comm) 674 ENDIF 675 ! 676 IF (flag==1) THEN 677 WRITE(stdout, '(/,6x,"Inverse dielectric constant at & 678 &frequency",f9.4," +",f9.4," i Ry")') fru(iu), fiu(iu) 679 WRITE(stdout, '(46x,f9.4," +",f9.4," i eV")') current_w * rytoev 680 WRITE(stdout,'(/,6x,"epsilon^-1(q,w) =",2f15.6)') epsm1(iu) 681 ! 682 WRITE( stdout, '(/,6x,"Charge-charge susceptibility:")') 683 ! 684 WRITE(stdout,'(/,6x,"chirr(q,w) =",2f15.6)') chirr(iu) 685 IF (lsda) THEN 686 WRITE(stdout,'(/,6x,"m_z-charge susceptibility:")') 687 WRITE(stdout,'(/,6x,"chizr(q,w) =",2f15.6)') chizr(iu) 688 ENDIF 689 ! 690 ELSEIF (lsda) THEN 691 WRITE( stdout, '(/,6x,"m_z - m_z susceptibility at & 692 &frequency",f9.4," +",f9.4," i Ry")') fru(iu), fiu(iu) 693 WRITE( stdout, '(43x,f9.4," +",f9.4," i eV")') current_w * rytoev 694 WRITE(stdout,'(/,6x,"chizz(q,w) =",2f15.6)') chizz(iu) 695 WRITE(stdout,'(/,6x,"chirz(q,w) =",2f15.6)') chirz(iu) 696 ENDIF 697 ! 698 deallocate (aux1) 699 IF (ldpsi1) THEN 700 deallocate (dpsi1) 701 deallocate (dvpsi1) 702 deallocate (h_diag) 703 deallocate (h_diag1) 704 deallocate (h_dia) 705 deallocate (s_dia) 706 ELSE 707 deallocate (h_diagr) 708 ENDIF 709 deallocate (dbecsum) 710 deallocate (dvscfout) 711 IF (okpaw) THEN 712 DEALLOCATE(mixin) 713 DEALLOCATE(mixout) 714 ENDIF 715 deallocate (drhoscfout) 716 if (doublegrid) deallocate (dvscfins) 717 deallocate (dvscfin) 718 if (noncolin) deallocate(dbecsum_nc) 719 deallocate(aux2) 720 IF ( dffts%has_task_groups ) THEN 721 ! 722 DEALLOCATE( tg_dv ) 723 DEALLOCATE( tg_psic) 724 ! 725 ENDIF 726 727! CLOSE( unit = iund0psi) 728! CLOSE( unit = iudwf) 729! CLOSE( unit = iu1dwf) 730 731 alpha_pv=alpha_pv0 732 ! 733 CALL stop_clock ('stern_step') 734 ! 735 RETURN 736 ! 737END SUBROUTINE one_sternheimer_step 738 739END MODULE lr_sternheimer 740