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_linter (irr, imode0, npe, drhoscf) 11 !----------------------------------------------------------------------- 12 ! 13 ! Driver routine for the solution of the linear system which 14 ! defines the change of the wavefunction due to a lattice distorsion 15 ! It performs the following tasks: 16 ! a) computes the bare potential term Delta V | psi > 17 ! and an additional term in the case of US pseudopotentials. 18 ! If lda_plus_u=.true. compute also the bare potential 19 ! term Delta V_hub | psi >. 20 ! b) adds to it the screening term Delta V_{SCF} | psi >. 21 ! If lda_plus_u=.true. compute also the SCF part 22 ! of the response Hubbard potential. 23 ! c) applies P_c^+ (orthogonalization to valence states) 24 ! d) calls cgsolve_all to solve the linear system 25 ! e) computes Delta rho, Delta V_{SCF} and symmetrizes them 26 ! f) If lda_plus_u=.true. compute also the response occupation 27 ! matrices dnsscf 28 ! g) (Introduced in February 2020) If noncolin=.true. and domag=.true. 29 ! the linear system is solved twice (nsolv = 2, the case 30 ! isolv = 2 needs the time-reversed wave functions). For the 31 ! theoretical background, please refer to Phys. Rev. B 100, 32 ! 045115 (2019) 33 ! 34 USE kinds, ONLY : DP 35 USE ions_base, ONLY : nat, ntyp => nsp, ityp 36 USE io_global, ONLY : stdout, ionode 37 USE io_files, ONLY : prefix, diropn 38 USE check_stop, ONLY : check_stop_now 39 USE wavefunctions, ONLY : evc 40 USE cell_base, ONLY : at 41 USE klist, ONLY : ltetra, lgauss, xk, wk, ngk, igk_k 42 USE gvect, ONLY : g 43 USE gvecs, ONLY : doublegrid 44 USE fft_base, ONLY : dfftp, dffts 45 USE lsda_mod, ONLY : lsda, nspin, current_spin, isk 46 USE spin_orb, ONLY : domag 47 USE wvfct, ONLY : nbnd, npwx, et 48 USE scf, ONLY : rho, vrs 49 USE uspp, ONLY : okvan, vkb, deeq_nc 50 USE uspp_param, ONLY : upf, nhm 51 USE noncollin_module, ONLY : noncolin, npol, nspin_mag 52 USE paw_variables, ONLY : okpaw 53 USE paw_onecenter, ONLY : paw_dpotential 54 USE paw_symmetry, ONLY : paw_dusymmetrize, paw_dumqsymmetrize 55 USE buffers, ONLY : save_buffer, get_buffer 56 USE control_ph, ONLY : rec_code, niter_ph, nmix_ph, tr2_ph, & 57 lgamma_gamma, convt, & 58 alpha_mix, rec_code_read, & 59 where_rec, flmixdpot, ext_recover 60 USE el_phon, ONLY : elph 61 USE uspp, ONLY : nlcc_any 62 USE units_ph, ONLY : iudrho, lrdrho, iudwf, lrdwf, iubar, lrbar, & 63 iudvscf, iuint3paw, lint3paw 64 USE units_lr, ONLY : iuwfc, lrwfc 65 USE output, ONLY : fildrho, fildvscf 66 USE phus, ONLY : becsumort, alphap, int1_nc 67 USE modes, ONLY : npertx, npert, u, t, tmq 68 USE recover_mod, ONLY : read_rec, write_rec 69 ! used to write fildrho: 70 USE dfile_autoname, ONLY : dfile_name 71 USE save_ph, ONLY : tmp_dir_save 72 ! used oly to write the restart file 73 USE mp_pools, ONLY : inter_pool_comm 74 USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp 75 USE mp, ONLY : mp_sum 76 USE efermi_shift, ONLY : ef_shift, ef_shift_paw, def 77 USE lrus, ONLY : int3_paw, becp1, int3_nc 78 USE lr_symm_base, ONLY : irotmq, minus_q, nsymq, rtau 79 USE eqv, ONLY : dvpsi, dpsi, evq 80 USE qpoint, ONLY : xq, nksq, ikks, ikqs 81 USE qpoint_aux, ONLY : ikmks, ikmkmqs, becpt, alphapt 82 USE control_lr, ONLY : nbnd_occ, lgamma 83 USE dv_of_drho_lr, ONLY : dv_of_drho 84 USE fft_helper_subroutines 85 USE fft_interfaces, ONLY : fft_interpolate 86 USE ldaU, ONLY : lda_plus_u 87 USE nc_mag_aux, ONLY : int1_nc_save, deeq_nc_save, int3_save 88 89 implicit none 90 91 integer :: irr, npe, imode0 92 ! input: the irreducible representation 93 ! input: the number of perturbation 94 ! input: the position of the modes 95 96 complex(DP) :: drhoscf (dfftp%nnr, nspin_mag, npe) 97 ! output: the change of the scf charge 98 99 real(DP) , allocatable :: h_diag (:,:) 100 ! h_diag: diagonal part of the Hamiltonian 101 real(DP) :: thresh, anorm, averlt, dr2, rsign 102 ! thresh: convergence threshold 103 ! anorm : the norm of the error 104 ! averlt: average number of iterations 105 ! dr2 : self-consistency error 106 ! rsign : sign or the term in the magnetization 107 real(DP) :: dos_ef, weight, aux_avg (2) 108 ! Misc variables for metals 109 ! dos_ef: density of states at Ef 110 111 complex(DP), allocatable, target :: dvscfin(:,:,:) 112 ! change of the scf potential 113 complex(DP), pointer :: dvscfins (:,:,:) 114 ! change of the scf potential (smooth part only) 115 complex(DP), allocatable :: drhoscfh (:,:,:), dvscfout (:,:,:) 116 ! change of rho / scf potential (output) 117 ! change of scf potential (output) 118 complex(DP), allocatable :: ldos (:,:), ldoss (:,:), mixin(:), mixout(:), & 119 dbecsum (:,:,:,:), dbecsum_nc(:,:,:,:,:,:), aux1 (:,:), tg_dv(:,:), & 120 tg_psic(:,:), aux2(:,:), drhoc(:), dbecsum_aux (:,:,:,:) 121 ! Misc work space 122 ! ldos : local density of states af Ef 123 ! ldoss: as above, without augmentation charges 124 ! dbecsum: the derivative of becsum 125 ! drhoc: response core charge density 126 REAL(DP), allocatable :: becsum1(:,:,:) 127 128 logical :: conv_root, & ! true if linear system is converged 129 exst, & ! used to open the recover file 130 lmetq0 ! true if xq=(0,0,0) in a metal 131 132 integer :: kter, & ! counter on iterations 133 iter0, & ! starting iteration 134 ipert, & ! counter on perturbations 135 ibnd, & ! counter on bands 136 iter, & ! counter on iterations 137 lter, & ! counter on iterations of linear system 138 ltaver, & ! average counter 139 lintercall, & ! average number of calls to cgsolve_all 140 ik, ikk, & ! counter on k points 141 ikq, & ! counter on k+q points 142 ig, & ! counter on G vectors 143 ndim, & 144 is, & ! counter on spin polarizations 145 nrec, & ! the record number for dvpsi and dpsi 146 ios, & ! integer variable for I/O control 147 incr, & ! used for tg 148 v_siz, & ! size of the potential 149 ipol, & ! counter on polarization 150 mode, & ! mode index 151 isolv, & ! counter on linear systems 152 nsolv, & ! number of linear systems 153 ikmk, & ! index of mk 154 ikmkmq ! index of mk-mq 155 156 integer :: npw, npwq 157 integer :: iq_dummy 158 real(DP) :: tcpu, get_clock ! timing variables 159 character(len=256) :: filename 160 161 external ch_psi_all, cg_psi 162 ! 163 IF (rec_code_read > 20 ) RETURN 164 165 call start_clock ('solve_linter') 166! 167! This routine is task group aware 168! 169 nsolv=1 170 IF (noncolin.AND.domag) nsolv=2 171 172 allocate (dvscfin ( dfftp%nnr , nspin_mag , npe)) 173 if (doublegrid) then 174 allocate (dvscfins (dffts%nnr , nspin_mag , npe)) 175 else 176 dvscfins => dvscfin 177 endif 178 allocate (drhoscfh ( dfftp%nnr, nspin_mag , npe)) 179 allocate (dvscfout ( dfftp%nnr, nspin_mag , npe)) 180 allocate (dbecsum ( (nhm * (nhm + 1))/2 , nat , nspin_mag , npe)) 181 IF (okpaw) THEN 182 allocate (mixin(dfftp%nnr*nspin_mag*npe+(nhm*(nhm+1)*nat*nspin_mag*npe)/2) ) 183 allocate (mixout(dfftp%nnr*nspin_mag*npe+(nhm*(nhm+1)*nat*nspin_mag*npe)/2) ) 184 mixin=(0.0_DP,0.0_DP) 185 ELSE 186 ALLOCATE(mixin(1)) 187 ENDIF 188 IF (noncolin) allocate (dbecsum_nc (nhm,nhm, nat , nspin , npe, nsolv)) 189 allocate (aux1 ( dffts%nnr, npol)) 190 allocate (h_diag ( npwx*npol, nbnd)) 191 allocate (aux2(npwx*npol, nbnd)) 192 allocate (drhoc(dfftp%nnr)) 193 IF (noncolin.AND.domag.AND.okvan) THEN 194 ALLOCATE (int3_save( nhm, nhm, nat, nspin_mag, npe, 2)) 195 ALLOCATE (dbecsum_aux ( (nhm * (nhm + 1))/2 , nat , nspin_mag , npe)) 196 ENDIF 197 incr=1 198 IF ( dffts%has_task_groups ) THEN 199 ! 200 v_siz = dffts%nnr_tg 201 ALLOCATE( tg_dv ( v_siz, nspin_mag ) ) 202 ALLOCATE( tg_psic( v_siz, npol ) ) 203 incr = fftx_ntgrp(dffts) 204 ! 205 ENDIF 206 ! 207 if (rec_code_read == 10.AND.ext_recover) then 208 ! restart from Phonon calculation 209 IF (okpaw) THEN 210 CALL read_rec(dr2, iter0, npe, dvscfin, dvscfins, drhoscfh, dbecsum) 211 IF (convt) THEN 212 CALL PAW_dpotential(dbecsum,rho%bec,int3_paw,npe) 213 ELSE 214 CALL setmixout(npe*dfftp%nnr*nspin_mag,& 215 (nhm*(nhm+1)*nat*nspin_mag*npe)/2,mixin,dvscfin,dbecsum,ndim,-1) 216 ENDIF 217 ELSE 218 CALL read_rec(dr2, iter0, npe, dvscfin, dvscfins, drhoscfh) 219 ENDIF 220 rec_code=0 221 else 222 iter0 = 0 223 convt =.FALSE. 224 where_rec='no_recover' 225 endif 226 227 IF (ionode .AND. fildrho /= ' ') THEN 228 INQUIRE (UNIT = iudrho, OPENED = exst) 229 IF (exst) CLOSE (UNIT = iudrho, STATUS='keep') 230 filename = dfile_name(xq, at, fildrho, TRIM(tmp_dir_save)//prefix, generate=.true., index_q=iq_dummy) 231 CALL diropn (iudrho, filename, lrdrho, exst) 232 END IF 233 234 IF (convt) GOTO 155 235 ! 236 ! if q=0 for a metal: allocate and compute local DOS at Ef 237 ! 238 239 lmetq0 = (lgauss .OR. ltetra) .AND. lgamma 240 if (lmetq0) then 241 allocate ( ldos ( dfftp%nnr , nspin_mag) ) 242 allocate ( ldoss( dffts%nnr , nspin_mag) ) 243 allocate (becsum1 ( (nhm * (nhm + 1))/2 , nat , nspin_mag)) 244 call localdos ( ldos , ldoss , becsum1, dos_ef ) 245 IF (.NOT.okpaw) deallocate(becsum1) 246 endif 247 ! 248 ! 249 ! In this case it has recovered after computing the contribution 250 ! to the dynamical matrix. This is a new iteration that has to 251 ! start from the beginning. 252 ! 253 IF (iter0==-1000) iter0=0 254 ! 255 ! The outside loop is over the iterations 256 ! 257 do kter = 1, niter_ph 258 ! 259 iter = kter + iter0 260 ltaver = 0 261 lintercall = 0 262 ! 263 drhoscf(:,:,:) = (0.d0, 0.d0) 264 dbecsum(:,:,:,:) = (0.d0, 0.d0) 265 IF (noncolin) dbecsum_nc = (0.d0, 0.d0) 266 ! 267 ! DFPT+U: at each ph iteration calculate dnsscf, 268 ! i.e. the scf variation of the occupation matrix ns. 269 ! 270 IF (lda_plus_u .AND. (iter.NE.1)) & 271 CALL dnsq_scf (npe, lmetq0, imode0, irr, .true.) 272 ! 273 do ik = 1, nksq 274 ! 275 ikk = ikks(ik) 276 ikq = ikqs(ik) 277 npw = ngk(ikk) 278 npwq= ngk(ikq) 279 ! 280 if (lsda) current_spin = isk (ikk) 281 ! 282 ! compute beta functions and kinetic energy for k-point ikq 283 ! needed by h_psi, called by ch_psi_all, called by cgsolve_all 284 ! 285 CALL init_us_2 (npwq, igk_k(1,ikq), xk (1, ikq), vkb) 286 CALL g2_kin (ikq) 287 ! 288 ! Start the loop on the two linear systems, one at B and one at 289 ! -B 290 ! 291 DO isolv=1, nsolv 292 IF (isolv==2) THEN 293 ikmk = ikmks(ik) 294 ikmkmq = ikmkmqs(ik) 295 rsign=-1.0_DP 296 ELSE 297 ikmk=ikk 298 ikmkmq=ikq 299 rsign=1.0_DP 300 ENDIF 301 ! 302 ! read unperturbed wavefunctions psi(k) and psi(k+q) 303 ! 304 if (nksq.gt.1.OR.nsolv==2) then 305 if (lgamma) then 306 call get_buffer (evc, lrwfc, iuwfc, ikmk) 307 else 308 call get_buffer (evc, lrwfc, iuwfc, ikmk) 309 call get_buffer (evq, lrwfc, iuwfc, ikmkmq) 310 end if 311 endif 312 ! 313 ! compute preconditioning matrix h_diag used by cgsolve_all 314 ! 315 CALL h_prec (ik, evq, h_diag) 316 ! 317 do ipert = 1, npe 318 mode = imode0 + ipert 319 nrec = (ipert - 1) * nksq + ik + (isolv-1) * npe * nksq 320 ! 321 ! and now adds the contribution of the self consistent term 322 ! 323 if (where_rec =='solve_lint'.or.iter>1) then 324 ! 325 ! After the first iteration dvbare_q*psi_kpoint is read from file 326 ! 327 call get_buffer (dvpsi, lrbar, iubar, nrec) 328 ! 329 ! calculates dvscf_q*psi_k in G_space, for all bands, k=kpoint 330 ! dvscf_q from previous iteration (mix_potential) 331 ! 332 call start_clock ('vpsifft') 333 ! 334 ! change the sign of the magnetic field if required 335 ! 336 IF (isolv==2) THEN 337 dvscfins(:,2:4,ipert)=-dvscfins(:,2:4,ipert) 338 IF (okvan) int3_nc(:,:,:,:,ipert)=int3_save(:,:,:,:,ipert,2) 339 ENDIF 340 ! 341 ! Set the potential for task groups 342 ! 343 IF( dffts%has_task_groups ) THEN 344 IF (noncolin) THEN 345 CALL tg_cgather( dffts, dvscfins(:,1,ipert), tg_dv(:,1)) 346 IF (domag) THEN 347 DO ipol=2,4 348 CALL tg_cgather( dffts, dvscfins(:,ipol,ipert), tg_dv(:,ipol)) 349 ENDDO 350 ENDIF 351 ELSE 352 CALL tg_cgather( dffts, dvscfins(:,current_spin,ipert), tg_dv(:,1)) 353 ENDIF 354 ENDIF 355 aux2=(0.0_DP,0.0_DP) 356 do ibnd = 1, nbnd_occ (ikk), incr 357 IF( dffts%has_task_groups ) THEN 358 call cft_wave_tg (ik, evc, tg_psic, 1, v_siz, ibnd, nbnd_occ (ikk) ) 359 call apply_dpot(v_siz, tg_psic, tg_dv, 1) 360 call cft_wave_tg (ik, aux2, tg_psic, -1, v_siz, ibnd, nbnd_occ (ikk)) 361 ELSE 362 call cft_wave (ik, evc (1, ibnd), aux1, +1) 363 call apply_dpot(dffts%nnr,aux1, dvscfins(1,1,ipert), current_spin) 364 call cft_wave (ik, aux2 (1, ibnd), aux1, -1) 365 ENDIF 366 enddo 367 dvpsi=dvpsi+aux2 368 call stop_clock ('vpsifft') 369 ! 370 ! In the case of US pseudopotentials there is an additional 371 ! selfconsist term which comes from the dependence of D on 372 ! V_{eff} on the bare change of the potential 373 ! 374 IF (isolv==1) THEN 375 call adddvscf_ph_mag (ipert, ik, becp1) 376 ! 377 ! DFPT+U: add to dvpsi the scf part of the response 378 ! Hubbard potential dV_hub 379 ! 380 if (lda_plus_u) call adddvhubscf (ipert, ik) 381 ELSE 382 call adddvscf_ph_mag (ipert, ik, becpt) 383 END IF 384 ! 385 ! reset the original magnetic field if it was changed 386 ! 387 IF (isolv==2) THEN 388 dvscfins(:,2:4,ipert)=-dvscfins(:,2:4,ipert) 389 IF (okvan) int3_nc(:,:,:,:,ipert)=int3_save(:,:,:,:,ipert,1) 390 ENDIF 391 ! 392 else 393 ! 394 ! At the first iteration dvbare_q*psi_kpoint is calculated 395 ! and written to file. 396 ! 397 IF (isolv==1) THEN 398 call dvqpsi_us (ik, u (1, mode),.false., becp1, alphap ) 399 ! 400 ! DFPT+U: At the first ph iteration the bare perturbed 401 ! Hubbard potential dvbare_hub_q * psi_kpoint 402 ! is calculated and added to dvpsi. 403 ! 404 if (lda_plus_u) call dvqhub_barepsi_us (ik, u(1,mode)) 405 ! 406 ELSE 407 IF (okvan) THEN 408 deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,2) 409 int1_nc(:,:,:,:,:)=int1_nc_save(:,:,:,:,:,2) 410 ENDIF 411 call dvqpsi_us (ik, u (1, mode),.false., becpt, alphapt) 412 IF (okvan) THEN 413 deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1) 414 int1_nc(:,:,:,:,:)=int1_nc_save(:,:,:,:,:,1) 415 ENDIF 416 ENDIF 417 call save_buffer (dvpsi, lrbar, iubar, nrec) 418 ! 419 endif 420 ! 421 ! Ortogonalize dvpsi to valence states: ps = <evq|dvpsi> 422 ! Apply -P_c^+. 423 ! 424 CALL orthogonalize(dvpsi, evq, ikmk, ikmkmq, dpsi, npwq, .false.) 425 ! 426 if (where_rec=='solve_lint'.or.iter > 1) then 427 ! 428 ! starting value for delta_psi is read from iudwf 429 ! 430 call get_buffer( dpsi, lrdwf, iudwf, nrec) 431 ! 432 ! threshold for iterative solution of the linear system 433 ! 434 thresh = min (1.d-1 * sqrt (dr2), 1.d-2) 435 else 436 ! 437 ! At the first iteration dpsi and dvscfin are set to zero 438 ! 439 dpsi(:,:) = (0.d0, 0.d0) 440 dvscfin (:, :, ipert) = (0.d0, 0.d0) 441 ! 442 ! starting threshold for iterative solution of the linear system 443 ! 444 thresh = 1.0d-2 445 endif 446 ! 447 ! iterative solution of the linear system (H-eS)*dpsi=dvpsi, 448 ! dvpsi=-P_c^+ (dvbare+dvscf)*psi , dvscf fixed. 449 ! 450 IF (isolv==2) THEN 451 vrs(:,2:4)=-vrs(:,2:4) 452 IF (okvan) deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,2) 453 ENDIF 454 conv_root = .true. 455 456 call cgsolve_all (ch_psi_all, cg_psi, et(1,ikmk), dvpsi, dpsi, & 457 h_diag, npwx, npwq, thresh, ik, lter, conv_root, & 458 anorm, nbnd_occ(ikk), npol ) 459 460 IF (isolv==2) THEN 461 vrs(:,2:4)=-vrs(:,2:4) 462 IF (okvan) deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1) 463 ENDIF 464 465 ltaver = ltaver + lter 466 lintercall = lintercall + 1 467 if (.not.conv_root) WRITE( stdout, '(5x,"kpoint",i4," ibnd",i4, & 468 & " solve_linter: root not converged ",es10.3)') & 469 & ik , ibnd, anorm 470 ! 471 ! writes delta_psi on iunit iudwf, k=kpoint, 472 ! 473 ! if (nksq.gt.1 .or. npert(irr).gt.1) 474 call save_buffer (dpsi, lrdwf, iudwf, nrec) 475 ! 476 ! calculates dvscf, sum over k => dvscf_q_ipert 477 ! 478 weight = wk (ikk) 479 IF (nsolv==2) weight=weight/2.0_DP 480 IF (noncolin) THEN 481 call incdrhoscf_nc(drhoscf(1,1,ipert),weight,ik, & 482 dbecsum_nc(1,1,1,1,ipert,isolv), dpsi, rsign) 483 ELSE 484 call incdrhoscf (drhoscf(1,current_spin,ipert), weight, ik, & 485 dbecsum(1,1,current_spin,ipert), dpsi) 486 END IF 487 ! on perturbations 488 enddo 489 ! on isolv 490 END DO 491 ! on k-points 492 enddo 493 ! 494 ! The calculation of dbecsum is distributed across processors (see addusdbec) 495 ! Sum over processors the contributions coming from each slice of bands 496 ! 497 IF (noncolin) THEN 498 call mp_sum ( dbecsum_nc, intra_bgrp_comm ) 499 ELSE 500 call mp_sum ( dbecsum, intra_bgrp_comm ) 501 ENDIF 502 503 if (doublegrid) then 504 do is = 1, nspin_mag 505 do ipert = 1, npe 506 call fft_interpolate (dffts, drhoscf(:,is,ipert), dfftp, drhoscfh(:,is,ipert)) 507 enddo 508 enddo 509 else 510 call zcopy (npe*nspin_mag*dfftp%nnr, drhoscf, 1, drhoscfh, 1) 511 endif 512 ! 513 ! In the noncolinear, spin-orbit case rotate dbecsum 514 ! 515 IF (noncolin.and.okvan) THEN 516 CALL set_dbecsum_nc(dbecsum_nc, dbecsum, npe) 517 IF (nsolv==2) THEN 518 dbecsum_aux=(0.0_DP,0.0_DP) 519 CALL set_dbecsum_nc(dbecsum_nc(1,1,1,1,1,2), dbecsum_aux, npe) 520 dbecsum(:,:,1,:)=dbecsum(:,:,1,:)+dbecsum_aux(:,:,1,:) 521 dbecsum(:,:,2:4,:)=dbecsum(:,:,2:4,:)-dbecsum_aux(:,:,2:4,:) 522 ENDIF 523 ENDIF 524 ! 525 ! Now we compute for all perturbations the total charge and potential 526 ! 527 call addusddens (drhoscfh, dbecsum, imode0, npe, 0) 528 ! 529 ! Reduce the delta rho across pools 530 ! 531 call mp_sum ( drhoscf, inter_pool_comm ) 532 call mp_sum ( drhoscfh, inter_pool_comm ) 533 IF (okpaw) call mp_sum ( dbecsum, inter_pool_comm ) 534 535 ! 536 ! q=0 in metallic case deserve special care (e_Fermi can shift) 537 ! 538 539 IF (okpaw) THEN 540 IF (lmetq0) & 541 call ef_shift_paw (drhoscfh, dbecsum, ldos, ldoss, becsum1, & 542 dos_ef, irr, npe, .false.) 543 DO ipert=1,npe 544 dbecsum(:,:,:,ipert)=2.0_DP *dbecsum(:,:,:,ipert) & 545 +becsumort(:,:,:,imode0+ipert) 546 ENDDO 547 ELSE 548 IF (lmetq0) call ef_shift(drhoscfh,ldos,ldoss,dos_ef,irr,npe,.false.) 549 ENDIF 550 ! 551 ! After the loop over the perturbations we have the linear change 552 ! in the charge density for each mode of this representation. 553 ! Here we symmetrize them ... 554 ! 555 IF (.not.lgamma_gamma) THEN 556 call psymdvscf (npe, irr, drhoscfh) 557 IF ( noncolin.and.domag ) CALL psym_dmag( npe, irr, drhoscfh) 558 IF (okpaw) THEN 559 IF (minus_q) CALL PAW_dumqsymmetrize(dbecsum,npe,irr, npertx,irotmq,rtau,xq,tmq) 560 CALL PAW_dusymmetrize(dbecsum,npe,irr,npertx,nsymq,rtau,xq,t) 561 END IF 562 ENDIF 563 ! 564 ! ... save them on disk and 565 ! compute the corresponding change in scf potential 566 ! 567 do ipert = 1, npe 568 if (fildrho.ne.' ') then 569 call davcio_drho (drhoscfh(1,1,ipert), lrdrho, iudrho, imode0+ipert, +1) 570! close(iudrho) 571 endif 572 ! 573 call zcopy (dfftp%nnr*nspin_mag,drhoscfh(1,1,ipert),1,dvscfout(1,1,ipert),1) 574 ! 575 ! Compute the response of the core charge density 576 ! IT: Should the condition "imode0+ipert > 0" be removed? 577 ! 578 if (imode0+ipert > 0) then 579 call addcore (imode0+ipert, drhoc) 580 else 581 drhoc(:) = (0.0_DP,0.0_DP) 582 endif 583 ! 584 ! Compute the response HXC potential 585 call dv_of_drho (dvscfout(1,1,ipert), .true., drhoc) 586 ! 587 enddo 588 ! 589 ! And we mix with the old potential 590 ! 591 IF (okpaw) THEN 592 ! 593 ! In this case we mix also dbecsum 594 ! 595 call setmixout(npe*dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag*npe)/2, & 596 mixout, dvscfout, dbecsum, ndim, -1 ) 597 call mix_potential (2*npe*dfftp%nnr*nspin_mag+2*ndim, mixout, mixin, & 598 alpha_mix(kter), dr2, npe*tr2_ph/npol, iter, & 599 nmix_ph, flmixdpot, convt) 600 call setmixout(npe*dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag*npe)/2, & 601 mixin, dvscfin, dbecsum, ndim, 1 ) 602 if (lmetq0.and.convt) & 603 call ef_shift_paw (drhoscf, dbecsum, ldos, ldoss, becsum1, & 604 dos_ef, irr, npe, .true.) 605 ELSE 606 call mix_potential (2*npe*dfftp%nnr*nspin_mag, dvscfout, dvscfin, & 607 alpha_mix(kter), dr2, npe*tr2_ph/npol, iter, & 608 nmix_ph, flmixdpot, convt) 609 if (lmetq0.and.convt) & 610 call ef_shift (drhoscf, ldos, ldoss, dos_ef, irr, npe, .true.) 611 ENDIF 612 ! check that convergent have been reached on ALL processors in this image 613 CALL check_all_convt(convt) 614 615 if (doublegrid) then 616 do ipert = 1, npe 617 do is = 1, nspin_mag 618 call fft_interpolate (dfftp, dvscfin(:,is,ipert), dffts, dvscfins(:,is,ipert)) 619 enddo 620 enddo 621 endif 622! 623! calculate here the change of the D1-~D1 coefficients due to the phonon 624! perturbation 625! 626 IF (okvan) THEN 627 IF (okpaw) CALL PAW_dpotential(dbecsum,rho%bec,int3_paw,npe) 628 ! 629 ! with the new change of the potential we compute the integrals 630 ! of the change of potential and Q 631 ! 632 call newdq (dvscfin, npe) 633 ! 634 ! In the noncollinear magnetic case computes the int3 coefficients with 635 ! the opposite sign of the magnetic field. They are saved in int3_save, 636 ! that must have been allocated by the calling routine 637 ! 638 IF (noncolin.AND.domag) THEN 639 int3_save(:,:,:,:,:,1)=int3_nc(:,:,:,:,:) 640 IF (okpaw) rho%bec(:,:,2:4)=-rho%bec(:,:,2:4) 641 DO ipert=1,npe 642 dvscfin(:,2:4,ipert)=-dvscfin(:,2:4,ipert) 643 IF (okpaw) dbecsum(:,:,2:4,ipert)=-dbecsum(:,:,2:4,ipert) 644 ENDDO 645 ! 646 ! if needed recompute the paw coeffients with the opposite sign of 647 ! the magnetic field 648 ! 649 IF (okpaw) CALL PAW_dpotential(dbecsum,rho%bec,int3_paw,npe) 650 ! 651 ! here compute the int3 integrals 652 ! 653 CALL newdq (dvscfin, npe) 654 int3_save(:,:,:,:,:,2)=int3_nc(:,:,:,:,:) 655 ! 656 ! restore the correct sign of the magnetic field. 657 ! 658 DO ipert=1,npe 659 dvscfin(:,2:4,ipert)=-dvscfin(:,2:4,ipert) 660 IF (okpaw) dbecsum(:,:,2:4,ipert)=-dbecsum(:,:,2:4,ipert) 661 ENDDO 662 IF (okpaw) rho%bec(:,:,2:4)=-rho%bec(:,:,2:4) 663 ! 664 ! put into int3_nc the coefficient with +B 665 ! 666 int3_nc(:,:,:,:,:)=int3_save(:,:,:,:,:,1) 667 ENDIF 668 END IF 669#if defined(__MPI) 670 aux_avg (1) = DBLE (ltaver) 671 aux_avg (2) = DBLE (lintercall) 672 call mp_sum ( aux_avg, inter_pool_comm ) 673 averlt = aux_avg (1) / aux_avg (2) 674#else 675 averlt = DBLE (ltaver) / lintercall 676#endif 677 tcpu = get_clock ('PHONON') 678 679 WRITE( stdout, '(/,5x," iter # ",i3," total cpu time :",f8.1, & 680 & " secs av.it.: ",f5.1)') iter, tcpu, averlt 681 dr2 = dr2 / npe 682 WRITE( stdout, '(5x," thresh=",es10.3, " alpha_mix = ",f6.3, & 683 & " |ddv_scf|^2 = ",es10.3 )') thresh, alpha_mix (kter) , dr2 684 ! 685 ! Here we save the information for recovering the run from this poin 686 ! 687 FLUSH( stdout ) 688 ! 689 rec_code=10 690 IF (okpaw) THEN 691 CALL write_rec('solve_lint', irr, dr2, iter, convt, npe, & 692 dvscfin, drhoscfh, dbecsum) 693 ELSE 694 CALL write_rec('solve_lint', irr, dr2, iter, convt, npe, & 695 dvscfin, drhoscfh) 696 ENDIF 697 698 if (check_stop_now()) call stop_smoothly_ph (.false.) 699 if (convt) goto 155 700 enddo 701155 iter0=0 702 ! 703 ! A part of the dynamical matrix requires the integral of 704 ! the self consistent change of the potential and the variation of 705 ! the charge due to the displacement of the atoms. 706 ! We compute it here. 707 ! 708 if (convt) then 709 call drhodvus (irr, imode0, dvscfin, npe) 710 if (fildvscf.ne.' ') then 711 do ipert = 1, npe 712 if(lmetq0) then 713 dvscfin(:,:,ipert) = dvscfin(:,:,ipert)-def(ipert) 714 if (doublegrid) dvscfins(:,:,ipert) = dvscfins(:,:,ipert)-def(ipert) 715 endif 716 call davcio_drho ( dvscfin(1,1,ipert), lrdrho, iudvscf, imode0 + ipert, +1 ) 717 IF (okpaw.AND.me_bgrp==0) CALL davcio( int3_paw(:,:,:,:,ipert), lint3paw, & 718 iuint3paw, imode0+ipert, + 1 ) 719 end do 720 if (elph) call elphel (irr, npe, imode0, dvscfins) 721 end if 722 endif 723 if (convt.and.nlcc_any) call addnlcc (imode0, drhoscfh, npe) 724 if (allocated(ldoss)) deallocate (ldoss) 725 if (allocated(ldos)) deallocate (ldos) 726 deallocate (h_diag) 727 deallocate (aux1) 728 deallocate (dbecsum) 729 IF (okpaw) THEN 730 if (allocated(becsum1)) deallocate (becsum1) 731 deallocate (mixin) 732 deallocate (mixout) 733 ENDIF 734 IF (noncolin) deallocate (dbecsum_nc) 735 deallocate (dvscfout) 736 deallocate (drhoscfh) 737 if (doublegrid) deallocate (dvscfins) 738 deallocate (dvscfin) 739 deallocate(aux2) 740 deallocate(drhoc) 741 IF ( dffts%has_task_groups ) THEN 742 DEALLOCATE( tg_dv ) 743 DEALLOCATE( tg_psic ) 744 ENDIF 745 IF (noncolin.AND.domag.AND.okvan) THEN 746 DEALLOCATE (int3_save) 747 DEALLOCATE (dbecsum_aux) 748 ENDIF 749 750 call stop_clock ('solve_linter') 751 752 RETURN 753 754END SUBROUTINE solve_linter 755 756SUBROUTINE check_all_convt(convt) 757 USE mp, ONLY : mp_sum 758 USE mp_images, ONLY : nproc_image, me_image, intra_image_comm 759 IMPLICIT NONE 760 LOGICAL,INTENT(in) :: convt 761 INTEGER :: tot_conv 762 ! 763 IF(nproc_image==1) RETURN 764 ! 765 ! Work out how many processes have converged 766 ! 767 tot_conv = 0 768 IF(convt) tot_conv = 1 769 CALL mp_sum(tot_conv, intra_image_comm) 770 ! 771 IF ((tot_conv > 0) .and. (tot_conv < nproc_image)) THEN 772 CALL errore('check_all_convt', 'Only some processors converged: '& 773 &' either something is wrong with solve_linter, or a different'& 774 &' parallelism scheme should be used.', 1) 775 ENDIF 776 ! 777 RETURN 778 ! 779END SUBROUTINE 780