1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Does all kind of post scf calculations for GPW/GAPW 8!> \par History 9!> Started as a copy from the relevant part of qs_scf 10!> Start to adapt for k-points [07.2015, JGH] 11!> \author Joost VandeVondele (10.2003) 12! ************************************************************************************************** 13MODULE qs_scf_post_gpw 14 USE admm_types, ONLY: admm_type 15 USE admm_utils, ONLY: admm_correct_for_eigenvalues,& 16 admm_uncorrect_for_eigenvalues 17 USE ai_onecenter, ONLY: sg_overlap 18 USE atom_kind_orbitals, ONLY: calculate_atomic_density 19 USE atomic_kind_types, ONLY: atomic_kind_type,& 20 get_atomic_kind 21 USE basis_set_types, ONLY: gto_basis_set_p_type,& 22 gto_basis_set_type 23 USE cell_types, ONLY: cell_type 24 USE cp_array_utils, ONLY: cp_1d_r_p_type 25 USE cp_blacs_env, ONLY: cp_blacs_env_type 26 USE cp_control_types, ONLY: dft_control_type 27 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 28 dbcsr_deallocate_matrix_set 29 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 30 USE cp_ddapc_util, ONLY: get_ddapc 31 USE cp_fm_diag, ONLY: choose_eigv_solver 32 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 33 cp_fm_struct_release,& 34 cp_fm_struct_type 35 USE cp_fm_types, ONLY: cp_fm_create,& 36 cp_fm_get_info,& 37 cp_fm_init_random,& 38 cp_fm_p_type,& 39 cp_fm_release,& 40 cp_fm_to_fm,& 41 cp_fm_type 42 USE cp_fm_vect, ONLY: cp_fm_vect_dealloc 43 USE cp_log_handling, ONLY: cp_get_default_logger,& 44 cp_logger_get_default_io_unit,& 45 cp_logger_type,& 46 cp_to_string 47 USE cp_output_handling, ONLY: cp_p_file,& 48 cp_print_key_finished_output,& 49 cp_print_key_should_output,& 50 cp_print_key_unit_nr 51 USE cp_para_types, ONLY: cp_para_env_type 52 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube 53 USE dbcsr_api, ONLY: & 54 dbcsr_convert_dbcsr_to_csr, dbcsr_copy, dbcsr_csr_create_from_dbcsr, & 55 dbcsr_csr_dbcsr_blkrow_dist, dbcsr_csr_destroy, dbcsr_csr_type, dbcsr_csr_write, & 56 dbcsr_desymmetrize, dbcsr_has_symmetry, dbcsr_p_type, dbcsr_release, dbcsr_type 57 USE dct, ONLY: pw_shrink 58 USE et_coupling_types, ONLY: set_et_coupling_type 59 USE hirshfeld_methods, ONLY: comp_hirshfeld_charges,& 60 comp_hirshfeld_i_charges,& 61 create_shape_function,& 62 write_hirshfeld_charges 63 USE hirshfeld_types, ONLY: create_hirshfeld_type,& 64 hirshfeld_type,& 65 release_hirshfeld_type,& 66 set_hirshfeld_info 67 USE input_constants, ONLY: do_loc_both,& 68 do_loc_homo,& 69 do_loc_lumo,& 70 ot_precond_full_all,& 71 radius_covalent,& 72 radius_user,& 73 ref_charge_atomic,& 74 ref_charge_mulliken 75 USE input_section_types, ONLY: section_get_ival,& 76 section_get_ivals,& 77 section_get_lval,& 78 section_get_rval,& 79 section_vals_get,& 80 section_vals_get_subs_vals,& 81 section_vals_type,& 82 section_vals_val_get 83 USE kinds, ONLY: default_path_length,& 84 default_string_length,& 85 dp 86 USE kpoint_types, ONLY: kpoint_type 87 USE lapack, ONLY: lapack_sgesv 88 USE mao_wfn_analysis, ONLY: mao_analysis 89 USE mathconstants, ONLY: pi 90 USE memory_utilities, ONLY: reallocate 91 USE message_passing, ONLY: mp_sum 92 USE minbas_wfn_analysis, ONLY: minbas_analysis 93 USE molden_utils, ONLY: write_mos_molden 94 USE molecule_types, ONLY: molecule_type 95 USE mulliken, ONLY: mulliken_charges 96 USE orbital_pointers, ONLY: indso 97 USE particle_list_types, ONLY: particle_list_type 98 USE particle_types, ONLY: particle_type 99 USE physcon, ONLY: angstrom,& 100 evolt 101 USE population_analyses, ONLY: lowdin_population_analysis,& 102 mulliken_population_analysis 103 USE preconditioner_types, ONLY: preconditioner_type 104 USE ps_implicit_types, ONLY: MIXED_BC,& 105 MIXED_PERIODIC_BC,& 106 NEUMANN_BC,& 107 PERIODIC_BC 108 USE pw_env_types, ONLY: pw_env_get,& 109 pw_env_type 110 USE pw_grids, ONLY: get_pw_grid_info 111 USE pw_methods, ONLY: pw_axpy,& 112 pw_copy,& 113 pw_derive,& 114 pw_integrate_function,& 115 pw_scale,& 116 pw_transfer,& 117 pw_zero 118 USE pw_poisson_methods, ONLY: pw_poisson_solve 119 USE pw_poisson_types, ONLY: pw_poisson_implicit,& 120 pw_poisson_type 121 USE pw_pool_types, ONLY: pw_pool_create_pw,& 122 pw_pool_give_back_pw,& 123 pw_pool_p_type,& 124 pw_pool_type 125 USE pw_types, ONLY: COMPLEXDATA1D,& 126 REALDATA3D,& 127 REALSPACE,& 128 RECIPROCALSPACE,& 129 pw_p_type,& 130 pw_type 131 USE qs_active_space_methods, ONLY: active_space_main 132 USE qs_collocate_density, ONLY: calculate_wavefunction 133 USE qs_commutators, ONLY: build_com_hr_matrix 134 USE qs_core_energies, ONLY: calculate_ptrace 135 USE qs_electric_field_gradient, ONLY: qs_efg_calc 136 USE qs_elf_methods, ONLY: qs_elf_calc 137 USE qs_energy_types, ONLY: qs_energy_type 138 USE qs_energy_window, ONLY: energy_windows 139 USE qs_environment_types, ONLY: get_qs_env,& 140 qs_environment_type,& 141 set_qs_env 142 USE qs_epr_hyp, ONLY: qs_epr_hyp_calc 143 USE qs_grid_atom, ONLY: grid_atom_type 144 USE qs_integral_utils, ONLY: basis_set_list_setup 145 USE qs_kind_types, ONLY: get_qs_kind,& 146 qs_kind_type 147 USE qs_ks_methods, ONLY: calc_rho_tot_gspace,& 148 qs_ks_update_qs_env 149 USE qs_ks_types, ONLY: qs_ks_did_change 150 USE qs_loc_dipole, ONLY: loc_dipole 151 USE qs_loc_states, ONLY: get_localization_info 152 USE qs_loc_types, ONLY: qs_loc_env_create,& 153 qs_loc_env_destroy,& 154 qs_loc_env_new_type 155 USE qs_loc_utils, ONLY: loc_write_restart,& 156 qs_loc_control_init,& 157 qs_loc_env_init,& 158 qs_loc_init,& 159 retain_history 160 USE qs_local_properties, ONLY: qs_local_energy,& 161 qs_local_stress 162 USE qs_mo_io, ONLY: write_dm_binary_restart,& 163 write_mo_set 164 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues,& 165 make_mo_eig 166 USE qs_mo_occupation, ONLY: set_mo_occupation 167 USE qs_mo_types, ONLY: get_mo_set,& 168 mo_set_p_type 169 USE qs_moments, ONLY: qs_moment_berry_phase,& 170 qs_moment_locop 171 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 172 get_neighbor_list_set_p,& 173 neighbor_list_iterate,& 174 neighbor_list_iterator_create,& 175 neighbor_list_iterator_p_type,& 176 neighbor_list_iterator_release,& 177 neighbor_list_set_p_type 178 USE qs_ot_eigensolver, ONLY: ot_eigensolver 179 USE qs_pdos, ONLY: calculate_projected_dos 180 USE qs_resp, ONLY: resp_fit 181 USE qs_rho0_types, ONLY: get_rho0_mpole,& 182 mpole_rho_atom,& 183 rho0_mpole_type 184 USE qs_rho_atom_types, ONLY: rho_atom_type 185 USE qs_rho_types, ONLY: qs_rho_get,& 186 qs_rho_type 187 USE qs_scf_types, ONLY: ot_method_nr,& 188 qs_scf_env_type 189 USE qs_scf_wfn_mix, ONLY: wfn_mix 190 USE qs_subsys_types, ONLY: qs_subsys_get,& 191 qs_subsys_type 192 USE qs_wannier90, ONLY: wannier90_interface 193 USE s_square_methods, ONLY: compute_s_square 194 USE scf_control_types, ONLY: scf_control_type 195 USE stm_images, ONLY: th_stm_image 196 USE transport, ONLY: qs_scf_post_transport 197 USE virial_types, ONLY: virial_type 198 USE xray_diffraction, ONLY: calculate_rhotot_elec_gspace,& 199 xray_diffraction_spectrum 200#include "./base/base_uses.f90" 201 202 IMPLICIT NONE 203 PRIVATE 204 205 ! Global parameters 206 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_post_gpw' 207 PUBLIC :: scf_post_calculation_gpw, & 208 qs_scf_post_moments, & 209 write_available_results, & 210 write_mo_free_results 211 PUBLIC :: make_lumo 212 213! ************************************************************************************************** 214 215CONTAINS 216 217! ************************************************************************************************** 218!> \brief collects possible post - scf calculations and prints info / computes properties. 219!> \param qs_env the qs_env in which the qs_env lives 220!> \param wf_type ... 221!> \par History 222!> 02.2003 created [fawzi] 223!> 10.2004 moved here from qs_scf [Joost VandeVondele] 224!> started splitting out different subroutines 225!> 10.2015 added header for wave-function correlated methods [Vladimir Rybkin] 226!> \author fawzi 227!> \note 228!> this function changes mo_eigenvectors and mo_eigenvalues, depending on the print keys. 229!> In particular, MO_CUBES causes the MOs to be rotated to make them eigenstates of the KS 230!> matrix, and mo_eigenvalues is updated accordingly. This can, for unconverged wavefunctions, 231!> change afterwards slightly the forces (hence small numerical differences between MD 232!> with and without the debug print level). Ideally this should not happen... 233! ************************************************************************************************** 234 SUBROUTINE scf_post_calculation_gpw(qs_env, wf_type) 235 236 TYPE(qs_environment_type), POINTER :: qs_env 237 CHARACTER(6), OPTIONAL :: wf_type 238 239 CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_gpw', & 240 routineP = moduleN//':'//routineN 241 242 INTEGER :: handle, homo, ispin, min_lumos, n_rep, & 243 nhomo, nlumo, nlumo_stm, nlumo_tddft, & 244 nlumos, nmo, output_unit, unit_nr 245 INTEGER, DIMENSION(:, :, :), POINTER :: marked_states 246 LOGICAL :: check_write, compute_lumos, do_homo, do_kpoints, do_mo_cubes, do_stm, & 247 do_wannier_cubes, has_homo, has_lumo, loc_explicit, loc_print_explicit, my_localized_wfn, & 248 p_loc, p_loc_homo, p_loc_lumo 249 REAL(dp) :: e_kin 250 REAL(KIND=dp) :: gap, homo_lumo(2, 2) 251 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues 252 TYPE(admm_type), POINTER :: admm_env 253 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 254 TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: occupied_evals, unoccupied_evals, & 255 unoccupied_evals_stm 256 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: homo_localized, lumo_localized, lumo_ptr, & 257 mo_loc_history, occupied_orbs, unoccupied_orbs, unoccupied_orbs_stm 258 TYPE(cp_fm_type), POINTER :: mo_coeff 259 TYPE(cp_logger_type), POINTER :: logger 260 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, matrix_s, mo_derivs 261 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: kinetic_m, rho_ao 262 TYPE(dft_control_type), POINTER :: dft_control 263 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 264 TYPE(molecule_type), POINTER :: molecule_set(:) 265 TYPE(particle_list_type), POINTER :: particles 266 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 267 TYPE(pw_env_type), POINTER :: pw_env 268 TYPE(pw_p_type) :: wf_g, wf_r 269 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 270 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 271 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 272 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env_homo, qs_loc_env_lumo 273 TYPE(qs_rho_type), POINTER :: rho 274 TYPE(qs_scf_env_type), POINTER :: scf_env 275 TYPE(qs_subsys_type), POINTER :: subsys 276 TYPE(scf_control_type), POINTER :: scf_control 277 TYPE(section_vals_type), POINTER :: dft_section, input, loc_print_section, & 278 localize_section, print_key, & 279 stm_section 280 281 CALL timeset(routineN, handle) 282 283 logger => cp_get_default_logger() 284 output_unit = cp_logger_get_default_io_unit(logger) 285 286 ! Print out the type of wavefunction to distinguish between SCF and post-SCF 287 IF (PRESENT(wf_type)) THEN 288 IF (output_unit > 0) THEN 289 WRITE (UNIT=output_unit, FMT='(/,(T1,A))') REPEAT("-", 40) 290 WRITE (UNIT=output_unit, FMT='(/,(T3,A,T19,A,T25,A))') "Properties from ", wf_type, " density" 291 WRITE (UNIT=output_unit, FMT='(/,(T1,A))') REPEAT("-", 40) 292 ENDIF 293 ENDIF 294 295 ! Writes the data that is already available in qs_env 296 CALL get_qs_env(qs_env, scf_env=scf_env) 297 CALL write_available_results(qs_env, scf_env) 298 299 my_localized_wfn = .FALSE. 300 NULLIFY (admm_env, dft_control, pw_env, auxbas_pw_pool, pw_pools, mos, rho, & 301 mo_coeff, ks_rmpv, matrix_s, qs_loc_env_homo, qs_loc_env_lumo, scf_control, & 302 unoccupied_orbs, unoccupied_orbs_stm, mo_eigenvalues, unoccupied_evals, & 303 unoccupied_evals_stm, molecule_set, mo_derivs, & 304 subsys, particles, input, print_key, kinetic_m, marked_states) 305 NULLIFY (homo_localized, lumo_localized, lumo_ptr, rho_ao) 306 307 has_homo = .FALSE. 308 has_lumo = .FALSE. 309 p_loc = .FALSE. 310 p_loc_homo = .FALSE. 311 p_loc_lumo = .FALSE. 312 313 CPASSERT(ASSOCIATED(scf_env)) 314 CPASSERT(scf_env%ref_count > 0) 315 CPASSERT(ASSOCIATED(qs_env)) 316 ! Here we start with data that needs a postprocessing... 317 CALL get_qs_env(qs_env, & 318 dft_control=dft_control, & 319 molecule_set=molecule_set, & 320 scf_control=scf_control, & 321 do_kpoints=do_kpoints, & 322 input=input, & 323 subsys=subsys, & 324 rho=rho, & 325 pw_env=pw_env, & 326 particle_set=particle_set, & 327 atomic_kind_set=atomic_kind_set, & 328 qs_kind_set=qs_kind_set) 329 CALL qs_subsys_get(subsys, particles=particles) 330 331 CALL qs_rho_get(rho, rho_ao_kp=rho_ao) 332 333 ! In MP2 case update the Hartree potential 334 IF (ASSOCIATED(qs_env%mp2_env)) CALL update_hartree_with_mp2(rho, qs_env) 335 336 ! **** the kinetic energy 337 IF (cp_print_key_should_output(logger%iter_info, input, & 338 "DFT%PRINT%KINETIC_ENERGY") /= 0) THEN 339 CALL get_qs_env(qs_env, kinetic_kp=kinetic_m) 340 CPASSERT(ASSOCIATED(kinetic_m)) 341 CPASSERT(ASSOCIATED(kinetic_m(1, 1)%matrix)) 342 CALL calculate_ptrace(kinetic_m, rho_ao, e_kin, dft_control%nspins) 343 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%KINETIC_ENERGY", & 344 extension=".Log") 345 IF (unit_nr > 0) THEN 346 WRITE (unit_nr, '(T3,A,T55,F25.14)') "Electronic kinetic energy:", e_kin 347 ENDIF 348 CALL cp_print_key_finished_output(unit_nr, logger, input, & 349 "DFT%PRINT%KINETIC_ENERGY") 350 END IF 351 352 ! Atomic Charges that require further computation 353 CALL qs_scf_post_charges(input, logger, qs_env) 354 355 ! Moments of charge distribution 356 CALL qs_scf_post_moments(input, logger, qs_env, output_unit) 357 358 ! Determine if we need to computer properties using the localized centers 359 dft_section => section_vals_get_subs_vals(input, "DFT") 360 localize_section => section_vals_get_subs_vals(dft_section, "LOCALIZE") 361 loc_print_section => section_vals_get_subs_vals(localize_section, "PRINT") 362 CALL section_vals_get(localize_section, explicit=loc_explicit) 363 CALL section_vals_get(loc_print_section, explicit=loc_print_explicit) 364 365 ! Print_keys controlled by localization 366 IF (loc_print_explicit) THEN 367 print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_DIPOLES") 368 p_loc = BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) 369 print_key => section_vals_get_subs_vals(loc_print_section, "TOTAL_DIPOLE") 370 p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) 371 print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_CENTERS") 372 p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) 373 print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_SPREADS") 374 p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) 375 print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_CUBES") 376 p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) 377 print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_STATES") 378 p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) 379 print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_MOMENTS") 380 p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) 381 ELSE 382 p_loc = .FALSE. 383 END IF 384 IF (loc_explicit) THEN 385 p_loc_homo = (section_get_ival(localize_section, "STATES") == do_loc_homo .OR. & 386 section_get_ival(localize_section, "STATES") == do_loc_both) .AND. p_loc 387 p_loc_lumo = (section_get_ival(localize_section, "STATES") == do_loc_lumo .OR. & 388 section_get_ival(localize_section, "STATES") == do_loc_both) .AND. p_loc 389 CALL section_vals_val_get(localize_section, "LIST_UNOCCUPIED", n_rep_val=n_rep) 390 ELSE 391 p_loc_homo = .FALSE. 392 p_loc_lumo = .FALSE. 393 n_rep = 0 394 END IF 395 396 IF (n_rep == 0 .AND. p_loc_lumo) THEN 397 CALL cp_abort(__LOCATION__, "No LIST_UNOCCUPIED was specified, "// & 398 "therefore localization of unoccupied states will be skipped!") 399 p_loc_lumo = .FALSE. 400 END IF 401 print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_STATES") 402 p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) 403 404 ! Control for STM 405 stm_section => section_vals_get_subs_vals(input, "DFT%PRINT%STM") 406 CALL section_vals_get(stm_section, explicit=do_stm) 407 nlumo_stm = 0 408 IF (do_stm) nlumo_stm = section_get_ival(stm_section, "NLUMO") 409 410 ! check for CUBES (MOs and WANNIERS) 411 do_mo_cubes = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_CUBES") & 412 , cp_p_file) 413 IF (loc_print_explicit) THEN 414 do_wannier_cubes = BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, & 415 "WANNIER_CUBES"), cp_p_file) 416 ELSE 417 do_wannier_cubes = .FALSE. 418 END IF 419 nlumo = section_get_ival(dft_section, "PRINT%MO_CUBES%NLUMO") 420 nhomo = section_get_ival(dft_section, "PRINT%MO_CUBES%NHOMO") 421 nlumo_tddft = 0 422 IF (dft_control%do_tddfpt_calculation) THEN 423 nlumo_tddft = section_get_ival(dft_section, "TDDFPT%NLUMO") 424 END IF 425 426 ! Setup the grids needed to compute a wavefunction given a vector.. 427 IF (((do_mo_cubes .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc) THEN 428 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 429 pw_pools=pw_pools) 430 CALL pw_pool_create_pw(auxbas_pw_pool, wf_r%pw, & 431 use_data=REALDATA3D, & 432 in_space=REALSPACE) 433 CALL pw_pool_create_pw(auxbas_pw_pool, wf_g%pw, & 434 use_data=COMPLEXDATA1D, & 435 in_space=RECIPROCALSPACE) 436 END IF 437 438 ! Makes the MOs eigenstates, computes eigenvalues, write cubes 439 IF (do_kpoints) THEN 440 IF (do_mo_cubes) THEN 441 CPWARN("Print MO Cubes not implemented for k-point calculations!!") 442 END IF 443 ELSE 444 CALL get_qs_env(qs_env, & 445 mos=mos, & 446 matrix_ks=ks_rmpv) 447 IF ((do_mo_cubes .AND. nhomo /= 0) .OR. do_stm .OR. dft_control%do_tddfpt_calculation) THEN 448 IF (dft_control%restricted) THEN 449 IF (output_unit > 0) WRITE (output_unit, *) & 450 " Unclear how we define MOs in the restricted case ... skipping" 451 ELSE 452 CALL get_qs_env(qs_env, mo_derivs=mo_derivs) 453 IF (dft_control%do_admm) THEN 454 CALL get_qs_env(qs_env, admm_env=admm_env) 455 CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs, admm_env=admm_env) 456 ELSE 457 CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs) 458 END IF 459 END IF 460 DO ispin = 1, dft_control%nspins 461 CALL get_mo_set(mo_set=mos(ispin)%mo_set, eigenvalues=mo_eigenvalues, homo=homo) 462 homo_lumo(ispin, 1) = mo_eigenvalues(homo) 463 END DO 464 has_homo = .TRUE. 465 END IF 466 IF (do_mo_cubes .AND. nhomo /= 0) THEN 467 DO ispin = 1, dft_control%nspins 468 ! Prints the cube files of OCCUPIED ORBITALS 469 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, & 470 eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo) 471 CALL qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env, & 472 mo_coeff, wf_g, wf_r, particles, homo, ispin) 473 END DO 474 ENDIF 475 ENDIF 476 477 ! Initialize the localization environment, needed e.g. for wannier functions and molecular states 478 ! Gets localization info for the occupied orbs 479 ! - Possibly gets wannier functions 480 ! - Possibly gets molecular states 481 IF (p_loc_homo) THEN 482 IF (do_kpoints) THEN 483 CPWARN("Localization not implemented for k-point calculations!!") 484 ELSEIF (dft_control%restricted) THEN 485 IF (output_unit > 0) WRITE (output_unit, *) & 486 " Unclear how we define MOs / localization in the restricted case ... skipping" 487 ELSE 488 ALLOCATE (occupied_orbs(dft_control%nspins)) 489 ALLOCATE (occupied_evals(dft_control%nspins)) 490 ALLOCATE (homo_localized(dft_control%nspins)) 491 DO ispin = 1, dft_control%nspins 492 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, & 493 eigenvalues=mo_eigenvalues) 494 occupied_orbs(ispin)%matrix => mo_coeff 495 occupied_evals(ispin)%array => mo_eigenvalues 496 CALL cp_fm_create(homo_localized(ispin)%matrix, occupied_orbs(ispin)%matrix%matrix_struct) 497 CALL cp_fm_to_fm(occupied_orbs(ispin)%matrix, homo_localized(ispin)%matrix) 498 END DO 499 500 CALL get_qs_env(qs_env, mo_loc_history=mo_loc_history) 501 do_homo = .TRUE. 502 503 CALL qs_loc_env_create(qs_loc_env_homo) 504 CALL qs_loc_control_init(qs_loc_env_homo, localize_section, do_homo=do_homo) 505 CALL qs_loc_init(qs_env, qs_loc_env_homo, localize_section, homo_localized, do_homo, & 506 do_mo_cubes, mo_loc_history=mo_loc_history) 507 CALL get_localization_info(qs_env, qs_loc_env_homo, localize_section, homo_localized, & 508 wf_r, wf_g, particles, occupied_orbs, occupied_evals, marked_states) 509 510 !retain the homo_localized for future use 511 IF (qs_loc_env_homo%localized_wfn_control%use_history) THEN 512 CALL retain_history(mo_loc_history, homo_localized) 513 CALL set_qs_env(qs_env, mo_loc_history=mo_loc_history) 514 ENDIF 515 516 !write restart for localization of occupied orbitals 517 CALL loc_write_restart(qs_loc_env_homo, loc_print_section, mos, & 518 homo_localized, do_homo) 519 CALL cp_fm_vect_dealloc(homo_localized) 520 DEALLOCATE (occupied_orbs) 521 DEALLOCATE (occupied_evals) 522 ! Print Total Dipole if the localization has been performed 523 IF (qs_loc_env_homo%do_localize) THEN 524 CALL loc_dipole(input, dft_control, qs_loc_env_homo, logger, qs_env) 525 END IF 526 END IF 527 ENDIF 528 529 ! Gets the lumos, and eigenvalues for the lumos, and localize them if requested 530 IF (do_kpoints) THEN 531 IF (do_mo_cubes .OR. p_loc_lumo) THEN 532 ! nothing at the moment, not implemented 533 CPWARN("Localization and MO related output not implemented for k-point calculations!") 534 END IF 535 ELSE 536 IF (nlumo .GT. -1) THEN 537 nlumo = MAX(nlumo, nlumo_tddft) 538 END IF 539 compute_lumos = (do_mo_cubes .OR. dft_control%do_tddfpt_calculation) .AND. nlumo .NE. 0 540 compute_lumos = compute_lumos .OR. p_loc_lumo 541 542 DO ispin = 1, dft_control%nspins 543 CALL get_mo_set(mo_set=mos(ispin)%mo_set, homo=homo, nmo=nmo) 544 compute_lumos = compute_lumos .AND. homo == nmo 545 END DO 546 547 IF (do_mo_cubes .AND. .NOT. compute_lumos) THEN 548 549 nlumo = section_get_ival(dft_section, "PRINT%MO_CUBES%NLUMO") 550 DO ispin = 1, dft_control%nspins 551 552 CALL get_mo_set(mo_set=mos(ispin)%mo_set, homo=homo, nmo=nmo, eigenvalues=mo_eigenvalues) 553 IF (nlumo > nmo - homo) THEN 554 ! this case not yet implemented 555 ELSE 556 IF (nlumo .EQ. -1) THEN 557 nlumo = nmo - homo 558 END IF 559 IF (output_unit > 0) WRITE (output_unit, *) " " 560 IF (output_unit > 0) WRITE (output_unit, *) " Lowest eigenvalues of the unoccupied subspace spin ", ispin 561 IF (output_unit > 0) WRITE (output_unit, *) "---------------------------------------------" 562 IF (output_unit > 0) WRITE (output_unit, '(4(1X,1F16.8))') mo_eigenvalues(homo + 1:homo + nlumo) 563 564 ! Prints the cube files of UNOCCUPIED ORBITALS 565 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 566 CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, & 567 mo_coeff, wf_g, wf_r, particles, nlumo, homo, ispin, lumo=homo + 1) 568 END IF 569 END DO 570 571 END IF 572 573 IF (compute_lumos) THEN 574 check_write = .TRUE. 575 min_lumos = nlumo 576 IF (nlumo == 0) check_write = .FALSE. 577 IF (p_loc_lumo) THEN 578 do_homo = .FALSE. 579 CALL qs_loc_env_create(qs_loc_env_lumo) 580 CALL qs_loc_control_init(qs_loc_env_lumo, localize_section, do_homo=do_homo) 581 min_lumos = MAX(MAXVAL(qs_loc_env_lumo%localized_wfn_control%loc_states(:, :)), nlumo) 582 END IF 583 584 ALLOCATE (unoccupied_orbs(dft_control%nspins)) 585 ALLOCATE (unoccupied_evals(dft_control%nspins)) 586 CALL make_lumo(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, min_lumos, nlumos) 587 lumo_ptr => unoccupied_orbs 588 DO ispin = 1, dft_control%nspins 589 has_lumo = .TRUE. 590 homo_lumo(ispin, 2) = unoccupied_evals(ispin)%array(1) 591 CALL get_mo_set(mo_set=mos(ispin)%mo_set, homo=homo) 592 IF (check_write) THEN 593 IF (p_loc_lumo .AND. nlumo .NE. -1) nlumos = MIN(nlumo, nlumos) 594 ! Prints the cube files of UNOCCUPIED ORBITALS 595 CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, & 596 unoccupied_orbs(ispin)%matrix, wf_g, wf_r, particles, nlumos, homo, ispin) 597 END IF 598 END DO 599 600 ! Save the info for tddfpt calculation 601 IF (dft_control%do_tddfpt_calculation) THEN 602 ALLOCATE (dft_control%tddfpt_control%lumos_eigenvalues(nlumos, dft_control%nspins)) 603 DO ispin = 1, dft_control%nspins 604 dft_control%tddfpt_control%lumos_eigenvalues(1:nlumos, ispin) = & 605 unoccupied_evals(ispin)%array(1:nlumos) 606 END DO 607 dft_control%tddfpt_control%lumos => unoccupied_orbs 608 END IF 609 610 IF (p_loc_lumo) THEN 611 ALLOCATE (lumo_localized(dft_control%nspins)) 612 DO ispin = 1, dft_control%nspins 613 CALL cp_fm_create(lumo_localized(ispin)%matrix, unoccupied_orbs(ispin)%matrix%matrix_struct) 614 CALL cp_fm_to_fm(unoccupied_orbs(ispin)%matrix, lumo_localized(ispin)%matrix) 615 END DO 616 CALL qs_loc_init(qs_env, qs_loc_env_lumo, localize_section, lumo_localized, do_homo, do_mo_cubes, & 617 evals=unoccupied_evals) 618 CALL qs_loc_env_init(qs_loc_env_lumo, qs_loc_env_lumo%localized_wfn_control, qs_env, & 619 loc_coeff=unoccupied_orbs) 620 CALL get_localization_info(qs_env, qs_loc_env_lumo, localize_section, & 621 lumo_localized, wf_r, wf_g, particles, & 622 unoccupied_orbs, unoccupied_evals, marked_states) 623 CALL loc_write_restart(qs_loc_env_lumo, loc_print_section, mos, homo_localized, do_homo, & 624 evals=unoccupied_evals) 625 lumo_ptr => lumo_localized 626 END IF 627 ENDIF 628 629 IF (has_homo .AND. has_lumo) THEN 630 IF (output_unit > 0) WRITE (output_unit, *) " " 631 DO ispin = 1, dft_control%nspins 632 IF (.NOT. scf_control%smear%do_smear) THEN 633 gap = homo_lumo(ispin, 2) - homo_lumo(ispin, 1) 634 IF (output_unit > 0) WRITE (output_unit, '(T2,A,F12.6)') & 635 "HOMO - LUMO gap [eV] :", gap*evolt 636 END IF 637 ENDDO 638 ENDIF 639 ENDIF 640 641 ! Deallocate grids needed to compute wavefunctions 642 IF (((do_mo_cubes .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc) THEN 643 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_r%pw) 644 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_g%pw) 645 END IF 646 647 ! Destroy the localization environment 648 IF (.NOT. do_kpoints) THEN 649 IF (p_loc_homo) CALL qs_loc_env_destroy(qs_loc_env_homo) 650 IF (p_loc_lumo) CALL qs_loc_env_destroy(qs_loc_env_lumo) 651 END IF 652 653 ! generate a mix of wfns, and write to a restart 654 IF (do_kpoints) THEN 655 ! nothing at the moment, not implemented 656 ELSE 657 CALL get_qs_env(qs_env, matrix_s=matrix_s) 658 CALL wfn_mix(mos, particle_set, dft_section, qs_kind_set, & 659 lumo_ptr, scf_env, matrix_s, output_unit, marked_states) 660 IF (p_loc_lumo) CALL cp_fm_vect_dealloc(lumo_localized) 661 END IF 662 IF (ASSOCIATED(marked_states)) THEN 663 DEALLOCATE (marked_states) 664 END IF 665 666 ! This is just a deallocation for printing MO_CUBES or TDDFPT 667 IF (.NOT. do_kpoints) THEN 668 IF (compute_lumos) THEN 669 DO ispin = 1, dft_control%nspins 670 DEALLOCATE (unoccupied_evals(ispin)%array) 671 IF (.NOT. dft_control%do_tddfpt_calculation) & 672 CALL cp_fm_release(unoccupied_orbs(ispin)%matrix) 673 ENDDO 674 DEALLOCATE (unoccupied_evals) 675 IF (.NOT. dft_control%do_tddfpt_calculation) THEN 676 DEALLOCATE (unoccupied_orbs) 677 END IF 678 ENDIF 679 ENDIF 680 681 !stm images 682 IF (do_stm) THEN 683 IF (do_kpoints) THEN 684 CPWARN("STM not implemented for k-point calculations!") 685 ELSE 686 IF (nlumo_stm > 0) THEN 687 ALLOCATE (unoccupied_orbs_stm(dft_control%nspins)) 688 ALLOCATE (unoccupied_evals_stm(dft_control%nspins)) 689 CALL make_lumo(qs_env, scf_env, unoccupied_orbs_stm, unoccupied_evals_stm, & 690 nlumo_stm, nlumos) 691 END IF 692 693 CALL th_stm_image(qs_env, stm_section, particles, unoccupied_orbs_stm, & 694 unoccupied_evals_stm) 695 696 IF (nlumo_stm > 0) THEN 697 DO ispin = 1, dft_control%nspins 698 DEALLOCATE (unoccupied_evals_stm(ispin)%array) 699 CALL cp_fm_release(unoccupied_orbs_stm(ispin)%matrix) 700 ENDDO 701 DEALLOCATE (unoccupied_evals_stm) 702 DEALLOCATE (unoccupied_orbs_stm) 703 END IF 704 END IF 705 END IF 706 707 ! Print coherent X-ray diffraction spectrum 708 CALL qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit) 709 710 ! Calculation of Electric Field Gradients 711 CALL qs_scf_post_efg(input, logger, qs_env) 712 713 ! Calculation of ET 714 CALL qs_scf_post_et(input, qs_env, dft_control) 715 716 ! Calculation of EPR Hyperfine Coupling Tensors 717 CALL qs_scf_post_epr(input, logger, qs_env) 718 719 ! Calculation of properties needed for BASIS_MOLOPT optimizations 720 CALL qs_scf_post_molopt(input, logger, qs_env) 721 722 ! Calculate ELF 723 CALL qs_scf_post_elf(input, logger, qs_env) 724 725 ! Use Wannier90 interface 726 CALL wannier90_interface(input, logger, qs_env) 727 728 ! Do active space calculation 729 CALL active_space_main(input, logger, qs_env) 730 731 CALL timestop(handle) 732 733 END SUBROUTINE scf_post_calculation_gpw 734 735! ************************************************************************************************** 736!> \brief Gets the lumos, and eigenvalues for the lumos 737!> \param qs_env ... 738!> \param scf_env ... 739!> \param unoccupied_orbs ... 740!> \param unoccupied_evals ... 741!> \param nlumo ... 742!> \param nlumos ... 743! ************************************************************************************************** 744 SUBROUTINE make_lumo(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, nlumo, nlumos) 745 746 TYPE(qs_environment_type), POINTER :: qs_env 747 TYPE(qs_scf_env_type), POINTER :: scf_env 748 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: unoccupied_orbs 749 TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: unoccupied_evals 750 INTEGER, INTENT(IN) :: nlumo 751 INTEGER, INTENT(OUT) :: nlumos 752 753 CHARACTER(len=*), PARAMETER :: routineN = 'make_lumo', routineP = moduleN//':'//routineN 754 755 INTEGER :: handle, homo, ispin, n, nao, nmo, & 756 output_unit 757 TYPE(admm_type), POINTER :: admm_env 758 TYPE(cp_blacs_env_type), POINTER :: blacs_env 759 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 760 TYPE(cp_fm_type), POINTER :: mo_coeff 761 TYPE(cp_logger_type), POINTER :: logger 762 TYPE(cp_para_env_type), POINTER :: para_env 763 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, matrix_s 764 TYPE(dft_control_type), POINTER :: dft_control 765 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 766 TYPE(preconditioner_type), POINTER :: local_preconditioner 767 TYPE(scf_control_type), POINTER :: scf_control 768 769 CALL timeset(routineN, handle) 770 771 NULLIFY (mos, ks_rmpv, scf_control, dft_control, admm_env, para_env, blacs_env) 772 CALL get_qs_env(qs_env, & 773 mos=mos, & 774 matrix_ks=ks_rmpv, & 775 scf_control=scf_control, & 776 dft_control=dft_control, & 777 matrix_s=matrix_s, & 778 admm_env=admm_env, & 779 para_env=para_env, & 780 blacs_env=blacs_env) 781 782 logger => cp_get_default_logger() 783 output_unit = cp_logger_get_default_io_unit(logger) 784 785 DO ispin = 1, dft_control%nspins 786 NULLIFY (unoccupied_orbs(ispin)%matrix) 787 NULLIFY (unoccupied_evals(ispin)%array) 788 ! Always write eigenvalues 789 IF (output_unit > 0) WRITE (output_unit, *) " " 790 IF (output_unit > 0) WRITE (output_unit, *) " Lowest Eigenvalues of the unoccupied subspace spin ", ispin 791 IF (output_unit > 0) WRITE (output_unit, FMT='(1X,A)') "-----------------------------------------------------" 792 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, homo=homo, nao=nao, nmo=nmo) 793 CALL cp_fm_get_info(mo_coeff, nrow_global=n) 794 nlumos = MAX(1, MIN(nlumo, nao - nmo)) 795 IF (nlumo == -1) nlumos = nao - nmo 796 ALLOCATE (unoccupied_evals(ispin)%array(nlumos)) 797 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, & 798 nrow_global=n, ncol_global=nlumos) 799 CALL cp_fm_create(unoccupied_orbs(ispin)%matrix, fm_struct_tmp, name="lumos") 800 CALL cp_fm_struct_release(fm_struct_tmp) 801 CALL cp_fm_init_random(unoccupied_orbs(ispin)%matrix, nlumos) 802 803 ! the full_all preconditioner makes not much sense for lumos search 804 NULLIFY (local_preconditioner) 805 IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN 806 local_preconditioner => scf_env%ot_preconditioner(1)%preconditioner 807 ! this one can for sure not be right (as it has to match a given C0) 808 IF (local_preconditioner%in_use == ot_precond_full_all) THEN 809 NULLIFY (local_preconditioner) 810 ENDIF 811 ENDIF 812 813 ! ** If we do ADMM, we add have to modify the kohn-sham matrix 814 IF (dft_control%do_admm) THEN 815 CALL admm_correct_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix) 816 END IF 817 818 CALL ot_eigensolver(matrix_h=ks_rmpv(ispin)%matrix, matrix_s=matrix_s(1)%matrix, & 819 matrix_c_fm=unoccupied_orbs(ispin)%matrix, & 820 matrix_orthogonal_space_fm=mo_coeff, & 821 eps_gradient=scf_control%eps_lumos, & 822 preconditioner=local_preconditioner, & 823 iter_max=scf_control%max_iter_lumos, & 824 size_ortho_space=nmo) 825 826 CALL calculate_subspace_eigenvalues(unoccupied_orbs(ispin)%matrix, ks_rmpv(ispin)%matrix, & 827 unoccupied_evals(ispin)%array, scr=output_unit, & 828 ionode=output_unit > 0) 829 830 ! ** If we do ADMM, we restore the original kohn-sham matrix 831 IF (dft_control%do_admm) THEN 832 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix) 833 END IF 834 835 END DO 836 837 CALL timestop(handle) 838 839 END SUBROUTINE make_lumo 840! ************************************************************************************************** 841!> \brief Computes and Prints Atomic Charges with several methods 842!> \param input ... 843!> \param logger ... 844!> \param qs_env the qs_env in which the qs_env lives 845! ************************************************************************************************** 846 SUBROUTINE qs_scf_post_charges(input, logger, qs_env) 847 TYPE(section_vals_type), POINTER :: input 848 TYPE(cp_logger_type), POINTER :: logger 849 TYPE(qs_environment_type), POINTER :: qs_env 850 851 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_charges', & 852 routineP = moduleN//':'//routineN 853 854 INTEGER :: handle, print_level, unit_nr 855 LOGICAL :: do_kpoints, print_it 856 TYPE(section_vals_type), POINTER :: density_fit_section, print_key 857 858 CALL timeset(routineN, handle) 859 860 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints) 861 862 ! Mulliken charges require no further computation and are printed from write_mo_free_results 863 864 ! Compute the Lowdin charges 865 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%LOWDIN") 866 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 867 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%LOWDIN", extension=".lowdin", & 868 log_filename=.FALSE.) 869 print_level = 1 870 CALL section_vals_val_get(print_key, "PRINT_GOP", l_val=print_it) 871 IF (print_it) print_level = 2 872 CALL section_vals_val_get(print_key, "PRINT_ALL", l_val=print_it) 873 IF (print_it) print_level = 3 874 IF (do_kpoints) THEN 875 CPWARN("Lowdin charges not implemented for k-point calculations!") 876 ELSE 877 CALL lowdin_population_analysis(qs_env, unit_nr, print_level) 878 END IF 879 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%LOWDIN") 880 END IF 881 882 ! Compute the RESP charges 883 CALL resp_fit(qs_env) 884 885 ! Compute the Density Derived Atomic Point charges with the Bloechl scheme 886 print_key => section_vals_get_subs_vals(input, "PROPERTIES%FIT_CHARGE") 887 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 888 unit_nr = cp_print_key_unit_nr(logger, input, "PROPERTIES%FIT_CHARGE", extension=".Fitcharge", & 889 log_filename=.FALSE.) 890 density_fit_section => section_vals_get_subs_vals(input, "DFT%DENSITY_FITTING") 891 CALL get_ddapc(qs_env, .FALSE., density_fit_section, iwc=unit_nr) 892 CALL cp_print_key_finished_output(unit_nr, logger, input, "PROPERTIES%FIT_CHARGE") 893 END IF 894 895 CALL timestop(handle) 896 897 END SUBROUTINE qs_scf_post_charges 898 899! ************************************************************************************************** 900!> \brief Computes and prints the Cube Files for MO 901!> \param input ... 902!> \param dft_section ... 903!> \param dft_control ... 904!> \param logger ... 905!> \param qs_env the qs_env in which the qs_env lives 906!> \param mo_coeff ... 907!> \param wf_g ... 908!> \param wf_r ... 909!> \param particles ... 910!> \param homo ... 911!> \param ispin ... 912! ************************************************************************************************** 913 SUBROUTINE qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env, & 914 mo_coeff, wf_g, wf_r, particles, homo, ispin) 915 TYPE(section_vals_type), POINTER :: input, dft_section 916 TYPE(dft_control_type), POINTER :: dft_control 917 TYPE(cp_logger_type), POINTER :: logger 918 TYPE(qs_environment_type), POINTER :: qs_env 919 TYPE(cp_fm_type), POINTER :: mo_coeff 920 TYPE(pw_p_type) :: wf_g, wf_r 921 TYPE(particle_list_type), POINTER :: particles 922 INTEGER, INTENT(IN) :: homo, ispin 923 924 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_occ_cubes', & 925 routineP = moduleN//':'//routineN 926 927 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title 928 INTEGER :: handle, i, ir, ivector, n_rep, nhomo, & 929 nlist, unit_nr 930 INTEGER, DIMENSION(:), POINTER :: list, list_index 931 LOGICAL :: append_cube, mpi_io 932 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 933 TYPE(cell_type), POINTER :: cell 934 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 935 TYPE(pw_env_type), POINTER :: pw_env 936 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 937 938 CALL timeset(routineN, handle) 939 940 NULLIFY (list_index) 941 942 IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_CUBES") & 943 , cp_p_file) .AND. section_get_lval(dft_section, "PRINT%MO_CUBES%WRITE_CUBE")) THEN 944 nhomo = section_get_ival(dft_section, "PRINT%MO_CUBES%NHOMO") 945 append_cube = section_get_lval(dft_section, "PRINT%MO_CUBES%APPEND") 946 my_pos_cube = "REWIND" 947 IF (append_cube) THEN 948 my_pos_cube = "APPEND" 949 END IF 950 CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%HOMO_LIST", n_rep_val=n_rep) 951 IF (n_rep > 0) THEN ! write the cubes of the list 952 nlist = 0 953 DO ir = 1, n_rep 954 NULLIFY (list) 955 CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%HOMO_LIST", i_rep_val=ir, & 956 i_vals=list) 957 IF (ASSOCIATED(list)) THEN 958 CALL reallocate(list_index, 1, nlist + SIZE(list)) 959 DO i = 1, SIZE(list) 960 list_index(i + nlist) = list(i) 961 END DO 962 nlist = nlist + SIZE(list) 963 END IF 964 END DO 965 ELSE 966 967 IF (nhomo == -1) nhomo = homo 968 nlist = homo - MAX(1, homo - nhomo + 1) + 1 969 ALLOCATE (list_index(nlist)) 970 DO i = 1, nlist 971 list_index(i) = MAX(1, homo - nhomo + 1) + i - 1 972 END DO 973 END IF 974 DO i = 1, nlist 975 ivector = list_index(i) 976 CALL get_qs_env(qs_env=qs_env, & 977 atomic_kind_set=atomic_kind_set, & 978 qs_kind_set=qs_kind_set, & 979 cell=cell, & 980 particle_set=particle_set, & 981 pw_env=pw_env) 982 CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, qs_kind_set, & 983 cell, dft_control, particle_set, pw_env) 984 WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", ivector, "_", ispin 985 mpi_io = .TRUE. 986 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MO_CUBES", extension=".cube", & 987 middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., & 988 mpi_io=mpi_io) 989 WRITE (title, *) "WAVEFUNCTION ", ivector, " spin ", ispin, " i.e. HOMO - ", ivector - homo 990 CALL cp_pw_to_cube(wf_r%pw, unit_nr, title, particles=particles, & 991 stride=section_get_ivals(dft_section, "PRINT%MO_CUBES%STRIDE"), mpi_io=mpi_io) 992 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MO_CUBES", mpi_io=mpi_io) 993 ENDDO 994 IF (ASSOCIATED(list_index)) DEALLOCATE (list_index) 995 END IF 996 997 CALL timestop(handle) 998 999 END SUBROUTINE qs_scf_post_occ_cubes 1000 1001! ************************************************************************************************** 1002!> \brief Computes and prints the Cube Files for MO 1003!> \param input ... 1004!> \param dft_section ... 1005!> \param dft_control ... 1006!> \param logger ... 1007!> \param qs_env the qs_env in which the qs_env lives 1008!> \param unoccupied_orbs ... 1009!> \param wf_g ... 1010!> \param wf_r ... 1011!> \param particles ... 1012!> \param nlumos ... 1013!> \param homo ... 1014!> \param ispin ... 1015!> \param lumo ... 1016! ************************************************************************************************** 1017 SUBROUTINE qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, & 1018 unoccupied_orbs, wf_g, wf_r, particles, nlumos, homo, ispin, lumo) 1019 TYPE(section_vals_type), POINTER :: input, dft_section 1020 TYPE(dft_control_type), POINTER :: dft_control 1021 TYPE(cp_logger_type), POINTER :: logger 1022 TYPE(qs_environment_type), POINTER :: qs_env 1023 TYPE(cp_fm_type), POINTER :: unoccupied_orbs 1024 TYPE(pw_p_type) :: wf_g, wf_r 1025 TYPE(particle_list_type), POINTER :: particles 1026 INTEGER, INTENT(IN) :: nlumos, homo, ispin 1027 INTEGER, INTENT(IN), OPTIONAL :: lumo 1028 1029 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_unocc_cubes', & 1030 routineP = moduleN//':'//routineN 1031 1032 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title 1033 INTEGER :: handle, ifirst, index_mo, ivector, & 1034 unit_nr 1035 LOGICAL :: append_cube, mpi_io 1036 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1037 TYPE(cell_type), POINTER :: cell 1038 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1039 TYPE(pw_env_type), POINTER :: pw_env 1040 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1041 1042 CALL timeset(routineN, handle) 1043 1044 IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_CUBES"), cp_p_file) & 1045 .AND. section_get_lval(dft_section, "PRINT%MO_CUBES%WRITE_CUBE")) THEN 1046 NULLIFY (qs_kind_set, particle_set, pw_env, cell) 1047 append_cube = section_get_lval(dft_section, "PRINT%MO_CUBES%APPEND") 1048 my_pos_cube = "REWIND" 1049 IF (append_cube) THEN 1050 my_pos_cube = "APPEND" 1051 END IF 1052 ifirst = 1 1053 IF (PRESENT(lumo)) ifirst = lumo 1054 DO ivector = ifirst, ifirst + nlumos - 1 1055 CALL get_qs_env(qs_env=qs_env, & 1056 atomic_kind_set=atomic_kind_set, & 1057 qs_kind_set=qs_kind_set, & 1058 cell=cell, & 1059 particle_set=particle_set, & 1060 pw_env=pw_env) 1061 CALL calculate_wavefunction(unoccupied_orbs, ivector, wf_r, wf_g, atomic_kind_set, & 1062 qs_kind_set, cell, dft_control, particle_set, pw_env) 1063 1064 IF (ifirst == 1) THEN 1065 index_mo = homo + ivector 1066 ELSE 1067 index_mo = ivector 1068 END IF 1069 WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", index_mo, "_", ispin 1070 mpi_io = .TRUE. 1071 1072 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MO_CUBES", extension=".cube", & 1073 middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., & 1074 mpi_io=mpi_io) 1075 WRITE (title, *) "WAVEFUNCTION ", index_mo, " spin ", ispin, " i.e. LUMO + ", ifirst + ivector - 2 1076 CALL cp_pw_to_cube(wf_r%pw, unit_nr, title, particles=particles, & 1077 stride=section_get_ivals(dft_section, "PRINT%MO_CUBES%STRIDE"), mpi_io=mpi_io) 1078 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MO_CUBES", mpi_io=mpi_io) 1079 ENDDO 1080 ENDIF 1081 1082 CALL timestop(handle) 1083 1084 END SUBROUTINE qs_scf_post_unocc_cubes 1085 1086! ************************************************************************************************** 1087!> \brief Computes and prints electric moments 1088!> \param input ... 1089!> \param logger ... 1090!> \param qs_env the qs_env in which the qs_env lives 1091!> \param output_unit ... 1092! ************************************************************************************************** 1093 SUBROUTINE qs_scf_post_moments(input, logger, qs_env, output_unit) 1094 TYPE(section_vals_type), POINTER :: input 1095 TYPE(cp_logger_type), POINTER :: logger 1096 TYPE(qs_environment_type), POINTER :: qs_env 1097 INTEGER, INTENT(IN) :: output_unit 1098 1099 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_moments', & 1100 routineP = moduleN//':'//routineN 1101 1102 CHARACTER(LEN=default_path_length) :: filename 1103 INTEGER :: handle, maxmom, reference, unit_nr 1104 LOGICAL :: do_kpoints, magnetic, periodic 1105 REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point 1106 TYPE(section_vals_type), POINTER :: print_key 1107 1108 CALL timeset(routineN, handle) 1109 1110 print_key => section_vals_get_subs_vals(section_vals=input, & 1111 subsection_name="DFT%PRINT%MOMENTS") 1112 1113 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 1114 1115 maxmom = section_get_ival(section_vals=input, & 1116 keyword_name="DFT%PRINT%MOMENTS%MAX_MOMENT") 1117 periodic = section_get_lval(section_vals=input, & 1118 keyword_name="DFT%PRINT%MOMENTS%PERIODIC") 1119 reference = section_get_ival(section_vals=input, & 1120 keyword_name="DFT%PRINT%MOMENTS%REFERENCE") 1121 magnetic = section_get_lval(section_vals=input, & 1122 keyword_name="DFT%PRINT%MOMENTS%MAGNETIC") 1123 NULLIFY (ref_point) 1124 CALL section_vals_val_get(input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point) 1125 unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=input, & 1126 print_key_path="DFT%PRINT%MOMENTS", extension=".dat", & 1127 middle_name="moments", log_filename=.FALSE.) 1128 1129 IF (output_unit > 0) THEN 1130 IF (unit_nr /= output_unit) THEN 1131 INQUIRE (UNIT=unit_nr, NAME=filename) 1132 WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") & 1133 "MOMENTS", "The electric/magnetic moments are written to file:", & 1134 TRIM(filename) 1135 ELSE 1136 WRITE (UNIT=output_unit, FMT="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS" 1137 END IF 1138 END IF 1139 1140 CALL get_qs_env(qs_env, do_kpoints=do_kpoints) 1141 1142 IF (do_kpoints) THEN 1143 CPWARN("Moments not implemented for k-point calculations!") 1144 ELSE 1145 IF (periodic) THEN 1146 CALL qs_moment_berry_phase(qs_env, magnetic, maxmom, reference, ref_point, unit_nr) 1147 ELSE 1148 CALL qs_moment_locop(qs_env, magnetic, maxmom, reference, ref_point, unit_nr) 1149 END IF 1150 END IF 1151 1152 CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, & 1153 basis_section=input, print_key_path="DFT%PRINT%MOMENTS") 1154 END IF 1155 1156 CALL timestop(handle) 1157 1158 END SUBROUTINE qs_scf_post_moments 1159 1160! ************************************************************************************************** 1161!> \brief Computes and prints the X-ray diffraction spectrum. 1162!> \param input ... 1163!> \param dft_section ... 1164!> \param logger ... 1165!> \param qs_env the qs_env in which the qs_env lives 1166!> \param output_unit ... 1167! ************************************************************************************************** 1168 SUBROUTINE qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit) 1169 1170 TYPE(section_vals_type), POINTER :: input, dft_section 1171 TYPE(cp_logger_type), POINTER :: logger 1172 TYPE(qs_environment_type), POINTER :: qs_env 1173 INTEGER, INTENT(IN) :: output_unit 1174 1175 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_post_xray', & 1176 routineP = moduleN//':'//routineN 1177 1178 CHARACTER(LEN=default_path_length) :: filename 1179 INTEGER :: handle, unit_nr 1180 REAL(KIND=dp) :: q_max 1181 TYPE(section_vals_type), POINTER :: print_key 1182 1183 CALL timeset(routineN, handle) 1184 1185 print_key => section_vals_get_subs_vals(section_vals=input, & 1186 subsection_name="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM") 1187 1188 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 1189 q_max = section_get_rval(section_vals=dft_section, & 1190 keyword_name="PRINT%XRAY_DIFFRACTION_SPECTRUM%Q_MAX") 1191 unit_nr = cp_print_key_unit_nr(logger=logger, & 1192 basis_section=input, & 1193 print_key_path="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM", & 1194 extension=".dat", & 1195 middle_name="xrd", & 1196 log_filename=.FALSE.) 1197 IF (output_unit > 0) THEN 1198 INQUIRE (UNIT=unit_nr, NAME=filename) 1199 WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") & 1200 "X-RAY DIFFRACTION SPECTRUM" 1201 IF (unit_nr /= output_unit) THEN 1202 WRITE (UNIT=output_unit, FMT="(/,T3,A,/,/,T3,A,/)") & 1203 "The coherent X-ray diffraction spectrum is written to the file:", & 1204 TRIM(filename) 1205 END IF 1206 END IF 1207 CALL xray_diffraction_spectrum(qs_env=qs_env, & 1208 unit_number=unit_nr, & 1209 q_max=q_max) 1210 CALL cp_print_key_finished_output(unit_nr=unit_nr, & 1211 logger=logger, & 1212 basis_section=input, & 1213 print_key_path="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM") 1214 END IF 1215 1216 CALL timestop(handle) 1217 1218 END SUBROUTINE qs_scf_post_xray 1219 1220! ************************************************************************************************** 1221!> \brief Computes and prints Electric Field Gradient 1222!> \param input ... 1223!> \param logger ... 1224!> \param qs_env the qs_env in which the qs_env lives 1225! ************************************************************************************************** 1226 SUBROUTINE qs_scf_post_efg(input, logger, qs_env) 1227 TYPE(section_vals_type), POINTER :: input 1228 TYPE(cp_logger_type), POINTER :: logger 1229 TYPE(qs_environment_type), POINTER :: qs_env 1230 1231 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_efg', & 1232 routineP = moduleN//':'//routineN 1233 1234 INTEGER :: handle 1235 TYPE(section_vals_type), POINTER :: print_key 1236 1237 CALL timeset(routineN, handle) 1238 1239 print_key => section_vals_get_subs_vals(section_vals=input, & 1240 subsection_name="DFT%PRINT%ELECTRIC_FIELD_GRADIENT") 1241 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), & 1242 cp_p_file)) THEN 1243 CALL qs_efg_calc(qs_env=qs_env) 1244 END IF 1245 1246 CALL timestop(handle) 1247 1248 END SUBROUTINE qs_scf_post_efg 1249 1250! ************************************************************************************************** 1251!> \brief Computes the Electron Transfer Coupling matrix element 1252!> \param input ... 1253!> \param qs_env the qs_env in which the qs_env lives 1254!> \param dft_control ... 1255! ************************************************************************************************** 1256 SUBROUTINE qs_scf_post_et(input, qs_env, dft_control) 1257 TYPE(section_vals_type), POINTER :: input 1258 TYPE(qs_environment_type), POINTER :: qs_env 1259 TYPE(dft_control_type), POINTER :: dft_control 1260 1261 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_et', routineP = moduleN//':'//routineN 1262 1263 INTEGER :: handle, ispin 1264 LOGICAL :: do_et 1265 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: my_mos 1266 TYPE(section_vals_type), POINTER :: et_section 1267 1268 CALL timeset(routineN, handle) 1269 1270 do_et = .FALSE. 1271 et_section => section_vals_get_subs_vals(input, "PROPERTIES%ET_COUPLING") 1272 CALL section_vals_get(et_section, explicit=do_et) 1273 IF (do_et) THEN 1274 IF (qs_env%et_coupling%first_run) THEN 1275 NULLIFY (my_mos) 1276 ALLOCATE (my_mos(dft_control%nspins)) 1277 ALLOCATE (qs_env%et_coupling%et_mo_coeff(dft_control%nspins)) 1278 DO ispin = 1, dft_control%nspins 1279 NULLIFY (my_mos(ispin)%matrix) 1280 CALL cp_fm_create(matrix=my_mos(ispin)%matrix, & 1281 matrix_struct=qs_env%mos(ispin)%mo_set%mo_coeff%matrix_struct, & 1282 name="FIRST_RUN_COEFF"//TRIM(ADJUSTL(cp_to_string(ispin)))//"MATRIX") 1283 CALL cp_fm_to_fm(qs_env%mos(ispin)%mo_set%mo_coeff, & 1284 my_mos(ispin)%matrix) 1285 END DO 1286 CALL set_et_coupling_type(qs_env%et_coupling, et_mo_coeff=my_mos) 1287 DEALLOCATE (my_mos) 1288 END IF 1289 END IF 1290 1291 CALL timestop(handle) 1292 1293 END SUBROUTINE qs_scf_post_et 1294 1295! ************************************************************************************************** 1296!> \brief compute the electron localization function 1297!> 1298!> \param input ... 1299!> \param logger ... 1300!> \param qs_env ... 1301!> \par History 1302!> 2012-07 Created [MI] 1303! ************************************************************************************************** 1304 SUBROUTINE qs_scf_post_elf(input, logger, qs_env) 1305 TYPE(section_vals_type), POINTER :: input 1306 TYPE(cp_logger_type), POINTER :: logger 1307 TYPE(qs_environment_type), POINTER :: qs_env 1308 1309 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_elf', & 1310 routineP = moduleN//':'//routineN 1311 1312 CHARACTER(LEN=default_path_length) :: filename, mpi_filename, my_pos_cube, & 1313 title 1314 INTEGER :: handle, ispin, output_unit, unit_nr 1315 LOGICAL :: append_cube, gapw, mpi_io 1316 REAL(dp) :: rho_cutoff 1317 TYPE(dft_control_type), POINTER :: dft_control 1318 TYPE(particle_list_type), POINTER :: particles 1319 TYPE(pw_env_type), POINTER :: pw_env 1320 TYPE(pw_p_type), DIMENSION(:), POINTER :: elf_r 1321 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 1322 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 1323 TYPE(qs_subsys_type), POINTER :: subsys 1324 TYPE(section_vals_type), POINTER :: elf_section 1325 1326 CALL timeset(routineN, handle) 1327 output_unit = cp_logger_get_default_io_unit(logger) 1328 1329 elf_section => section_vals_get_subs_vals(input, "DFT%PRINT%ELF_CUBE") 1330 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 1331 "DFT%PRINT%ELF_CUBE"), cp_p_file)) THEN 1332 1333 NULLIFY (dft_control, pw_env, auxbas_pw_pool, pw_pools, particles, subsys, elf_r) 1334 CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, subsys=subsys) 1335 CALL qs_subsys_get(subsys, particles=particles) 1336 1337 gapw = dft_control%qs_control%gapw 1338 IF (.NOT. gapw) THEN 1339 ! allocate 1340 ALLOCATE (elf_r(dft_control%nspins)) 1341 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 1342 pw_pools=pw_pools) 1343 DO ispin = 1, dft_control%nspins 1344 CALL pw_pool_create_pw(auxbas_pw_pool, elf_r(ispin)%pw, & 1345 use_data=REALDATA3D, & 1346 in_space=REALSPACE) 1347 CALL pw_zero(elf_r(ispin)%pw) 1348 END DO 1349 1350 IF (output_unit > 0) THEN 1351 WRITE (UNIT=output_unit, FMT="(/,T15,A,/)") & 1352 " ----- ELF is computed on the real space grid -----" 1353 END IF 1354 rho_cutoff = section_get_rval(elf_section, "density_cutoff") 1355 CALL qs_elf_calc(qs_env, elf_r, rho_cutoff) 1356 1357 ! write ELF into cube file 1358 append_cube = section_get_lval(elf_section, "APPEND") 1359 my_pos_cube = "REWIND" 1360 IF (append_cube) THEN 1361 my_pos_cube = "APPEND" 1362 END IF 1363 1364 DO ispin = 1, dft_control%nspins 1365 WRITE (filename, '(a5,I1.1)') "ELF_S", ispin 1366 WRITE (title, *) "ELF spin ", ispin 1367 mpi_io = .TRUE. 1368 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%ELF_CUBE", extension=".cube", & 1369 middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., & 1370 mpi_io=mpi_io, fout=mpi_filename) 1371 IF (output_unit > 0) THEN 1372 IF (.NOT. mpi_io) THEN 1373 INQUIRE (UNIT=unit_nr, NAME=filename) 1374 ELSE 1375 filename = mpi_filename 1376 END IF 1377 WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & 1378 "ELF is written in cube file format to the file:", & 1379 TRIM(filename) 1380 END IF 1381 1382 CALL cp_pw_to_cube(elf_r(ispin)%pw, unit_nr, title, particles=particles, & 1383 stride=section_get_ivals(elf_section, "STRIDE"), mpi_io=mpi_io) 1384 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%ELF_CUBE", mpi_io=mpi_io) 1385 1386 CALL pw_pool_give_back_pw(auxbas_pw_pool, elf_r(ispin)%pw) 1387 END DO 1388 1389 ! deallocate 1390 DEALLOCATE (elf_r) 1391 1392 ELSE 1393 ! not implemented 1394 CPWARN("ELF not implemented for GAPW calculations!!") 1395 1396 END IF 1397 1398 END IF ! print key 1399 1400 CALL timestop(handle) 1401 1402 END SUBROUTINE qs_scf_post_elf 1403 1404! ************************************************************************************************** 1405!> \brief computes the condition number of the overlap matrix and 1406!> prints the value of the total energy. This is needed 1407!> for BASIS_MOLOPT optimizations 1408!> \param input ... 1409!> \param logger ... 1410!> \param qs_env the qs_env in which the qs_env lives 1411!> \par History 1412!> 2007-07 Created [Joost VandeVondele] 1413! ************************************************************************************************** 1414 SUBROUTINE qs_scf_post_molopt(input, logger, qs_env) 1415 TYPE(section_vals_type), POINTER :: input 1416 TYPE(cp_logger_type), POINTER :: logger 1417 TYPE(qs_environment_type), POINTER :: qs_env 1418 1419 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_molopt', & 1420 routineP = moduleN//':'//routineN 1421 1422 INTEGER :: handle, nao, unit_nr 1423 REAL(KIND=dp) :: S_cond_number 1424 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues 1425 TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct 1426 TYPE(cp_fm_type), POINTER :: fm_s, fm_work, mo_coeff 1427 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 1428 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 1429 TYPE(qs_energy_type), POINTER :: energy 1430 TYPE(section_vals_type), POINTER :: print_key 1431 1432 CALL timeset(routineN, handle) 1433 1434 print_key => section_vals_get_subs_vals(section_vals=input, & 1435 subsection_name="DFT%PRINT%BASIS_MOLOPT_QUANTITIES") 1436 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), & 1437 cp_p_file)) THEN 1438 1439 CALL get_qs_env(qs_env, energy=energy, matrix_s=matrix_s, mos=mos) 1440 1441 ! set up the two needed full matrices, using mo_coeff as a template 1442 CALL get_mo_set(mo_set=mos(1)%mo_set, mo_coeff=mo_coeff, nao=nao) 1443 CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, & 1444 nrow_global=nao, ncol_global=nao, & 1445 template_fmstruct=mo_coeff%matrix_struct) 1446 CALL cp_fm_create(fm_s, matrix_struct=ao_ao_fmstruct, & 1447 name="fm_s") 1448 CALL cp_fm_create(fm_work, matrix_struct=ao_ao_fmstruct, & 1449 name="fm_work") 1450 CALL cp_fm_struct_release(ao_ao_fmstruct) 1451 ALLOCATE (eigenvalues(nao)) 1452 1453 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, fm_s) 1454 CALL choose_eigv_solver(fm_s, fm_work, eigenvalues) 1455 1456 CALL cp_fm_release(fm_s) 1457 CALL cp_fm_release(fm_work) 1458 1459 S_cond_number = MAXVAL(ABS(eigenvalues))/MAX(MINVAL(ABS(eigenvalues)), EPSILON(0.0_dp)) 1460 1461 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%BASIS_MOLOPT_QUANTITIES", & 1462 extension=".molopt") 1463 1464 IF (unit_nr > 0) THEN 1465 ! please keep this format fixed, needs to be grepable for molopt 1466 ! optimizations 1467 WRITE (unit_nr, '(T2,A28,2A25)') "", "Tot. Ener.", "S Cond. Numb." 1468 WRITE (unit_nr, '(T2,A28,2E25.17)') "BASIS_MOLOPT_QUANTITIES", energy%total, S_cond_number 1469 ENDIF 1470 1471 CALL cp_print_key_finished_output(unit_nr, logger, input, & 1472 "DFT%PRINT%BASIS_MOLOPT_QUANTITIES") 1473 1474 END IF 1475 1476 CALL timestop(handle) 1477 1478 END SUBROUTINE qs_scf_post_molopt 1479 1480! ************************************************************************************************** 1481!> \brief Dumps EPR 1482!> \param input ... 1483!> \param logger ... 1484!> \param qs_env the qs_env in which the qs_env lives 1485! ************************************************************************************************** 1486 SUBROUTINE qs_scf_post_epr(input, logger, qs_env) 1487 TYPE(section_vals_type), POINTER :: input 1488 TYPE(cp_logger_type), POINTER :: logger 1489 TYPE(qs_environment_type), POINTER :: qs_env 1490 1491 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_epr', & 1492 routineP = moduleN//':'//routineN 1493 1494 INTEGER :: handle 1495 TYPE(section_vals_type), POINTER :: print_key 1496 1497 CALL timeset(routineN, handle) 1498 1499 print_key => section_vals_get_subs_vals(section_vals=input, & 1500 subsection_name="DFT%PRINT%HYPERFINE_COUPLING_TENSOR") 1501 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), & 1502 cp_p_file)) THEN 1503 CALL qs_epr_hyp_calc(qs_env=qs_env) 1504 END IF 1505 1506 CALL timestop(handle) 1507 1508 END SUBROUTINE qs_scf_post_epr 1509 1510! ************************************************************************************************** 1511!> \brief Interface routine to trigger writing of results available from normal 1512!> SCF. Can write MO-dependent and MO free results (needed for call from 1513!> the linear scaling code) 1514!> \param qs_env the qs_env in which the qs_env lives 1515!> \param scf_env ... 1516! ************************************************************************************************** 1517 SUBROUTINE write_available_results(qs_env, scf_env) 1518 TYPE(qs_environment_type), POINTER :: qs_env 1519 TYPE(qs_scf_env_type), OPTIONAL, POINTER :: scf_env 1520 1521 CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results', & 1522 routineP = moduleN//':'//routineN 1523 1524 INTEGER :: handle 1525 1526 CALL timeset(routineN, handle) 1527 1528 ! those properties that require MOs (not suitable density matrix based methods) 1529 CALL write_mo_dependent_results(qs_env, scf_env) 1530 1531 ! those that depend only on the density matrix, they should be linear scaling in their implementation 1532 CALL write_mo_free_results(qs_env) 1533 1534 CALL timestop(handle) 1535 1536 END SUBROUTINE write_available_results 1537 1538! ************************************************************************************************** 1539!> \brief Write QS results available if MO's are present (if switched on through the print_keys) 1540!> Writes only MO dependent results. Split is necessary as ls_scf does not 1541!> provide MO's 1542!> \param qs_env the qs_env in which the qs_env lives 1543!> \param scf_env ... 1544! ************************************************************************************************** 1545 SUBROUTINE write_mo_dependent_results(qs_env, scf_env) 1546 TYPE(qs_environment_type), POINTER :: qs_env 1547 TYPE(qs_scf_env_type), OPTIONAL, POINTER :: scf_env 1548 1549 CHARACTER(len=*), PARAMETER :: routineN = 'write_mo_dependent_results', & 1550 routineP = moduleN//':'//routineN 1551 1552 INTEGER :: handle, homo, ik, ispin, nmo, output_unit 1553 LOGICAL :: all_equal, do_kpoints 1554 REAL(KIND=dp) :: maxocc, s_square, s_square_ideal, & 1555 total_abs_spin_dens 1556 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, occupation_numbers 1557 TYPE(admm_type), POINTER :: admm_env 1558 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1559 TYPE(cell_type), POINTER :: cell 1560 TYPE(cp_fm_type), POINTER :: mo_coeff 1561 TYPE(cp_logger_type), POINTER :: logger 1562 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, matrix_s 1563 TYPE(dbcsr_type), POINTER :: mo_coeff_deriv 1564 TYPE(dft_control_type), POINTER :: dft_control 1565 TYPE(kpoint_type), POINTER :: kpoints 1566 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 1567 TYPE(mo_set_p_type), DIMENSION(:, :), POINTER :: mos_k 1568 TYPE(molecule_type), POINTER :: molecule_set(:) 1569 TYPE(particle_list_type), POINTER :: particles 1570 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1571 TYPE(pw_env_type), POINTER :: pw_env 1572 TYPE(pw_p_type) :: wf_r 1573 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r 1574 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 1575 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 1576 TYPE(qs_energy_type), POINTER :: energy 1577 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1578 TYPE(qs_rho_type), POINTER :: rho 1579 TYPE(qs_subsys_type), POINTER :: subsys 1580 TYPE(scf_control_type), POINTER :: scf_control 1581 TYPE(section_vals_type), POINTER :: dft_section, input, sprint_section 1582 1583 CALL timeset(routineN, handle) 1584 NULLIFY (cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, mo_coeff, & 1585 mo_coeff_deriv, mo_eigenvalues, mos, atomic_kind_set, qs_kind_set, & 1586 particle_set, rho, ks_rmpv, matrix_s, scf_control, dft_section, & 1587 molecule_set, input, particles, subsys, rho_r) 1588 1589 logger => cp_get_default_logger() 1590 output_unit = cp_logger_get_default_io_unit(logger) 1591 1592 CPASSERT(ASSOCIATED(qs_env)) 1593 CALL get_qs_env(qs_env, & 1594 dft_control=dft_control, & 1595 molecule_set=molecule_set, & 1596 atomic_kind_set=atomic_kind_set, & 1597 particle_set=particle_set, & 1598 qs_kind_set=qs_kind_set, & 1599 admm_env=admm_env, & 1600 scf_control=scf_control, & 1601 input=input, & 1602 cell=cell, & 1603 subsys=subsys) 1604 CALL qs_subsys_get(subsys, particles=particles) 1605 CALL get_qs_env(qs_env, rho=rho) 1606 CALL qs_rho_get(rho, rho_r=rho_r) 1607 1608 ! kpoints 1609 CALL get_qs_env(qs_env, do_kpoints=do_kpoints) 1610 1611 ! *** if the dft_section tells you to do so, write last wavefunction to screen 1612 dft_section => section_vals_get_subs_vals(input, "DFT") 1613 IF (.NOT. qs_env%run_rtp) THEN 1614 IF (.NOT. do_kpoints) THEN 1615 CALL get_qs_env(qs_env, mos=mos) 1616 IF (dft_control%nspins == 2) THEN 1617 CALL write_mo_set(mos(1)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 1618 dft_section, spin="ALPHA", last=.TRUE.) 1619 CALL write_mo_set(mos(2)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 1620 dft_section, spin="BETA", last=.TRUE.) 1621 ELSE 1622 CALL write_mo_set(mos(1)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 1623 dft_section, last=.TRUE.) 1624 END IF 1625 CALL get_qs_env(qs_env, matrix_ks=ks_rmpv) 1626 CALL write_dm_binary_restart(mos, dft_section, ks_rmpv) 1627 sprint_section => section_vals_get_subs_vals(dft_section, "PRINT%MO_MOLDEN") 1628 CALL write_mos_molden(mos, qs_kind_set, particle_set, sprint_section) 1629 ELSE ! *** if we have k points, let's print the eigenvalues at each of them 1630 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints) 1631 1632 IF (kpoints%nkp /= 0) THEN 1633 DO ik = 1, SIZE(kpoints%kp_env) 1634 mos_k => kpoints%kp_env(ik)%kpoint_env%mos 1635 CPASSERT(ASSOCIATED(mos_k)) 1636 1637 IF (dft_control%nspins == 2) THEN 1638 CALL write_mo_set(mos_k(1, 1)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 1639 dft_section, spin="ALPHA", last=.TRUE., kpt=ik) 1640 CALL write_mo_set(mos_k(1, 2)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 1641 dft_section, spin="BETA", last=.TRUE., kpt=ik) 1642 ELSE 1643 CALL write_mo_set(mos_k(1, 1)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 1644 dft_section, last=.TRUE., kpt=ik) 1645 END IF 1646 END DO 1647 END IF 1648 END IF 1649 1650 ! *** at the end of scf print out the projected dos per kind 1651 IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%PDOS") & 1652 , cp_p_file)) THEN 1653 IF (do_kpoints) THEN 1654 CPWARN("Projected density of states not implemented for k-points.") 1655 ELSE 1656 CALL get_qs_env(qs_env, & 1657 mos=mos, & 1658 matrix_ks=ks_rmpv) 1659 DO ispin = 1, dft_control%nspins 1660 1661 ! ** If we do ADMM, we add have to modify the kohn-sham matrix 1662 IF (dft_control%do_admm) THEN 1663 CALL admm_correct_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix) 1664 END IF 1665 IF (PRESENT(scf_env)) THEN 1666 IF (scf_env%method == ot_method_nr) THEN 1667 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, & 1668 eigenvalues=mo_eigenvalues) 1669 IF (ASSOCIATED(qs_env%mo_derivs)) THEN 1670 mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix 1671 ELSE 1672 mo_coeff_deriv => NULL() 1673 ENDIF 1674 1675 CALL calculate_subspace_eigenvalues(mo_coeff, ks_rmpv(ispin)%matrix, mo_eigenvalues, & 1676 do_rotation=.TRUE., & 1677 co_rotate_dbcsr=mo_coeff_deriv) 1678 CALL set_mo_occupation(mo_set=mos(ispin)%mo_set) 1679 1680 END IF 1681 END IF 1682 IF (dft_control%nspins == 2) THEN 1683 CALL calculate_projected_dos(mos(ispin)%mo_set, atomic_kind_set, & 1684 qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin) 1685 ELSE 1686 CALL calculate_projected_dos(mos(ispin)%mo_set, atomic_kind_set, & 1687 qs_kind_set, particle_set, qs_env, dft_section) 1688 END IF 1689 1690 ! ** If we do ADMM, we add have to modify the kohn-sham matrix 1691 IF (dft_control%do_admm) THEN 1692 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix) 1693 END IF 1694 1695 END DO 1696 ENDIF 1697 ENDIF 1698 END IF 1699 1700 ! *** Integrated absolute spin density and spin contamination *** 1701 IF (dft_control%nspins .EQ. 2) THEN 1702 CALL get_qs_env(qs_env, mos=mos) 1703 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 1704 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 1705 pw_pools=pw_pools) 1706 CALL pw_pool_create_pw(auxbas_pw_pool, wf_r%pw, & 1707 use_data=REALDATA3D, & 1708 in_space=REALSPACE) 1709 CALL pw_copy(rho_r(1)%pw, wf_r%pw) 1710 CALL pw_axpy(rho_r(2)%pw, wf_r%pw, alpha=-1._dp) 1711 total_abs_spin_dens = pw_integrate_function(wf_r%pw, oprt="ABS") 1712 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,F20.10))') & 1713 "Integrated absolute spin density : ", total_abs_spin_dens 1714 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_r%pw) 1715 ! 1716 ! XXX Fix Me XXX 1717 ! should be extended to the case where added MOs are present 1718 ! should be extended to the k-point case 1719 ! 1720 IF (do_kpoints) THEN 1721 CPWARN("Spin contamination estimate not implemented for k-points.") 1722 ELSE 1723 all_equal = .TRUE. 1724 DO ispin = 1, dft_control%nspins 1725 CALL get_mo_set(mo_set=mos(ispin)%mo_set, & 1726 occupation_numbers=occupation_numbers, & 1727 homo=homo, & 1728 nmo=nmo, & 1729 maxocc=maxocc) 1730 IF (nmo > 0) THEN 1731 all_equal = all_equal .AND. & 1732 (ALL(occupation_numbers(1:homo) == maxocc) .AND. & 1733 ALL(occupation_numbers(homo + 1:nmo) == 0.0_dp)) 1734 END IF 1735 END DO 1736 IF (.NOT. all_equal) THEN 1737 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A)") & 1738 "WARNING: S**2 computation does not yet treat fractional occupied orbitals" 1739 ELSE 1740 CALL get_qs_env(qs_env=qs_env, & 1741 matrix_s=matrix_s, & 1742 energy=energy) 1743 CALL compute_s_square(mos=mos, matrix_s=matrix_s, s_square=s_square, & 1744 s_square_ideal=s_square_ideal) 1745 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(T3,A,T51,2F15.6)') & 1746 "Ideal and single determinant S**2 : ", s_square_ideal, s_square 1747 energy%s_square = s_square 1748 END IF 1749 ENDIF 1750 ENDIF 1751 1752 CALL timestop(handle) 1753 1754 END SUBROUTINE write_mo_dependent_results 1755 1756! ************************************************************************************************** 1757!> \brief Write QS results always available (if switched on through the print_keys) 1758!> Can be called from ls_scf 1759!> \param qs_env the qs_env in which the qs_env lives 1760! ************************************************************************************************** 1761 SUBROUTINE write_mo_free_results(qs_env) 1762 TYPE(qs_environment_type), POINTER :: qs_env 1763 1764 CHARACTER(len=*), PARAMETER :: routineN = 'write_mo_free_results', & 1765 routineP = moduleN//':'//routineN 1766 CHARACTER(len=1), DIMENSION(3), PARAMETER :: cdir = (/"x", "y", "z"/) 1767 1768 CHARACTER(LEN=2) :: element_symbol 1769 CHARACTER(LEN=default_path_length) :: filename, mpi_filename, my_pos_cube 1770 CHARACTER(LEN=default_string_length) :: name 1771 INTEGER :: after, handle, i, iat, id, ikind, img, & 1772 iso, ispin, iw, l, nd(3), ngto, niso, & 1773 nkind, np, nr, output_unit, & 1774 print_level, unit_nr 1775 LOGICAL :: append_cube, do_kpoints, mpi_io, omit_headers, print_it, print_total_density, & 1776 rho_r_valid, write_ks, write_xc, xrd_interface 1777 REAL(KIND=dp) :: q_max, rho_hard, rho_soft, rho_total, & 1778 rho_total_rspace, udvol, volume 1779 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: bfun 1780 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: aedens, ccdens, ppdens 1781 REAL(KIND=dp), DIMENSION(3) :: dr 1782 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1783 TYPE(atomic_kind_type), POINTER :: atomic_kind 1784 TYPE(cell_type), POINTER :: cell 1785 TYPE(cp_logger_type), POINTER :: logger 1786 TYPE(cp_para_env_type), POINTER :: para_env 1787 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hr 1788 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_rmpv, matrix_vxc, rho_ao 1789 TYPE(dft_control_type), POINTER :: dft_control 1790 TYPE(grid_atom_type), POINTER :: grid_atom 1791 TYPE(particle_list_type), POINTER :: particles 1792 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1793 TYPE(pw_env_type), POINTER :: pw_env 1794 TYPE(pw_p_type) :: aux_g, aux_r, rho_elec_gspace, & 1795 rho_elec_rspace, wf_r 1796 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r 1797 TYPE(pw_p_type), POINTER :: rho0_s_gs, rho_core, vee 1798 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 1799 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 1800 TYPE(pw_type), POINTER :: v_hartree_rspace 1801 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1802 TYPE(qs_kind_type), POINTER :: qs_kind 1803 TYPE(qs_rho_type), POINTER :: rho 1804 TYPE(qs_subsys_type), POINTER :: subsys 1805 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set 1806 TYPE(rho_atom_type), POINTER :: rho_atom 1807 TYPE(section_vals_type), POINTER :: dft_section, input, print_key, xc_section 1808 1809 CALL timeset(routineN, handle) 1810 NULLIFY (cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, & 1811 atomic_kind_set, qs_kind_set, particle_set, rho, ks_rmpv, rho_ao, rho_r, & 1812 dft_section, xc_section, input, particles, subsys, matrix_vxc, v_hartree_rspace, vee) 1813 1814 logger => cp_get_default_logger() 1815 output_unit = cp_logger_get_default_io_unit(logger) 1816 1817 CPASSERT(ASSOCIATED(qs_env)) 1818 CALL get_qs_env(qs_env, & 1819 atomic_kind_set=atomic_kind_set, & 1820 qs_kind_set=qs_kind_set, & 1821 particle_set=particle_set, & 1822 cell=cell, & 1823 para_env=para_env, & 1824 dft_control=dft_control, & 1825 input=input, & 1826 do_kpoints=do_kpoints, & 1827 subsys=subsys) 1828 dft_section => section_vals_get_subs_vals(input, "DFT") 1829 CALL qs_subsys_get(subsys, particles=particles) 1830 1831 ! Print the total density (electronic + core charge) 1832 CALL get_qs_env(qs_env, rho=rho) 1833 CALL qs_rho_get(rho, rho_r=rho_r) 1834 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 1835 "DFT%PRINT%TOT_DENSITY_CUBE"), cp_p_file)) THEN 1836 NULLIFY (rho_core, rho0_s_gs) 1837 append_cube = section_get_lval(input, "DFT%PRINT%TOT_DENSITY_CUBE%APPEND") 1838 my_pos_cube = "REWIND" 1839 IF (append_cube) THEN 1840 my_pos_cube = "APPEND" 1841 END IF 1842 1843 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho_core=rho_core, & 1844 rho0_s_gs=rho0_s_gs) 1845 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 1846 pw_pools=pw_pools) 1847 CALL pw_pool_create_pw(auxbas_pw_pool, wf_r%pw, & 1848 use_data=REALDATA3D, & 1849 in_space=REALSPACE) 1850 IF (dft_control%qs_control%gapw) THEN 1851 CALL pw_transfer(rho0_s_gs%pw, wf_r%pw) 1852 IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN 1853 CALL pw_axpy(rho_core%pw, wf_r%pw) 1854 END IF 1855 ELSE 1856 CALL pw_transfer(rho_core%pw, wf_r%pw) 1857 END IF 1858 DO ispin = 1, dft_control%nspins 1859 CALL pw_axpy(rho_r(ispin)%pw, wf_r%pw) 1860 END DO 1861 filename = "TOTAL_DENSITY" 1862 mpi_io = .TRUE. 1863 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%TOT_DENSITY_CUBE", & 1864 extension=".cube", middle_name=TRIM(filename), file_position=my_pos_cube, & 1865 log_filename=.FALSE., mpi_io=mpi_io) 1866 CALL cp_pw_to_cube(wf_r%pw, unit_nr, "TOTAL DENSITY", & 1867 particles=particles, & 1868 stride=section_get_ivals(dft_section, "PRINT%TOT_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) 1869 CALL cp_print_key_finished_output(unit_nr, logger, input, & 1870 "DFT%PRINT%TOT_DENSITY_CUBE", mpi_io=mpi_io) 1871 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_r%pw) 1872 END IF 1873 1874 ! Write cube file with electron density 1875 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 1876 "DFT%PRINT%E_DENSITY_CUBE"), cp_p_file)) THEN 1877 CALL section_vals_val_get(dft_section, & 1878 keyword_name="PRINT%E_DENSITY_CUBE%TOTAL_DENSITY", & 1879 l_val=print_total_density) 1880 append_cube = section_get_lval(input, "DFT%PRINT%E_DENSITY_CUBE%APPEND") 1881 my_pos_cube = "REWIND" 1882 IF (append_cube) THEN 1883 my_pos_cube = "APPEND" 1884 END IF 1885 ! Write the info on core densities for the interface between cp2k and the XRD code 1886 ! together with the valence density they are used to compute the form factor (Fourier transform) 1887 xrd_interface = section_get_lval(input, "DFT%PRINT%E_DENSITY_CUBE%XRD_INTERFACE") 1888 IF (xrd_interface) THEN 1889 !cube file only contains soft density (GAPW) 1890 IF (dft_control%qs_control%gapw) print_total_density = .FALSE. 1891 1892 filename = "ELECTRON_DENSITY" 1893 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & 1894 extension=".xrd", middle_name=TRIM(filename), & 1895 file_position=my_pos_cube, log_filename=.FALSE.) 1896 ngto = section_get_ival(input, "DFT%PRINT%E_DENSITY_CUBE%NGAUSS") 1897 IF (output_unit > 0) THEN 1898 INQUIRE (UNIT=unit_nr, NAME=filename) 1899 WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & 1900 "The electron density (atomic part) is written to the file:", & 1901 TRIM(filename) 1902 END IF 1903 1904 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 1905 nkind = SIZE(atomic_kind_set) 1906 IF (unit_nr > 0) THEN 1907 WRITE (unit_nr, *) "Atomic (core) densities" 1908 WRITE (unit_nr, *) "Unit cell" 1909 WRITE (unit_nr, FMT="(3F20.12)") cell%hmat(1, 1), cell%hmat(1, 2), cell%hmat(1, 3) 1910 WRITE (unit_nr, FMT="(3F20.12)") cell%hmat(2, 1), cell%hmat(2, 2), cell%hmat(2, 3) 1911 WRITE (unit_nr, FMT="(3F20.12)") cell%hmat(3, 1), cell%hmat(3, 2), cell%hmat(3, 3) 1912 WRITE (unit_nr, *) "Atomic types" 1913 WRITE (unit_nr, *) nkind 1914 END IF 1915 ! calculate atomic density and core density 1916 ALLOCATE (ppdens(ngto, 2, nkind), aedens(ngto, 2, nkind), ccdens(ngto, 2, nkind)) 1917 DO ikind = 1, nkind 1918 atomic_kind => atomic_kind_set(ikind) 1919 qs_kind => qs_kind_set(ikind) 1920 CALL get_atomic_kind(atomic_kind, name=name, element_symbol=element_symbol) 1921 CALL calculate_atomic_density(ppdens(:, :, ikind), atomic_kind, qs_kind, ngto, & 1922 iunit=output_unit, confine=.TRUE.) 1923 CALL calculate_atomic_density(aedens(:, :, ikind), atomic_kind, qs_kind, ngto, & 1924 iunit=output_unit, allelectron=.TRUE., confine=.TRUE.) 1925 ccdens(:, 1, ikind) = aedens(:, 1, ikind) 1926 ccdens(:, 2, ikind) = 0._dp 1927 CALL project_function_a(ccdens(1:ngto, 2, ikind), ccdens(1:ngto, 1, ikind), & 1928 ppdens(1:ngto, 2, ikind), ppdens(1:ngto, 1, ikind), 0) 1929 ccdens(:, 2, ikind) = aedens(:, 2, ikind) - ccdens(:, 2, ikind) 1930 IF (unit_nr > 0) THEN 1931 WRITE (unit_nr, FMT="(I6,A10,A20)") ikind, TRIM(element_symbol), TRIM(name) 1932 WRITE (unit_nr, FMT="(I6)") ngto 1933 WRITE (unit_nr, *) " Total density" 1934 WRITE (unit_nr, FMT="(2G24.12)") (aedens(i, 1, ikind), aedens(i, 2, ikind), i=1, ngto) 1935 WRITE (unit_nr, *) " Core density" 1936 WRITE (unit_nr, FMT="(2G24.12)") (ccdens(i, 1, ikind), ccdens(i, 2, ikind), i=1, ngto) 1937 END IF 1938 NULLIFY (atomic_kind) 1939 END DO 1940 1941 IF (dft_control%qs_control%gapw) THEN 1942 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set) 1943 1944 IF (unit_nr > 0) THEN 1945 WRITE (unit_nr, *) "Coordinates and GAPW density" 1946 END IF 1947 np = particles%n_els 1948 DO iat = 1, np 1949 CALL get_atomic_kind(particles%els(iat)%atomic_kind, kind_number=ikind) 1950 CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid_atom) 1951 rho_atom => rho_atom_set(iat) 1952 IF (ASSOCIATED(rho_atom%rho_rad_h(1)%r_coef)) THEN 1953 nr = SIZE(rho_atom%rho_rad_h(1)%r_coef, 1) 1954 niso = SIZE(rho_atom%rho_rad_h(1)%r_coef, 2) 1955 ELSE 1956 nr = 0 1957 niso = 0 1958 ENDIF 1959 CALL mp_sum(nr, para_env%group) 1960 CALL mp_sum(niso, para_env%group) 1961 1962 ALLOCATE (bfun(nr, niso)) 1963 bfun = 0._dp 1964 DO ispin = 1, dft_control%nspins 1965 IF (ASSOCIATED(rho_atom%rho_rad_h(1)%r_coef)) THEN 1966 bfun(:, :) = bfun + rho_atom%rho_rad_h(ispin)%r_coef - rho_atom%rho_rad_s(ispin)%r_coef 1967 END IF 1968 END DO 1969 CALL mp_sum(bfun, para_env%group) 1970 ccdens(:, 1, ikind) = ppdens(:, 1, ikind) 1971 ccdens(:, 2, ikind) = 0._dp 1972 IF (unit_nr > 0) THEN 1973 WRITE (unit_nr, '(I10,I5,3f12.6)') iat, ikind, particles%els(iat)%r 1974 END IF 1975 DO iso = 1, niso 1976 l = indso(1, iso) 1977 CALL project_function_b(ccdens(:, 2, ikind), ccdens(:, 1, ikind), bfun(:, iso), grid_atom, l) 1978 IF (unit_nr > 0) THEN 1979 WRITE (unit_nr, FMT="(3I6)") iso, l, ngto 1980 WRITE (unit_nr, FMT="(2G24.12)") (ccdens(i, 1, ikind), ccdens(i, 2, ikind), i=1, ngto) 1981 END IF 1982 END DO 1983 DEALLOCATE (bfun) 1984 END DO 1985 ELSE 1986 IF (unit_nr > 0) THEN 1987 WRITE (unit_nr, *) "Coordinates" 1988 np = particles%n_els 1989 DO iat = 1, np 1990 CALL get_atomic_kind(particles%els(iat)%atomic_kind, kind_number=ikind) 1991 WRITE (unit_nr, '(I10,I5,3f12.6)') iat, ikind, particles%els(iat)%r 1992 END DO 1993 END IF 1994 END IF 1995 1996 DEALLOCATE (ppdens, aedens, ccdens) 1997 1998 CALL cp_print_key_finished_output(unit_nr, logger, input, & 1999 "DFT%PRINT%E_DENSITY_CUBE") 2000 2001 END IF 2002 IF (dft_control%qs_control%gapw .AND. print_total_density) THEN 2003 ! total density in g-space not implemented for k-points 2004 CPASSERT(.NOT. do_kpoints) 2005 ! Print total electronic density 2006 CALL get_qs_env(qs_env=qs_env, & 2007 pw_env=pw_env) 2008 CALL pw_env_get(pw_env=pw_env, & 2009 auxbas_pw_pool=auxbas_pw_pool, & 2010 pw_pools=pw_pools) 2011 CALL pw_pool_create_pw(pool=auxbas_pw_pool, & 2012 pw=rho_elec_rspace%pw, & 2013 use_data=REALDATA3D, & 2014 in_space=REALSPACE) 2015 CALL pw_zero(rho_elec_rspace%pw) 2016 CALL pw_pool_create_pw(pool=auxbas_pw_pool, & 2017 pw=rho_elec_gspace%pw, & 2018 use_data=COMPLEXDATA1D, & 2019 in_space=RECIPROCALSPACE) 2020 CALL pw_zero(rho_elec_gspace%pw) 2021 CALL get_pw_grid_info(pw_grid=rho_elec_gspace%pw%pw_grid, & 2022 dr=dr, & 2023 vol=volume) 2024 q_max = SQRT(SUM((pi/dr(:))**2)) 2025 CALL calculate_rhotot_elec_gspace(qs_env=qs_env, & 2026 auxbas_pw_pool=auxbas_pw_pool, & 2027 rhotot_elec_gspace=rho_elec_gspace, & 2028 q_max=q_max, & 2029 rho_hard=rho_hard, & 2030 rho_soft=rho_soft) 2031 rho_total = rho_hard + rho_soft 2032 CALL get_pw_grid_info(pw_grid=rho_elec_gspace%pw%pw_grid, & 2033 vol=volume) 2034 CALL pw_transfer(rho_elec_gspace%pw, rho_elec_rspace%pw, debug=.FALSE.) 2035 rho_total_rspace = pw_integrate_function(rho_elec_rspace%pw, isign=-1)/volume 2036 filename = "TOTAL_ELECTRON_DENSITY" 2037 mpi_io = .TRUE. 2038 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & 2039 extension=".cube", middle_name=TRIM(filename), & 2040 file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & 2041 fout=mpi_filename) 2042 IF (output_unit > 0) THEN 2043 IF (.NOT. mpi_io) THEN 2044 INQUIRE (UNIT=unit_nr, NAME=filename) 2045 ELSE 2046 filename = mpi_filename 2047 END IF 2048 WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & 2049 "The total electron density is written in cube file format to the file:", & 2050 TRIM(filename) 2051 WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") & 2052 "q(max) [1/Angstrom] :", q_max/angstrom, & 2053 "Soft electronic charge (G-space) :", rho_soft, & 2054 "Hard electronic charge (G-space) :", rho_hard, & 2055 "Total electronic charge (G-space):", rho_total, & 2056 "Total electronic charge (R-space):", rho_total_rspace 2057 END IF 2058 CALL cp_pw_to_cube(rho_elec_rspace%pw, unit_nr, "TOTAL ELECTRON DENSITY", & 2059 particles=particles, & 2060 stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) 2061 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2062 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) 2063 ! Print total spin density for spin-polarized systems 2064 IF (dft_control%nspins > 1) THEN 2065 CALL pw_zero(rho_elec_gspace%pw) 2066 CALL pw_zero(rho_elec_rspace%pw) 2067 CALL calculate_rhotot_elec_gspace(qs_env=qs_env, & 2068 auxbas_pw_pool=auxbas_pw_pool, & 2069 rhotot_elec_gspace=rho_elec_gspace, & 2070 q_max=q_max, & 2071 rho_hard=rho_hard, & 2072 rho_soft=rho_soft, & 2073 fsign=-1.0_dp) 2074 rho_total = rho_hard + rho_soft 2075 CALL pw_transfer(rho_elec_gspace%pw, rho_elec_rspace%pw, debug=.FALSE.) 2076 rho_total_rspace = pw_integrate_function(rho_elec_rspace%pw, isign=-1)/volume 2077 filename = "TOTAL_SPIN_DENSITY" 2078 mpi_io = .TRUE. 2079 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & 2080 extension=".cube", middle_name=TRIM(filename), & 2081 file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & 2082 fout=mpi_filename) 2083 IF (output_unit > 0) THEN 2084 IF (.NOT. mpi_io) THEN 2085 INQUIRE (UNIT=unit_nr, NAME=filename) 2086 ELSE 2087 filename = mpi_filename 2088 END IF 2089 WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & 2090 "The total spin density is written in cube file format to the file:", & 2091 TRIM(filename) 2092 WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") & 2093 "q(max) [1/Angstrom] :", q_max/angstrom, & 2094 "Soft part of the spin density (G-space):", rho_soft, & 2095 "Hard part of the spin density (G-space):", rho_hard, & 2096 "Total spin density (G-space) :", rho_total, & 2097 "Total spin density (R-space) :", rho_total_rspace 2098 END IF 2099 CALL cp_pw_to_cube(rho_elec_rspace%pw, unit_nr, "TOTAL SPIN DENSITY", & 2100 particles=particles, & 2101 stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) 2102 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2103 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) 2104 END IF 2105 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_elec_gspace%pw) 2106 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_elec_rspace%pw) 2107 ELSE 2108 IF (dft_control%nspins > 1) THEN 2109 CALL get_qs_env(qs_env=qs_env, & 2110 pw_env=pw_env) 2111 CALL pw_env_get(pw_env=pw_env, & 2112 auxbas_pw_pool=auxbas_pw_pool, & 2113 pw_pools=pw_pools) 2114 CALL pw_pool_create_pw(pool=auxbas_pw_pool, & 2115 pw=rho_elec_rspace%pw, & 2116 use_data=REALDATA3D, & 2117 in_space=REALSPACE) 2118 CALL pw_copy(rho_r(1)%pw, rho_elec_rspace%pw) 2119 CALL pw_axpy(rho_r(2)%pw, rho_elec_rspace%pw) 2120 filename = "ELECTRON_DENSITY" 2121 mpi_io = .TRUE. 2122 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & 2123 extension=".cube", middle_name=TRIM(filename), & 2124 file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & 2125 fout=mpi_filename) 2126 IF (output_unit > 0) THEN 2127 IF (.NOT. mpi_io) THEN 2128 INQUIRE (UNIT=unit_nr, NAME=filename) 2129 ELSE 2130 filename = mpi_filename 2131 END IF 2132 WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & 2133 "The sum of alpha and beta density is written in cube file format to the file:", & 2134 TRIM(filename) 2135 END IF 2136 CALL cp_pw_to_cube(rho_elec_rspace%pw, unit_nr, "SUM OF ALPHA AND BETA DENSITY", & 2137 particles=particles, stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), & 2138 mpi_io=mpi_io) 2139 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2140 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) 2141 CALL pw_copy(rho_r(1)%pw, rho_elec_rspace%pw) 2142 CALL pw_axpy(rho_r(2)%pw, rho_elec_rspace%pw, alpha=-1.0_dp) 2143 filename = "SPIN_DENSITY" 2144 mpi_io = .TRUE. 2145 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & 2146 extension=".cube", middle_name=TRIM(filename), & 2147 file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & 2148 fout=mpi_filename) 2149 IF (output_unit > 0) THEN 2150 IF (.NOT. mpi_io) THEN 2151 INQUIRE (UNIT=unit_nr, NAME=filename) 2152 ELSE 2153 filename = mpi_filename 2154 END IF 2155 WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & 2156 "The spin density is written in cube file format to the file:", & 2157 TRIM(filename) 2158 END IF 2159 CALL cp_pw_to_cube(rho_elec_rspace%pw, unit_nr, "SPIN DENSITY", & 2160 particles=particles, & 2161 stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) 2162 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2163 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) 2164 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_elec_rspace%pw) 2165 ELSE 2166 filename = "ELECTRON_DENSITY" 2167 mpi_io = .TRUE. 2168 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%E_DENSITY_CUBE", & 2169 extension=".cube", middle_name=TRIM(filename), & 2170 file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, & 2171 fout=mpi_filename) 2172 IF (output_unit > 0) THEN 2173 IF (.NOT. mpi_io) THEN 2174 INQUIRE (UNIT=unit_nr, NAME=filename) 2175 ELSE 2176 filename = mpi_filename 2177 END IF 2178 WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & 2179 "The electron density is written in cube file format to the file:", & 2180 TRIM(filename) 2181 END IF 2182 CALL cp_pw_to_cube(rho_r(1)%pw, unit_nr, "ELECTRON DENSITY", & 2183 particles=particles, & 2184 stride=section_get_ivals(dft_section, "PRINT%E_DENSITY_CUBE%STRIDE"), mpi_io=mpi_io) 2185 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2186 "DFT%PRINT%E_DENSITY_CUBE", mpi_io=mpi_io) 2187 END IF ! nspins 2188 END IF ! total density for GAPW 2189 END IF ! print key 2190 2191 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 2192 dft_section, "PRINT%ENERGY_WINDOWS"), cp_p_file) .AND. .NOT. do_kpoints) THEN 2193 CALL energy_windows(qs_env) 2194 ENDIF 2195 2196 ! Print the hartree potential 2197 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 2198 "DFT%PRINT%V_HARTREE_CUBE"), cp_p_file)) THEN 2199 2200 CALL get_qs_env(qs_env=qs_env, & 2201 pw_env=pw_env, & 2202 v_hartree_rspace=v_hartree_rspace) 2203 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 2204 CALL pw_pool_create_pw(auxbas_pw_pool, aux_r%pw, & 2205 use_data=REALDATA3D, & 2206 in_space=REALSPACE) 2207 2208 append_cube = section_get_lval(input, "DFT%PRINT%V_HARTREE_CUBE%APPEND") 2209 my_pos_cube = "REWIND" 2210 IF (append_cube) THEN 2211 my_pos_cube = "APPEND" 2212 END IF 2213 mpi_io = .TRUE. 2214 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 2215 CALL pw_env_get(pw_env) 2216 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%V_HARTREE_CUBE", & 2217 extension=".cube", middle_name="v_hartree", file_position=my_pos_cube, mpi_io=mpi_io) 2218 udvol = 1.0_dp/v_hartree_rspace%pw_grid%dvol 2219 2220 CALL pw_copy(v_hartree_rspace, aux_r%pw) 2221 CALL pw_scale(aux_r%pw, udvol) 2222 2223 CALL cp_pw_to_cube(aux_r%pw, unit_nr, "HARTREE POTENTIAL", particles=particles, & 2224 stride=section_get_ivals(dft_section, & 2225 "PRINT%V_HARTREE_CUBE%STRIDE"), mpi_io=mpi_io) 2226 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2227 "DFT%PRINT%V_HARTREE_CUBE", mpi_io=mpi_io) 2228 2229 CALL pw_pool_give_back_pw(auxbas_pw_pool, aux_r%pw) 2230 ENDIF 2231 2232 ! Print the external potential 2233 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 2234 "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE"), cp_p_file)) THEN 2235 IF (dft_control%apply_external_potential) THEN 2236 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, vee=vee) 2237 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 2238 CALL pw_pool_create_pw(auxbas_pw_pool, aux_r%pw, & 2239 use_data=REALDATA3D, & 2240 in_space=REALSPACE) 2241 2242 append_cube = section_get_lval(input, "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE%APPEND") 2243 my_pos_cube = "REWIND" 2244 IF (append_cube) THEN 2245 my_pos_cube = "APPEND" 2246 END IF 2247 mpi_io = .TRUE. 2248 CALL pw_env_get(pw_env) 2249 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE", & 2250 extension=".cube", middle_name="ext_pot", file_position=my_pos_cube, mpi_io=mpi_io) 2251 2252 CALL pw_copy(vee%pw, aux_r%pw) 2253 2254 CALL cp_pw_to_cube(aux_r%pw, unit_nr, "EXTERNAL POTENTIAL", particles=particles, & 2255 stride=section_get_ivals(dft_section, & 2256 "PRINT%EXTERNAL_POTENTIAL_CUBE%STRIDE"), mpi_io=mpi_io) 2257 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2258 "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE", mpi_io=mpi_io) 2259 2260 CALL pw_pool_give_back_pw(auxbas_pw_pool, aux_r%pw) 2261 ENDIF 2262 ENDIF 2263 2264 ! Print the Electrical Field Components 2265 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 2266 "DFT%PRINT%EFIELD_CUBE"), cp_p_file)) THEN 2267 2268 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 2269 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 2270 CALL pw_pool_create_pw(auxbas_pw_pool, aux_r%pw, & 2271 use_data=REALDATA3D, & 2272 in_space=REALSPACE) 2273 CALL pw_pool_create_pw(auxbas_pw_pool, aux_g%pw, & 2274 use_data=COMPLEXDATA1D, & 2275 in_space=RECIPROCALSPACE) 2276 2277 append_cube = section_get_lval(input, "DFT%PRINT%EFIELD_CUBE%APPEND") 2278 my_pos_cube = "REWIND" 2279 IF (append_cube) THEN 2280 my_pos_cube = "APPEND" 2281 END IF 2282 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, & 2283 v_hartree_rspace=v_hartree_rspace) 2284 CALL pw_env_get(pw_env) 2285 udvol = 1.0_dp/v_hartree_rspace%pw_grid%dvol 2286 DO id = 1, 3 2287 mpi_io = .TRUE. 2288 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EFIELD_CUBE", & 2289 extension=".cube", middle_name="efield_"//cdir(id), file_position=my_pos_cube, & 2290 mpi_io=mpi_io) 2291 2292 CALL pw_transfer(v_hartree_rspace, aux_g%pw) 2293 nd = 0 2294 nd(id) = 1 2295 CALL pw_derive(aux_g%pw, nd) 2296 CALL pw_transfer(aux_g%pw, aux_r%pw) 2297 CALL pw_scale(aux_r%pw, udvol) 2298 2299 CALL cp_pw_to_cube(aux_r%pw, & 2300 unit_nr, "ELECTRIC FIELD", particles=particles, & 2301 stride=section_get_ivals(dft_section, & 2302 "PRINT%EFIELD_CUBE%STRIDE"), mpi_io=mpi_io) 2303 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2304 "DFT%PRINT%EFIELD_CUBE", mpi_io=mpi_io) 2305 END DO 2306 2307 CALL pw_pool_give_back_pw(auxbas_pw_pool, aux_r%pw) 2308 CALL pw_pool_give_back_pw(auxbas_pw_pool, aux_g%pw) 2309 END IF 2310 2311 ! Write cube files from the local energy 2312 CALL qs_scf_post_local_energy(input, logger, qs_env) 2313 2314 ! Write cube files from the local stress tensor 2315 CALL qs_scf_post_local_stress(input, logger, qs_env) 2316 2317 ! Write cube files from the implicit Poisson solver 2318 CALL qs_scf_post_ps_implicit(input, logger, qs_env) 2319 2320 ! post SCF Transport 2321 CALL qs_scf_post_transport(qs_env) 2322 2323 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers) 2324 ! Write the density matrices 2325 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 2326 "DFT%PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN 2327 iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/DENSITY", & 2328 extension=".Log") 2329 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 2330 CALL qs_rho_get(rho, rho_ao_kp=rho_ao) 2331 after = MIN(MAX(after, 1), 16) 2332 DO ispin = 1, dft_control%nspins 2333 DO img = 1, dft_control%nimages 2334 CALL cp_dbcsr_write_sparse_matrix(rho_ao(ispin, img)%matrix, 4, after, qs_env, & 2335 para_env, output_unit=iw, omit_headers=omit_headers) 2336 END DO 2337 END DO 2338 CALL cp_print_key_finished_output(iw, logger, input, & 2339 "DFT%PRINT%AO_MATRICES/DENSITY") 2340 END IF 2341 2342 ! Write the Kohn-Sham matrices 2343 write_ks = BTEST(cp_print_key_should_output(logger%iter_info, input, & 2344 "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file) 2345 write_xc = BTEST(cp_print_key_should_output(logger%iter_info, input, & 2346 "DFT%PRINT%AO_MATRICES/MATRIX_VXC"), cp_p_file) 2347 ! we need to update stuff before writing, potentially computing the matrix_vxc 2348 IF (write_ks .OR. write_xc) THEN 2349 IF (write_xc) qs_env%requires_matrix_vxc = .TRUE. 2350 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 2351 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., & 2352 just_energy=.FALSE.) 2353 IF (write_xc) qs_env%requires_matrix_vxc = .FALSE. 2354 END IF 2355 2356 ! Write the Kohn-Sham matrices 2357 IF (write_ks) THEN 2358 iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", & 2359 extension=".Log") 2360 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=ks_rmpv) 2361 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 2362 after = MIN(MAX(after, 1), 16) 2363 DO ispin = 1, dft_control%nspins 2364 DO img = 1, dft_control%nimages 2365 CALL cp_dbcsr_write_sparse_matrix(ks_rmpv(ispin, img)%matrix, 4, after, qs_env, & 2366 para_env, output_unit=iw, omit_headers=omit_headers) 2367 END DO 2368 END DO 2369 CALL cp_print_key_finished_output(iw, logger, input, & 2370 "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX") 2371 END IF 2372 2373 ! write csr matrices 2374 ! matrices in terms of the PAO basis will be taken care of in pao_post_scf. 2375 IF (.NOT. dft_control%qs_control%pao) THEN 2376 CALL write_ks_matrix_csr(qs_env, input) 2377 CALL write_s_matrix_csr(qs_env, input) 2378 END IF 2379 2380 ! write adjacency matrix 2381 CALL write_adjacency_matrix(qs_env, input) 2382 2383 ! Write the xc matrix 2384 IF (write_xc) THEN 2385 CALL get_qs_env(qs_env=qs_env, matrix_vxc_kp=matrix_vxc) 2386 CPASSERT(ASSOCIATED(matrix_vxc)) 2387 iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/MATRIX_VXC", & 2388 extension=".Log") 2389 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 2390 after = MIN(MAX(after, 1), 16) 2391 DO ispin = 1, dft_control%nspins 2392 DO img = 1, dft_control%nimages 2393 CALL cp_dbcsr_write_sparse_matrix(matrix_vxc(ispin, img)%matrix, 4, after, qs_env, & 2394 para_env, output_unit=iw, omit_headers=omit_headers) 2395 END DO 2396 END DO 2397 CALL cp_print_key_finished_output(iw, logger, input, & 2398 "DFT%PRINT%AO_MATRICES/MATRIX_VXC") 2399 END IF 2400 2401 ! Write the [H,r] commutator matrices 2402 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 2403 "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR"), cp_p_file)) THEN 2404 iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR", & 2405 extension=".Log") 2406 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 2407 NULLIFY (matrix_hr) 2408 CALL build_com_hr_matrix(qs_env, matrix_hr) 2409 after = MIN(MAX(after, 1), 16) 2410 DO img = 1, 3 2411 CALL cp_dbcsr_write_sparse_matrix(matrix_hr(img)%matrix, 4, after, qs_env, & 2412 para_env, output_unit=iw, omit_headers=omit_headers) 2413 END DO 2414 CALL dbcsr_deallocate_matrix_set(matrix_hr) 2415 CALL cp_print_key_finished_output(iw, logger, input, & 2416 "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR") 2417 END IF 2418 2419 ! Compute the Mulliken charges 2420 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MULLIKEN") 2421 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 2422 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MULLIKEN", extension=".mulliken", log_filename=.FALSE.) 2423 print_level = 1 2424 CALL section_vals_val_get(print_key, "PRINT_GOP", l_val=print_it) 2425 IF (print_it) print_level = 2 2426 CALL section_vals_val_get(print_key, "PRINT_ALL", l_val=print_it) 2427 IF (print_it) print_level = 3 2428 CALL mulliken_population_analysis(qs_env, unit_nr, print_level) 2429 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MULLIKEN") 2430 END IF 2431 2432 ! Compute the Hirshfeld charges 2433 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%HIRSHFELD") 2434 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 2435 ! we check if real space density is available 2436 NULLIFY (rho) 2437 CALL get_qs_env(qs_env=qs_env, rho=rho) 2438 CALL qs_rho_get(rho, rho_r_valid=rho_r_valid) 2439 IF (rho_r_valid) THEN 2440 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%HIRSHFELD", extension=".hirshfeld", log_filename=.FALSE.) 2441 CALL hirshfeld_charges(qs_env, print_key, unit_nr) 2442 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%HIRSHFELD") 2443 END IF 2444 END IF 2445 2446 ! MAO analysis 2447 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MAO_ANALYSIS") 2448 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 2449 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MAO_ANALYSIS", extension=".mao", log_filename=.FALSE.) 2450 CALL mao_analysis(qs_env, print_key, unit_nr) 2451 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MAO_ANALYSIS") 2452 END IF 2453 2454 ! MINBAS analysis 2455 print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MINBAS_ANALYSIS") 2456 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN 2457 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MINBAS_ANALYSIS", extension=".mao", log_filename=.FALSE.) 2458 CALL minbas_analysis(qs_env, print_key, unit_nr) 2459 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MINBAS_ANALYSIS") 2460 END IF 2461 2462 CALL timestop(handle) 2463 2464 END SUBROUTINE write_mo_free_results 2465 2466! ************************************************************************************************** 2467!> \brief Calculates Hirshfeld charges 2468!> \param qs_env the qs_env where to calculate the charges 2469!> \param input_section the input section for Hirshfeld charges 2470!> \param unit_nr the output unit number 2471! ************************************************************************************************** 2472 SUBROUTINE hirshfeld_charges(qs_env, input_section, unit_nr) 2473 TYPE(qs_environment_type), POINTER :: qs_env 2474 TYPE(section_vals_type), POINTER :: input_section 2475 INTEGER, INTENT(IN) :: unit_nr 2476 2477 CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_charges', & 2478 routineP = moduleN//':'//routineN 2479 2480 INTEGER :: i, iat, ikind, natom, nkind, nspin, & 2481 radius_type, refc, shapef 2482 INTEGER, DIMENSION(:), POINTER :: atom_list 2483 LOGICAL :: do_radius, do_sc, paw_atom 2484 REAL(KIND=dp) :: zeff 2485 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 2486 REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges 2487 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2488 TYPE(atomic_kind_type), POINTER :: atomic_kind 2489 TYPE(cp_para_env_type), POINTER :: para_env 2490 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s 2491 TYPE(dft_control_type), POINTER :: dft_control 2492 TYPE(hirshfeld_type), POINTER :: hirshfeld_env 2493 TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho 2494 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 2495 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 2496 TYPE(qs_rho_type), POINTER :: rho 2497 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 2498 2499 NULLIFY (hirshfeld_env) 2500 NULLIFY (radii) 2501 CALL create_hirshfeld_type(hirshfeld_env) 2502 ! 2503 CALL get_qs_env(qs_env, nkind=nkind, natom=natom) 2504 ALLOCATE (hirshfeld_env%charges(natom)) 2505 ! input options 2506 CALL section_vals_val_get(input_section, "SELF_CONSISTENT", l_val=do_sc) 2507 CALL section_vals_val_get(input_section, "USER_RADIUS", l_val=do_radius) 2508 CALL section_vals_val_get(input_section, "SHAPE_FUNCTION", i_val=shapef) 2509 CALL section_vals_val_get(input_section, "REFERENCE_CHARGE", i_val=refc) 2510 IF (do_radius) THEN 2511 radius_type = radius_user 2512 CALL section_vals_val_get(input_section, "ATOMIC_RADII", r_vals=radii) 2513 IF (.NOT. SIZE(radii) == nkind) & 2514 CALL cp_abort(__LOCATION__, & 2515 "Length of keyword HIRSHFELD\ATOMIC_RADII does not "// & 2516 "match number of atomic kinds in the input coordinate file.") 2517 ELSE 2518 radius_type = radius_covalent 2519 END IF 2520 CALL set_hirshfeld_info(hirshfeld_env, shape_function_type=shapef, & 2521 iterative=do_sc, ref_charge=refc, & 2522 radius_type=radius_type) 2523 ! shape function 2524 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set) 2525 CALL create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, & 2526 radii_list=radii) 2527 ! reference charges 2528 CALL get_qs_env(qs_env, rho=rho) 2529 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 2530 nspin = SIZE(matrix_p, 1) 2531 ALLOCATE (charges(natom, nspin)) 2532 SELECT CASE (refc) 2533 CASE (ref_charge_atomic) 2534 DO ikind = 1, nkind 2535 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff) 2536 atomic_kind => atomic_kind_set(ikind) 2537 CALL get_atomic_kind(atomic_kind, atom_list=atom_list) 2538 DO iat = 1, SIZE(atom_list) 2539 i = atom_list(iat) 2540 hirshfeld_env%charges(i) = zeff 2541 END DO 2542 END DO 2543 CASE (ref_charge_mulliken) 2544 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env) 2545 CALL mulliken_charges(matrix_p, matrix_s, para_env, charges) 2546 DO iat = 1, natom 2547 hirshfeld_env%charges(iat) = SUM(charges(iat, :)) 2548 END DO 2549 CASE DEFAULT 2550 CPABORT("Unknown type of reference charge for Hirshfeld partitioning.") 2551 END SELECT 2552 ! 2553 charges = 0.0_dp 2554 IF (hirshfeld_env%iterative) THEN 2555 ! Hirshfeld-I charges 2556 CALL comp_hirshfeld_i_charges(qs_env, hirshfeld_env, charges, unit_nr) 2557 ELSE 2558 ! Hirshfeld charges 2559 CALL comp_hirshfeld_charges(qs_env, hirshfeld_env, charges) 2560 END IF 2561 CALL get_qs_env(qs_env, particle_set=particle_set, dft_control=dft_control) 2562 IF (dft_control%qs_control%gapw) THEN 2563 ! GAPW: add core charges (rho_hard - rho_soft) 2564 CALL get_qs_env(qs_env, rho0_mpole=rho0_mpole) 2565 CALL get_rho0_mpole(rho0_mpole, mp_rho=mp_rho) 2566 DO iat = 1, natom 2567 atomic_kind => particle_set(iat)%atomic_kind 2568 CALL get_atomic_kind(atomic_kind, kind_number=ikind) 2569 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom) 2570 IF (paw_atom) THEN 2571 charges(iat, 1:nspin) = charges(iat, 1:nspin) + mp_rho(iat)%q0(1:nspin) 2572 END IF 2573 END DO 2574 END IF 2575 ! 2576 IF (unit_nr > 0) THEN 2577 CALL write_hirshfeld_charges(charges, hirshfeld_env, particle_set, & 2578 qs_kind_set, unit_nr) 2579 END IF 2580 ! 2581 CALL release_hirshfeld_type(hirshfeld_env) 2582 DEALLOCATE (charges) 2583 2584 END SUBROUTINE hirshfeld_charges 2585 2586! ************************************************************************************************** 2587!> \brief ... 2588!> \param ca ... 2589!> \param a ... 2590!> \param cb ... 2591!> \param b ... 2592!> \param l ... 2593! ************************************************************************************************** 2594 SUBROUTINE project_function_a(ca, a, cb, b, l) 2595 ! project function cb on ca 2596 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ca 2597 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, cb, b 2598 INTEGER, INTENT(IN) :: l 2599 2600 CHARACTER(len=*), PARAMETER :: routineN = 'project_function_a', & 2601 routineP = moduleN//':'//routineN 2602 2603 INTEGER :: info, n 2604 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv 2605 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smat, tmat, v 2606 2607 n = SIZE(ca) 2608 ALLOCATE (smat(n, n), tmat(n, n), v(n, 1), ipiv(n)) 2609 2610 CALL sg_overlap(smat, l, a, a) 2611 CALL sg_overlap(tmat, l, a, b) 2612 v(:, 1) = MATMUL(tmat, cb) 2613 CALL lapack_sgesv(n, 1, smat, n, ipiv, v, n, info) 2614 CPASSERT(info == 0) 2615 ca(:) = v(:, 1) 2616 2617 DEALLOCATE (smat, tmat, v, ipiv) 2618 2619 END SUBROUTINE project_function_a 2620 2621! ************************************************************************************************** 2622!> \brief ... 2623!> \param ca ... 2624!> \param a ... 2625!> \param bfun ... 2626!> \param grid_atom ... 2627!> \param l ... 2628! ************************************************************************************************** 2629 SUBROUTINE project_function_b(ca, a, bfun, grid_atom, l) 2630 ! project function f on ca 2631 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ca 2632 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, bfun 2633 TYPE(grid_atom_type), POINTER :: grid_atom 2634 INTEGER, INTENT(IN) :: l 2635 2636 CHARACTER(len=*), PARAMETER :: routineN = 'project_function_b', & 2637 routineP = moduleN//':'//routineN 2638 2639 INTEGER :: i, info, n, nr 2640 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv 2641 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: afun 2642 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smat, v 2643 2644 n = SIZE(ca) 2645 nr = grid_atom%nr 2646 ALLOCATE (smat(n, n), v(n, 1), ipiv(n), afun(nr)) 2647 2648 CALL sg_overlap(smat, l, a, a) 2649 DO i = 1, n 2650 afun(:) = grid_atom%rad(:)**l*EXP(-a(i)*grid_atom%rad2(:)) 2651 v(i, 1) = SUM(afun(:)*bfun(:)*grid_atom%wr(:)) 2652 END DO 2653 CALL lapack_sgesv(n, 1, smat, n, ipiv, v, n, info) 2654 CPASSERT(info == 0) 2655 ca(:) = v(:, 1) 2656 2657 DEALLOCATE (smat, v, ipiv, afun) 2658 2659 END SUBROUTINE project_function_b 2660 2661! ************************************************************************************************** 2662!> \brief Performs printing of cube files from local energy 2663!> \param input input 2664!> \param logger the logger 2665!> \param qs_env the qs_env in which the qs_env lives 2666!> \par History 2667!> 07.2019 created 2668!> \author JGH 2669! ************************************************************************************************** 2670 SUBROUTINE qs_scf_post_local_energy(input, logger, qs_env) 2671 TYPE(section_vals_type), POINTER :: input 2672 TYPE(cp_logger_type), POINTER :: logger 2673 TYPE(qs_environment_type), POINTER :: qs_env 2674 2675 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_local_energy', & 2676 routineP = moduleN//':'//routineN 2677 2678 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube 2679 INTEGER :: handle, io_unit, natom, unit_nr 2680 LOGICAL :: append_cube, gapw, gapw_xc, mpi_io 2681 TYPE(dft_control_type), POINTER :: dft_control 2682 TYPE(particle_list_type), POINTER :: particles 2683 TYPE(pw_env_type), POINTER :: pw_env 2684 TYPE(pw_p_type) :: eden 2685 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 2686 TYPE(qs_subsys_type), POINTER :: subsys 2687 TYPE(section_vals_type), POINTER :: dft_section 2688 2689 CALL timeset(routineN, handle) 2690 io_unit = cp_logger_get_default_io_unit(logger) 2691 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 2692 "DFT%PRINT%LOCAL_ENERGY_CUBE"), cp_p_file)) THEN 2693 dft_section => section_vals_get_subs_vals(input, "DFT") 2694 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, natom=natom) 2695 gapw = dft_control%qs_control%gapw 2696 gapw_xc = dft_control%qs_control%gapw_xc 2697 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys) 2698 CALL qs_subsys_get(subsys, particles=particles) 2699 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 2700 CALL pw_pool_create_pw(auxbas_pw_pool, eden%pw, use_data=REALDATA3D, in_space=REALSPACE) 2701 ! 2702 CALL qs_local_energy(qs_env, eden) 2703 ! 2704 append_cube = section_get_lval(input, "DFT%PRINT%LOCAL_ENERGY_CUBE%APPEND") 2705 IF (append_cube) THEN 2706 my_pos_cube = "APPEND" 2707 ELSE 2708 my_pos_cube = "REWIND" 2709 END IF 2710 mpi_io = .TRUE. 2711 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%LOCAL_ENERGY_CUBE", & 2712 extension=".cube", middle_name="local_energy", & 2713 file_position=my_pos_cube, mpi_io=mpi_io) 2714 CALL cp_pw_to_cube(eden%pw, & 2715 unit_nr, "LOCAL ENERGY", particles=particles, & 2716 stride=section_get_ivals(dft_section, & 2717 "PRINT%LOCAL_ENERGY_CUBE%STRIDE"), mpi_io=mpi_io) 2718 IF (io_unit > 0) THEN 2719 INQUIRE (UNIT=unit_nr, NAME=filename) 2720 IF (gapw .OR. gapw_xc) THEN 2721 WRITE (UNIT=io_unit, FMT="(/,T3,A,A)") & 2722 "The soft part of the local energy is written to the file: ", TRIM(ADJUSTL(filename)) 2723 ELSE 2724 WRITE (UNIT=io_unit, FMT="(/,T3,A,A)") & 2725 "The local energy is written to the file: ", TRIM(ADJUSTL(filename)) 2726 END IF 2727 END IF 2728 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2729 "DFT%PRINT%LOCAL_ENERGY_CUBE", mpi_io=mpi_io) 2730 ! 2731 CALL pw_pool_give_back_pw(auxbas_pw_pool, eden%pw) 2732 END IF 2733 CALL timestop(handle) 2734 2735 END SUBROUTINE qs_scf_post_local_energy 2736 2737! ************************************************************************************************** 2738!> \brief Performs printing of cube files from local energy 2739!> \param input input 2740!> \param logger the logger 2741!> \param qs_env the qs_env in which the qs_env lives 2742!> \par History 2743!> 07.2019 created 2744!> \author JGH 2745! ************************************************************************************************** 2746 SUBROUTINE qs_scf_post_local_stress(input, logger, qs_env) 2747 TYPE(section_vals_type), POINTER :: input 2748 TYPE(cp_logger_type), POINTER :: logger 2749 TYPE(qs_environment_type), POINTER :: qs_env 2750 2751 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_local_stress', & 2752 routineP = moduleN//':'//routineN 2753 2754 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube 2755 INTEGER :: handle, io_unit, natom, unit_nr 2756 LOGICAL :: append_cube, gapw, gapw_xc, mpi_io 2757 REAL(KIND=dp) :: beta 2758 TYPE(dft_control_type), POINTER :: dft_control 2759 TYPE(particle_list_type), POINTER :: particles 2760 TYPE(pw_env_type), POINTER :: pw_env 2761 TYPE(pw_p_type) :: stress 2762 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 2763 TYPE(qs_subsys_type), POINTER :: subsys 2764 TYPE(section_vals_type), POINTER :: dft_section 2765 2766 CALL timeset(routineN, handle) 2767 io_unit = cp_logger_get_default_io_unit(logger) 2768 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, & 2769 "DFT%PRINT%LOCAL_STRESS_CUBE"), cp_p_file)) THEN 2770 dft_section => section_vals_get_subs_vals(input, "DFT") 2771 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, natom=natom) 2772 gapw = dft_control%qs_control%gapw 2773 gapw_xc = dft_control%qs_control%gapw_xc 2774 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys) 2775 CALL qs_subsys_get(subsys, particles=particles) 2776 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 2777 CALL pw_pool_create_pw(auxbas_pw_pool, stress%pw, use_data=REALDATA3D, in_space=REALSPACE) 2778 ! 2779 ! use beta=0: kinetic energy density in symmetric form 2780 beta = 0.0_dp 2781 CALL qs_local_stress(qs_env, pressure=stress, beta=beta) 2782 ! 2783 append_cube = section_get_lval(input, "DFT%PRINT%LOCAL_STRESS_CUBE%APPEND") 2784 IF (append_cube) THEN 2785 my_pos_cube = "APPEND" 2786 ELSE 2787 my_pos_cube = "REWIND" 2788 END IF 2789 mpi_io = .TRUE. 2790 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%LOCAL_STRESS_CUBE", & 2791 extension=".cube", middle_name="local_stress", & 2792 file_position=my_pos_cube, mpi_io=mpi_io) 2793 CALL cp_pw_to_cube(stress%pw, & 2794 unit_nr, "LOCAL STRESS", particles=particles, & 2795 stride=section_get_ivals(dft_section, & 2796 "PRINT%LOCAL_STRESS_CUBE%STRIDE"), mpi_io=mpi_io) 2797 IF (io_unit > 0) THEN 2798 INQUIRE (UNIT=unit_nr, NAME=filename) 2799 WRITE (UNIT=io_unit, FMT="(/,T3,A)") "Write 1/3*Tr(sigma) to cube file" 2800 IF (gapw .OR. gapw_xc) THEN 2801 WRITE (UNIT=io_unit, FMT="(T3,A,A)") & 2802 "The soft part of the local stress is written to the file: ", TRIM(ADJUSTL(filename)) 2803 ELSE 2804 WRITE (UNIT=io_unit, FMT="(T3,A,A)") & 2805 "The local stress is written to the file: ", TRIM(ADJUSTL(filename)) 2806 END IF 2807 END IF 2808 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2809 "DFT%PRINT%LOCAL_STRESS_CUBE", mpi_io=mpi_io) 2810 ! 2811 CALL pw_pool_give_back_pw(auxbas_pw_pool, stress%pw) 2812 END IF 2813 2814 CALL timestop(handle) 2815 2816 END SUBROUTINE qs_scf_post_local_stress 2817 2818! ************************************************************************************************** 2819!> \brief Performs printing of cube files related to the implicit Poisson solver 2820!> \param input input 2821!> \param logger the logger 2822!> \param qs_env the qs_env in which the qs_env lives 2823!> \par History 2824!> 03.2016 refactored from write_mo_free_results [Hossein Bani-Hashemian] 2825!> \author Mohammad Hossein Bani-Hashemian 2826! ************************************************************************************************** 2827 SUBROUTINE qs_scf_post_ps_implicit(input, logger, qs_env) 2828 TYPE(section_vals_type), POINTER :: input 2829 TYPE(cp_logger_type), POINTER :: logger 2830 TYPE(qs_environment_type), POINTER :: qs_env 2831 2832 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_ps_implicit', & 2833 routineP = moduleN//':'//routineN 2834 2835 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube 2836 INTEGER :: boundary_condition, handle, i, j, & 2837 n_cstr, n_tiles, unit_nr 2838 LOGICAL :: append_cube, do_cstr_charge_cube, do_dielectric_cube, do_dirichlet_bc_cube, & 2839 has_dirichlet_bc, has_implicit_ps, mpi_io, tile_cubes 2840 TYPE(particle_list_type), POINTER :: particles 2841 TYPE(pw_env_type), POINTER :: pw_env 2842 TYPE(pw_p_type) :: aux_r 2843 TYPE(pw_poisson_type), POINTER :: poisson_env 2844 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 2845 TYPE(pw_type), POINTER :: dirichlet_tile 2846 TYPE(qs_subsys_type), POINTER :: subsys 2847 TYPE(section_vals_type), POINTER :: dft_section 2848 2849 CALL timeset(routineN, handle) 2850 2851 NULLIFY (pw_env, auxbas_pw_pool, dft_section, particles) 2852 2853 dft_section => section_vals_get_subs_vals(input, "DFT") 2854 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys) 2855 CALL qs_subsys_get(subsys, particles=particles) 2856 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 2857 2858 has_implicit_ps = .FALSE. 2859 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 2860 IF (pw_env%poisson_env%parameters%solver .EQ. pw_poisson_implicit) has_implicit_ps = .TRUE. 2861 2862 ! Write the dielectric constant into a cube file 2863 do_dielectric_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, & 2864 "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE"), cp_p_file) 2865 IF (has_implicit_ps .AND. do_dielectric_cube) THEN 2866 append_cube = section_get_lval(input, "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%APPEND") 2867 my_pos_cube = "REWIND" 2868 IF (append_cube) THEN 2869 my_pos_cube = "APPEND" 2870 END IF 2871 mpi_io = .TRUE. 2872 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE", & 2873 extension=".cube", middle_name="DIELECTRIC_CONSTANT", file_position=my_pos_cube, & 2874 mpi_io=mpi_io) 2875 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool) 2876 CALL pw_pool_create_pw(auxbas_pw_pool, aux_r%pw, use_data=REALDATA3D, in_space=REALSPACE) 2877 2878 boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition 2879 SELECT CASE (boundary_condition) 2880 CASE (PERIODIC_BC, MIXED_PERIODIC_BC) 2881 CALL pw_copy(poisson_env%implicit_env%dielectric%eps, aux_r%pw) 2882 CASE (MIXED_BC, NEUMANN_BC) 2883 CALL pw_shrink(pw_env%poisson_env%parameters%ps_implicit_params%neumann_directions, & 2884 pw_env%poisson_env%implicit_env%dct_env%dests_shrink, & 2885 pw_env%poisson_env%implicit_env%dct_env%srcs_shrink, & 2886 pw_env%poisson_env%implicit_env%dct_env%bounds_local_shftd, & 2887 poisson_env%implicit_env%dielectric%eps, aux_r%pw) 2888 END SELECT 2889 2890 CALL cp_pw_to_cube(aux_r%pw, unit_nr, "DIELECTRIC CONSTANT", particles=particles, & 2891 stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%STRIDE"), & 2892 mpi_io=mpi_io) 2893 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2894 "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE", mpi_io=mpi_io) 2895 2896 CALL pw_pool_give_back_pw(auxbas_pw_pool, aux_r%pw) 2897 ENDIF 2898 2899 ! Write Dirichlet constraint charges into a cube file 2900 do_cstr_charge_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, & 2901 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE"), cp_p_file) 2902 2903 has_dirichlet_bc = .FALSE. 2904 IF (has_implicit_ps) THEN 2905 boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition 2906 IF (boundary_condition .EQ. MIXED_PERIODIC_BC .OR. boundary_condition .EQ. MIXED_BC) THEN 2907 has_dirichlet_bc = .TRUE. 2908 END IF 2909 END IF 2910 2911 IF (has_implicit_ps .AND. do_cstr_charge_cube .AND. has_dirichlet_bc) THEN 2912 append_cube = section_get_lval(input, & 2913 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%APPEND") 2914 my_pos_cube = "REWIND" 2915 IF (append_cube) THEN 2916 my_pos_cube = "APPEND" 2917 END IF 2918 mpi_io = .TRUE. 2919 unit_nr = cp_print_key_unit_nr(logger, input, & 2920 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE", & 2921 extension=".cube", middle_name="dirichlet_cstr_charge", file_position=my_pos_cube, & 2922 mpi_io=mpi_io) 2923 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool) 2924 CALL pw_pool_create_pw(auxbas_pw_pool, aux_r%pw, use_data=REALDATA3D, in_space=REALSPACE) 2925 2926 boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition 2927 SELECT CASE (boundary_condition) 2928 CASE (MIXED_PERIODIC_BC) 2929 CALL pw_copy(poisson_env%implicit_env%cstr_charge, aux_r%pw) 2930 CASE (MIXED_BC) 2931 CALL pw_shrink(pw_env%poisson_env%parameters%ps_implicit_params%neumann_directions, & 2932 pw_env%poisson_env%implicit_env%dct_env%dests_shrink, & 2933 pw_env%poisson_env%implicit_env%dct_env%srcs_shrink, & 2934 pw_env%poisson_env%implicit_env%dct_env%bounds_local_shftd, & 2935 poisson_env%implicit_env%cstr_charge, aux_r%pw) 2936 END SELECT 2937 2938 CALL cp_pw_to_cube(aux_r%pw, unit_nr, "DIRICHLET CONSTRAINT CHARGE", particles=particles, & 2939 stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%STRIDE"), & 2940 mpi_io=mpi_io) 2941 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2942 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE", mpi_io=mpi_io) 2943 2944 CALL pw_pool_give_back_pw(auxbas_pw_pool, aux_r%pw) 2945 ENDIF 2946 2947 ! Write Dirichlet type constranits into cube files 2948 do_dirichlet_bc_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, & 2949 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE"), cp_p_file) 2950 has_dirichlet_bc = .FALSE. 2951 IF (has_implicit_ps) THEN 2952 boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition 2953 IF (boundary_condition .EQ. MIXED_PERIODIC_BC .OR. boundary_condition .EQ. MIXED_BC) THEN 2954 has_dirichlet_bc = .TRUE. 2955 END IF 2956 END IF 2957 2958 IF (has_implicit_ps .AND. has_dirichlet_bc .AND. do_dirichlet_bc_cube) THEN 2959 append_cube = section_get_lval(input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%APPEND") 2960 my_pos_cube = "REWIND" 2961 IF (append_cube) THEN 2962 my_pos_cube = "APPEND" 2963 END IF 2964 tile_cubes = section_get_lval(input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%TILE_CUBES") 2965 2966 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool) 2967 CALL pw_pool_create_pw(auxbas_pw_pool, aux_r%pw, use_data=REALDATA3D, in_space=REALSPACE) 2968 CALL pw_zero(aux_r%pw) 2969 2970 IF (tile_cubes) THEN 2971 ! one cube file per tile 2972 n_cstr = SIZE(poisson_env%implicit_env%contacts) 2973 DO j = 1, n_cstr 2974 n_tiles = poisson_env%implicit_env%contacts(j)%dirichlet_bc%n_tiles 2975 DO i = 1, n_tiles 2976 filename = "dirichlet_cstr_"//TRIM(ADJUSTL(cp_to_string(j)))// & 2977 "_tile_"//TRIM(ADJUSTL(cp_to_string(i))) 2978 mpi_io = .TRUE. 2979 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", & 2980 extension=".cube", middle_name=filename, file_position=my_pos_cube, & 2981 mpi_io=mpi_io) 2982 2983 CALL pw_copy(poisson_env%implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, aux_r%pw) 2984 2985 CALL cp_pw_to_cube(aux_r%pw, unit_nr, "DIRICHLET TYPE CONSTRAINT", particles=particles, & 2986 stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%STRIDE"), & 2987 mpi_io=mpi_io) 2988 CALL cp_print_key_finished_output(unit_nr, logger, input, & 2989 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", mpi_io=mpi_io) 2990 END DO 2991 END DO 2992 ELSE 2993 ! a single cube file 2994 NULLIFY (dirichlet_tile) 2995 CALL pw_pool_create_pw(auxbas_pw_pool, dirichlet_tile, use_data=REALDATA3D, in_space=REALSPACE) 2996 CALL pw_zero(dirichlet_tile) 2997 mpi_io = .TRUE. 2998 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", & 2999 extension=".cube", middle_name="DIRICHLET_CSTR", file_position=my_pos_cube, & 3000 mpi_io=mpi_io) 3001 3002 n_cstr = SIZE(poisson_env%implicit_env%contacts) 3003 DO j = 1, n_cstr 3004 n_tiles = poisson_env%implicit_env%contacts(j)%dirichlet_bc%n_tiles 3005 DO i = 1, n_tiles 3006 CALL pw_copy(poisson_env%implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, dirichlet_tile) 3007 CALL pw_axpy(dirichlet_tile, aux_r%pw) 3008 END DO 3009 END DO 3010 3011 CALL cp_pw_to_cube(aux_r%pw, unit_nr, "DIRICHLET TYPE CONSTRAINT", particles=particles, & 3012 stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%STRIDE"), & 3013 mpi_io=mpi_io) 3014 CALL cp_print_key_finished_output(unit_nr, logger, input, & 3015 "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", mpi_io=mpi_io) 3016 CALL pw_pool_give_back_pw(auxbas_pw_pool, dirichlet_tile) 3017 END IF 3018 3019 CALL pw_pool_give_back_pw(auxbas_pw_pool, aux_r%pw) 3020 ENDIF 3021 3022 CALL timestop(handle) 3023 3024 END SUBROUTINE qs_scf_post_ps_implicit 3025 3026!************************************************************************************************** 3027!> \brief writing the KS matrix in csr format into a file 3028!> \param qs_env qs environment 3029!> \param input the input 3030!> \author Mohammad Hossein Bani-Hashemian 3031! ************************************************************************************************** 3032 SUBROUTINE write_ks_matrix_csr(qs_env, input) 3033 TYPE(qs_environment_type), POINTER :: qs_env 3034 TYPE(section_vals_type), POINTER :: input 3035 3036 CHARACTER(len=*), PARAMETER :: routineN = 'write_ks_matrix_csr', & 3037 routineP = moduleN//':'//routineN 3038 3039 CHARACTER(LEN=default_path_length) :: file_name, fileformat 3040 INTEGER :: handle, ispin, output_unit, unit_nr 3041 LOGICAL :: bin, do_kpoints, do_ks_csr_write, uptr 3042 REAL(KIND=dp) :: thld 3043 TYPE(cp_logger_type), POINTER :: logger 3044 TYPE(dbcsr_csr_type) :: ks_mat_csr 3045 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks 3046 TYPE(dbcsr_type) :: matrix_ks_nosym 3047 TYPE(section_vals_type), POINTER :: dft_section 3048 3049 CALL timeset(routineN, handle) 3050 3051 NULLIFY (dft_section) 3052 3053 logger => cp_get_default_logger() 3054 output_unit = cp_logger_get_default_io_unit(logger) 3055 3056 dft_section => section_vals_get_subs_vals(input, "DFT") 3057 do_ks_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, & 3058 "PRINT%KS_CSR_WRITE"), cp_p_file) 3059 3060 ! NOTE: k-points has to be treated differently later. k-points has KS matrix as double pointer. 3061 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints) 3062 3063 IF (do_ks_csr_write .AND. (.NOT. do_kpoints)) THEN 3064 CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%THRESHOLD", r_val=thld) 3065 CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%UPPER_TRIANGULAR", l_val=uptr) 3066 CALL section_vals_val_get(dft_section, "PRINT%KS_CSR_WRITE%BINARY", l_val=bin) 3067 CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks) 3068 3069 IF (bin) THEN 3070 fileformat = "UNFORMATTED" 3071 ELSE 3072 fileformat = "FORMATTED" 3073 END IF 3074 3075 DO ispin = 1, SIZE(matrix_ks) 3076 3077 IF (dbcsr_has_symmetry(matrix_ks(ispin)%matrix)) THEN 3078 CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, matrix_ks_nosym) 3079 ELSE 3080 CALL dbcsr_copy(matrix_ks_nosym, matrix_ks(ispin)%matrix) 3081 END IF 3082 3083 CALL dbcsr_csr_create_from_dbcsr(matrix_ks_nosym, ks_mat_csr, dbcsr_csr_dbcsr_blkrow_dist) 3084 CALL dbcsr_convert_dbcsr_to_csr(matrix_ks_nosym, ks_mat_csr) 3085 3086 WRITE (file_name, '(A,I0)') "KS_SPIN_", ispin 3087 unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%KS_CSR_WRITE", & 3088 extension=".csr", middle_name=TRIM(file_name), & 3089 file_status="REPLACE", file_form=fileformat) 3090 CALL dbcsr_csr_write(ks_mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin) 3091 3092 CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%KS_CSR_WRITE") 3093 3094 CALL dbcsr_csr_destroy(ks_mat_csr) 3095 CALL dbcsr_release(matrix_ks_nosym) 3096 END DO 3097 END IF 3098 3099 CALL timestop(handle) 3100 3101 END SUBROUTINE write_ks_matrix_csr 3102 3103!************************************************************************************************** 3104!> \brief writing the KS matrix in csr format into a file 3105!> \param qs_env qs environment 3106!> \param input the input 3107!> \author Mohammad Hossein Bani-Hashemian 3108! ************************************************************************************************** 3109 SUBROUTINE write_s_matrix_csr(qs_env, input) 3110 TYPE(qs_environment_type), POINTER :: qs_env 3111 TYPE(section_vals_type), POINTER :: input 3112 3113 CHARACTER(len=*), PARAMETER :: routineN = 'write_s_matrix_csr', & 3114 routineP = moduleN//':'//routineN 3115 3116 CHARACTER(LEN=default_path_length) :: file_name, fileformat 3117 INTEGER :: handle, ispin, output_unit, unit_nr 3118 LOGICAL :: bin, do_kpoints, do_s_csr_write, uptr 3119 REAL(KIND=dp) :: thld 3120 TYPE(cp_logger_type), POINTER :: logger 3121 TYPE(dbcsr_csr_type) :: s_mat_csr 3122 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 3123 TYPE(dbcsr_type) :: matrix_s_nosym 3124 TYPE(section_vals_type), POINTER :: dft_section 3125 3126 CALL timeset(routineN, handle) 3127 3128 NULLIFY (dft_section) 3129 3130 logger => cp_get_default_logger() 3131 output_unit = cp_logger_get_default_io_unit(logger) 3132 3133 dft_section => section_vals_get_subs_vals(input, "DFT") 3134 do_s_csr_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, & 3135 "PRINT%S_CSR_WRITE"), cp_p_file) 3136 3137 ! NOTE: k-points has to be treated differently later. k-points has overlap matrix as double pointer. 3138 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints) 3139 3140 IF (do_s_csr_write .AND. (.NOT. do_kpoints)) THEN 3141 CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%THRESHOLD", r_val=thld) 3142 CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%UPPER_TRIANGULAR", l_val=uptr) 3143 CALL section_vals_val_get(dft_section, "PRINT%S_CSR_WRITE%BINARY", l_val=bin) 3144 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s) 3145 3146 IF (bin) THEN 3147 fileformat = "UNFORMATTED" 3148 ELSE 3149 fileformat = "FORMATTED" 3150 END IF 3151 3152 DO ispin = 1, SIZE(matrix_s) 3153 3154 IF (dbcsr_has_symmetry(matrix_s(ispin)%matrix)) THEN 3155 CALL dbcsr_desymmetrize(matrix_s(ispin)%matrix, matrix_s_nosym) 3156 ELSE 3157 CALL dbcsr_copy(matrix_s_nosym, matrix_s(ispin)%matrix) 3158 END IF 3159 3160 CALL dbcsr_csr_create_from_dbcsr(matrix_s_nosym, s_mat_csr, dbcsr_csr_dbcsr_blkrow_dist) 3161 CALL dbcsr_convert_dbcsr_to_csr(matrix_s_nosym, s_mat_csr) 3162 3163 WRITE (file_name, '(A,I0)') "S_SPIN_", ispin 3164 unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%S_CSR_WRITE", & 3165 extension=".csr", middle_name=TRIM(file_name), & 3166 file_status="REPLACE", file_form=fileformat) 3167 CALL dbcsr_csr_write(s_mat_csr, unit_nr, upper_triangle=uptr, threshold=thld, binary=bin) 3168 3169 CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%S_CSR_WRITE") 3170 3171 CALL dbcsr_csr_destroy(s_mat_csr) 3172 CALL dbcsr_release(matrix_s_nosym) 3173 END DO 3174 END IF 3175 3176 CALL timestop(handle) 3177 3178 END SUBROUTINE write_s_matrix_csr 3179 3180!************************************************************************************************** 3181!> \brief write an adjacency (interaction) matrix 3182!> \param qs_env qs environment 3183!> \param input the input 3184!> \author Mohammad Hossein Bani-Hashemian 3185! ************************************************************************************************** 3186 SUBROUTINE write_adjacency_matrix(qs_env, input) 3187 TYPE(qs_environment_type), POINTER :: qs_env 3188 TYPE(section_vals_type), POINTER :: input 3189 3190 CHARACTER(len=*), PARAMETER :: routineN = 'write_adjacency_matrix', & 3191 routineP = moduleN//':'//routineN 3192 3193 INTEGER :: adjm_size, colind, handle, iatom, ikind, & 3194 ind, jatom, jkind, k, natom, nkind, & 3195 output_unit, rowind, unit_nr 3196 INTEGER, ALLOCATABLE, DIMENSION(:) :: interact_adjm 3197 LOGICAL :: do_adjm_write, do_symmetric 3198 TYPE(cp_logger_type), POINTER :: logger 3199 TYPE(cp_para_env_type), POINTER :: para_env 3200 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b 3201 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 3202 TYPE(neighbor_list_iterator_p_type), & 3203 DIMENSION(:), POINTER :: nl_iterator 3204 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 3205 POINTER :: nl 3206 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 3207 TYPE(section_vals_type), POINTER :: dft_section 3208 3209 CALL timeset(routineN, handle) 3210 3211 NULLIFY (dft_section) 3212 3213 logger => cp_get_default_logger() 3214 output_unit = cp_logger_get_default_io_unit(logger) 3215 3216 dft_section => section_vals_get_subs_vals(input, "DFT") 3217 do_adjm_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, & 3218 "PRINT%ADJMAT_WRITE"), cp_p_file) 3219 3220 IF (do_adjm_write) THEN 3221 NULLIFY (qs_kind_set, nl_iterator) 3222 NULLIFY (basis_set_list_a, basis_set_list_b, basis_set_a, basis_set_b) 3223 3224 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, sab_orb=nl, natom=natom, para_env=para_env) 3225 3226 nkind = SIZE(qs_kind_set) 3227 CPASSERT(SIZE(nl) .GT. 0) 3228 CALL get_neighbor_list_set_p(neighbor_list_sets=nl, symmetric=do_symmetric) 3229 CPASSERT(do_symmetric) 3230 ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind)) 3231 CALL basis_set_list_setup(basis_set_list_a, "ORB", qs_kind_set) 3232 CALL basis_set_list_setup(basis_set_list_b, "ORB", qs_kind_set) 3233 3234 adjm_size = ((natom + 1)*natom)/2 3235 ALLOCATE (interact_adjm(4*adjm_size)) 3236 interact_adjm = 0 3237 3238 NULLIFY (nl_iterator) 3239 CALL neighbor_list_iterator_create(nl_iterator, nl) 3240 DO WHILE (neighbor_list_iterate(nl_iterator) .EQ. 0) 3241 CALL get_iterator_info(nl_iterator, & 3242 ikind=ikind, jkind=jkind, & 3243 iatom=iatom, jatom=jatom) 3244 3245 basis_set_a => basis_set_list_a(ikind)%gto_basis_set 3246 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 3247 basis_set_b => basis_set_list_b(jkind)%gto_basis_set 3248 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 3249 3250 ! move everything to the upper triangular part 3251 IF (iatom .LE. jatom) THEN 3252 rowind = iatom 3253 colind = jatom 3254 ELSE 3255 rowind = jatom 3256 colind = iatom 3257 ! swap the kinds too 3258 ikind = ikind + jkind 3259 jkind = ikind - jkind 3260 ikind = ikind - jkind 3261 END IF 3262 3263 ! indexing upper triangular matrix 3264 ind = adjm_size - (natom - rowind + 1)*((natom - rowind + 1) + 1)/2 + colind - rowind + 1 3265 ! convert the upper triangular matrix into a adjm_size x 4 matrix 3266 ! columns are: iatom, jatom, ikind, jkind 3267 interact_adjm((ind - 1)*4 + 1) = rowind 3268 interact_adjm((ind - 1)*4 + 2) = colind 3269 interact_adjm((ind - 1)*4 + 3) = ikind 3270 interact_adjm((ind - 1)*4 + 4) = jkind 3271 ENDDO 3272 3273 CALL mp_sum(interact_adjm, para_env%group) 3274 3275 unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%ADJMAT_WRITE", & 3276 extension=".adjmat", file_form="FORMATTED", & 3277 file_status="REPLACE") 3278 IF (unit_nr .GT. 0) THEN 3279 WRITE (unit_nr, "(1A,2X,1A,5X,1A,4X,A5,3X,A5)") "#", "iatom", "jatom", "ikind", "jkind" 3280 DO k = 1, 4*adjm_size, 4 3281 ! print only the interacting atoms 3282 IF (interact_adjm(k) .GT. 0 .AND. interact_adjm(k + 1) .GT. 0) THEN 3283 WRITE (unit_nr, "(I8,2X,I8,3X,I6,2X,I6)") interact_adjm(k:k + 3) 3284 END IF 3285 END DO 3286 END IF 3287 3288 CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%ADJMAT_WRITE") 3289 3290 CALL neighbor_list_iterator_release(nl_iterator) 3291 DEALLOCATE (basis_set_list_a, basis_set_list_b) 3292 END IF 3293 3294 CALL timestop(handle) 3295 3296 END SUBROUTINE write_adjacency_matrix 3297 3298! ************************************************************************************************** 3299!> \brief Updates Hartree potential with MP2 density. Important for REPEAT charges 3300!> \param rho ... 3301!> \param qs_env ... 3302!> \author Vladimir Rybkin 3303! ************************************************************************************************** 3304 SUBROUTINE update_hartree_with_mp2(rho, qs_env) 3305 TYPE(qs_rho_type), POINTER :: rho 3306 TYPE(qs_environment_type), POINTER :: qs_env 3307 3308 LOGICAL :: use_virial 3309 TYPE(pw_env_type), POINTER :: pw_env 3310 TYPE(pw_p_type) :: rho_tot_gspace, v_hartree_gspace, & 3311 v_hartree_rspace 3312 TYPE(pw_p_type), POINTER :: rho_core 3313 TYPE(pw_poisson_type), POINTER :: poisson_env 3314 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 3315 TYPE(qs_energy_type), POINTER :: energy 3316 TYPE(virial_type), POINTER :: virial 3317 3318 NULLIFY (auxbas_pw_pool, pw_env, poisson_env, energy, rho_core, v_hartree_rspace%pw, virial) 3319 CALL get_qs_env(qs_env, pw_env=pw_env, energy=energy, & 3320 rho_core=rho_core, virial=virial, & 3321 v_hartree_rspace=v_hartree_rspace%pw) 3322 3323 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 3324 3325 IF (.NOT. use_virial) THEN 3326 3327 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 3328 poisson_env=poisson_env) 3329 CALL pw_pool_create_pw(auxbas_pw_pool, & 3330 v_hartree_gspace%pw, & 3331 use_data=COMPLEXDATA1D, & 3332 in_space=RECIPROCALSPACE) 3333 CALL pw_pool_create_pw(auxbas_pw_pool, & 3334 rho_tot_gspace%pw, & 3335 use_data=COMPLEXDATA1D, & 3336 in_space=RECIPROCALSPACE) 3337 3338 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho) 3339 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, energy%hartree, & 3340 v_hartree_gspace%pw, rho_core=rho_core) 3341 3342 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw) 3343 CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol) 3344 3345 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw) 3346 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw) 3347 ENDIF 3348 3349 END SUBROUTINE update_hartree_with_mp2 3350 3351END MODULE qs_scf_post_gpw 3352