1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7MODULE qs_gspace_mixing 8 9 USE cp_control_types, ONLY: dft_control_type 10 USE cp_para_types, ONLY: cp_para_env_type 11 USE dbcsr_api, ONLY: dbcsr_add,& 12 dbcsr_copy,& 13 dbcsr_p_type,& 14 dbcsr_set,& 15 dbcsr_type 16 USE kinds, ONLY: dp 17 USE mathlib, ONLY: invert_matrix 18 USE message_passing, ONLY: mp_sum 19 USE pw_env_types, ONLY: pw_env_get,& 20 pw_env_type 21 USE pw_methods, ONLY: pw_axpy,& 22 pw_copy,& 23 pw_integrate_function,& 24 pw_scale,& 25 pw_transfer,& 26 pw_zero 27 USE pw_pool_types, ONLY: pw_pool_create_pw,& 28 pw_pool_give_back_pw,& 29 pw_pool_type 30 USE pw_types, ONLY: COMPLEXDATA1D,& 31 RECIPROCALSPACE,& 32 pw_p_type 33 USE qs_density_mixing_types, ONLY: broyden_mixing_new_nr,& 34 broyden_mixing_nr,& 35 cp_1d_z_p_type,& 36 gspace_mixing_nr,& 37 mixing_storage_type,& 38 multisecant_mixing_nr,& 39 pulay_mixing_nr 40 USE qs_environment_types, ONLY: get_qs_env,& 41 qs_environment_type 42 USE qs_rho_atom_types, ONLY: rho_atom_type 43 USE qs_rho_types, ONLY: qs_rho_get,& 44 qs_rho_type 45 USE qs_scf_methods, ONLY: cp_sm_mix 46#include "./base/base_uses.f90" 47 48 IMPLICIT NONE 49 50 PRIVATE 51 52 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gspace_mixing' 53 54 PUBLIC :: gspace_mixing 55 56CONTAINS 57 58! ************************************************************************************************** 59!> \brief Driver for the g-space mixing, calls the proper routine given the 60!>requested method 61!> \param qs_env ... 62!> \param mixing_method ... 63!> \param mixing_store ... 64!> \param rho ... 65!> \param para_env ... 66!> \param iter_count ... 67!> \par History 68!> 05.2009 69!> 02.2015 changed input to make g-space mixing available in linear scaling 70!> SCF [Patrick Seewald] 71!> \author MI 72! ************************************************************************************************** 73 SUBROUTINE gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count) 74 TYPE(qs_environment_type), POINTER :: qs_env 75 INTEGER :: mixing_method 76 TYPE(mixing_storage_type), POINTER :: mixing_store 77 TYPE(qs_rho_type), POINTER :: rho 78 TYPE(cp_para_env_type), POINTER :: para_env 79 INTEGER :: iter_count 80 81 CHARACTER(len=*), PARAMETER :: routineN = 'gspace_mixing', routineP = moduleN//':'//routineN 82 83 INTEGER :: handle, i, iatom, ig, ispin, natom, ng, & 84 nimg, nspin 85 LOGICAL :: gapw 86 REAL(dp) :: alpha 87 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r 88 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp 89 TYPE(dft_control_type), POINTER :: dft_control 90 TYPE(pw_env_type), POINTER :: pw_env 91 TYPE(pw_p_type) :: rho_tmp 92 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r 93 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 94 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom 95 96 CALL timeset(routineN, handle) 97 98 CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env) 99 IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb .OR. & 100 dft_control%qs_control%semi_empirical)) THEN 101 102 CPASSERT(ASSOCIATED(mixing_store)) 103 NULLIFY (auxbas_pw_pool, rho_ao_kp, rho_atom, rho_g, rho_r, tot_rho_r) 104 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g, rho_r=rho_r, tot_rho_r=tot_rho_r) 105 106 nspin = SIZE(rho_g, 1) 107 nimg = dft_control%nimages 108 ng = SIZE(rho_g(1)%pw%pw_grid%gsq) 109 CPASSERT((ng == SIZE(mixing_store%rhoin(1)%cc))) 110 111 alpha = mixing_store%alpha 112 gapw = dft_control%qs_control%gapw 113 114 IF (nspin == 2) THEN 115 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool) 116 CALL pw_pool_create_pw(auxbas_pw_pool, & 117 rho_tmp%pw, & 118 use_data=COMPLEXDATA1D, & 119 in_space=RECIPROCALSPACE) 120 CALL pw_zero(rho_tmp%pw) 121 CALL pw_copy(rho_g(1)%pw, rho_tmp%pw) 122 CALL pw_axpy(rho_g(2)%pw, rho_g(1)%pw, 1.0_dp) 123 CALL pw_axpy(rho_tmp%pw, rho_g(2)%pw, -1.0_dp) 124 CALL pw_scale(rho_g(2)%pw, -1.0_dp) 125 END IF 126 127 IF (iter_count + 1 <= mixing_store%nskip_mixing) THEN 128 ! skip mixing 129 DO ispin = 1, nspin 130 DO ig = 1, ng 131 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 132 END DO 133 IF (mixing_store%gmix_p) THEN 134 DO i = 1, nimg 135 CALL dbcsr_copy(mixing_store%rho_ao_in(ispin, i)%matrix, rho_ao_kp(ispin, i)%matrix) 136 END DO 137 END IF 138 END DO 139 IF (mixing_store%gmix_p .AND. gapw) THEN 140 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom) 141 natom = SIZE(rho_atom) 142 DO ispin = 1, nspin 143 DO iatom = 1, natom 144 IF (mixing_store%paw(iatom)) THEN 145 mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef 146 mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef 147 END IF 148 END DO 149 END DO 150 END IF 151 152 mixing_store%iter_method = "NoMix" 153 IF (nspin == 2) THEN 154 CALL pw_axpy(rho_g(2)%pw, rho_g(1)%pw, 1.0_dp) 155 CALL pw_scale(rho_g(1)%pw, 0.5_dp) 156 CALL pw_axpy(rho_g(1)%pw, rho_g(2)%pw, -1.0_dp) 157 CALL pw_scale(rho_g(2)%pw, -1.0_dp) 158 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tmp%pw) 159 END IF 160 CALL timestop(handle) 161 RETURN 162 END IF 163 164 IF ((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) THEN 165 CALL gmix_potential_only(qs_env, mixing_store, rho, para_env) 166 mixing_store%iter_method = "Kerker" 167 ELSEIF (mixing_method == gspace_mixing_nr) THEN 168 CALL gmix_potential_only(qs_env, mixing_store, rho, para_env) 169 mixing_store%iter_method = "Kerker" 170 ELSEIF (mixing_method == pulay_mixing_nr) THEN 171 CALL pulay_mixing(qs_env, mixing_store, rho, para_env) 172 mixing_store%iter_method = "Pulay" 173 ELSEIF (mixing_method == broyden_mixing_nr) THEN 174 CALL broyden_mixing(qs_env, mixing_store, rho, para_env) 175 mixing_store%iter_method = "Broy." 176 ELSEIF (mixing_method == broyden_mixing_new_nr) THEN 177 CPASSERT(.NOT. gapw) 178 CALL broyden_mixing_new(mixing_store, rho, para_env) 179 mixing_store%iter_method = "Broy." 180 ELSEIF (mixing_method == multisecant_mixing_nr) THEN 181 CPASSERT(.NOT. gapw) 182 CALL multisecant_mixing(mixing_store, rho, para_env) 183 mixing_store%iter_method = "MSec." 184 END IF 185 186 IF (nspin == 2) THEN 187 CALL pw_axpy(rho_g(2)%pw, rho_g(1)%pw, 1.0_dp) 188 CALL pw_scale(rho_g(1)%pw, 0.5_dp) 189 CALL pw_axpy(rho_g(1)%pw, rho_g(2)%pw, -1.0_dp) 190 CALL pw_scale(rho_g(2)%pw, -1.0_dp) 191 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tmp%pw) 192 END IF 193 194 DO ispin = 1, nspin 195 CALL pw_transfer(rho_g(ispin)%pw, rho_r(ispin)%pw) 196 tot_rho_r(ispin) = pw_integrate_function(rho_r(ispin)%pw, isign=-1) 197 END DO 198 END IF 199 200 CALL timestop(handle) 201 202 END SUBROUTINE gspace_mixing 203 204! ************************************************************************************************** 205!> \brief G-space mixing performed via the Kerker damping on the density on the grid 206!> thus affecting only the caluclation of the potential, but not the denisity matrix 207!> \param qs_env ... 208!> \param mixing_store ... 209!> \param rho ... 210!> \param para_env ... 211!> \par History 212!> 02.2009 213!> \author MI 214! ************************************************************************************************** 215 216 SUBROUTINE gmix_potential_only(qs_env, mixing_store, rho, para_env) 217 218 TYPE(qs_environment_type), POINTER :: qs_env 219 TYPE(mixing_storage_type), POINTER :: mixing_store 220 TYPE(qs_rho_type), POINTER :: rho 221 TYPE(cp_para_env_type), POINTER :: para_env 222 223 CHARACTER(len=*), PARAMETER :: routineN = 'gmix_potential_only', & 224 routineP = moduleN//':'//routineN 225 226 COMPLEX(dp), DIMENSION(:), POINTER :: cc_new 227 INTEGER :: handle, iatom, ic, ig, ispin, natom, ng, & 228 nspin 229 LOGICAL :: gapw 230 REAL(dp) :: alpha, f_mix, tmp 231 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, rho_ao_kp 232 TYPE(dft_control_type), POINTER :: dft_control 233 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g 234 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom 235 236 CPASSERT(ASSOCIATED(mixing_store%rhoin)) 237 CPASSERT(ASSOCIATED(mixing_store%kerker_factor)) 238 239 CALL timeset(routineN, handle) 240 NULLIFY (cc_new, dft_control, rho_ao_kp, rho_atom, rho_g) 241 242 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s_kp=matrix_s) 243 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g) 244 245 nspin = SIZE(rho_g, 1) 246 ng = SIZE(rho_g(1)%pw%pw_grid%gsq) 247 248 gapw = dft_control%qs_control%gapw 249 alpha = mixing_store%alpha 250 251 DO ispin = 1, nspin 252 cc_new => rho_g(ispin)%pw%cc 253 DO ig = 1, mixing_store%ig_max ! ng 254 f_mix = mixing_store%alpha*mixing_store%kerker_factor(ig) 255 cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig) 256 mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig) 257 END DO 258 DO ig = mixing_store%ig_max + 1, ng 259 f_mix = mixing_store%alpha 260 cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig) 261 mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig) 262 END DO 263 264 IF (mixing_store%gmix_p) THEN 265 DO ic = 1, SIZE(rho_ao_kp, 2) 266 CALL cp_sm_mix(m1=rho_ao_kp(ispin, ic)%matrix, & 267 m2=mixing_store%rho_ao_in(ispin, ic)%matrix, & 268 p_mix=alpha, & 269 delta=tmp, & 270 para_env=para_env) 271 CALL dbcsr_copy(mixing_store%rho_ao_in(ispin, ic)%matrix, rho_ao_kp(ispin, ic)%matrix) 272 END DO 273 END IF 274 END DO 275 276 IF (mixing_store%gmix_p .AND. gapw) THEN 277 CALL get_qs_env(qs_env=qs_env, & 278 rho_atom_set=rho_atom) 279 natom = SIZE(rho_atom) 280 DO ispin = 1, nspin 281 DO iatom = 1, natom 282 IF (mixing_store%paw(iatom)) THEN 283 rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef + & 284 mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha) 285 rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef + & 286 mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha) 287 mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef 288 mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef 289 END IF 290 END DO 291 END DO 292 END IF 293 294 CALL timestop(handle) 295 296 END SUBROUTINE gmix_potential_only 297 298! ************************************************************************************************** 299!> \brief Pulay Mixing using as metrics for the residual the Kerer damping factor 300!> The mixing is applied directly on the density in reciprocal space, 301!> therefore it affects the potentials 302!> on the grid but not the density matrix 303!> \param qs_env ... 304!> \param mixing_store ... 305!> \param rho ... 306!> \param para_env ... 307!> \par History 308!> 03.2009 309!> \author MI 310! ************************************************************************************************** 311 312 SUBROUTINE pulay_mixing(qs_env, mixing_store, rho, para_env) 313 314 TYPE(qs_environment_type), POINTER :: qs_env 315 TYPE(mixing_storage_type), POINTER :: mixing_store 316 TYPE(qs_rho_type), POINTER :: rho 317 TYPE(cp_para_env_type), POINTER :: para_env 318 319 CHARACTER(len=*), PARAMETER :: routineN = 'pulay_mixing', routineP = moduleN//':'//routineN 320 321 COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: cc_mix 322 INTEGER :: handle, i, iatom, ib, ibb, ic, ig, & 323 ispin, jb, n1, n2, natom, nb, nb1, & 324 nbuffer, ng, nspin 325 LOGICAL :: gapw 326 REAL(dp) :: alpha_kerker, alpha_pulay, beta, f_mix, & 327 inv_err, norm, norm_c_inv, vol 328 REAL(dp), ALLOCATABLE, DIMENSION(:) :: alpha_c, ev 329 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, c, c_inv, cpc_h_mix, cpc_s_mix 330 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp 331 TYPE(dbcsr_type), POINTER :: rho_ao_mix 332 TYPE(dft_control_type), POINTER :: dft_control 333 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g 334 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom 335 336 CPASSERT(ASSOCIATED(mixing_store%res_buffer)) 337 CPASSERT(ASSOCIATED(mixing_store%rhoin_buffer)) 338 NULLIFY (dft_control, rho_ao_mix, rho_ao_kp, rho_atom, rho_g) 339 CALL timeset(routineN, handle) 340 341 CALL get_qs_env(qs_env, dft_control=dft_control) 342 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g) 343 nspin = SIZE(rho_g, 1) 344 ng = SIZE(mixing_store%res_buffer(1, 1)%cc) 345 vol = rho_g(1)%pw%pw_grid%vol 346 347 alpha_kerker = mixing_store%alpha 348 beta = mixing_store%pulay_beta 349 alpha_pulay = mixing_store%pulay_alpha 350 nbuffer = mixing_store%nbuffer 351 gapw = dft_control%qs_control%gapw 352 IF (gapw) THEN 353 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom) 354 natom = SIZE(rho_atom) 355 END IF 356 357 ib = MODULO(mixing_store%ncall, nbuffer) + 1 358 mixing_store%ncall = mixing_store%ncall + 1 359 nb = MIN(mixing_store%ncall, nbuffer) 360 ibb = MODULO(mixing_store%ncall, nbuffer) + 1 361 362 nb1 = nb + 1 363 ALLOCATE (a(nb1, nb1)) 364 a = 0.0_dp 365 ALLOCATE (b(nb1, nb1)) 366 b = 0.0_dp 367 ALLOCATE (c(nb, nb)) 368 c = 0.0_dp 369 ALLOCATE (c_inv(nb, nb)) 370 c_inv = 0.0_dp 371 ALLOCATE (alpha_c(nb)) 372 alpha_c = 0.0_dp 373 norm_c_inv = 0.0_dp 374 ALLOCATE (ev(nb1)) 375 ev = 0.0_dp 376 ALLOCATE (cc_mix(ng)) 377 378 DO ispin = 1, nspin 379 mixing_store%res_buffer(ib, ispin)%cc(:) = CMPLX(0._dp, 0._dp, KIND=dp) 380 norm = 0.0_dp 381 IF (nb == 1) mixing_store%rhoin_buffer(1, ispin)%cc = mixing_store%rhoin(ispin)%cc 382 DO ig = 1, ng 383 f_mix = mixing_store%kerker_factor(ig) 384 mixing_store%res_buffer(ib, ispin)%cc(ig) = f_mix*(rho_g(ispin)%pw%cc(ig) - & 385 mixing_store%rhoin_buffer(ib, ispin)%cc(ig)) 386 END DO 387 388 DO jb = 1, nb 389 mixing_store%pulay_matrix(jb, ib) = 0.0_dp 390 DO ig = 1, ng 391 f_mix = mixing_store%special_metric(ig) 392 mixing_store%pulay_matrix(jb, ib) = mixing_store%pulay_matrix(jb, ib) + & 393 f_mix*(REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp) & 394 *REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + & 395 AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))* & 396 AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))) 397 END DO 398 CALL mp_sum(mixing_store%pulay_matrix(jb, ib), para_env%group) 399 mixing_store%pulay_matrix(jb, ib) = mixing_store%pulay_matrix(jb, ib)*vol*vol 400 mixing_store%pulay_matrix(ib, jb) = mixing_store%pulay_matrix(jb, ib) 401 END DO 402 403 IF (nb == 1) THEN 404 DO ig = 1, ng 405 f_mix = alpha_kerker*mixing_store%kerker_factor(ig) 406 cc_mix(ig) = rho_g(ispin)%pw%cc(ig) - & 407 mixing_store%rhoin_buffer(ib, ispin)%cc(ig) 408 rho_g(ispin)%pw%cc(ig) = f_mix*cc_mix(ig) + & 409 mixing_store%rhoin_buffer(ib, ispin)%cc(ig) 410 mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 411 END DO 412 IF (mixing_store%gmix_p) THEN 413 DO ic = 1, SIZE(rho_ao_kp, 2) 414 CALL dbcsr_copy(mixing_store%rho_ao_res_buffer(ispin, ic, ib)%matrix, rho_ao_kp(ispin, ic)%matrix) 415 CALL dbcsr_add(matrix_a=mixing_store%rho_ao_res_buffer(ispin, ic, ib)%matrix, alpha_scalar=1.0_dp, & 416 matrix_b=mixing_store%rho_ao_in(ispin, ic)%matrix, beta_scalar=-1.0_dp) 417 418 CALL dbcsr_add(matrix_a=rho_ao_kp(ispin, ic)%matrix, alpha_scalar=alpha_kerker, & 419 matrix_b=mixing_store%rho_ao_in(ispin, ic)%matrix, & 420 beta_scalar=(1.0_dp - alpha_kerker)) 421 422 CALL dbcsr_copy(mixing_store%rho_ao_in_buffer(ispin, ic, ib)%matrix, mixing_store%rho_ao_in(ispin, ic)%matrix) 423 CALL dbcsr_copy(mixing_store%rho_ao_in_buffer(ispin, ic, ibb)%matrix, rho_ao_kp(ispin, ic)%matrix) 424 END DO 425 IF (gapw) THEN 426 DO iatom = 1, natom 427 IF (mixing_store%paw(iatom)) THEN 428 mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - & 429 mixing_store%cpc_h_in(iatom, ispin)%r_coef 430 mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - & 431 mixing_store%cpc_s_in(iatom, ispin)%r_coef 432 433 rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_h(ispin)%r_coef + & 434 (1.0_dp - alpha_kerker)*mixing_store%cpc_h_in(iatom, ispin)%r_coef 435 rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_s(ispin)%r_coef + & 436 (1.0_dp - alpha_kerker)*mixing_store%cpc_s_in(iatom, ispin)%r_coef 437 438 mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef 439 mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef 440 441 mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef 442 mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef 443 END IF 444 END DO 445 END IF 446 END IF 447 ELSE 448 b(1:nb, 1:nb) = mixing_store%pulay_matrix(1:nb, 1:nb) 449 c(1:nb, 1:nb) = b(1:nb, 1:nb) 450 b(nb1, 1:nb) = -1.0_dp 451 b(1:nb, nb1) = -1.0_dp 452 b(nb1, nb1) = 0.0_dp 453 454 CALL invert_matrix(c, c_inv, inv_err) 455 alpha_c = 0.0_dp 456 DO i = 1, nb 457 DO jb = 1, nb 458 alpha_c(i) = alpha_c(i) + c_inv(jb, i) 459 norm_c_inv = norm_c_inv + c_inv(jb, i) 460 END DO 461 END DO 462 alpha_c(1:nb) = alpha_c(1:nb)/norm_c_inv 463 cc_mix = CMPLX(0._dp, 0._dp, KIND=dp) 464 DO jb = 1, nb 465 DO ig = 1, ng 466 cc_mix(ig) = cc_mix(ig) + & 467 alpha_c(jb)*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) + & 468 mixing_store%pulay_beta*mixing_store%res_buffer(jb, ispin)%cc(ig)) 469 END DO 470 END DO 471 mixing_store%rhoin_buffer(ibb, ispin)%cc = CMPLX(0._dp, 0._dp, KIND=dp) 472 IF (alpha_pulay > 0.0_dp) THEN 473 DO ig = 1, ng 474 f_mix = alpha_pulay*mixing_store%kerker_factor(ig) 475 rho_g(ispin)%pw%cc(ig) = f_mix*rho_g(ispin)%pw%cc(ig) + & 476 (1.0_dp - f_mix)*cc_mix(ig) 477 mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 478 END DO 479 ELSE 480 DO ig = 1, ng 481 rho_g(ispin)%pw%cc(ig) = cc_mix(ig) 482 mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 483 END DO 484 END IF 485 486 IF (mixing_store%gmix_p) THEN 487 CALL dbcsr_set(mixing_store%rho_ao_mix, 0.0_dp) 488 rho_ao_mix => mixing_store%rho_ao_mix 489 DO ic = 1, SIZE(rho_ao_kp, 2) 490 CALL dbcsr_copy(mixing_store%rho_ao_res_buffer(ispin, ic, ib)%matrix, rho_ao_kp(ispin, ic)%matrix) 491 CALL dbcsr_add(matrix_a=mixing_store%rho_ao_res_buffer(ispin, ic, ib)%matrix, alpha_scalar=1.0_dp, & 492 matrix_b=mixing_store%rho_ao_in_buffer(ispin, ic, ib)%matrix, beta_scalar=-1.0_dp) 493 494 DO jb = 1, nb 495 CALL dbcsr_add(matrix_a=rho_ao_mix, alpha_scalar=1.0_dp, & 496 matrix_b=mixing_store%rho_ao_in_buffer(ispin, ic, jb)%matrix, & 497 beta_scalar=alpha_c(jb)) 498 CALL dbcsr_add(matrix_a=rho_ao_mix, alpha_scalar=1.0_dp, & 499 matrix_b=mixing_store%rho_ao_res_buffer(ispin, ic, jb)%matrix, & 500 beta_scalar=alpha_c(jb)*beta) 501 END DO 502 CALL dbcsr_add(matrix_a=rho_ao_kp(ispin, ic)%matrix, alpha_scalar=alpha_pulay, & 503 matrix_b=rho_ao_mix, beta_scalar=(1.0_dp - alpha_pulay)) 504 CALL dbcsr_copy(mixing_store%rho_ao_in_buffer(ispin, ic, ibb)%matrix, rho_ao_kp(ispin, ic)%matrix) 505 END DO 506 END IF 507 IF (mixing_store%gmix_p .AND. gapw) THEN 508 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom) 509 DO iatom = 1, natom 510 IF (mixing_store%paw(iatom)) THEN 511 n1 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 1) 512 n2 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 2) 513 ALLOCATE (cpc_h_mix(n1, n2)) 514 ALLOCATE (cpc_s_mix(n1, n2)) 515 516 mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - & 517 mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef 518 mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - & 519 mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef 520 cpc_h_mix = 0.0_dp 521 cpc_s_mix = 0.0_dp 522 DO jb = 1, nb 523 cpc_h_mix(:, :) = cpc_h_mix(:, :) + & 524 alpha_c(jb)*mixing_store%cpc_h_in_buffer(jb, iatom, ispin)%r_coef(:, :) + & 525 alpha_c(jb)*beta*mixing_store%cpc_h_res_buffer(jb, iatom, ispin)%r_coef(:, :) 526 cpc_s_mix(:, :) = cpc_s_mix(:, :) + & 527 alpha_c(jb)*mixing_store%cpc_s_in_buffer(jb, iatom, ispin)%r_coef(:, :) + & 528 alpha_c(jb)*beta*mixing_store%cpc_s_res_buffer(jb, iatom, ispin)%r_coef(:, :) 529 END DO 530 rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_h(ispin)%r_coef + & 531 (1.0_dp - alpha_pulay)*cpc_h_mix 532 rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_s(ispin)%r_coef + & 533 (1.0_dp - alpha_pulay)*cpc_s_mix 534 mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef 535 mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef 536 DEALLOCATE (cpc_h_mix) 537 DEALLOCATE (cpc_s_mix) 538 END IF 539 END DO 540 END IF 541 END IF ! nb==1 542 END DO ! ispin 543 544 DEALLOCATE (a) 545 DEALLOCATE (b) 546 DEALLOCATE (ev) 547 DEALLOCATE (cc_mix) 548 DEALLOCATE (c, c_inv, alpha_c) 549 550 CALL timestop(handle) 551 552 END SUBROUTINE pulay_mixing 553 554! ************************************************************************************************** 555!> \brief Broyden Mixing using as metrics for the residual the Kerer damping factor 556!> The mixing is applied directly on the density expansion in reciprocal space 557!> \param qs_env ... 558!> \param mixing_store ... 559!> \param rho ... 560!> \param para_env ... 561!> \par History 562!> 03.2009 563!> \author MI 564! ************************************************************************************************** 565 566 SUBROUTINE broyden_mixing(qs_env, mixing_store, rho, para_env) 567 568 TYPE(qs_environment_type), POINTER :: qs_env 569 TYPE(mixing_storage_type), POINTER :: mixing_store 570 TYPE(qs_rho_type), POINTER :: rho 571 TYPE(cp_para_env_type), POINTER :: para_env 572 573 CHARACTER(len=*), PARAMETER :: routineN = 'broyden_mixing', routineP = moduleN//':'//routineN 574 575 COMPLEX(dp) :: cc_mix 576 COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: res_rho 577 INTEGER :: handle, i, iatom, ib, ic, ig, ispin, j, & 578 jb, kb, n1, n2, natom, nb, nbuffer, & 579 ng, nspin 580 LOGICAL :: gapw, only_kerker 581 REAL(dp) :: alpha, dcpc_h_res, dcpc_s_res, & 582 delta_norm, f_mix, inv_err, norm, & 583 res_norm, valh, vals, w0 584 REAL(dp), ALLOCATABLE, DIMENSION(:) :: c, g 585 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b 586 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp 587 TYPE(dbcsr_type), POINTER :: rho_ao_mix, rho_ao_res 588 TYPE(dft_control_type), POINTER :: dft_control 589 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g 590 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom 591 592 CPASSERT(ASSOCIATED(mixing_store%res_buffer)) 593 CPASSERT(ASSOCIATED(mixing_store%rhoin)) 594 CPASSERT(ASSOCIATED(mixing_store%rhoin_old)) 595 CPASSERT(ASSOCIATED(mixing_store%drho_buffer)) 596 NULLIFY (dft_control, rho_ao_kp, rho_ao_mix, rho_ao_res, rho_atom, rho_g) 597 CALL timeset(routineN, handle) 598 599 CALL get_qs_env(qs_env, dft_control=dft_control) 600 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g) 601 602 nspin = SIZE(rho_g, 1) 603 ng = SIZE(mixing_store%res_buffer(1, 1)%cc) 604 605 alpha = mixing_store%alpha 606 w0 = mixing_store%broy_w0 607 nbuffer = mixing_store%nbuffer 608 gapw = dft_control%qs_control%gapw 609 610 ALLOCATE (res_rho(ng)) 611 612 mixing_store%ncall = mixing_store%ncall + 1 613 IF (mixing_store%ncall == 1) THEN 614 only_kerker = .TRUE. 615 ELSE 616 only_kerker = .FALSE. 617 nb = MIN(mixing_store%ncall - 1, nbuffer) 618 ib = MODULO(mixing_store%ncall - 2, nbuffer) + 1 619 END IF 620 621 IF (gapw) THEN 622 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom) 623 natom = SIZE(rho_atom) 624 ELSE 625 natom = 0 626 END IF 627 628 DO ispin = 1, nspin 629 res_rho = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 630 DO ig = 1, ng 631 res_rho(ig) = rho_g(ispin)%pw%cc(ig) - mixing_store%rhoin(ispin)%cc(ig) 632 END DO 633 634 IF (only_kerker) THEN 635 DO ig = 1, ng 636 mixing_store%last_res(ispin)%cc(ig) = res_rho(ig) 637 f_mix = alpha*mixing_store%kerker_factor(ig) 638 rho_g(ispin)%pw%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig) 639 mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) 640 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 641 END DO 642 643 IF (mixing_store%gmix_p) THEN 644 DO ic = 1, SIZE(rho_ao_kp, 2) 645 CALL dbcsr_copy(mixing_store%rho_ao_lastres(ispin, ic)%matrix, rho_ao_kp(ispin, ic)%matrix) 646 CALL dbcsr_add(matrix_a=mixing_store%rho_ao_lastres(ispin, ic)%matrix, alpha_scalar=1.0_dp, & 647 matrix_b=mixing_store%rho_ao_in(ispin, ic)%matrix, beta_scalar=-1.0_dp) 648 649 CALL dbcsr_add(matrix_a=rho_ao_kp(ispin, ic)%matrix, alpha_scalar=alpha, & 650 matrix_b=mixing_store%rho_ao_in(ispin, ic)%matrix, & 651 beta_scalar=(1.0_dp - alpha)) 652 653 CALL dbcsr_copy(mixing_store%rho_ao_in_old(ispin, ic)%matrix, mixing_store%rho_ao_in(ispin, ic)%matrix) 654 CALL dbcsr_copy(mixing_store%rho_ao_in(ispin, ic)%matrix, rho_ao_kp(ispin, ic)%matrix) 655 END DO 656 657 IF (gapw) THEN 658 DO iatom = 1, natom 659 IF (mixing_store%paw(iatom)) THEN 660 mixing_store%cpc_h_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - & 661 mixing_store%cpc_h_in(iatom, ispin)%r_coef 662 mixing_store%cpc_s_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - & 663 mixing_store%cpc_s_in(iatom, ispin)%r_coef 664 665 rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef + & 666 mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha) 667 rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef + & 668 mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha) 669 670 mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef 671 mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef 672 mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef 673 mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef 674 END IF 675 END DO 676 END IF 677 END IF 678 ELSE 679 680 ALLOCATE (a(nb, nb)) 681 a = 0.0_dp 682 ALLOCATE (b(nb, nb)) 683 b = 0.0_dp 684 ALLOCATE (c(nb)) 685 c = 0.0_dp 686 ALLOCATE (g(nb)) 687 g = 0.0_dp 688 689 delta_norm = 0.0_dp 690 res_norm = 0.0_dp 691 DO ig = 1, ng 692 mixing_store%res_buffer(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig) 693 mixing_store%last_res(ispin)%cc(ig) = res_rho(ig) 694 res_norm = res_norm + & 695 REAL(res_rho(ig), dp)*REAL(res_rho(ig), dp) + & 696 AIMAG(res_rho(ig))*AIMAG(res_rho(ig)) 697 delta_norm = delta_norm + & 698 REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp)* & 699 REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + & 700 AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))* & 701 AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig)) 702 END DO 703 DO ig = 1, ng 704 mixing_store%drho_buffer(ib, ispin)%cc(ig) = & 705 mixing_store%rhoin(ispin)%cc(ig) - & 706 mixing_store%rhoin_old(ispin)%cc(ig) 707 END DO 708 CALL mp_sum(delta_norm, para_env%group) 709 delta_norm = SQRT(delta_norm) 710 CALL mp_sum(res_norm, para_env%group) 711 res_norm = SQRT(res_norm) 712 mixing_store%res_buffer(ib, ispin)%cc(:) = mixing_store%res_buffer(ib, ispin)%cc(:)/delta_norm 713 mixing_store%drho_buffer(ib, ispin)%cc(:) = mixing_store%drho_buffer(ib, ispin)%cc(:)/delta_norm 714 715 IF (mixing_store%gmix_p .AND. gapw) THEN 716 DO iatom = 1, natom 717 IF (mixing_store%paw(iatom)) THEN 718 n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1) 719 n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2) 720 DO i = 1, n2 721 DO j = 1, n1 722 mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = & 723 (mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i) - & 724 mixing_store%cpc_h_old(iatom, ispin)%r_coef(j, i))/delta_norm 725 dcpc_h_res = ((rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - & 726 mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)) - & 727 mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i))/delta_norm 728 mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - & 729 mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i) 730 731 mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = & 732 (mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i) - & 733 mixing_store%cpc_s_old(iatom, ispin)%r_coef(j, i))/delta_norm 734 dcpc_s_res = ((rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - & 735 mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)) - & 736 mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i))/delta_norm 737 mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - & 738 mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i) 739 740 mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = & 741 alpha*dcpc_h_res + & 742 mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) 743 mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = & 744 alpha*dcpc_s_res + & 745 mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) 746 END DO 747 END DO 748 END IF 749 END DO 750 END IF 751 752 a(:, :) = 0.0_dp 753 DO ig = 1, ng 754 f_mix = alpha*mixing_store%kerker_factor(ig) 755 mixing_store%drho_buffer(ib, ispin)%cc(ig) = & 756 f_mix*mixing_store%res_buffer(ib, ispin)%cc(ig) + & 757 mixing_store%drho_buffer(ib, ispin)%cc(ig) 758 END DO 759 DO jb = 1, nb 760 DO kb = jb, nb 761 DO ig = 1, mixing_store%ig_max !ng 762 a(kb, jb) = a(kb, jb) + mixing_store%p_metric(ig)*( & 763 REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)* & 764 REAL(mixing_store%res_buffer(kb, ispin)%cc(ig), dp) + & 765 AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))* & 766 AIMAG(mixing_store%res_buffer(kb, ispin)%cc(ig))) 767 END DO 768 a(jb, kb) = a(kb, jb) 769 END DO 770 END DO 771 CALL mp_sum(a, para_env%group) 772 773 C = 0.0_dp 774 DO jb = 1, nb 775 a(jb, jb) = w0 + a(jb, jb) 776 DO ig = 1, mixing_store%ig_max !ng 777 c(jb) = c(jb) + mixing_store%p_metric(ig)*( & 778 REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)*REAL(res_rho(ig), dp) + & 779 AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))*AIMAG(res_rho(ig))) 780 END DO 781 END DO 782 CALL mp_sum(c, para_env%group) 783 CALL invert_matrix(a, b, inv_err) 784 785 CALL dgemv('T', nb, nb, 1.0_dp, B, nb, C, 1, 0.0_dp, G, 1) 786 787 DO ig = 1, ng 788 cc_mix = CMPLX(0.0_dp, 0.0_dp, kind=dp) 789 DO jb = 1, nb 790 cc_mix = cc_mix - G(jb)*mixing_store%drho_buffer(jb, ispin)%cc(ig) 791 END DO 792 f_mix = alpha*mixing_store%kerker_factor(ig) 793 rho_g(ispin)%pw%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) + & 794 f_mix*res_rho(ig) + cc_mix 795 mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) 796 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 797 END DO 798 799 IF (mixing_store%gmix_p) THEN 800 norm = 1.0_dp/delta_norm 801 f_mix = alpha*norm 802 rho_ao_mix => mixing_store%rho_ao_mix 803 rho_ao_res => mixing_store%rho_ao_res 804 DO ic = 1, SIZE(rho_ao_kp, 2) 805 CALL dbcsr_copy(rho_ao_res, rho_ao_kp(ispin, ic)%matrix) 806 CALL dbcsr_add(rho_ao_res, alpha_scalar=1.0_dp, & 807 matrix_b=mixing_store%rho_ao_in(ispin, ic)%matrix, beta_scalar=-1.0_dp) 808 CALL dbcsr_copy(rho_ao_mix, mixing_store%rho_ao_lastres(ispin, ic)%matrix) 809 CALL dbcsr_copy(mixing_store%rho_ao_lastres(ispin, ic)%matrix, rho_ao_res) 810 811 CALL dbcsr_add(rho_ao_res, alpha_scalar=f_mix, & 812 matrix_b=rho_ao_mix, beta_scalar=-f_mix) 813 814 CALL dbcsr_add(rho_ao_res, alpha_scalar=1.0_dp, & 815 matrix_b=mixing_store%rho_ao_in(ispin, ic)%matrix, beta_scalar=norm) 816 CALL dbcsr_add(rho_ao_res, alpha_scalar=1.0_dp, & 817 matrix_b=mixing_store%rho_ao_in_old(ispin, ic)%matrix, beta_scalar=-norm) 818 819 CALL dbcsr_copy(mixing_store%rho_ao_res_buffer(ispin, ic, ib)%matrix, rho_ao_res) 820 821 CALL dbcsr_set(mixing_store%rho_ao_mix, 0.0_dp) 822 DO jb = 1, nb 823 CALL dbcsr_add(matrix_a=rho_ao_mix, alpha_scalar=1.0_dp, & 824 matrix_b=mixing_store%rho_ao_res_buffer(ispin, ic, jb)%matrix, & 825 beta_scalar=g(jb)) 826 END DO 827 CALL dbcsr_add(matrix_a=rho_ao_kp(ispin, ic)%matrix, alpha_scalar=alpha, & 828 matrix_b=mixing_store%rho_ao_in(ispin, ic)%matrix, & 829 beta_scalar=(1.0_dp - alpha)) 830 831 CALL dbcsr_add(matrix_a=rho_ao_kp(ispin, ic)%matrix, alpha_scalar=1.0_dp, & 832 matrix_b=rho_ao_mix, beta_scalar=-1.0_dp) 833 834 CALL dbcsr_copy(mixing_store%rho_ao_in_old(ispin, ic)%matrix, mixing_store%rho_ao_in(ispin, ic)%matrix) 835 CALL dbcsr_copy(mixing_store%rho_ao_in(ispin, ic)%matrix, rho_ao_kp(ispin, ic)%matrix) 836 END DO 837 838 IF (gapw) THEN 839 DO iatom = 1, natom 840 IF (mixing_store%paw(iatom)) THEN 841 n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1) 842 n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2) 843 DO i = 1, n2 844 DO j = 1, n1 845 valh = 0.0_dp 846 vals = 0.0_dp 847 DO jb = 1, nb 848 valh = valh - G(jb)*mixing_store%dcpc_h_in(jb, iatom, ispin)%r_coef(j, i) 849 vals = vals - G(jb)*mixing_store%dcpc_s_in(jb, iatom, ispin)%r_coef(j, i) 850 END DO 851 rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) = & 852 alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) + & 853 mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + valh 854 rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) = & 855 alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) + & 856 mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + vals 857 END DO 858 END DO 859 860 mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef 861 mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef 862 mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef 863 mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef 864 END IF 865 END DO 866 END IF 867 END IF 868 869 DEALLOCATE (a, b, c, g) 870 END IF 871 872 END DO ! ispin 873 874 DEALLOCATE (res_rho) 875 876 CALL timestop(handle) 877 878 END SUBROUTINE broyden_mixing 879 880! ************************************************************************************************** 881!> \brief ... 882!> \param mixing_store ... 883!> \param rho ... 884!> \param para_env ... 885! ************************************************************************************************** 886 SUBROUTINE broyden_mixing_new(mixing_store, rho, para_env) 887 888 TYPE(mixing_storage_type), POINTER :: mixing_store 889 TYPE(qs_rho_type), POINTER :: rho 890 TYPE(cp_para_env_type), POINTER :: para_env 891 892 CHARACTER(len=*), PARAMETER :: routineN = 'broyden_mixing_new', & 893 routineP = moduleN//':'//routineN 894 895 COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: delta_res_p, res_rho, res_rho_p, tmp 896 INTEGER :: handle, ib, ibb, ig, info, ispin, jb, & 897 kb, kkb, lwork, nb, nbuffer, ng, nspin 898 INTEGER, ALLOCATABLE, DIMENSION(:) :: IWORK 899 LOGICAL :: only_kerker, skip_bq 900 REAL(dp) :: alpha, beta, delta, delta_p, delta_rhog, delta_rhog_p, f_mix, imp, imp_j, & 901 inv_err, norm, norm_ig, rep, rep_j, sqt_uvol, sqt_vol, uvol, vol, wc, wmax 902 REAL(dp), ALLOCATABLE, DIMENSION(:) :: aval, work_dgesdd 903 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, au, av, b, bq, work 904 REAL(dp), DIMENSION(:), POINTER :: p_metric 905 REAL(dp), DIMENSION(:, :), POINTER :: fmat, weight 906 TYPE(cp_1d_z_p_type), DIMENSION(:), POINTER :: tmp_z 907 TYPE(cp_1d_z_p_type), DIMENSION(:, :), POINTER :: delta_res, u_vec, z_vec 908 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g 909 910 CPASSERT(ASSOCIATED(mixing_store%rhoin_buffer)) 911 912 CALL timeset(routineN, handle) 913 914 NULLIFY (delta_res, u_vec, z_vec) 915 NULLIFY (fmat, rho_g) 916 CALL qs_rho_get(rho, rho_g=rho_g) 917 918 nspin = SIZE(rho_g, 1) 919 ng = SIZE(mixing_store%rhoin_buffer(1, 1)%cc) 920 vol = rho_g(1)%pw%pw_grid%vol 921 sqt_vol = SQRT(vol) 922 uvol = 1.0_dp/vol 923 sqt_uvol = SQRT(uvol) 924 alpha = mixing_store%alpha 925 beta = mixing_store%beta 926 927 wc = mixing_store%wc 928 wmax = mixing_store%wmax 929 nbuffer = mixing_store%nbuffer 930 931 mixing_store%ncall = mixing_store%ncall + 1 932 IF (mixing_store%ncall == 1) THEN 933 only_kerker = .TRUE. 934 ELSE 935 only_kerker = .FALSE. 936 nb = MIN(mixing_store%ncall - 1, nbuffer) 937 ib = MODULO(mixing_store%ncall - 2, nbuffer) + 1 938 939 ALLOCATE (a(nb, nb)) 940 a = 0.0_dp 941 ALLOCATE (b(nb, nb)) 942 b = 0.0_dp 943 ALLOCATE (bq(nb, nb)) 944 bq = 0.0_dp 945 946 ALLOCATE (tmp(ng)) 947 ALLOCATE (delta_res_p(ng)) 948 END IF 949 950 ibb = 0 951 952 ALLOCATE (res_rho(ng)) 953 ALLOCATE (res_rho_p(ng)) 954 955 p_metric => mixing_store%p_metric 956 weight => mixing_store%weight 957 CPASSERT(ASSOCIATED(mixing_store%delta_res)) 958 delta_res => mixing_store%delta_res 959 CPASSERT(ASSOCIATED(mixing_store%u_vec)) 960 u_vec => mixing_store%u_vec 961 CPASSERT(ASSOCIATED(mixing_store%z_vec)) 962 z_vec => mixing_store%z_vec 963 964 delta_rhog = 0.0_dp 965 delta_rhog_p = 0.0_dp 966 967 DO ispin = 1, nspin 968 969 fmat => mixing_store%fmat(:, :, ispin) 970 971 delta = 0.0_dp 972 delta_p = 0.0_dp 973 ! Residual at this step R_i(G) (rho_out(G)-rho_in(G)) 974 ! Residual multiplied by the metrics RP_i(G) = (rho_out(G)-rho_in(G)) * P(G) 975 ! Delta is the norm of the residual, measures how far we are from convergence 976 DO ig = 1, ng 977 res_rho(ig) = rho_g(ispin)%pw%cc(ig) - mixing_store%rhoin(ispin)%cc(ig) 978 res_rho_p(ig) = res_rho(ig)*p_metric(ig) !*sqt_uvol 979 norm_ig = REAL(res_rho(ig), dp)*REAL(res_rho(ig), dp) + AIMAG(res_rho(ig))*AIMAG(res_rho(ig)) 980 delta = delta + norm_ig 981 delta_p = delta_p + norm_ig*p_metric(ig) !*p_metric(ig) 982 END DO 983 CALL mp_sum(delta, para_env%group) 984 delta = SQRT(delta) 985 CALL mp_sum(delta_p, para_env%group) 986 delta_p = SQRT(delta_p) 987 delta_rhog = delta_rhog + delta 988 delta_rhog_p = delta_rhog_p + delta_p 989 990 weight(ib, ispin) = 1.0_dp ! wc 991 IF (wc < 0.0_dp) weight(ib, ispin) = 0.01_dp*ABS(wc)/(delta_p*delta_p) 992 IF (weight(ib, ispin) == 0.0_dp) weight(ib, ispin) = 100.0_dp 993 IF (weight(ib, ispin) < 1.0_dp) weight(ib, ispin) = 1.0_dp 994 995 IF (only_kerker) THEN 996 ! Simple Kerker damping : linear mixing rho(G) = rho_in(G) - alpha k(G)*(rho_out(G)-rho_in(G)) 997 DO ig = 1, ng 998 f_mix = alpha*mixing_store%kerker_factor(ig) 999 rho_g(ispin)%pw%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig) 1000 mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) 1001 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 1002 mixing_store%last_res(ispin)%cc(ig) = res_rho(ig) 1003 END DO 1004 ELSE 1005 norm = 0.0_dp 1006 ! Difference of residuals DR_{i-1)} (G) = R_i(G) - R_{i-1}(G) 1007 DO ig = 1, ng 1008 delta_res(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig) 1009 delta_res_p(ig) = p_metric(ig)*delta_res(ib, ispin)%cc(ig) 1010 norm_ig = REAL(delta_res(ib, ispin)%cc(ig), dp)*REAL(delta_res_p(ig), dp) + & 1011 AIMAG(delta_res(ib, ispin)%cc(ig))*AIMAG(delta_res_p(ig)) 1012 norm = norm + norm_ig 1013 END DO 1014 CALL mp_sum(norm, para_env%group) 1015 norm = 1._dp/SQRT(norm) 1016 delta_res(ib, ispin)%cc(:) = delta_res(ib, ispin)%cc(:)*norm 1017 delta_res_p(:) = delta_res_p(:)*norm 1018 1019 IF (para_env%ionode) WRITE (*, *) 'norm ', norm 1020 ! Vector U_{i-1}(G) = Drho_{i-1} + k(G) * DR_{i-1}(G) 1021 DO ig = 1, ng 1022 tmp(ig) = (mixing_store%rhoin(ispin)%cc(ig) - & 1023 mixing_store%rhoin_old(ispin)%cc(ig))*norm 1024 u_vec(ib, ispin)%cc(ig) = (tmp(ig) + & 1025 mixing_store%kerker_factor(ig)*delta_res(ib, ispin)%cc(ig)) 1026 END DO 1027 1028 DO jb = 1, nb 1029 fmat(jb, ib) = 0.0_dp 1030 1031 DO ig = 1, ng 1032 rep_j = REAL(delta_res(jb, ispin)%cc(ig), dp) 1033 imp_j = AIMAG(delta_res(jb, ispin)%cc(ig)) 1034 ! < DR_{j} | DR_{i-1} > 1035 rep = REAL(delta_res_p(ig), dp) 1036 imp = AIMAG(delta_res_p(ig)) 1037 fmat(jb, ib) = fmat(jb, ib) + rep_j*rep + imp_j*imp 1038 END DO 1039 1040 END DO 1041 CALL mp_sum(fmat(1:nb, ib), para_env%group) 1042 1043 fmat(ib, ib) = 1.0_dp 1044 1045 DO jb = 1, nb 1046 fmat(ib, jb) = fmat(jb, ib) 1047 ENDDO 1048 1049 DO jb = 1, nb 1050 a(jb, jb) = weight(jb, ispin)*weight(jb, ispin)*fmat(jb, jb) 1051 DO kb = 1, jb - 1 1052 a(jb, kb) = weight(jb, ispin)*weight(kb, ispin)*fmat(jb, kb) 1053 a(kb, jb) = weight(jb, ispin)*weight(kb, ispin)*fmat(kb, jb) 1054 ENDDO 1055 END DO 1056 IF (.TRUE.) THEN 1057 b = 0.0_dp 1058 CALL invert_matrix(a, b, inv_err) 1059 ELSE 1060 b = 0.0_dp 1061 ALLOCATE (work(ib - 1, ib - 1), aval(ib - 1), av(ib - 1, ib - 1), au(ib - 1, ib - 1)) 1062 work(:, :) = a 1063 ALLOCATE (iwork(8*(ib - 1)), work_dgesdd(1)) 1064 lwork = -1 1065 CALL DGESDD('S', ib - 1, ib - 1, work, ib - 1, aval, au, ib - 1, av, ib - 1, work_dgesdd, lwork, iwork, info) 1066 lwork = INT(work_dgesdd(1)) 1067 DEALLOCATE (work_dgesdd); ALLOCATE (work_dgesdd(lwork)) 1068 CALL DGESDD('S', ib - 1, ib - 1, work, ib - 1, aval, au, ib - 1, av, ib - 1, work_dgesdd, lwork, iwork, info) 1069 ! construct the inverse 1070 DO kb = 1, ib - 1 1071 ! invert SV 1072 IF (aval(kb) < 1.E-6_dp) THEN 1073 aval(kb) = 0.0_dp 1074 ELSE 1075 aval(kb) = 1.0_dp/aval(kb) 1076 ENDIF 1077 av(kb, :) = av(kb, :)*aval(kb) 1078 ENDDO 1079 CALL DGEMM('T', 'T', ib - 1, ib - 1, ib - 1, 1.0_dp, av, ib - 1, au, ib - 1, 0.0_dp, b, ib - 1) 1080 DEALLOCATE (iwork, aval, au, av, work_dgesdd, work) 1081 END IF 1082 1083 bq = 0.0_dp 1084 ! Broyden second method requires also bq (NYI) 1085 skip_bq = .TRUE. 1086 DO jb = 1, nb 1087 DO kb = 1, nb 1088 bq(jb, kb) = 0.0_dp 1089 DO kkb = 1, nb 1090 bq(jb, kb) = bq(jb, kb) - weight(jb, ispin)*weight(kkb, ispin)*b(jb, kkb)*fmat(kkb, kb) 1091 END DO 1092 END DO 1093 bq(jb, jb) = 1.0_dp + bq(jb, jb) 1094 END DO 1095 1096 IF (.NOT. skip_bq) THEN 1097 ! in this case the old z_vec is needed 1098 ! a temporary array is needed to store the new one 1099 ALLOCATE (tmp_z(nb)) 1100 DO jb = 1, nb 1101 ALLOCATE (tmp_z(jb)%cc(ng)) 1102 END DO 1103 END IF 1104 DO jb = 1, nb 1105 tmp(:) = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1106 IF (.NOT. skip_bq) THEN 1107 ! sum_{kb} bq(jb,kb) * z_vec_{kb,iter-2} 1108 ! added only for small weights 1109 DO kb = 1, nb 1110 IF (weight(kb, ispin) >= (10.0_dp*wmax)) CYCLE 1111 DO ig = 1, ng 1112 tmp(ig) = tmp(ig) + bq(jb, kb)*z_vec(kb, ispin)%cc(ig) 1113 END DO 1114 END DO ! kb 1115 END IF 1116 1117 ! sum_{kb} w(jb)*w(kb)*b(jb,kb) * u_vec_{kb} 1118 DO kb = 1, ib - 1 1119 DO ig = 1, ng 1120 tmp(ig) = tmp(ig) + weight(kb, ispin)*weight(jb, ispin)*b(jb, kb)*u_vec(kb, ispin)%cc(ig) 1121 END DO 1122 END DO 1123 1124 ! store the new z_vec(jb) 1125 IF (skip_bq .OR. (weight(jb, ispin) >= (10._dp*wmax))) THEN 1126 z_vec(jb, ispin)%cc(:) = tmp(:) 1127 ELSE 1128 ! temporary array: old z_vec may still be needed 1129 tmp_z(jb)%cc(:) = tmp(:) 1130 END IF 1131 END DO !jb 1132 1133 IF (.NOT. skip_bq) THEN 1134 DO jb = 1, ib - 1 1135 IF (weight(jb, ispin) < (10._dp*wmax)) z_vec(jb, ispin)%cc(:) = tmp_z(jb)%cc(:) 1136 DEALLOCATE (tmp_z(jb)%cc) 1137 END DO 1138 DEALLOCATE (tmp_z) 1139 END IF 1140 1141 ! Overwrite the density i reciprocal space 1142 rho_g(ispin)%pw%cc(:) = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1143 DO jb = 1, ib - 1 1144 norm = 0.0_dp 1145 DO ig = 1, ng 1146 rep_j = REAL(delta_res(jb, ispin)%cc(ig), dp) 1147 imp_j = AIMAG(delta_res(jb, ispin)%cc(ig)) 1148 rep = REAL(res_rho_p(ig), dp) 1149 imp = AIMAG(res_rho_p(ig)) 1150 norm = norm + rep_j*rep + imp_j*imp 1151 END DO 1152 CALL mp_sum(norm, para_env%group) 1153 ! Subtract |Z_jb)><DR_jb|P|R_{iter}> 1154 DO ig = 1, ng 1155 rho_g(ispin)%pw%cc(ig) = rho_g(ispin)%pw%cc(ig) - norm*z_vec(jb, ispin)%cc(ig)*sqt_vol 1156 1157 END DO 1158 END DO 1159 1160 DO ig = 1, ng 1161 f_mix = alpha*mixing_store%kerker_factor(ig) 1162 rho_g(ispin)%pw%cc(ig) = rho_g(ispin)%pw%cc(ig) + & 1163 mixing_store%rhoin_buffer(ib, ispin)%cc(ig) + f_mix*res_rho(ig) 1164 mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 1165 END DO 1166 1167 mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) 1168 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 1169 mixing_store%last_res(ispin)%cc(:) = res_rho(:) 1170 END IF ! ib 1171 1172 END DO ! ispin 1173 IF (.NOT. only_kerker) THEN 1174 DEALLOCATE (a, b, bq) 1175 DEALLOCATE (delta_res_p, tmp) 1176 END IF 1177 DEALLOCATE (res_rho, res_rho_p) 1178 1179 CALL timestop(handle) 1180 1181 END SUBROUTINE broyden_mixing_new 1182 1183! ************************************************************************************************** 1184!> \brief Multisecant scheme to perform charge density Mixing 1185!> as preconditioner we use the Kerker damping factor 1186!> The mixing is applied directly on the density in reciprocal space, 1187!> therefore it affects the potentials 1188!> on the grid but not the density matrix 1189!> \param mixing_store ... 1190!> \param rho ... 1191!> \param para_env ... 1192!> \par History 1193!> 05.2009 1194!> \author MI 1195! ************************************************************************************************** 1196 SUBROUTINE multisecant_mixing(mixing_store, rho, para_env) 1197 1198 TYPE(mixing_storage_type), POINTER :: mixing_store 1199 TYPE(qs_rho_type), POINTER :: rho 1200 TYPE(cp_para_env_type), POINTER :: para_env 1201 1202 CHARACTER(len=*), PARAMETER :: routineN = 'multisecant_mixing', & 1203 routineP = moduleN//':'//routineN 1204 COMPLEX(KIND=dp), PARAMETER :: cmone = (-1.0_dp, 0.0_dp), & 1205 cone = (1.0_dp, 0.0_dp), & 1206 czero = (0.0_dp, 0.0_dp) 1207 1208 COMPLEX(dp) :: saa, yaa 1209 COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: gn_global, pgn, res_matrix_global, & 1210 tmp_vec, ugn 1211 COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, res_matrix, sa, step_matrix, ya 1212 COMPLEX(dp), DIMENSION(:), POINTER :: gn 1213 INTEGER :: handle, ib, ib_next, ib_prev, ig, & 1214 ig_global, iig, ispin, jb, kb, nb, & 1215 nbuffer, ng, ng_global, nspin 1216 LOGICAL :: use_zgemm, use_zgemm_rev 1217 REAL(dp) :: alpha, f_mix, gn_norm, gn_norm_old, inv_err, n_low, n_up, pgn_norm, prec, & 1218 r_step, reg_par, sigma_max, sigma_tilde, step_size 1219 REAL(dp), ALLOCATABLE, DIMENSION(:) :: norm_res, norm_res_low, norm_res_up 1220 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b_matrix, binv_matrix 1221 REAL(dp), DIMENSION(:), POINTER :: g2 1222 REAL(dp), SAVE :: sigma_old = 1.0_dp 1223 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g 1224 1225 CPASSERT(ASSOCIATED(mixing_store)) 1226 1227 CALL timeset(routineN, handle) 1228 1229 NULLIFY (gn, rho_g) 1230 1231 use_zgemm = .FALSE. 1232 use_zgemm_rev = .TRUE. 1233 1234 ! prepare the parameters 1235 CALL qs_rho_get(rho, rho_g=rho_g) 1236 1237 nspin = SIZE(rho_g, 1) 1238 ! not implemented for large grids. 1239 CPASSERT(rho_g(1)%pw%pw_grid%ngpts < HUGE(ng_global)) 1240 ng_global = INT(rho_g(1)%pw%pw_grid%ngpts) 1241 ng = SIZE(mixing_store%rhoin_buffer(1, 1)%cc) 1242 alpha = mixing_store%alpha 1243 1244 sigma_max = mixing_store%sigma_max 1245 reg_par = mixing_store%reg_par 1246 r_step = mixing_store%r_step 1247 nbuffer = mixing_store%nbuffer 1248 1249 ! determine the step number, and multisecant iteration 1250 nb = MIN(mixing_store%ncall, nbuffer - 1) 1251 ib = MODULO(mixing_store%ncall, nbuffer) + 1 1252 IF (mixing_store%ncall > 0) THEN 1253 ib_prev = MODULO(mixing_store%ncall - 1, nbuffer) + 1 1254 ELSE 1255 ib_prev = 0 1256 END IF 1257 mixing_store%ncall = mixing_store%ncall + 1 1258 ib_next = MODULO(mixing_store%ncall, nbuffer) + 1 1259 1260 ! compute the residual gn and its norm gn_norm 1261 DO ispin = 1, nspin 1262 gn => mixing_store%res_buffer(ib, ispin)%cc 1263 gn_norm = 0.0_dp 1264 DO ig = 1, ng 1265 gn(ig) = (rho_g(ispin)%pw%cc(ig) - mixing_store%rhoin_buffer(ib, ispin)%cc(ig)) 1266 gn_norm = gn_norm + & 1267 REAL(gn(ig), dp)*REAL(gn(ig), dp) + AIMAG(gn(ig))*AIMAG(gn(ig)) 1268 END DO 1269 CALL mp_sum(gn_norm, para_env%group) 1270 gn_norm = SQRT(gn_norm) 1271 mixing_store%norm_res_buffer(ib, ispin) = gn_norm 1272 END DO 1273 1274 IF (nb == 0) THEN 1275 !simple mixing 1276 DO ispin = 1, nspin 1277 DO ig = 1, ng 1278 f_mix = alpha*mixing_store%kerker_factor(ig) 1279 rho_g(ispin)%pw%cc(ig) = mixing_store%rhoin_buffer(1, ispin)%cc(ig) + & 1280 f_mix*mixing_store%res_buffer(1, ispin)%cc(ig) 1281 mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 1282 END DO 1283 END DO 1284 CALL timestop(handle) 1285 RETURN 1286 END IF 1287 1288 ! allocate temporary arrays 1289 ! step_matrix S ngxnb 1290 ALLOCATE (step_matrix(ng, nb)) 1291 ! res_matrix Y ngxnb 1292 ALLOCATE (res_matrix(ng, nb)) 1293 ! matrix A nbxnb 1294 ALLOCATE (a_matrix(nb, ng_global)) 1295 ! PSI nb vector of norms 1296 ALLOCATE (norm_res(nb)) 1297 ALLOCATE (norm_res_low(nb)) 1298 ALLOCATE (norm_res_up(nb)) 1299 1300 ! matrix B nbxnb 1301 ALLOCATE (b_matrix(nb, nb)) 1302 ! matrix B_inv nbxnb 1303 ALLOCATE (binv_matrix(nb, nb)) 1304 1305 ALLOCATE (gn_global(ng_global)) 1306 ALLOCATE (res_matrix_global(ng_global)) 1307 IF (use_zgemm) THEN 1308 ALLOCATE (sa(ng, ng_global)) 1309 ALLOCATE (ya(ng, ng_global)) 1310 END IF 1311 IF (use_zgemm_rev) THEN 1312 ALLOCATE (tmp_vec(nb)) 1313 END IF 1314 ALLOCATE (pgn(ng)) 1315 ALLOCATE (ugn(ng)) 1316 1317 DO ispin = 1, nspin 1318 ! generate the global vector with the present residual 1319 gn => mixing_store%res_buffer(ib, ispin)%cc 1320 gn_global = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1321 DO ig = 1, ng 1322 ig_global = mixing_store%ig_global_index(ig) 1323 gn_global(ig_global) = gn(ig) 1324 END DO 1325 CALL mp_sum(gn_global, para_env%group) 1326 1327 ! compute steps (matrix S) and residual differences (matrix Y) 1328 ! with respect to the present 1329 ! step and the present residual (use stored rho_in and res_buffer) 1330 1331 ! These quantities are pre-conditioned by means of Kerker factor multipliccation 1332 g2 => rho_g(1)%pw%pw_grid%gsq 1333 1334 DO jb = 1, nb + 1 1335 IF (jb < ib) THEN 1336 kb = jb 1337 ELSEIF (jb > ib) THEN 1338 kb = jb - 1 1339 ELSE 1340 CYCLE 1341 END IF 1342 norm_res(kb) = 0.0_dp 1343 norm_res_low(kb) = 0.0_dp 1344 norm_res_up(kb) = 0.0_dp 1345 DO ig = 1, ng 1346 1347 prec = 1.0_dp 1348 1349 step_matrix(ig, kb) = prec*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) - & 1350 mixing_store%rhoin_buffer(ib, ispin)%cc(ig)) 1351 res_matrix(ig, kb) = (mixing_store%res_buffer(jb, ispin)%cc(ig) - & 1352 mixing_store%res_buffer(ib, ispin)%cc(ig)) 1353 norm_res(kb) = norm_res(kb) + REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + & 1354 AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb)) 1355 IF (g2(ig) < 4.0_dp) THEN 1356 norm_res_low(kb) = norm_res_low(kb) + & 1357 REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + & 1358 AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb)) 1359 ELSE 1360 norm_res_up(kb) = norm_res_up(kb) + & 1361 REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + & 1362 AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb)) 1363 END IF 1364 res_matrix(ig, kb) = prec*res_matrix(ig, kb) 1365 END DO 1366 END DO !jb 1367 1368 ! normalize each column of S and Y => Snorm Ynorm 1369 CALL mp_sum(norm_res, para_env%group) 1370 CALL mp_sum(norm_res_up, para_env%group) 1371 CALL mp_sum(norm_res_low, para_env%group) 1372 norm_res(1:nb) = 1.0_dp/SQRT(norm_res(1:nb)) 1373 n_low = 0.0_dp 1374 n_up = 0.0_dp 1375 DO jb = 1, nb 1376 n_low = n_low + norm_res_low(jb)/norm_res(jb) 1377 n_up = n_up + norm_res_up(jb)/norm_res(jb) 1378 END DO 1379 DO ig = 1, ng 1380 IF (g2(ig) > 4.0_dp) THEN 1381 step_matrix(ig, 1:nb) = step_matrix(ig, 1:nb)*SQRT(n_low/n_up) 1382 res_matrix(ig, 1:nb) = res_matrix(ig, 1:nb)*SQRT(n_low/n_up) 1383 END IF 1384 END DO 1385 DO kb = 1, nb 1386 step_matrix(1:ng, kb) = step_matrix(1:ng, kb)*norm_res(kb) 1387 res_matrix(1:ng, kb) = res_matrix(1:ng, kb)*norm_res(kb) 1388 END DO 1389 1390 ! compute A as [(Ynorm^t Ynorm) + (alpha I)]^(-1) Ynorm^t 1391 ! compute B 1392 DO jb = 1, nb 1393 DO kb = 1, nb 1394 b_matrix(kb, jb) = 0.0_dp 1395 DO ig = 1, ng 1396 ! it is assumed that summing over all G vector gives a real, because 1397 ! y(-G,kb) = (y(G,kb))* 1398 b_matrix(kb, jb) = b_matrix(kb, jb) + REAL(res_matrix(ig, kb)*res_matrix(ig, jb), dp) 1399 END DO 1400 END DO 1401 END DO 1402 1403 CALL mp_sum(b_matrix, para_env%group) 1404 DO jb = 1, nb 1405 b_matrix(jb, jb) = b_matrix(jb, jb) + reg_par 1406 END DO 1407 ! invert B 1408 CALL invert_matrix(b_matrix, binv_matrix, inv_err) 1409 1410 ! A = Binv Ynorm^t 1411 a_matrix = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1412 DO kb = 1, nb 1413 res_matrix_global = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1414 DO ig = 1, ng 1415 ig_global = mixing_store%ig_global_index(ig) 1416 res_matrix_global(ig_global) = res_matrix(ig, kb) 1417 END DO 1418 CALL mp_sum(res_matrix_global, para_env%group) 1419 1420 DO jb = 1, nb 1421 DO ig = 1, ng 1422 ig_global = mixing_store%ig_global_index(ig) 1423 a_matrix(jb, ig_global) = a_matrix(jb, ig_global) + & 1424 binv_matrix(jb, kb)*res_matrix_global(ig_global) 1425 END DO 1426 END DO 1427 END DO 1428 CALL mp_sum(a_matrix, para_env%group) 1429 1430 ! compute the two components of gn that will be used to update rho 1431 gn => mixing_store%res_buffer(ib, ispin)%cc 1432 pgn_norm = 0.0_dp 1433 1434 IF (use_zgemm) THEN 1435 1436 CALL zgemm("N", "N", ng, ng_global, nb, cmone, step_matrix(1, 1), ng, & 1437 a_matrix(1, 1), nb, czero, sa(1, 1), ng) 1438 CALL zgemm("N", "N", ng, ng_global, nb, cmone, res_matrix(1, 1), ng, & 1439 a_matrix(1, 1), nb, czero, ya(1, 1), ng) 1440 DO ig = 1, ng 1441 ig_global = mixing_store%ig_global_index(ig) 1442 ya(ig, ig_global) = ya(ig, ig_global) + CMPLX(1.0_dp, 0.0_dp, KIND=dp) 1443 END DO 1444 1445 CALL zgemv("N", ng, ng_global, cone, sa(1, 1), & 1446 ng, gn_global(1), 1, czero, pgn(1), 1) 1447 CALL zgemv("N", ng, ng_global, cone, ya(1, 1), & 1448 ng, gn_global(1), 1, czero, ugn(1), 1) 1449 1450 DO ig = 1, ng 1451 pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + & 1452 AIMAG(pgn(ig))*AIMAG(pgn(ig)) 1453 END DO 1454 CALL mp_sum(pgn_norm, para_env%group) 1455 ELSEIF (use_zgemm_rev) THEN 1456 1457 CALL zgemv("N", nb, ng_global, cone, a_matrix(1, 1), & 1458 nb, gn_global(1), 1, czero, tmp_vec(1), 1) 1459 1460 CALL zgemv("N", ng, nb, cmone, step_matrix(1, 1), ng, & 1461 tmp_vec(1), 1, czero, pgn(1), 1) 1462 1463 CALL zgemv("N", ng, nb, cmone, res_matrix(1, 1), ng, & 1464 tmp_vec(1), 1, czero, ugn(1), 1) 1465 1466 DO ig = 1, ng 1467 pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + & 1468 AIMAG(pgn(ig))*AIMAG(pgn(ig)) 1469 ugn(ig) = ugn(ig) + gn(ig) 1470 END DO 1471 CALL mp_sum(pgn_norm, para_env%group) 1472 1473 ELSE 1474 DO ig = 1, ng 1475 pgn(ig) = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1476 ugn(ig) = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1477 ig_global = mixing_store%ig_global_index(ig) 1478 DO iig = 1, ng_global 1479 saa = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1480 yaa = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1481 1482 IF (ig_global == iig) yaa = CMPLX(1.0_dp, 0.0_dp, KIND=dp) 1483 1484 DO jb = 1, nb 1485 saa = saa - step_matrix(ig, jb)*a_matrix(jb, iig) 1486 yaa = yaa - res_matrix(ig, jb)*a_matrix(jb, iig) 1487 END DO 1488 pgn(ig) = pgn(ig) + saa*gn_global(iig) 1489 ugn(ig) = ugn(ig) + yaa*gn_global(iig) 1490 END DO 1491 END DO 1492 DO ig = 1, ng 1493 pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + & 1494 AIMAG(pgn(ig))*AIMAG(pgn(ig)) 1495 END DO 1496 CALL mp_sum(pgn_norm, para_env%group) 1497 END IF 1498 1499 gn_norm = mixing_store%norm_res_buffer(ib, ispin) 1500 gn_norm_old = mixing_store%norm_res_buffer(ib_prev, ispin) 1501 IF (ib_prev /= 0) THEN 1502 sigma_tilde = sigma_old*MAX(0.5_dp, MIN(2.0_dp, gn_norm_old/gn_norm)) 1503 ELSE 1504 sigma_tilde = 0.5_dp 1505 END IF 1506 sigma_tilde = 0.1_dp 1507 ! Step size for the unpredicted component 1508 step_size = MIN(sigma_tilde, r_step*pgn_norm/gn_norm, sigma_max) 1509 sigma_old = step_size 1510 1511 ! update the density 1512 DO ig = 1, ng 1513 prec = mixing_store%kerker_factor(ig) 1514 rho_g(ispin)%pw%cc(ig) = mixing_store%rhoin_buffer(ib, ispin)%cc(ig) & 1515 - prec*step_size*ugn(ig) + prec*pgn(ig) ! - 0.1_dp * prec* gn(ig) 1516 mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig) 1517 END DO 1518 1519 END DO ! ispin 1520 1521 ! Deallocate temporary arrays 1522 DEALLOCATE (step_matrix, res_matrix) 1523 DEALLOCATE (norm_res) 1524 DEALLOCATE (norm_res_up) 1525 DEALLOCATE (norm_res_low) 1526 DEALLOCATE (a_matrix, b_matrix, binv_matrix) 1527 DEALLOCATE (ugn, pgn) 1528 IF (use_zgemm) THEN 1529 DEALLOCATE (sa, ya) 1530 END IF 1531 IF (use_zgemm_rev) THEN 1532 DEALLOCATE (tmp_vec) 1533 END IF 1534 DEALLOCATE (gn_global, res_matrix_global) 1535 1536 CALL timestop(handle) 1537 1538 END SUBROUTINE multisecant_mixing 1539 1540END MODULE qs_gspace_mixing 1541