1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for a Harris type energy correction on top of a 8!> Kohn-Sham calculation 9!> \par History 10!> 03.2014 created 11!> 09.2019 Moved from KG to Kohn-Sham 12!> \author JGH 13! ************************************************************************************************** 14MODULE energy_corrections 15 USE atomic_kind_types, ONLY: atomic_kind_type,& 16 get_atomic_kind,& 17 get_atomic_kind_set 18 USE basis_set_types, ONLY: get_gto_basis_set,& 19 gto_basis_set_type 20 USE cell_types, ONLY: cell_type,& 21 pbc 22 USE core_ppl, ONLY: build_core_ppl 23 USE core_ppnl, ONLY: build_core_ppnl 24 USE cp_blacs_env, ONLY: cp_blacs_env_type 25 USE cp_control_types, ONLY: dft_control_type 26 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 27 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 28 copy_fm_to_dbcsr,& 29 cp_dbcsr_sm_fm_multiply,& 30 dbcsr_allocate_matrix_set,& 31 dbcsr_deallocate_matrix_set 32 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 33 cp_fm_triangular_invert 34 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose 35 USE cp_fm_diag, ONLY: choose_eigv_solver 36 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 37 cp_fm_struct_release,& 38 cp_fm_struct_type 39 USE cp_fm_types, ONLY: cp_fm_create,& 40 cp_fm_get_info,& 41 cp_fm_p_type,& 42 cp_fm_release,& 43 cp_fm_set_all,& 44 cp_fm_to_fm_triangular,& 45 cp_fm_type 46 USE cp_log_handling, ONLY: cp_get_default_logger,& 47 cp_logger_get_default_unit_nr,& 48 cp_logger_type 49 USE cp_output_handling, ONLY: cp_p_file,& 50 cp_print_key_finished_output,& 51 cp_print_key_should_output,& 52 cp_print_key_unit_nr 53 USE cp_para_types, ONLY: cp_para_env_type 54 USE dbcsr_api, ONLY: & 55 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, & 56 dbcsr_dot, dbcsr_filter, dbcsr_get_info, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, & 57 dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, & 58 dbcsr_type_symmetric 59 USE distribution_1d_types, ONLY: distribution_1d_type 60 USE distribution_2d_types, ONLY: distribution_2d_type 61 USE dm_ls_scf_methods, ONLY: density_matrix_sign,& 62 density_matrix_tc2,& 63 density_matrix_trs4 64 USE ec_efield_local, ONLY: ec_efield_integrals,& 65 ec_efield_local_operator 66 USE ec_env_types, ONLY: energy_correction_type 67 USE external_potential_types, ONLY: get_potential,& 68 gth_potential_type,& 69 sgp_potential_type 70 USE input_constants, ONLY: ec_curvy_steps,& 71 ec_diagonalization,& 72 ec_functional_harris,& 73 ec_matrix_sign,& 74 ec_matrix_tc2,& 75 ec_matrix_trs4,& 76 ls_s_sqrt_ns,& 77 ls_s_sqrt_proot 78 USE input_section_types, ONLY: section_get_ival,& 79 section_get_lval,& 80 section_get_rval,& 81 section_vals_get_subs_vals,& 82 section_vals_type,& 83 section_vals_val_get 84 USE iterate_matrix, ONLY: invert_Hotelling,& 85 matrix_sqrt_Newton_Schulz,& 86 matrix_sqrt_proot 87 USE kinds, ONLY: default_path_length,& 88 default_string_length,& 89 dp 90 USE mao_basis, ONLY: mao_generate_basis 91 USE message_passing, ONLY: mp_sum 92 USE molecule_types, ONLY: molecule_type 93 USE moments_utils, ONLY: get_reference_point 94 USE particle_types, ONLY: particle_type 95 USE physcon, ONLY: debye 96 USE pw_env_types, ONLY: pw_env_get,& 97 pw_env_type 98 USE pw_methods, ONLY: pw_axpy,& 99 pw_integral_ab,& 100 pw_scale,& 101 pw_transfer,& 102 pw_zero 103 USE pw_poisson_methods, ONLY: pw_poisson_solve 104 USE pw_poisson_types, ONLY: pw_poisson_type 105 USE pw_pool_types, ONLY: pw_pool_create_pw,& 106 pw_pool_give_back_pw,& 107 pw_pool_type 108 USE pw_types, ONLY: COMPLEXDATA1D,& 109 REALDATA3D,& 110 REALSPACE,& 111 RECIPROCALSPACE,& 112 pw_p_type 113 USE qs_collocate_density, ONLY: calculate_rho_elec 114 USE qs_core_energies, ONLY: calculate_ecore_overlap 115 USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot 116 USE qs_energy_types, ONLY: qs_energy_type 117 USE qs_environment_types, ONLY: get_qs_env,& 118 qs_environment_type,& 119 set_qs_env 120 USE qs_force_types, ONLY: allocate_qs_force,& 121 qs_force_type,& 122 sum_qs_force,& 123 total_qs_force,& 124 zero_qs_force 125 USE qs_integrate_potential, ONLY: integrate_v_core_rspace,& 126 integrate_v_rspace 127 USE qs_kind_types, ONLY: get_qs_kind,& 128 get_qs_kind_set,& 129 qs_kind_type 130 USE qs_kinetic, ONLY: build_kinetic_matrix 131 USE qs_ks_methods, ONLY: calc_rho_tot_gspace 132 USE qs_ks_types, ONLY: qs_ks_env_type 133 USE qs_mo_types, ONLY: get_mo_set,& 134 mo_set_p_type 135 USE qs_moments, ONLY: build_local_moment_matrix 136 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 137 USE qs_neighbor_lists, ONLY: atom2d_build,& 138 atom2d_cleanup,& 139 build_neighbor_lists,& 140 local_atoms_type,& 141 pair_radius_setup 142 USE qs_overlap, ONLY: build_overlap_matrix 143 USE qs_rho_types, ONLY: qs_rho_get,& 144 qs_rho_type 145 USE qs_vxc, ONLY: qs_vxc_create 146 USE response_solver, ONLY: ks_ref_potential,& 147 response_equation,& 148 response_force 149 USE task_list_methods, ONLY: generate_qs_task_list 150 USE task_list_types, ONLY: allocate_task_list,& 151 deallocate_task_list 152 USE virial_types, ONLY: virial_type 153 USE xc, ONLY: xc_calc_2nd_deriv,& 154 xc_prep_2nd_deriv 155 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 156 xc_dset_release 157 USE xc_derivatives, ONLY: xc_functionals_get_needs 158 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 159 USE xc_rho_set_types, ONLY: xc_rho_set_create,& 160 xc_rho_set_release,& 161 xc_rho_set_type,& 162 xc_rho_set_update 163#include "./base/base_uses.f90" 164 165 IMPLICIT NONE 166 167 PRIVATE 168 169! *** Global parameters *** 170 171 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'energy_corrections' 172 173 LOGICAL, PARAMETER :: debug_forces = .TRUE. 174 175 PUBLIC :: energy_correction 176 177CONTAINS 178 179! ************************************************************************************************** 180!> \brief Energy correction to a KG simulation 181!> 182!> \param qs_env ... 183!> \param ec_init ... 184!> \param calculate_forces ... 185!> \par History 186!> 03.2014 created 187!> \author JGH 188! ************************************************************************************************** 189 SUBROUTINE energy_correction(qs_env, ec_init, calculate_forces) 190 TYPE(qs_environment_type), POINTER :: qs_env 191 LOGICAL, INTENT(IN), OPTIONAL :: ec_init, calculate_forces 192 193 CHARACTER(len=*), PARAMETER :: routineN = 'energy_correction', & 194 routineP = moduleN//':'//routineN 195 196 INTEGER :: handle, nkind, unit_nr 197 INTEGER, ALLOCATABLE, DIMENSION(:) :: natom_of_kind 198 LOGICAL :: my_calc_forces 199 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 200 TYPE(cp_logger_type), POINTER :: logger 201 TYPE(energy_correction_type), POINTER :: ec_env 202 TYPE(qs_energy_type), POINTER :: energy 203 TYPE(qs_force_type), DIMENSION(:), POINTER :: ec_force, ks_force 204 205 CALL timeset(routineN, handle) 206 207 ! Check for energy correction 208 IF (qs_env%energy_correction) THEN 209 NULLIFY (ec_env) 210 CALL get_qs_env(qs_env, ec_env=ec_env) 211 212 ec_env%should_update = .TRUE. 213 IF (PRESENT(ec_init)) ec_env%should_update = ec_init 214 215 my_calc_forces = .FALSE. 216 IF (PRESENT(calculate_forces)) my_calc_forces = calculate_forces 217 218 logger => cp_get_default_logger() 219 IF (logger%para_env%ionode) THEN 220 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 221 ELSE 222 unit_nr = -1 223 ENDIF 224 225 IF (ec_env%should_update) THEN 226 ec_env%etotal = 0.0_dp 227 ec_env%eband = 0.0_dp 228 ec_env%ehartree = 0.0_dp 229 ec_env%exc = 0.0_dp 230 ec_env%vhxc = 0.0_dp 231 ec_env%edispersion = 0.0_dp 232 END IF 233 234 IF (my_calc_forces) THEN 235 IF (unit_nr > 0) THEN 236 WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 25), & 237 " Energy Correction Forces ", REPEAT("-", 26), "!" 238 END IF 239 CALL get_qs_env(qs_env, force=ks_force, atomic_kind_set=atomic_kind_set) 240 nkind = SIZE(atomic_kind_set) 241 ALLOCATE (natom_of_kind(nkind)) 242 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind) 243 NULLIFY (ec_force) 244 CALL allocate_qs_force(ec_force, natom_of_kind) 245 DEALLOCATE (natom_of_kind) 246 CALL zero_qs_force(ec_force) 247 CALL set_qs_env(qs_env, force=ec_force) 248 ELSE 249 IF (unit_nr > 0) THEN 250 WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 29), & 251 " Energy Correction ", REPEAT("-", 29), "!" 252 END IF 253 END IF 254 255 CALL get_qs_env(qs_env, energy=energy) 256 SELECT CASE (ec_env%energy_functional) 257 CASE (ec_functional_harris) 258 ! 259 CALL harris_energy(qs_env, ec_env, my_calc_forces, unit_nr) 260 ! 261 IF (ec_env%should_update) THEN 262 energy%nonscf_correction = ec_env%etotal - energy%total 263 END IF 264 IF (my_calc_forces) THEN 265 CALL get_qs_env(qs_env, force=ec_force) 266 ec_env%force => ec_force 267 CALL zero_qs_force(ks_force) 268 CALL sum_qs_force(ks_force, ec_force) 269 CALL set_qs_env(qs_env, force=ks_force) 270 END IF 271 energy%total = energy%total + energy%nonscf_correction 272 CASE DEFAULT 273 CPABORT("unknown energy correction") 274 END SELECT 275 276 IF (.NOT. my_calc_forces .AND. unit_nr > 0) THEN 277 WRITE (unit_nr, '(T3,A,T65,F16.10)') "Energy Correction ", energy%nonscf_correction 278 END IF 279 IF (unit_nr > 0) THEN 280 WRITE (unit_nr, '(T2,A,A,A)') "!", REPEAT("-", 77), "!" 281 END IF 282 283 END IF 284 285 CALL timestop(handle) 286 287 END SUBROUTINE energy_correction 288 289! ************************************************************************************************** 290!> \brief Harris Energy Correction to a Kohn-Sham simulation 291!> 292!> \param qs_env ... 293!> \param ec_env ... 294!> \param calculate_forces ... 295!> \param unit_nr ... 296!> \par History 297!> 03.2014 created 298!> \author JGH 299! ************************************************************************************************** 300 SUBROUTINE harris_energy(qs_env, ec_env, calculate_forces, unit_nr) 301 TYPE(qs_environment_type), POINTER :: qs_env 302 TYPE(energy_correction_type), POINTER :: ec_env 303 LOGICAL, INTENT(IN) :: calculate_forces 304 INTEGER, INTENT(IN) :: unit_nr 305 306 CHARACTER(len=*), PARAMETER :: routineN = 'harris_energy', routineP = moduleN//':'//routineN 307 308 REAL(KIND=dp) :: exc 309 310 IF (ec_env%should_update) THEN 311 CALL ec_build_neighborlist(qs_env, ec_env) 312 ! 313 CALL ec_build_core_hamiltonian(qs_env, ec_env) 314 CALL ks_ref_potential(qs_env, ec_env%vh_rspace, ec_env%vxc_rspace, ec_env%vtau_rspace, & 315 ec_env%ehartree, exc) 316 CALL ec_build_ks_matrix(qs_env, ec_env) 317 ! 318 IF (ec_env%mao) THEN 319 ! MAO basis 320 IF (ASSOCIATED(ec_env%mao_coef)) CALL dbcsr_deallocate_matrix_set(ec_env%mao_coef) 321 NULLIFY (ec_env%mao_coef) 322 CALL mao_generate_basis(qs_env, ec_env%mao_coef, ref_basis_set="HARRIS", molecular=.TRUE., & 323 max_iter=ec_env%mao_max_iter, eps_grad=ec_env%mao_eps_grad, unit_nr=unit_nr) 324 END IF 325 ! 326 CALL ec_ks_solver(qs_env, ec_env) 327 ! 328 ! dispersion through pairpotentials 329 CALL ec_disp(qs_env, ec_env, calculate_forces=.FALSE.) 330 CALL ec_energy(ec_env, unit_nr) 331 END IF 332 IF (calculate_forces) THEN 333 CALL ec_disp(qs_env, ec_env, calculate_forces=.TRUE.) 334 CALL ec_build_core_hamiltonian_force(qs_env, ec_env) 335 CALL ec_build_ks_matrix_force(qs_env, ec_env) 336 ! 337 CALL response_equation(qs_env, ec_env%p_env, ec_env%cpmos) 338 ! 339 CALL response_force(qs_env, ec_env%p_env, & 340 ec_env%vh_rspace, ec_env%vxc_rspace, ec_env%vtau_rspace, & 341 ec_env%matrix_hz) 342 ! 343 CALL ec_properties(qs_env, ec_env) 344 END IF 345 346 END SUBROUTINE harris_energy 347 348! ************************************************************************************************** 349!> \brief ... 350!> \param qs_env ... 351!> \param ec_env ... 352!> \param calculate_forces ... 353! ************************************************************************************************** 354 SUBROUTINE ec_disp(qs_env, ec_env, calculate_forces) 355 TYPE(qs_environment_type), POINTER :: qs_env 356 TYPE(energy_correction_type), POINTER :: ec_env 357 LOGICAL, INTENT(IN) :: calculate_forces 358 359 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_disp', routineP = moduleN//':'//routineN 360 361 REAL(KIND=dp) :: edisp 362 363 CALL calculate_dispersion_pairpot(qs_env, ec_env%dispersion_env, edisp, calculate_forces) 364 IF (.NOT. calculate_forces) ec_env%edispersion = ec_env%edispersion + edisp 365 366 END SUBROUTINE ec_disp 367 368! ************************************************************************************************** 369!> \brief Construction of the Core Hamiltonian Matrix 370!> Short version of qs_core_hamiltonian 371!> \param qs_env ... 372!> \param ec_env ... 373!> \author Creation (03.2014,JGH) 374! ************************************************************************************************** 375 SUBROUTINE ec_build_core_hamiltonian(qs_env, ec_env) 376 TYPE(qs_environment_type), POINTER :: qs_env 377 TYPE(energy_correction_type), POINTER :: ec_env 378 379 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_core_hamiltonian', & 380 routineP = moduleN//':'//routineN 381 382 INTEGER :: handle, iounit, nder, nimages 383 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 384 LOGICAL :: calculate_forces, use_virial 385 REAL(KIND=dp) :: eps_ppnl 386 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 387 TYPE(cp_logger_type), POINTER :: logger 388 TYPE(cp_para_env_type), POINTER :: para_env 389 TYPE(dft_control_type), POINTER :: dft_control 390 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 391 POINTER :: sab_orb, sac_ppl, sap_ppnl 392 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 393 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 394 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 395 TYPE(qs_ks_env_type), POINTER :: ks_env 396 TYPE(virial_type), POINTER :: virial 397 398 CALL timeset(routineN, handle) 399 400 logger => cp_get_default_logger() 401 iounit = cp_logger_get_default_unit_nr(logger) 402 403 ! no k-points possible 404 CALL get_qs_env(qs_env=qs_env, para_env=para_env, dft_control=dft_control) 405 nimages = dft_control%nimages 406 IF (nimages /= 1) THEN 407 CPABORT("K-points for Harris functional not implemented") 408 END IF 409 410 ! check for GAPW/GAPW_XC 411 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN 412 CPABORT("Harris functional for GAPW not implemented") 413 END IF 414 415 ! get neighbor lists, we need the full sab_orb list from the ec_env 416 NULLIFY (sab_orb, sac_ppl, sap_ppnl) 417 sab_orb => ec_env%sab_orb 418 sac_ppl => ec_env%sac_ppl 419 sap_ppnl => ec_env%sap_ppnl 420 421 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env) 422 nder = 0 423 ! Overlap and kinetic energy matrices 424 CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, & 425 matrix_name="OVERLAP MATRIX", & 426 basis_type_a="HARRIS", & 427 basis_type_b="HARRIS", & 428 sab_nl=sab_orb) 429 CALL build_kinetic_matrix(ks_env, matrixkp_t=ec_env%matrix_t, & 430 matrix_name="KINETIC ENERGY MATRIX", & 431 basis_type="HARRIS", & 432 sab_nl=sab_orb) 433 434 ! initialize H matrix 435 CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1) 436 ALLOCATE (ec_env%matrix_h(1, 1)%matrix) 437 CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix) 438 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb) 439 440 ! add kinetic energy 441 CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, ec_env%matrix_t(1, 1)%matrix, & 442 keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX") 443 444 ! compute the ppl contribution to the core hamiltonian 445 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, & 446 atomic_kind_set=atomic_kind_set) 447 NULLIFY (cell_to_index, virial) 448 use_virial = .FALSE. 449 calculate_forces = .FALSE. 450 IF (ASSOCIATED(sac_ppl)) THEN 451 CALL build_core_ppl(ec_env%matrix_h, ec_env%matrix_p, force, & 452 virial, calculate_forces, use_virial, nder, & 453 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & 454 nimages, cell_to_index, "HARRIS") 455 END IF 456 457 ! compute the ppnl contribution to the core hamiltonian *** 458 eps_ppnl = dft_control%qs_control%eps_ppnl 459 IF (ASSOCIATED(sap_ppnl)) THEN 460 CALL build_core_ppnl(ec_env%matrix_h, ec_env%matrix_p, force, & 461 virial, calculate_forces, use_virial, nder, & 462 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, & 463 nimages, cell_to_index, "HARRIS") 464 END IF 465 466 ! External field (nonperiodic case) 467 ec_env%efield_nuclear = 0.0_dp 468 CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces) 469 470 CALL timestop(handle) 471 472 END SUBROUTINE ec_build_core_hamiltonian 473 474! ************************************************************************************************** 475!> \brief Solve KS equation for a given matrix 476!> \brief calculate the complete KS matrix 477!> \param qs_env ... 478!> \param ec_env ... 479!> \par History 480!> 03.2014 adapted from qs_ks_build_kohn_sham_matrix [JGH] 481!> \author JGH 482! ************************************************************************************************** 483 SUBROUTINE ec_build_ks_matrix(qs_env, ec_env) 484 TYPE(qs_environment_type), POINTER :: qs_env 485 TYPE(energy_correction_type), POINTER :: ec_env 486 487 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_ks_matrix', & 488 routineP = moduleN//':'//routineN 489 490 CHARACTER(LEN=default_string_length) :: headline 491 INTEGER :: handle, iounit, ispin, nspins 492 LOGICAL :: calculate_forces, use_virial 493 REAL(dp) :: eexc, evhxc 494 TYPE(cp_blacs_env_type), POINTER :: blacs_env 495 TYPE(cp_logger_type), POINTER :: logger 496 TYPE(cp_para_env_type), POINTER :: para_env 497 TYPE(dft_control_type), POINTER :: dft_control 498 TYPE(pw_env_type), POINTER :: pw_env 499 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, tau_r, v_rspace, v_tau_rspace 500 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 501 TYPE(qs_ks_env_type), POINTER :: ks_env 502 TYPE(qs_rho_type), POINTER :: rho 503 TYPE(virial_type), POINTER :: virial 504 505 CALL timeset(routineN, handle) 506 507 logger => cp_get_default_logger() 508 iounit = cp_logger_get_default_unit_nr(logger) 509 510 calculate_forces = .FALSE. 511 512 ! get all information on the electronic density 513 NULLIFY (rho, ks_env) 514 CALL get_qs_env(qs_env=qs_env, rho=rho, virial=virial, dft_control=dft_control, & 515 para_env=para_env, blacs_env=blacs_env, ks_env=ks_env) 516 517 nspins = dft_control%nspins 518 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 519 CPASSERT(.NOT. use_virial) 520 521 ! Kohn-Sham matrix 522 IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks) 523 CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1) 524 DO ispin = 1, nspins 525 headline = "KOHN-SHAM MATRIX" 526 ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix) 527 CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=TRIM(headline), & 528 template=ec_env%matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) 529 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%sab_orb) 530 CALL dbcsr_set(ec_env%matrix_ks(ispin, 1)%matrix, 0.0_dp) 531 ENDDO 532 533 NULLIFY (pw_env) 534 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 535 CPASSERT(ASSOCIATED(pw_env)) 536 537 ! v_rspace and v_tau_rspace are generated from the auxbas pool 538 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 539 NULLIFY (v_rspace, v_tau_rspace) 540 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, & 541 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.FALSE.) 542 IF (.NOT. ASSOCIATED(v_rspace)) THEN 543 ALLOCATE (v_rspace(nspins)) 544 DO ispin = 1, nspins 545 NULLIFY (v_rspace(ispin)%pw) 546 CALL pw_pool_create_pw(auxbas_pw_pool, v_rspace(ispin)%pw, & 547 use_data=REALDATA3D, in_space=REALSPACE) 548 CALL pw_zero(v_rspace(ispin)%pw) 549 ENDDO 550 END IF 551 552 evhxc = 0.0_dp 553 CALL qs_rho_get(rho, rho_r=rho_r) 554 IF (ASSOCIATED(v_tau_rspace)) THEN 555 CALL qs_rho_get(rho, tau_r=tau_r) 556 END IF 557 DO ispin = 1, nspins 558 ! Add v_hartree + v_xc = v_rspace 559 CALL pw_scale(v_rspace(ispin)%pw, v_rspace(ispin)%pw%pw_grid%dvol) 560 CALL pw_axpy(ec_env%vh_rspace%pw, v_rspace(ispin)%pw) 561 ! integrate over potential <a|V|b> 562 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), & 563 hmat=ec_env%matrix_ks(ispin, 1), & 564 qs_env=qs_env, & 565 calculate_forces=.FALSE., & 566 basis_type="HARRIS", & 567 task_list_external=ec_env%task_list) 568 569 IF (ASSOCIATED(v_tau_rspace)) THEN 570 ! integrate over Tau-potential <nabla.a|V|nabla.b> 571 CALL pw_scale(v_tau_rspace(ispin)%pw, v_tau_rspace(ispin)%pw%pw_grid%dvol) 572 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), hmat=ec_env%matrix_ks(ispin, 1), & 573 qs_env=qs_env, calculate_forces=.FALSE., compute_tau=.TRUE., & 574 basis_type="HARRIS", & 575 task_list_external=ec_env%task_list) 576 END IF 577 578 ! calclulate Int(vhxc*rho)dr and Int(vtau*tau)dr 579 evhxc = evhxc + pw_integral_ab(rho_r(ispin)%pw, v_rspace(ispin)%pw)/v_rspace(1)%pw%pw_grid%dvol 580 IF (ASSOCIATED(v_tau_rspace)) THEN 581 evhxc = evhxc + pw_integral_ab(tau_r(ispin)%pw, v_tau_rspace(ispin)%pw)/ & 582 v_tau_rspace(ispin)%pw%pw_grid%dvol 583 END IF 584 585 END DO 586 587 ! return pw grids 588 DO ispin = 1, nspins 589 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace(ispin)%pw) 590 IF (ASSOCIATED(v_tau_rspace)) THEN 591 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw) 592 END IF 593 ENDDO 594 DEALLOCATE (v_rspace) 595 596 ! energies 597 ec_env%exc = eexc 598 ec_env%vhxc = evhxc 599 600 ! add the core matrix 601 DO ispin = 1, nspins 602 CALL dbcsr_add(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_h(1, 1)%matrix, & 603 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 604 CALL dbcsr_filter(ec_env%matrix_ks(ispin, 1)%matrix, & 605 dft_control%qs_control%eps_filter_matrix) 606 END DO 607 608 CALL timestop(handle) 609 610 END SUBROUTINE ec_build_ks_matrix 611 612! ************************************************************************************************** 613!> \brief Construction of the Core Hamiltonian Matrix 614!> Short version of qs_core_hamiltonian 615!> \param qs_env ... 616!> \param ec_env ... 617!> \author Creation (03.2014,JGH) 618! ************************************************************************************************** 619 SUBROUTINE ec_build_core_hamiltonian_force(qs_env, ec_env) 620 TYPE(qs_environment_type), POINTER :: qs_env 621 TYPE(energy_correction_type), POINTER :: ec_env 622 623 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_core_hamiltonian_force', & 624 routineP = moduleN//':'//routineN 625 626 INTEGER :: handle, iounit, nder, nimages 627 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 628 LOGICAL :: calculate_forces, use_virial 629 REAL(KIND=dp) :: eps_ppnl 630 REAL(KIND=dp), DIMENSION(3) :: fodeb 631 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 632 TYPE(cp_logger_type), POINTER :: logger 633 TYPE(cp_para_env_type), POINTER :: para_env 634 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm 635 TYPE(dft_control_type), POINTER :: dft_control 636 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 637 POINTER :: sab_orb, sac_ppl, sap_ppnl 638 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 639 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 640 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 641 TYPE(qs_ks_env_type), POINTER :: ks_env 642 TYPE(virial_type), POINTER :: virial 643 644 CALL timeset(routineN, handle) 645 646 logger => cp_get_default_logger() 647 iounit = cp_logger_get_default_unit_nr(logger) 648 649 calculate_forces = .TRUE. 650 651 ! no k-points possible 652 CALL get_qs_env(qs_env=qs_env, para_env=para_env, dft_control=dft_control) 653 nimages = dft_control%nimages 654 IF (nimages /= 1) THEN 655 CPABORT("K-points for Harris functional not implemented") 656 END IF 657 658 ! check for virial (currently no stress tensor available) 659 CALL get_qs_env(qs_env=qs_env, virial=virial) 660 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 661 IF (use_virial) THEN 662 CPABORT("Stress Tensor for Harris functional not implemented") 663 END IF 664 665 ! check for GAPW/GAPW_XC 666 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN 667 CPABORT("Harris functional for GAPW not implemented") 668 END IF 669 670 ! get neighbor lists, we need the full sab_orb list from the ec_env 671 NULLIFY (sab_orb, sac_ppl, sap_ppnl) 672 sab_orb => ec_env%sab_orb 673 sac_ppl => ec_env%sac_ppl 674 sap_ppnl => ec_env%sap_ppnl 675 676 ! initialize src matrix 677 NULLIFY (scrm) 678 CALL dbcsr_allocate_matrix_set(scrm, 1, 1) 679 ALLOCATE (scrm(1, 1)%matrix) 680 CALL dbcsr_create(scrm(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix) 681 CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb) 682 683 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, force=force) 684 nder = 1 685 IF (SIZE(ec_env%matrix_p, 1) == 2) THEN 686 CALL dbcsr_add(ec_env%matrix_p(1, 1)%matrix, ec_env%matrix_p(2, 1)%matrix, & 687 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 688 CALL dbcsr_add(ec_env%matrix_w(1, 1)%matrix, ec_env%matrix_w(2, 1)%matrix, & 689 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 690 END IF 691 ! Overlap and kinetic energy matrices 692 IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1) 693 CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, & 694 matrix_name="OVERLAP MATRIX", & 695 basis_type_a="HARRIS", & 696 basis_type_b="HARRIS", & 697 sab_nl=sab_orb, calculate_forces=.TRUE., & 698 matrixkp_p=ec_env%matrix_w) 699 IF (debug_forces) THEN 700 fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3) 701 CALL mp_sum(fodeb, para_env%group) 702 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wout*dS ", fodeb 703 fodeb(1:3) = force(1)%kinetic(1:3, 1) 704 END IF 705 CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, & 706 matrix_name="KINETIC ENERGY MATRIX", & 707 basis_type="HARRIS", & 708 sab_nl=sab_orb, calculate_forces=.TRUE., & 709 matrixkp_p=ec_env%matrix_p) 710 IF (debug_forces) THEN 711 fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3) 712 CALL mp_sum(fodeb, para_env%group) 713 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dT ", fodeb 714 END IF 715 IF (SIZE(ec_env%matrix_p, 1) == 2) THEN 716 CALL dbcsr_add(ec_env%matrix_p(1, 1)%matrix, ec_env%matrix_p(2, 1)%matrix, & 717 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) 718 END IF 719 720 ! compute the ppl contribution to the core hamiltonian 721 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, & 722 atomic_kind_set=atomic_kind_set) 723 NULLIFY (cell_to_index, virial) 724 use_virial = .FALSE. 725 IF (ASSOCIATED(sac_ppl)) THEN 726 IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1) 727 CALL build_core_ppl(scrm, ec_env%matrix_p, force, & 728 virial, calculate_forces, use_virial, nder, & 729 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & 730 nimages, cell_to_index, "HARRIS") 731 IF (calculate_forces .AND. debug_forces) THEN 732 fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3) 733 CALL mp_sum(fodeb, para_env%group) 734 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPL ", fodeb 735 END IF 736 END IF 737 738 ! compute the ppnl contribution to the core hamiltonian *** 739 eps_ppnl = dft_control%qs_control%eps_ppnl 740 IF (ASSOCIATED(sap_ppnl)) THEN 741 IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) 742 CALL build_core_ppnl(scrm, ec_env%matrix_p, force, & 743 virial, calculate_forces, use_virial, nder, & 744 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, & 745 nimages, cell_to_index, "HARRIS") 746 IF (calculate_forces .AND. debug_forces) THEN 747 fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3) 748 CALL mp_sum(fodeb, para_env%group) 749 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPNL", fodeb 750 END IF 751 END IF 752 753 ! External field (nonperiodic case) 754 ec_env%efield_nuclear = 0.0_dp 755 IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%efield(1:3, 1) 756 CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces) 757 IF (calculate_forces .AND. debug_forces) THEN 758 fodeb(1:3) = force(1)%efield(1:3, 1) - fodeb(1:3) 759 CALL mp_sum(fodeb, para_env%group) 760 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dEfield", fodeb 761 END IF 762 763 ! delete scr matrix 764 CALL dbcsr_deallocate_matrix_set(scrm) 765 766 CALL timestop(handle) 767 768 END SUBROUTINE ec_build_core_hamiltonian_force 769 770! ************************************************************************************************** 771!> \brief Solve KS equation for a given matrix 772!> \brief calculate the complete KS matrix 773!> \param qs_env ... 774!> \param ec_env ... 775!> \par History 776!> 03.2014 adapted from qs_ks_build_kohn_sham_matrix [JGH] 777!> \author JGH 778! ************************************************************************************************** 779 SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env) 780 TYPE(qs_environment_type), POINTER :: qs_env 781 TYPE(energy_correction_type), POINTER :: ec_env 782 783 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_ks_matrix_force', & 784 routineP = moduleN//':'//routineN 785 786 INTEGER :: handle, i, iounit, ispin, nao, natom, & 787 norb, nspins 788 INTEGER, DIMENSION(2, 3) :: bo 789 LOGICAL :: lsd, use_virial 790 REAL(dp) :: dehartree, eexc, eovrl, focc, & 791 total_rhoout 792 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftot 793 REAL(dp), DIMENSION(3) :: fodeb 794 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 795 TYPE(cp_blacs_env_type), POINTER :: blacs_env 796 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: vmos 797 TYPE(cp_fm_struct_type), POINTER :: fm_struct 798 TYPE(cp_fm_type), POINTER :: mo_coeff 799 TYPE(cp_logger_type), POINTER :: logger 800 TYPE(cp_para_env_type), POINTER :: para_env 801 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s, matrix_v, scrm 802 TYPE(dft_control_type), POINTER :: dft_control 803 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 804 TYPE(pw_env_type), POINTER :: pw_env 805 TYPE(pw_p_type) :: dv_hartree_rspace, rho_tot_gspace, & 806 v_hartree_gspace, v_hartree_rspace, & 807 vtot_rspace 808 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r, rhoout_g, rhoout_r, tau_r, & 809 v_rspace, v_tau_rspace, v_xc 810 TYPE(pw_poisson_type), POINTER :: poisson_env 811 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 812 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 813 TYPE(qs_ks_env_type), POINTER :: ks_env 814 TYPE(qs_rho_type), POINTER :: rho 815 TYPE(section_vals_type), POINTER :: input, xc_fun_section, xc_section 816 TYPE(virial_type), POINTER :: virial 817 TYPE(xc_derivative_set_type), POINTER :: deriv_set 818 TYPE(xc_rho_cflags_type) :: needs 819 TYPE(xc_rho_set_type), POINTER :: rho1_set, rho_set 820 821 CALL timeset(routineN, handle) 822 823 logger => cp_get_default_logger() 824 iounit = cp_logger_get_default_unit_nr(logger) 825 826 ! get all information on the electronic density 827 NULLIFY (rho, ks_env) 828 CALL get_qs_env(qs_env=qs_env, rho=rho, virial=virial, dft_control=dft_control, & 829 para_env=para_env, blacs_env=blacs_env, ks_env=ks_env) 830 831 nspins = dft_control%nspins 832 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 833 CPASSERT(.NOT. use_virial) 834 835 NULLIFY (pw_env) 836 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 837 CPASSERT(ASSOCIATED(pw_env)) 838 839 NULLIFY (auxbas_pw_pool, poisson_env) 840 ! gets the tmp grids 841 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 842 poisson_env=poisson_env) 843 844 ! Calculate the Hartree potential 845 NULLIFY (v_hartree_gspace%pw, rho_tot_gspace%pw, v_hartree_rspace%pw) 846 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 847 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 848 CALL pw_pool_create_pw(auxbas_pw_pool, rho_tot_gspace%pw, & 849 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 850 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_rspace%pw, & 851 use_data=REALDATA3D, in_space=REALSPACE) 852 853 ! Get the total density in g-space [ions + electrons] 854 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho) 855 856 CALL pw_transfer(ec_env%vh_rspace%pw, v_hartree_rspace%pw) 857 858 CALL get_qs_env(qs_env=qs_env, force=force) 859 ! calculate output density on grid 860 ! rho_in(R): CALL qs_rho_get(rho, rho_r=rho_r) 861 ! rho_in(G): CALL qs_rho_get(rho, rho_g=rho_g) 862 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g, tau_r=tau_r) 863 ALLOCATE (rhoout_r(nspins), rhoout_g(nspins)) 864 DO ispin = 1, nspins 865 NULLIFY (rhoout_r(ispin)%pw, rhoout_g(ispin)%pw) 866 CALL pw_pool_create_pw(auxbas_pw_pool, rhoout_r(ispin)%pw, & 867 use_data=REALDATA3D, in_space=REALSPACE) 868 CALL pw_pool_create_pw(auxbas_pw_pool, rhoout_g(ispin)%pw, & 869 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 870 ENDDO 871 CALL pw_pool_create_pw(auxbas_pw_pool, dv_hartree_rspace%pw, & 872 use_data=REALDATA3D, in_space=REALSPACE) 873 CALL pw_pool_create_pw(auxbas_pw_pool, vtot_rspace%pw, & 874 use_data=REALDATA3D, in_space=REALSPACE) 875 876 CALL pw_zero(rho_tot_gspace%pw) 877 DO ispin = 1, nspins 878 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=ec_env%matrix_p(ispin, 1)%matrix, & 879 rho=rhoout_r(ispin), & 880 rho_gspace=rhoout_g(ispin), & 881 total_rho=total_rhoout, & 882 basis_type="HARRIS", & 883 task_list_external=ec_env%task_list) 884 ! rho_out - rho_in 885 CALL pw_axpy(rho_r(ispin)%pw, rhoout_r(ispin)%pw, -1.0_dp) 886 CALL pw_axpy(rho_g(ispin)%pw, rhoout_g(ispin)%pw, -1.0_dp) 887 CALL pw_axpy(rhoout_g(ispin)%pw, rho_tot_gspace%pw) 888 END DO 889 ! calculate associated hartree potential 890 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, dehartree, & 891 v_hartree_gspace%pw) 892 CALL pw_transfer(v_hartree_gspace%pw, dv_hartree_rspace%pw) 893 CALL pw_scale(dv_hartree_rspace%pw, dv_hartree_rspace%pw%pw_grid%dvol) 894 ! Getting nuclear force contribution from the core charge density 895 ! Vh(rho_in + rho_c) + Vh(rho_out - rho_in) 896 CALL pw_transfer(v_hartree_rspace%pw, vtot_rspace%pw) 897 CALL pw_axpy(dv_hartree_rspace%pw, vtot_rspace%pw) 898 IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1) 899 CALL integrate_v_core_rspace(vtot_rspace, qs_env) 900 IF (debug_forces) THEN 901 fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3) 902 CALL mp_sum(fodeb, para_env%group) 903 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vtot*dncore", fodeb 904 END IF 905 ! 906 ! Pulay force from Tr P_in (V_H(drho)+ Fxc(rho_in)*drho) 907 ! RHS of CPKS equations: (V_H(drho)+ Fxc(rho_in)*drho)*C0 908 ! Fxc*drho term 909 ALLOCATE (v_xc(nspins)) 910 DO ispin = 1, nspins 911 NULLIFY (v_xc(ispin)%pw) 912 CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, & 913 use_data=REALDATA3D, in_space=REALSPACE) 914 CALL pw_zero(v_xc(ispin)%pw) 915 END DO 916 CALL get_qs_env(qs_env, input=input) 917 xc_section => ec_env%xc_section 918 lsd = (nspins == 2) 919 NULLIFY (deriv_set, rho_set, rho1_set) 920 CALL xc_prep_2nd_deriv(deriv_set, rho_set, rho_r, auxbas_pw_pool, xc_section=xc_section) 921 bo = rhoout_r(1)%pw%pw_grid%bounds_local 922 CALL xc_rho_set_create(rho1_set, bo, & 923 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), & 924 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), & 925 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF")) 926 927 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL") 928 needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.) 929 930 ! calculate the arguments needed by the functionals 931 CALL xc_rho_set_update(rho1_set, rhoout_r, rhoout_g, tau_r, needs, & 932 section_get_ival(xc_section, "XC_GRID%XC_DERIV"), & 933 section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), & 934 auxbas_pw_pool) 935 CALL xc_calc_2nd_deriv(v_xc, deriv_set, rho_set, & 936 rho1_set, auxbas_pw_pool, xc_section=xc_section) 937 CALL xc_dset_release(deriv_set) 938 CALL xc_rho_set_release(rho_set) 939 CALL xc_rho_set_release(rho1_set) 940 ! 941 CALL get_qs_env(qs_env=qs_env, rho=rho, matrix_s_kp=matrix_s) 942 NULLIFY (matrix_v) 943 CALL dbcsr_allocate_matrix_set(matrix_v, nspins, 1) 944 DO ispin = 1, nspins 945 ALLOCATE (matrix_v(ispin, 1)%matrix) 946 CALL dbcsr_create(matrix_v(ispin, 1)%matrix, template=matrix_s(1, 1)%matrix) 947 CALL dbcsr_copy(matrix_v(ispin, 1)%matrix, matrix_s(1, 1)%matrix) 948 CALL dbcsr_set(matrix_v(ispin, 1)%matrix, 0.0_dp) 949 END DO 950 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 951 ! vtot = v_xc(ispin) + dv_hartree 952 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) 953 DO ispin = 1, nspins 954 CALL pw_scale(v_xc(ispin)%pw, v_xc(ispin)%pw%pw_grid%dvol) 955 CALL pw_axpy(dv_hartree_rspace%pw, v_xc(ispin)%pw) 956 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), & 957 hmat=matrix_v(ispin, 1), & 958 pmat=matrix_p(ispin, 1), & 959 calculate_forces=.TRUE.) 960 END DO 961 IF (debug_forces) THEN 962 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3) 963 CALL mp_sum(fodeb, para_env%group) 964 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKdrho", fodeb 965 END IF 966 ! CPKS vector 967 ! 968 ! sign is changed in linres_solver! 969 ! Projector Q applied in linres_solver! 970 focc = 2.0_dp 971 IF (nspins == 1) focc = 4.0_dp 972 CALL get_qs_env(qs_env=qs_env, mos=mos) 973 ALLOCATE (vmos(nspins)) 974 DO ispin = 1, nspins 975 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 976 CALL cp_fm_get_info(mo_coeff, matrix_struct=fm_struct, nrow_global=nao, ncol_global=norb) 977 NULLIFY (vmos(ispin)%matrix) 978 CALL cp_fm_create(vmos(ispin)%matrix, fm_struct) 979 CALL cp_dbcsr_sm_fm_multiply(matrix_v(ispin, 1)%matrix, mo_coeff, vmos(ispin)%matrix, & 980 norb, alpha=focc, beta=0.0_dp) 981 END DO 982 !MK The Intel compiler has trouble with that pointer assignment 983 !MK ec_env%matrix_hz => matrix_v(:, 1) 984 !MK Begin explicit copy 985 CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins) 986 DO ispin = 1, nspins 987 ALLOCATE (ec_env%matrix_hz(ispin)%matrix) 988 CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_v(ispin, 1)%matrix) 989 CALL dbcsr_copy(matrix_v(ispin, 1)%matrix, ec_env%matrix_hz(ispin)%matrix) 990 END DO 991 CALL dbcsr_deallocate_matrix_set(matrix_v) 992 !MK End 993 IF (ASSOCIATED(ec_env%cpmos)) THEN 994 DO i = 1, SIZE(ec_env%cpmos) 995 CALL cp_fm_release(ec_env%cpmos(i)%matrix) 996 END DO 997 DEALLOCATE (ec_env%cpmos) 998 NULLIFY (ec_env%cpmos) 999 END IF 1000 ec_env%cpmos => vmos 1001 ! 1002 CALL pw_pool_give_back_pw(auxbas_pw_pool, dv_hartree_rspace%pw) 1003 CALL pw_pool_give_back_pw(auxbas_pw_pool, vtot_rspace%pw) 1004 DO ispin = 1, nspins 1005 CALL pw_pool_give_back_pw(auxbas_pw_pool, rhoout_r(ispin)%pw) 1006 CALL pw_pool_give_back_pw(auxbas_pw_pool, rhoout_g(ispin)%pw) 1007 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_xc(ispin)%pw) 1008 END DO 1009 DEALLOCATE (rhoout_r, rhoout_g, v_xc) 1010 1011 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw) 1012 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw) 1013 1014 ! initialize src matrix 1015 NULLIFY (scrm) 1016 CALL dbcsr_allocate_matrix_set(scrm, 1, 1) 1017 ALLOCATE (scrm(1, 1)%matrix) 1018 CALL dbcsr_create(scrm(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix) 1019 CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, ec_env%sab_orb) 1020 1021 ! v_rspace and v_tau_rspace are generated from the auxbas pool 1022 NULLIFY (v_rspace, v_tau_rspace) 1023 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, & 1024 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.FALSE.) 1025 IF (.NOT. ASSOCIATED(v_rspace)) THEN 1026 ALLOCATE (v_rspace(nspins)) 1027 DO ispin = 1, nspins 1028 NULLIFY (v_rspace(ispin)%pw) 1029 CALL pw_pool_create_pw(auxbas_pw_pool, v_rspace(ispin)%pw, & 1030 use_data=REALDATA3D, in_space=REALSPACE) 1031 CALL pw_zero(v_rspace(ispin)%pw) 1032 ENDDO 1033 END IF 1034 1035 CALL qs_rho_get(rho, rho_r=rho_r) 1036 IF (ASSOCIATED(v_tau_rspace)) THEN 1037 CALL qs_rho_get(rho, tau_r=tau_r) 1038 END IF 1039 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) 1040 DO ispin = 1, nspins 1041 ! Add v_hartree + v_xc = v_rspace 1042 CALL pw_scale(v_rspace(ispin)%pw, v_rspace(ispin)%pw%pw_grid%dvol) 1043 CALL pw_axpy(v_hartree_rspace%pw, v_rspace(ispin)%pw) 1044 ! integrate over potential <a|V|b> 1045 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), & 1046 hmat=scrm(1, 1), & 1047 pmat=ec_env%matrix_p(ispin, 1), & 1048 qs_env=qs_env, & 1049 calculate_forces=.TRUE., & 1050 basis_type="HARRIS", & 1051 task_list_external=ec_env%task_list) 1052 1053 IF (ASSOCIATED(v_tau_rspace)) THEN 1054 ! integrate over Tau-potential <nabla.a|V|nabla.b> 1055 CALL pw_scale(v_tau_rspace(ispin)%pw, v_tau_rspace(ispin)%pw%pw_grid%dvol) 1056 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), hmat=scrm(1, 1), & 1057 pmat=ec_env%matrix_p(ispin, 1), & 1058 qs_env=qs_env, calculate_forces=.TRUE., compute_tau=.TRUE., & 1059 basis_type="HARRIS", & 1060 task_list_external=ec_env%task_list) 1061 END IF 1062 1063 END DO 1064 IF (debug_forces) THEN 1065 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3) 1066 CALL mp_sum(fodeb, para_env%group) 1067 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc ", fodeb 1068 END IF 1069 1070 ! delete scr matrix 1071 CALL dbcsr_deallocate_matrix_set(scrm) 1072 1073 ! return pw grids 1074 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw) 1075 DO ispin = 1, nspins 1076 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace(ispin)%pw) 1077 IF (ASSOCIATED(v_tau_rspace)) THEN 1078 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw) 1079 END IF 1080 ENDDO 1081 1082 ! Core overlap 1083 IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1) 1084 CALL calculate_ecore_overlap(qs_env, para_env, .TRUE., E_overlap_core=eovrl) 1085 IF (debug_forces) THEN 1086 fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3) 1087 CALL mp_sum(fodeb, para_env%group) 1088 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: CoreOverlap", fodeb 1089 END IF 1090 1091 IF (debug_forces) THEN 1092 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set) 1093 ALLOCATE (ftot(3, natom)) 1094 CALL total_qs_force(ftot, force, atomic_kind_set) 1095 fodeb(1:3) = ftot(1:3, 1) 1096 DEALLOCATE (ftot) 1097 CALL mp_sum(fodeb, para_env%group) 1098 IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Force Explicit", fodeb 1099 END IF 1100 1101 DEALLOCATE (v_rspace) 1102 1103 CALL timestop(handle) 1104 1105 END SUBROUTINE ec_build_ks_matrix_force 1106 1107! ************************************************************************************************** 1108!> \brief Solve KS equation for a given matrix 1109!> \param qs_env ... 1110!> \param ec_env ... 1111!> \par History 1112!> 03.2014 created [JGH] 1113!> \author JGH 1114! ************************************************************************************************** 1115 1116 SUBROUTINE ec_ks_solver(qs_env, ec_env) 1117 1118 TYPE(qs_environment_type), POINTER :: qs_env 1119 TYPE(energy_correction_type), POINTER :: ec_env 1120 1121 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ks_solver', routineP = moduleN//':'//routineN 1122 1123 CHARACTER(LEN=default_string_length) :: headline 1124 INTEGER :: handle, ispin, nspins 1125 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, pmat, smat, wmat 1126 TYPE(dft_control_type), POINTER :: dft_control 1127 1128 CALL timeset(routineN, handle) 1129 1130 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 1131 nspins = dft_control%nspins 1132 1133 ! create density matrix 1134 IF (.NOT. ASSOCIATED(ec_env%matrix_p)) THEN 1135 headline = "DENSITY MATRIX" 1136 CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1) 1137 DO ispin = 1, nspins 1138 ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix) 1139 CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=TRIM(headline), & 1140 template=ec_env%matrix_s(1, 1)%matrix) 1141 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, ec_env%sab_orb) 1142 END DO 1143 END IF 1144 ! create energy weighted density matrix 1145 IF (.NOT. ASSOCIATED(ec_env%matrix_w)) THEN 1146 headline = "ENERGY WEIGHTED DENSITY MATRIX" 1147 CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1) 1148 DO ispin = 1, nspins 1149 ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix) 1150 CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=TRIM(headline), & 1151 template=ec_env%matrix_s(1, 1)%matrix) 1152 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, ec_env%sab_orb) 1153 END DO 1154 END IF 1155 1156 IF (ec_env%mao) THEN 1157 CALL mao_create_matrices(ec_env, ksmat, smat, pmat, wmat) 1158 ELSE 1159 ksmat => ec_env%matrix_ks 1160 smat => ec_env%matrix_s 1161 pmat => ec_env%matrix_p 1162 wmat => ec_env%matrix_w 1163 ENDIF 1164 1165 SELECT CASE (ec_env%ks_solver) 1166 CASE (ec_diagonalization) 1167 CALL ec_diag_solver(qs_env, ksmat, smat, pmat, wmat) 1168 CASE (ec_curvy_steps, ec_matrix_sign, ec_matrix_trs4, ec_matrix_tc2) 1169 CALL ec_ls_solver(qs_env, ec_env%ks_solver, ksmat, smat, pmat, wmat) 1170 CASE DEFAULT 1171 CPASSERT(.FALSE.) 1172 END SELECT 1173 1174 IF (ec_env%mao) THEN 1175 CALL mao_release_matrices(ec_env, ksmat, smat, pmat, wmat) 1176 ENDIF 1177 1178 CALL timestop(handle) 1179 1180 END SUBROUTINE ec_ks_solver 1181 1182! ************************************************************************************************** 1183!> \brief Create matrices with MAO sizes 1184!> \param ec_env ... 1185!> \param ksmat ... 1186!> \param smat ... 1187!> \param pmat ... 1188!> \param wmat ... 1189!> \par History 1190!> 08.2016 created [JGH] 1191!> \author JGH 1192! ************************************************************************************************** 1193 1194 SUBROUTINE mao_create_matrices(ec_env, ksmat, smat, pmat, wmat) 1195 1196 TYPE(energy_correction_type), POINTER :: ec_env 1197 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, smat, pmat, wmat 1198 1199 CHARACTER(LEN=*), PARAMETER :: routineN = 'mao_create_matrices', & 1200 routineP = moduleN//':'//routineN 1201 1202 INTEGER :: handle, ispin, nspins 1203 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes 1204 TYPE(dbcsr_distribution_type) :: dbcsr_dist 1205 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef 1206 TYPE(dbcsr_type) :: cgmat 1207 1208 CALL timeset(routineN, handle) 1209 1210 mao_coef => ec_env%mao_coef 1211 1212 NULLIFY (ksmat, smat, pmat, wmat) 1213 nspins = SIZE(ec_env%matrix_ks, 1) 1214 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist) 1215 CALL dbcsr_allocate_matrix_set(ksmat, nspins, 1) 1216 CALL dbcsr_allocate_matrix_set(smat, nspins, 1) 1217 DO ispin = 1, nspins 1218 ALLOCATE (ksmat(ispin, 1)%matrix) 1219 CALL dbcsr_create(ksmat(ispin, 1)%matrix, dist=dbcsr_dist, name="MAO KS mat", & 1220 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, & 1221 col_blk_size=col_blk_sizes, nze=0) 1222 ALLOCATE (smat(ispin, 1)%matrix) 1223 CALL dbcsr_create(smat(ispin, 1)%matrix, dist=dbcsr_dist, name="MAO S mat", & 1224 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, & 1225 col_blk_size=col_blk_sizes, nze=0) 1226 END DO 1227 ! 1228 CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix) 1229 DO ispin = 1, nspins 1230 CALL dbcsr_multiply("N", "N", 1.0_dp, ec_env%matrix_s(1, 1)%matrix, mao_coef(ispin)%matrix, & 1231 0.0_dp, cgmat) 1232 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, smat(ispin, 1)%matrix) 1233 CALL dbcsr_multiply("N", "N", 1.0_dp, ec_env%matrix_ks(1, 1)%matrix, mao_coef(ispin)%matrix, & 1234 0.0_dp, cgmat) 1235 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, ksmat(ispin, 1)%matrix) 1236 END DO 1237 CALL dbcsr_release(cgmat) 1238 1239 CALL dbcsr_allocate_matrix_set(pmat, nspins, 1) 1240 DO ispin = 1, nspins 1241 ALLOCATE (pmat(ispin, 1)%matrix) 1242 CALL dbcsr_create(pmat(ispin, 1)%matrix, template=smat(1, 1)%matrix) 1243 CALL cp_dbcsr_alloc_block_from_nbl(pmat(ispin, 1)%matrix, ec_env%sab_orb) 1244 END DO 1245 1246 CALL dbcsr_allocate_matrix_set(wmat, nspins, 1) 1247 DO ispin = 1, nspins 1248 ALLOCATE (wmat(ispin, 1)%matrix) 1249 CALL dbcsr_create(wmat(ispin, 1)%matrix, template=smat(1, 1)%matrix) 1250 CALL cp_dbcsr_alloc_block_from_nbl(wmat(ispin, 1)%matrix, ec_env%sab_orb) 1251 END DO 1252 1253 CALL timestop(handle) 1254 1255 END SUBROUTINE mao_create_matrices 1256 1257! ************************************************************************************************** 1258!> \brief Release matrices with MAO sizes 1259!> \param ec_env ... 1260!> \param ksmat ... 1261!> \param smat ... 1262!> \param pmat ... 1263!> \param wmat ... 1264!> \par History 1265!> 08.2016 created [JGH] 1266!> \author JGH 1267! ************************************************************************************************** 1268 1269 SUBROUTINE mao_release_matrices(ec_env, ksmat, smat, pmat, wmat) 1270 1271 TYPE(energy_correction_type), POINTER :: ec_env 1272 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, smat, pmat, wmat 1273 1274 CHARACTER(LEN=*), PARAMETER :: routineN = 'mao_release_matrices', & 1275 routineP = moduleN//':'//routineN 1276 1277 INTEGER :: handle, ispin, nspins 1278 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef 1279 TYPE(dbcsr_type) :: cgmat 1280 1281 CALL timeset(routineN, handle) 1282 1283 mao_coef => ec_env%mao_coef 1284 nspins = SIZE(mao_coef, 1) 1285 1286 ! save pmat/wmat in full basis format 1287 CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix) 1288 DO ispin = 1, nspins 1289 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef(ispin)%matrix, pmat(ispin, 1)%matrix, 0.0_dp, cgmat) 1290 CALL dbcsr_multiply("N", "T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, & 1291 ec_env%matrix_p(ispin, 1)%matrix, retain_sparsity=.TRUE.) 1292 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef(ispin)%matrix, wmat(ispin, 1)%matrix, 0.0_dp, cgmat) 1293 CALL dbcsr_multiply("N", "T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, & 1294 ec_env%matrix_w(ispin, 1)%matrix, retain_sparsity=.TRUE.) 1295 END DO 1296 CALL dbcsr_release(cgmat) 1297 1298 CALL dbcsr_deallocate_matrix_set(ksmat) 1299 CALL dbcsr_deallocate_matrix_set(smat) 1300 CALL dbcsr_deallocate_matrix_set(pmat) 1301 CALL dbcsr_deallocate_matrix_set(wmat) 1302 1303 CALL timestop(handle) 1304 1305 END SUBROUTINE mao_release_matrices 1306 1307! ************************************************************************************************** 1308!> \brief Solve KS equation using diagonalization 1309!> \param qs_env ... 1310!> \param matrix_ks ... 1311!> \param matrix_s ... 1312!> \param matrix_p ... 1313!> \param matrix_w ... 1314!> \par History 1315!> 03.2014 created [JGH] 1316!> \author JGH 1317! ************************************************************************************************** 1318 1319 SUBROUTINE ec_diag_solver(qs_env, matrix_ks, matrix_s, matrix_p, matrix_w) 1320 1321 TYPE(qs_environment_type), POINTER :: qs_env 1322 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_p, matrix_w 1323 1324 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_diag_solver', routineP = moduleN//':'//routineN 1325 1326 INTEGER :: handle, info, ispin, nmo(2), nsize, & 1327 nspins 1328 REAL(KIND=dp) :: eps_filter, focc(2) 1329 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues 1330 TYPE(cp_blacs_env_type), POINTER :: blacs_env 1331 TYPE(cp_fm_struct_type), POINTER :: fm_struct 1332 TYPE(cp_fm_type), POINTER :: fm_ks, fm_mo, fm_ortho 1333 TYPE(cp_para_env_type), POINTER :: para_env 1334 TYPE(dbcsr_type), POINTER :: buf1_dbcsr, buf2_dbcsr, buf3_dbcsr, & 1335 ortho_dbcsr, ref_matrix 1336 TYPE(dft_control_type), POINTER :: dft_control 1337 1338 CALL timeset(routineN, handle) 1339 1340 NULLIFY (blacs_env, para_env) 1341 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, para_env=para_env) 1342 1343 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 1344 eps_filter = dft_control%qs_control%eps_filter_matrix 1345 nspins = dft_control%nspins 1346 1347 nmo = 0 1348 CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo) 1349 focc = 1._dp 1350 IF (nspins == 1) THEN 1351 focc = 2._dp 1352 nmo(1) = nmo(1)/2 1353 END IF 1354 1355 CALL dbcsr_get_info(matrix_ks(1, 1)%matrix, nfullrows_total=nsize) 1356 ALLOCATE (eigenvalues(nsize)) 1357 1358 NULLIFY (fm_ortho, fm_ks, fm_mo, fm_struct, ref_matrix) 1359 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nsize, & 1360 ncol_global=nsize, para_env=para_env) 1361 CALL cp_fm_create(fm_ortho, fm_struct) 1362 CALL cp_fm_create(fm_ks, fm_struct) 1363 CALL cp_fm_create(fm_mo, fm_struct) 1364 CALL cp_fm_struct_release(fm_struct) 1365 1366 ! factorization 1367 ref_matrix => matrix_s(1, 1)%matrix 1368 NULLIFY (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr) 1369 CALL dbcsr_init_p(ortho_dbcsr) 1370 CALL dbcsr_create(ortho_dbcsr, template=ref_matrix, & 1371 matrix_type=dbcsr_type_no_symmetry) 1372 CALL dbcsr_init_p(buf1_dbcsr) 1373 CALL dbcsr_create(buf1_dbcsr, template=ref_matrix, & 1374 matrix_type=dbcsr_type_no_symmetry) 1375 CALL dbcsr_init_p(buf2_dbcsr) 1376 CALL dbcsr_create(buf2_dbcsr, template=ref_matrix, & 1377 matrix_type=dbcsr_type_no_symmetry) 1378 CALL dbcsr_init_p(buf3_dbcsr) 1379 CALL dbcsr_create(buf3_dbcsr, template=ref_matrix, & 1380 matrix_type=dbcsr_type_no_symmetry) 1381 1382 ref_matrix => matrix_s(1, 1)%matrix 1383 CALL copy_dbcsr_to_fm(ref_matrix, fm_ortho) 1384 CALL cp_fm_cholesky_decompose(fm_ortho) 1385 CALL cp_fm_triangular_invert(fm_ortho) 1386 CALL cp_fm_set_all(fm_ks, 0.0_dp) 1387 CALL cp_fm_to_fm_triangular(fm_ortho, fm_ks, "U") 1388 CALL copy_fm_to_dbcsr(fm_ks, ortho_dbcsr) 1389 DO ispin = 1, nspins 1390 ! calculate ZHZ(T) 1391 CALL dbcsr_desymmetrize(matrix_ks(ispin, 1)%matrix, buf1_dbcsr) 1392 CALL dbcsr_multiply("N", "N", 1.0_dp, buf1_dbcsr, ortho_dbcsr, & 1393 0.0_dp, buf2_dbcsr, filter_eps=eps_filter) 1394 CALL dbcsr_multiply("T", "N", 1.0_dp, ortho_dbcsr, buf2_dbcsr, & 1395 0.0_dp, buf1_dbcsr, filter_eps=eps_filter) 1396 ! copy to fm format 1397 CALL copy_dbcsr_to_fm(buf1_dbcsr, fm_ks) 1398 CALL choose_eigv_solver(fm_ks, fm_mo, eigenvalues, info) 1399 CPASSERT(info == 0) 1400 ! back transform of mos c = Z(T)*c 1401 CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr) 1402 CALL dbcsr_multiply("N", "N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, & 1403 0.0_dp, buf2_dbcsr, filter_eps=eps_filter) 1404 ! density matrix 1405 CALL dbcsr_set(matrix_p(ispin, 1)%matrix, 0.0_dp) 1406 CALL dbcsr_multiply("N", "T", focc(ispin), buf2_dbcsr, buf2_dbcsr, & 1407 1.0_dp, matrix_p(ispin, 1)%matrix, retain_sparsity=.TRUE., last_k=nmo(ispin)) 1408 ! energy weighted density matrix 1409 CALL dbcsr_set(matrix_w(ispin, 1)%matrix, 0.0_dp) 1410 CALL cp_fm_column_scale(fm_mo, eigenvalues) 1411 CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr) 1412 CALL dbcsr_multiply("N", "N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, & 1413 0.0_dp, buf3_dbcsr, filter_eps=eps_filter) 1414 CALL dbcsr_multiply("N", "T", focc(ispin), buf2_dbcsr, buf3_dbcsr, & 1415 1.0_dp, matrix_w(ispin, 1)%matrix, retain_sparsity=.TRUE., last_k=nmo(ispin)) 1416 END DO 1417 1418 CALL cp_fm_release(fm_ks) 1419 CALL cp_fm_release(fm_mo) 1420 CALL cp_fm_release(fm_ortho) 1421 CALL dbcsr_release(ortho_dbcsr) 1422 CALL dbcsr_release(buf1_dbcsr) 1423 CALL dbcsr_release(buf2_dbcsr) 1424 CALL dbcsr_release(buf3_dbcsr) 1425 DEALLOCATE (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr) 1426 DEALLOCATE (eigenvalues) 1427 1428 CALL timestop(handle) 1429 1430 END SUBROUTINE ec_diag_solver 1431 1432! ************************************************************************************************** 1433!> \brief Solve KS equation using linear scaling methods 1434!> \param qs_env ... 1435!> \param ls_method ... 1436!> \param matrix_ks ... 1437!> \param matrix_s ... 1438!> \param matrix_p ... 1439!> \param matrix_w ... 1440!> \par History 1441!> 12.2019 created [JGH] 1442!> \author JGH 1443! ************************************************************************************************** 1444 1445 SUBROUTINE ec_ls_solver(qs_env, ls_method, matrix_ks, matrix_s, matrix_p, matrix_w) 1446 1447 TYPE(qs_environment_type), POINTER :: qs_env 1448 INTEGER, INTENT(IN) :: ls_method 1449 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_p, matrix_w 1450 1451 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ls_solver', routineP = moduleN//':'//routineN 1452 1453 INTEGER :: cluster_type, extrapolation_order, handle, ispin, max_iter_lanczos, nelectron, & 1454 nspins, s_inversion_type, s_preconditioner_type, s_sqrt_method, s_sqrt_order, & 1455 sign_method, sign_order 1456 INTEGER, DIMENSION(2) :: nmo 1457 LOGICAL :: check_s_inv, converged, & 1458 dynamic_threshold, fixed_mu, & 1459 non_monotonic, report_all_sparsities, & 1460 single_precision 1461 REAL(KIND=dp) :: eps_filter, eps_lanczos, mu 1462 REAL(KIND=dp), DIMENSION(2) :: e_homo, e_lumo, mu_spin 1463 TYPE(dbcsr_type) :: matrix_s_inv, matrix_s_sqrt, & 1464 matrix_s_sqrt_inv, matrix_tmp 1465 TYPE(dft_control_type), POINTER :: dft_control 1466 TYPE(section_vals_type), POINTER :: input, solver_section 1467 1468 CALL timeset(routineN, handle) 1469 1470 CPASSERT(ls_method > 0) 1471 CPASSERT(ASSOCIATED(qs_env)) 1472 CPASSERT(ASSOCIATED(matrix_ks)) 1473 CPASSERT(ASSOCIATED(matrix_s)) 1474 CPASSERT(ASSOCIATED(matrix_p)) 1475 CPASSERT(ASSOCIATED(matrix_w)) 1476 1477 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, input=input) 1478 nspins = dft_control%nspins 1479 ! 1480 solver_section => section_vals_get_subs_vals(input, "DFT%ENERGY_CORRECTION%LS_SOLVER") 1481 CALL section_vals_val_get(solver_section, "EPS_FILTER", r_val=eps_filter) 1482 CALL section_vals_val_get(solver_section, "MU", r_val=mu) 1483 CALL section_vals_val_get(solver_section, "FIXED_MU", l_val=fixed_mu) 1484 mu_spin = mu 1485 CALL section_vals_val_get(solver_section, "S_PRECONDITIONER", i_val=s_preconditioner_type) 1486 CALL section_vals_val_get(solver_section, "MATRIX_CLUSTER_TYPE", i_val=cluster_type) 1487 CALL section_vals_val_get(solver_section, "SINGLE_PRECISION_MATRICES", l_val=single_precision) 1488 CALL section_vals_val_get(solver_section, "S_INVERSION", i_val=s_inversion_type) 1489 CALL section_vals_val_get(solver_section, "CHECK_S_INV", l_val=check_s_inv) 1490 CALL section_vals_val_get(solver_section, "REPORT_ALL_SPARSITIES", l_val=report_all_sparsities) 1491 CALL section_vals_val_get(solver_section, "SIGN_METHOD", i_val=sign_method) 1492 CALL section_vals_val_get(solver_section, "SIGN_ORDER", i_val=sign_order) 1493 CALL section_vals_val_get(solver_section, "DYNAMIC_THRESHOLD", l_val=dynamic_threshold) 1494 CALL section_vals_val_get(solver_section, "NON_MONOTONIC", l_val=non_monotonic) 1495 CALL section_vals_val_get(solver_section, "S_SQRT_METHOD", i_val=s_sqrt_method) 1496 CALL section_vals_val_get(solver_section, "S_SQRT_ORDER", i_val=s_sqrt_order) 1497 CALL section_vals_val_get(solver_section, "EXTRAPOLATION_ORDER", i_val=extrapolation_order) 1498 CALL section_vals_val_get(solver_section, "EPS_LANCZOS", r_val=eps_lanczos) 1499 CALL section_vals_val_get(solver_section, "MAX_ITER_LANCZOS", i_val=max_iter_lanczos) 1500 ! 1501 nmo = 0 1502 CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo) 1503 IF (nspins == 1) nmo(1) = nmo(1)/2 1504 e_homo = -0.1_dp 1505 e_lumo = 0.1_dp 1506 1507 SELECT CASE (ls_method) 1508 CASE (ec_curvy_steps) 1509 CPABORT("Curvy Step Method not available") 1510 CASE (ec_matrix_sign) 1511 CALL dbcsr_create(matrix_s_inv, template=matrix_s(1, 1)%matrix, & 1512 matrix_type=dbcsr_type_no_symmetry) 1513 CALL invert_Hotelling(matrix_s_inv, matrix_s(1, 1)%matrix, eps_filter) 1514 CASE (ec_matrix_trs4, ec_matrix_tc2) 1515 CALL dbcsr_create(matrix_s_sqrt, template=matrix_s(1, 1)%matrix, & 1516 matrix_type=dbcsr_type_no_symmetry) 1517 CALL dbcsr_create(matrix_s_sqrt_inv, template=matrix_s(1, 1)%matrix, & 1518 matrix_type=dbcsr_type_no_symmetry) 1519 SELECT CASE (s_sqrt_method) 1520 CASE (ls_s_sqrt_proot) 1521 CALL matrix_sqrt_proot(matrix_s_sqrt, matrix_s_sqrt_inv, & 1522 matrix_s(1, 1)%matrix, eps_filter, & 1523 s_sqrt_order, eps_lanczos, max_iter_lanczos, symmetrize=.TRUE.) 1524 CASE (ls_s_sqrt_ns) 1525 CALL matrix_sqrt_Newton_Schulz(matrix_s_sqrt, matrix_s_sqrt_inv, & 1526 matrix_s(1, 1)%matrix, eps_filter, & 1527 s_sqrt_order, eps_lanczos, max_iter_lanczos) 1528 CASE DEFAULT 1529 CPABORT("Unknown sqrt method.") 1530 END SELECT 1531 CALL dbcsr_release(matrix_s_sqrt) 1532 CASE DEFAULT 1533 CPABORT("Wrong ls_method") 1534 END SELECT 1535 1536 DO ispin = 1, nspins 1537 nelectron = nmo(ispin) 1538 1539 SELECT CASE (ls_method) 1540 CASE (ec_curvy_steps) 1541 CPABORT("Curvy Step Method not available") 1542 CASE (ec_matrix_sign) 1543 CALL density_matrix_sign(matrix_p(ispin, 1)%matrix, mu_spin(ispin), fixed_mu, & 1544 sign_method, sign_order, matrix_ks(ispin, 1)%matrix, & 1545 matrix_s(1, 1)%matrix, matrix_s_inv, nelectron, eps_filter) 1546 CASE (ec_matrix_trs4) 1547 CALL density_matrix_trs4(matrix_p(ispin, 1)%matrix, matrix_ks(ispin, 1)%matrix, & 1548 matrix_s_sqrt_inv, nelectron, eps_filter, & 1549 e_homo(ispin), e_lumo(ispin), mu_spin(ispin), & 1550 max_iter_lanczos=max_iter_lanczos, & 1551 eps_lanczos=eps_lanczos, converged=converged) 1552 CASE (ec_matrix_tc2) 1553 CPABORT("TC2 Method not available") 1554 CALL density_matrix_tc2(matrix_p(ispin, 1)%matrix, matrix_ks(ispin, 1)%matrix, & 1555 matrix_s_sqrt_inv, nelectron, eps_filter, & 1556 e_homo(ispin), e_lumo(ispin), & 1557 non_monotonic, eps_lanczos, max_iter_lanczos) 1558 CASE DEFAULT 1559 CPABORT("Wrong ls_method") 1560 END SELECT 1561 END DO 1562 1563 SELECT CASE (ls_method) 1564 CASE (ec_curvy_steps) 1565 CPABORT("Curvy Step Method not available") 1566 CASE (ec_matrix_sign) 1567 CALL dbcsr_release(matrix_s_inv) 1568 CASE (ec_matrix_trs4, ec_matrix_tc2) 1569 CALL dbcsr_release(matrix_s_sqrt_inv) 1570 CASE DEFAULT 1571 CPABORT("Wrong ls_method") 1572 END SELECT 1573 1574 ! W matrix 1575 CALL dbcsr_create(matrix_tmp, template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry) 1576 DO ispin = 1, nspins 1577 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p(ispin, 1)%matrix, matrix_ks(ispin, 1)%matrix, & 1578 0.0_dp, matrix_tmp, filter_eps=eps_filter) 1579 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_p(ispin, 1)%matrix, & 1580 0.0_dp, matrix_w(ispin, 1)%matrix, filter_eps=eps_filter) 1581 ENDDO 1582 CALL dbcsr_release(matrix_tmp) 1583 1584 IF (nspins == 1) THEN 1585 CALL dbcsr_scale(matrix_p(1, 1)%matrix, 2.0_dp) 1586 CALL dbcsr_scale(matrix_w(1, 1)%matrix, 2.0_dp) 1587 END IF 1588 1589 CALL timestop(handle) 1590 1591 END SUBROUTINE ec_ls_solver 1592 1593! ************************************************************************************************** 1594!> \brief Calculate the energy correction 1595!> \param ec_env ... 1596!> \param unit_nr ... 1597!> \author Creation (03.2014,JGH) 1598! ************************************************************************************************** 1599 SUBROUTINE ec_energy(ec_env, unit_nr) 1600 1601 TYPE(energy_correction_type) :: ec_env 1602 INTEGER, INTENT(IN) :: unit_nr 1603 1604 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_energy', routineP = moduleN//':'//routineN 1605 1606 INTEGER :: handle, ispin, nspins 1607 REAL(KIND=dp) :: eband, trace 1608 1609 CALL timeset(routineN, handle) 1610 1611 SELECT CASE (ec_env%energy_functional) 1612 CASE (ec_functional_harris) 1613 nspins = SIZE(ec_env%matrix_ks, 1) 1614 eband = 0.0_dp 1615 DO ispin = 1, nspins 1616 CALL dbcsr_dot(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_p(ispin, 1)%matrix, trace) 1617 eband = eband + trace 1618 END DO 1619 ec_env%eband = eband + ec_env%efield_nuclear 1620 ec_env%etotal = ec_env%eband + ec_env%ehartree + ec_env%exc - ec_env%vhxc + ec_env%edispersion 1621 IF (unit_nr > 0) THEN 1622 WRITE (unit_nr, '(T3,A,T65,F16.10)') "Eband ", ec_env%eband 1623 WRITE (unit_nr, '(T3,A,T65,F16.10)') "Ehartree ", ec_env%ehartree 1624 WRITE (unit_nr, '(T3,A,T65,F16.10)') "Exc ", ec_env%exc 1625 WRITE (unit_nr, '(T3,A,T65,F16.10)') "Evhxc ", ec_env%vhxc 1626 WRITE (unit_nr, '(T3,A,T65,F16.10)') "Edisp ", ec_env%edispersion 1627 WRITE (unit_nr, '(T3,A,T65,F16.10)') "Etotal Harris Functional ", ec_env%etotal 1628 END IF 1629 1630 CASE DEFAULT 1631 1632 CPASSERT(.FALSE.) 1633 1634 END SELECT 1635 1636 CALL timestop(handle) 1637 1638 END SUBROUTINE ec_energy 1639 1640! ************************************************************************************************** 1641!> \brief builds either the full neighborlist or neighborlists of molecular 1642!> \brief subsets, depending on parameter values 1643!> \param qs_env ... 1644!> \param ec_env ... 1645!> \par History 1646!> 2012.07 created [Martin Haeufel] 1647!> 2016.07 Adapted for Harris functional {JGH] 1648!> \author Martin Haeufel 1649! ************************************************************************************************** 1650 SUBROUTINE ec_build_neighborlist(qs_env, ec_env) 1651 TYPE(qs_environment_type), POINTER :: qs_env 1652 TYPE(energy_correction_type), POINTER :: ec_env 1653 1654 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_neighborlist', & 1655 routineP = moduleN//':'//routineN 1656 1657 INTEGER :: handle, ikind, nkind 1658 LOGICAL :: gth_potential_present, & 1659 sgp_potential_present, & 1660 skip_load_balance_distributed 1661 LOGICAL, ALLOCATABLE, DIMENSION(:) :: orb_present, ppl_present, ppnl_present 1662 REAL(dp) :: subcells 1663 REAL(dp), ALLOCATABLE, DIMENSION(:) :: orb_radius, ppl_radius, ppnl_radius 1664 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius 1665 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1666 TYPE(cell_type), POINTER :: cell 1667 TYPE(dft_control_type), POINTER :: dft_control 1668 TYPE(distribution_1d_type), POINTER :: distribution_1d 1669 TYPE(distribution_2d_type), POINTER :: distribution_2d 1670 TYPE(gth_potential_type), POINTER :: gth_potential 1671 TYPE(gto_basis_set_type), POINTER :: basis_set 1672 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d 1673 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1674 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1675 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1676 TYPE(qs_kind_type), POINTER :: qs_kind 1677 TYPE(qs_ks_env_type), POINTER :: ks_env 1678 TYPE(sgp_potential_type), POINTER :: sgp_potential 1679 1680 CALL timeset(routineN, handle) 1681 1682 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set) 1683 CALL get_qs_kind_set(qs_kind_set, gth_potential_present=gth_potential_present, & 1684 sgp_potential_present=sgp_potential_present) 1685 nkind = SIZE(qs_kind_set) 1686 ALLOCATE (orb_radius(nkind), ppl_radius(nkind), ppnl_radius(nkind)) 1687 ALLOCATE (orb_present(nkind), ppl_present(nkind), ppnl_present(nkind)) 1688 ALLOCATE (pair_radius(nkind, nkind)) 1689 ALLOCATE (atom2d(nkind)) 1690 1691 CALL get_qs_env(qs_env, & 1692 atomic_kind_set=atomic_kind_set, & 1693 cell=cell, & 1694 distribution_2d=distribution_2d, & 1695 local_particles=distribution_1d, & 1696 particle_set=particle_set, & 1697 molecule_set=molecule_set) 1698 1699 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, & 1700 molecule_set, .FALSE., particle_set) 1701 1702 DO ikind = 1, nkind 1703 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list) 1704 qs_kind => qs_kind_set(ikind) 1705 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="HARRIS") 1706 IF (ASSOCIATED(basis_set)) THEN 1707 orb_present(ikind) = .TRUE. 1708 CALL get_gto_basis_set(gto_basis_set=basis_set, kind_radius=orb_radius(ikind)) 1709 ELSE 1710 orb_present(ikind) = .FALSE. 1711 orb_radius(ikind) = 0.0_dp 1712 END IF 1713 CALL get_qs_kind(qs_kind, gth_potential=gth_potential, sgp_potential=sgp_potential) 1714 IF (gth_potential_present .OR. sgp_potential_present) THEN 1715 IF (ASSOCIATED(gth_potential)) THEN 1716 CALL get_potential(potential=gth_potential, & 1717 ppl_present=ppl_present(ikind), & 1718 ppl_radius=ppl_radius(ikind), & 1719 ppnl_present=ppnl_present(ikind), & 1720 ppnl_radius=ppnl_radius(ikind)) 1721 ELSE IF (ASSOCIATED(sgp_potential)) THEN 1722 CALL get_potential(potential=sgp_potential, & 1723 ppl_present=ppl_present(ikind), & 1724 ppl_radius=ppl_radius(ikind), & 1725 ppnl_present=ppnl_present(ikind), & 1726 ppnl_radius=ppnl_radius(ikind)) 1727 ELSE 1728 ppl_present(ikind) = .FALSE. 1729 ppl_radius(ikind) = 0.0_dp 1730 ppnl_present(ikind) = .FALSE. 1731 ppnl_radius(ikind) = 0.0_dp 1732 END IF 1733 END IF 1734 END DO 1735 1736 CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells) 1737 1738 ! overlap 1739 CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius) 1740 CALL build_neighbor_lists(ec_env%sab_orb, particle_set, atom2d, cell, pair_radius, & 1741 subcells=subcells, nlname="sab_orb") 1742 ! pseudopotential 1743 IF (gth_potential_present .OR. sgp_potential_present) THEN 1744 IF (ANY(ppl_present)) THEN 1745 CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius) 1746 CALL build_neighbor_lists(ec_env%sac_ppl, particle_set, atom2d, cell, pair_radius, & 1747 subcells=subcells, operator_type="ABC", nlname="sac_ppl") 1748 END IF 1749 1750 IF (ANY(ppnl_present)) THEN 1751 CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius) 1752 CALL build_neighbor_lists(ec_env%sap_ppnl, particle_set, atom2d, cell, pair_radius, & 1753 subcells=subcells, operator_type="ABBA", nlname="sap_ppnl") 1754 END IF 1755 END IF 1756 1757 ! Release work storage 1758 CALL atom2d_cleanup(atom2d) 1759 DEALLOCATE (atom2d) 1760 DEALLOCATE (orb_present, ppl_present, ppnl_present) 1761 DEALLOCATE (orb_radius, ppl_radius, ppnl_radius) 1762 DEALLOCATE (pair_radius) 1763 1764 ! Task list 1765 CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control) 1766 skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed 1767 IF (ASSOCIATED(ec_env%task_list)) CALL deallocate_task_list(ec_env%task_list) 1768 CALL allocate_task_list(ec_env%task_list) 1769 CALL generate_qs_task_list(ks_env, ec_env%task_list, & 1770 reorder_rs_grid_ranks=.FALSE., soft_valid=.FALSE., & 1771 skip_load_balance_distributed=skip_load_balance_distributed, & 1772 basis_type="HARRIS", sab_orb_external=ec_env%sab_orb) 1773 1774 CALL timestop(handle) 1775 1776 END SUBROUTINE ec_build_neighborlist 1777 1778! ************************************************************************************************** 1779!> \brief ... 1780!> \param qs_env ... 1781!> \param ec_env ... 1782! ************************************************************************************************** 1783 SUBROUTINE ec_properties(qs_env, ec_env) 1784 TYPE(qs_environment_type), POINTER :: qs_env 1785 TYPE(energy_correction_type), POINTER :: ec_env 1786 1787 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_properties', routineP = moduleN//':'//routineN 1788 1789 CHARACTER(LEN=8), DIMENSION(3) :: rlab 1790 CHARACTER(LEN=default_path_length) :: filename 1791 INTEGER :: akind, handle, i, ia, iatom, idir, & 1792 ikind, iounit, ispin, maxmom, nspins, & 1793 reference, unit_nr 1794 LOGICAL :: magnetic, periodic 1795 REAL(KIND=dp) :: charge, dd, focc, tmp 1796 REAL(KIND=dp), DIMENSION(3) :: cdip, pdip, rcc, rdip, ria, tdip 1797 REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point 1798 TYPE(atomic_kind_type), POINTER :: atomic_kind 1799 TYPE(cell_type), POINTER :: cell 1800 TYPE(cp_logger_type), POINTER :: logger 1801 TYPE(cp_para_env_type), POINTER :: para_env 1802 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments 1803 TYPE(distribution_1d_type), POINTER :: local_particles 1804 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1805 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1806 TYPE(section_vals_type), POINTER :: print_key 1807 1808 CALL timeset(routineN, handle) 1809 1810 rlab(1) = "X" 1811 rlab(2) = "Y" 1812 rlab(3) = "Z" 1813 1814 logger => cp_get_default_logger() 1815 iounit = cp_logger_get_default_unit_nr(logger) 1816 1817 print_key => section_vals_get_subs_vals(section_vals=qs_env%input, & 1818 subsection_name="DFT%PRINT%MOMENTS") 1819 1820 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 1821 1822 maxmom = section_get_ival(section_vals=qs_env%input, & 1823 keyword_name="DFT%PRINT%MOMENTS%MAX_MOMENT") 1824 periodic = section_get_lval(section_vals=qs_env%input, & 1825 keyword_name="DFT%PRINT%MOMENTS%PERIODIC") 1826 reference = section_get_ival(section_vals=qs_env%input, & 1827 keyword_name="DFT%PRINT%MOMENTS%REFERENCE") 1828 magnetic = section_get_lval(section_vals=qs_env%input, & 1829 keyword_name="DFT%PRINT%MOMENTS%MAGNETIC") 1830 NULLIFY (ref_point) 1831 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point) 1832 unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=qs_env%input, & 1833 print_key_path="DFT%PRINT%MOMENTS", extension=".dat", & 1834 middle_name="moments", log_filename=.FALSE.) 1835 1836 IF (iounit > 0) THEN 1837 IF (unit_nr /= iounit .AND. unit_nr > 0) THEN 1838 INQUIRE (UNIT=unit_nr, NAME=filename) 1839 WRITE (UNIT=iounit, FMT="(/,T2,A,2(/,T3,A),/)") & 1840 "MOMENTS", "The electric/magnetic moments are written to file:", & 1841 TRIM(filename) 1842 ELSE 1843 WRITE (UNIT=iounit, FMT="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS" 1844 END IF 1845 END IF 1846 1847 IF (periodic) THEN 1848 CPABORT("Periodic moments not implemented with EC") 1849 ELSE 1850 CPASSERT(maxmom < 2) 1851 CPASSERT(.NOT. magnetic) 1852 IF (maxmom == 1) THEN 1853 CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env) 1854 ! reference point 1855 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point) 1856 ! nuclear contribution 1857 cdip = 0.0_dp 1858 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, & 1859 qs_kind_set=qs_kind_set, local_particles=local_particles) 1860 DO ikind = 1, SIZE(local_particles%n_el) 1861 DO ia = 1, local_particles%n_el(ikind) 1862 iatom = local_particles%list(ikind)%array(ia) 1863 ! fold atomic positions back into unit cell 1864 ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc 1865 ria = ria - rcc 1866 atomic_kind => particle_set(iatom)%atomic_kind 1867 CALL get_atomic_kind(atomic_kind, kind_number=akind) 1868 CALL get_qs_kind(qs_kind_set(akind), core_charge=charge) 1869 cdip(1:3) = cdip(1:3) - charge*ria(1:3) 1870 END DO 1871 END DO 1872 CALL mp_sum(cdip, para_env%group) 1873 ! 1874 ! direct density contribution 1875 CALL ec_efield_integrals(qs_env, ec_env, rcc) 1876 ! 1877 nspins = SIZE(ec_env%matrix_p, 1) 1878 pdip = 0.0_dp 1879 DO ispin = 1, nspins 1880 DO idir = 1, 3 1881 CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, & 1882 ec_env%efield%dipmat(idir)%matrix, tmp) 1883 pdip(idir) = pdip(idir) + tmp 1884 END DO 1885 END DO 1886 ! 1887 ! response contribution 1888 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s) 1889 NULLIFY (moments) 1890 CALL dbcsr_allocate_matrix_set(moments, 4) 1891 DO i = 1, 4 1892 ALLOCATE (moments(i)%matrix) 1893 CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments") 1894 CALL dbcsr_set(moments(i)%matrix, 0.0_dp) 1895 END DO 1896 CALL build_local_moment_matrix(qs_env, moments, 1, ref_point=rcc) 1897 ! 1898 focc = 2.0_dp 1899 IF (nspins == 2) focc = 1.0_dp 1900 rdip = 0.0_dp 1901 DO ispin = 1, nspins 1902 DO idir = 1, 3 1903 CALL dbcsr_dot(ec_env%p_env%p1(ispin)%matrix, moments(idir)%matrix, tmp) 1904 rdip(idir) = rdip(idir) + tmp 1905 END DO 1906 END DO 1907 CALL dbcsr_deallocate_matrix_set(moments) 1908 ! 1909 IF (unit_nr > 0) THEN 1910 tdip = -(rdip + pdip + cdip) 1911 WRITE (unit_nr, "(T3,A)") "Dipoles are based on the traditional operator." 1912 dd = SQRT(SUM(tdip(1:3)**2))*debye 1913 WRITE (unit_nr, "(T3,A)") "Dipole moment [Debye]" 1914 WRITE (unit_nr, "(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") & 1915 (TRIM(rlab(i)), "=", tdip(i)*debye, i=1, 3), "Total=", dd 1916 END IF 1917 ENDIF 1918 END IF 1919 1920 CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, & 1921 basis_section=qs_env%input, print_key_path="DFT%PRINT%MOMENTS") 1922 END IF 1923 1924 CALL timestop(handle) 1925 1926 END SUBROUTINE ec_properties 1927 1928! ************************************************************************************************** 1929 1930END MODULE energy_corrections 1931