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