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! 8! 9#define ZERO ( 0._dp, 0._dp ) 10! 11! This macro force the normalization of betamix matrix, usually not necessary 12!#define __NORMALIZE_BETAMIX 13! 14#if defined(__GFORTRAN__) 15#if (__GNUC__<4) || ((__GNUC__==4) && (__GNUC_MINOR__<8)) 16#define __GFORTRAN_HACK 17#endif 18#endif 19 20MODULE mix_save 21#if defined(__GFORTRAN_HACK) 22! gfortran hack - for some mysterious reason gfortran doesn't save 23! derived-type variables even with the SAVE attribute 24 USE scf, ONLY : mix_type 25 TYPE(mix_type), ALLOCATABLE, SAVE :: & 26 df(:), &! information from preceding iterations 27 dv(:) ! " " " " " " 28#endif 29END MODULE mix_save 30 31!---------------------------------------------------------------------------- 32SUBROUTINE mix_rho( input_rhout, rhoin, alphamix, dr2, tr2_min, iter, n_iter,& 33 iunmix, conv ) 34 !---------------------------------------------------------------------------- 35 !! * Modified Broyden's method for charge density mixing: D.D. Johnson, 36 !! PRB 38, 12807 (1988) ; 37 !! * Thomas-Fermi preconditioning described in: Raczkowski, Canning, Wang, 38 !! PRB 64,121101 (2001) ; 39 !! * Extended to mix also quantities needed for PAW, meta-GGA, DFT+U(+V) ; 40 !! * Electric field (all these are included into \(\text{mix_type}\)) ; 41 !! * On output: the mixed density is in \(\text{rhoin}\), 42 !! \(\text{input_rhout}\) is unchanged. 43 ! 44 USE kinds, ONLY : DP 45 USE ions_base, ONLY : nat, ntyp => nsp 46 USE gvect, ONLY : ngm 47 USE gvecs, ONLY : ngms 48 USE lsda_mod, ONLY : nspin 49 USE control_flags, ONLY : imix, ngm0, tr2, io_level 50 ! ... for PAW: 51 USE uspp_param, ONLY : nhm 52 USE scf, ONLY : scf_type, create_scf_type, destroy_scf_type, & 53 mix_type, create_mix_type, destroy_mix_type, & 54 assign_scf_to_mix_type, assign_mix_to_scf_type, & 55 mix_type_AXPY, davcio_mix_type, rho_ddot, & 56 high_frequency_mixing, nsg_ddot, & 57 mix_type_COPY, mix_type_SCAL 58 USE io_global, ONLY : stdout 59 USE ldaU, ONLY : lda_plus_u, lda_plus_u_kind, ldim_u, & 60 max_num_neighbors, nsg, nsgnew 61 USE io_files, ONLY : diropn 62#if defined(__GFORTRAN_HACK) 63 USE mix_save 64#endif 65 ! 66 IMPLICIT NONE 67 ! 68 ! ... First the I/O variable 69 ! 70 INTEGER, INTENT(IN) :: iter 71 !! counter of the number of iterations 72 INTEGER, INTENT(IN) :: n_iter 73 !! number of iterations used in mixing 74 INTEGER, INTENT(IN) :: iunmix 75 !! I/O unit where data from previous iterations is stored 76 REAL(DP), INTENT(IN) :: alphamix 77 !! mixing factor 78 REAL(DP), INTENT(IN) :: tr2_min 79 !! estimated error in diagonalization. If the estimated 80 !! scf error is smaller than this, exit: a more accurate 81 !! diagonalization is needed 82 REAL(DP), INTENT(OUT) :: dr2 83 !! the estimated error on the energy 84 LOGICAL, INTENT(OUT) :: conv 85 !! .TRUE. if the convergence has been reached 86 ! 87 TYPE(scf_type), INTENT(INOUT) :: input_rhout 88 TYPE(scf_type), INTENT(INOUT) :: rhoin 89 ! 90 ! ... local variables 91 ! 92 TYPE(mix_type) :: rhout_m, rhoin_m 93 INTEGER, PARAMETER :: & 94 maxmix = 25 ! max number of iterations for charge mixing 95 INTEGER :: & 96 iter_used, &! actual number of iterations used 97 ipos, &! index of the present iteration 98 inext, &! index of the next iteration 99 i, j, &! counters on number of iterations 100 info, &! flag saying if the exec. of libr. routines was ok 101 ldim, &! 2 * Hubbard_lmax + 1 102 iunmix_nsg, &! the unit for Hubbard mixing within DFT+U+V 103 nt ! index of the atomic type 104 REAL(DP),ALLOCATABLE :: betamix(:,:), work(:) 105 INTEGER, ALLOCATABLE :: iwork(:) 106 COMPLEX(DP), ALLOCATABLE :: nsginsave(:,:,:,:,:), nsgoutsave(:,:,:,:,:) 107 COMPLEX(DP), ALLOCATABLE :: deltansg(:,:,:,:,:) 108 LOGICAL :: exst 109 REAL(DP) :: gamma0 110#if defined(__NORMALIZE_BETAMIX) 111 REAL(DP) :: norm2, obn 112#endif 113 ! 114 ! ... saved variables and arrays 115 ! 116 INTEGER, SAVE :: & 117 mixrho_iter = 0 ! history of mixing 118#if !defined(__GFORTRAN_HACK) 119 TYPE(mix_type), ALLOCATABLE, SAVE :: & 120 df(:), &! information from preceding iterations 121 dv(:) ! " " " " " " 122#endif 123 REAL(DP) :: norm 124 INTEGER, PARAMETER :: read_ = -1, write_ = +1 125 ! 126 ! ... external functions 127 ! 128 INTEGER, EXTERNAL :: find_free_unit 129 ! 130 COMPLEX(DP), ALLOCATABLE, SAVE :: df_nsg(:,:,:,:,:,:), dv_nsg(:,:,:,:,:,:) 131 ! 132 CALL start_clock( 'mix_rho' ) 133 ! 134 ngm0 = ngms 135 ! 136 mixrho_iter = iter 137 ! 138 IF ( n_iter > maxmix ) CALL errore( 'mix_rho', 'n_iter too big', 1 ) 139 ! 140 ! define mix_type variables and copy scf_type variables there 141 ! 142 call create_mix_type(rhout_m) 143 call create_mix_type(rhoin_m) 144 ! 145 ! DFT+U+V case 146 ! 147 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 148 ldim = 0 149 DO nt = 1, ntyp 150 ldim = MAX(ldim, ldim_u(nt)) 151 ENDDO 152 ALLOCATE ( deltansg (ldim, ldim, max_num_neighbors, nat, nspin) ) 153 deltansg(:,:,:,:,:) = nsgnew(:,:,:,:,:) - nsg(:,:,:,:,:) 154 ENDIF 155 ! 156 call assign_scf_to_mix_type(rhoin, rhoin_m) 157 call assign_scf_to_mix_type(input_rhout, rhout_m) 158 159 call mix_type_AXPY ( -1.d0, rhoin_m, rhout_m ) 160 ! 161 dr2 = rho_ddot( rhout_m, rhout_m, ngms ) !!!! this used to be ngm NOT ngms 162 ! 163 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) & 164 dr2 = dr2 + nsg_ddot( deltansg, deltansg, nspin ) 165 ! 166 IF (dr2 < 0.0_DP) CALL errore('mix_rho','negative dr2',1) 167 ! 168 conv = ( dr2 < tr2 ) 169 ! 170 IF ( conv .OR. dr2 < tr2_min ) THEN 171 ! 172 ! ... if convergence is achieved or if the self-consistency error (dr2) is 173 ! ... smaller than the estimated error due to diagonalization (tr2_min), 174 ! ... exit and leave rhoin and rhocout unchanged 175 ! 176 IF ( ALLOCATED( df ) ) THEN 177 DO i=1, n_iter 178 call destroy_mix_type(df(i)) 179 END DO 180 DEALLOCATE( df ) 181 END IF 182 IF ( ALLOCATED( dv ) ) THEN 183 DO i=1, n_iter 184 call destroy_mix_type(dv(i)) 185 END DO 186 DEALLOCATE( dv ) 187 END IF 188 ! 189 call destroy_mix_type(rhoin_m) 190 call destroy_mix_type(rhout_m) 191 ! 192 IF ( ALLOCATED( dv_nsg ) ) DEALLOCATE( dv_nsg ) 193 IF ( ALLOCATED( df_nsg ) ) DEALLOCATE( df_nsg ) 194 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 195 nsgnew(:,:,:,:,:) = deltansg(:,:,:,:,:) + nsg(:,:,:,:,:) 196 DEALLOCATE(deltansg) 197 ENDIF 198 ! 199 CALL stop_clock( 'mix_rho' ) 200 ! 201 RETURN 202 ! 203 END IF 204 ! 205 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 206 iunmix_nsg = find_free_unit() 207 CALL diropn( iunmix_nsg, 'mix.nsg', ldim*ldim*nspin*nat*max_num_neighbors, exst) 208 ENDIF 209 ! 210 IF ( .NOT. ALLOCATED( df ) ) THEN 211 ALLOCATE( df( n_iter ) ) 212 DO i=1,n_iter 213 CALL create_mix_type( df(i) ) 214 END DO 215 END IF 216 IF ( .NOT. ALLOCATED( dv ) ) THEN 217 ALLOCATE( dv( n_iter ) ) 218 DO i=1,n_iter 219 CALL create_mix_type( dv(i) ) 220 END DO 221 END IF 222 ! 223 IF (lda_plus_u .AND. lda_plus_u_kind .EQ. 2) THEN 224 IF ( .NOT. ALLOCATED( df_nsg ) ) & 225 ALLOCATE( df_nsg(ldim,ldim,max_num_neighbors,nat,nspin,n_iter) ) 226 IF ( .NOT. ALLOCATED( dv_nsg ) ) & 227 ALLOCATE( dv_nsg(ldim,ldim,max_num_neighbors,nat,nspin,n_iter) ) 228 ENDIF 229 ! 230 ! ... iter_used = mixrho_iter-1 if mixrho_iter <= n_iter 231 ! ... iter_used = n_iter if mixrho_iter > n_iter 232 ! 233 iter_used = MIN( ( mixrho_iter - 1 ), n_iter ) 234 ! 235 ! ... ipos is the position in which results from the present iteration 236 ! ... are stored. ipos=mixrho_iter-1 until ipos=n_iter, then back to 1,2,... 237 ! 238 ipos = mixrho_iter - 1 - ( ( mixrho_iter - 2 ) / n_iter ) * n_iter 239 ! 240 IF ( mixrho_iter > 1 ) THEN 241 ! 242 CALL davcio_mix_type( df(ipos), iunmix, 1, read_ ) 243 CALL davcio_mix_type( dv(ipos), iunmix, 2, read_ ) 244 ! 245 call mix_type_AXPY ( -1.d0, rhout_m, df(ipos) ) 246 call mix_type_AXPY ( -1.d0, rhoin_m, dv(ipos) ) 247 ! 248 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 249 df_nsg(:,:,:,:,:,ipos) = df_nsg(:,:,:,:,:,ipos) - deltansg !nsgnew 250 dv_nsg(:,:,:,:,:,ipos) = dv_nsg(:,:,:,:,:,ipos) - nsg 251 ENDIF 252 ! 253#if defined(__NORMALIZE_BETAMIX) 254 ! NORMALIZE 255 norm2 = rho_ddot( df(ipos), df(ipos), ngm0 ) 256 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 257 norm2 = norm2 + nsg_ddot( df_nsg(1,1,1,1,1,ipos), df_nsg(1,1,1,1,1,ipos), nspin ) 258 ENDIF 259 obn = 1.d0/sqrt(norm2) 260 call mix_type_SCAL (obn,df(ipos)) 261 call mix_type_SCAL (obn,dv(ipos)) 262 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 263 df_nsg(:,:,:,:,:,ipos) = df_nsg(:,:,:,:,:,ipos) * obn 264 dv_nsg(:,:,:,:,:,ipos) = dv_nsg(:,:,:,:,:,ipos) * obn 265 ENDIF 266#endif 267 ! 268 END IF 269 ! 270 DO i = 1, iter_used 271 ! 272 IF ( i /= ipos ) THEN 273 ! 274 CALL davcio_mix_type( df(i), iunmix, 2*i+1, read_ ) 275 CALL davcio_mix_type( dv(i), iunmix, 2*i+2, read_ ) 276 END IF 277 ! 278 END DO 279 ! 280 CALL davcio_mix_type( rhout_m, iunmix, 1, write_ ) 281 CALL davcio_mix_type( rhoin_m, iunmix, 2, write_ ) 282 ! 283 IF ( mixrho_iter > 1 ) THEN 284 CALL davcio_mix_type( df(ipos), iunmix, 2*ipos+1, write_ ) 285 CALL davcio_mix_type( dv(ipos), iunmix, 2*ipos+2, write_ ) 286 END IF 287 ! 288 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 289 ! 290 ALLOCATE( nsginsave( ldim,ldim,max_num_neighbors,nat,nspin ), & 291 nsgoutsave( ldim,ldim,max_num_neighbors,nat,nspin ) ) 292 nsginsave = (0.d0,0.d0) 293 nsgoutsave = (0.d0,0.d0) 294 nsginsave = nsg 295 nsgoutsave = deltansg !nsgnew 296 ! 297 ENDIF 298 ! 299 ! Nothing else to do on first iteration 300 skip_on_first: & 301 IF (iter_used > 0) THEN 302 ! 303 ALLOCATE(betamix(iter_used, iter_used)) !iter_used)) 304 betamix = 0._dp 305 ! 306 DO i = 1, iter_used 307 ! 308 DO j = i, iter_used 309 ! 310 betamix(i,j) = rho_ddot( df(j), df(i), ngm0 ) 311 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) & 312 betamix(i,j) = betamix(i,j) + & 313 nsg_ddot( df_nsg(1,1,1,1,1,j), df_nsg(1,1,1,1,1,i), nspin ) 314 betamix(j,i) = betamix(i,j) 315 ! 316 END DO 317 ! 318 END DO 319 ! 320 ! allocate(e(iter_used), v(iter_used, iter_used)) 321 ! CALL rdiagh(iter_used, betamix, iter_used, e, v) 322 ! write(*,'(1e11.3)') e(:) 323 ! write(*,*) 324 ! deallocate(e,v) 325 allocate(work(iter_used), iwork(iter_used)) 326 !write(*,*) betamix(:,:) 327 CALL DSYTRF( 'U', iter_used, betamix, iter_used, iwork, work, iter_used, info ) 328 CALL errore( 'broyden', 'factorization', abs(info) ) 329 ! 330 CALL DSYTRI( 'U', iter_used, betamix, iter_used, iwork, work, info ) 331 CALL errore( 'broyden', 'DSYTRI', abs(info) ) ! 332 deallocate(iwork) 333 ! 334 FORALL( i = 1:iter_used, & 335 j = 1:iter_used, j > i ) betamix(j,i) = betamix(i,j) 336 ! 337 DO i = 1, iter_used 338 work(i) = rho_ddot( df(i), rhout_m, ngm0 ) 339 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) & 340 work(i) = work(i) + nsg_ddot( df_nsg(1,1,1,1,1,i), deltansg, nspin ) 341 END DO 342 ! 343 DO i = 1, iter_used 344 ! 345 gamma0 = DOT_PRODUCT( betamix(1:iter_used,i), work(1:iter_used) ) 346 ! 347 call mix_type_AXPY ( -gamma0, dv(i), rhoin_m ) 348 call mix_type_AXPY ( -gamma0, df(i), rhout_m ) 349 ! 350 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 351 nsg(:,:,:,:,:) = nsg(:,:,:,:,:) - gamma0*dv_nsg(:,:,:,:,:,i) 352 deltansg(:,:,:,:,:) = deltansg(:,:,:,:,:) - gamma0*df_nsg(:,:,:,:,:,i) 353 ENDIF 354 ! 355 END DO 356 DEALLOCATE(betamix, work) 357 ! 358 ! ... auxiliary vectors dv and df not needed anymore 359 ! 360 ENDIF skip_on_first 361 ! 362 IF ( ALLOCATED( df ) ) THEN 363 DO i=1, n_iter 364 call destroy_mix_type(df(i)) 365 END DO 366 DEALLOCATE( df ) 367 END IF 368 IF ( ALLOCATED( dv ) ) THEN 369 DO i=1, n_iter 370 call destroy_mix_type(dv(i)) 371 END DO 372 DEALLOCATE( dv ) 373 END IF 374 ! 375 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 376 inext = mixrho_iter - ( ( mixrho_iter - 1 ) / n_iter ) * n_iter 377 df_nsg(:,:,:,:,:,inext) = nsgoutsave 378 dv_nsg(:,:,:,:,:,inext) = nsginsave 379 IF (ALLOCATED(nsginsave)) DEALLOCATE(nsginsave) 380 IF (ALLOCATED(nsgoutsave)) DEALLOCATE(nsgoutsave) 381 ENDIF 382 ! 383 ! ... preconditioning the new search direction 384 ! 385 IF ( imix == 1 ) THEN 386 ! 387 CALL approx_screening( rhout_m ) 388 ! 389 ELSE IF ( imix == 2 ) THEN 390 ! 391 CALL approx_screening2( rhout_m, rhoin_m ) 392 ! 393 END IF 394 ! 395 ! ... set new trial density 396 ! 397 call mix_type_AXPY ( alphamix, rhout_m, rhoin_m ) 398 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 399 nsg = nsg + alphamix * deltansg 400 IF (ALLOCATED(deltansg)) DEALLOCATE(deltansg) 401 ENDIF 402 ! ... simple mixing for high_frequencies (and set to zero the smooth ones) 403 call high_frequency_mixing ( rhoin, input_rhout, alphamix ) 404 ! ... add the mixed rho for the smooth frequencies 405 call assign_mix_to_scf_type(rhoin_m,rhoin) 406 ! 407 call destroy_mix_type(rhout_m) 408 call destroy_mix_type(rhoin_m) 409 ! 410 IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) THEN 411 CLOSE( iunmix_nsg, STATUS = 'KEEP' ) 412 ENDIF 413 ! 414 CALL stop_clock( 'mix_rho' ) 415 ! 416 RETURN 417 ! 418END SUBROUTINE mix_rho 419! 420!---------------------------------------------------------------------------- 421SUBROUTINE approx_screening( drho ) 422 !---------------------------------------------------------------------------- 423 !! Apply an average TF preconditioning to \(\text{drho}\). 424 ! 425 USE kinds, ONLY : DP 426 USE constants, ONLY : e2, pi, fpi 427 USE cell_base, ONLY : omega, tpiba2 428 USE gvect, ONLY : gg, ngm 429 USE klist, ONLY : nelec 430 USE control_flags, ONLY : ngm0 431 USE scf, ONLY : mix_type 432 USE wavefunctions, ONLY : psic 433 ! 434 IMPLICIT NONE 435 ! 436 type (mix_type), intent(INOUT) :: drho ! (in/out) 437 ! 438 REAL(DP) :: rs, agg0 439 ! 440 rs = ( 3.D0 * omega / fpi / nelec )**( 1.D0 / 3.D0 ) 441 ! 442 agg0 = ( 12.D0 / pi )**( 2.D0 / 3.D0 ) / tpiba2 / rs 443 ! 444 drho%of_g(:ngm0,1) = drho%of_g(:ngm0,1) * gg(:ngm0) / (gg(:ngm0)+agg0) 445 ! 446 RETURN 447 ! 448END SUBROUTINE approx_screening 449! 450!---------------------------------------------------------------------------- 451SUBROUTINE approx_screening2( drho, rhobest ) 452 !---------------------------------------------------------------------------- 453 !! Apply a local-density dependent TF preconditioning to \(\text{drho}\). 454 ! 455 USE kinds, ONLY : DP 456 USE constants, ONLY : e2, pi, tpi, fpi, eps8, eps32 457 USE cell_base, ONLY : omega, tpiba2 458 USE gvect, ONLY : gg, ngm 459 USE wavefunctions, ONLY : psic 460 USE klist, ONLY : nelec 461 USE control_flags, ONLY : ngm0, gamma_only 462 USE scf, ONLY : mix_type, local_tf_ddot 463 USE mp, ONLY : mp_sum 464 USE mp_bands, ONLY : intra_bgrp_comm 465 USE fft_base, ONLY : dffts 466 USE fft_interfaces, ONLY : fwfft, invfft 467 ! 468 IMPLICIT NONE 469 ! 470 type(mix_type), intent(inout) :: drho 471 type(mix_type), intent(in) :: rhobest 472 ! 473 INTEGER, PARAMETER :: mmx = 12 474 ! 475 INTEGER :: & 476 iwork(mmx), i, j, m, info, is 477 REAL(DP) :: & 478 avg_rsm1, target, dr2_best 479 REAL(DP) :: & 480 aa(mmx,mmx), invaa(mmx,mmx), bb(mmx), work(mmx), vec(mmx), agg0 481 COMPLEX(DP), ALLOCATABLE :: & 482 v(:,:), &! v(ngm0,mmx) 483 w(:,:), &! w(ngm0,mmx) 484 dv(:), &! dv(ngm0) 485 vbest(:), &! vbest(ngm0) 486 wbest(:) ! wbest(ngm0) 487 REAL(DP), ALLOCATABLE :: & 488 alpha(:) ! alpha(dffts%nnr) 489 ! 490 INTEGER :: ir, ig 491 REAL(DP), PARAMETER :: one_third = 1.D0 / 3.D0 492 ! 493 target = 0.D0 494 ! 495 IF ( gg(1) < eps8 ) drho%of_g(1,1) = ZERO 496 ! 497 ALLOCATE( alpha( dffts%nnr ) ) 498 ALLOCATE( v( ngm0, mmx ), & 499 w( ngm0, mmx ), dv( ngm0 ), vbest( ngm0 ), wbest( ngm0 ) ) 500 ! 501 !$omp parallel 502 ! 503 CALL threaded_barrier_memset(psic, 0.0_DP, dffts%nnr*2) 504 !$omp do 505 DO ig = 1, ngm0 506 psic(dffts%nl(ig)) = rhobest%of_g(ig,1) 507 ENDDO 508 !$omp end do nowait 509 ! 510 !$omp end parallel 511 ! 512 ! ... calculate alpha from density 513 ! 514 CALL invfft ('Rho', psic, dffts) 515 ! 516 avg_rsm1 = 0.D0 517 ! 518 !$omp parallel do reduction(+:avg_rsm1) 519 DO ir = 1, dffts%nnr 520 alpha(ir) = ABS( REAL( psic(ir) ) ) 521 ! 522 IF ( alpha(ir) > eps32 ) THEN 523 ! 524 alpha(ir) = ( 3.D0 / fpi / alpha(ir) )**one_third 525 avg_rsm1 = avg_rsm1 + 1.D0 / alpha(ir) 526 ! 527 END IF 528 ! 529 alpha(ir) = 3.D0 * ( tpi / 3.D0 )**( 5.D0 / 3.D0 ) * alpha(ir) 530 ! 531 END DO 532 !$omp end parallel do 533 ! 534 CALL mp_sum( avg_rsm1 , intra_bgrp_comm ) 535 avg_rsm1 = ( dffts%nr1*dffts%nr2*dffts%nr3 ) / avg_rsm1 536 agg0 = ( 12.D0 / pi )**( 2.D0 / 3.D0 ) / tpiba2 / avg_rsm1 537 ! 538 ! ... calculate deltaV and the first correction vector 539 ! 540 !$omp parallel 541 CALL threaded_barrier_memset(psic, 0.0_DP, dffts%nnr*2) 542 !$omp do 543 DO ig = 1, ngm0 544 psic(dffts%nl(ig)) = drho%of_g(ig,1) 545 ENDDO 546 !$omp end do nowait 547 ! 548 IF ( gamma_only ) THEN 549 !$omp do 550 DO ig = 1, ngm0 551 psic(dffts%nlm(ig)) = CONJG( psic(dffts%nl(ig)) ) 552 ENDDO 553 !$omp end do nowait 554 ENDIF 555 !$omp end parallel 556 ! 557 CALL invfft ('Rho', psic, dffts) 558 ! 559 !$omp parallel do 560 DO ir = 1, dffts%nnr 561 psic(ir) = psic(ir) * alpha(ir) 562 ENDDO 563 !$omp end parallel do 564 ! 565 CALL fwfft ('Rho', psic, dffts) 566 ! 567 !$omp parallel do 568 DO ig = 1, ngm0 569 dv(ig) = psic(dffts%nl(ig)) * gg(ig) * tpiba2 570 v(ig,1)= psic(dffts%nl(ig)) * gg(ig) / ( gg(ig) + agg0 ) 571 ENDDO 572 !$omp end parallel do 573 ! 574 m = 1 575 aa(:,:) = 0.D0 576 bb(:) = 0.D0 577 ! 578 repeat_loop: DO 579 ! 580 ! ... generate the vector w 581 ! 582 !$omp parallel 583 CALL threaded_barrier_memset(psic, 0.0_DP, dffts%nnr*2) 584 !$omp do 585 DO ig = 1, ngm0 586 ! 587 w(ig,m) = fpi * e2 * v(ig,m) 588 ! 589 psic(dffts%nl(ig)) = v(ig,m) 590 ENDDO 591 !$omp end do nowait 592 ! 593 IF ( gamma_only ) THEN 594 !$omp do 595 DO ig = 1, ngm0 596 psic(dffts%nlm(ig)) = CONJG( psic(dffts%nl(ig)) ) 597 ENDDO 598 !$omp end do nowait 599 ENDIF 600 !$omp end parallel 601 ! 602 CALL invfft ('Rho', psic, dffts) 603 ! 604 !$omp parallel do 605 DO ir = 1, dffts%nnr 606 psic(ir) = psic(ir) * alpha(ir) 607 ENDDO 608 !$omp end parallel do 609 ! 610 CALL fwfft ('Rho', psic, dffts) 611 ! 612 !$omp parallel do 613 DO ig = 1, ngm0 614 w(ig,m) = w(ig,m) + gg(ig) * tpiba2 * psic(dffts%nl(ig)) 615 ENDDO 616 !$omp end parallel do 617 ! 618 ! ... build the linear system 619 ! 620 DO i = 1, m 621 ! 622 aa(i,m) = local_tf_ddot( w(1,i), w(1,m), ngm0) 623 ! 624 aa(m,i) = aa(i,m) 625 ! 626 END DO 627 ! 628 bb(m) = local_tf_ddot( w(1,m), dv, ngm0) 629 ! 630 ! ... solve it -> vec 631 ! 632 invaa = aa 633 ! 634 CALL DSYTRF( 'U', m, invaa, mmx, iwork, work, mmx, info ) 635 CALL errore( 'broyden', 'factorization', info ) 636 ! 637 CALL DSYTRI( 'U', m, invaa, mmx, iwork, work, info ) 638 CALL errore( 'broyden', 'DSYTRI', info ) 639 ! 640 FORALL( i = 1:m, j = 1:m, j > i ) invaa(j,i) = invaa(i,j) 641 ! 642 FORALL( i = 1:m ) vec(i) = SUM( invaa(i,:)*bb(:) ) 643 ! 644 !$omp parallel 645 !$omp do 646 DO ig = 1, ngm0 647 vbest(ig) = ZERO 648 wbest(ig) = dv(ig) 649 ENDDO 650 !$omp end do nowait 651 ! 652 DO i = 1, m 653 !$omp do 654 DO ig = 1, ngm0 655 vbest(ig) = vbest(ig) + vec(i) * v(ig,i) 656 wbest(ig) = wbest(ig) - vec(i) * w(ig,i) 657 ENDDO 658 !$omp end do nowait 659 END DO 660 !$omp end parallel 661 ! 662 dr2_best = local_tf_ddot( wbest, wbest, ngm0 ) 663 ! 664 IF ( target == 0.D0 ) target = MAX( 1.D-12, 1.D-6*dr2_best ) 665 ! 666 IF ( dr2_best < target ) THEN 667 ! 668 !$omp parallel 669 !$omp do 670 DO ig = 1, ngm0 671 drho%of_g(ig,1) = vbest(ig) 672 ENDDO 673 !$omp end do nowait 674 ! 675 !$omp end parallel 676 ! 677 DEALLOCATE( alpha, v, w, dv, vbest, wbest ) 678 ! 679 EXIT repeat_loop 680 ! 681 ELSE IF ( m >= mmx ) THEN 682 ! 683 m = 1 684 ! 685 !$omp parallel do 686 DO ig = 1, ngm0 687 v(ig,m) = vbest(ig) 688 ENDDO 689 !$omp end parallel do 690 aa(:,:) = 0.D0 691 bb(:) = 0.D0 692 ! 693 CYCLE repeat_loop 694 ! 695 END IF 696 ! 697 m = m + 1 698 ! 699 !$omp parallel do 700 DO ig = 1, ngm0 701 v(ig,m) = wbest(ig) / ( gg(ig) + agg0 ) 702 ENDDO 703 !$omp end parallel do 704 ! 705 END DO repeat_loop 706 ! 707 RETURN 708 ! 709END SUBROUTINE approx_screening2 710