1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief localize wavefunctions 8!> linear response scf 9!> \par History 10!> created 07-2005 [MI] 11!> \author MI 12! ************************************************************************************************** 13MODULE qs_linres_methods 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind 16 USE cp_control_types, ONLY: dft_control_type 17 USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,& 18 cp_dbcsr_sm_fm_multiply,& 19 dbcsr_allocate_matrix_set 20 USE cp_external_control, ONLY: external_control 21 USE cp_files, ONLY: close_file,& 22 open_file 23 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,& 24 cp_fm_trace 25 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 26 cp_fm_struct_release,& 27 cp_fm_struct_type 28 USE cp_fm_types, ONLY: cp_fm_create,& 29 cp_fm_get_info,& 30 cp_fm_get_submatrix,& 31 cp_fm_p_type,& 32 cp_fm_release,& 33 cp_fm_set_submatrix,& 34 cp_fm_to_fm,& 35 cp_fm_type 36 USE cp_fm_vect, ONLY: cp_fm_vect_dealloc 37 USE cp_gemm_interface, ONLY: cp_gemm 38 USE cp_log_handling, ONLY: cp_get_default_logger,& 39 cp_logger_type,& 40 cp_to_string 41 USE cp_output_handling, ONLY: cp_p_file,& 42 cp_print_key_finished_output,& 43 cp_print_key_generate_filename,& 44 cp_print_key_should_output,& 45 cp_print_key_unit_nr 46 USE cp_para_types, ONLY: cp_para_env_type 47 USE dbcsr_api, ONLY: dbcsr_checksum,& 48 dbcsr_copy,& 49 dbcsr_filter,& 50 dbcsr_p_type,& 51 dbcsr_set,& 52 dbcsr_type 53 USE hartree_local_methods, ONLY: Vh_1c_gg_integrals 54 USE input_constants, ONLY: do_loc_none,& 55 op_loc_berry,& 56 ot_precond_none,& 57 ot_precond_solver_default,& 58 state_loc_all 59 USE input_section_types, ONLY: section_get_ival,& 60 section_get_lval,& 61 section_get_rval,& 62 section_vals_get_subs_vals,& 63 section_vals_type,& 64 section_vals_val_get 65 USE kinds, ONLY: default_path_length,& 66 default_string_length,& 67 dp 68 USE lri_environment_types, ONLY: lri_density_type,& 69 lri_environment_type,& 70 lri_kind_type 71 USE lri_ks_methods, ONLY: calculate_lri_ks_matrix 72 USE machine, ONLY: m_flush,& 73 m_walltime 74 USE message_passing, ONLY: mp_bcast,& 75 mp_sum 76 USE mulliken, ONLY: ao_charges 77 USE particle_types, ONLY: particle_type 78 USE preconditioner, ONLY: apply_preconditioner,& 79 make_preconditioner 80 USE pw_env_types, ONLY: pw_env_get,& 81 pw_env_type 82 USE pw_methods, ONLY: pw_axpy,& 83 pw_copy,& 84 pw_transfer,& 85 pw_zero 86 USE pw_poisson_methods, ONLY: pw_poisson_solve 87 USE pw_poisson_types, ONLY: pw_poisson_type 88 USE pw_pool_types, ONLY: pw_pool_create_pw,& 89 pw_pool_give_back_pw,& 90 pw_pool_p_type,& 91 pw_pool_type 92 USE pw_types, ONLY: COMPLEXDATA1D,& 93 REALDATA3D,& 94 REALSPACE,& 95 RECIPROCALSPACE,& 96 pw_create,& 97 pw_p_type,& 98 pw_release,& 99 pw_retain 100 USE qs_environment_types, ONLY: get_qs_env,& 101 qs_environment_type 102 USE qs_gapw_densities, ONLY: prepare_gapw_den 103 USE qs_integrate_potential, ONLY: integrate_v_rspace,& 104 integrate_v_rspace_diagonal,& 105 integrate_v_rspace_one_center 106 USE qs_kind_types, ONLY: get_qs_kind,& 107 get_qs_kind_set,& 108 qs_kind_type 109 USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type 110 USE qs_ks_atom, ONLY: update_ks_atom 111 USE qs_linres_types, ONLY: linres_control_type 112 USE qs_loc_methods, ONLY: qs_loc_driver 113 USE qs_loc_types, ONLY: get_qs_loc_env,& 114 localized_wfn_control_type,& 115 qs_loc_env_create,& 116 qs_loc_env_new_type,& 117 qs_loc_env_release,& 118 qs_loc_env_retain 119 USE qs_loc_utils, ONLY: loc_write_restart,& 120 qs_loc_control_init,& 121 qs_loc_init 122 USE qs_mo_types, ONLY: get_mo_set,& 123 mo_set_p_type 124 USE qs_p_env_types, ONLY: qs_p_env_type 125 USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace 126 USE qs_rho_methods, ONLY: qs_rho_rebuild,& 127 qs_rho_update_rho 128 USE qs_rho_types, ONLY: qs_rho_get,& 129 qs_rho_type 130 USE qs_vxc_atom, ONLY: calculate_xc_2nd_deriv_atom 131 USE string_utilities, ONLY: xstring 132 USE xc, ONLY: xc_calc_2nd_deriv,& 133 xc_prep_2nd_deriv,& 134 xc_vxc_pw_create 135 USE xc_derivatives, ONLY: xc_functionals_get_needs 136 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 137 USE xc_rho_set_types, ONLY: xc_rho_set_create,& 138 xc_rho_set_get,& 139 xc_rho_set_release,& 140 xc_rho_set_type,& 141 xc_rho_set_update 142 USE xtb_ehess, ONLY: xtb_coulomb_hessian 143 USE xtb_types, ONLY: get_xtb_atom_param,& 144 xtb_atom_type 145#include "./base/base_uses.f90" 146 147 IMPLICIT NONE 148 149 PRIVATE 150 151 ! *** Public subroutines *** 152 PUBLIC :: linres_localize, linres_solver 153 PUBLIC :: linres_write_restart, linres_read_restart 154 155 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_methods' 156 157! ************************************************************************************************** 158 159CONTAINS 160 161! ************************************************************************************************** 162!> \brief Find the centers and spreads of the wfn, 163!> if required apply a localization algorithm 164!> \param qs_env ... 165!> \param linres_control ... 166!> \param nspins ... 167!> \param centers_only ... 168!> \par History 169!> 07.2005 created [MI] 170!> \author MI 171! ************************************************************************************************** 172 SUBROUTINE linres_localize(qs_env, linres_control, nspins, centers_only) 173 174 TYPE(qs_environment_type), POINTER :: qs_env 175 TYPE(linres_control_type), POINTER :: linres_control 176 INTEGER, INTENT(IN) :: nspins 177 LOGICAL, INTENT(IN), OPTIONAL :: centers_only 178 179 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_localize', & 180 routineP = moduleN//':'//routineN 181 182 INTEGER :: iounit, ispin, istate, nmoloc(2) 183 LOGICAL :: my_centers_only 184 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_localized 185 TYPE(cp_fm_type), POINTER :: mo_coeff 186 TYPE(cp_logger_type), POINTER :: logger 187 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control 188 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 189 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 190 TYPE(section_vals_type), POINTER :: loc_print_section, loc_section, & 191 lr_section 192 193 NULLIFY (logger, lr_section, loc_section, loc_print_section, localized_wfn_control) 194 logger => cp_get_default_logger() 195 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") 196 loc_section => section_vals_get_subs_vals(lr_section, "LOCALIZE") 197 loc_print_section => section_vals_get_subs_vals(lr_section, "LOCALIZE%PRINT") 198 iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 199 extension=".linresLog") 200 my_centers_only = .FALSE. 201 IF (PRESENT(centers_only)) my_centers_only = centers_only 202 203 NULLIFY (mos, mo_coeff, qs_loc_env, mos_localized) 204 CALL get_qs_env(qs_env=qs_env, mos=mos) 205 CALL qs_loc_env_create(qs_loc_env) 206 CALL qs_loc_env_retain(qs_loc_env) 207 linres_control%qs_loc_env => qs_loc_env 208 CALL qs_loc_env_release(qs_loc_env) 209 qs_loc_env => linres_control%qs_loc_env 210 CALL qs_loc_control_init(qs_loc_env, loc_section, do_homo=.TRUE.) 211 CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control) 212 213 ALLOCATE (mos_localized(nspins)) 214 DO ispin = 1, nspins 215 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 216 CALL cp_fm_create(mos_localized(ispin)%matrix, mo_coeff%matrix_struct) 217 CALL cp_fm_to_fm(mo_coeff, mos_localized(ispin)%matrix) 218 END DO 219 220 nmoloc(1:2) = 0 221 IF (my_centers_only) THEN 222 localized_wfn_control%set_of_states = state_loc_all 223 localized_wfn_control%localization_method = do_loc_none 224 localized_wfn_control%operator_type = op_loc_berry 225 ENDIF 226 227 CALL qs_loc_init(qs_env, qs_loc_env, loc_section, mos_localized=mos_localized, & 228 do_homo=.TRUE.) 229 230 ! The orbital centers are stored in linres_control%localized_wfn_control 231 DO ispin = 1, nspins 232 CALL qs_loc_driver(qs_env, qs_loc_env, loc_print_section, myspin=ispin, & 233 ext_mo_coeff=mos_localized(ispin)%matrix) 234 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 235 CALL cp_fm_to_fm(mos_localized(ispin)%matrix, mo_coeff) 236 END DO 237 238 CALL loc_write_restart(qs_loc_env, loc_print_section, mos, & 239 mos_localized, do_homo=.TRUE.) 240 CALL cp_fm_vect_dealloc(mos_localized) 241 242 ! Write Centers and Spreads on std out 243 IF (iounit > 0) THEN 244 DO ispin = 1, nspins 245 WRITE (iounit, "(/,T2,A,I2)") & 246 "WANNIER CENTERS for spin ", ispin 247 WRITE (iounit, "(/,T18,A,3X,A)") & 248 "--------------- Centers --------------- ", & 249 "--- Spreads ---" 250 DO istate = 1, SIZE(localized_wfn_control%centers_set(ispin)%array, 2) 251 WRITE (iounit, "(T5,A6,I6,2X,3f12.6,5X,f12.6)") & 252 'state ', istate, localized_wfn_control%centers_set(ispin)%array(1:3, istate), & 253 localized_wfn_control%centers_set(ispin)%array(4, istate) 254 END DO 255 END DO 256 CALL m_flush(iounit) 257 END IF 258 259 END SUBROUTINE linres_localize 260 261! ************************************************************************************************** 262!> \brief scf loop to optimize the first order wavefunctions (psi1) 263!> given a perturbation as an operator applied to the ground 264!> state orbitals (h1_psi0) 265!> \param p_env ... 266!> \param qs_env ... 267!> \param psi1 ... 268!> \param h1_psi0 ... 269!> \param psi0_order ... 270!> \param iounit ... 271!> \param should_stop ... 272!> \par History 273!> 07.2005 created [MI] 274!> \author MI 275! ************************************************************************************************** 276 SUBROUTINE linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, should_stop) 277 TYPE(qs_p_env_type), POINTER :: p_env 278 TYPE(qs_environment_type), POINTER :: qs_env 279 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: psi1, h1_psi0, psi0_order 280 INTEGER, INTENT(IN) :: iounit 281 LOGICAL, INTENT(OUT) :: should_stop 282 283 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_solver', routineP = moduleN//':'//routineN 284 285 INTEGER :: handle, ispin, iter, maxnmo, maxnmo_o, & 286 nao, ncol, nmo, nspins 287 LOGICAL :: restart 288 REAL(dp) :: norm_res, t1, t2 289 REAL(dp), DIMENSION(:), POINTER :: alpha, beta, tr_pAp, tr_rz0, tr_rz00, & 290 tr_rz1 291 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: Ap, chc, mo_coeff_array, p, r, Sc, z 292 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 293 TYPE(cp_fm_type), POINTER :: buf, mo_coeff 294 TYPE(cp_para_env_type), POINTER :: para_env 295 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_t 296 TYPE(dbcsr_type), POINTER :: matrix_x 297 TYPE(dft_control_type), POINTER :: dft_control 298 TYPE(linres_control_type), POINTER :: linres_control 299 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 300 301 CALL timeset(routineN, handle) 302 303 NULLIFY (dft_control, linres_control, matrix_s, matrix_t, matrix_ks, para_env) 304 NULLIFY (Ap, r, p, z, buf, mos, tmp_fm_struct, mo_coeff) 305 NULLIFY (Sc, chc) 306 307 t1 = m_walltime() 308 309 CALL get_qs_env(qs_env=qs_env, & 310 matrix_ks=matrix_ks, & 311 matrix_s=matrix_s, & 312 kinetic=matrix_t, & 313 dft_control=dft_control, & 314 linres_control=linres_control, & 315 para_env=para_env, & 316 mos=mos) 317 318 nspins = dft_control%nspins 319 CALL get_mo_set(mos(1)%mo_set, nao=nao) 320 maxnmo = 0 321 maxnmo_o = 0 322 DO ispin = 1, nspins 323 CALL get_mo_set(mos(ispin)%mo_set, nmo=ncol) 324 maxnmo = MAX(maxnmo, ncol) 325 CALL cp_fm_get_info(psi0_order(ispin)%matrix, ncol_global=ncol) 326 maxnmo_o = MAX(maxnmo_o, ncol) 327 ENDDO 328 ! 329 CALL check_p_env_init(p_env, linres_control, nspins) 330 ! 331 ! allocate the vectors 332 ALLOCATE (alpha(nspins), beta(nspins), tr_pAp(nspins), tr_rz0(nspins), tr_rz00(nspins), tr_rz1(nspins), & 333 r(nspins), p(nspins), z(nspins), Ap(nspins), mo_coeff_array(nspins)) 334 DO ispin = 1, nspins 335 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff) 336 mo_coeff_array(ispin)%matrix => mo_coeff 337 ENDDO 338 ! 339 DO ispin = 1, nspins 340 NULLIFY (r(ispin)%matrix, p(ispin)%matrix, z(ispin)%matrix, Ap(ispin)%matrix) 341 CALL cp_fm_create(r(ispin)%matrix, psi1(ispin)%matrix%matrix_struct) 342 CALL cp_fm_create(p(ispin)%matrix, psi1(ispin)%matrix%matrix_struct) 343 CALL cp_fm_create(z(ispin)%matrix, psi1(ispin)%matrix%matrix_struct) 344 CALL cp_fm_create(Ap(ispin)%matrix, psi1(ispin)%matrix%matrix_struct) 345 ENDDO 346 ! 347 ! compute S*C0, C0_order'*H*C0_order (this should be done once for all) 348 ALLOCATE (chc(nspins), Sc(nspins)) 349 DO ispin = 1, nspins 350 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff, nmo=nmo) 351 CALL cp_fm_create(Sc(ispin)%matrix, mo_coeff%matrix_struct) 352 NULLIFY (tmp_fm_struct) 353 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, & 354 ncol_global=nmo, para_env=para_env, & 355 context=mo_coeff%matrix_struct%context) 356 CALL cp_fm_create(chc(ispin)%matrix, tmp_fm_struct) 357 CALL cp_fm_struct_release(tmp_fm_struct) 358 ENDDO 359 ! 360 DO ispin = 1, nspins 361 ! 362 ! C0_order' * H * C0_order 363 mo_coeff => psi0_order(ispin)%matrix 364 CALL cp_fm_create(buf, mo_coeff%matrix_struct) 365 CALL cp_fm_get_info(mo_coeff, ncol_global=ncol) 366 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, mo_coeff, buf, ncol) 367 CALL cp_gemm('T', 'N', ncol, ncol, nao, -1.0_dp, mo_coeff, buf, 0.0_dp, chc(ispin)%matrix) 368 CALL cp_fm_release(buf) 369 ! 370 ! S * C0 371 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff) 372 CALL cp_fm_get_info(mo_coeff, ncol_global=ncol) 373 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, Sc(ispin)%matrix, ncol) 374 ENDDO 375 ! 376 ! 377 ! 378 ! header 379 IF (iounit > 0) THEN 380 WRITE (iounit, "(/,T3,A,T16,A,T25,A,T38,A,T52,A,T72,A,/,T3,A)") & 381 "Iteration", "Method", "Restart", "Stepsize", "Convergence", "Time", & 382 REPEAT("-", 80) 383 ENDIF 384 ! 385 ! orthogonalize x with respect to the psi0 386 CALL preortho(psi1, mo_coeff_array, Sc) 387 ! 388 ! build the preconditioner 389 IF (linres_control%preconditioner_type /= ot_precond_none) THEN 390 IF (p_env%new_preconditioner) THEN 391 p_env%os_valid = .FALSE. 392 DO ispin = 1, nspins 393 IF (ASSOCIATED(matrix_t)) THEN 394 CALL make_preconditioner(p_env%preconditioner(ispin), & 395 linres_control%preconditioner_type, ot_precond_solver_default, & 396 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, matrix_t(1)%matrix, & 397 mos(ispin)%mo_set, linres_control%energy_gap) 398 ELSE 399 NULLIFY (matrix_x) 400 CALL make_preconditioner(p_env%preconditioner(ispin), & 401 linres_control%preconditioner_type, ot_precond_solver_default, & 402 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, matrix_x, & 403 mos(ispin)%mo_set, linres_control%energy_gap) 404 END IF 405 ENDDO 406 p_env%new_preconditioner = .FALSE. 407 ENDIF 408 ENDIF 409 ! 410 ! initalization of the linear solver 411 ! 412 ! A * x0 413 CALL apply_op(qs_env, p_env, psi0_order, psi1, Ap, chc) 414 ! 415 ! 416 ! r_0 = b - Ax0 417 DO ispin = 1, nspins 418 CALL cp_fm_to_fm(h1_psi0(ispin)%matrix, r(ispin)%matrix) 419 CALL cp_fm_scale_and_add(-1.0_dp, r(ispin)%matrix, -1.0_dp, Ap(ispin)%matrix) 420 ENDDO 421 ! 422 ! proj r 423 CALL postortho(r, mo_coeff_array, Sc) 424 ! 425 ! preconditioner 426 linres_control%flag = "" 427 IF (linres_control%preconditioner_type .EQ. ot_precond_none) THEN 428 ! 429 ! z_0 = r_0 430 DO ispin = 1, nspins 431 CALL cp_fm_to_fm(r(ispin)%matrix, z(ispin)%matrix) 432 ENDDO 433 linres_control%flag = "CG" 434 ELSE 435 ! 436 ! z_0 = M * r_0 437 DO ispin = 1, nspins 438 CALL apply_preconditioner(p_env%preconditioner(ispin), r(ispin)%matrix, & 439 z(ispin)%matrix) 440 ENDDO 441 linres_control%flag = "PCG" 442 ENDIF 443 ! 444 norm_res = 0.0_dp 445 DO ispin = 1, nspins 446 ! 447 ! p_0 = z_0 448 CALL cp_fm_to_fm(z(ispin)%matrix, p(ispin)%matrix) 449 ! 450 ! trace(r_0 * z_0) 451 CALL cp_fm_trace(r(ispin)%matrix, z(ispin)%matrix, tr_rz0(ispin)) 452 IF (tr_rz0(ispin) .LT. 0.0_dp) CPABORT("tr(r_j*z_j) < 0") 453 norm_res = MAX(norm_res, ABS(tr_rz0(ispin))/SQRT(REAL(nao*maxnmo_o, dp))) 454 ENDDO 455 ! 456 alpha(:) = 0.0_dp 457 restart = .FALSE. 458 should_stop = .FALSE. 459 iteration: DO iter = 1, linres_control%max_iter 460 ! 461 ! check convergence 462 linres_control%converged = .FALSE. 463 IF (norm_res .LT. linres_control%eps) THEN 464 linres_control%converged = .TRUE. 465 ENDIF 466 ! 467 t2 = m_walltime() 468 IF (iter .EQ. 1 .OR. MOD(iter, 1) .EQ. 0 .OR. linres_control%converged .OR. restart .OR. should_stop) THEN 469 IF (iounit > 0) THEN 470 WRITE (iounit, "(T5,I5,T18,A3,T28,L1,T38,1E8.2,T48,F16.10,T68,F8.2)") & 471 iter, linres_control%flag, restart, MAXVAL(alpha), norm_res, t2 - t1 472 CALL m_flush(iounit) 473 ENDIF 474 ENDIF 475 ! 476 IF (linres_control%converged) THEN 477 IF (iounit > 0) THEN 478 WRITE (iounit, "(/,T2,A,I4,A,/)") "The linear solver converged in ", iter, " iterations." 479 CALL m_flush(iounit) 480 ENDIF 481 EXIT iteration 482 ELSE IF (should_stop) THEN 483 IF (iounit > 0) THEN 484 WRITE (iounit, "(/,T2,A,I4,A,/)") "The linear solver did NOT converge! External stop" 485 CALL m_flush(iounit) 486 END IF 487 EXIT iteration 488 ENDIF 489 ! 490 ! Max number of iteration reached 491 IF (iter == linres_control%max_iter) THEN 492 IF (iounit > 0) THEN 493 WRITE (iounit, "(/,T2,A/)") & 494 "The linear solver didnt converge! Maximum number of iterations reached." 495 CALL m_flush(iounit) 496 ENDIF 497 linres_control%converged = .FALSE. 498 ENDIF 499 ! 500 ! Apply the operators that do not depend on the perturbation 501 CALL apply_op(qs_env, p_env, psi0_order, p, Ap, chc) 502 ! 503 ! proj Ap onto the virtual subspace 504 CALL postortho(Ap, mo_coeff_array, Sc) 505 ! 506 DO ispin = 1, nspins 507 ! 508 ! tr(Ap_j*p_j) 509 CALL cp_fm_trace(Ap(ispin)%matrix, p(ispin)%matrix, tr_pAp(ispin)) 510 IF (tr_pAp(ispin) .LT. 0.0_dp) THEN 511 ! try to fix it by getting rid of the preconditioner 512 IF (iter > 1) THEN 513 CALL cp_fm_scale_and_add(beta(ispin), p(ispin)%matrix, -1.0_dp, z(ispin)%matrix) 514 CALL cp_fm_trace(r(ispin)%matrix, r(ispin)%matrix, tr_rz1(ispin)) 515 beta(ispin) = tr_rz1(ispin)/tr_rz00(ispin) 516 CALL cp_fm_scale_and_add(beta(ispin), p(ispin)%matrix, 1.0_dp, r(ispin)%matrix) 517 tr_rz0(ispin) = tr_rz1(ispin) 518 ELSE 519 CALL cp_fm_to_fm(r(ispin)%matrix, p(ispin)%matrix) 520 CALL cp_fm_trace(r(ispin)%matrix, r(ispin)%matrix, tr_rz0(ispin)) 521 END IF 522 linres_control%flag = "CG" 523 524 CALL apply_op(qs_env, p_env, psi0_order, p, Ap, chc) 525 CALL postortho(Ap, mo_coeff_array, Sc) 526 CALL cp_fm_trace(Ap(ispin)%matrix, p(ispin)%matrix, tr_pAp(ispin)) 527 CPABORT("tr(Ap_j*p_j) < 0") 528 END IF 529 ! 530 ! alpha = tr(r_j*z_j) / tr(Ap_j*p_j) 531 IF (tr_pAp(ispin) .LT. 1.0e-10_dp) THEN 532 alpha(ispin) = 1.0_dp 533 ELSE 534 alpha(ispin) = tr_rz0(ispin)/tr_pAp(ispin) 535 ENDIF 536 ! 537 ! x_j+1 = x_j + alpha * p_j 538 CALL cp_fm_scale_and_add(1.0_dp, psi1(ispin)%matrix, alpha(ispin), p(ispin)%matrix) 539 ENDDO 540 ! 541 ! need to recompute the residue 542 restart = .FALSE. 543 IF (MOD(iter, linres_control%restart_every) .EQ. 0) THEN 544 ! 545 ! r_j+1 = b - A * x_j+1 546 CALL apply_op(qs_env, p_env, psi0_order, psi1, Ap, chc) 547 ! 548 DO ispin = 1, nspins 549 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff) 550 CALL cp_fm_to_fm(h1_psi0(ispin)%matrix, r(ispin)%matrix) 551 CALL cp_fm_scale_and_add(-1.0_dp, r(ispin)%matrix, -1.0_dp, Ap(ispin)%matrix) 552 ENDDO 553 CALL postortho(r, mo_coeff_array, Sc) 554 ! 555 restart = .TRUE. 556 ELSE 557 ! proj Ap onto the virtual subspace 558 CALL postortho(Ap, mo_coeff_array, Sc) 559 ! 560 ! r_j+1 = r_j - alpha * Ap_j 561 DO ispin = 1, nspins 562 CALL cp_fm_scale_and_add(1.0_dp, r(ispin)%matrix, -alpha(ispin), Ap(ispin)%matrix) 563 ENDDO 564 restart = .FALSE. 565 ENDIF 566 ! 567 ! preconditioner 568 linres_control%flag = "" 569 IF (linres_control%preconditioner_type .EQ. ot_precond_none) THEN 570 ! 571 ! z_j+1 = r_j+1 572 DO ispin = 1, nspins 573 CALL cp_fm_to_fm(r(ispin)%matrix, z(ispin)%matrix) 574 ENDDO 575 linres_control%flag = "CG" 576 ELSE 577 ! 578 ! z_j+1 = M * r_j+1 579 DO ispin = 1, nspins 580 CALL apply_preconditioner(p_env%preconditioner(ispin), r(ispin)%matrix, & 581 z(ispin)%matrix) 582 ENDDO 583 linres_control%flag = "PCG" 584 ENDIF 585 ! 586 norm_res = 0.0_dp 587 DO ispin = 1, nspins 588 ! 589 ! tr(r_j+1*z_j+1) 590 CALL cp_fm_trace(r(ispin)%matrix, z(ispin)%matrix, tr_rz1(ispin)) 591 IF (tr_rz1(ispin) .LT. 0.0_dp) CPABORT("tr(r_j+1*z_j+1) < 0") 592 norm_res = MAX(norm_res, tr_rz1(ispin)/SQRT(REAL(nao*maxnmo_o, dp))) 593 ! 594 ! beta = tr(r_j+1*z_j+1) / tr(r_j*z_j) 595 IF (tr_rz0(ispin) .LT. 1.0e-10_dp) THEN 596 beta(ispin) = 0.0_dp 597 ELSE 598 beta(ispin) = tr_rz1(ispin)/tr_rz0(ispin) 599 ENDIF 600 ! 601 ! p_j+1 = z_j+1 + beta * p_j 602 CALL cp_fm_scale_and_add(beta(ispin), p(ispin)%matrix, 1.0_dp, z(ispin)%matrix) 603 tr_rz00(ispin) = tr_rz0(ispin) 604 tr_rz0(ispin) = tr_rz1(ispin) 605 ENDDO 606 607 ! Can we exit the SCF loop? 608 CALL external_control(should_stop, "LINRES", target_time=qs_env%target_time, & 609 start_time=qs_env%start_time) 610 611 ENDDO iteration 612 ! 613 ! proj psi1 614 CALL preortho(psi1, mo_coeff_array, Sc) 615 ! 616 ! clean up 617 DO ispin = 1, nspins 618 CALL cp_fm_release(r(ispin)%matrix) 619 CALL cp_fm_release(p(ispin)%matrix) 620 CALL cp_fm_release(z(ispin)%matrix) 621 CALL cp_fm_release(Ap(ispin)%matrix) 622 ! 623 CALL cp_fm_release(Sc(ispin)%matrix) 624 CALL cp_fm_release(chc(ispin)%matrix) 625 ENDDO 626 DEALLOCATE (alpha, beta, tr_pAp, tr_rz0, tr_rz00, tr_rz1, r, p, z, Ap, Sc, chc, mo_coeff_array) 627 ! 628 CALL timestop(handle) 629 ! 630 END SUBROUTINE linres_solver 631 632! ************************************************************************************************** 633!> \brief ... 634!> \param qs_env ... 635!> \param p_env ... 636!> \param c0 ... 637!> \param v ... 638!> \param Av ... 639!> \param chc ... 640! ************************************************************************************************** 641 SUBROUTINE apply_op(qs_env, p_env, c0, v, Av, chc) 642 ! 643 TYPE(qs_environment_type), POINTER :: qs_env 644 TYPE(qs_p_env_type), POINTER :: p_env 645 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av, chc 646 647 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op', routineP = moduleN//':'//routineN 648 649 INTEGER :: handle, ispin, nspins 650 REAL(dp) :: chksum 651 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho1_ao 652 TYPE(dft_control_type), POINTER :: dft_control 653 TYPE(linres_control_type), POINTER :: linres_control 654 TYPE(qs_rho_type), POINTER :: rho 655 656 CALL timeset(routineN, handle) 657 658 CPASSERT(ASSOCIATED(v)) 659 CPASSERT(ASSOCIATED(Av)) 660 CPASSERT(ASSOCIATED(chc)) 661 662 NULLIFY (dft_control, matrix_ks, matrix_s, linres_control, rho1_ao) 663 CALL get_qs_env(qs_env=qs_env, & 664 matrix_ks=matrix_ks, & 665 matrix_s=matrix_s, & 666 dft_control=dft_control, & 667 linres_control=linres_control) 668 669 nspins = dft_control%nspins 670 671 ! apply the uncoupled operator 672 DO ispin = 1, nspins 673 CALL apply_op_1(v(ispin)%matrix, Av(ispin)%matrix, matrix_ks(ispin)%matrix, & 674 matrix_s(1)%matrix, chc(ispin)%matrix) 675 ENDDO 676 677 IF (linres_control%do_kernel) THEN 678 679 ! build DM, refill p1, build_dm_response keeps sparse structure 680 DO ispin = 1, nspins 681 CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix) 682 ENDDO 683 CALL build_dm_response(c0, v, p_env%p1) 684 DO ispin = 1, nspins 685 CALL dbcsr_filter(p_env%p1(ispin)%matrix, linres_control%eps_filter) 686 ENDDO 687 688 chksum = 0.0_dp 689 DO ispin = 1, nspins 690 chksum = chksum + dbcsr_checksum(p_env%p1(ispin)%matrix) 691 ENDDO 692 693 ! skip the kernel if the DM is very small 694 IF (chksum .GT. 1.0E-14_dp) THEN 695 696 CALL p_env_check_i_alloc(p_env, qs_env) 697 698 CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao) 699 DO ispin = 1, nspins 700 CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1(ispin)%matrix) 701 ENDDO 702 703 CALL qs_rho_update_rho(rho_struct=p_env%rho1, & 704 local_rho_set=p_env%local_rho_set, & 705 qs_env=qs_env) 706 707 CALL get_qs_env(qs_env, rho=rho) ! that could be called before 708 CALL qs_rho_update_rho(rho, qs_env=qs_env) ! that could be called before 709 710 CALL apply_op_2(qs_env, p_env, c0, v, Av, chc) 711 712 ENDIF 713 714 ENDIF 715 716 CALL timestop(handle) 717 718 END SUBROUTINE apply_op 719 720! ************************************************************************************************** 721!> \brief ... 722!> \param c0 ... 723!> \param c1 ... 724!> \param dm ... 725! ************************************************************************************************** 726 SUBROUTINE build_dm_response(c0, c1, dm) 727 ! 728 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, c1 729 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dm 730 731 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_dm_response', & 732 routineP = moduleN//':'//routineN 733 734 INTEGER :: ispin, ncol, nspins 735 736 CPASSERT(ASSOCIATED(c0)) 737 CPASSERT(ASSOCIATED(c1)) 738 CPASSERT(ASSOCIATED(dm)) 739 740 nspins = SIZE(dm, 1) 741 742 DO ispin = 1, nspins 743 CALL dbcsr_set(dm(ispin)%matrix, 0.0_dp) 744 CALL cp_fm_get_info(c0(ispin)%matrix, ncol_global=ncol) 745 CALL cp_dbcsr_plus_fm_fm_t(dm(ispin)%matrix, & 746 matrix_v=c0(ispin)%matrix, & 747 matrix_g=c1(ispin)%matrix, & 748 ncol=ncol, alpha=1.0_dp) 749 CALL cp_dbcsr_plus_fm_fm_t(dm(ispin)%matrix, & 750 matrix_v=c1(ispin)%matrix, & 751 matrix_g=c0(ispin)%matrix, & 752 ncol=ncol, alpha=1.0_dp) 753 ENDDO 754 755 END SUBROUTINE build_dm_response 756 757! ************************************************************************************************** 758!> \brief ... 759!> \param v ... 760!> \param Av ... 761!> \param matrix_ks ... 762!> \param matrix_s ... 763!> \param chc ... 764! ************************************************************************************************** 765 SUBROUTINE apply_op_1(v, Av, matrix_ks, matrix_s, chc) 766 ! 767 TYPE(cp_fm_type), POINTER :: v, Av 768 TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_s 769 TYPE(cp_fm_type), POINTER :: chc 770 771 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op_1', routineP = moduleN//':'//routineN 772 773 INTEGER :: handle, ncol, nrow 774 TYPE(cp_fm_type), POINTER :: buf 775 776 CALL timeset(routineN, handle) 777 ! 778 CPASSERT(ASSOCIATED(v)) 779 CPASSERT(ASSOCIATED(Av)) 780 CPASSERT(ASSOCIATED(matrix_ks)) 781 CPASSERT(ASSOCIATED(matrix_s)) 782 CPASSERT(ASSOCIATED(chc)) 783 NULLIFY (buf) 784 ! 785 CALL cp_fm_create(buf, v%matrix_struct) 786 ! 787 CALL cp_fm_get_info(v, ncol_global=ncol, nrow_global=nrow) 788 ! H * v 789 CALL cp_dbcsr_sm_fm_multiply(matrix_ks, v, Av, ncol) 790 ! v * e (chc already multiplied by -1) 791 CALL cp_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, v, chc, 0.0_dp, buf) 792 ! S * ve 793 CALL cp_dbcsr_sm_fm_multiply(matrix_s, buf, Av, ncol, alpha=1.0_dp, beta=1.0_dp) 794 !Results is H*C1 - S*<iHj>*C1 795 ! 796 CALL cp_fm_release(buf) 797 ! 798 CALL timestop(handle) 799 ! 800 END SUBROUTINE apply_op_1 801 802! ************************************************************************************************** 803!> \brief ... 804!> \param qs_env ... 805!> \param p_env ... 806!> \param c0 ... 807!> \param v ... 808!> \param Av ... 809!> \param chc ... 810! ************************************************************************************************** 811 SUBROUTINE apply_op_2(qs_env, p_env, c0, v, Av, chc) 812 ! 813 TYPE(qs_environment_type), POINTER :: qs_env 814 TYPE(qs_p_env_type), POINTER :: p_env 815 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av, chc 816 817 CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2', routineP = moduleN//':'//routineN 818 819 TYPE(dft_control_type), POINTER :: dft_control 820 821 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 822 IF (dft_control%qs_control%semi_empirical) THEN 823 CPABORT("Linear response not available with SE methods") 824 ELSEIF (dft_control%qs_control%dftb) THEN 825 CPABORT("Linear response not available with DFTB") 826 ELSEIF (dft_control%qs_control%xtb) THEN 827 CALL apply_op_2_xtb(qs_env, p_env, c0, v, Av, chc) 828 ELSE 829 CALL apply_op_2_dft(qs_env, p_env, c0, v, Av, chc) 830 END IF 831 832 END SUBROUTINE apply_op_2 833 834! ************************************************************************************************** 835!> \brief ... 836!> \param qs_env ... 837!> \param p_env ... 838!> \param c0 ... 839!> \param v ... 840!> \param Av ... 841!> \param chc ... 842! ************************************************************************************************** 843 SUBROUTINE apply_op_2_dft(qs_env, p_env, c0, v, Av, chc) 844 TYPE(qs_environment_type), POINTER :: qs_env 845 TYPE(qs_p_env_type), POINTER :: p_env 846 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av, chc 847 848 CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2_dft', routineP = moduleN//':'//routineN 849 REAL(KIND=dp), PARAMETER :: h = 0.001_dp 850 851 INTEGER :: handle, ikind, ispin, ncol, nkind, ns, & 852 nspins 853 INTEGER, DIMENSION(2, 3) :: bo 854 LOGICAL :: deriv2_analytic, gapw, gapw_xc, & 855 lr_triplet, lrigpw, lsd 856 REAL(KIND=dp) :: energy_hartree, energy_hartree_1c, exc, & 857 fac 858 REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc 859 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho3, rhoa, rhob 860 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 861 TYPE(cp_logger_type), POINTER :: logger 862 TYPE(cp_para_env_type), POINTER :: para_env 863 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: k1mat, rho1_ao, rho_ao 864 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, psmat 865 TYPE(dft_control_type), POINTER :: dft_control 866 TYPE(linres_control_type), POINTER :: linres_control 867 TYPE(lri_density_type), POINTER :: lri_density 868 TYPE(lri_environment_type), POINTER :: lri_env 869 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int 870 TYPE(pw_env_type), POINTER :: pw_env 871 TYPE(pw_p_type) :: rho1_tot_gspace, v_hartree_gspace, & 872 v_hartree_rspace 873 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g, rho1_g_pw, rho1_r, rho1_r_pw, rho_g, & 874 rho_r, tau, tau_pw, v_rspace_new, v_xc, vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4 875 TYPE(pw_poisson_type), POINTER :: poisson_env 876 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 877 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 878 TYPE(qs_rho_type), POINTER :: rho, rho1, rho1_xc 879 TYPE(section_vals_type), POINTER :: input, scf_section, xc_fun_section, & 880 xc_section 881 TYPE(xc_rho_cflags_type) :: needs 882 TYPE(xc_rho_set_type), POINTER :: rho1_set 883 884 CALL timeset(routineN, handle) 885 886 NULLIFY (auxbas_pw_pool, pw_pools, pw_env, v_rspace_new, & 887 rho1_r, rho1_g_pw, tau_pw, v_xc, rho1_set, rho1_ao, rho_ao, & 888 poisson_env, input, scf_section, rho, dft_control, logger, rho1_g) 889 logger => cp_get_default_logger() 890 891 energy_hartree = 0.0_dp 892 energy_hartree_1c = 0.0_dp 893 894 CPASSERT(ASSOCIATED(c0)) 895 CPASSERT(ASSOCIATED(v)) 896 CPASSERT(ASSOCIATED(Av)) 897 CPASSERT(ASSOCIATED(chc)) 898 899 CPASSERT(ASSOCIATED(p_env%kpp1_env)) 900 CPASSERT(ASSOCIATED(p_env%kpp1)) 901 rho1 => p_env%rho1 902 rho1_xc => p_env%rho1_xc 903 CPASSERT(ASSOCIATED(rho1)) 904 905 CPASSERT(p_env%kpp1_env%ref_count > 0) 906 907 CALL get_qs_env(qs_env=qs_env, & 908 pw_env=pw_env, & 909 input=input, & 910 rho=rho, & 911 linres_control=linres_control, & 912 dft_control=dft_control) 913 914 lrigpw = dft_control%qs_control%lrigpw 915 IF (lrigpw) THEN 916 CALL get_qs_env(qs_env, & 917 lri_env=lri_env, & 918 lri_density=lri_density, & 919 atomic_kind_set=atomic_kind_set) 920 ENDIF 921 922 CALL qs_rho_get(rho, rho_ao=rho_ao) 923 924 lr_triplet = linres_control%lr_triplet 925 CALL kpp1_check_i_alloc(p_env%kpp1_env, qs_env, lr_triplet) 926 gapw = dft_control%qs_control%gapw 927 gapw_xc = dft_control%qs_control%gapw_xc 928 IF (gapw_xc) THEN 929 CPASSERT(ASSOCIATED(rho1_xc)) 930 END IF 931 932 nspins = SIZE(p_env%kpp1) 933 lsd = (nspins == 2) 934 935 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 936 scf_section => section_vals_get_subs_vals(input, "DFT%SCF") 937 938 p_env%kpp1_env%iter = p_env%kpp1_env%iter + 1 939 940 ! gets the tmp grids 941 CPASSERT(ASSOCIATED(pw_env)) 942 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 943 pw_pools=pw_pools, poisson_env=poisson_env) 944 ALLOCATE (v_rspace_new(nspins)) 945 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 946 use_data=COMPLEXDATA1D, & 947 in_space=RECIPROCALSPACE) 948 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_rspace%pw, & 949 use_data=REALDATA3D, & 950 in_space=REALSPACE) 951 952 IF (gapw .OR. gapw_xc) & 953 CALL prepare_gapw_den(qs_env, p_env%local_rho_set, do_rho0=(.NOT. gapw_xc)) 954 955 ! *** calculate the hartree potential on the total density *** 956 CALL pw_pool_create_pw(auxbas_pw_pool, rho1_tot_gspace%pw, & 957 use_data=COMPLEXDATA1D, & 958 in_space=RECIPROCALSPACE) 959 960 CALL qs_rho_get(rho1, rho_g=rho1_g) 961 CALL pw_copy(rho1_g(1)%pw, rho1_tot_gspace%pw) 962 DO ispin = 2, nspins 963 CALL pw_axpy(rho1_g(ispin)%pw, rho1_tot_gspace%pw) 964 END DO 965 IF (gapw) & 966 CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs%pw, rho1_tot_gspace%pw) 967 968 IF (.NOT. (nspins == 1 .AND. lr_triplet)) THEN 969 CALL pw_poisson_solve(poisson_env, rho1_tot_gspace%pw, & 970 energy_hartree, & 971 v_hartree_gspace%pw) 972 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw) 973 ENDIF 974 975 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho1_tot_gspace%pw) 976 977 ! *** calculate the xc potential *** 978 IF (gapw_xc) THEN 979 CALL qs_rho_get(rho1_xc, rho_r=rho1_r) 980 ELSE 981 CALL qs_rho_get(rho1, rho_r=rho1_r) 982 END IF 983 984 IF (nspins == 1 .AND. (lr_triplet)) THEN 985 986 lsd = .TRUE. 987 ALLOCATE (rho1_r_pw(2)) 988 DO ispin = 1, 2 989 NULLIFY (rho1_r_pw(ispin)%pw) 990 CALL pw_create(rho1_r_pw(ispin)%pw, rho1_r(1)%pw%pw_grid, & 991 rho1_r(1)%pw%in_use, rho1_r(1)%pw%in_space) 992 CALL pw_transfer(rho1_r(1)%pw, rho1_r_pw(ispin)%pw) 993 END DO 994 995 ELSE 996 997 ALLOCATE (rho1_r_pw(nspins)) 998 DO ispin = 1, nspins 999 rho1_r_pw(ispin)%pw => rho1_r(ispin)%pw 1000 CALL pw_retain(rho1_r_pw(ispin)%pw) 1001 END DO 1002 1003 END IF 1004 1005 NULLIFY (tau_pw) 1006 1007 deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL") 1008 1009 IF (deriv2_analytic) THEN 1010 !------! 1011 ! rho1 ! 1012 !------! 1013 bo = rho1_r(1)%pw%pw_grid%bounds_local 1014 ! create the place where to store the argument for the functionals 1015 CALL xc_rho_set_create(rho1_set, bo, & 1016 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), & 1017 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), & 1018 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF")) 1019 1020 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL") 1021 needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.) 1022 1023 ! calculate the arguments needed by the functionals 1024 CALL xc_rho_set_update(rho1_set, rho1_r_pw, rho1_g_pw, tau_pw, needs, & 1025 section_get_ival(xc_section, "XC_GRID%XC_DERIV"), & 1026 section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), & 1027 auxbas_pw_pool) 1028 1029 ALLOCATE (v_xc(nspins)) 1030 DO ispin = 1, nspins 1031 NULLIFY (v_xc(ispin)%pw) 1032 CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, & 1033 use_data=REALDATA3D, & 1034 in_space=REALSPACE) 1035 CALL pw_zero(v_xc(ispin)%pw) 1036 END DO 1037 1038 fac = 0._dp 1039 IF (nspins == 1) THEN 1040 IF (lr_triplet) fac = -1.0_dp 1041 END IF 1042 1043 CALL xc_calc_2nd_deriv(v_xc, p_env%kpp1_env%deriv_set, p_env%kpp1_env%rho_set, & 1044 rho1_set, auxbas_pw_pool, xc_section=xc_section, & 1045 tddfpt_fac=fac) 1046 1047 DO ispin = 1, nspins 1048 v_rspace_new(ispin)%pw => v_xc(ispin)%pw 1049 END DO 1050 DEALLOCATE (v_xc) 1051 1052 IF (gapw) CALL calculate_xc_2nd_deriv_atom(p_env, qs_env, xc_section, do_tddft=.FALSE., & 1053 do_triplet=lr_triplet) 1054 1055 CALL xc_rho_set_release(rho1_set) 1056 1057 ELSE 1058 NULLIFY (tau) 1059 ALLOCATE (v_xc(nspins)) 1060 DO ispin = 1, nspins 1061 NULLIFY (v_xc(ispin)%pw) 1062 CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, & 1063 use_data=REALDATA3D, & 1064 in_space=REALSPACE) 1065 CALL pw_zero(v_xc(ispin)%pw) 1066 END DO 1067 ! finite differences 1068 CPASSERT(.NOT. gapw) 1069 IF (nspins == 2) THEN 1070 NULLIFY (vxc_rho_1, vxc_rho_2, rho_g) 1071 ALLOCATE (rho_r(2)) 1072 DO ispin = 1, nspins 1073 NULLIFY (rho_r(ispin)%pw) 1074 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(ispin)%pw, in_space=REALSPACE, use_data=REALDATA3D) 1075 END DO 1076 CALL xc_rho_set_get(p_env%kpp1_env%rho_set, rhoa=rhoa, rhob=rhob) 1077 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) + 0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :) 1078 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) + 0.5_dp*h*rho1_r(2)%pw%cr3d(:, :, :) 1079 CALL xc_vxc_pw_create(vxc_rho_1, tau_pw, exc, rho_r, rho_g, tau, xc_section, & 1080 auxbas_pw_pool, .FALSE., virial_xc) 1081 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) - 0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :) 1082 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) - 0.5_dp*h*rho1_r(2)%pw%cr3d(:, :, :) 1083 CALL xc_vxc_pw_create(vxc_rho_2, tau_pw, exc, rho_r, rho_g, tau, xc_section, & 1084 auxbas_pw_pool, .FALSE., virial_xc) 1085 v_xc(1)%pw%cr3d(:, :, :) = (vxc_rho_1(1)%pw%cr3d(:, :, :) - vxc_rho_2(1)%pw%cr3d(:, :, :))/h 1086 v_xc(2)%pw%cr3d(:, :, :) = (vxc_rho_1(2)%pw%cr3d(:, :, :) - vxc_rho_2(2)%pw%cr3d(:, :, :))/h 1087 DO ispin = 1, nspins 1088 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(ispin)%pw) 1089 END DO 1090 DEALLOCATE (rho_r) 1091 DO ispin = 1, nspins 1092 CALL pw_release(vxc_rho_1(ispin)%pw) 1093 CALL pw_release(vxc_rho_2(ispin)%pw) 1094 END DO 1095 DEALLOCATE (vxc_rho_1, vxc_rho_2) 1096 ELSE IF (nspins == 1 .AND. (lr_triplet)) THEN 1097 NULLIFY (vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4, rho_g) 1098 ALLOCATE (rho_r(2)) 1099 DO ispin = 1, 2 1100 NULLIFY (rho_r(ispin)%pw) 1101 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(ispin)%pw, in_space=REALSPACE, use_data=REALDATA3D) 1102 END DO 1103 CALL xc_rho_set_get(p_env%kpp1_env%rho_set, rhoa=rhoa, rhob=rhob) 1104 ! K(alpha,alpha) 1105 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) + 0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :) 1106 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) 1107 CALL xc_vxc_pw_create(vxc_rho_1, tau_pw, exc, rho_r, rho_g, tau, xc_section, & 1108 auxbas_pw_pool, .FALSE., virial_xc) 1109 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) - 0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :) 1110 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) 1111 CALL xc_vxc_pw_create(vxc_rho_2, tau_pw, exc, rho_r, rho_g, tau, xc_section, & 1112 auxbas_pw_pool, .FALSE., virial_xc) 1113 v_xc(1)%pw%cr3d(:, :, :) = (vxc_rho_1(1)%pw%cr3d(:, :, :) - vxc_rho_2(1)%pw%cr3d(:, :, :))/h 1114 ! K(alpha,beta) 1115 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) 1116 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) + 0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :) 1117 CALL xc_vxc_pw_create(vxc_rho_3, tau_pw, exc, rho_r, rho_g, tau, xc_section, & 1118 auxbas_pw_pool, .FALSE., virial_xc) 1119 rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :) 1120 rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :) - 0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :) 1121 CALL xc_vxc_pw_create(vxc_rho_4, tau_pw, exc, rho_r, rho_g, tau, xc_section, & 1122 auxbas_pw_pool, .FALSE., virial_xc) 1123 v_xc(1)%pw%cr3d(:, :, :) = v_xc(1)%pw%cr3d(:, :, :) - & 1124 (vxc_rho_3(1)%pw%cr3d(:, :, :) - vxc_rho_4(1)%pw%cr3d(:, :, :))/h 1125 DO ispin = 1, 2 1126 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(ispin)%pw) 1127 END DO 1128 DEALLOCATE (rho_r) 1129 DO ispin = 1, 2 1130 CALL pw_release(vxc_rho_1(ispin)%pw) 1131 CALL pw_release(vxc_rho_2(ispin)%pw) 1132 CALL pw_release(vxc_rho_3(ispin)%pw) 1133 CALL pw_release(vxc_rho_4(ispin)%pw) 1134 END DO 1135 DEALLOCATE (vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4) 1136 ELSE 1137 NULLIFY (vxc_rho_1, vxc_rho_2, rho_r, rho_g) 1138 ALLOCATE (rho_r(1)) 1139 NULLIFY (rho_r(1)%pw) 1140 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(1)%pw, in_space=REALSPACE, use_data=REALDATA3D) 1141 CALL xc_rho_set_get(p_env%kpp1_env%rho_set, rho=rho3) 1142 rho_r(1)%pw%cr3d(:, :, :) = rho3(:, :, :) + 0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :) 1143 CALL xc_vxc_pw_create(vxc_rho_1, tau_pw, exc, rho_r, rho_g, tau, xc_section, & 1144 auxbas_pw_pool, .FALSE., virial_xc) 1145 rho_r(1)%pw%cr3d(:, :, :) = rho3(:, :, :) - 0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :) 1146 CALL xc_vxc_pw_create(vxc_rho_2, tau_pw, exc, rho_r, rho_g, tau, xc_section, & 1147 auxbas_pw_pool, .FALSE., virial_xc) 1148 v_xc(1)%pw%cr3d(:, :, :) = (vxc_rho_1(1)%pw%cr3d(:, :, :) - vxc_rho_2(1)%pw%cr3d(:, :, :))/h 1149 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(1)%pw) 1150 DEALLOCATE (rho_r) 1151 CALL pw_release(vxc_rho_1(1)%pw) 1152 CALL pw_release(vxc_rho_2(1)%pw) 1153 DEALLOCATE (vxc_rho_1, vxc_rho_2) 1154 END IF 1155 DO ispin = 1, nspins 1156 v_rspace_new(ispin)%pw => v_xc(ispin)%pw 1157 END DO 1158 DEALLOCATE (v_xc) 1159 END IF 1160 1161 DO ispin = 1, SIZE(rho1_r_pw) 1162 CALL pw_release(rho1_r_pw(ispin)%pw) 1163 END DO 1164 DEALLOCATE (rho1_r_pw) 1165 1166 !-------------------------------! 1167 ! Add both hartree and xc terms ! 1168 !-------------------------------! 1169 DO ispin = 1, nspins 1170 1171 IF (gapw_xc) THEN 1172 ! XC and Hartree are integrated separatedly 1173 ! XC uses the sofft basis set only 1174 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d* & 1175 v_rspace_new(ispin)%pw%pw_grid%dvol 1176 1177 IF (nspins == 1) THEN 1178 1179 IF (.NOT. (lr_triplet)) THEN 1180 1181 v_rspace_new(1)%pw%cr3d = 2.0_dp*v_rspace_new(1)%pw%cr3d 1182 1183 END IF 1184 CALL qs_rho_get(rho1, rho_ao=rho1_ao) 1185 ! remove kpp1_env%v_ao and work directly on k_p_p1 ? 1186 CALL dbcsr_set(p_env%kpp1_env%v_ao(ispin)%matrix, 0.0_dp) 1187 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), & 1188 pmat=rho1_ao(ispin), & 1189 hmat=p_env%kpp1_env%v_ao(ispin), & 1190 qs_env=qs_env, & 1191 calculate_forces=.FALSE., gapw=gapw_xc) 1192 1193 ! add hartree only for SINGLETS 1194 IF (.NOT. lr_triplet) THEN 1195 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d* & 1196 v_hartree_rspace%pw%pw_grid%dvol 1197 v_rspace_new(1)%pw%cr3d = 2.0_dp*v_hartree_rspace%pw%cr3d 1198 1199 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), & 1200 pmat=rho_ao(ispin), & 1201 hmat=p_env%kpp1_env%v_ao(ispin), & 1202 qs_env=qs_env, & 1203 calculate_forces=.FALSE., gapw=gapw) 1204 END IF 1205 ELSE 1206 ! remove kpp1_env%v_ao and work directly on k_p_p1 ? 1207 CALL dbcsr_set(p_env%kpp1_env%v_ao(ispin)%matrix, 0.0_dp) 1208 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), & 1209 pmat=rho_ao(ispin), & 1210 hmat=p_env%kpp1_env%v_ao(ispin), & 1211 qs_env=qs_env, & 1212 calculate_forces=.FALSE., gapw=gapw_xc) 1213 1214 IF (ispin == 1) THEN 1215 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d* & 1216 v_hartree_rspace%pw%pw_grid%dvol 1217 END IF 1218 v_rspace_new(ispin)%pw%cr3d = v_hartree_rspace%pw%cr3d 1219 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), & 1220 pmat=rho_ao(ispin), & 1221 hmat=p_env%kpp1_env%v_ao(ispin), & 1222 qs_env=qs_env, & 1223 calculate_forces=.FALSE., gapw=gapw) 1224 END IF 1225 1226 ELSE 1227 1228 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d* & 1229 v_rspace_new(ispin)%pw%pw_grid%dvol 1230 1231 IF (nspins == 1) THEN 1232 1233 IF (.NOT. (lr_triplet)) THEN 1234 v_rspace_new(1)%pw%cr3d = 2.0_dp*v_rspace_new(1)%pw%cr3d 1235 1236 ENDIF 1237 1238 ! add hartree only for SINGLETS 1239 !IF (res_etype == tddfpt_singlet) THEN 1240 IF (.NOT. lr_triplet) THEN 1241 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d* & 1242 v_hartree_rspace%pw%pw_grid%dvol 1243 v_rspace_new(1)%pw%cr3d = v_rspace_new(1)%pw%cr3d+ & 1244 2.0_dp*v_hartree_rspace%pw%cr3d 1245 END IF 1246 ELSE 1247 IF (ispin == 1) THEN 1248 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d* & 1249 v_hartree_rspace%pw%pw_grid%dvol 1250 END IF 1251 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d+ & 1252 v_hartree_rspace%pw%cr3d 1253 END IF 1254 1255 ! remove kpp1_env%v_ao and work directly on k_p_p1 ? 1256 CALL dbcsr_set(p_env%kpp1_env%v_ao(ispin)%matrix, 0.0_dp) 1257 1258 IF (lrigpw) THEN 1259 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 1260 CALL get_qs_env(qs_env, para_env=para_env, nkind=nkind) 1261 DO ikind = 1, nkind 1262 lri_v_int(ikind)%v_int = 0.0_dp 1263 END DO 1264 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, & 1265 lri_v_int, .FALSE., "LRI_AUX") 1266 DO ikind = 1, nkind 1267 CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group) 1268 END DO 1269 ALLOCATE (k1mat(1)) 1270 k1mat(1)%matrix => p_env%kpp1_env%v_ao(ispin)%matrix 1271 IF (lri_env%exact_1c_terms) THEN 1272 CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), k1mat(1)%matrix, & 1273 rho_ao(ispin)%matrix, qs_env, .FALSE., "ORB") 1274 END IF 1275 CALL calculate_lri_ks_matrix(lri_env, lri_v_int, k1mat, atomic_kind_set) 1276 DEALLOCATE (k1mat) 1277 ELSE 1278 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), & 1279 pmat=rho_ao(ispin), & 1280 hmat=p_env%kpp1_env%v_ao(ispin), & 1281 qs_env=qs_env, & 1282 calculate_forces=.FALSE., gapw=gapw) 1283 END IF 1284 1285 END IF 1286 1287 CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, p_env%kpp1_env%v_ao(ispin)%matrix) 1288 END DO 1289 1290 IF (gapw) THEN 1291 1292 IF (.NOT. ((nspins == 1 .AND. lr_triplet))) THEN 1293 CALL Vh_1c_gg_integrals(qs_env, energy_hartree_1c, tddft=.TRUE., do_triplet=lr_triplet, & 1294 p_env=p_env) 1295 1296 CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, & 1297 .FALSE., tddft=.TRUE., do_triplet=lr_triplet, & 1298 p_env=p_env) 1299 END IF 1300 1301 ! *** Add single atom contributions to the KS matrix *** 1302 ! remap pointer 1303 ns = SIZE(p_env%kpp1) 1304 ksmat(1:ns, 1:1) => p_env%kpp1(1:ns) 1305 ns = SIZE(rho_ao) 1306 psmat(1:ns, 1:1) => rho_ao(1:ns) 1307 1308 CALL update_ks_atom(qs_env, ksmat, psmat, .FALSE., .TRUE., p_env) 1309 1310 END IF 1311 1312 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw) 1313 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw) 1314 DO ispin = 1, nspins 1315 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_new(ispin)%pw) 1316 END DO 1317 DEALLOCATE (v_rspace_new) 1318 1319 DO ispin = 1, nspins 1320 CALL cp_fm_get_info(c0(ispin)%matrix, ncol_global=ncol) 1321 CALL cp_dbcsr_sm_fm_multiply(p_env%kpp1(ispin)%matrix, & 1322 c0(ispin)%matrix, & 1323 Av(ispin)%matrix, & 1324 ncol=ncol, alpha=1.0_dp, beta=1.0_dp) 1325 ENDDO 1326 ! 1327 CALL timestop(handle) 1328 ! 1329 END SUBROUTINE apply_op_2_dft 1330 1331! ************************************************************************************************** 1332!> \brief ... 1333!> \param qs_env ... 1334!> \param p_env ... 1335!> \param c0 ... 1336!> \param v ... 1337!> \param Av ... 1338!> \param chc ... 1339! ************************************************************************************************** 1340 SUBROUTINE apply_op_2_xtb(qs_env, p_env, c0, v, Av, chc) 1341 TYPE(qs_environment_type), POINTER :: qs_env 1342 TYPE(qs_p_env_type), POINTER :: p_env 1343 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av, chc 1344 1345 CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2_xtb', routineP = moduleN//':'//routineN 1346 1347 INTEGER :: atom_a, handle, iatom, ikind, is, ispin, & 1348 na, natom, natorb, ncol, nkind, ns, & 1349 nsgf, nspins 1350 INTEGER, DIMENSION(25) :: lao 1351 INTEGER, DIMENSION(5) :: occ 1352 LOGICAL :: lr_triplet 1353 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1 1354 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1 1355 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1356 TYPE(cp_para_env_type), POINTER :: para_env 1357 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix, rho_ao 1358 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_p1, matrix_s 1359 TYPE(dbcsr_type), POINTER :: s_matrix 1360 TYPE(dft_control_type), POINTER :: dft_control 1361 TYPE(linres_control_type), POINTER :: linres_control 1362 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1363 TYPE(pw_env_type), POINTER :: pw_env 1364 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1365 TYPE(qs_rho_type), POINTER :: rho, rho1 1366 TYPE(xtb_atom_type), POINTER :: xtb_kind 1367 1368 CALL timeset(routineN, handle) 1369 1370 CPASSERT(ASSOCIATED(c0)) 1371 CPASSERT(ASSOCIATED(v)) 1372 CPASSERT(ASSOCIATED(Av)) 1373 CPASSERT(ASSOCIATED(chc)) 1374 1375 CPASSERT(ASSOCIATED(p_env%kpp1_env)) 1376 CPASSERT(ASSOCIATED(p_env%kpp1)) 1377 1378 rho1 => p_env%rho1 1379 CPASSERT(ASSOCIATED(rho1)) 1380 1381 CPASSERT(p_env%kpp1_env%ref_count > 0) 1382 1383 CALL get_qs_env(qs_env=qs_env, & 1384 pw_env=pw_env, & 1385 para_env=para_env, & 1386 rho=rho, & 1387 linres_control=linres_control, & 1388 dft_control=dft_control) 1389 1390 CALL qs_rho_get(rho, rho_ao=rho_ao) 1391 1392 lr_triplet = linres_control%lr_triplet 1393 CPASSERT(.NOT. lr_triplet) 1394 1395 nspins = SIZE(p_env%kpp1) 1396 p_env%kpp1_env%iter = p_env%kpp1_env%iter + 1 1397 1398 DO ispin = 1, nspins 1399 CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp) 1400 ENDDO 1401 1402 IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN 1403 ! Mulliken charges 1404 CALL get_qs_env(qs_env, particle_set=particle_set, matrix_s_kp=matrix_s) 1405 natom = SIZE(particle_set) 1406 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 1407 CALL qs_rho_get(rho1, rho_ao_kp=matrix_p1) 1408 ALLOCATE (mcharge(natom), charges(natom, 5)) 1409 ALLOCATE (mcharge1(natom), charges1(natom, 5)) 1410 charges = 0.0_dp 1411 charges1 = 0.0_dp 1412 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set) 1413 nkind = SIZE(atomic_kind_set) 1414 CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf) 1415 ALLOCATE (aocg(nsgf, natom)) 1416 aocg = 0.0_dp 1417 ALLOCATE (aocg1(nsgf, natom)) 1418 aocg1 = 0.0_dp 1419 p_matrix => matrix_p(:, 1) 1420 s_matrix => matrix_s(1, 1)%matrix 1421 CALL ao_charges(matrix_p, matrix_s, aocg, para_env) 1422 CALL ao_charges(matrix_p1, matrix_s, aocg1, para_env) 1423 DO ikind = 1, nkind 1424 CALL get_atomic_kind(atomic_kind_set(ikind), natom=na) 1425 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind) 1426 CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ) 1427 DO iatom = 1, na 1428 atom_a = atomic_kind_set(ikind)%atom_list(iatom) 1429 charges(atom_a, :) = REAL(occ(:), KIND=dp) 1430 DO is = 1, natorb 1431 ns = lao(is) + 1 1432 charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a) 1433 charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a) 1434 END DO 1435 END DO 1436 END DO 1437 DEALLOCATE (aocg, aocg1) 1438 DO iatom = 1, natom 1439 mcharge(iatom) = SUM(charges(iatom, :)) 1440 mcharge1(iatom) = SUM(charges1(iatom, :)) 1441 END DO 1442 ! Coulomb Kernel 1443 CALL xtb_coulomb_hessian(qs_env, p_env%kpp1, charges1, mcharge1, mcharge) 1444 ! 1445 DEALLOCATE (charges, mcharge, charges1, mcharge1) 1446 END IF 1447 1448 DO ispin = 1, nspins 1449 CALL cp_fm_get_info(c0(ispin)%matrix, ncol_global=ncol) 1450 CALL cp_dbcsr_sm_fm_multiply(p_env%kpp1(ispin)%matrix, & 1451 c0(ispin)%matrix, & 1452 Av(ispin)%matrix, & 1453 ncol=ncol, alpha=1.0_dp, beta=1.0_dp) 1454 ENDDO 1455 1456 CALL timestop(handle) 1457 1458 END SUBROUTINE apply_op_2_xtb 1459 1460! ************************************************************************************************** 1461!> \brief ... 1462!> \param p_env ... 1463!> \param qs_env ... 1464! ************************************************************************************************** 1465 SUBROUTINE p_env_check_i_alloc(p_env, qs_env) 1466 TYPE(qs_p_env_type), POINTER :: p_env 1467 TYPE(qs_environment_type), POINTER :: qs_env 1468 1469 CHARACTER(len=*), PARAMETER :: routineN = 'p_env_check_i_alloc', & 1470 routineP = moduleN//':'//routineN 1471 1472 CHARACTER(len=25) :: name 1473 INTEGER :: handle, ispin, nspins 1474 LOGICAL :: gapw_xc 1475 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 1476 TYPE(dft_control_type), POINTER :: dft_control 1477 1478 CALL timeset(routineN, handle) 1479 1480 NULLIFY (dft_control, matrix_s) 1481 1482 CPASSERT(ASSOCIATED(p_env)) 1483 CPASSERT(p_env%ref_count > 0) 1484 CALL get_qs_env(qs_env, dft_control=dft_control) 1485 gapw_xc = dft_control%qs_control%gapw_xc 1486 IF (.NOT. ASSOCIATED(p_env%kpp1)) THEN 1487 CALL get_qs_env(qs_env, matrix_s=matrix_s) 1488 nspins = dft_control%nspins 1489 1490 CALL dbcsr_allocate_matrix_set(p_env%kpp1, nspins) 1491 name = "p_env"//cp_to_string(p_env%id_nr)//"%kpp1-" 1492 !CALL compress(name,full=.TRUE.) 1493 DO ispin = 1, nspins 1494 ALLOCATE (p_env%kpp1(ispin)%matrix) 1495 CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, matrix_s(1)%matrix, & 1496 name=TRIM(name)//ADJUSTL(cp_to_string(ispin))) 1497 CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp) 1498 END DO 1499 1500 CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env) 1501 IF (gapw_xc) THEN 1502 CALL qs_rho_rebuild(p_env%rho1_xc, qs_env=qs_env) 1503 END IF 1504 END IF 1505 1506 IF (.NOT. ASSOCIATED(p_env%rho1)) THEN 1507 CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env) 1508 IF (gapw_xc) THEN 1509 CALL qs_rho_rebuild(p_env%rho1_xc, qs_env=qs_env) 1510 END IF 1511 END IF 1512 1513 CALL timestop(handle) 1514 1515 END SUBROUTINE p_env_check_i_alloc 1516 1517! ************************************************************************************************** 1518!> \brief ... 1519!> \param kpp1_env ... 1520!> \param qs_env ... 1521!> \param lr_triplet ... 1522! ************************************************************************************************** 1523 SUBROUTINE kpp1_check_i_alloc(kpp1_env, qs_env, lr_triplet) 1524 1525 TYPE(qs_kpp1_env_type), POINTER :: kpp1_env 1526 TYPE(qs_environment_type), POINTER :: qs_env 1527 LOGICAL, INTENT(IN) :: lr_triplet 1528 1529 CHARACTER(len=*), PARAMETER :: routineN = 'kpp1_check_i_alloc', & 1530 routineP = moduleN//':'//routineN 1531 1532 INTEGER :: ispin, nspins 1533 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 1534 TYPE(pw_env_type), POINTER :: pw_env 1535 TYPE(pw_p_type), DIMENSION(:), POINTER :: my_rho_r, rho_r 1536 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 1537 TYPE(qs_rho_type), POINTER :: rho 1538 TYPE(section_vals_type), POINTER :: input, xc_section 1539 1540 NULLIFY (pw_env, auxbas_pw_pool, matrix_s, rho, rho_r, input) 1541 1542 CPASSERT(ASSOCIATED(kpp1_env)) 1543 CPASSERT(kpp1_env%ref_count > 0) 1544 1545 CALL get_qs_env(qs_env, pw_env=pw_env, & 1546 matrix_s=matrix_s, input=input, rho=rho) 1547 1548 CALL qs_rho_get(rho, rho_r=rho_r) 1549 nspins = SIZE(rho_r) 1550 1551 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 1552 1553 IF (.NOT. ASSOCIATED(kpp1_env%v_rspace)) THEN 1554 ALLOCATE (kpp1_env%v_rspace(nspins)) 1555 DO ispin = 1, nspins 1556 CALL pw_pool_create_pw(auxbas_pw_pool, & 1557 kpp1_env%v_rspace(ispin)%pw, & 1558 use_data=REALDATA3D, in_space=REALSPACE) 1559 END DO 1560 END IF 1561 1562 IF (.NOT. ASSOCIATED(kpp1_env%v_ao)) THEN 1563 CALL dbcsr_allocate_matrix_set(kpp1_env%v_ao, nspins) 1564 DO ispin = 1, nspins 1565 ALLOCATE (kpp1_env%v_ao(ispin)%matrix) 1566 CALL dbcsr_copy(kpp1_env%v_ao(ispin)%matrix, matrix_s(1)%matrix, & 1567 name="kpp1%v_ao-"//ADJUSTL(cp_to_string(ispin))) 1568 END DO 1569 END IF 1570 1571 IF (.NOT. ASSOCIATED(kpp1_env%deriv_set)) THEN 1572 1573 IF (nspins == 1 .AND. lr_triplet) THEN 1574 ALLOCATE (my_rho_r(2)) 1575 DO ispin = 1, 2 1576 CALL pw_pool_create_pw(auxbas_pw_pool, my_rho_r(ispin)%pw, & 1577 use_data=rho_r(1)%pw%in_use, in_space=rho_r(1)%pw%in_space) 1578 my_rho_r(ispin)%pw%cr3d = 0.5_dp*rho_r(1)%pw%cr3d 1579 END DO 1580 ELSE 1581 ALLOCATE (my_rho_r(SIZE(rho_r))) 1582 DO ispin = 1, SIZE(rho_r) 1583 my_rho_r(ispin)%pw => rho_r(ispin)%pw 1584 CALL pw_retain(my_rho_r(ispin)%pw) 1585 END DO 1586 END IF 1587 1588 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 1589 1590 CALL xc_prep_2nd_deriv(kpp1_env%deriv_set, kpp1_env%rho_set, & 1591 my_rho_r, auxbas_pw_pool, & 1592 xc_section=xc_section) 1593 1594 DO ispin = 1, SIZE(my_rho_r) 1595 CALL pw_release(my_rho_r(ispin)%pw) 1596 ENDDO 1597 DEALLOCATE (my_rho_r) 1598 ENDIF 1599 1600 END SUBROUTINE kpp1_check_i_alloc 1601 1602! ************************************************************************************************** 1603!> \brief ... 1604!> \param v ... 1605!> \param psi0 ... 1606!> \param S_psi0 ... 1607! ************************************************************************************************** 1608 SUBROUTINE preortho(v, psi0, S_psi0) 1609 !v = (I-PS)v 1610 ! 1611 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: v, psi0, S_psi0 1612 1613 CHARACTER(LEN=*), PARAMETER :: routineN = 'preortho', routineP = moduleN//':'//routineN 1614 1615 INTEGER :: handle, ispin, mp, mt, mv, np, nspins, & 1616 nt, nv 1617 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 1618 TYPE(cp_fm_type), POINTER :: buf 1619 1620 CALL timeset(routineN, handle) 1621 ! 1622 CPASSERT(ASSOCIATED(v)) 1623 CPASSERT(ASSOCIATED(S_psi0)) 1624 CPASSERT(ASSOCIATED(psi0)) 1625 NULLIFY (buf, tmp_fm_struct) 1626 ! 1627 nspins = SIZE(v, 1) 1628 ! 1629 DO ispin = 1, nspins 1630 CALL cp_fm_get_info(v(ispin)%matrix, ncol_global=mv, nrow_global=nv) 1631 CALL cp_fm_get_info(psi0(ispin)%matrix, ncol_global=mp, nrow_global=np) 1632 ! 1633 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nv, ncol_global=mp, & 1634 para_env=v(ispin)%matrix%matrix_struct%para_env, & 1635 context=v(ispin)%matrix%matrix_struct%context) 1636 CALL cp_fm_create(buf, tmp_fm_struct) 1637 CALL cp_fm_struct_release(tmp_fm_struct) 1638 ! 1639 CALL cp_fm_get_info(buf, ncol_global=mt, nrow_global=nt) 1640 CPASSERT(nv == np) 1641 CPASSERT(mt >= mv) 1642 CPASSERT(mt >= mp) 1643 CPASSERT(nt == nv) 1644 ! 1645 ! buf = v' * S_psi0 1646 CALL cp_gemm('T', 'N', mv, mp, nv, 1.0_dp, v(ispin)%matrix, S_psi0(ispin)%matrix, 0.0_dp, buf) 1647 ! v = v - psi0 * buf' 1648 CALL cp_gemm('N', 'T', nv, mv, mp, -1.0_dp, psi0(ispin)%matrix, buf, 1.0_dp, v(ispin)%matrix) 1649 ! 1650 CALL cp_fm_release(buf) 1651 ENDDO 1652 ! 1653 CALL timestop(handle) 1654 ! 1655 END SUBROUTINE preortho 1656 1657! ************************************************************************************************** 1658!> \brief ... 1659!> \param v ... 1660!> \param psi0 ... 1661!> \param S_psi0 ... 1662! ************************************************************************************************** 1663 SUBROUTINE postortho(v, psi0, S_psi0) 1664 !v = (I-SP)v 1665 ! 1666 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: v, psi0, S_psi0 1667 1668 CHARACTER(LEN=*), PARAMETER :: routineN = 'postortho', routineP = moduleN//':'//routineN 1669 1670 INTEGER :: handle, ispin, mp, mt, mv, np, nspins, & 1671 nt, nv 1672 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 1673 TYPE(cp_fm_type), POINTER :: buf 1674 1675 CALL timeset(routineN, handle) 1676 ! 1677 CPASSERT(ASSOCIATED(v)) 1678 CPASSERT(ASSOCIATED(S_psi0)) 1679 CPASSERT(ASSOCIATED(psi0)) 1680 NULLIFY (buf, tmp_fm_struct) 1681 ! 1682 nspins = SIZE(v, 1) 1683 ! 1684 DO ispin = 1, nspins 1685 CALL cp_fm_get_info(v(ispin)%matrix, ncol_global=mv, nrow_global=nv) 1686 CALL cp_fm_get_info(psi0(ispin)%matrix, ncol_global=mp, nrow_global=np) 1687 ! 1688 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nv, ncol_global=mp, & 1689 para_env=v(ispin)%matrix%matrix_struct%para_env, & 1690 context=v(ispin)%matrix%matrix_struct%context) 1691 CALL cp_fm_create(buf, tmp_fm_struct) 1692 CALL cp_fm_struct_release(tmp_fm_struct) 1693 ! 1694 CALL cp_fm_get_info(buf, ncol_global=mt, nrow_global=nt) 1695 CPASSERT(nv == np) 1696 CPASSERT(mt >= mv) 1697 CPASSERT(mt >= mp) 1698 CPASSERT(nt == nv) 1699 ! 1700 ! buf = v' * psi0 1701 CALL cp_gemm('T', 'N', mv, mp, nv, 1.0_dp, v(ispin)%matrix, psi0(ispin)%matrix, 0.0_dp, buf) 1702 ! v = v - S_psi0 * buf' 1703 CALL cp_gemm('N', 'T', nv, mv, mp, -1.0_dp, S_psi0(ispin)%matrix, buf, 1.0_dp, v(ispin)%matrix) 1704 ! 1705 CALL cp_fm_release(buf) 1706 ENDDO 1707 ! 1708 CALL timestop(handle) 1709 ! 1710 END SUBROUTINE postortho 1711 1712! ************************************************************************************************** 1713!> \brief ... 1714!> \param qs_env ... 1715!> \param linres_section ... 1716!> \param vec ... 1717!> \param ivec ... 1718!> \param tag ... 1719!> \param ind ... 1720! ************************************************************************************************** 1721 SUBROUTINE linres_write_restart(qs_env, linres_section, vec, ivec, tag, ind) 1722 TYPE(qs_environment_type), POINTER :: qs_env 1723 TYPE(section_vals_type), POINTER :: linres_section 1724 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: vec 1725 INTEGER, INTENT(IN) :: ivec 1726 CHARACTER(LEN=*) :: tag 1727 INTEGER, INTENT(IN), OPTIONAL :: ind 1728 1729 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_write_restart', & 1730 routineP = moduleN//':'//routineN 1731 1732 CHARACTER(LEN=default_path_length) :: filename 1733 CHARACTER(LEN=default_string_length) :: my_middle, my_pos, my_status 1734 INTEGER :: handle, i, i_block, ia, ie, iounit, & 1735 ispin, j, max_block, nao, nmo, nspins, & 1736 rst_unit 1737 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer 1738 TYPE(cp_fm_type), POINTER :: mo_coeff 1739 TYPE(cp_logger_type), POINTER :: logger 1740 TYPE(cp_para_env_type), POINTER :: para_env 1741 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 1742 TYPE(section_vals_type), POINTER :: print_key 1743 1744 NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer) 1745 1746 CALL timeset(routineN, handle) 1747 1748 logger => cp_get_default_logger() 1749 1750 IF (BTEST(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", & 1751 used_print_key=print_key), & 1752 cp_p_file)) THEN 1753 1754 iounit = cp_print_key_unit_nr(logger, linres_section, & 1755 "PRINT%PROGRAM_RUN_INFO", extension=".Log") 1756 1757 CALL get_qs_env(qs_env=qs_env, & 1758 mos=mos, & 1759 para_env=para_env) 1760 1761 nspins = SIZE(mos) 1762 1763 my_status = "REPLACE" 1764 my_pos = "REWIND" 1765 CALL XSTRING(tag, ia, ie) 1766 IF (PRESENT(ind)) THEN 1767 my_middle = "RESTART-"//tag(ia:ie)//TRIM(ADJUSTL(cp_to_string(ivec))) 1768 ELSE 1769 my_middle = "RESTART-"//tag(ia:ie) 1770 IF (ivec > 1) THEN 1771 my_status = "OLD" 1772 my_pos = "APPEND" 1773 END IF 1774 END IF 1775 rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", & 1776 extension=".lr", middle_name=TRIM(my_middle), file_status=TRIM(my_status), & 1777 file_position=TRIM(my_pos), file_action="WRITE", file_form="UNFORMATTED") 1778 1779 filename = cp_print_key_generate_filename(logger, print_key, & 1780 extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.) 1781 1782 IF (iounit > 0) THEN 1783 WRITE (UNIT=iounit, FMT="(T2,A)") & 1784 "LINRES| Writing response functions to the restart file <"//TRIM(ADJUSTL(filename))//">" 1785 END IF 1786 1787 ! 1788 ! write data to file 1789 ! use the scalapack block size as a default for buffering columns 1790 CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff) 1791 CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block) 1792 ALLOCATE (vecbuffer(nao, max_block)) 1793 1794 IF (PRESENT(ind)) THEN 1795 IF (rst_unit > 0) WRITE (rst_unit) ind, ivec, nspins, nao 1796 ELSE 1797 IF (rst_unit > 0) WRITE (rst_unit) ivec, nspins, nao 1798 END IF 1799 1800 DO ispin = 1, nspins 1801 CALL cp_fm_get_info(vec(ispin)%matrix, ncol_global=nmo) 1802 1803 IF (rst_unit > 0) WRITE (rst_unit) nmo 1804 1805 DO i = 1, nmo, MAX(max_block, 1) 1806 i_block = MIN(max_block, nmo - i + 1) 1807 CALL cp_fm_get_submatrix(vec(ispin)%matrix, vecbuffer, 1, i, nao, i_block) 1808 ! doing this in one write would increase efficiency, but breaks RESTART compatibility. 1809 ! to old ones, and in cases where max_block is different between runs, as might happen during 1810 ! restarts with a different number of CPUs 1811 DO j = 1, i_block 1812 IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j) 1813 ENDDO 1814 ENDDO 1815 ENDDO 1816 1817 DEALLOCATE (vecbuffer) 1818 1819 CALL cp_print_key_finished_output(rst_unit, logger, linres_section, & 1820 "PRINT%RESTART") 1821 ENDIF 1822 1823 CALL timestop(handle) 1824 1825 END SUBROUTINE linres_write_restart 1826 1827! ************************************************************************************************** 1828!> \brief ... 1829!> \param qs_env ... 1830!> \param linres_section ... 1831!> \param vec ... 1832!> \param ivec ... 1833!> \param tag ... 1834!> \param ind ... 1835! ************************************************************************************************** 1836 SUBROUTINE linres_read_restart(qs_env, linres_section, vec, ivec, tag, ind) 1837 TYPE(qs_environment_type), POINTER :: qs_env 1838 TYPE(section_vals_type), POINTER :: linres_section 1839 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: vec 1840 INTEGER, INTENT(IN) :: ivec 1841 CHARACTER(LEN=*) :: tag 1842 INTEGER, INTENT(INOUT), OPTIONAL :: ind 1843 1844 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_read_restart', & 1845 routineP = moduleN//':'//routineN 1846 1847 CHARACTER(LEN=default_path_length) :: filename 1848 CHARACTER(LEN=default_string_length) :: my_middle 1849 INTEGER :: group, handle, i, i_block, ia, ie, iostat, iounit, ispin, iv, iv1, ivec_tmp, j, & 1850 max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, nspins, nspins_tmp, rst_unit, source 1851 LOGICAL :: file_exists 1852 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer 1853 TYPE(cp_fm_type), POINTER :: mo_coeff 1854 TYPE(cp_logger_type), POINTER :: logger 1855 TYPE(cp_para_env_type), POINTER :: para_env 1856 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 1857 TYPE(section_vals_type), POINTER :: print_key 1858 1859 file_exists = .FALSE. 1860 1861 CALL timeset(routineN, handle) 1862 1863 NULLIFY (mos, para_env, logger, print_key, vecbuffer) 1864 logger => cp_get_default_logger() 1865 1866 iounit = cp_print_key_unit_nr(logger, linres_section, & 1867 "PRINT%PROGRAM_RUN_INFO", extension=".Log") 1868 1869 CALL get_qs_env(qs_env=qs_env, & 1870 para_env=para_env, & 1871 mos=mos) 1872 1873 nspins = SIZE(mos) 1874 group = para_env%group 1875 source = para_env%source !ionode??? 1876 1877 rst_unit = -1 1878 IF (para_env%ionode) THEN 1879 CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", & 1880 n_rep_val=n_rep_val) 1881 1882 CALL XSTRING(tag, ia, ie) 1883 IF (PRESENT(ind)) THEN 1884 my_middle = "RESTART-"//tag(ia:ie)//TRIM(ADJUSTL(cp_to_string(ivec))) 1885 ELSE 1886 my_middle = "RESTART-"//tag(ia:ie) 1887 END IF 1888 1889 IF (n_rep_val > 0) THEN 1890 CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename) 1891 CALL xstring(filename, ia, ie) 1892 filename = filename(ia:ie)//TRIM(my_middle)//".lr" 1893 ELSE 1894 ! try to read from the filename that is generated automatically from the printkey 1895 print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART") 1896 filename = cp_print_key_generate_filename(logger, print_key, & 1897 extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.) 1898 ENDIF 1899 INQUIRE (FILE=filename, exist=file_exists) 1900 ! 1901 ! open file 1902 IF (file_exists) THEN 1903 CALL open_file(file_name=TRIM(filename), & 1904 file_action="READ", & 1905 file_form="UNFORMATTED", & 1906 file_position="REWIND", & 1907 file_status="OLD", & 1908 unit_number=rst_unit) 1909 1910 IF (iounit > 0) WRITE (iounit, "(T2,A)") & 1911 "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">" 1912 ELSE 1913 IF (iounit > 0) WRITE (iounit, "(T2,A)") & 1914 "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found" 1915 ENDIF 1916 ENDIF 1917 1918 CALL mp_bcast(file_exists, source, group) 1919 1920 IF (file_exists) THEN 1921 1922 CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff) 1923 CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block) 1924 1925 ALLOCATE (vecbuffer(nao, max_block)) 1926 ! 1927 ! read headers 1928 IF (PRESENT(ind)) THEN 1929 iv1 = ivec 1930 ELSE 1931 iv1 = 1 1932 END IF 1933 DO iv = iv1, ivec 1934 1935 IF (PRESENT(ind)) THEN 1936 IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) ind, ivec_tmp, nspins_tmp, nao_tmp 1937 CALL mp_bcast(iostat, source, group) 1938 CALL mp_bcast(ind, source, group) 1939 ELSE 1940 IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) ivec_tmp, nspins_tmp, nao_tmp 1941 CALL mp_bcast(iostat, source, group) 1942 END IF 1943 1944 IF (iostat .NE. 0) EXIT 1945 CALL mp_bcast(ivec_tmp, source, group) 1946 CALL mp_bcast(nspins_tmp, source, group) 1947 CALL mp_bcast(nao_tmp, source, group) 1948 1949 ! check that the number nao, nmo and nspins are 1950 ! the same as in the current mos 1951 IF (nspins_tmp .NE. nspins) CPABORT("nspins not consistent") 1952 IF (nao_tmp .NE. nao) CPABORT("nao not consistent") 1953 ! 1954 DO ispin = 1, nspins 1955 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff) 1956 CALL cp_fm_get_info(mo_coeff, ncol_global=nmo) 1957 ! 1958 IF (rst_unit > 0) READ (rst_unit) nmo_tmp 1959 CALL mp_bcast(nmo_tmp, source, group) 1960 IF (nmo_tmp .NE. nmo) CPABORT("nmo not consistent") 1961 ! 1962 ! read the response 1963 DO i = 1, nmo, MAX(max_block, 1) 1964 i_block = MIN(max_block, nmo - i + 1) 1965 DO j = 1, i_block 1966 IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j) 1967 ENDDO 1968 IF (iv .EQ. ivec_tmp) THEN 1969 CALL mp_bcast(vecbuffer, source, group) 1970 CALL cp_fm_set_submatrix(vec(ispin)%matrix, vecbuffer, 1, i, nao, i_block) 1971 ENDIF 1972 ENDDO 1973 ENDDO 1974 IF (ivec .EQ. ivec_tmp) EXIT 1975 ENDDO 1976 1977 IF (iostat /= 0) THEN 1978 IF (iounit > 0) WRITE (iounit, "(T2,A)") & 1979 "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found" 1980 ENDIF 1981 1982 DEALLOCATE (vecbuffer) 1983 1984 ENDIF 1985 1986 IF (para_env%ionode) THEN 1987 IF (file_exists) CALL close_file(unit_number=rst_unit) 1988 ENDIF 1989 1990 CALL timestop(handle) 1991 1992 END SUBROUTINE linres_read_restart 1993 1994! ************************************************************************************************** 1995 1996! ************************************************************************************************** 1997!> \brief ... 1998!> \param p_env ... 1999!> \param linres_control ... 2000!> \param nspins ... 2001! ************************************************************************************************** 2002 SUBROUTINE check_p_env_init(p_env, linres_control, nspins) 2003 ! 2004 TYPE(qs_p_env_type), POINTER :: p_env 2005 TYPE(linres_control_type), POINTER :: linres_control 2006 INTEGER, INTENT(IN) :: nspins 2007 2008 CHARACTER(LEN=*), PARAMETER :: routineN = 'check_p_env_init', & 2009 routineP = moduleN//':'//routineN 2010 2011 INTEGER :: ispin, ncol, nrow 2012 2013 p_env%iter = 0 2014 p_env%only_energy = .FALSE. 2015 p_env%ls_count = 0 2016 2017 p_env%ls_pos = 0.0_dp 2018 p_env%ls_energy = 0.0_dp 2019 p_env%ls_grad = 0.0_dp 2020 p_env%gnorm_old = 1.0_dp 2021 2022 IF (linres_control%preconditioner_type /= ot_precond_none) THEN 2023 CPASSERT(ASSOCIATED(p_env%preconditioner)) 2024 DO ispin = 1, nspins 2025 CALL cp_fm_get_info(p_env%PS_psi0(ispin)%matrix, nrow_global=nrow, ncol_global=ncol) 2026 CPASSERT(nrow == p_env%n_ao(ispin)) 2027 CPASSERT(ncol == p_env%n_mo(ispin)) 2028 ENDDO 2029 ENDIF 2030 2031 END SUBROUTINE check_p_env_init 2032 2033END MODULE qs_linres_methods 2034