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 hp_solve_linear_system (na, iq) 11 !----------------------------------------------------------------------- 12 ! 13 ! This is a driver routine for the solution of the linear-response Kohn-Sham 14 ! equations (45) in Ref. [1]. The solution defines the change of Kohn-Sham 15 ! wavefunctions due change of occupations. 16 ! [1] Phys. Rev. B 98, 085127 (2018) 17 ! 18 USE kinds, ONLY : DP 19 USE ions_base, ONLY : nat 20 USE io_global, ONLY : stdout 21 USE check_stop, ONLY : check_stop_now 22 USE wavefunctions, ONLY : evc 23 USE cell_base, ONLY : tpiba2 24 USE klist, ONLY : lgauss, ltetra, xk, wk, nelec, ngk, igk_k 25 USE gvect, ONLY : g 26 USE gvecs, ONLY : doublegrid 27 USE scf, ONLY : rho 28 USE fft_base, ONLY : dfftp, dffts 29 USE lsda_mod, ONLY : lsda, current_spin, isk 30 USE wvfct, ONLY : nbnd, npwx, g2kin, et 31 USE uspp, ONLY : okvan, vkb, nkb 32 USE uspp_param, ONLY : nhm 33 USE becmod, ONLY : allocate_bec_type, deallocate_bec_type, becp 34 USE buffers, ONLY : save_buffer, get_buffer 35 USE noncollin_module, ONLY : npol, nspin_mag 36 USE paw_variables, ONLY : okpaw 37 USE paw_onecenter, ONLY : paw_dpotential 38 USE paw_symmetry, ONLY : paw_dusymmetrize, paw_dumqsymmetrize 39 USE mp_pools, ONLY : inter_pool_comm, intra_pool_comm 40 USE mp, ONLY : mp_sum 41 USE hp_efermi_shift, ONLY : hp_ef_shift, def 42 USE eqv, ONLY : dvpsi, dpsi, evq 43 USE qpoint, ONLY : nksq, ikks, ikqs, xq 44 USE control_lr, ONLY : lgamma, nbnd_occ 45 USE units_lr, ONLY : iuwfc, lrwfc 46 USE lrus, ONLY : int3, int3_paw 47 USE dv_of_drho_lr, ONLY : dv_of_drho 48 USE fft_helper_subroutines 49 USE fft_interfaces, ONLY : fft_interpolate 50 USE lr_symm_base, ONLY : irotmq, minus_q, nsymq, rtau 51 USE ldaU_hp, ONLY : thresh_init, dnsscf, dns0, trace_dns_tot_old, & 52 conv_thr_chi_best, iter_best, niter_max, nmix, & 53 alpha_mix, iudwfc, lrdwfc, code 54 ! 55 IMPLICIT NONE 56 ! 57 INTEGER, INTENT(IN) :: na ! number of the perturbed atom 58 INTEGER, INTENT(IN) :: iq ! number of the q point 59 ! 60 REAL(DP), ALLOCATABLE :: h_diag (:,:) ! diagonal part of the Hamiltonian 61 ! 62 REAL(DP) :: thresh, & ! convergence threshold 63 anorm, & ! the norm of the error 64 averlt, & ! average number of iterations 65 dr2 ! self-consistency error 66 ! 67 REAL(DP) :: dos_ef, & ! density of states at the Fermi level 68 weight, & ! Misc variables for metals 69 aux_avg(2) ! Misc variables for metals 70 ! 71 REAL(DP), ALLOCATABLE :: becsum1(:,:,:) 72 ! 73 COMPLEX(DP), ALLOCATABLE, TARGET :: dvscfin(:,:) 74 ! change of the scf potential (input) 75 ! 76 COMPLEX(DP), POINTER :: dvscfins(:,:) 77 ! change of the scf potential (smooth part only) 78 ! 79 COMPLEX(DP), ALLOCATABLE :: drhoscf (:,:), & 80 drhoscfh (:,:), & 81 dvscfout (:,:) 82 ! change of rho / scf potential (output) 83 ! 84 COMPLEX(DP), ALLOCATABLE :: & 85 ldos (:,:), & ! local density of states at Ef 86 ldoss (:,:), & ! as above, without augmentation charges 87 dbecsum (:,:,:,:), & ! the derivative of becsum 88 aux1 (:,:), aux2 (:,:), & ! auxiliary arrays 89 mixin(:), mixout(:), & ! auxiliary arrays for mixing of the response potential 90 tg_dv(:,:), & ! Task groups: auxiliary array for potential * wfct 91 tg_psic(:,:) ! Task groups: auxiliary array for wavefunctions 92 93 COMPLEX(DP), ALLOCATABLE :: t(:,:,:,:), tmq(:,:,:) 94 ! PAW: auxiliary arrays 95 96 LOGICAL :: conv_root, & ! true if linear system is converged 97 exst, & ! used to open the recover file 98 lmetq0, & ! true if xq=(0,0,0) in a metal 99 convt, & ! not needed for HP 100 convt_chi ! used instead of convt to control the convergence 101 102 REAL(DP), PARAMETER :: tr2 = 1.D-30 ! threshold parameter 103 104 INTEGER :: ibnd, & ! counter on bands 105 iter, & ! counter on iterations 106 lter, & ! counter on iterations of linear system 107 ltaver, & ! average counter 108 lintercall, & ! average number of calls to cgsolve_all 109 ik, ikk, & ! counter on k points 110 ikq, & ! counter on k+q points 111 ig, & ! counter on G vectors 112 ndim, & 113 is, & ! counter on spin polarizations 114 nt, & ! counter on types 115 ios, & ! integer variable for I/O control 116 incr, & ! used for task groups 117 v_siz, & ! size of the potential 118 npw, & ! number of plane waves at k 119 npwq ! number of plane waves at k+q 120 121 REAL(DP) :: tcpu, get_clock ! timing variables 122 CHARACTER(LEN=256) :: filename, & 123 flmixdpot = 'mixd' 124 EXTERNAL ch_psi_all, cg_psi 125 ! 126 CALL start_clock ('hp_solve_linear_system') 127 ! 128 WRITE( stdout,*) " =--------------------------------------------=" 129 WRITE( stdout, '(13x,"START SOLVING THE LINEAR SYSTEM")') 130 WRITE( stdout,*) " =--------------------------------------------=" 131 ! 132 ! Allocate arrays for the SCF density/potential 133 ! 134 ALLOCATE (drhoscf (dffts%nnr, nspin_mag)) 135 ALLOCATE (drhoscfh(dfftp%nnr, nspin_mag)) 136 ALLOCATE (dvscfin (dfftp%nnr, nspin_mag)) 137 ALLOCATE (dvscfout(dfftp%nnr, nspin_mag)) 138 ! 139 IF (doublegrid) THEN 140 ALLOCATE (dvscfins(dffts%nnr, nspin_mag)) 141 ELSE 142 dvscfins => dvscfin 143 ENDIF 144 ! 145 ! USPP-specific allocations 146 ! 147 IF (okvan) ALLOCATE (int3 ( nhm, nhm, nat, nspin_mag, 1)) 148 IF (okpaw) ALLOCATE (int3_paw ( nhm, nhm, nat, nspin_mag, 1)) 149 CALL allocate_bec_type (nkb, nbnd, becp) 150 ! 151 ALLOCATE (dbecsum((nhm*(nhm+1))/2, nat, nspin_mag, 1)) 152 ! 153 IF (okpaw) THEN 154 ! 155 ALLOCATE (mixin(dfftp%nnr*nspin_mag+(nhm*(nhm+1)*nat*nspin_mag)/2) ) 156 ALLOCATE (mixout(dfftp%nnr*nspin_mag+(nhm*(nhm+1)*nat*nspin_mag)/2) ) 157 mixin = (0.0_DP,0.0_DP) 158 ! 159 ! Auxiliary unitary arrays 160 ALLOCATE ( tmq(1,1,3*nat) ) 161 ALLOCATE ( t(1,1,48,3*nat) ) 162 t(:,:,:,:) = (1.0_DP, 0.0_DP) 163 tmq(:,:,:) = (1.0_DP, 0.0_DP) 164 ! 165 ENDIF 166 ! 167 ALLOCATE (aux1(dffts%nnr, npol)) 168 ALLOCATE (aux2(npwx*npol, nbnd)) 169 ALLOCATE (h_diag(npwx*npol, nbnd)) 170 ALLOCATE (trace_dns_tot_old(nat)) 171 trace_dns_tot_old(:) = (0.d0, 0.d0) 172 ! 173 convt = .FALSE. 174 convt_chi = .FALSE. 175 ! 176 incr = 1 177 IF ( dffts%has_task_groups ) THEN 178 ! 179 v_siz = dffts%nnr_tg 180 ALLOCATE( tg_dv ( v_siz, nspin_mag ) ) 181 ALLOCATE( tg_psic( v_siz, npol ) ) 182 incr = fftx_ntgrp(dffts) 183 ! 184 ENDIF 185 ! 186 ! If q=0 for a metal: allocate and compute local DOS and DOS at Ef 187 ! 188 lmetq0 = (lgauss .OR. ltetra) .AND. lgamma 189 ! 190 IF (lmetq0) THEN 191 ALLOCATE (ldos (dfftp%nnr, nspin_mag)) 192 ALLOCATE (ldoss(dffts%nnr, nspin_mag)) 193 ALLOCATE (becsum1 ( (nhm * (nhm + 1))/2, nat, nspin_mag)) 194 CALL localdos (ldos, ldoss, becsum1, dos_ef) 195 IF (.NOT.okpaw) DEALLOCATE (becsum1) 196 ENDIF 197 ! 198 ! The loop of the linear-response calculation 199 ! 200 DO iter = 1, niter_max 201 ! 202 WRITE(stdout,'(/6x,"atom #",i3,3x,"q point #",i4,3x,"iter # ",i3)') na, iq, iter 203 ! 204 ltaver = 0 205 lintercall = 0 206 ! 207 drhoscf(:,:) = (0.d0, 0.d0) 208 dvscfout(:,:) = (0.d0, 0.d0) 209 dbecsum(:,:,:,:) = (0.d0, 0.d0) 210 ! 211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 212 !!!!!!!!!!!!!!!!!!! START OF THE K LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 213 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 214 ! 215 DO ik = 1, nksq 216 ! 217 ikk = ikks(ik) 218 ikq = ikqs(ik) 219 npw = ngk(ikk) 220 npwq = ngk(ikq) 221 ! 222 IF (lsda) current_spin = isk(ikk) 223 ! 224 ! Read unperturbed KS wavefuctions psi(k) and psi(k+q) 225 ! 226 IF (nksq.gt.1) THEN 227 IF (lgamma) THEN 228 CALL get_buffer (evc, lrwfc, iuwfc, ikk) 229 ELSE 230 CALL get_buffer (evc, lrwfc, iuwfc, ikk) 231 CALL get_buffer (evq, lrwfc, iuwfc, ikq) 232 ENDIF 233 ENDIF 234 ! 235 ! USPP: Compute the projectors vkb at k+q 236 ! 237 CALL init_us_2 (npwq, igk_k(1,ikq), xk(1,ikq), vkb) 238 ! 239 ! Compute the kinetic energy at k+q 240 ! 241 CALL g2_kin (ikq) 242 ! 243 ! Compute preconditioning matrix h_diag used by cgsolve_all 244 ! 245 CALL h_prec (ik, evq, h_diag) 246 ! 247 ! Computes (iter=1) or reads (iter>1) the action of the perturbing 248 ! potential on the unperturbed KS wavefunctions: |dvpsi> = dV_pert * |evc> 249 ! See Eq. (46) in Ref. [1] 250 ! 251 CALL hp_dvpsi_pert(ik) 252 ! 253 IF ( iter > 1 ) THEN 254 ! 255 ! Add the contribution of the self consistent term. 256 ! Calculates dvscf_q*psi(k) in G-space, for all bands, k=ik 257 ! dvscf_q from previous iteration (mix_potential) 258 ! 259 CALL start_clock ('hp_vpsifft') 260 IF ( dffts%has_task_groups ) & 261 & CALL tg_cgather( dffts, dvscfins(:,current_spin), tg_dv(:,1)) 262 aux2 = (0.0_DP, 0.0_DP) 263 DO ibnd = 1, nbnd_occ(ikk), incr 264 IF ( dffts%has_task_groups ) THEN 265 CALL cft_wave_tg (ik, evc, tg_psic, 1, v_siz, ibnd, nbnd_occ(ikk) ) 266 CALL apply_dpot(v_siz, tg_psic, tg_dv, 1) 267 CALL cft_wave_tg (ik, aux2, tg_psic, -1, v_siz, ibnd, nbnd_occ(ikk)) 268 ELSE 269 CALL cft_wave (ik, evc(:,ibnd), aux1, +1) 270 CALL apply_dpot(dffts%nnr, aux1, dvscfins, current_spin) 271 CALL cft_wave (ik, aux2(:,ibnd), aux1, -1) 272 ENDIF 273 ENDDO 274 dvpsi = dvpsi + aux2 275 CALL stop_clock ('hp_vpsifft') 276 ! 277 ! USPP: there is an additional self-consistent term proportional to int3 278 ! |dvpsi> = |dvpsi> + dV_HXC*|evc> + int3 * |beta><beta|evc> 279 ! 280 IF (okvan) CALL adddvscf(1, ik) 281 ! 282 ENDIF 283 ! 284 ! Ortogonalize dvpsi to valence states: ps = <evq|dvpsi> 285 ! Apply -P_c^+. See Eq. (A21) in Ref. [1] 286 ! 287 CALL orthogonalize(dvpsi, evq, ikk, ikq, dpsi, npwq, .FALSE.) 288 ! 289 IF ( iter == 1 ) THEN 290 ! 291 ! At the first iteration dpsi and dvscfin are set to zero 292 ! 293 dpsi(:,:) = (0.d0, 0.d0) 294 dvscfin(:,:) = (0.d0, 0.d0) 295 ! 296 ! Starting threshold for iterative solution of the linear system. 297 ! A strickt threshold for the first iteration is needed, 298 ! because we need dns0 to very high precision. 299 ! 300 thresh = thresh_init * nelec 301 ! 302 ELSE 303 ! 304 ! Starting value for dpsi is read from file 305 ! 306 CALL get_buffer( dpsi, lrdwfc, iudwfc, ik) 307 ! 308 ! Threshold for iterative solution of the linear system. 309 ! We start with not a strict threshold for iter=2, and then 310 ! it decreases with iterations. 311 ! 312 thresh = MIN (1.D-1 * SQRT(dr2), 1.D-2) 313 ! 314 ENDIF 315 ! 316 ! Iterative solution of the linear system: 317 ! (H + Q - eS) * |dpsi> = |dvpsi>, 318 ! where |dvpsi> = - P_c^+ (dV_HXC + dV_pert) * |evc> 319 ! See Eq. (43) in Ref. [1] 320 ! 321 CALL cgsolve_all (ch_psi_all, cg_psi, et(1,ikk), dvpsi, dpsi, h_diag, & 322 & npwx, npwq, thresh, ik, lter, conv_root, anorm, nbnd_occ(ikk), npol ) 323 ! 324 ltaver = ltaver + lter 325 ! 326 lintercall = lintercall + 1 327 ! 328 IF (.NOT.conv_root) THEN 329 WRITE( stdout, '(6x,"kpoint",i4, & 330 & " hp_solve_linear_system: root not converged, thresh < ",e10.3)') ik , anorm 331 IF (iter == 1) WRITE( stdout, '(6x,"Try to increase thresh_init...")') 332 ENDIF 333 ! 334 ! Writes dpsi on file for a given k 335 ! 336 CALL save_buffer (dpsi, lrdwfc, iudwfc, ik) 337 ! 338 ! Setup the weight at point k (normalized by the number of k points) 339 ! 340 weight = wk(ikk) 341 ! 342 ! Calculates the response charge density (sum over k) 343 ! See Eq. (48) in Ref. [1] 344 ! 345 CALL incdrhoscf (drhoscf(:,current_spin), weight, ik, & 346 & dbecsum(:,:,current_spin,1), dpsi) 347 ! 348 ENDDO ! k points 349 ! 350 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 351 !!!!!!!!!!!!!!!!!!!!! END OF THE K LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 352 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 353 ! 354#if defined (__MPI) 355 ! 356 ! USPP: The calculation of dbecsum is distributed across processors (see addusdbec) 357 ! Sum over processors the contributions coming from each slice of bands 358 ! 359 CALL mp_sum ( dbecsum, intra_pool_comm ) 360 ! 361#endif 362 ! 363 ! Copy/interpolate the response density drhoscf -> drhoscfh 364 ! 365 IF (doublegrid) THEN 366 do is = 1, nspin_mag 367 CALL fft_interpolate (dffts, drhoscf(:,is), dfftp, drhoscfh(:,is)) 368 enddo 369 ELSE 370 CALL zcopy (nspin_mag*dfftp%nnr, drhoscf, 1, drhoscfh, 1) 371 ENDIF 372 ! 373 ! USPP: Compute the total response charge density (standard term + US term) 374 ! 375 IF (okvan) CALL lr_addusddens (drhoscfh, dbecsum) 376 ! 377#if defined (__MPI) 378 CALL mp_sum ( drhoscfh, inter_pool_comm ) 379 IF (okpaw) CALL mp_sum ( dbecsum, inter_pool_comm ) 380#endif 381 ! 382 ! PAW: the factor of 2 is due to the presence of the CC term 383 ! (see first two terms in Eq.(9) in PRB 81, 075123 (2010)) 384 ! 385 IF (okpaw) dbecsum = 2.0_DP * dbecsum 386 ! 387 ! Metallic case and q=0: add a correction to the response charge density 388 ! due to the shift of the Fermi energy (see Eq.(75) in Rev. Mod. Phys. 73, 515 (2001)). 389 ! This term is added to the response charge density (in order to obtain correct 390 ! response HXC potential) and to the response occupation matrices. 391 ! 392 IF (lmetq0) THEN 393 ! 394 IF (okpaw) THEN 395 CALL hp_ef_shift (drhoscfh, ldos, ldoss, dos_ef, dbecsum, becsum1) 396 ELSE 397 CALL hp_ef_shift (drhoscfh, ldos, ldoss, dos_ef) 398 ENDIF 399 ! 400 ! Check that def is not too large (it is in Ry). 401 ! 402 IF ( ABS(DBLE(def)) > 5.0d0 ) THEN 403 ! 404 WRITE( stdout, '(/6x,"WARNING: The Fermi energy shift is too big!")') 405 WRITE( stdout, '(6x, "This may happen in two cases:")') 406 WRITE( stdout, '(6x, "1. The DOS at the Fermi level is too small:")') 407 WRITE( stdout, '(6x, " DOS(E_Fermi) = ",1x,2e12.4)') dos_ef 408 WRITE( stdout, '(6x, " This means that most likely the system has a gap,")') 409 WRITE( stdout, '(6x, " and hence it should NOT be treated as a metal")') 410 WRITE( stdout, '(6x, " (otherwise numerical instabilities will appear).")') 411 WRITE( stdout, '(6x, "2. Numerical instabilities due to too low cutoff")') 412 WRITE( stdout, '(6x, " for hard pseudopotentials.")') 413 WRITE( stdout, '(/6x,"Stopping...")') 414 ! 415 CALL hp_stop_smoothly (.FALSE.) 416 ! 417 ENDIF 418 ! 419 ENDIF 420 ! 421 ! Symmetrization of the response charge density. 422 ! 423 CALL hp_psymdvscf (drhoscfh) 424 ! 425 ! Symmetrize dbecsum 426 ! 427 IF (okpaw) THEN 428 IF (minus_q) CALL PAW_dumqsymmetrize(dbecsum,1,1,1,irotmq,rtau,xq,tmq) 429 CALL PAW_dusymmetrize(dbecsum,1,1,1,nsymq,rtau,xq,t) 430 ENDIF 431 ! 432 ! Copy drhoscfh to dvscfout 433 ! 434 CALL zcopy (nspin_mag*dfftp%nnr, drhoscfh, 1, dvscfout, 1) 435 ! 436 ! Compute the response potential dV_HXC from the response charge density. 437 ! See Eq. (47) in Ref. [1] 438 ! 439 CALL dv_of_drho (dvscfout, .FALSE.) 440 ! 441 ! Mix the new HXC response potential (dvscfout) with the old one (dvscfin). 442 ! Note: dvscfin = 0 for iter = 1 (so there is no mixing). 443 ! Output: dvscfin becomes a mixed potential 444 ! dvscfout contains the difference vout-vin (it is not needed) 445 ! PAW: mix also dbecsum 446 ! 447 IF (okpaw) THEN 448 CALL setmixout(dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag)/2, & 449 & mixout, dvscfout, dbecsum, ndim, -1 ) 450 CALL mix_potential (2*dfftp%nnr*nspin_mag+2*ndim, mixout, mixin, & 451 & alpha_mix(iter), dr2, tr2/npol, iter, nmix, flmixdpot, convt) 452 CALL setmixout(dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag)/2, & 453 & mixin, dvscfin, dbecsum, ndim, 1 ) 454 ELSE 455 CALL mix_potential (2*dfftp%nnr*nspin_mag, dvscfout, dvscfin, & 456 & alpha_mix(iter), dr2, tr2/npol, iter, nmix, flmixdpot, convt) 457 ENDIF 458 ! 459 ! NCPP case: dvscfins is a pointer to dvscfin. 460 ! USPP case: Interpolate dvscfin from dense to smooth grid 461 ! and put the result in dvscfins (needed for the next iteration). 462 ! 463 IF (doublegrid) THEN 464 DO is = 1, nspin_mag 465 CALL fft_interpolate (dfftp, dvscfin(:,is), dffts, dvscfins(:,is)) 466 ENDDO 467 ENDIF 468 ! 469 ! PAW: compute the response PAW potential 470 ! (see the last term in Eq.(12) in PRB 81, 075123 (2010)) 471 ! 472 IF (okpaw) CALL PAW_dpotential(dbecsum,rho%bec,int3_paw,1) 473 ! 474 ! USPP: With the new change of the HXC response potential dV_HXC 475 ! we compute the integral of the change of this potential and Q. 476 ! int3 = \int Q(r) dV_HXC(r) dr 477 ! PAW: int3_paw is added to int3 inside of the routine newdq 478 ! 479 IF (okvan) CALL newdq (dvscfin, 1) 480 ! 481 ! Calculate the response occupation matrix 482 ! See Eq. (43) in Ref. [1] 483 ! 484 CALL hp_dnsq (lmetq0, iter, convt_chi, dnsscf(:,:,:,:,iq)) 485 ! 486 ! Save the response occupation matrix after the first iteration, 487 ! which was computed from dpsi corresponding to the perturbing 488 ! potential only (dV_HXC=0). This is needed for the calculation of chi0. 489 ! 490 IF ( iter == 1 ) dns0(:,:,:,:,iq) = dnsscf(:,:,:,:,iq) 491 ! 492 ! Compute the average number of iterations 493 ! 494#if defined (__MPI) 495 aux_avg(1) = DBLE(ltaver) 496 aux_avg(2) = DBLE(lintercall) 497 CALL mp_sum ( aux_avg, inter_pool_comm ) 498 averlt = aux_avg(1) / aux_avg(2) 499#else 500 averlt = DBLE(ltaver) / DBLE(lintercall) 501#endif 502 ! 503 tcpu = get_clock(code) 504 ! 505 WRITE( stdout, '(6x,"Average number of iter. to solve lin. system:",2x,f5.1)') averlt 506 WRITE( stdout, '(6x,"Total CPU time :",f8.1,1x,"s")') tcpu 507 ! 508 IF ( check_stop_now() ) CALL hp_stop_smoothly (.FALSE.) 509 ! 510 IF ( iter==niter_max .AND. .NOT.convt_chi) THEN 511 WRITE( stdout, '(/6x,"Convergence has not been reached after",1x,i3,1x,"iterations!")') niter_max 512 IF ( iter > 1 ) THEN 513 WRITE( stdout, '(6x,"The best overall accuracy which was reached :")') 514 WRITE( stdout, '(6x,"diff = ",1x,f14.10,1x," iter =",1x,i3)') conv_thr_chi_best, iter_best 515 ENDIF 516 WRITE( stdout, '(/6x,"Stopping...")') 517 CALL hp_stop_smoothly (.TRUE.) 518 ENDIF 519 ! 520 IF (convt_chi) goto 155 521 ! 522 ENDDO ! loop over the iterations iter 523 ! 524155 CONTINUE 525 ! 526 DEALLOCATE (h_diag) 527 DEALLOCATE (aux1) 528 DEALLOCATE (aux2) 529 DEALLOCATE (dbecsum) 530 DEALLOCATE (drhoscf) 531 DEALLOCATE (drhoscfh) 532 DEALLOCATE (dvscfin) 533 DEALLOCATE (dvscfout) 534 DEALLOCATE (trace_dns_tot_old) 535 IF (doublegrid) DEALLOCATE (dvscfins) 536 IF (ALLOCATED(ldoss)) DEALLOCATE (ldoss) 537 IF (ALLOCATED(ldos)) DEALLOCATE (ldos) 538 IF (ALLOCATED(becsum1)) DEALLOCATE (becsum1) 539 IF (okvan) DEALLOCATE (int3) 540 IF (okpaw) THEN 541 DEALLOCATE (int3_paw) 542 DEALLOCATE (mixin, mixout) 543 DEALLOCATE (t) 544 DEALLOCATE (tmq) 545 ENDIF 546 CALL deallocate_bec_type (becp) 547 IF ( dffts%has_task_groups ) THEN 548 DEALLOCATE( tg_dv ) 549 DEALLOCATE( tg_psic ) 550 ENDIF 551 ! 552 WRITE( stdout,*) " " 553 WRITE( stdout,*) " =--------------------------------------------=" 554 WRITE( stdout, '(13x,"CONVERGENCE HAS BEEN REACHED")') 555 WRITE( stdout,*) " =--------------------------------------------=" 556 ! 557 CALL stop_clock ('hp_solve_linear_system') 558 ! 559 RETURN 560 ! 561END SUBROUTINE hp_solve_linear_system 562