1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief routines that build the Kohn-Sham matrix (i.e calculate the coulomb 8!> and xc parts 9!> \author Fawzi Mohamed 10!> \par History 11!> - 05.2002 moved from qs_scf (see there the history) [fawzi] 12!> - JGH [30.08.02] multi-grid arrays independent from density and potential 13!> - 10.2002 introduced pools, uses updated rho as input, 14!> removed most temporary variables, renamed may vars, 15!> began conversion to LSD [fawzi] 16!> - 10.2004 moved calculate_w_matrix here [Joost VandeVondele] 17!> introduced energy derivative wrt MOs [Joost VandeVondele] 18!> - SCCS implementation (16.10.2013,MK) 19! ************************************************************************************************** 20MODULE qs_ks_methods 21 USE admm_dm_methods, ONLY: admm_dm_calc_rho_aux,& 22 admm_dm_merge_ks_matrix 23 USE admm_methods, ONLY: admm_mo_calc_rho_aux,& 24 admm_mo_merge_derivs,& 25 admm_mo_merge_ks_matrix 26 USE admm_types, ONLY: admm_type 27 USE cell_types, ONLY: cell_type 28 USE cp_blacs_env, ONLY: cp_blacs_env_type 29 USE cp_control_types, ONLY: dft_control_type 30 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 31 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 32 copy_fm_to_dbcsr,& 33 cp_dbcsr_plus_fm_fm_t,& 34 cp_dbcsr_sm_fm_multiply,& 35 dbcsr_allocate_matrix_set,& 36 dbcsr_copy_columns_hack 37 USE cp_ddapc, ONLY: qs_ks_ddapc 38 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 39 cp_fm_scale_and_add,& 40 cp_fm_symm,& 41 cp_fm_transpose,& 42 cp_fm_upper_to_full 43 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 44 cp_fm_struct_release,& 45 cp_fm_struct_type 46 USE cp_fm_types, ONLY: cp_fm_create,& 47 cp_fm_get_info,& 48 cp_fm_p_type,& 49 cp_fm_release,& 50 cp_fm_to_fm,& 51 cp_fm_type 52 USE cp_gemm_interface, ONLY: cp_gemm 53 USE cp_log_handling, ONLY: cp_get_default_logger,& 54 cp_logger_get_default_io_unit,& 55 cp_logger_type 56 USE cp_output_handling, ONLY: cp_p_file,& 57 cp_print_key_should_output 58 USE cp_para_types, ONLY: cp_para_env_type 59 USE dbcsr_api, ONLY: & 60 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_filter, dbcsr_get_info, dbcsr_multiply, & 61 dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_symmetric 62 USE dft_plus_u, ONLY: plus_u 63 USE efield_utils, ONLY: efield_potential 64 USE hartree_local_methods, ONLY: Vh_1c_gg_integrals 65 USE hfx_admm_utils, ONLY: hfx_admm_init,& 66 hfx_ks_matrix 67 USE input_constants, ONLY: do_ppl_grid,& 68 outer_scf_becke_constraint,& 69 outer_scf_hirshfeld_constraint 70 USE input_section_types, ONLY: section_vals_get,& 71 section_vals_get_subs_vals,& 72 section_vals_type,& 73 section_vals_val_get 74 USE kg_correction, ONLY: kg_ekin_subset 75 USE kinds, ONLY: default_string_length,& 76 dp 77 USE lri_environment_methods, ONLY: v_int_ppl_energy 78 USE lri_environment_types, ONLY: lri_density_type,& 79 lri_environment_type,& 80 lri_kind_type 81 USE mathlib, ONLY: abnormal_value 82 USE message_passing, ONLY: mp_sum 83 USE pw_env_types, ONLY: pw_env_get,& 84 pw_env_type 85 USE pw_methods, ONLY: pw_axpy,& 86 pw_copy,& 87 pw_integral_ab,& 88 pw_integrate_function,& 89 pw_scale,& 90 pw_transfer,& 91 pw_zero 92 USE pw_poisson_methods, ONLY: pw_poisson_solve 93 USE pw_poisson_types, ONLY: pw_poisson_implicit,& 94 pw_poisson_type 95 USE pw_pool_types, ONLY: pw_pool_create_pw,& 96 pw_pool_give_back_pw,& 97 pw_pool_type 98 USE pw_types, ONLY: COMPLEXDATA1D,& 99 REALDATA3D,& 100 REALSPACE,& 101 RECIPROCALSPACE,& 102 pw_p_type,& 103 pw_release 104 USE qmmm_image_charge, ONLY: add_image_pot_to_hartree_pot,& 105 calculate_image_pot,& 106 integrate_potential_devga_rspace 107 USE qs_cdft_types, ONLY: cdft_control_type 108 USE qs_charges_types, ONLY: qs_charges_type 109 USE qs_core_energies, ONLY: calculate_ptrace 110 USE qs_dftb_matrices, ONLY: build_dftb_ks_matrix 111 USE qs_efield_berry, ONLY: qs_efield_berry_phase 112 USE qs_efield_local, ONLY: qs_efield_local_operator 113 USE qs_energy_types, ONLY: qs_energy_type 114 USE qs_environment_types, ONLY: get_qs_env,& 115 qs_environment_type 116 USE qs_gapw_densities, ONLY: prepare_gapw_den 117 USE qs_integrate_potential, ONLY: integrate_ppl_rspace,& 118 integrate_rho_nlcc,& 119 integrate_v_core_rspace 120 USE qs_ks_apply_restraints, ONLY: qs_ks_cdft_constraint,& 121 qs_ks_mulliken_restraint,& 122 qs_ks_s2_restraint 123 USE qs_ks_atom, ONLY: update_ks_atom 124 USE qs_ks_qmmm_methods, ONLY: qmmm_calculate_energy,& 125 qmmm_modify_hartree_pot 126 USE qs_ks_types, ONLY: qs_ks_env_type,& 127 set_ks_env 128 USE qs_ks_utils, ONLY: & 129 calc_v_sic_rspace, calculate_zmp_potential, compute_matrix_vxc, & 130 get_embed_potential_energy, low_spin_roks, print_densities, print_detailed_energy, & 131 sic_explicit_orbitals, sum_up_and_integrate 132 USE qs_mo_types, ONLY: get_mo_set,& 133 mo_set_p_type,& 134 mo_set_type 135 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 136 USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace 137 USE qs_rho_types, ONLY: qs_rho_get,& 138 qs_rho_type 139 USE qs_sccs, ONLY: sccs 140 USE qs_vxc, ONLY: qs_vxc_create 141 USE qs_vxc_atom, ONLY: calculate_vxc_atom 142 USE rtp_admm_methods, ONLY: rtp_admm_calc_rho_aux,& 143 rtp_admm_merge_ks_matrix 144 USE se_fock_matrix, ONLY: build_se_fock_matrix 145 USE surface_dipole, ONLY: calc_dipsurf_potential 146 USE virial_types, ONLY: virial_type 147 USE xtb_matrices, ONLY: build_xtb_ks_matrix 148#include "./base/base_uses.f90" 149 150 IMPLICIT NONE 151 152 PRIVATE 153 154 INTERFACE calculate_w_matrix 155 MODULE PROCEDURE calculate_w_matrix_1, & 156 calculate_w_matrix_roks 157 158 END INTERFACE 159 160 LOGICAL, PARAMETER :: debug_this_module = .TRUE. 161 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_methods' 162 163 PUBLIC :: calc_rho_tot_gspace, qs_ks_update_qs_env, qs_ks_build_kohn_sham_matrix, & 164 calculate_w_matrix, calculate_w_matrix_ot, qs_ks_allocate_basics 165 166 TYPE qs_potential_type 167 TYPE(pw_p_type), POINTER :: v_efield_rspace, v_hartree_gspace, & 168 vppl_rspace, v_sic_rspace, v_spin_ddapc_rest_r 169 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_new, & 170 v_rspace_new_aux_fit, & 171 v_tau_rspace, & 172 v_tau_rspace_aux_fit, & 173 v_rspace_embed 174 END TYPE qs_potential_type 175 176CONTAINS 177 178! ************************************************************************************************** 179!> \brief routine where the real calculations are made: the 180!> KS matrix is calculated 181!> \param qs_env the qs_env to update 182!> \param calculate_forces if true calculate the quantities needed 183!> to calculate the forces. Defaults to false. 184!> \param just_energy if true updates the energies but not the 185!> ks matrix. Defaults to false 186!> \param print_active ... 187!> \param ext_ks_matrix ... 188!> \par History 189!> 06.2002 moved from qs_scf to qs_ks_methods, use of ks_env 190!> new did_change scheme [fawzi] 191!> 10.2002 introduced pools, uses updated rho as input, LSD [fawzi] 192!> 10.2004 build_kohn_sham matrix now also computes the derivatives 193!> of the total energy wrt to the MO coefs, if instructed to 194!> do so. This appears useful for orbital dependent functionals 195!> where the KS matrix alone (however this might be defined) 196!> does not contain the info to construct this derivative. 197!> \author Matthias Krack 198!> \note 199!> make rho, energy and qs_charges optional, defaulting 200!> to qs_env components? 201! ************************************************************************************************** 202 SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, & 203 print_active, ext_ks_matrix) 204 TYPE(qs_environment_type), POINTER :: qs_env 205 LOGICAL, INTENT(in) :: calculate_forces, just_energy 206 LOGICAL, INTENT(IN), OPTIONAL :: print_active 207 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & 208 POINTER :: ext_ks_matrix 209 210 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_build_kohn_sham_matrix', & 211 routineP = moduleN//':'//routineN 212 213 CHARACTER(len=default_string_length) :: name 214 INTEGER :: handle, iatom, idir, img, ispin, & 215 nimages, ns, nspins 216 LOGICAL :: do_adiabatic_rescaling, do_ddapc, do_hfx, do_ppl, gapw, gapw_xc, & 217 hfx_treat_lsd_in_core, just_energy_xc, lrigpw, my_print, rigpw, use_virial 218 REAL(KIND=dp) :: ecore_ppl, edisp, ee_ener, ekin_mol, & 219 mulliken_order_p 220 REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_loc 221 TYPE(admm_type), POINTER :: admm_env 222 TYPE(cdft_control_type), POINTER :: cdft_control 223 TYPE(cell_type), POINTER :: cell 224 TYPE(cp_logger_type), POINTER :: logger 225 TYPE(cp_para_env_type), POINTER :: para_env 226 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_vxc, mo_derivs 227 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h, matrix_s, my_rho, & 228 rho_ao 229 TYPE(dft_control_type), POINTER :: dft_control 230 TYPE(lri_density_type), POINTER :: lri_density 231 TYPE(lri_environment_type), POINTER :: lri_env 232 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int 233 TYPE(pw_env_type), POINTER :: pw_env 234 TYPE(pw_p_type) :: rho_tot_gspace, v_efield_rspace, v_hartree_gspace, v_hartree_rspace, & 235 v_minus_veps, v_sccs_rspace, v_sic_rspace, v_spin_ddapc_rest_r 236 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, v_rspace_embed, v_rspace_new, & 237 v_rspace_new_aux_fit, v_tau_rspace, & 238 v_tau_rspace_aux_fit 239 TYPE(pw_p_type), POINTER :: rho0_s_rs, rho_core, rho_nlcc, vee, & 240 vppl_rspace 241 TYPE(pw_poisson_type), POINTER :: poisson_env 242 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 243 TYPE(qs_energy_type), POINTER :: energy 244 TYPE(qs_ks_env_type), POINTER :: ks_env 245 TYPE(qs_rho_type), POINTER :: rho, rho_struct, rho_xc 246 TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, & 247 hfx_sections, input, scf_section, & 248 xc_section 249 TYPE(virial_type), POINTER :: virial 250 251 CALL timeset(routineN, handle) 252 NULLIFY (admm_env, cell, dft_control, logger, & 253 mo_derivs, my_rho, rho_struct, para_env, pw_env, virial, & 254 vppl_rspace, v_hartree_rspace%pw, adiabatic_rescaling_section, & 255 hfx_sections, input, scf_section, xc_section, matrix_h, matrix_s, & 256 auxbas_pw_pool, poisson_env, v_rspace_new, & 257 v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, & 258 matrix_vxc, vee, rho_nlcc, ks_env, & 259 ks_matrix, rho, energy, rho_xc, rho_r, rho_ao, rho_core) 260 261 CPASSERT(ASSOCIATED(qs_env)) 262 263 logger => cp_get_default_logger() 264 my_print = .TRUE. 265 IF (PRESENT(print_active)) my_print = print_active 266 267 CALL get_qs_env(qs_env, & 268 ks_env=ks_env, & 269 dft_control=dft_control, & 270 matrix_h_kp=matrix_h, & 271 matrix_s_kp=matrix_s, & 272 matrix_ks_kp=ks_matrix, & 273 matrix_vxc=matrix_vxc, & 274 pw_env=pw_env, & 275 cell=cell, & 276 para_env=para_env, & 277 input=input, & 278 virial=virial, & 279 v_hartree_rspace=v_hartree_rspace%pw, & 280 vee=vee, & 281 rho_nlcc=rho_nlcc, & 282 rho=rho, & 283 rho_core=rho_core, & 284 rho_xc=rho_xc, & 285 energy=energy) 286 287 CALL qs_rho_get(rho, rho_r=rho_r, rho_ao_kp=rho_ao) 288 289 IF (PRESENT(ext_ks_matrix)) THEN 290 ! remap pointer to allow for non-kpoint external ks matrix 291 ! ext_ks_matrix is used in linear response code 292 ns = SIZE(ext_ks_matrix) 293 ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns) 294 END IF 295 296 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 297 298 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF") 299 CALL section_vals_get(hfx_sections, explicit=do_hfx) 300 IF (do_hfx) THEN 301 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, & 302 i_rep_section=1) 303 END IF 304 adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING") 305 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling) 306 just_energy_xc = just_energy 307 IF (do_adiabatic_rescaling) THEN 308 !! If we perform adiabatic rescaling, the xc potential has to be scaled by the xc- and 309 !! HFX-energy. Thus, let us first calculate the energy 310 just_energy_xc = .TRUE. 311 END IF 312 313 nimages = dft_control%nimages 314 nspins = dft_control%nspins 315 CPASSERT(ASSOCIATED(matrix_h)) 316 CPASSERT(ASSOCIATED(matrix_s)) 317 CPASSERT(ASSOCIATED(rho)) 318 CPASSERT(ASSOCIATED(pw_env)) 319 CPASSERT(SIZE(ks_matrix, 1) > 0) 320 321 ! Setup the possible usage of DDAPC charges 322 do_ddapc = dft_control%qs_control%ddapc_restraint .OR. & 323 qs_env%cp_ddapc_ewald%do_decoupling .OR. & 324 qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. & 325 qs_env%cp_ddapc_ewald%do_solvation 326 327 ! Check if LRIGPW is used 328 lrigpw = dft_control%qs_control%lrigpw 329 rigpw = dft_control%qs_control%rigpw 330 IF (rigpw) THEN 331 CPASSERT(nimages == 1) 332 ENDIF 333 IF (lrigpw .AND. rigpw) THEN 334 CPABORT(" LRI and RI are not compatible") 335 ENDIF 336 337 ! Check for GAPW method : additional terms for local densities 338 gapw = dft_control%qs_control%gapw 339 gapw_xc = dft_control%qs_control%gapw_xc 340 IF (gapw_xc .AND. gapw) THEN 341 CPABORT(" GAPW and GAPW_XC are not compatible") 342 END IF 343 IF ((gapw .AND. lrigpw) .OR. (gapw_xc .AND. lrigpw)) THEN 344 CPABORT(" GAPW/GAPW_XC and LRIGPW are not compatible") 345 ENDIF 346 IF ((gapw .AND. rigpw) .OR. (gapw_xc .AND. rigpw)) THEN 347 CPABORT(" GAPW/GAPW_XC and RIGPW are not compatible") 348 ENDIF 349 350 do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid 351 IF (do_ppl) THEN 352 CPASSERT(.NOT. gapw) 353 CALL get_qs_env(qs_env=qs_env, vppl=vppl_rspace) 354 END IF 355 356 IF (gapw_xc) THEN 357 CPASSERT(ASSOCIATED(rho_xc)) 358 END IF 359 360 ! gets the tmp grids 361 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env) 362 363 IF (gapw .AND. (poisson_env%parameters%solver .EQ. pw_poisson_implicit)) THEN 364 CPABORT("The implicit Poisson solver cannot be used in conjunction with GAPW.") 365 END IF 366 367 ! *** Prepare densities for gapw *** 368 IF (gapw .OR. gapw_xc) THEN 369 CALL prepare_gapw_den(qs_env, do_rho0=(.NOT. gapw_xc)) 370 ENDIF 371 372 ! Calculate the Hartree potential 373 CALL pw_pool_create_pw(auxbas_pw_pool, & 374 v_hartree_gspace%pw, & 375 use_data=COMPLEXDATA1D, & 376 in_space=RECIPROCALSPACE) 377 CALL pw_pool_create_pw(auxbas_pw_pool, & 378 rho_tot_gspace%pw, & 379 use_data=COMPLEXDATA1D, & 380 in_space=RECIPROCALSPACE) 381 382 scf_section => section_vals_get_subs_vals(input, "DFT%SCF") 383 IF (BTEST(cp_print_key_should_output(logger%iter_info, scf_section, & 384 "PRINT%DETAILED_ENERGY"), & 385 cp_p_file) .AND. & 386 (.NOT. gapw) .AND. (.NOT. gapw_xc) .AND. & 387 (.NOT. (poisson_env%parameters%solver .EQ. pw_poisson_implicit))) THEN 388 CALL pw_zero(rho_tot_gspace%pw) 389 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density=.TRUE.) 390 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, energy%e_hartree, & 391 v_hartree_gspace%pw) 392 CALL pw_zero(rho_tot_gspace%pw) 393 CALL pw_zero(v_hartree_gspace%pw) 394 END IF 395 396 ! Get the total density in g-space [ions + electrons] 397 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho) 398 399 IF (my_print) THEN 400 CALL print_densities(qs_env, rho) 401 END IF 402 403 IF (dft_control%do_sccs) THEN 404 ! Self-consistent continuum solvation (SCCS) model 405 CALL pw_pool_create_pw(auxbas_pw_pool, & 406 v_sccs_rspace%pw, & 407 use_data=REALDATA3D, & 408 in_space=REALSPACE) 409 410 IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN 411 CPABORT("The implicit Poisson solver cannot be used together with SCCS.") 412 END IF 413 414 IF (use_virial .AND. calculate_forces) THEN 415 CALL sccs(qs_env, rho_tot_gspace%pw, v_hartree_gspace%pw, v_sccs_rspace%pw, & 416 h_stress=h_stress) 417 virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp) 418 virial%pv_hartree = virial%pv_hartree + h_stress/REAL(para_env%num_pe, dp) 419 ELSE 420 CALL sccs(qs_env, rho_tot_gspace%pw, v_hartree_gspace%pw, v_sccs_rspace%pw) 421 END IF 422 ELSE 423 ! Getting the Hartree energy and Hartree potential. Also getting the stress tensor 424 ! from the Hartree term if needed. No nuclear force information here 425 IF (use_virial .AND. calculate_forces) THEN 426 h_stress(:, :) = 0.0_dp 427 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, energy%hartree, & 428 v_hartree_gspace%pw, h_stress=h_stress, & 429 rho_core=rho_core) 430 virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp) 431 virial%pv_hartree = virial%pv_hartree + h_stress/REAL(para_env%num_pe, dp) 432 ELSE 433 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, energy%hartree, & 434 v_hartree_gspace%pw, rho_core=rho_core) 435 END IF 436 END IF 437 438 ! In case decouple periodic images and/or apply restraints to charges 439 IF (do_ddapc) THEN 440 CALL qs_ks_ddapc(qs_env, auxbas_pw_pool, rho_tot_gspace, v_hartree_gspace, & 441 v_spin_ddapc_rest_r, energy, calculate_forces, ks_matrix, & 442 just_energy) 443 ELSE 444 dft_control%qs_control%ddapc_explicit_potential = .FALSE. 445 dft_control%qs_control%ddapc_restraint_is_spin = .FALSE. 446 IF (.NOT. just_energy) THEN 447 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw) 448 CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol) 449 END IF 450 END IF 451 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw) 452 453 IF (dft_control%apply_efield_field) THEN 454 CALL pw_pool_create_pw(auxbas_pw_pool, & 455 v_efield_rspace%pw, & 456 use_data=REALDATA3D, & 457 in_space=REALSPACE) 458 CALL efield_potential(qs_env, v_efield_rspace) 459 CALL pw_scale(v_efield_rspace%pw, v_efield_rspace%pw%pw_grid%dvol) 460 END IF 461 462 IF (dft_control%correct_surf_dip) THEN 463 CALL calc_dipsurf_potential(qs_env, energy) 464 energy%hartree = energy%hartree + energy%surf_dipole 465 END IF 466 467 ! SIC 468 CALL calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, & 469 just_energy, calculate_forces, auxbas_pw_pool) 470 471 IF (gapw) CALL Vh_1c_gg_integrals(qs_env, energy%hartree_1c) 472 473 ! Check if CDFT constraint is needed 474 CALL qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control) 475 476 ! Adds the External Potential if requested 477 IF (dft_control%apply_external_potential) THEN 478 ! Compute the energy due to the external potential 479 ee_ener = 0.0_dp 480 DO ispin = 1, nspins 481 ee_ener = ee_ener + pw_integral_ab(rho_r(ispin)%pw, vee%pw) 482 END DO 483 IF (.NOT. just_energy) THEN 484 IF (gapw) THEN 485 CALL get_qs_env(qs_env=qs_env, & 486 rho0_s_rs=rho0_s_rs) 487 CPASSERT(ASSOCIATED(rho0_s_rs)) 488 ee_ener = ee_ener + pw_integral_ab(rho0_s_rs%pw, vee%pw) 489 END IF 490 END IF 491 ! the sign accounts for the charge of the electrons 492 energy%ee = -ee_ener 493 END IF 494 495 ! Adds the QM/MM potential 496 IF (qs_env%qmmm) THEN 497 CALL qmmm_calculate_energy(qs_env=qs_env, & 498 rho=rho_r, & 499 v_qmmm=qs_env%ks_qmmm_env%v_qmmm_rspace, & 500 qmmm_energy=energy%qmmm_el) 501 IF (qs_env%qmmm_env_qm%image_charge) THEN 502 CALL calculate_image_pot(v_hartree_rspace=v_hartree_rspace, & 503 rho_hartree_gspace=rho_tot_gspace, & 504 energy=energy, & 505 qmmm_env=qs_env%qmmm_env_qm, & 506 qs_env=qs_env) 507 IF (.NOT. just_energy) THEN 508 CALL add_image_pot_to_hartree_pot(v_hartree=v_hartree_rspace, & 509 v_metal=qs_env%ks_qmmm_env%v_metal_rspace, & 510 qs_env=qs_env) 511 IF (calculate_forces) THEN 512 CALL integrate_potential_devga_rspace( & 513 potential=v_hartree_rspace, coeff=qs_env%image_coeff, & 514 forces=qs_env%qmmm_env_qm%image_charge_pot%image_forcesMM, & 515 qmmm_env=qs_env%qmmm_env_qm, qs_env=qs_env) 516 ENDIF 517 ENDIF 518 CALL pw_release(qs_env%ks_qmmm_env%v_metal_rspace%pw) 519 END IF 520 IF (.NOT. just_energy) THEN 521 CALL qmmm_modify_hartree_pot(v_hartree=v_hartree_rspace, & 522 v_qmmm=qs_env%ks_qmmm_env%v_qmmm_rspace, scale=1.0_dp) 523 END IF 524 END IF 525 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw) 526 527 ! calculate the density matrix for the fitted mo_coeffs 528 IF (dft_control%do_admm) THEN 529 CALL hfx_admm_init(qs_env) 530 IF (dft_control%do_admm_mo) THEN 531 IF (qs_env%run_rtp) THEN 532 CALL rtp_admm_calc_rho_aux(qs_env) 533 ELSE 534 CALL admm_mo_calc_rho_aux(qs_env) 535 END IF 536 ELSEIF (dft_control%do_admm_dm) THEN 537 CALL admm_dm_calc_rho_aux(ks_env) 538 ENDIF 539 ENDIF 540 541 ! only activate stress calculation if 542 IF (use_virial .AND. calculate_forces) virial%pv_calculate = .TRUE. 543 544 ! *** calculate the xc potential on the pw density *** 545 ! *** associates v_rspace_new if the xc potential needs to be computed. 546 ! If we do wavefunction fitting, we need the vxc_potential in the auxiliary basis set 547 IF (dft_control%do_admm) THEN 548 CALL get_qs_env(qs_env, admm_env=admm_env) 549 xc_section => admm_env%xc_section_aux 550 IF (gapw_xc) THEN 551 CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct) 552 ELSE 553 CALL get_qs_env(qs_env=qs_env, rho_aux_fit=rho_struct) 554 END IF 555 556 ! here we ignore a possible vdW section in admm_env%xc_section_aux 557 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, & 558 vxc_rho=v_rspace_new_aux_fit, vxc_tau=v_tau_rspace_aux_fit, exc=energy%exc_aux_fit, & 559 just_energy=just_energy_xc) 560 561 NULLIFY (rho_struct) 562 563 IF (use_virial .AND. calculate_forces) THEN 564 virial%pv_virial = virial%pv_virial - virial%pv_xc 565 ! ** virial%pv_xc will be zeroed in the xc routines 566 END IF 567 xc_section => admm_env%xc_section_primary 568 ELSE 569 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 570 END IF 571 572 IF (gapw_xc) THEN 573 CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct) 574 ELSE 575 CALL get_qs_env(qs_env=qs_env, rho=rho_struct) 576 END IF 577 578 ! zmp 579 IF (dft_control%apply_external_density .OR. dft_control%apply_external_vxc) THEN 580 energy%exc = 0.0_dp 581 CALL calculate_zmp_potential(qs_env, v_rspace_new, rho, exc=energy%exc) 582 ELSE 583 ! Embedding potential 584 IF (dft_control%apply_embed_pot) THEN 585 NULLIFY (v_rspace_embed) 586 energy%embed_corr = 0.0_dp 587 CALL get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, & 588 energy%embed_corr, just_energy) 589 ENDIF 590 ! Everything else 591 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, & 592 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, & 593 edisp=edisp, dispersion_env=qs_env%dispersion_env, & 594 just_energy=just_energy_xc) 595 IF (edisp /= 0.0_dp) energy%dispersion = edisp 596 IF (qs_env%requires_matrix_vxc .AND. ASSOCIATED(v_rspace_new)) THEN 597 CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace_new, matrix_vxc=matrix_vxc) 598 CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc) 599 ENDIF 600 601 IF (gapw .OR. gapw_xc) THEN 602 CALL calculate_vxc_atom(qs_env, just_energy_xc) 603 ENDIF 604 605 ENDIF 606 607 NULLIFY (rho_struct) 608 IF (use_virial .AND. calculate_forces) THEN 609 virial%pv_virial = virial%pv_virial - virial%pv_xc 610 virial%pv_exc = -virial%pv_xc 611 DO idir = 1, 3 612 virial%pv_exc(idir, idir) = virial%pv_exc(idir, idir) - energy%exc 613 virial%pv_hartree(idir, idir) = virial%pv_hartree(idir, idir) - 2.0_dp*energy%hartree 614 END DO 615 IF (dft_control%do_admm) THEN 616 DO idir = 1, 3 617 virial%pv_exc(idir, idir) = virial%pv_exc(idir, idir) - energy%exc_aux_fit 618 END DO 619 END IF 620 ENDIF 621 622 ! *** Add Hartree-Fock contribution if required *** 623 IF (do_hfx) THEN 624 CALL hfx_ks_matrix(qs_env, ks_matrix, rho, energy, calculate_forces, & 625 just_energy, v_rspace_new, v_tau_rspace) 626!! Adiabatic rescaling only if do_hfx; right????? 627 END IF !do_hfx 628 629 IF (do_ppl .AND. calculate_forces) THEN 630 CPASSERT(.NOT. gapw) 631 DO ispin = 1, nspins 632 CALL integrate_ppl_rspace(rho_r(ispin), qs_env) 633 END DO 634 END IF 635 636 IF (ASSOCIATED(rho_nlcc) .AND. calculate_forces) THEN 637 DO ispin = 1, nspins 638 CALL integrate_rho_nlcc(v_rspace_new(ispin), qs_env) 639 IF (dft_control%do_admm) CALL integrate_rho_nlcc(v_rspace_new_aux_fit(ispin), qs_env) 640 END DO 641 ENDIF 642 643 ! calculate KG correction 644 IF (dft_control%qs_control%do_kg .AND. just_energy) THEN 645 646 CPASSERT(.NOT. (gapw .OR. gapw_xc)) 647 CPASSERT(nimages == 1) 648 ksmat => ks_matrix(:, 1) 649 CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces) 650 651 ! substract kg corr from the total energy 652 energy%exc = energy%exc - ekin_mol 653 654 END IF 655 656 ! *** Single atom contributions *** 657 IF (.NOT. just_energy) THEN 658 IF (calculate_forces) THEN 659 ! Getting nuclear force contribution from the core charge density 660 IF ((poisson_env%parameters%solver .EQ. pw_poisson_implicit) .AND. & 661 (poisson_env%parameters%dielectric_params%dielec_core_correction)) THEN 662 CALL pw_pool_create_pw(auxbas_pw_pool, v_minus_veps%pw, use_data=REALDATA3D, in_space=REALSPACE) 663 CALL pw_copy(v_hartree_rspace%pw, v_minus_veps%pw) 664 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_minus_veps%pw, -v_hartree_rspace%pw%pw_grid%dvol) 665 CALL integrate_v_core_rspace(v_minus_veps, qs_env) 666 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_minus_veps%pw) 667 ELSE 668 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env) 669 END IF 670 END IF 671 672 IF (.NOT. do_hfx) THEN 673 ! Initialize the Kohn-Sham matrix with the core Hamiltonian matrix 674 ! (sets ks sparsity equal to matrix_h sparsity) 675 DO ispin = 1, nspins 676 DO img = 1, nimages 677 CALL dbcsr_get_info(ks_matrix(ispin, img)%matrix, name=name) ! keep the name 678 CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix, name=name) 679 END DO 680 END DO 681 END IF 682 683 IF (use_virial .AND. calculate_forces) THEN 684 pv_loc = virial%pv_virial 685 END IF 686 ! sum up potentials and integrate 687 ! Pointing my_rho to the density matrix rho_ao 688 my_rho => rho_ao 689 CALL sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, & 690 v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, & 691 v_efield_rspace, v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, & 692 cdft_control, calculate_forces) 693 694 IF (use_virial .AND. calculate_forces) THEN 695 virial%pv_hartree = virial%pv_hartree + (virial%pv_virial - pv_loc) 696 END IF 697 IF (dft_control%qs_control%do_kg) THEN 698 CPASSERT(.NOT. (gapw .OR. gapw_xc)) 699 CPASSERT(nimages == 1) 700 ksmat => ks_matrix(:, 1) 701 CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces) 702 ! substract kg corr from the total energy 703 energy%exc = energy%exc - ekin_mol 704 ! virial corrections 705 IF (use_virial .AND. calculate_forces) THEN 706 virial%pv_virial = virial%pv_virial + virial%pv_xc 707 virial%pv_xc = 0.0_dp 708 END IF 709 END IF 710 711 END IF ! .NOT. just energy 712 713 IF (dft_control%qs_control%ddapc_explicit_potential) THEN 714 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_spin_ddapc_rest_r%pw) 715 END IF 716 717 IF (dft_control%apply_efield_field) THEN 718 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_efield_rspace%pw) 719 END IF 720 721 IF (calculate_forces .AND. dft_control%qs_control%cdft) THEN 722 IF (.NOT. cdft_control%transfer_pot) THEN 723 DO iatom = 1, SIZE(cdft_control%group) 724 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%group(iatom)%weight%pw) 725 END DO 726 IF (cdft_control%atomic_charges) THEN 727 DO iatom = 1, cdft_control%natoms 728 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%charge(iatom)%pw) 729 END DO 730 DEALLOCATE (cdft_control%charge) 731 END IF 732 IF (cdft_control%type == outer_scf_becke_constraint .AND. & 733 cdft_control%becke_control%cavity_confine) THEN 734 IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN 735 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%becke_control%cavity%pw) 736 ELSE 737 DEALLOCATE (cdft_control%becke_control%cavity_mat) 738 END IF 739 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN 740 IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) & 741 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%hirshfeld_control%hirshfeld_env%fnorm%pw) 742 END IF 743 IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment) 744 cdft_control%save_pot = .FALSE. 745 cdft_control%need_pot = .TRUE. 746 cdft_control%external_control = .FALSE. 747 END IF 748 END IF 749 750 IF (dft_control%do_sccs) THEN 751 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_sccs_rspace%pw) 752 END IF 753 754 IF (gapw) THEN 755 IF (dft_control%apply_external_potential) THEN 756 ! Integrals of the Hartree potential with g0_soft 757 CALL qmmm_modify_hartree_pot(v_hartree=v_hartree_rspace, & 758 v_qmmm=vee, scale=-1.0_dp) 759 ENDIF 760 CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, calculate_forces) 761 END IF 762 763 IF (gapw .OR. gapw_xc) THEN 764 ! Single atom contributions in the KS matrix *** 765 CALL update_ks_atom(qs_env, ks_matrix, rho_ao, calculate_forces) 766 END IF 767 768 !Calculation of Mulliken restraint, if requested 769 CALL qs_ks_mulliken_restraint(energy, dft_control, just_energy, para_env, & 770 ks_matrix, matrix_s, rho, mulliken_order_p) 771 772 ! Add DFT+U contribution, if requested 773 IF (dft_control%dft_plus_u) THEN 774 CPASSERT(nimages == 1) 775 IF (just_energy) THEN 776 CALL plus_u(qs_env=qs_env) 777 ELSE 778 ksmat => ks_matrix(:, 1) 779 CALL plus_u(qs_env=qs_env, matrix_h=ksmat) 780 END IF 781 ELSE 782 energy%dft_plus_u = 0.0_dp 783 END IF 784 785 ! At this point the ks matrix should be up to date, filter it if requested 786 DO ispin = 1, nspins 787 DO img = 1, nimages 788 CALL dbcsr_filter(ks_matrix(ispin, img)%matrix, & 789 dft_control%qs_control%eps_filter_matrix) 790 END DO 791 ENDDO 792 793 !** merge the auxiliary KS matrix and the primary one 794 IF (dft_control%do_admm_mo) THEN 795 IF (qs_env%run_rtp) THEN 796 CALL rtp_admm_merge_ks_matrix(qs_env) 797 ELSE 798 CALL admm_mo_merge_ks_matrix(qs_env) 799 ENDIF 800 ELSEIF (dft_control%do_admm_dm) THEN 801 CALL admm_dm_merge_ks_matrix(ks_env) 802 ENDIF 803 804 ! External field (nonperiodic case) 805 CALL qs_efield_local_operator(qs_env, just_energy, calculate_forces) 806 807 ! Right now we can compute the orbital derivative here, as it depends currently only on the available 808 ! Kohn-Sham matrix. This might change in the future, in which case more pieces might need to be assembled 809 ! from this routine, notice that this part of the calculation in not linear scaling 810 ! right now this operation is only non-trivial because of occupation numbers and the restricted keyword 811 IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy .AND. .NOT. qs_env%run_rtp) THEN 812 CALL get_qs_env(qs_env, mo_derivs=mo_derivs) 813 CPASSERT(nimages == 1) 814 ksmat => ks_matrix(:, 1) 815 CALL calc_mo_derivatives(qs_env, ksmat, mo_derivs) 816 ENDIF 817 818 ! deal with low spin roks 819 CALL low_spin_roks(energy, qs_env, dft_control, just_energy, & 820 calculate_forces, auxbas_pw_pool) 821 822 ! deal with sic on explicit orbitals 823 CALL sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, & 824 calculate_forces, auxbas_pw_pool) 825 826 ! Periodic external field 827 CALL qs_efield_berry_phase(qs_env, just_energy, calculate_forces) 828 829 ! adds s2_restraint energy and orbital derivatives 830 CALL qs_ks_s2_restraint(dft_control, qs_env, matrix_s, & 831 energy, calculate_forces, just_energy) 832 833 IF (do_ppl) THEN 834 ! update core energy for grid based local pseudopotential 835 ecore_ppl = 0._dp 836 DO ispin = 1, nspins 837 ecore_ppl = ecore_ppl + & 838 SUM(vppl_rspace%pw%cr3d*rho_r(ispin)%pw%cr3d)*vppl_rspace%pw%pw_grid%dvol 839 END DO 840 CALL mp_sum(ecore_ppl, para_env%group) 841 energy%core = energy%core + ecore_ppl 842 END IF 843 844 IF (lrigpw) THEN 845 ! update core energy for ppl_ri method 846 CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density) 847 IF (lri_env%ppl_ri) THEN 848 ecore_ppl = 0._dp 849 DO ispin = 1, nspins 850 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 851 CALL v_int_ppl_energy(qs_env, lri_v_int, ecore_ppl) 852 END DO 853 energy%core = energy%core + ecore_ppl 854 END IF 855 END IF 856 857 ! Sum all energy terms to obtain the total energy 858 energy%total = energy%core_overlap + energy%core_self + energy%core + energy%hartree + & 859 energy%hartree_1c + energy%exc + energy%exc1 + energy%ex + & 860 energy%dispersion + energy%gcp + energy%qmmm_el + energy%mulliken + & 861 SUM(energy%ddapc_restraint) + energy%s2_restraint + & 862 energy%dft_plus_u + energy%kTS + & 863 energy%efield + energy%efield_core + energy%ee + & 864 energy%ee_core + energy%exc_aux_fit + energy%image_charge + & 865 energy%sccs_pol + energy%sccs_mpc + energy%cdft 866 867 IF (dft_control%apply_embed_pot) energy%total = energy%total + energy%embed_corr 868 869 IF (abnormal_value(energy%total)) & 870 CPABORT("KS energy is an abnormal value (NaN/Inf).") 871 872 ! Print detailed energy 873 IF (my_print) THEN 874 CALL print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p) 875 END IF 876 877 CALL timestop(handle) 878 879 END SUBROUTINE qs_ks_build_kohn_sham_matrix 880 881! ************************************************************************************************** 882!> \brief ... 883!> \param rho_tot_gspace ... 884!> \param qs_env ... 885!> \param rho ... 886!> \param skip_nuclear_density ... 887! ************************************************************************************************** 888 SUBROUTINE calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density) 889 TYPE(pw_p_type), INTENT(INOUT) :: rho_tot_gspace 890 TYPE(qs_environment_type), POINTER :: qs_env 891 TYPE(qs_rho_type), POINTER :: rho 892 LOGICAL, INTENT(IN), OPTIONAL :: skip_nuclear_density 893 894 CHARACTER(*), PARAMETER :: routineN = 'calc_rho_tot_gspace', & 895 routineP = moduleN//':'//routineN 896 897 INTEGER :: ispin 898 LOGICAL :: my_skip 899 TYPE(dft_control_type), POINTER :: dft_control 900 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g 901 TYPE(pw_p_type), POINTER :: rho0_s_gs, rho_core 902 TYPE(qs_charges_type), POINTER :: qs_charges 903 904 my_skip = .FALSE. 905 IF (PRESENT(skip_nuclear_density)) my_skip = skip_nuclear_density 906 907 CALL qs_rho_get(rho, rho_g=rho_g) 908 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 909 910 IF (.NOT. my_skip) THEN 911 NULLIFY (rho_core) 912 CALL get_qs_env(qs_env=qs_env, rho_core=rho_core) 913 IF (dft_control%qs_control%gapw) THEN 914 NULLIFY (rho0_s_gs) 915 CALL get_qs_env(qs_env=qs_env, rho0_s_gs=rho0_s_gs) 916 CPASSERT(ASSOCIATED(rho0_s_gs)) 917 CALL pw_copy(rho0_s_gs%pw, rho_tot_gspace%pw) 918 IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN 919 CALL pw_axpy(rho_core%pw, rho_tot_gspace%pw) 920 END IF 921 ELSE 922 CALL pw_copy(rho_core%pw, rho_tot_gspace%pw) 923 END IF 924 DO ispin = 1, dft_control%nspins 925 CALL pw_axpy(rho_g(ispin)%pw, rho_tot_gspace%pw) 926 END DO 927 CALL get_qs_env(qs_env=qs_env, qs_charges=qs_charges) 928 qs_charges%total_rho_gspace = pw_integrate_function(rho_tot_gspace%pw, isign=-1) 929 ELSE 930 DO ispin = 1, dft_control%nspins 931 CALL pw_axpy(rho_g(ispin)%pw, rho_tot_gspace%pw) 932 END DO 933 END IF 934 935 END SUBROUTINE calc_rho_tot_gspace 936 937! ************************************************************************************************** 938!> \brief compute MO derivatives 939!> \param qs_env the qs_env to update 940!> \param ks_matrix ... 941!> \param mo_derivs ... 942!> \par History 943!> 01.2014 created, transferred from qs_ks_build_kohn_sham_matrix in 944!> separate subroutine 945!> \author Dorothea Golze 946! ************************************************************************************************** 947 SUBROUTINE calc_mo_derivatives(qs_env, ks_matrix, mo_derivs) 948 TYPE(qs_environment_type), POINTER :: qs_env 949 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, mo_derivs 950 951 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mo_derivatives', & 952 routineP = moduleN//':'//routineN 953 954 INTEGER :: ispin 955 LOGICAL :: uniform_occupation 956 REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers 957 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mo_derivs_aux_fit, mo_derivs_tmp 958 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 959 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit 960 TYPE(dbcsr_type) :: mo_derivs2_tmp1, mo_derivs2_tmp2 961 TYPE(dbcsr_type), POINTER :: mo_coeff_b 962 TYPE(dft_control_type), POINTER :: dft_control 963 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array, mos_aux_fit 964 965 NULLIFY (dft_control, mo_array, mo_coeff, mo_coeff_aux_fit, mo_coeff_b, & 966 mo_derivs_aux_fit, mo_derivs_tmp, mos_aux_fit, & 967 occupation_numbers, matrix_ks_aux_fit) 968 969 CALL get_qs_env(qs_env, & 970 dft_control=dft_control, & 971 mos=mo_array) 972 973 IF (dft_control%do_admm_mo) THEN !fm->dbcsr 974 NULLIFY (mo_derivs_tmp) !fm->dbcsr 975 ALLOCATE (mo_derivs_tmp(SIZE(mo_derivs))) 976 DO ispin = 1, SIZE(mo_derivs) !fm->dbcsr 977 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, mo_coeff=mo_coeff) !fm->dbcsr 978 NULLIFY (mo_derivs_tmp(ispin)%matrix) 979 CALL cp_fm_create(mo_derivs_tmp(ispin)%matrix, mo_coeff%matrix_struct) !fm->dbcsr 980 ENDDO !fm->dbcsr 981 ENDIF !fm->dbcsr 982 983 DO ispin = 1, SIZE(mo_derivs) 984 985 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, mo_coeff=mo_coeff, & 986 mo_coeff_b=mo_coeff_b, occupation_numbers=occupation_numbers) 987 CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin)%matrix, mo_coeff_b, & 988 0.0_dp, mo_derivs(ispin)%matrix) 989 990 IF (dft_control%do_admm_mo) THEN 991 CALL get_qs_env(qs_env, & 992 mos_aux_fit=mos_aux_fit, & 993 mo_derivs_aux_fit=mo_derivs_aux_fit, & 994 matrix_ks_aux_fit=matrix_ks_aux_fit) 995 CALL get_mo_set(mo_set=mos_aux_fit(ispin)%mo_set, mo_coeff=mo_coeff_aux_fit) 996 997 CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, & 998 mo_derivs_tmp(ispin)%matrix) !fm->dbcsr 999 CALL admm_mo_merge_derivs(ispin, qs_env%admm_env, mo_array(ispin)%mo_set, & 1000 mo_coeff, mo_coeff_aux_fit, mo_derivs_tmp, mo_derivs_aux_fit, & 1001 matrix_ks_aux_fit) 1002 CALL copy_fm_to_dbcsr(mo_derivs_tmp(ispin)%matrix, mo_derivs(ispin)%matrix) !fm->dbcsr 1003 END IF 1004 1005 IF (dft_control%restricted) THEN 1006 ! only the first mo_set are actual variables, but we still need both 1007 CPASSERT(ispin == 1) 1008 CPASSERT(SIZE(mo_array) == 2) 1009 ! use a temporary array with the same size as the first spin for the second spin 1010 1011 ! uniform_occupation is needed for this case, otherwise we can no 1012 ! reconstruct things in ot, since we irreversibly sum 1013 CALL get_mo_set(mo_set=mo_array(1)%mo_set, uniform_occupation=uniform_occupation) 1014 CPASSERT(uniform_occupation) 1015 CALL get_mo_set(mo_set=mo_array(2)%mo_set, uniform_occupation=uniform_occupation) 1016 CPASSERT(uniform_occupation) 1017 1018 ! The beta-spin might have fewer orbitals than alpa-spin... 1019 ! create tempoary matrices with beta_nmo columns 1020 CALL get_mo_set(mo_set=mo_array(2)%mo_set, mo_coeff_b=mo_coeff_b) 1021 CALL dbcsr_create(mo_derivs2_tmp1, template=mo_coeff_b) 1022 1023 ! calculate beta derivatives 1024 CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(2)%matrix, mo_coeff_b, 0.0_dp, mo_derivs2_tmp1) 1025 1026 ! create larger matrix with alpha_nmo columns 1027 CALL dbcsr_create(mo_derivs2_tmp2, template=mo_derivs(1)%matrix) 1028 CALL dbcsr_set(mo_derivs2_tmp2, 0.0_dp) 1029 1030 ! copy into larger matrix, fills the first beta_nmo columns 1031 CALL dbcsr_copy_columns_hack(mo_derivs2_tmp2, mo_derivs2_tmp1, & 1032 mo_array(2)%mo_set%nmo, 1, 1, & 1033 para_env=mo_array(1)%mo_set%mo_coeff%matrix_struct%para_env, & 1034 blacs_env=mo_array(1)%mo_set%mo_coeff%matrix_struct%context) 1035 1036 ! add beta contribution to alpa mo_derivs 1037 CALL dbcsr_add(mo_derivs(1)%matrix, mo_derivs2_tmp2, 1.0_dp, 1.0_dp) 1038 CALL dbcsr_release(mo_derivs2_tmp1) 1039 CALL dbcsr_release(mo_derivs2_tmp2) 1040 ENDIF 1041 1042 ENDDO 1043 1044 IF (dft_control%do_admm_mo) THEN !fm->dbcsr 1045 DO ispin = 1, SIZE(mo_derivs) !fm->dbcsr 1046 CALL cp_fm_release(mo_derivs_tmp(ispin)%matrix) !fm->dbcsr 1047 ENDDO !fm->dbcsr 1048 DEALLOCATE (mo_derivs_tmp) !fm->dbcsr 1049 ENDIF !fm->dbcsr 1050 1051 END SUBROUTINE calc_mo_derivatives 1052 1053! ************************************************************************************************** 1054!> \brief updates the Kohn Sham matrix of the given qs_env (facility method) 1055!> \param qs_env the qs_env to update 1056!> \param calculate_forces if true calculate the quantities needed 1057!> to calculate the forces. Defaults to false. 1058!> \param just_energy if true updates the energies but not the 1059!> ks matrix. Defaults to false 1060!> \param print_active ... 1061!> \par History 1062!> 4.2002 created [fawzi] 1063!> 8.2014 kpoints [JGH] 1064!> 10.2014 refractored [Ole Schuett] 1065!> \author Fawzi Mohamed 1066! ************************************************************************************************** 1067 SUBROUTINE qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, & 1068 print_active) 1069 TYPE(qs_environment_type), POINTER :: qs_env 1070 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces, just_energy, & 1071 print_active 1072 1073 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_update_qs_env', & 1074 routineP = moduleN//':'//routineN 1075 1076 INTEGER :: handle, unit_nr 1077 LOGICAL :: c_forces, do_rebuild, energy_only, & 1078 forces_up_to_date, potential_changed, & 1079 rho_changed, s_mstruct_changed 1080 TYPE(qs_ks_env_type), POINTER :: ks_env 1081 1082 NULLIFY (ks_env) 1083 unit_nr = cp_logger_get_default_io_unit() 1084 1085 c_forces = .FALSE. 1086 energy_only = .FALSE. 1087 IF (PRESENT(just_energy)) energy_only = just_energy 1088 IF (PRESENT(calculate_forces)) c_forces = calculate_forces 1089 1090 IF (c_forces) THEN 1091 CALL timeset(routineN//'_forces', handle) 1092 ELSE 1093 CALL timeset(routineN, handle) 1094 ENDIF 1095 1096 CPASSERT(ASSOCIATED(qs_env)) 1097 1098 CALL get_qs_env(qs_env, & 1099 ks_env=ks_env, & 1100 rho_changed=rho_changed, & 1101 s_mstruct_changed=s_mstruct_changed, & 1102 potential_changed=potential_changed, & 1103 forces_up_to_date=forces_up_to_date) 1104 1105 do_rebuild = .FALSE. 1106 do_rebuild = do_rebuild .OR. rho_changed 1107 do_rebuild = do_rebuild .OR. s_mstruct_changed 1108 do_rebuild = do_rebuild .OR. potential_changed 1109 do_rebuild = do_rebuild .OR. (c_forces .AND. .NOT. forces_up_to_date) 1110 1111 IF (do_rebuild) THEN 1112 CALL evaluate_core_matrix_traces(qs_env) 1113 1114 ! the ks matrix will be rebuilt so this is fine now 1115 CALL set_ks_env(ks_env, potential_changed=.FALSE.) 1116 1117 CALL rebuild_ks_matrix(qs_env, & 1118 calculate_forces=c_forces, & 1119 just_energy=energy_only, & 1120 print_active=print_active) 1121 1122 IF (.NOT. energy_only) THEN 1123 CALL set_ks_env(ks_env, & 1124 rho_changed=.FALSE., & 1125 s_mstruct_changed=.FALSE., & 1126 forces_up_to_date=forces_up_to_date .OR. c_forces) 1127 END IF 1128 END IF 1129 1130 CALL timestop(handle) 1131 1132 END SUBROUTINE qs_ks_update_qs_env 1133 1134! ************************************************************************************************** 1135!> \brief Calculates the traces of the core matrices and the density matrix. 1136!> \param qs_env ... 1137!> \author Ole Schuett 1138! ************************************************************************************************** 1139 SUBROUTINE evaluate_core_matrix_traces(qs_env) 1140 TYPE(qs_environment_type), POINTER :: qs_env 1141 1142 CHARACTER(LEN=*), PARAMETER :: routineN = 'evaluate_core_matrix_traces', & 1143 routineP = moduleN//':'//routineN 1144 1145 INTEGER :: handle 1146 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_t, rho_ao_kp 1147 TYPE(dft_control_type), POINTER :: dft_control 1148 TYPE(qs_energy_type), POINTER :: energy 1149 TYPE(qs_rho_type), POINTER :: rho 1150 1151 CALL timeset(routineN, handle) 1152 NULLIFY (energy, rho, dft_control, rho_ao_kp, matrixkp_t, matrixkp_h) 1153 1154 CALL get_qs_env(qs_env, & 1155 rho=rho, & 1156 energy=energy, & 1157 dft_control=dft_control, & 1158 kinetic_kp=matrixkp_t, & 1159 matrix_h_kp=matrixkp_h) 1160 1161 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 1162 1163 CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy%core, dft_control%nspins) 1164 1165 ! kinetic energy 1166 IF (ASSOCIATED(matrixkp_t)) & 1167 CALL calculate_ptrace(matrixkp_t, rho_ao_kp, energy%kinetic, dft_control%nspins) 1168 1169 CALL timestop(handle) 1170 END SUBROUTINE evaluate_core_matrix_traces 1171 1172! ************************************************************************************************** 1173!> \brief Constructs a new Khon-Sham matrix 1174!> \param qs_env ... 1175!> \param calculate_forces ... 1176!> \param just_energy ... 1177!> \param print_active ... 1178!> \author Ole Schuett 1179! ************************************************************************************************** 1180 SUBROUTINE rebuild_ks_matrix(qs_env, calculate_forces, just_energy, print_active) 1181 TYPE(qs_environment_type), POINTER :: qs_env 1182 LOGICAL, INTENT(IN) :: calculate_forces, just_energy 1183 LOGICAL, INTENT(IN), OPTIONAL :: print_active 1184 1185 CHARACTER(LEN=*), PARAMETER :: routineN = 'rebuild_ks_matrix', & 1186 routineP = moduleN//':'//routineN 1187 1188 INTEGER :: handle 1189 TYPE(dft_control_type), POINTER :: dft_control 1190 1191 CALL timeset(routineN, handle) 1192 NULLIFY (dft_control) 1193 1194 CALL get_qs_env(qs_env, dft_control=dft_control) 1195 1196 IF (dft_control%qs_control%semi_empirical) THEN 1197 CALL build_se_fock_matrix(qs_env, & 1198 calculate_forces=calculate_forces, & 1199 just_energy=just_energy) 1200 1201 ELSEIF (dft_control%qs_control%dftb) THEN 1202 CALL build_dftb_ks_matrix(qs_env, & 1203 calculate_forces=calculate_forces, & 1204 just_energy=just_energy) 1205 1206 ELSEIF (dft_control%qs_control%xtb) THEN 1207 CALL build_xtb_ks_matrix(qs_env, & 1208 calculate_forces=calculate_forces, & 1209 just_energy=just_energy) 1210 1211 ELSE 1212 CALL qs_ks_build_kohn_sham_matrix(qs_env, & 1213 calculate_forces=calculate_forces, & 1214 just_energy=just_energy, & 1215 print_active=print_active) 1216 END IF 1217 1218 CALL timestop(handle) 1219 1220 END SUBROUTINE rebuild_ks_matrix 1221 1222! ************************************************************************************************** 1223!> \brief Calculate the W matrix from the MO eigenvectors, MO eigenvalues, 1224!> and the MO occupation numbers. Only works if they are eigenstates 1225!> \param mo_set type containing the full matrix of the MO and the eigenvalues 1226!> \param w_matrix sparse matrix 1227!> error 1228!> \par History 1229!> Creation (03.03.03,MK) 1230!> Modification that computes it as a full block, several times (e.g. 20) 1231!> faster at the cost of some additional memory 1232!> \author MK 1233! ************************************************************************************************** 1234 SUBROUTINE calculate_w_matrix_1(mo_set, w_matrix) 1235 1236 TYPE(mo_set_type), POINTER :: mo_set 1237 TYPE(dbcsr_type), POINTER :: w_matrix 1238 1239 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_1', & 1240 routineP = moduleN//':'//routineN 1241 1242 INTEGER :: handle, imo 1243 REAL(KIND=dp), DIMENSION(:), POINTER :: eigocc 1244 TYPE(cp_fm_type), POINTER :: weighted_vectors 1245 1246 CALL timeset(routineN, handle) 1247 NULLIFY (weighted_vectors) 1248 1249 CALL dbcsr_set(w_matrix, 0.0_dp) 1250 CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors") 1251 CALL cp_fm_to_fm(mo_set%mo_coeff, weighted_vectors) 1252 1253 ! scale every column with the occupation 1254 ALLOCATE (eigocc(mo_set%homo)) 1255 1256 DO imo = 1, mo_set%homo 1257 eigocc(imo) = mo_set%eigenvalues(imo)*mo_set%occupation_numbers(imo) 1258 ENDDO 1259 CALL cp_fm_column_scale(weighted_vectors, eigocc) 1260 DEALLOCATE (eigocc) 1261 1262 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=w_matrix, & 1263 matrix_v=mo_set%mo_coeff, & 1264 matrix_g=weighted_vectors, & 1265 ncol=mo_set%homo) 1266 1267 CALL cp_fm_release(weighted_vectors) 1268 1269 CALL timestop(handle) 1270 1271 END SUBROUTINE calculate_w_matrix_1 1272 1273! ************************************************************************************************** 1274!> \brief Calculate the W matrix from the MO coefs, MO derivs 1275!> could overwrite the mo_derivs for increased memory efficiency 1276!> \param mo_set type containing the full matrix of the MO coefs 1277!> mo_deriv: 1278!> \param mo_deriv ... 1279!> \param w_matrix sparse matrix 1280!> \param s_matrix sparse matrix for the overlap 1281!> error 1282!> \par History 1283!> Creation (JV) 1284!> \author MK 1285! ************************************************************************************************** 1286 SUBROUTINE calculate_w_matrix_ot(mo_set, mo_deriv, w_matrix, s_matrix) 1287 1288 TYPE(mo_set_type), POINTER :: mo_set 1289 TYPE(dbcsr_type), POINTER :: mo_deriv, w_matrix, s_matrix 1290 1291 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_ot', & 1292 routineP = moduleN//':'//routineN 1293 LOGICAL, PARAMETER :: check_gradient = .FALSE., & 1294 do_symm = .FALSE. 1295 1296 INTEGER :: handle, ncol_block, ncol_global, & 1297 nrow_block, nrow_global 1298 REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers, scaling_factor 1299 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 1300 TYPE(cp_fm_type), POINTER :: gradient, h_block, h_block_t, & 1301 weighted_vectors 1302 1303 CALL timeset(routineN, handle) 1304 NULLIFY (weighted_vectors, h_block, fm_struct_tmp) 1305 1306 CALL cp_fm_get_info(matrix=mo_set%mo_coeff, & 1307 ncol_global=ncol_global, & 1308 nrow_global=nrow_global, & 1309 nrow_block=nrow_block, & 1310 ncol_block=ncol_block) 1311 1312 CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors") 1313 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_global, ncol_global=ncol_global, & 1314 para_env=mo_set%mo_coeff%matrix_struct%para_env, & 1315 context=mo_set%mo_coeff%matrix_struct%context) 1316 CALL cp_fm_create(h_block, fm_struct_tmp, name="h block") 1317 IF (do_symm) CALL cp_fm_create(h_block_t, fm_struct_tmp, name="h block t") 1318 CALL cp_fm_struct_release(fm_struct_tmp) 1319 1320 CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers) 1321 ALLOCATE (scaling_factor(SIZE(occupation_numbers))) 1322 scaling_factor = 2.0_dp*occupation_numbers 1323 CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors) 1324 CALL cp_fm_column_scale(weighted_vectors, scaling_factor) 1325 DEALLOCATE (scaling_factor) 1326 1327 ! the convention seems to require the half here, the factor of two is presumably taken care of 1328 ! internally in qs_core_hamiltonian 1329 CALL cp_gemm('T', 'N', ncol_global, ncol_global, nrow_global, 0.5_dp, & 1330 mo_set%mo_coeff, weighted_vectors, 0.0_dp, h_block) 1331 1332 IF (do_symm) THEN 1333 ! at the minimum things are anyway symmetric, but numerically it might not be the case 1334 ! needs some investigation to find out if using this is better 1335 CALL cp_fm_transpose(h_block, h_block_t) 1336 CALL cp_fm_scale_and_add(0.5_dp, h_block, 0.5_dp, h_block_t) 1337 ENDIF 1338 1339 ! this could overwrite the mo_derivs to save the weighted_vectors 1340 CALL cp_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, & 1341 mo_set%mo_coeff, h_block, 0.0_dp, weighted_vectors) 1342 1343 CALL dbcsr_set(w_matrix, 0.0_dp) 1344 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=w_matrix, & 1345 matrix_v=mo_set%mo_coeff, & 1346 matrix_g=weighted_vectors, & 1347 ncol=mo_set%homo) 1348 1349 IF (check_gradient) THEN 1350 CALL cp_fm_create(gradient, mo_set%mo_coeff%matrix_struct, "gradient") 1351 CALL cp_dbcsr_sm_fm_multiply(s_matrix, weighted_vectors, & 1352 gradient, ncol_global) 1353 1354 ALLOCATE (scaling_factor(SIZE(occupation_numbers))) 1355 scaling_factor = 2.0_dp*occupation_numbers 1356 CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors) 1357 CALL cp_fm_column_scale(weighted_vectors, scaling_factor) 1358 DEALLOCATE (scaling_factor) 1359 1360 WRITE (*, *) " maxabs difference ", MAXVAL(ABS(weighted_vectors%local_data - 2.0_dp*gradient%local_data)) 1361 CALL cp_fm_release(gradient) 1362 ENDIF 1363 1364 IF (do_symm) CALL cp_fm_release(h_block_t) 1365 CALL cp_fm_release(weighted_vectors) 1366 CALL cp_fm_release(h_block) 1367 1368 CALL timestop(handle) 1369 1370 END SUBROUTINE calculate_w_matrix_ot 1371 1372! ************************************************************************************************** 1373!> \brief Calculate the energy-weighted density matrix W if ROKS is active. 1374!> The W matrix is returned in matrix_w. 1375!> \param mo_set ... 1376!> \param matrix_ks ... 1377!> \param matrix_p ... 1378!> \param matrix_w ... 1379!> \author 04.05.06,MK 1380! ************************************************************************************************** 1381 SUBROUTINE calculate_w_matrix_roks(mo_set, matrix_ks, matrix_p, matrix_w) 1382 TYPE(mo_set_type), POINTER :: mo_set 1383 TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_p, matrix_w 1384 1385 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_roks', & 1386 routineP = moduleN//':'//routineN 1387 1388 INTEGER :: handle, nao 1389 TYPE(cp_blacs_env_type), POINTER :: context 1390 TYPE(cp_fm_struct_type), POINTER :: fm_struct 1391 TYPE(cp_fm_type), POINTER :: c, ks, p, work 1392 TYPE(cp_para_env_type), POINTER :: para_env 1393 1394 CALL timeset(routineN, handle) 1395 1396 NULLIFY (c) 1397 NULLIFY (context) 1398 NULLIFY (fm_struct) 1399 NULLIFY (ks) 1400 NULLIFY (p) 1401 NULLIFY (para_env) 1402 NULLIFY (work) 1403 1404 CALL get_mo_set(mo_set=mo_set, mo_coeff=c) 1405 CALL cp_fm_get_info(c, context=context, nrow_global=nao, para_env=para_env) 1406 CALL cp_fm_struct_create(fm_struct, context=context, nrow_global=nao, & 1407 ncol_global=nao, para_env=para_env) 1408 CALL cp_fm_create(ks, fm_struct, name="Kohn-Sham matrix") 1409 CALL cp_fm_create(p, fm_struct, name="Density matrix") 1410 CALL cp_fm_create(work, fm_struct, name="Work matrix") 1411 CALL cp_fm_struct_release(fm_struct) 1412 CALL copy_dbcsr_to_fm(matrix_ks, ks) 1413 CALL copy_dbcsr_to_fm(matrix_p, p) 1414 CALL cp_fm_upper_to_full(p, work) 1415 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ks, p, 0.0_dp, work) 1416 CALL cp_gemm("T", "N", nao, nao, nao, 1.0_dp, p, work, 0.0_dp, ks) 1417 CALL dbcsr_set(matrix_w, 0.0_dp) 1418 CALL copy_fm_to_dbcsr(ks, matrix_w, keep_sparsity=.TRUE.) 1419 1420 CALL cp_fm_release(work) 1421 CALL cp_fm_release(p) 1422 CALL cp_fm_release(ks) 1423 1424 CALL timestop(handle) 1425 1426 END SUBROUTINE calculate_w_matrix_roks 1427 1428! ************************************************************************************************** 1429!> \brief Allocate ks_matrix and ks_env if necessary 1430!> \param qs_env ... 1431!> \par History 1432!> refactoring 04.03.2011 [MI] 1433!> \author 1434! ************************************************************************************************** 1435 1436 SUBROUTINE qs_ks_allocate_basics(qs_env) 1437 TYPE(qs_environment_type), POINTER :: qs_env 1438 1439 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_allocate_basics', & 1440 routineP = moduleN//':'//routineN 1441 1442 CHARACTER(LEN=default_string_length) :: headline 1443 INTEGER :: ic, ispin, nimg 1444 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit, & 1445 matrix_ks_aux_fit_dft, & 1446 matrix_ks_aux_fit_hfx, matrix_s_aux_fit 1447 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, matrixkp_ks 1448 TYPE(dbcsr_type), POINTER :: refmatrix 1449 TYPE(dft_control_type), POINTER :: dft_control 1450 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1451 POINTER :: sab_aux_fit, sab_orb 1452 TYPE(qs_ks_env_type), POINTER :: ks_env 1453 1454 NULLIFY (dft_control, ks_env, matrix_s_kp, matrix_s_aux_fit, sab_orb, sab_aux_fit, matrixkp_ks) 1455 1456 CALL get_qs_env(qs_env, & 1457 dft_control=dft_control, & 1458 matrix_s_kp=matrix_s_kp, & 1459 ks_env=ks_env, & 1460 sab_orb=sab_orb, & 1461 sab_aux_fit=sab_aux_fit, & 1462 matrix_ks_kp=matrixkp_ks) 1463 1464 IF (.NOT. ASSOCIATED(matrixkp_ks)) THEN 1465 nimg = dft_control%nimages 1466 CALL dbcsr_allocate_matrix_set(matrixkp_ks, dft_control%nspins, nimg) 1467 refmatrix => matrix_s_kp(1, 1)%matrix 1468 DO ispin = 1, dft_control%nspins 1469 DO ic = 1, nimg 1470 IF (dft_control%nspins > 1) THEN 1471 IF (ispin == 1) THEN 1472 headline = "KOHN-SHAM MATRIX FOR ALPHA SPIN" 1473 ELSE 1474 headline = "KOHN-SHAM MATRIX FOR BETA SPIN" 1475 END IF 1476 ELSE 1477 headline = "KOHN-SHAM MATRIX" 1478 END IF 1479 ALLOCATE (matrixkp_ks(ispin, ic)%matrix) 1480 CALL dbcsr_create(matrix=matrixkp_ks(ispin, ic)%matrix, template=refmatrix, & 1481 name=TRIM(headline), matrix_type=dbcsr_type_symmetric, nze=0) 1482 CALL cp_dbcsr_alloc_block_from_nbl(matrixkp_ks(ispin, ic)%matrix, sab_orb) 1483 CALL dbcsr_set(matrixkp_ks(ispin, ic)%matrix, 0.0_dp) 1484 ENDDO 1485 ENDDO 1486 CALL set_ks_env(ks_env, matrix_ks_kp=matrixkp_ks) 1487 END IF 1488 1489 ! Allocate matrix_ks_aux_fit if requested and put it in the QS environment 1490 ! Also matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx for ADMMS 1491 IF (dft_control%do_admm) THEN 1492 NULLIFY (matrix_ks_aux_fit, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx) 1493 CALL get_qs_env(qs_env, & 1494 matrix_ks_aux_fit=matrix_ks_aux_fit, & 1495 matrix_ks_aux_fit_dft=matrix_ks_aux_fit_dft, & 1496 matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx, & 1497 matrix_s_aux_fit=matrix_s_aux_fit) 1498 IF (.NOT. ASSOCIATED(matrix_ks_aux_fit)) THEN 1499 refmatrix => matrix_s_aux_fit(1)%matrix 1500 CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit, dft_control%nspins) 1501 DO ispin = 1, dft_control%nspins 1502 ALLOCATE (matrix_ks_aux_fit(ispin)%matrix) 1503 CALL dbcsr_create(matrix=matrix_ks_aux_fit(ispin)%matrix, template=refmatrix, & 1504 name="KOHN-SHAM_MATRIX for ADMM", & 1505 matrix_type=dbcsr_type_symmetric, nze=0) 1506 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit(ispin)%matrix, sab_aux_fit) 1507 CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp) 1508 ENDDO 1509 CALL set_ks_env(ks_env, matrix_ks_aux_fit=matrix_ks_aux_fit) 1510 END IF 1511 ! same allocation procedure for matrix_ks_aux_fit_dft 1512 IF (.NOT. ASSOCIATED(matrix_ks_aux_fit_dft)) THEN 1513 refmatrix => matrix_s_aux_fit(1)%matrix 1514 CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_dft, dft_control%nspins) 1515 DO ispin = 1, dft_control%nspins 1516 ALLOCATE (matrix_ks_aux_fit_dft(ispin)%matrix) 1517 CALL dbcsr_create(matrix=matrix_ks_aux_fit_dft(ispin)%matrix, template=refmatrix, & 1518 name="AUX. KOHN-SHAM MATRIX FOR ADMM ONLY DFT EXCHANGE", & 1519 matrix_type=dbcsr_type_symmetric, nze=0) 1520 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_dft(ispin)%matrix, sab_aux_fit) 1521 CALL dbcsr_set(matrix_ks_aux_fit_dft(ispin)%matrix, 0.0_dp) 1522 ENDDO 1523 CALL set_ks_env(ks_env, matrix_ks_aux_fit_dft=matrix_ks_aux_fit_dft) 1524 END IF 1525 ! same allocation procedure for matrix_ks_aux_fit_hfx 1526 IF (.NOT. ASSOCIATED(matrix_ks_aux_fit_hfx)) THEN 1527 refmatrix => matrix_s_aux_fit(1)%matrix 1528 CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_hfx, dft_control%nspins) 1529 DO ispin = 1, dft_control%nspins 1530 ALLOCATE (matrix_ks_aux_fit_hfx(ispin)%matrix) 1531 CALL dbcsr_create(matrix=matrix_ks_aux_fit_hfx(ispin)%matrix, template=refmatrix, & 1532 name="AUX. KOHN-SHAM MATRIX FOR ADMM ONLY HF EXCHANGE", & 1533 matrix_type=dbcsr_type_symmetric, nze=0) 1534 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_hfx(ispin)%matrix, sab_aux_fit) 1535 CALL dbcsr_set(matrix_ks_aux_fit_hfx(ispin)%matrix, 0.0_dp) 1536 ENDDO 1537 CALL set_ks_env(ks_env, matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx) 1538 END IF 1539 1540 END IF 1541 1542 END SUBROUTINE qs_ks_allocate_basics 1543 1544END MODULE qs_ks_methods 1545