LCOV - code coverage report
Current view: top level - src - qs_scf_post_gpw.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:6d276e9) Lines: 90.5 % 1745 1579
Test Date: 2026-09-10 07:29:18 Functions: 97.2 % 36 35

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Does all kind of post scf calculations for GPW/GAPW
      10              : !> \par History
      11              : !>      Started as a copy from the relevant part of qs_scf
      12              : !>      Start to adapt for k-points [07.2015, JGH]
      13              : !> \author Joost VandeVondele (10.2003)
      14              : ! **************************************************************************************************
      15              : MODULE qs_scf_post_gpw
      16              :    USE admm_types,                      ONLY: admm_type
      17              :    USE admm_utils,                      ONLY: admm_correct_for_eigenvalues,&
      18              :                                               admm_uncorrect_for_eigenvalues
      19              :    USE ai_onecenter,                    ONLY: sg_overlap
      20              :    USE atom_kind_orbitals,              ONLY: calculate_atomic_density
      21              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      22              :                                               get_atomic_kind
      23              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      24              :                                               gto_basis_set_type
      25              :    USE casino_utils,                    ONLY: write_casino
      26              :    USE cell_types,                      ONLY: cell_type
      27              :    USE cp_array_utils,                  ONLY: cp_1d_r_p_type
      28              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      29              :    USE cp_control_types,                ONLY: dft_control_type,&
      30              :                                               rtp_control_type
      31              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      32              :                                               dbcsr_p_type,&
      33              :                                               dbcsr_type
      34              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_checksum
      35              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      36              :                                               dbcsr_deallocate_matrix_set
      37              :    USE cp_dbcsr_output,                 ONLY: cp_dbcsr_write_sparse_matrix
      38              :    USE cp_ddapc_util,                   ONLY: get_ddapc
      39              :    USE cp_fm_diag,                      ONLY: choose_eigv_solver
      40              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      41              :                                               cp_fm_struct_release,&
      42              :                                               cp_fm_struct_type
      43              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      44              :                                               cp_fm_get_info,&
      45              :                                               cp_fm_init_random,&
      46              :                                               cp_fm_release,&
      47              :                                               cp_fm_to_fm,&
      48              :                                               cp_fm_type
      49              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      50              :                                               cp_logger_get_default_io_unit,&
      51              :                                               cp_logger_type,&
      52              :                                               cp_to_string
      53              :    USE cp_output_handling,              ONLY: cp_iter_string,&
      54              :                                               cp_p_file,&
      55              :                                               cp_print_key_finished_output,&
      56              :                                               cp_print_key_should_output,&
      57              :                                               cp_print_key_unit_nr
      58              :    USE cp_output_handling_openpmd,      ONLY: cp_openpmd_close_iterations,&
      59              :                                               cp_openpmd_print_key_finished_output,&
      60              :                                               cp_openpmd_print_key_unit_nr
      61              :    USE cp_realspace_grid_cube,          ONLY: cp_pw_to_cube
      62              :    USE cp_realspace_grid_openpmd,       ONLY: cp_pw_to_openpmd
      63              :    USE cp_result_methods,               ONLY: get_results
      64              :    USE cp_result_types,                 ONLY: cp_result_type
      65              :    USE dct,                             ONLY: pw_shrink
      66              :    USE ed_analysis,                     ONLY: edmf_analysis
      67              :    USE eeq_method,                      ONLY: eeq_print
      68              :    USE et_coupling_types,               ONLY: set_et_coupling_type
      69              :    USE gapw_gspace_reconstruction,      ONLY: calculate_rhotot_elec_gspace
      70              :    USE hfx_ri,                          ONLY: print_ri_hfx
      71              :    USE hirshfeld_methods,               ONLY: comp_hirshfeld_charges,&
      72              :                                               comp_hirshfeld_i_charges,&
      73              :                                               create_shape_function,&
      74              :                                               save_hirshfeld_charges,&
      75              :                                               write_hirshfeld_charges
      76              :    USE hirshfeld_types,                 ONLY: create_hirshfeld_type,&
      77              :                                               hirshfeld_type,&
      78              :                                               release_hirshfeld_type,&
      79              :                                               set_hirshfeld_info
      80              :    USE iao_analysis,                    ONLY: iao_wfn_analysis
      81              :    USE iao_types,                       ONLY: iao_env_type,&
      82              :                                               iao_read_input
      83              :    USE input_constants,                 ONLY: &
      84              :         do_loc_both, do_loc_homo, do_loc_jacobi, do_loc_lumo, do_loc_mixed, do_loc_none, &
      85              :         moments_format_trajectory, ot_precond_full_all, radius_covalent, radius_user, &
      86              :         ref_charge_atomic, ref_charge_mulliken
      87              :    USE input_section_types,             ONLY: section_get_ival,&
      88              :                                               section_get_ivals,&
      89              :                                               section_get_lval,&
      90              :                                               section_get_rval,&
      91              :                                               section_vals_get,&
      92              :                                               section_vals_get_subs_vals,&
      93              :                                               section_vals_type,&
      94              :                                               section_vals_val_get
      95              :    USE kinds,                           ONLY: default_path_length,&
      96              :                                               default_string_length,&
      97              :                                               dp
      98              :    USE kpoint_mo_dump,                  ONLY: write_kpoint_mo_data
      99              :    USE kpoint_types,                    ONLY: kpoint_type
     100              :    USE localized_moments,               ONLY: calculate_kg_moments
     101              :    USE mao_wfn_analysis,                ONLY: mao_analysis
     102              :    USE mathconstants,                   ONLY: pi
     103              :    USE memory_utilities,                ONLY: reallocate
     104              :    USE message_passing,                 ONLY: mp_para_env_type
     105              :    USE minbas_wfn_analysis,             ONLY: minbas_analysis
     106              :    USE molden_utils,                    ONLY: write_mos_molden
     107              :    USE molecule_types,                  ONLY: molecule_type
     108              :    USE mulliken,                        ONLY: mulliken_charges
     109              :    USE orbital_pointers,                ONLY: indso
     110              :    USE particle_list_types,             ONLY: particle_list_type
     111              :    USE particle_types,                  ONLY: particle_type
     112              :    USE physcon,                         ONLY: a_bohr,&
     113              :                                               angstrom,&
     114              :                                               debye,&
     115              :                                               evolt
     116              :    USE population_analyses,             ONLY: lowdin_population_analysis,&
     117              :                                               mulliken_population_analysis
     118              :    USE preconditioner_types,            ONLY: preconditioner_type
     119              :    USE ps_implicit_types,               ONLY: MIXED_BC,&
     120              :                                               MIXED_PERIODIC_BC,&
     121              :                                               NEUMANN_BC,&
     122              :                                               PERIODIC_BC
     123              :    USE pw_env_types,                    ONLY: pw_env_get,&
     124              :                                               pw_env_type
     125              :    USE pw_grids,                        ONLY: get_pw_grid_info
     126              :    USE pw_methods,                      ONLY: pw_axpy,&
     127              :                                               pw_copy,&
     128              :                                               pw_derive,&
     129              :                                               pw_integrate_function,&
     130              :                                               pw_scale,&
     131              :                                               pw_transfer,&
     132              :                                               pw_zero
     133              :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
     134              :    USE pw_poisson_types,                ONLY: pw_poisson_implicit,&
     135              :                                               pw_poisson_type
     136              :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
     137              :                                               pw_pool_type
     138              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
     139              :                                               pw_r3d_rs_type
     140              :    USE qs_chargemol,                    ONLY: write_wfx
     141              :    USE qs_charges_types,                ONLY: qs_charges_type
     142              :    USE qs_collocate_density,            ONLY: calculate_rho_resp_all,&
     143              :                                               calculate_wavefunction
     144              :    USE qs_commutators,                  ONLY: build_com_hr_matrix
     145              :    USE qs_core_energies,                ONLY: calculate_ptrace
     146              :    USE qs_dos,                          ONLY: calculate_dos,&
     147              :                                               calculate_dos_kp
     148              :    USE qs_dos_utils,                    ONLY: get_dos_pdos_flags
     149              :    USE qs_electric_field_gradient,      ONLY: qs_efg_calc
     150              :    USE qs_elf_methods,                  ONLY: qs_elf_calc
     151              :    USE qs_energy_types,                 ONLY: qs_energy_type
     152              :    USE qs_energy_window,                ONLY: energy_windows
     153              :    USE qs_environment_types,            ONLY: get_qs_env,&
     154              :                                               qs_environment_type,&
     155              :                                               set_qs_env
     156              :    USE qs_epr_hyp,                      ONLY: qs_epr_hyp_calc
     157              :    USE qs_grid_atom,                    ONLY: grid_atom_type
     158              :    USE qs_integral_utils,               ONLY: basis_set_list_setup
     159              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
     160              :                                               qs_kind_type
     161              :    USE qs_ks_methods,                   ONLY: calc_rho_tot_gspace,&
     162              :                                               qs_ks_update_qs_env
     163              :    USE qs_ks_types,                     ONLY: qs_ks_did_change
     164              :    USE qs_kubo_transport,               ONLY: qs_scf_post_kubo_transport
     165              :    USE qs_loc_dipole,                   ONLY: loc_dipole
     166              :    USE qs_loc_states,                   ONLY: get_localization_info
     167              :    USE qs_loc_types,                    ONLY: qs_loc_env_create,&
     168              :                                               qs_loc_env_release,&
     169              :                                               qs_loc_env_type
     170              :    USE qs_loc_utils,                    ONLY: loc_write_restart,&
     171              :                                               qs_loc_control_init,&
     172              :                                               qs_loc_env_init,&
     173              :                                               qs_loc_init,&
     174              :                                               retain_history
     175              :    USE qs_local_properties,             ONLY: qs_local_energy,&
     176              :                                               qs_local_stress
     177              :    USE qs_mo_io,                        ONLY: write_dm_binary_restart
     178              :    USE qs_mo_methods,                   ONLY: calculate_subspace_eigenvalues,&
     179              :                                               make_mo_eig
     180              :    USE qs_mo_occupation,                ONLY: set_mo_occupation
     181              :    USE qs_mo_types,                     ONLY: get_mo_set,&
     182              :                                               mo_set_type
     183              :    USE qs_moments,                      ONLY: qs_moment_berry_phase,&
     184              :                                               qs_moment_kpoints,&
     185              :                                               qs_moment_locop
     186              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
     187              :                                               get_neighbor_list_set_p,&
     188              :                                               neighbor_list_iterate,&
     189              :                                               neighbor_list_iterator_create,&
     190              :                                               neighbor_list_iterator_p_type,&
     191              :                                               neighbor_list_iterator_release,&
     192              :                                               neighbor_list_set_p_type
     193              :    USE qs_ot_eigensolver,               ONLY: ot_eigensolver
     194              :    USE qs_pdos,                         ONLY: calculate_projected_dos,&
     195              :                                               calculate_projected_dos_kp
     196              :    USE qs_resp,                         ONLY: resp_fit
     197              :    USE qs_rho0_types,                   ONLY: get_rho0_mpole,&
     198              :                                               mpole_rho_atom,&
     199              :                                               rho0_mpole_type
     200              :    USE qs_rho_atom_types,               ONLY: rho_atom_type
     201              :    USE qs_rho_methods,                  ONLY: qs_rho_update_rho
     202              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
     203              :                                               qs_rho_type
     204              :    USE qs_scf_csr_write,                ONLY: write_hcore_matrix_csr,&
     205              :                                               write_ks_matrix_csr,&
     206              :                                               write_p_matrix_csr,&
     207              :                                               write_s_matrix_csr
     208              :    USE qs_scf_output,                   ONLY: qs_scf_write_mos
     209              :    USE qs_scf_types,                    ONLY: ot_method_nr,&
     210              :                                               qs_scf_env_type
     211              :    USE qs_scf_wfn_mix,                  ONLY: wfn_mix
     212              :    USE qs_subsys_types,                 ONLY: qs_subsys_get,&
     213              :                                               qs_subsys_type
     214              :    USE qs_wannier90,                    ONLY: wannier90_interface
     215              :    USE s_square_methods,                ONLY: compute_s_square
     216              :    USE scf_control_types,               ONLY: scf_control_type
     217              :    USE stm_images,                      ONLY: th_stm_image
     218              :    USE transport,                       ONLY: qs_scf_post_transport
     219              :    USE trexio_utils,                    ONLY: write_trexio
     220              :    USE virial_types,                    ONLY: virial_type
     221              :    USE voronoi_interface,               ONLY: entry_voronoi_or_bqb
     222              :    USE xray_diffraction,                ONLY: xray_diffraction_spectrum
     223              : #include "./base/base_uses.f90"
     224              : 
     225              :    IMPLICIT NONE
     226              :    PRIVATE
     227              : 
     228              :    ! Global parameters
     229              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_post_gpw'
     230              :    PUBLIC :: make_lumo_gpw, &
     231              :              scf_post_calculation_gpw, &
     232              :              qs_scf_post_moments, &
     233              :              write_mo_dependent_results, &
     234              :              write_mo_free_results
     235              : 
     236              :    CHARACTER(len=*), PARAMETER :: &
     237              :       str_mo_cubes = "PRINT%MO_CUBES", &
     238              :       str_mo_openpmd = "PRINT%MO_OPENPMD", &
     239              :       str_elf_cubes = "PRINT%ELF_CUBE", &
     240              :       str_elf_openpmd = "PRINT%ELF_OPENPMD", &
     241              :       str_e_density_cubes = "PRINT%E_DENSITY_CUBE", &
     242              :       str_e_density_openpmd = "PRINT%E_DENSITY_OPENPMD"
     243              : 
     244              :    INTEGER, PARAMETER :: grid_output_cubes = 1, grid_output_openpmd = 2
     245              : 
     246              :    REAL(kind=dp), DIMENSION(7), PARAMETER :: openpmd_unit_dimension_density = &
     247              :                                              [-3, 0, 0, 0, 0, 0, 0]
     248              :    REAL(kind=dp), DIMENSION(7), PARAMETER :: openpmd_unit_dimension_dimensionless = &
     249              :                                              [0, 0, 0, 0, 0, 0, 0]
     250              :    REAL(kind=dp), DIMENSION(7), PARAMETER :: openpmd_unit_dimension_wavefunction = &
     251              :                                              [-1.5_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp]
     252              :    REAL(kind=dp), PARAMETER :: openpmd_unit_si_density = a_bohr**(-3)
     253              :    REAL(kind=dp), PARAMETER :: openpmd_unit_si_dimensionless = 1.0_dp
     254              :    REAL(kind=dp), PARAMETER :: openpmd_unit_si_wavefunction = a_bohr**(-1.5_dp)
     255              : 
     256              :    ! Generic information on whether a certain output section has been activated
     257              :    ! or not, and on whether it has been activated in the Cube or openPMD variant.
     258              :    ! Create with function cube_or_openpmd(), see there for further details.
     259              :    TYPE cp_section_key
     260              :       CHARACTER(len=default_string_length) :: relative_section_key = "" ! e.g. PRINT%MO_CUBES
     261              :       CHARACTER(len=default_string_length) :: absolute_section_key = "" ! e.g. DFT%PRINT%MO_CUBES
     262              :       CHARACTER(len=7) :: format_name = "" ! 'openPMD' or 'Cube', for logging
     263              :       INTEGER :: grid_output = -1 ! either 1 for grid_output_cubes or 2 for grid_output_openpmd
     264              :       LOGICAL :: do_output = .FALSE.
     265              :    CONTAINS
     266              :       ! Open a file as either Cube or openPMD
     267              :       PROCEDURE, PUBLIC :: print_key_unit_nr => cp_forward_print_key_unit_nr
     268              :       ! Write either to the Cube or openPMD file
     269              :       PROCEDURE, PUBLIC :: write_pw => cp_forward_write_pw
     270              :       ! Close either the Cube or openPMD file
     271              :       PROCEDURE, PUBLIC :: print_key_finished_output => cp_forward_print_key_finished_output
     272              :       ! Helpers
     273              :       PROCEDURE, PUBLIC :: do_openpmd => cp_section_key_do_openpmd
     274              :       PROCEDURE, PUBLIC :: do_cubes => cp_section_key_do_cubes
     275              :       PROCEDURE, PUBLIC :: concat_to_relative => cp_section_key_concat_to_relative
     276              :       PROCEDURE, PUBLIC :: concat_to_absolute => cp_section_key_concat_to_absolute
     277              :    END TYPE cp_section_key
     278              : 
     279              : CONTAINS
     280              : 
     281              : ! **************************************************************************************************
     282              : !> \brief Collects the effective core charge for every atom in a QS environment
     283              : !> \param qs_env the QS environment
     284              : !> \param zcharge effective core charges ordered by atom index
     285              : ! **************************************************************************************************
     286        13809 :    SUBROUTINE get_effective_core_charges(qs_env, zcharge)
     287              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     288              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     289              :          INTENT(OUT)                                     :: zcharge
     290              : 
     291              :       INTEGER                                            :: iat, iatom, ikind, nat, natom, nkind
     292              :       REAL(KIND=dp)                                      :: zeff
     293        13809 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     294        13809 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     295              : 
     296              :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
     297        13809 :                       nkind=nkind, natom=natom)
     298        41427 :       ALLOCATE (zcharge(natom))
     299        38087 :       DO ikind = 1, nkind
     300        24278 :          CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
     301        24278 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
     302        89994 :          DO iatom = 1, nat
     303        51907 :             iat = atomic_kind_set(ikind)%atom_list(iatom)
     304        76185 :             zcharge(iat) = zeff
     305              :          END DO
     306              :       END DO
     307        13809 :    END SUBROUTINE get_effective_core_charges
     308              : 
     309              : ! **************************************************************************************************
     310              : !> \brief Append `extend_by` to the absolute path of the base section.
     311              : !> \param self ...
     312              : !> \param extend_by ...
     313              : !> \return ...
     314              : ! **************************************************************************************************
     315          310 :    FUNCTION cp_section_key_concat_to_absolute(self, extend_by) RESULT(res)
     316              :       CLASS(cp_section_key), INTENT(IN) :: self
     317              :       CHARACTER(*), INTENT(IN) :: extend_by
     318              :       CHARACTER(len=default_string_length) :: res
     319              : 
     320          310 :       IF (LEN(TRIM(extend_by)) > 0 .AND. extend_by(1:1) == "%") THEN
     321          310 :          res = TRIM(self%absolute_section_key)//TRIM(extend_by)
     322              :       ELSE
     323            0 :          res = TRIM(self%absolute_section_key)//"%"//TRIM(extend_by)
     324              :       END IF
     325          310 :    END FUNCTION cp_section_key_concat_to_absolute
     326              : 
     327              : ! **************************************************************************************************
     328              : !> \brief Append `extend_by` to the relative path (e.g. without DFT%) of the base section.
     329              : !> \param self ...
     330              : !> \param extend_by ...
     331              : !> \return ...
     332              : ! **************************************************************************************************
     333        25610 :    FUNCTION cp_section_key_concat_to_relative(self, extend_by) RESULT(res)
     334              :       CLASS(cp_section_key), INTENT(IN) :: self
     335              :       CHARACTER(*), INTENT(IN) :: extend_by
     336              :       CHARACTER(len=default_string_length) :: res
     337              : 
     338        25610 :       IF (LEN(TRIM(extend_by)) > 0 .AND. extend_by(1:1) == "%") THEN
     339        25610 :          res = TRIM(self%relative_section_key)//TRIM(extend_by)
     340              :       ELSE
     341            0 :          res = TRIM(self%relative_section_key)//"%"//TRIM(extend_by)
     342              :       END IF
     343        25610 :    END FUNCTION cp_section_key_concat_to_relative
     344              : 
     345              : ! **************************************************************************************************
     346              : !> \brief Is Cube output active for the current base section?
     347              : !> \param self ...
     348              : !> \return ...
     349              : ! **************************************************************************************************
     350          668 :    FUNCTION cp_section_key_do_cubes(self) RESULT(res)
     351              :       CLASS(cp_section_key) :: self
     352              :       LOGICAL :: res
     353              : 
     354          668 :       res = self%do_output .AND. self%grid_output == grid_output_cubes
     355          668 :    END FUNCTION cp_section_key_do_cubes
     356              : 
     357              : ! **************************************************************************************************
     358              : !> \brief Is openPMD output active for the current base section?
     359              : !> \param self ...
     360              : !> \return ...
     361              : ! **************************************************************************************************
     362          668 :    FUNCTION cp_section_key_do_openpmd(self) RESULT(res)
     363              :       CLASS(cp_section_key) :: self
     364              :       LOGICAL :: res
     365              : 
     366          668 :       res = self%do_output .AND. self%grid_output == grid_output_openpmd
     367          668 :    END FUNCTION cp_section_key_do_openpmd
     368              : 
     369              : ! **************************************************************************************************
     370              : !> \brief Forwards to either `cp_print_key_unit_nr` or `cp_openpmd_print_key_unit_nr`,
     371              : !>        depending on the configuration of the current base section.
     372              : !>        Opens either a Cube or openPMD output file
     373              : !> \param self ...
     374              : !> \param logger ...
     375              : !> \param basis_section ...
     376              : !> \param print_key_path ...
     377              : !> \param extension ...
     378              : !> \param middle_name ...
     379              : !> \param local ...
     380              : !> \param log_filename ...
     381              : !> \param ignore_should_output ...
     382              : !> \param file_form ...
     383              : !> \param file_position ...
     384              : !> \param file_action ...
     385              : !> \param file_status ...
     386              : !> \param do_backup ...
     387              : !> \param on_file ...
     388              : !> \param is_new_file ...
     389              : !> \param mpi_io ...
     390              : !> \param fout ...
     391              : !> \param openpmd_basename ...
     392              : !> \param openpmd_unit_dimension ...
     393              : !> \param openpmd_unit_si ...
     394              : !> \param sim_time ...
     395              : !> \return ...
     396              : ! **************************************************************************************************
     397          552 :    FUNCTION cp_forward_print_key_unit_nr( &
     398              :       self, &
     399              :       logger, &
     400              :       basis_section, &
     401              :       print_key_path, &
     402              :       extension, &
     403              :       middle_name, &
     404              :       local, &
     405              :       log_filename, &
     406              :       ignore_should_output, &
     407              :       file_form, &
     408              :       file_position, &
     409              :       file_action, &
     410              :       file_status, &
     411              :       do_backup, &
     412              :       on_file, &
     413              :       is_new_file, &
     414              :       mpi_io, &
     415              :       fout, &
     416              :       openpmd_basename, &
     417              :       openpmd_unit_dimension, &
     418              :       openpmd_unit_si, &
     419              :       sim_time) RESULT(res)
     420              : 
     421              :       CLASS(cp_section_key), INTENT(IN)                  :: self
     422              :       TYPE(cp_logger_type), POINTER                      :: logger
     423              :       TYPE(section_vals_type), INTENT(IN)                :: basis_section
     424              :       CHARACTER(len=*), INTENT(IN), OPTIONAL             :: print_key_path
     425              :       CHARACTER(len=*), INTENT(IN)                       :: extension
     426              :       CHARACTER(len=*), INTENT(IN), OPTIONAL             :: middle_name
     427              :       LOGICAL, INTENT(IN), OPTIONAL                      :: local, log_filename, ignore_should_output
     428              :       CHARACTER(len=*), INTENT(IN), OPTIONAL             :: file_form, file_position, file_action, &
     429              :                                                             file_status
     430              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_backup, on_file
     431              :       LOGICAL, INTENT(OUT), OPTIONAL                     :: is_new_file
     432              :       LOGICAL, INTENT(INOUT), OPTIONAL                   :: mpi_io
     433              :       CHARACTER(len=default_path_length), INTENT(OUT), &
     434              :          OPTIONAL                                        :: fout
     435              :       CHARACTER(len=*), INTENT(IN), OPTIONAL             :: openpmd_basename
     436              :       REAL(kind=dp), DIMENSION(7), OPTIONAL, INTENT(IN)  :: openpmd_unit_dimension
     437              :       REAL(kind=dp), OPTIONAL, INTENT(IN)                :: openpmd_unit_si
     438              :       REAL(kind=dp), OPTIONAL, INTENT(IN)                :: sim_time
     439              :       INTEGER                                            :: res
     440              : 
     441          552 :       IF (self%grid_output == grid_output_cubes) THEN
     442              :          res = cp_print_key_unit_nr( &
     443              :                logger, basis_section, print_key_path, extension=extension, &
     444              :                middle_name=middle_name, local=local, log_filename=log_filename, &
     445              :                ignore_should_output=ignore_should_output, file_form=file_form, &
     446              :                file_position=file_position, file_action=file_action, &
     447              :                file_status=file_status, do_backup=do_backup, on_file=on_file, &
     448         2466 :                is_new_file=is_new_file, mpi_io=mpi_io, fout=fout)
     449              :       ELSE
     450              :          res = cp_openpmd_print_key_unit_nr( &
     451              :                logger, &
     452              :                basis_section, &
     453              :                print_key_path, &
     454              :                middle_name=middle_name, &
     455              :                ignore_should_output=ignore_should_output, &
     456              :                mpi_io=mpi_io, &
     457              :                fout=fout, &
     458              :                openpmd_basename=openpmd_basename, &
     459              :                openpmd_unit_dimension=openpmd_unit_dimension, &
     460              :                openpmd_unit_si=openpmd_unit_si, &
     461            0 :                sim_time=sim_time)
     462              :       END IF
     463          552 :    END FUNCTION cp_forward_print_key_unit_nr
     464              : 
     465              : ! **************************************************************************************************
     466              : !> \brief Forwards to either `cp_pw_to_cube` or `cp_pw_to_openpmd`,
     467              : !>        depending on the configuration of the current base section.
     468              : !>        Writes data to either a Cube or an openPMD file.
     469              : !> \param self ...
     470              : !> \param pw ...
     471              : !> \param unit_nr ...
     472              : !> \param title ...
     473              : !> \param particles ...
     474              : !> \param zeff ...
     475              : !> \param stride ...
     476              : !> \param max_file_size_mb ...
     477              : !> \param zero_tails ...
     478              : !> \param silent ...
     479              : !> \param mpi_io ...
     480              : ! **************************************************************************************************
     481          552 :    SUBROUTINE cp_forward_write_pw( &
     482              :       self, &
     483              :       pw, &
     484              :       unit_nr, &
     485              :       title, &
     486              :       particles, &
     487          552 :       zeff, &
     488              :       stride, &
     489              :       max_file_size_mb, &
     490              :       zero_tails, &
     491              :       silent, &
     492              :       mpi_io &
     493              :       )
     494              :       CLASS(cp_section_key), INTENT(IN)                  :: self
     495              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: pw
     496              :       INTEGER, INTENT(IN)                                :: unit_nr
     497              :       CHARACTER(*), INTENT(IN), OPTIONAL                 :: title
     498              :       TYPE(particle_list_type), POINTER                  :: particles
     499              :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: stride
     500              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: max_file_size_mb
     501              :       LOGICAL, INTENT(IN), OPTIONAL                      :: zero_tails, silent, mpi_io
     502              :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: zeff
     503              : 
     504          552 :       IF (self%grid_output == grid_output_cubes) THEN
     505          552 :          CALL cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
     506              :       ELSE
     507            0 :          CALL cp_pw_to_openpmd(pw, unit_nr, title, particles, zeff, stride, zero_tails, silent, mpi_io)
     508              :       END IF
     509          552 :    END SUBROUTINE cp_forward_write_pw
     510              : 
     511              : ! **************************************************************************************************
     512              : !> \brief Forwards to either `cp_print_key_finished_output` or `cp_openpmd_print_key_finished_output`,
     513              : !>        depending on the configuration of the current base section.
     514              : !>        Closes either a Cube file or a reference to a section within an openPMD file.
     515              : !> \param self ...
     516              : !> \param unit_nr ...
     517              : !> \param logger ...
     518              : !> \param basis_section ...
     519              : !> \param print_key_path ...
     520              : !> \param local ...
     521              : !> \param ignore_should_output ...
     522              : !> \param on_file ...
     523              : !> \param mpi_io ...
     524              : ! **************************************************************************************************
     525          552 :    SUBROUTINE cp_forward_print_key_finished_output(self, unit_nr, logger, basis_section, &
     526              :                                                    print_key_path, local, ignore_should_output, on_file, &
     527              :                                                    mpi_io)
     528              :       CLASS(cp_section_key), INTENT(IN)                  :: self
     529              :       INTEGER, INTENT(INOUT)                             :: unit_nr
     530              :       TYPE(cp_logger_type), POINTER                      :: logger
     531              :       TYPE(section_vals_type), INTENT(IN)                :: basis_section
     532              :       CHARACTER(len=*), INTENT(IN), OPTIONAL             :: print_key_path
     533              :       LOGICAL, INTENT(IN), OPTIONAL                      :: local, ignore_should_output, on_file, &
     534              :                                                             mpi_io
     535              : 
     536          552 :       IF (self%grid_output == grid_output_cubes) THEN
     537          552 :      CALL cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
     538              :       ELSE
     539            0 :       CALL cp_openpmd_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, mpi_io)
     540              :       END IF
     541          552 :    END SUBROUTINE cp_forward_print_key_finished_output
     542              : 
     543              :    !
     544              : ! **************************************************************************************************
     545              : !> \brief Decides if a particular output routine will write to openPMD, to Cube or to none.
     546              : !>        Writing to both is not supported.
     547              : !>        The distinction between Cube and openPMD output works such that the output configuration
     548              : !>        sections exist as duplicates: E.g. for DFT%PRINT%MO_CUBES,
     549              : !>        there additionally exists DFT%PRINT%MO_OPENPMD.
     550              : !>        The internal base configuration for such sections is identical; additionally there
     551              : !>        exist format-specific options such as APPEND for Cube or OPENPMD_CFG_FILE for openPMD.
     552              : !>        The routines in this file alternate between using relative section paths without the
     553              : !>        %DFT prefix (e.g. PRINT%MO_CUBES) or absolute section paths with the %DF% prefix
     554              : !>        (e.g. DFT%PRINT%MO_CUBES). Call this routine with the relative paths.
     555              : !> \param input ...
     556              : !> \param str_cubes ...
     557              : !> \param str_openpmd ...
     558              : !> \param logger ...
     559              : !> \return ...
     560              : ! **************************************************************************************************
     561        37595 :    FUNCTION cube_or_openpmd(input, str_cubes, str_openpmd, logger) RESULT(res)
     562              :       TYPE(section_vals_type), POINTER                   :: input
     563              :       CHARACTER(len=*), INTENT(IN)                       :: str_cubes, str_openpmd
     564              :       TYPE(cp_logger_type), POINTER                      :: logger
     565              :       TYPE(cp_section_key)                               :: res
     566              : 
     567              :       LOGICAL                                            :: do_cubes, do_openpmd
     568              : 
     569              :       do_cubes = BTEST(cp_print_key_should_output( &
     570              :                        logger%iter_info, input, &
     571        37595 :                        "DFT%"//TRIM(ADJUSTL(str_cubes))), cp_p_file)
     572              :       do_openpmd = BTEST(cp_print_key_should_output( &
     573              :                          logger%iter_info, input, &
     574        37595 :                          "DFT%"//TRIM(ADJUSTL(str_openpmd))), cp_p_file)
     575              :       ! Having Cube and openPMD output both active should be theoretically possible.
     576              :       ! It would require some extra handling for the unit_nr return values.
     577              :       ! (e.g. returning the Cube unit_nr and internally storing the associated openPMD unit_nr).
     578        37595 :       CPASSERT(.NOT. (do_cubes .AND. do_openpmd))
     579        37595 :       res%do_output = do_cubes .OR. do_openpmd
     580        37595 :       IF (do_openpmd) THEN
     581            0 :          res%grid_output = grid_output_openpmd
     582            0 :          res%relative_section_key = TRIM(ADJUSTL(str_openpmd))
     583            0 :          res%format_name = "openPMD"
     584              :       ELSE
     585        37595 :          res%grid_output = grid_output_cubes
     586        37595 :          res%relative_section_key = TRIM(ADJUSTL(str_cubes))
     587        37595 :          res%format_name = "Cube"
     588              :       END IF
     589        37595 :       res%absolute_section_key = "DFT%"//TRIM(ADJUSTL(res%relative_section_key))
     590        37595 :    END FUNCTION cube_or_openpmd
     591              : 
     592              : ! **************************************************************************************************
     593              : !> \brief This section key is named WRITE_CUBE for Cube which does not make much sense
     594              : !>        for openPMD, so this key name has to be distinguished.
     595              : !> \param grid_output ...
     596              : !> \return ...
     597              : ! **************************************************************************************************
     598          296 :    FUNCTION section_key_do_write(grid_output) RESULT(res)
     599              :       INTEGER, INTENT(IN)                                :: grid_output
     600              :       CHARACTER(len=32)                                  :: res
     601              : 
     602          296 :       IF (grid_output == grid_output_cubes) THEN
     603          296 :          res = "%WRITE_CUBE"
     604            0 :       ELSE IF (grid_output == grid_output_openpmd) THEN
     605            0 :          res = "%WRITE_OPENPMD"
     606              :       END IF
     607          296 :    END FUNCTION section_key_do_write
     608              : 
     609              : ! **************************************************************************************************
     610              : !> \brief Prints the output message for density file writing
     611              : !> \param output_unit Unit number for output
     612              : !> \param prefix The message prefix (e.g., "The total electron density")
     613              : !> \param e_density_section Section key containing grid_output and format_name
     614              : !> \param filename The actual filename or pattern used
     615              : ! **************************************************************************************************
     616          105 :    SUBROUTINE print_density_output_message(output_unit, prefix, e_density_section, filename)
     617              :       INTEGER, INTENT(IN)                                :: output_unit
     618              :       CHARACTER(len=*), INTENT(IN)                       :: prefix
     619              :       TYPE(cp_section_key), INTENT(IN)                   :: e_density_section
     620              :       CHARACTER(len=*), INTENT(IN)                       :: filename
     621              : 
     622          105 :       IF (e_density_section%grid_output == grid_output_openpmd) THEN
     623              :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
     624              :             TRIM(prefix)//" is written in " &
     625              :             //e_density_section%format_name &
     626            0 :             //" file format to the file / file pattern:", &
     627            0 :             TRIM(filename)
     628              :       ELSE
     629              :          WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") &
     630              :             TRIM(prefix)//" is written in " &
     631              :             //e_density_section%format_name &
     632          105 :             //" file format to the file:", &
     633          210 :             TRIM(filename)
     634              :       END IF
     635          105 :    END SUBROUTINE print_density_output_message
     636              : 
     637              : ! **************************************************************************************************
     638              : !> \brief collects possible post - scf calculations and prints info / computes properties.
     639              : !> \param qs_env the qs_env in which the qs_env lives
     640              : !> \param wf_type ...
     641              : !> \param do_mp2 ...
     642              : !> \par History
     643              : !>      02.2003 created [fawzi]
     644              : !>      10.2004 moved here from qs_scf [Joost VandeVondele]
     645              : !>              started splitting out different subroutines
     646              : !>      10.2015 added header for wave-function correlated methods [Vladimir Rybkin]
     647              : !> \author fawzi
     648              : !> \note
     649              : !>      this function changes mo_eigenvectors and mo_eigenvalues, depending on the print keys.
     650              : !>      In particular, MO_CUBES causes the MOs to be rotated to make them eigenstates of the KS
     651              : !>      matrix, and mo_eigenvalues is updated accordingly. This can, for unconverged wavefunctions,
     652              : !>      change afterwards slightly the forces (hence small numerical differences between MD
     653              : !>      with and without the debug print level). Ideally this should not happen...
     654              : ! **************************************************************************************************
     655        12077 :    SUBROUTINE scf_post_calculation_gpw(qs_env, wf_type, do_mp2)
     656              : 
     657              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     658              :       CHARACTER(6), OPTIONAL                             :: wf_type
     659              :       LOGICAL, OPTIONAL                                  :: do_mp2
     660              : 
     661              :       CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_gpw', &
     662              :          warning_cube_kpoint = "Print MO cubes not implemented for k-point calculations", &
     663              :          warning_openpmd_kpoint = "Writing to openPMD not implemented for k-point calculations"
     664              : 
     665              :       INTEGER                                            :: handle, homo, ispin, min_lumos, n_rep, &
     666              :                                                             nchk_nmoloc, nhomo, nlumo, nlumo_stm, &
     667              :                                                             nlumos, nmo, nspins, output_unit, &
     668              :                                                             unit_nr
     669        12077 :       INTEGER, DIMENSION(:, :, :), POINTER               :: marked_states
     670              :       LOGICAL :: check_write, compute_lumos, do_homo, do_kpoints, do_mixed, do_stm, &
     671              :          do_wannier_cubes, has_homo, has_lumo, loc_explicit, loc_print_explicit, my_do_mp2, &
     672              :          my_localized_wfn, p_loc, p_loc_homo, p_loc_lumo, p_loc_mixed
     673              :       REAL(dp)                                           :: e_kin
     674              :       REAL(KIND=dp)                                      :: gap, homo_lumo(2, 2), total_zeff_corr
     675        12077 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues
     676              :       TYPE(admm_type), POINTER                           :: admm_env
     677        12077 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     678        12077 :       TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER        :: mixed_evals, occupied_evals, &
     679        12077 :                                                             unoccupied_evals, unoccupied_evals_stm
     680        12077 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: mixed_orbs, occupied_orbs
     681              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
     682        12077 :          TARGET                                          :: homo_localized, lumo_localized, &
     683        12077 :                                                             mixed_localized
     684        12077 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: lumo_ptr, mo_loc_history, &
     685        12077 :                                                             unoccupied_orbs, unoccupied_orbs_stm
     686              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     687              :       TYPE(cp_logger_type), POINTER                      :: logger
     688              :       TYPE(cp_section_key)                               :: mo_section
     689        12077 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_rmpv, matrix_p_mp2, matrix_s, &
     690        12077 :                                                             mo_derivs
     691        12077 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: kinetic_m, rho_ao
     692              :       TYPE(dft_control_type), POINTER                    :: dft_control
     693        12077 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     694        12077 :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     695              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     696              :       TYPE(particle_list_type), POINTER                  :: particles
     697        12077 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     698              :       TYPE(pw_c1d_gs_type)                               :: wf_g
     699              :       TYPE(pw_env_type), POINTER                         :: pw_env
     700        12077 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     701              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     702              :       TYPE(pw_r3d_rs_type)                               :: wf_r
     703        12077 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     704              :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env_homo, qs_loc_env_lumo, &
     705              :                                                             qs_loc_env_mixed
     706              :       TYPE(qs_rho_type), POINTER                         :: rho
     707              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     708              :       TYPE(qs_subsys_type), POINTER                      :: subsys
     709              :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     710              :       TYPE(scf_control_type), POINTER                    :: scf_control
     711              :       TYPE(section_vals_type), POINTER                   :: dft_section, input, loc_print_section, &
     712              :                                                             localize_section, print_key, &
     713              :                                                             stm_section
     714              : 
     715        12077 :       CALL timeset(routineN, handle)
     716              : 
     717        12077 :       logger => cp_get_default_logger()
     718        12077 :       output_unit = cp_logger_get_default_io_unit(logger)
     719              : 
     720              :       ! Print out the type of wavefunction to distinguish between SCF and post-SCF
     721        12077 :       my_do_mp2 = .FALSE.
     722        12077 :       IF (PRESENT(do_mp2)) my_do_mp2 = do_mp2
     723        12077 :       IF (PRESENT(wf_type)) THEN
     724          328 :          IF (output_unit > 0) THEN
     725          164 :             WRITE (UNIT=output_unit, FMT='(/,(T1,A))') REPEAT("-", 40)
     726          164 :             WRITE (UNIT=output_unit, FMT='(/,(T3,A,T19,A,T25,A))') "Properties from ", wf_type, " density"
     727          164 :             WRITE (UNIT=output_unit, FMT='(/,(T1,A))') REPEAT("-", 40)
     728              :          END IF
     729              :       END IF
     730              : 
     731              :       ! Writes the data that is already available in qs_env
     732        12077 :       CALL get_qs_env(qs_env, scf_env=scf_env)
     733              : 
     734        12077 :       my_localized_wfn = .FALSE.
     735        12077 :       NULLIFY (admm_env, dft_control, pw_env, auxbas_pw_pool, pw_pools, mos, rho, &
     736        12077 :                mo_coeff, ks_rmpv, matrix_s, qs_loc_env_homo, qs_loc_env_lumo, scf_control, &
     737        12077 :                unoccupied_orbs, mo_eigenvalues, unoccupied_evals, &
     738        12077 :                unoccupied_evals_stm, molecule_set, mo_derivs, &
     739        12077 :                subsys, particles, input, print_key, kinetic_m, marked_states, &
     740        12077 :                mixed_evals, qs_loc_env_mixed)
     741        12077 :       NULLIFY (lumo_ptr, rho_ao)
     742              : 
     743        12077 :       has_homo = .FALSE.
     744        12077 :       has_lumo = .FALSE.
     745        12077 :       p_loc = .FALSE.
     746        12077 :       p_loc_homo = .FALSE.
     747        12077 :       p_loc_lumo = .FALSE.
     748        12077 :       p_loc_mixed = .FALSE.
     749              : 
     750        12077 :       CPASSERT(ASSOCIATED(scf_env))
     751        12077 :       CPASSERT(ASSOCIATED(qs_env))
     752              :       ! Here we start with data that needs a postprocessing...
     753              :       CALL get_qs_env(qs_env, &
     754              :                       dft_control=dft_control, &
     755              :                       molecule_set=molecule_set, &
     756              :                       scf_control=scf_control, &
     757              :                       do_kpoints=do_kpoints, &
     758              :                       input=input, &
     759              :                       subsys=subsys, &
     760              :                       rho=rho, &
     761              :                       pw_env=pw_env, &
     762              :                       particle_set=particle_set, &
     763              :                       atomic_kind_set=atomic_kind_set, &
     764        12077 :                       qs_kind_set=qs_kind_set)
     765        12077 :       rtp_control => dft_control%rtp_control
     766        12077 :       CALL qs_subsys_get(subsys, particles=particles)
     767              : 
     768        12077 :       CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
     769              : 
     770        12077 :       IF (my_do_mp2) THEN
     771              :          ! Get the HF+MP2 density
     772          322 :          CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2)
     773          742 :          DO ispin = 1, dft_control%nspins
     774          742 :             CALL dbcsr_add(rho_ao(ispin, 1)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
     775              :          END DO
     776          322 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
     777          322 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
     778              :          ! In MP2 case update the Hartree potential
     779          322 :          CALL update_hartree_with_mp2(rho, qs_env)
     780              :       END IF
     781              : 
     782        12077 :       CALL write_available_results(qs_env, scf_env)
     783              : 
     784              :       !    **** the kinetic energy
     785        12077 :       IF (cp_print_key_should_output(logger%iter_info, input, &
     786              :                                      "DFT%PRINT%KINETIC_ENERGY") /= 0) THEN
     787           80 :          CALL get_qs_env(qs_env, kinetic_kp=kinetic_m)
     788           80 :          CPASSERT(ASSOCIATED(kinetic_m))
     789           80 :          CPASSERT(ASSOCIATED(kinetic_m(1, 1)%matrix))
     790           80 :          CALL calculate_ptrace(kinetic_m, rho_ao, e_kin, dft_control%nspins)
     791              :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%KINETIC_ENERGY", &
     792           80 :                                         extension=".Log")
     793           80 :          IF (unit_nr > 0) THEN
     794           40 :             WRITE (unit_nr, '(T3,A,T55,F25.14)') "Electronic kinetic energy:", e_kin
     795              :          END IF
     796              :          CALL cp_print_key_finished_output(unit_nr, logger, input, &
     797           80 :                                            "DFT%PRINT%KINETIC_ENERGY")
     798              :       END IF
     799              : 
     800              :       ! Atomic Charges that require further computation
     801        12077 :       CALL qs_scf_post_charges(input, logger, qs_env)
     802              : 
     803              :       ! Moments of charge distribution
     804        12077 :       CALL qs_scf_post_moments(input, logger, qs_env, output_unit)
     805              : 
     806              :       ! Determine if we need to computer properties using the localized centers
     807        12077 :       dft_section => section_vals_get_subs_vals(input, "DFT")
     808        12077 :       localize_section => section_vals_get_subs_vals(dft_section, "LOCALIZE")
     809        12077 :       loc_print_section => section_vals_get_subs_vals(localize_section, "PRINT")
     810        12077 :       CALL section_vals_get(localize_section, explicit=loc_explicit)
     811        12077 :       CALL section_vals_get(loc_print_section, explicit=loc_print_explicit)
     812              : 
     813              :       ! Print_keys controlled by localization
     814        12077 :       IF (loc_print_explicit) THEN
     815          100 :          print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_DIPOLES")
     816          100 :          p_loc = BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
     817          100 :          print_key => section_vals_get_subs_vals(loc_print_section, "TOTAL_DIPOLE")
     818          100 :          p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
     819          100 :          print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_CENTERS")
     820          100 :          p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
     821          100 :          print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_SPREADS")
     822          100 :          p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
     823          100 :          print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_CUBES")
     824          100 :          p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
     825          100 :          print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_STATES")
     826          100 :          p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
     827          100 :          print_key => section_vals_get_subs_vals(loc_print_section, "MOLECULAR_MOMENTS")
     828          100 :          p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
     829          100 :          print_key => section_vals_get_subs_vals(loc_print_section, "LOCALIZED_MOMENTS")
     830          100 :          p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
     831          100 :          print_key => section_vals_get_subs_vals(loc_print_section, "WANNIER_STATES")
     832          100 :          p_loc = p_loc .OR. BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
     833              :       ELSE
     834              :          p_loc = .FALSE.
     835              :       END IF
     836        12077 :       IF (loc_explicit) THEN
     837              :          p_loc_homo = (section_get_ival(localize_section, "STATES") == do_loc_homo .OR. &
     838          100 :                        section_get_ival(localize_section, "STATES") == do_loc_both) .AND. p_loc
     839              :          p_loc_lumo = (section_get_ival(localize_section, "STATES") == do_loc_lumo .OR. &
     840          100 :                        section_get_ival(localize_section, "STATES") == do_loc_both) .AND. p_loc
     841          100 :          p_loc_mixed = (section_get_ival(localize_section, "STATES") == do_loc_mixed) .AND. p_loc
     842          100 :          CALL section_vals_val_get(localize_section, "LIST_UNOCCUPIED", n_rep_val=n_rep)
     843              :       ELSE
     844        11977 :          p_loc_homo = .FALSE.
     845        11977 :          p_loc_lumo = .FALSE.
     846        11977 :          p_loc_mixed = .FALSE.
     847        11977 :          n_rep = 0
     848              :       END IF
     849              : 
     850        12077 :       IF (n_rep == 0 .AND. p_loc_lumo) THEN
     851              :          CALL cp_abort(__LOCATION__, "No LIST_UNOCCUPIED was specified, "// &
     852            0 :                        "therefore localization of unoccupied states will be skipped!")
     853            0 :          p_loc_lumo = .FALSE.
     854              :       END IF
     855              : 
     856              :       ! Control for STM
     857        12077 :       stm_section => section_vals_get_subs_vals(input, "DFT%PRINT%STM")
     858        12077 :       CALL section_vals_get(stm_section, explicit=do_stm)
     859        12077 :       nlumo_stm = 0
     860        12077 :       IF (do_stm) nlumo_stm = section_get_ival(stm_section, "NLUMO")
     861              : 
     862              :       ! check for CUBES or openPMD (MOs and WANNIERS)
     863        12077 :       mo_section = cube_or_openpmd(input, str_mo_cubes, str_mo_openpmd, logger)
     864              : 
     865        12077 :       IF (loc_print_explicit) THEN
     866              :          do_wannier_cubes = BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
     867          100 :                                                              "WANNIER_CUBES"), cp_p_file)
     868              :       ELSE
     869              :          do_wannier_cubes = .FALSE.
     870              :       END IF
     871        12077 :       nlumo = section_get_ival(dft_section, mo_section%concat_to_relative("%NLUMO"))
     872        12077 :       nhomo = section_get_ival(dft_section, mo_section%concat_to_relative("%NHOMO"))
     873              : 
     874              :       ! Setup the grids needed to compute a wavefunction given a vector..
     875        12077 :       IF (((mo_section%do_output .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc) THEN
     876              :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     877          216 :                          pw_pools=pw_pools)
     878          216 :          CALL auxbas_pw_pool%create_pw(wf_r)
     879          216 :          CALL auxbas_pw_pool%create_pw(wf_g)
     880              :       END IF
     881              : 
     882        12077 :       IF (dft_control%restricted) THEN
     883              :          !For ROKS useful only first term
     884           80 :          nspins = 1
     885              :       ELSE
     886        11997 :          nspins = dft_control%nspins
     887              :       END IF
     888              :       !Some info about ROKS
     889        12077 :       IF (dft_control%restricted .AND. (mo_section%do_output .OR. p_loc_homo)) THEN
     890            0 :          CALL cp_abort(__LOCATION__, "Unclear how we define MOs / localization in the restricted case ... ")
     891              :          ! It is possible to obtain Wannier centers for ROKS without rotations for SINGLE OCCUPIED ORBITALS
     892              :       END IF
     893              :       ! Makes the MOs eigenstates, computes eigenvalues, write cubes
     894        12077 :       IF (do_kpoints) THEN
     895          668 :          CPWARN_IF(mo_section%do_cubes(), warning_cube_kpoint)
     896          668 :          CPWARN_IF(mo_section%do_openpmd(), warning_openpmd_kpoint)
     897              :       ELSE
     898              :          CALL get_qs_env(qs_env, &
     899              :                          mos=mos, &
     900        11409 :                          matrix_ks=ks_rmpv)
     901        11409 :          IF ((mo_section%do_output .AND. nhomo /= 0) .OR. do_stm) THEN
     902          136 :             CALL get_qs_env(qs_env, mo_derivs=mo_derivs)
     903          136 :             IF (dft_control%do_admm) THEN
     904            0 :                CALL get_qs_env(qs_env, admm_env=admm_env)
     905            0 :                CALL make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env=admm_env)
     906              :             ELSE
     907          136 :                IF (dft_control%hairy_probes) THEN
     908            0 :                   scf_control%smear%do_smear = .FALSE.
     909              :                   CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs, &
     910              :                                    hairy_probes=dft_control%hairy_probes, &
     911            0 :                                    probe=dft_control%probe)
     912              :                ELSE
     913          136 :                   CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs)
     914              :                END IF
     915              :             END IF
     916          292 :             DO ispin = 1, dft_control%nspins
     917          156 :                CALL get_mo_set(mo_set=mos(ispin), eigenvalues=mo_eigenvalues, homo=homo)
     918          292 :                homo_lumo(ispin, 1) = mo_eigenvalues(homo)
     919              :             END DO
     920              :             has_homo = .TRUE.
     921              :          END IF
     922        11409 :          IF (mo_section%do_output .AND. nhomo /= 0) THEN
     923          278 :             DO ispin = 1, nspins
     924              :                ! Prints the cube files of OCCUPIED ORBITALS
     925              :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
     926          148 :                                eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
     927              :                CALL qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env, &
     928          278 :                                           mo_coeff, wf_g, wf_r, particles, homo, ispin, mo_section)
     929              :             END DO
     930              :          END IF
     931              :       END IF
     932              : 
     933              :       ! Initialize the localization environment, needed e.g. for wannier functions and molecular states
     934              :       ! Gets localization info for the occupied orbs
     935              :       !  - Possibly gets wannier functions
     936              :       !  - Possibly gets molecular states
     937        12077 :       IF (p_loc_homo) THEN
     938           94 :          IF (do_kpoints) THEN
     939            0 :             CPWARN("Localization not implemented for k-point calculations!")
     940              :          ELSE IF (dft_control%restricted &
     941              :                   .AND. (section_get_ival(localize_section, "METHOD") /= do_loc_none) &
     942           94 :                   .AND. (section_get_ival(localize_section, "METHOD") /= do_loc_jacobi)) THEN
     943            0 :             CPABORT("ROKS works only with LOCALIZE METHOD NONE or JACOBI")
     944              :          ELSE
     945          392 :             ALLOCATE (occupied_orbs(dft_control%nspins))
     946          392 :             ALLOCATE (occupied_evals(dft_control%nspins))
     947          392 :             ALLOCATE (homo_localized(dft_control%nspins))
     948          204 :             DO ispin = 1, dft_control%nspins
     949              :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
     950          110 :                                eigenvalues=mo_eigenvalues)
     951          110 :                occupied_orbs(ispin) = mo_coeff
     952          110 :                occupied_evals(ispin)%array => mo_eigenvalues
     953          110 :                CALL cp_fm_create(homo_localized(ispin), occupied_orbs(ispin)%matrix_struct)
     954          204 :                CALL cp_fm_to_fm(occupied_orbs(ispin), homo_localized(ispin))
     955              :             END DO
     956              : 
     957           94 :             CALL get_qs_env(qs_env, mo_loc_history=mo_loc_history)
     958           94 :             do_homo = .TRUE.
     959              : 
     960          752 :             ALLOCATE (qs_loc_env_homo)
     961           94 :             CALL qs_loc_env_create(qs_loc_env_homo)
     962           94 :             CALL qs_loc_control_init(qs_loc_env_homo, localize_section, do_homo=do_homo)
     963              :             CALL qs_loc_init(qs_env, qs_loc_env_homo, localize_section, homo_localized, do_homo, &
     964           94 :                              mo_section%do_output, mo_loc_history=mo_loc_history)
     965              :             CALL get_localization_info(qs_env, qs_loc_env_homo, localize_section, homo_localized, &
     966           94 :                                        wf_r, wf_g, particles, occupied_orbs, occupied_evals, marked_states)
     967              : 
     968              :             !retain the homo_localized for future use
     969           94 :             IF (qs_loc_env_homo%localized_wfn_control%use_history) THEN
     970           10 :                CALL retain_history(mo_loc_history, homo_localized)
     971           10 :                CALL set_qs_env(qs_env, mo_loc_history=mo_loc_history)
     972              :             END IF
     973              : 
     974              :             !write restart for localization of occupied orbitals
     975              :             CALL loc_write_restart(qs_loc_env_homo, loc_print_section, mos, &
     976           94 :                                    homo_localized, do_homo)
     977           94 :             CALL cp_fm_release(homo_localized)
     978           94 :             DEALLOCATE (occupied_orbs)
     979           94 :             DEALLOCATE (occupied_evals)
     980              :             ! Print Total Dipole if the localization has been performed
     981          188 :             IF (qs_loc_env_homo%do_localize) THEN
     982           78 :                CALL loc_dipole(input, dft_control, qs_loc_env_homo, logger, qs_env)
     983              :             END IF
     984              :          END IF
     985              :       END IF
     986              : 
     987              :       ! Gets the lumos, and eigenvalues for the lumos, and localize them if requested
     988        12077 :       IF (do_kpoints) THEN
     989          668 :          IF (mo_section%do_output .OR. p_loc_lumo) THEN
     990              :             ! nothing at the moment, not implemented
     991            2 :             CPWARN("Localization and MO related output not implemented for k-point calculations!")
     992              :          END IF
     993              :       ELSE
     994        11409 :          compute_lumos = mo_section%do_output .AND. nlumo /= 0
     995        11409 :          compute_lumos = compute_lumos .OR. p_loc_lumo
     996              : 
     997        25012 :          DO ispin = 1, dft_control%nspins
     998        13603 :             CALL get_mo_set(mo_set=mos(ispin), homo=homo, nmo=nmo)
     999        38567 :             compute_lumos = compute_lumos .AND. homo == nmo
    1000              :          END DO
    1001              : 
    1002        11409 :          IF (mo_section%do_output .AND. .NOT. compute_lumos) THEN
    1003              : 
    1004           98 :             nlumo = section_get_ival(dft_section, mo_section%concat_to_relative("%NLUMO"))
    1005          198 :             DO ispin = 1, dft_control%nspins
    1006              : 
    1007          100 :                CALL get_mo_set(mo_set=mos(ispin), homo=homo, nmo=nmo, eigenvalues=mo_eigenvalues)
    1008          198 :                IF (nlumo > nmo - homo) THEN
    1009              :                   ! this case not yet implemented
    1010              :                ELSE
    1011          100 :                   IF (nlumo == -1) THEN
    1012            0 :                      nlumo = nmo - homo
    1013              :                   END IF
    1014          100 :                   IF (output_unit > 0) WRITE (output_unit, *) " "
    1015          100 :                   IF (output_unit > 0) WRITE (output_unit, *) " Lowest eigenvalues of the unoccupied subspace spin ", ispin
    1016          100 :                   IF (output_unit > 0) WRITE (output_unit, *) "---------------------------------------------"
    1017          107 :                   IF (output_unit > 0) WRITE (output_unit, '(4(1X,1F16.8))') mo_eigenvalues(homo + 1:homo + nlumo)
    1018              : 
    1019              :                   ! Prints the cube files of UNOCCUPIED ORBITALS
    1020          100 :                   CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
    1021              :                   CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, &
    1022          100 :                                           mo_coeff, wf_g, wf_r, particles, nlumo, homo, ispin, lumo=homo + 1, mo_section=mo_section)
    1023              :                END IF
    1024              :             END DO
    1025              : 
    1026              :          END IF
    1027              : 
    1028        11377 :          IF (compute_lumos) THEN
    1029           32 :             check_write = .TRUE.
    1030           32 :             min_lumos = nlumo
    1031           32 :             IF (nlumo == 0) check_write = .FALSE.
    1032           32 :             IF (p_loc_lumo) THEN
    1033            6 :                do_homo = .FALSE.
    1034           48 :                ALLOCATE (qs_loc_env_lumo)
    1035            6 :                CALL qs_loc_env_create(qs_loc_env_lumo)
    1036            6 :                CALL qs_loc_control_init(qs_loc_env_lumo, localize_section, do_homo=do_homo)
    1037           98 :                min_lumos = MAX(MAXVAL(qs_loc_env_lumo%localized_wfn_control%loc_states(:, :)), nlumo)
    1038              :             END IF
    1039              : 
    1040          144 :             ALLOCATE (unoccupied_orbs(dft_control%nspins))
    1041          144 :             ALLOCATE (unoccupied_evals(dft_control%nspins))
    1042           32 :             CALL make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, min_lumos, nlumos)
    1043           32 :             lumo_ptr => unoccupied_orbs
    1044           80 :             DO ispin = 1, dft_control%nspins
    1045           48 :                has_lumo = .TRUE.
    1046           48 :                homo_lumo(ispin, 2) = unoccupied_evals(ispin)%array(1)
    1047           48 :                CALL get_mo_set(mo_set=mos(ispin), homo=homo)
    1048           80 :                IF (check_write) THEN
    1049           48 :                   IF (p_loc_lumo .AND. nlumo /= -1) nlumos = MIN(nlumo, nlumos)
    1050              :                   ! Prints the cube files of UNOCCUPIED ORBITALS
    1051              :                   CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, &
    1052           48 :                                           unoccupied_orbs(ispin), wf_g, wf_r, particles, nlumos, homo, ispin, mo_section=mo_section)
    1053              :                END IF
    1054              :             END DO
    1055              : 
    1056           64 :             IF (p_loc_lumo) THEN
    1057           30 :                ALLOCATE (lumo_localized(dft_control%nspins))
    1058           18 :                DO ispin = 1, dft_control%nspins
    1059           12 :                   CALL cp_fm_create(lumo_localized(ispin), unoccupied_orbs(ispin)%matrix_struct)
    1060           18 :                   CALL cp_fm_to_fm(unoccupied_orbs(ispin), lumo_localized(ispin))
    1061              :                END DO
    1062              :                CALL qs_loc_init(qs_env, qs_loc_env_lumo, localize_section, lumo_localized, do_homo, mo_section%do_output, &
    1063            6 :                                 evals=unoccupied_evals)
    1064              :                CALL qs_loc_env_init(qs_loc_env_lumo, qs_loc_env_lumo%localized_wfn_control, qs_env, &
    1065            6 :                                     loc_coeff=unoccupied_orbs)
    1066              :                CALL get_localization_info(qs_env, qs_loc_env_lumo, localize_section, &
    1067              :                                           lumo_localized, wf_r, wf_g, particles, &
    1068            6 :                                           unoccupied_orbs, unoccupied_evals, marked_states)
    1069              :                CALL loc_write_restart(qs_loc_env_lumo, loc_print_section, mos, homo_localized, do_homo, &
    1070            6 :                                       evals=unoccupied_evals)
    1071            6 :                lumo_ptr => lumo_localized
    1072              :             END IF
    1073              :          END IF
    1074              : 
    1075           32 :          IF (has_homo .AND. has_lumo) THEN
    1076           32 :             IF (output_unit > 0) WRITE (output_unit, *) " "
    1077           80 :             DO ispin = 1, dft_control%nspins
    1078           80 :                IF (.NOT. scf_control%smear%do_smear) THEN
    1079           48 :                   gap = homo_lumo(ispin, 2) - homo_lumo(ispin, 1)
    1080           48 :                   IF (output_unit > 0) WRITE (output_unit, '(T2,A,F12.6)') &
    1081           24 :                      "HOMO - LUMO gap [eV] :", gap*evolt
    1082              :                END IF
    1083              :             END DO
    1084              :          END IF
    1085              :       END IF
    1086              : 
    1087        12077 :       IF (p_loc_mixed) THEN
    1088            2 :          IF (do_kpoints) THEN
    1089            0 :             CPWARN("Localization not implemented for k-point calculations!")
    1090            2 :          ELSE IF (dft_control%restricted) THEN
    1091            0 :             IF (output_unit > 0) WRITE (output_unit, *) &
    1092            0 :                " Unclear how we define MOs / localization in the restricted case... skipping"
    1093              :          ELSE
    1094              : 
    1095            8 :             ALLOCATE (mixed_orbs(dft_control%nspins))
    1096            8 :             ALLOCATE (mixed_evals(dft_control%nspins))
    1097            8 :             ALLOCATE (mixed_localized(dft_control%nspins))
    1098            4 :             DO ispin = 1, dft_control%nspins
    1099              :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
    1100            2 :                                eigenvalues=mo_eigenvalues)
    1101            2 :                mixed_orbs(ispin) = mo_coeff
    1102            2 :                mixed_evals(ispin)%array => mo_eigenvalues
    1103            2 :                CALL cp_fm_create(mixed_localized(ispin), mixed_orbs(ispin)%matrix_struct)
    1104            4 :                CALL cp_fm_to_fm(mixed_orbs(ispin), mixed_localized(ispin))
    1105              :             END DO
    1106              : 
    1107            2 :             CALL get_qs_env(qs_env, mo_loc_history=mo_loc_history)
    1108            2 :             do_homo = .FALSE.
    1109            2 :             do_mixed = .TRUE.
    1110            2 :             total_zeff_corr = scf_env%sum_zeff_corr
    1111           16 :             ALLOCATE (qs_loc_env_mixed)
    1112            2 :             CALL qs_loc_env_create(qs_loc_env_mixed)
    1113            2 :             CALL qs_loc_control_init(qs_loc_env_mixed, localize_section, do_homo=do_homo, do_mixed=do_mixed)
    1114              :             CALL qs_loc_init(qs_env, qs_loc_env_mixed, localize_section, mixed_localized, do_homo, &
    1115              :                              mo_section%do_output, mo_loc_history=mo_loc_history, tot_zeff_corr=total_zeff_corr, &
    1116            2 :                              do_mixed=do_mixed)
    1117              : 
    1118            4 :             DO ispin = 1, dft_control%nspins
    1119            4 :                CALL cp_fm_get_info(mixed_localized(ispin), ncol_global=nchk_nmoloc)
    1120              :             END DO
    1121              : 
    1122              :             CALL get_localization_info(qs_env, qs_loc_env_mixed, localize_section, mixed_localized, &
    1123            2 :                                        wf_r, wf_g, particles, mixed_orbs, mixed_evals, marked_states)
    1124              : 
    1125              :             !retain the homo_localized for future use
    1126            2 :             IF (qs_loc_env_mixed%localized_wfn_control%use_history) THEN
    1127            0 :                CALL retain_history(mo_loc_history, mixed_localized)
    1128            0 :                CALL set_qs_env(qs_env, mo_loc_history=mo_loc_history)
    1129              :             END IF
    1130              : 
    1131              :             !write restart for localization of occupied orbitals
    1132              :             CALL loc_write_restart(qs_loc_env_mixed, loc_print_section, mos, &
    1133            2 :                                    mixed_localized, do_homo, do_mixed=do_mixed)
    1134            2 :             CALL cp_fm_release(mixed_localized)
    1135            2 :             DEALLOCATE (mixed_orbs)
    1136            4 :             DEALLOCATE (mixed_evals)
    1137              :          END IF
    1138              :       END IF
    1139              : 
    1140              :       ! Deallocate grids needed to compute wavefunctions
    1141        12077 :       IF (((mo_section%do_output .OR. do_wannier_cubes) .AND. (nlumo /= 0 .OR. nhomo /= 0)) .OR. p_loc) THEN
    1142          216 :          CALL auxbas_pw_pool%give_back_pw(wf_r)
    1143          216 :          CALL auxbas_pw_pool%give_back_pw(wf_g)
    1144              :       END IF
    1145              : 
    1146              :       ! Destroy the localization environment
    1147        12077 :       IF (.NOT. do_kpoints) THEN
    1148        11409 :          IF (p_loc_homo) THEN
    1149           94 :             CALL qs_loc_env_release(qs_loc_env_homo)
    1150           94 :             DEALLOCATE (qs_loc_env_homo)
    1151              :          END IF
    1152        11409 :          IF (p_loc_lumo) THEN
    1153            6 :             CALL qs_loc_env_release(qs_loc_env_lumo)
    1154            6 :             DEALLOCATE (qs_loc_env_lumo)
    1155              :          END IF
    1156        11409 :          IF (p_loc_mixed) THEN
    1157            2 :             CALL qs_loc_env_release(qs_loc_env_mixed)
    1158            2 :             DEALLOCATE (qs_loc_env_mixed)
    1159              :          END IF
    1160              :       END IF
    1161              : 
    1162              :       ! generate a mix of wfns, and write to a restart
    1163        12077 :       IF (do_kpoints) THEN
    1164              :          ! nothing at the moment, not implemented
    1165              :       ELSE
    1166        11409 :          CALL get_qs_env(qs_env, matrix_s=matrix_s, para_env=para_env)
    1167              :          CALL wfn_mix(mos, particle_set, dft_section, qs_kind_set, para_env, &
    1168              :                       output_unit, unoccupied_orbs=lumo_ptr, scf_env=scf_env, &
    1169        11409 :                       matrix_s=matrix_s, marked_states=marked_states)
    1170              : 
    1171        11409 :          IF (p_loc_lumo) CALL cp_fm_release(lumo_localized)
    1172              :       END IF
    1173        12077 :       IF (ASSOCIATED(marked_states)) THEN
    1174           16 :          DEALLOCATE (marked_states)
    1175              :       END IF
    1176              : 
    1177              :       ! This is just a deallocation for printing MO_CUBES or TDDFPT
    1178        12077 :       IF (.NOT. do_kpoints) THEN
    1179        11409 :          IF (compute_lumos) THEN
    1180           80 :             DO ispin = 1, dft_control%nspins
    1181           48 :                DEALLOCATE (unoccupied_evals(ispin)%array)
    1182           80 :                CALL cp_fm_release(unoccupied_orbs(ispin))
    1183              :             END DO
    1184           32 :             DEALLOCATE (unoccupied_evals)
    1185           32 :             DEALLOCATE (unoccupied_orbs)
    1186              :          END IF
    1187              :       END IF
    1188              : 
    1189              :       !stm images
    1190        12077 :       IF (do_stm) THEN
    1191            6 :          IF (do_kpoints) THEN
    1192            0 :             CPWARN("STM not implemented for k-point calculations!")
    1193              :          ELSE
    1194            6 :             NULLIFY (unoccupied_orbs_stm, unoccupied_evals_stm)
    1195            6 :             IF (nlumo_stm > 0) THEN
    1196            8 :                ALLOCATE (unoccupied_orbs_stm(dft_control%nspins))
    1197            8 :                ALLOCATE (unoccupied_evals_stm(dft_control%nspins))
    1198              :                CALL make_lumo_gpw(qs_env, scf_env, unoccupied_orbs_stm, unoccupied_evals_stm, &
    1199            2 :                                   nlumo_stm, nlumos)
    1200              :             END IF
    1201              : 
    1202              :             CALL th_stm_image(qs_env, stm_section, particles, unoccupied_orbs_stm, &
    1203            6 :                               unoccupied_evals_stm)
    1204              : 
    1205            6 :             IF (nlumo_stm > 0) THEN
    1206            4 :                DO ispin = 1, dft_control%nspins
    1207            4 :                   DEALLOCATE (unoccupied_evals_stm(ispin)%array)
    1208              :                END DO
    1209            2 :                DEALLOCATE (unoccupied_evals_stm)
    1210            2 :                CALL cp_fm_release(unoccupied_orbs_stm)
    1211              :             END IF
    1212              :          END IF
    1213              :       END IF
    1214              : 
    1215              :       ! Print coherent X-ray diffraction spectrum
    1216        12077 :       CALL qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit)
    1217              : 
    1218              :       ! Calculation of Electric Field Gradients
    1219        12077 :       CALL qs_scf_post_efg(input, logger, qs_env)
    1220              : 
    1221              :       ! Calculation of ET
    1222        12077 :       CALL qs_scf_post_et(input, qs_env, dft_control)
    1223              : 
    1224              :       ! Calculation of EPR Hyperfine Coupling Tensors
    1225        12077 :       CALL qs_scf_post_epr(input, logger, qs_env)
    1226              : 
    1227              :       ! Calculation of properties needed for BASIS_MOLOPT optimizations
    1228        12077 :       CALL qs_scf_post_molopt(input, logger, qs_env)
    1229              : 
    1230              :       ! Calculate ELF
    1231        12077 :       CALL qs_scf_post_elf(input, logger, qs_env)
    1232              : 
    1233              :       ! Use Wannier90 interface
    1234        12077 :       CALL wannier90_interface(input, logger, qs_env)
    1235              : 
    1236        12077 :       IF (my_do_mp2) THEN
    1237              :          ! Get everything back
    1238          742 :          DO ispin = 1, dft_control%nspins
    1239          742 :             CALL dbcsr_add(rho_ao(ispin, 1)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
    1240              :          END DO
    1241          322 :          CALL qs_rho_update_rho(rho, qs_env=qs_env)
    1242          322 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
    1243              :       END IF
    1244              : 
    1245        12077 :       CALL cp_openpmd_close_iterations()
    1246              : 
    1247        12077 :       CALL timestop(handle)
    1248              : 
    1249        24154 :    END SUBROUTINE scf_post_calculation_gpw
    1250              : 
    1251              : ! **************************************************************************************************
    1252              : !> \brief Gets the LUMOs and their eigenvalues for all spin channels.
    1253              : !> \param qs_env ...
    1254              : !> \param scf_env ...
    1255              : !> \param unoccupied_orbs ...
    1256              : !> \param unoccupied_evals ...
    1257              : !> \param nlumo ...
    1258              : !> \param nlumos ...
    1259              : ! **************************************************************************************************
    1260           36 :    SUBROUTINE make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, nlumo, nlumos)
    1261              : 
    1262              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1263              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
    1264              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT)      :: unoccupied_orbs
    1265              :       TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER        :: unoccupied_evals
    1266              :       INTEGER, INTENT(IN)                                :: nlumo
    1267              :       INTEGER, INTENT(OUT)                               :: nlumos
    1268              : 
    1269              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'make_lumo_gpw'
    1270              : 
    1271              :       INTEGER                                            :: handle, homo, ispin, n, nao, nmo, &
    1272              :                                                             output_unit
    1273              :       TYPE(admm_type), POINTER                           :: admm_env
    1274              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    1275              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
    1276              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    1277              :       TYPE(cp_logger_type), POINTER                      :: logger
    1278           36 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_rmpv, matrix_s
    1279              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1280           36 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1281              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1282              :       TYPE(preconditioner_type), POINTER                 :: local_preconditioner
    1283              :       TYPE(scf_control_type), POINTER                    :: scf_control
    1284              : 
    1285           36 :       CALL timeset(routineN, handle)
    1286              : 
    1287           36 :       NULLIFY (ks_rmpv, matrix_s, scf_control, dft_control, admm_env, para_env, blacs_env, mos)
    1288              :       CALL get_qs_env(qs_env, &
    1289              :                       matrix_ks=ks_rmpv, &
    1290              :                       matrix_s=matrix_s, &
    1291              :                       scf_control=scf_control, &
    1292              :                       dft_control=dft_control, &
    1293              :                       admm_env=admm_env, &
    1294              :                       para_env=para_env, &
    1295              :                       blacs_env=blacs_env, &
    1296           36 :                       mos=mos)
    1297              : 
    1298           36 :       logger => cp_get_default_logger()
    1299           36 :       output_unit = cp_logger_get_default_io_unit(logger)
    1300              : 
    1301           88 :       DO ispin = 1, dft_control%nspins
    1302           52 :          NULLIFY (unoccupied_evals(ispin)%array)
    1303           52 :          IF (output_unit > 0) WRITE (output_unit, *) " "
    1304           52 :          IF (output_unit > 0) WRITE (output_unit, *) &
    1305           26 :             " Using OT eigensolver for additional unoccupied orbitals spin ", ispin
    1306           52 :          IF (output_unit > 0) WRITE (output_unit, *) &
    1307           26 :             " Lowest Eigenvalues of the unoccupied subspace spin ", ispin
    1308           52 :          IF (output_unit > 0) WRITE (output_unit, FMT='(1X,A)') "-----------------------------------------------------"
    1309           52 :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=homo, nao=nao, nmo=nmo)
    1310           52 :          CALL cp_fm_get_info(mo_coeff, nrow_global=n)
    1311           52 :          nlumos = MAX(1, MIN(nlumo, nao - nmo))
    1312           52 :          IF (nlumo == -1) nlumos = nao - nmo
    1313          156 :          ALLOCATE (unoccupied_evals(ispin)%array(nlumos))
    1314              :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
    1315           52 :                                   nrow_global=n, ncol_global=nlumos)
    1316           52 :          CALL cp_fm_create(unoccupied_orbs(ispin), fm_struct_tmp, name="lumos")
    1317           52 :          CALL cp_fm_struct_release(fm_struct_tmp)
    1318           52 :          CALL cp_fm_init_random(unoccupied_orbs(ispin), nlumos)
    1319              : 
    1320              :          ! FULL_ALL has column-dependent occupied-state weights and cannot be reused here.
    1321           52 :          NULLIFY (local_preconditioner)
    1322           52 :          IF (ASSOCIATED(scf_env)) THEN
    1323           52 :             IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
    1324           28 :                local_preconditioner => scf_env%ot_preconditioner(1)%preconditioner
    1325           28 :                IF (local_preconditioner%in_use == ot_precond_full_all) THEN
    1326            4 :                   NULLIFY (local_preconditioner)
    1327              :                END IF
    1328              :             END IF
    1329              :          END IF
    1330              : 
    1331              :          ! If we do ADMM, we add have to modify the Kohn-Sham matrix
    1332           52 :          IF (dft_control%do_admm) THEN
    1333            0 :             CALL admm_correct_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
    1334              :          END IF
    1335              : 
    1336              :          CALL ot_eigensolver(matrix_h=ks_rmpv(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
    1337              :                              matrix_c_fm=unoccupied_orbs(ispin), &
    1338              :                              matrix_orthogonal_space_fm=mo_coeff, &
    1339              :                              eps_gradient=scf_control%eps_lumos, &
    1340              :                              preconditioner=local_preconditioner, &
    1341              :                              iter_max=scf_control%max_iter_lumos, &
    1342           52 :                              size_ortho_space=nmo)
    1343              : 
    1344              :          CALL calculate_subspace_eigenvalues(unoccupied_orbs(ispin), ks_rmpv(ispin)%matrix, &
    1345              :                                              unoccupied_evals(ispin)%array, scr=output_unit, &
    1346           52 :                                              ionode=output_unit > 0)
    1347              : 
    1348              :          ! If we do ADMM, we restore the original Kohn-Sham matrix
    1349          140 :          IF (dft_control%do_admm) THEN
    1350            0 :             CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
    1351              :          END IF
    1352              : 
    1353              :       END DO
    1354              : 
    1355           36 :       CALL timestop(handle)
    1356              : 
    1357           36 :    END SUBROUTINE make_lumo_gpw
    1358              : 
    1359              : ! **************************************************************************************************
    1360              : !> \brief Computes and Prints Atomic Charges with several methods
    1361              : !> \param input ...
    1362              : !> \param logger ...
    1363              : !> \param qs_env the qs_env in which the qs_env lives
    1364              : ! **************************************************************************************************
    1365        12077 :    SUBROUTINE qs_scf_post_charges(input, logger, qs_env)
    1366              :       TYPE(section_vals_type), POINTER                   :: input
    1367              :       TYPE(cp_logger_type), POINTER                      :: logger
    1368              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1369              : 
    1370              :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_charges'
    1371              : 
    1372              :       INTEGER                                            :: handle, print_level, unit_nr
    1373              :       LOGICAL                                            :: do_kpoints, print_it
    1374              :       TYPE(section_vals_type), POINTER                   :: density_fit_section, print_key
    1375              : 
    1376        12077 :       CALL timeset(routineN, handle)
    1377              : 
    1378        12077 :       CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
    1379              : 
    1380              :       ! Mulliken charges require no further computation and are printed from write_mo_free_results
    1381              : 
    1382              :       ! Compute the Lowdin charges
    1383        12077 :       print_key => section_vals_get_subs_vals(input, "DFT%PRINT%LOWDIN")
    1384        12077 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    1385              :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%LOWDIN", extension=".lowdin", &
    1386           92 :                                         log_filename=.FALSE.)
    1387           92 :          print_level = 1
    1388           92 :          CALL section_vals_val_get(print_key, "PRINT_GOP", l_val=print_it)
    1389           92 :          IF (print_it) print_level = 2
    1390           92 :          CALL section_vals_val_get(print_key, "PRINT_ALL", l_val=print_it)
    1391           92 :          IF (print_it) print_level = 3
    1392           92 :          CALL lowdin_population_analysis(qs_env, unit_nr, print_level)
    1393           92 :          CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%LOWDIN")
    1394              :       END IF
    1395              : 
    1396              :       ! Compute the RESP charges
    1397        12077 :       CALL resp_fit(qs_env)
    1398              : 
    1399              :       ! Compute the Density Derived Atomic Point charges with the Bloechl scheme
    1400        12077 :       print_key => section_vals_get_subs_vals(input, "PROPERTIES%FIT_CHARGE")
    1401        12077 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    1402              :          unit_nr = cp_print_key_unit_nr(logger, input, "PROPERTIES%FIT_CHARGE", extension=".Fitcharge", &
    1403          102 :                                         log_filename=.FALSE.)
    1404          102 :          density_fit_section => section_vals_get_subs_vals(input, "DFT%DENSITY_FITTING")
    1405          102 :          CALL get_ddapc(qs_env, .FALSE., density_fit_section, iwc=unit_nr)
    1406          102 :          CALL cp_print_key_finished_output(unit_nr, logger, input, "PROPERTIES%FIT_CHARGE")
    1407              :       END IF
    1408              : 
    1409        12077 :       CALL timestop(handle)
    1410              : 
    1411        12077 :    END SUBROUTINE qs_scf_post_charges
    1412              : 
    1413              : ! **************************************************************************************************
    1414              : !> \brief Computes and prints the Cube Files for MO
    1415              : !> \param input ...
    1416              : !> \param dft_section ...
    1417              : !> \param dft_control ...
    1418              : !> \param logger ...
    1419              : !> \param qs_env the qs_env in which the qs_env lives
    1420              : !> \param mo_coeff ...
    1421              : !> \param wf_g ...
    1422              : !> \param wf_r ...
    1423              : !> \param particles ...
    1424              : !> \param homo ...
    1425              : !> \param ispin ...
    1426              : !> \param mo_section ...
    1427              : ! **************************************************************************************************
    1428          148 :    SUBROUTINE qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env, &
    1429              :                                     mo_coeff, wf_g, wf_r, particles, homo, ispin, mo_section)
    1430              :       TYPE(section_vals_type), POINTER                   :: input, dft_section
    1431              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1432              :       TYPE(cp_logger_type), POINTER                      :: logger
    1433              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1434              :       TYPE(cp_fm_type), INTENT(IN)                       :: mo_coeff
    1435              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                :: wf_g
    1436              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: wf_r
    1437              :       TYPE(particle_list_type), POINTER                  :: particles
    1438              :       INTEGER, INTENT(IN)                                :: homo, ispin
    1439              :       TYPE(cp_section_key)                               :: mo_section
    1440              : 
    1441              :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_occ_cubes'
    1442              : 
    1443              :       CHARACTER(LEN=default_path_length)                 :: filename, my_pos_cube, title
    1444              :       INTEGER                                            :: handle, i, ir, ivector, n_rep, nhomo, &
    1445              :                                                             nlist, unit_nr
    1446          148 :       INTEGER, DIMENSION(:), POINTER                     :: list, list_index
    1447              :       LOGICAL                                            :: append_cube, mpi_io
    1448          148 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: zcharge
    1449          148 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1450              :       TYPE(cell_type), POINTER                           :: cell
    1451          148 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1452              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1453          148 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1454              : 
    1455          148 :       CALL timeset(routineN, handle)
    1456              : 
    1457              : #ifndef __OPENPMD
    1458              :       ! Error should usually be caught earlier as PRINT%MO_OPENPMD is not added to the input section
    1459              :       ! if openPMD is not activated
    1460          148 :       CPASSERT(mo_section%grid_output /= grid_output_openpmd)
    1461              : #endif
    1462              : 
    1463          148 :       NULLIFY (list_index)
    1464              : 
    1465              :       IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, mo_section%relative_section_key) &
    1466          148 :  , cp_p_file) .AND. section_get_lval(dft_section, mo_section%concat_to_relative(section_key_do_write(mo_section%grid_output)))) THEN
    1467          110 :          CALL get_effective_core_charges(qs_env, zcharge)
    1468          110 :          nhomo = section_get_ival(dft_section, mo_section%concat_to_relative("%NHOMO"))
    1469              :          ! For openPMD, refer to access modes instead of APPEND key
    1470          110 :          IF (mo_section%grid_output == grid_output_cubes) THEN
    1471          110 :             append_cube = section_get_lval(dft_section, mo_section%concat_to_relative("%APPEND"))
    1472              :          END IF
    1473          110 :          my_pos_cube = "REWIND"
    1474          110 :          IF (append_cube) THEN
    1475            0 :             my_pos_cube = "APPEND"
    1476              :          END IF
    1477          110 :          CALL section_vals_val_get(dft_section, mo_section%concat_to_relative("%HOMO_LIST"), n_rep_val=n_rep)
    1478          110 :          IF (n_rep > 0) THEN ! write the cubes of the list
    1479            0 :             nlist = 0
    1480            0 :             DO ir = 1, n_rep
    1481            0 :                NULLIFY (list)
    1482              :                CALL section_vals_val_get(dft_section, mo_section%concat_to_relative("%HOMO_LIST"), i_rep_val=ir, &
    1483            0 :                                          i_vals=list)
    1484            0 :                IF (ASSOCIATED(list)) THEN
    1485            0 :                   CALL reallocate(list_index, 1, nlist + SIZE(list))
    1486            0 :                   DO i = 1, SIZE(list)
    1487            0 :                      list_index(i + nlist) = list(i)
    1488              :                   END DO
    1489            0 :                   nlist = nlist + SIZE(list)
    1490              :                END IF
    1491              :             END DO
    1492              :          ELSE
    1493              : 
    1494          110 :             IF (nhomo == -1) nhomo = homo
    1495          110 :             nlist = homo - MAX(1, homo - nhomo + 1) + 1
    1496          330 :             ALLOCATE (list_index(nlist))
    1497          224 :             DO i = 1, nlist
    1498          224 :                list_index(i) = MAX(1, homo - nhomo + 1) + i - 1
    1499              :             END DO
    1500              :          END IF
    1501          224 :          DO i = 1, nlist
    1502          114 :             ivector = list_index(i)
    1503              :             CALL get_qs_env(qs_env=qs_env, &
    1504              :                             atomic_kind_set=atomic_kind_set, &
    1505              :                             qs_kind_set=qs_kind_set, &
    1506              :                             cell=cell, &
    1507              :                             particle_set=particle_set, &
    1508          114 :                             pw_env=pw_env)
    1509              :             CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
    1510          114 :                                         cell, dft_control, particle_set, pw_env)
    1511          114 :             WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", ivector, "_", ispin
    1512          114 :             mpi_io = .TRUE.
    1513              : 
    1514              :             unit_nr = mo_section%print_key_unit_nr( &
    1515              :                       logger, &
    1516              :                       input, &
    1517              :                       mo_section%absolute_section_key, &
    1518              :                       extension=".cube", &
    1519              :                       middle_name=TRIM(filename), &
    1520              :                       file_position=my_pos_cube, &
    1521              :                       log_filename=.FALSE., &
    1522              :                       mpi_io=mpi_io, &
    1523              :                       openpmd_basename="dft-mo", &
    1524              :                       openpmd_unit_dimension=openpmd_unit_dimension_wavefunction, &
    1525              :                       openpmd_unit_si=openpmd_unit_si_wavefunction, &
    1526          114 :                       sim_time=qs_env%sim_time)
    1527          114 :             WRITE (title, *) "WAVEFUNCTION ", ivector, " spin ", ispin, " i.e. HOMO - ", ivector - homo
    1528              :             CALL mo_section%write_pw(wf_r, unit_nr, title, particles=particles, zeff=zcharge, &
    1529              :                                      stride=section_get_ivals(dft_section, mo_section%concat_to_relative("%STRIDE")), &
    1530              :                                      max_file_size_mb=section_get_rval(dft_section, "PRINT%MO_CUBES%MAX_FILE_SIZE_MB"), &
    1531          114 :                                      mpi_io=mpi_io)
    1532          224 :             CALL mo_section%print_key_finished_output(unit_nr, logger, input, mo_section%absolute_section_key, mpi_io=mpi_io)
    1533              :          END DO
    1534          110 :          IF (ASSOCIATED(list_index)) DEALLOCATE (list_index)
    1535          258 :          DEALLOCATE (zcharge)
    1536              :       END IF
    1537              : 
    1538          148 :       CALL timestop(handle)
    1539              : 
    1540          296 :    END SUBROUTINE qs_scf_post_occ_cubes
    1541              : 
    1542              : ! **************************************************************************************************
    1543              : !> \brief Computes and prints the Cube Files for MO
    1544              : !> \param input ...
    1545              : !> \param dft_section ...
    1546              : !> \param dft_control ...
    1547              : !> \param logger ...
    1548              : !> \param qs_env the qs_env in which the qs_env lives
    1549              : !> \param unoccupied_orbs ...
    1550              : !> \param wf_g ...
    1551              : !> \param wf_r ...
    1552              : !> \param particles ...
    1553              : !> \param nlumos ...
    1554              : !> \param homo ...
    1555              : !> \param ispin ...
    1556              : !> \param lumo ...
    1557              : !> \param mo_section ...
    1558              : ! **************************************************************************************************
    1559          148 :    SUBROUTINE qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env, &
    1560              :                                       unoccupied_orbs, wf_g, wf_r, particles, nlumos, homo, ispin, lumo, mo_section)
    1561              : 
    1562              :       TYPE(section_vals_type), POINTER                   :: input, dft_section
    1563              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1564              :       TYPE(cp_logger_type), POINTER                      :: logger
    1565              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1566              :       TYPE(cp_fm_type), INTENT(IN)                       :: unoccupied_orbs
    1567              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                :: wf_g
    1568              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: wf_r
    1569              :       TYPE(particle_list_type), POINTER                  :: particles
    1570              :       INTEGER, INTENT(IN)                                :: nlumos, homo, ispin
    1571              :       INTEGER, INTENT(IN), OPTIONAL                      :: lumo
    1572              :       TYPE(cp_section_key)                               :: mo_section
    1573              : 
    1574              :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_unocc_cubes'
    1575              : 
    1576              :       CHARACTER(LEN=default_path_length)                 :: filename, my_pos_cube, title
    1577              :       INTEGER                                            :: handle, ifirst, index_mo, ivector, &
    1578              :                                                             unit_nr
    1579              :       LOGICAL                                            :: append_cube, mpi_io
    1580          148 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: zcharge
    1581          148 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1582              :       TYPE(cell_type), POINTER                           :: cell
    1583          148 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1584              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1585          148 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1586              : 
    1587          148 :       CALL timeset(routineN, handle)
    1588              : 
    1589              : #ifndef __OPENPMD
    1590              :       ! Error should usually be caught earlier as PRINT%MO_OPENPMD is not added to the input section
    1591              :       ! if openPMD is not activated
    1592          148 :       CPASSERT(mo_section%grid_output /= grid_output_openpmd)
    1593              : #endif
    1594              : 
    1595              :       IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, mo_section%relative_section_key), cp_p_file) &
    1596          148 :           .AND. section_get_lval(dft_section, mo_section%concat_to_relative(section_key_do_write(mo_section%grid_output)))) THEN
    1597          110 :          CALL get_effective_core_charges(qs_env, zcharge)
    1598          110 :          NULLIFY (qs_kind_set, particle_set, pw_env, cell)
    1599              :          ! For openPMD, refer to access modes instead of APPEND key
    1600          110 :          IF (mo_section%grid_output == grid_output_cubes) THEN
    1601          110 :             append_cube = section_get_lval(dft_section, mo_section%concat_to_relative("%APPEND"))
    1602              :          END IF
    1603          110 :          my_pos_cube = "REWIND"
    1604          110 :          IF (append_cube) THEN
    1605            0 :             my_pos_cube = "APPEND"
    1606              :          END IF
    1607          110 :          ifirst = 1
    1608          110 :          IF (PRESENT(lumo)) ifirst = lumo
    1609          254 :          DO ivector = ifirst, ifirst + nlumos - 1
    1610              :             CALL get_qs_env(qs_env=qs_env, &
    1611              :                             atomic_kind_set=atomic_kind_set, &
    1612              :                             qs_kind_set=qs_kind_set, &
    1613              :                             cell=cell, &
    1614              :                             particle_set=particle_set, &
    1615          144 :                             pw_env=pw_env)
    1616              :             CALL calculate_wavefunction(unoccupied_orbs, ivector, wf_r, wf_g, atomic_kind_set, &
    1617          144 :                                         qs_kind_set, cell, dft_control, particle_set, pw_env)
    1618              : 
    1619          144 :             IF (ifirst == 1) THEN
    1620          130 :                index_mo = homo + ivector
    1621              :             ELSE
    1622           14 :                index_mo = ivector
    1623              :             END IF
    1624          144 :             WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", index_mo, "_", ispin
    1625          144 :             mpi_io = .TRUE.
    1626              : 
    1627              :             unit_nr = mo_section%print_key_unit_nr( &
    1628              :                       logger, &
    1629              :                       input, &
    1630              :                       mo_section%absolute_section_key, &
    1631              :                       extension=".cube", &
    1632              :                       middle_name=TRIM(filename), &
    1633              :                       file_position=my_pos_cube, &
    1634              :                       log_filename=.FALSE., &
    1635              :                       mpi_io=mpi_io, &
    1636              :                       openpmd_basename="dft-mo", &
    1637              :                       openpmd_unit_dimension=openpmd_unit_dimension_wavefunction, &
    1638              :                       openpmd_unit_si=openpmd_unit_si_wavefunction, &
    1639          144 :                       sim_time=qs_env%sim_time)
    1640          144 :             WRITE (title, *) "WAVEFUNCTION ", index_mo, " spin ", ispin, " i.e. LUMO + ", ifirst + ivector - 2
    1641              :             CALL mo_section%write_pw(wf_r, unit_nr, title, particles=particles, zeff=zcharge, &
    1642              :                                      stride=section_get_ivals(dft_section, mo_section%concat_to_relative("%STRIDE")), &
    1643              :                                      max_file_size_mb=section_get_rval(dft_section, "PRINT%MO_CUBES%MAX_FILE_SIZE_MB"), &
    1644          144 :                                      mpi_io=mpi_io)
    1645          254 :             CALL mo_section%print_key_finished_output(unit_nr, logger, input, mo_section%absolute_section_key, mpi_io=mpi_io)
    1646              : 
    1647              :          END DO
    1648          258 :          DEALLOCATE (zcharge)
    1649              :       END IF
    1650              : 
    1651          148 :       CALL timestop(handle)
    1652              : 
    1653          296 :    END SUBROUTINE qs_scf_post_unocc_cubes
    1654              : 
    1655              : ! **************************************************************************************************
    1656              : !> \brief Computes and prints electric moments
    1657              : !> \param input ...
    1658              : !> \param logger ...
    1659              : !> \param qs_env the qs_env in which the qs_env lives
    1660              : !> \param output_unit ...
    1661              : ! **************************************************************************************************
    1662        13381 :    SUBROUTINE qs_scf_post_moments(input, logger, qs_env, output_unit)
    1663              :       TYPE(section_vals_type), POINTER                   :: input
    1664              :       TYPE(cp_logger_type), POINTER                      :: logger
    1665              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1666              :       INTEGER, INTENT(IN)                                :: output_unit
    1667              : 
    1668              :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_moments'
    1669              : 
    1670              :       CHARACTER(LEN=default_path_length)                 :: filename
    1671              :       INTEGER                                            :: handle, max_nmo, maxmom, moments_format, &
    1672              :                                                             moments_unit_nr, reference, unit_nr
    1673              :       LOGICAL                                            :: com_nl, do_kg, do_kpoints, magnetic, &
    1674              :                                                             new_file, periodic, second_ref_point, &
    1675              :                                                             vel_reprs
    1676        13381 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ref_point
    1677              :       TYPE(section_vals_type), POINTER                   :: print_key
    1678              : 
    1679        13381 :       CALL timeset(routineN, handle)
    1680              : 
    1681              :       print_key => section_vals_get_subs_vals(section_vals=input, &
    1682        13381 :                                               subsection_name="DFT%PRINT%MOMENTS")
    1683              : 
    1684        13381 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    1685              : 
    1686              :          maxmom = section_get_ival(section_vals=input, &
    1687         1570 :                                    keyword_name="DFT%PRINT%MOMENTS%MAX_MOMENT")
    1688              :          moments_format = section_get_ival(section_vals=input, &
    1689         1570 :                                            keyword_name="DFT%PRINT%MOMENTS%FORMAT")
    1690              :          periodic = section_get_lval(section_vals=input, &
    1691         1570 :                                      keyword_name="DFT%PRINT%MOMENTS%PERIODIC")
    1692              :          reference = section_get_ival(section_vals=input, &
    1693         1570 :                                       keyword_name="DFT%PRINT%MOMENTS%REFERENCE")
    1694              :          magnetic = section_get_lval(section_vals=input, &
    1695         1570 :                                      keyword_name="DFT%PRINT%MOMENTS%MAGNETIC")
    1696              :          vel_reprs = section_get_lval(section_vals=input, &
    1697         1570 :                                       keyword_name="DFT%PRINT%MOMENTS%VEL_REPRS")
    1698              :          com_nl = section_get_lval(section_vals=input, &
    1699         1570 :                                    keyword_name="DFT%PRINT%MOMENTS%COM_NL")
    1700              :          second_ref_point = section_get_lval(section_vals=input, &
    1701         1570 :                                              keyword_name="DFT%PRINT%MOMENTS%SECOND_REFERENCE_POINT")
    1702              :          do_kg = section_get_lval(section_vals=input, &
    1703         1570 :                                   keyword_name="DFT%PRINT%MOMENTS%KG")
    1704              :          max_nmo = section_get_ival(section_vals=input, &
    1705         1570 :                                     keyword_name="DFT%PRINT%MOMENTS%MAX_NMO")
    1706              : 
    1707         1570 :          NULLIFY (ref_point)
    1708         1570 :          CALL section_vals_val_get(input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
    1709              :          unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=input, &
    1710              :                                         print_key_path="DFT%PRINT%MOMENTS", extension=".dat", &
    1711              :                                         middle_name="moments", log_filename=.FALSE., &
    1712         1570 :                                         is_new_file=new_file)
    1713              : 
    1714         1570 :          IF (output_unit > 0) THEN
    1715          795 :             IF (unit_nr /= output_unit) THEN
    1716           49 :                INQUIRE (UNIT=unit_nr, NAME=filename)
    1717              :                WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
    1718           49 :                   "MOMENTS", "The electric/magnetic moments are written to file:", &
    1719           98 :                   TRIM(filename)
    1720          746 :             ELSE IF (moments_format /= moments_format_trajectory) THEN
    1721          743 :                WRITE (UNIT=output_unit, FMT="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS"
    1722              :             END IF
    1723              :          END IF
    1724              : 
    1725         1570 :          CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
    1726              : 
    1727         1570 :          IF (moments_format == moments_format_trajectory) THEN
    1728            6 :             IF (do_kpoints) THEN
    1729            0 :                CPABORT("MOMENTS FORMAT TRAJECTORY is not available for k-point calculations.")
    1730              :             END IF
    1731            6 :             IF (maxmom /= 1) CPABORT("MOMENTS FORMAT TRAJECTORY requires MAX_MOMENT 1.")
    1732            6 :             IF (magnetic) CPABORT("MOMENTS FORMAT TRAJECTORY does not support MAGNETIC moments.")
    1733            6 :             IF (vel_reprs) CPABORT("MOMENTS FORMAT TRAJECTORY does not support VEL_REPRS.")
    1734            6 :             IF (do_kg) CPABORT("MOMENTS FORMAT TRAJECTORY does not support KG moments.")
    1735            6 :             moments_unit_nr = -1
    1736              :          ELSE
    1737         1564 :             moments_unit_nr = unit_nr
    1738              :          END IF
    1739              : 
    1740         1570 :          IF (do_kpoints) THEN
    1741           10 :             CALL qs_moment_kpoints(qs_env, maxmom, reference, ref_point, max_nmo, moments_unit_nr)
    1742              :          ELSE
    1743         1560 :             IF (periodic) THEN
    1744          478 :                CALL qs_moment_berry_phase(qs_env, magnetic, maxmom, reference, ref_point, moments_unit_nr)
    1745              :             ELSE
    1746         1082 :                CALL qs_moment_locop(qs_env, magnetic, maxmom, reference, ref_point, moments_unit_nr, vel_reprs, com_nl)
    1747              :             END IF
    1748         1560 :             IF (do_kg) THEN
    1749            0 :                CALL calculate_kg_moments(qs_env, moments_unit_nr, maxmom, magnetic, vel_reprs, com_nl)
    1750              :             END IF
    1751              :          END IF
    1752         1570 :          IF (moments_format == moments_format_trajectory) THEN
    1753            6 :             CALL write_moments_trajectory(unit_nr, logger, qs_env, periodic, new_file, "MOMENTS|")
    1754              :          END IF
    1755              : 
    1756              :          CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, &
    1757         1570 :                                            basis_section=input, print_key_path="DFT%PRINT%MOMENTS")
    1758              : 
    1759         1570 :          IF (second_ref_point) THEN
    1760              :             reference = section_get_ival(section_vals=input, &
    1761            0 :                                          keyword_name="DFT%PRINT%MOMENTS%REFERENCE_2")
    1762              : 
    1763            0 :             NULLIFY (ref_point)
    1764            0 :             CALL section_vals_val_get(input, "DFT%PRINT%MOMENTS%REF_POINT_2", r_vals=ref_point)
    1765              :             unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=input, &
    1766              :                                            print_key_path="DFT%PRINT%MOMENTS", extension=".dat", &
    1767              :                                            middle_name="moments_refpoint_2", log_filename=.FALSE., &
    1768            0 :                                            is_new_file=new_file)
    1769              : 
    1770            0 :             IF (output_unit > 0) THEN
    1771            0 :                IF (unit_nr /= output_unit) THEN
    1772            0 :                   INQUIRE (UNIT=unit_nr, NAME=filename)
    1773              :                   WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
    1774            0 :                      "MOMENTS", "The electric/magnetic moments for the second reference point are written to file:", &
    1775            0 :                      TRIM(filename)
    1776            0 :                ELSE IF (moments_format /= moments_format_trajectory) THEN
    1777            0 :                   WRITE (UNIT=output_unit, FMT="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS"
    1778              :                END IF
    1779              :             END IF
    1780            0 :             IF (moments_format /= moments_format_trajectory) moments_unit_nr = unit_nr
    1781            0 :             IF (do_kpoints) THEN
    1782            0 :                CALL qs_moment_kpoints(qs_env, maxmom, reference, ref_point, max_nmo, moments_unit_nr)
    1783              :             ELSE
    1784            0 :                IF (periodic) THEN
    1785            0 :                   CALL qs_moment_berry_phase(qs_env, magnetic, maxmom, reference, ref_point, moments_unit_nr)
    1786              :                ELSE
    1787              :                   CALL qs_moment_locop(qs_env, magnetic, maxmom, reference, ref_point, &
    1788            0 :                                        moments_unit_nr, vel_reprs, com_nl)
    1789              :                END IF
    1790              :             END IF
    1791            0 :             IF (moments_format == moments_format_trajectory) THEN
    1792            0 :                CALL write_moments_trajectory(unit_nr, logger, qs_env, periodic, new_file, "MOMENTS_REF2|")
    1793              :             END IF
    1794              :             CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, &
    1795            0 :                                               basis_section=input, print_key_path="DFT%PRINT%MOMENTS")
    1796              :          END IF
    1797              : 
    1798              :       END IF
    1799              : 
    1800        13381 :       CALL timestop(handle)
    1801              : 
    1802        13381 :    END SUBROUTINE qs_scf_post_moments
    1803              : 
    1804              : ! **************************************************************************************************
    1805              : !> \brief Writes one machine-readable electric-dipole record for the current iteration.
    1806              : !> \param unit_nr output unit
    1807              : !> \param logger logger carrying the current iteration levels
    1808              : !> \param qs_env quantum environment containing the dipole result and cell
    1809              : !> \param periodic whether the Berry-phase operator was used
    1810              : !> \param new_file whether the output file has just been created
    1811              : !> \param label record label
    1812              : ! **************************************************************************************************
    1813            6 :    SUBROUTINE write_moments_trajectory(unit_nr, logger, qs_env, periodic, new_file, label)
    1814              :       INTEGER, INTENT(IN)                                :: unit_nr
    1815              :       TYPE(cp_logger_type), POINTER                      :: logger
    1816              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1817              :       LOGICAL, INTENT(IN)                                :: periodic, new_file
    1818              :       CHARACTER(LEN=*), INTENT(IN)                       :: label
    1819              : 
    1820              :       CHARACTER(LEN=default_string_length)               :: description, iter
    1821              :       REAL(KIND=dp), DIMENSION(3)                        :: dipole
    1822              :       TYPE(cell_type), POINTER                           :: cell
    1823              :       TYPE(cp_result_type), POINTER                      :: results
    1824              : 
    1825            6 :       IF (unit_nr <= 0) RETURN
    1826              : 
    1827            3 :       NULLIFY (cell, results)
    1828            3 :       CALL get_qs_env(qs_env, cell=cell, results=results)
    1829            3 :       description = "[DIPOLE]"
    1830            3 :       CALL get_results(results=results, description=description, values=dipole)
    1831              : 
    1832            3 :       IF (new_file) THEN
    1833            0 :          IF (periodic) THEN
    1834              :             WRITE (unit_nr, "(A)") "# "//TRIM(label)// &
    1835              :                " iter_level dipole_x dipole_y dipole_z dipole_norm cell_xx cell_xy cell_xz"// &
    1836            0 :                " cell_yx cell_yy cell_yz cell_zx cell_zy cell_zz [Debye]"
    1837              :          ELSE
    1838              :             WRITE (unit_nr, "(A)") "# "//TRIM(label)// &
    1839            0 :                " iter_level dipole_x dipole_y dipole_z dipole_norm [Debye]"
    1840              :          END IF
    1841              :       END IF
    1842              : 
    1843            3 :       iter = cp_iter_string(logger%iter_info)
    1844            3 :       IF (periodic) THEN
    1845            3 :          WRITE (unit_nr, "(1X,A,1X,A15,13(1X,ES18.10))") TRIM(label), iter(1:15), &
    1846           21 :             dipole*debye, SQRT(SUM(dipole**2))*debye, &
    1847           33 :             cell%hmat(1, :)*debye, cell%hmat(2, :)*debye, cell%hmat(3, :)*debye
    1848              :       ELSE
    1849            0 :          WRITE (unit_nr, "(1X,A,1X,A15,4(1X,ES18.10))") TRIM(label), iter(1:15), &
    1850            0 :             dipole*debye, SQRT(SUM(dipole**2))*debye
    1851              :       END IF
    1852              : 
    1853              :    END SUBROUTINE write_moments_trajectory
    1854              : 
    1855              : ! **************************************************************************************************
    1856              : !> \brief Computes and prints the X-ray diffraction spectrum.
    1857              : !> \param input ...
    1858              : !> \param dft_section ...
    1859              : !> \param logger ...
    1860              : !> \param qs_env the qs_env in which the qs_env lives
    1861              : !> \param output_unit ...
    1862              : ! **************************************************************************************************
    1863        12077 :    SUBROUTINE qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit)
    1864              : 
    1865              :       TYPE(section_vals_type), POINTER                   :: input, dft_section
    1866              :       TYPE(cp_logger_type), POINTER                      :: logger
    1867              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1868              :       INTEGER, INTENT(IN)                                :: output_unit
    1869              : 
    1870              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_scf_post_xray'
    1871              : 
    1872              :       CHARACTER(LEN=default_path_length)                 :: filename
    1873              :       INTEGER                                            :: handle, unit_nr
    1874              :       REAL(KIND=dp)                                      :: q_max
    1875              :       TYPE(section_vals_type), POINTER                   :: print_key
    1876              : 
    1877        12077 :       CALL timeset(routineN, handle)
    1878              : 
    1879              :       print_key => section_vals_get_subs_vals(section_vals=input, &
    1880        12077 :                                               subsection_name="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM")
    1881              : 
    1882        12077 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    1883              :          q_max = section_get_rval(section_vals=dft_section, &
    1884           30 :                                   keyword_name="PRINT%XRAY_DIFFRACTION_SPECTRUM%Q_MAX")
    1885              :          unit_nr = cp_print_key_unit_nr(logger=logger, &
    1886              :                                         basis_section=input, &
    1887              :                                         print_key_path="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM", &
    1888              :                                         extension=".dat", &
    1889              :                                         middle_name="xrd", &
    1890           30 :                                         log_filename=.FALSE.)
    1891           30 :          IF (output_unit > 0) THEN
    1892           15 :             INQUIRE (UNIT=unit_nr, NAME=filename)
    1893              :             WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
    1894           15 :                "X-RAY DIFFRACTION SPECTRUM"
    1895           15 :             IF (unit_nr /= output_unit) THEN
    1896              :                WRITE (UNIT=output_unit, FMT="(/,T3,A,/,/,T3,A,/)") &
    1897           14 :                   "The coherent X-ray diffraction spectrum is written to the file:", &
    1898           28 :                   TRIM(filename)
    1899              :             END IF
    1900              :          END IF
    1901              :          CALL xray_diffraction_spectrum(qs_env=qs_env, &
    1902              :                                         unit_number=unit_nr, &
    1903           30 :                                         q_max=q_max)
    1904              :          CALL cp_print_key_finished_output(unit_nr=unit_nr, &
    1905              :                                            logger=logger, &
    1906              :                                            basis_section=input, &
    1907           30 :                                            print_key_path="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM")
    1908              :       END IF
    1909              : 
    1910        12077 :       CALL timestop(handle)
    1911              : 
    1912        12077 :    END SUBROUTINE qs_scf_post_xray
    1913              : 
    1914              : ! **************************************************************************************************
    1915              : !> \brief Computes and prints Electric Field Gradient
    1916              : !> \param input ...
    1917              : !> \param logger ...
    1918              : !> \param qs_env the qs_env in which the qs_env lives
    1919              : ! **************************************************************************************************
    1920        12077 :    SUBROUTINE qs_scf_post_efg(input, logger, qs_env)
    1921              :       TYPE(section_vals_type), POINTER                   :: input
    1922              :       TYPE(cp_logger_type), POINTER                      :: logger
    1923              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1924              : 
    1925              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_scf_post_efg'
    1926              : 
    1927              :       INTEGER                                            :: handle
    1928              :       TYPE(section_vals_type), POINTER                   :: print_key
    1929              : 
    1930        12077 :       CALL timeset(routineN, handle)
    1931              : 
    1932              :       print_key => section_vals_get_subs_vals(section_vals=input, &
    1933        12077 :                                               subsection_name="DFT%PRINT%ELECTRIC_FIELD_GRADIENT")
    1934        12077 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
    1935              :                 cp_p_file)) THEN
    1936           30 :          CALL qs_efg_calc(qs_env=qs_env)
    1937              :       END IF
    1938              : 
    1939        12077 :       CALL timestop(handle)
    1940              : 
    1941        12077 :    END SUBROUTINE qs_scf_post_efg
    1942              : 
    1943              : ! **************************************************************************************************
    1944              : !> \brief Computes the Electron Transfer Coupling matrix element
    1945              : !> \param input ...
    1946              : !> \param qs_env the qs_env in which the qs_env lives
    1947              : !> \param dft_control ...
    1948              : ! **************************************************************************************************
    1949        24154 :    SUBROUTINE qs_scf_post_et(input, qs_env, dft_control)
    1950              :       TYPE(section_vals_type), POINTER                   :: input
    1951              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1952              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1953              : 
    1954              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_scf_post_et'
    1955              : 
    1956              :       INTEGER                                            :: handle, ispin
    1957              :       LOGICAL                                            :: do_et
    1958        12077 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: my_mos
    1959              :       TYPE(section_vals_type), POINTER                   :: et_section
    1960              : 
    1961        12077 :       CALL timeset(routineN, handle)
    1962              : 
    1963              :       do_et = .FALSE.
    1964        12077 :       et_section => section_vals_get_subs_vals(input, "PROPERTIES%ET_COUPLING")
    1965        12077 :       CALL section_vals_get(et_section, explicit=do_et)
    1966        12077 :       IF (do_et) THEN
    1967           10 :          IF (qs_env%et_coupling%first_run) THEN
    1968           10 :             NULLIFY (my_mos)
    1969           50 :             ALLOCATE (my_mos(dft_control%nspins))
    1970           50 :             ALLOCATE (qs_env%et_coupling%et_mo_coeff(dft_control%nspins))
    1971           30 :             DO ispin = 1, dft_control%nspins
    1972              :                CALL cp_fm_create(matrix=my_mos(ispin), &
    1973              :                                  matrix_struct=qs_env%mos(ispin)%mo_coeff%matrix_struct, &
    1974           20 :                                  name="FIRST_RUN_COEFF"//TRIM(ADJUSTL(cp_to_string(ispin)))//"MATRIX")
    1975              :                CALL cp_fm_to_fm(qs_env%mos(ispin)%mo_coeff, &
    1976           30 :                                 my_mos(ispin))
    1977              :             END DO
    1978           10 :             CALL set_et_coupling_type(qs_env%et_coupling, et_mo_coeff=my_mos)
    1979           10 :             DEALLOCATE (my_mos)
    1980              :          END IF
    1981              :       END IF
    1982              : 
    1983        12077 :       CALL timestop(handle)
    1984              : 
    1985        12077 :    END SUBROUTINE qs_scf_post_et
    1986              : 
    1987              : ! **************************************************************************************************
    1988              : !> \brief compute the electron localization function
    1989              : !>
    1990              : !> \param input ...
    1991              : !> \param logger ...
    1992              : !> \param qs_env ...
    1993              : !> \par History
    1994              : !>      2012-07 Created [MI]
    1995              : ! **************************************************************************************************
    1996        12077 :    SUBROUTINE qs_scf_post_elf(input, logger, qs_env)
    1997              :       TYPE(section_vals_type), POINTER                   :: input
    1998              :       TYPE(cp_logger_type), POINTER                      :: logger
    1999              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2000              : 
    2001              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_scf_post_elf'
    2002              : 
    2003              :       CHARACTER(LEN=default_path_length)                 :: filename, mpi_filename, my_pos_cube, &
    2004              :                                                             title
    2005              :       INTEGER                                            :: handle, ispin, output_unit, unit_nr
    2006              :       LOGICAL                                            :: append_cube, gapw, mpi_io
    2007              :       REAL(dp)                                           :: rho_cutoff
    2008        12077 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: zcharge
    2009              :       TYPE(cp_section_key)                               :: elf_section_key
    2010              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2011              :       TYPE(particle_list_type), POINTER                  :: particles
    2012              :       TYPE(pw_env_type), POINTER                         :: pw_env
    2013        12077 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
    2014              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    2015        12077 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)    :: elf_r
    2016              :       TYPE(qs_subsys_type), POINTER                      :: subsys
    2017              :       TYPE(section_vals_type), POINTER                   :: elf_section
    2018              : 
    2019        12077 :       CALL timeset(routineN, handle)
    2020        12077 :       output_unit = cp_logger_get_default_io_unit(logger)
    2021              : 
    2022        12077 :       elf_section_key = cube_or_openpmd(input, str_elf_cubes, str_elf_openpmd, logger)
    2023              : 
    2024        12077 :       elf_section => section_vals_get_subs_vals(input, elf_section_key%absolute_section_key)
    2025        12077 :       IF (elf_section_key%do_output) THEN
    2026              : 
    2027           82 :          NULLIFY (dft_control, pw_env, auxbas_pw_pool, pw_pools, particles, subsys)
    2028           82 :          CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, subsys=subsys)
    2029           82 :          CALL qs_subsys_get(subsys, particles=particles)
    2030              : 
    2031           82 :          gapw = dft_control%qs_control%gapw
    2032           82 :          IF (.NOT. gapw) THEN
    2033              :             ! allocate
    2034          330 :             ALLOCATE (elf_r(dft_control%nspins))
    2035              :             CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
    2036           82 :                             pw_pools=pw_pools)
    2037          166 :             DO ispin = 1, dft_control%nspins
    2038           84 :                CALL auxbas_pw_pool%create_pw(elf_r(ispin))
    2039          166 :                CALL pw_zero(elf_r(ispin))
    2040              :             END DO
    2041              : 
    2042           82 :             IF (output_unit > 0) THEN
    2043              :                WRITE (UNIT=output_unit, FMT="(/,T15,A,/)") &
    2044           41 :                   " ----- ELF is computed on the real space grid -----"
    2045              :             END IF
    2046           82 :             rho_cutoff = section_get_rval(elf_section, "density_cutoff")
    2047           82 :             CALL qs_elf_calc(qs_env, elf_r, rho_cutoff)
    2048           82 :             CALL get_effective_core_charges(qs_env, zcharge)
    2049              : 
    2050              :             ! write ELF into cube file
    2051              : 
    2052              :             ! For openPMD, refer to access modes instead of APPEND key
    2053           82 :             IF (elf_section_key%grid_output == grid_output_cubes) THEN
    2054           82 :                append_cube = section_get_lval(elf_section, "APPEND")
    2055              :             END IF
    2056           82 :             my_pos_cube = "REWIND"
    2057           82 :             IF (append_cube) THEN
    2058            0 :                my_pos_cube = "APPEND"
    2059              :             END IF
    2060              : 
    2061          166 :             DO ispin = 1, dft_control%nspins
    2062           84 :                WRITE (filename, '(a5,I1.1)') "ELF_S", ispin
    2063           84 :                WRITE (title, *) "ELF spin ", ispin
    2064           84 :                mpi_io = .TRUE.
    2065              :                unit_nr = elf_section_key%print_key_unit_nr( &
    2066              :                          logger, &
    2067              :                          input, &
    2068              :                          elf_section_key%absolute_section_key, &
    2069              :                          extension=".cube", &
    2070              :                          middle_name=TRIM(filename), &
    2071              :                          file_position=my_pos_cube, &
    2072              :                          log_filename=.FALSE., &
    2073              :                          mpi_io=mpi_io, &
    2074              :                          fout=mpi_filename, &
    2075              :                          openpmd_basename="dft-elf", &
    2076              :                          openpmd_unit_dimension=openpmd_unit_dimension_dimensionless, &
    2077              :                          openpmd_unit_si=openpmd_unit_si_dimensionless, &
    2078           84 :                          sim_time=qs_env%sim_time)
    2079           84 :                IF (output_unit > 0) THEN
    2080           42 :                   IF (.NOT. mpi_io) THEN
    2081            0 :                      INQUIRE (UNIT=unit_nr, NAME=filename)
    2082              :                   ELSE
    2083           42 :                      filename = mpi_filename
    2084              :                   END IF
    2085              :                   WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") &
    2086           42 :                      "ELF is written in "//elf_section_key%format_name//" file format to the file:", &
    2087           84 :                      TRIM(filename)
    2088              :                END IF
    2089              : 
    2090              :                CALL elf_section_key%write_pw(elf_r(ispin), unit_nr, title, particles=particles, zeff=zcharge, &
    2091           84 :                                              stride=section_get_ivals(elf_section, "STRIDE"), mpi_io=mpi_io)
    2092              :                CALL elf_section_key%print_key_finished_output( &
    2093              :                   unit_nr, &
    2094              :                   logger, &
    2095              :                   input, &
    2096              :                   elf_section_key%absolute_section_key, &
    2097           84 :                   mpi_io=mpi_io)
    2098              : 
    2099          166 :                CALL auxbas_pw_pool%give_back_pw(elf_r(ispin))
    2100              :             END DO
    2101              : 
    2102              :             ! deallocate
    2103           82 :             DEALLOCATE (elf_r, zcharge)
    2104              : 
    2105              :          ELSE
    2106              :             ! not implemented
    2107            0 :             CPWARN("ELF not implemented for GAPW calculations!")
    2108              :          END IF
    2109              : 
    2110              :       END IF ! print key
    2111              : 
    2112        12077 :       CALL timestop(handle)
    2113              : 
    2114        24154 :    END SUBROUTINE qs_scf_post_elf
    2115              : 
    2116              : ! **************************************************************************************************
    2117              : !> \brief computes the condition number of the overlap matrix and
    2118              : !>      prints the value of the total energy. This is needed
    2119              : !>      for BASIS_MOLOPT optimizations
    2120              : !> \param input ...
    2121              : !> \param logger ...
    2122              : !> \param qs_env the qs_env in which the qs_env lives
    2123              : !> \par History
    2124              : !>      2007-07 Created [Joost VandeVondele]
    2125              : ! **************************************************************************************************
    2126        12077 :    SUBROUTINE qs_scf_post_molopt(input, logger, qs_env)
    2127              :       TYPE(section_vals_type), POINTER                   :: input
    2128              :       TYPE(cp_logger_type), POINTER                      :: logger
    2129              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2130              : 
    2131              :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_molopt'
    2132              : 
    2133              :       INTEGER                                            :: handle, nao, unit_nr
    2134              :       REAL(KIND=dp)                                      :: S_cond_number
    2135        12077 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues
    2136              :       TYPE(cp_fm_struct_type), POINTER                   :: ao_ao_fmstruct
    2137              :       TYPE(cp_fm_type)                                   :: fm_s, fm_work
    2138              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    2139        12077 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
    2140        12077 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    2141              :       TYPE(qs_energy_type), POINTER                      :: energy
    2142              :       TYPE(section_vals_type), POINTER                   :: print_key
    2143              : 
    2144        12077 :       CALL timeset(routineN, handle)
    2145              : 
    2146              :       print_key => section_vals_get_subs_vals(section_vals=input, &
    2147        12077 :                                               subsection_name="DFT%PRINT%BASIS_MOLOPT_QUANTITIES")
    2148        12077 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
    2149              :                 cp_p_file)) THEN
    2150              : 
    2151           28 :          CALL get_qs_env(qs_env, energy=energy, matrix_s=matrix_s, mos=mos)
    2152              : 
    2153              :          ! set up the two needed full matrices, using mo_coeff as a template
    2154           28 :          CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, nao=nao)
    2155              :          CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, &
    2156              :                                   nrow_global=nao, ncol_global=nao, &
    2157           28 :                                   template_fmstruct=mo_coeff%matrix_struct)
    2158              :          CALL cp_fm_create(fm_s, matrix_struct=ao_ao_fmstruct, &
    2159           28 :                            name="fm_s")
    2160              :          CALL cp_fm_create(fm_work, matrix_struct=ao_ao_fmstruct, &
    2161           28 :                            name="fm_work")
    2162           28 :          CALL cp_fm_struct_release(ao_ao_fmstruct)
    2163           84 :          ALLOCATE (eigenvalues(nao))
    2164              : 
    2165           28 :          CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, fm_s)
    2166           28 :          CALL choose_eigv_solver(fm_s, fm_work, eigenvalues)
    2167              : 
    2168           28 :          CALL cp_fm_release(fm_s)
    2169           28 :          CALL cp_fm_release(fm_work)
    2170              : 
    2171          992 :          S_cond_number = MAXVAL(ABS(eigenvalues))/MAX(MINVAL(ABS(eigenvalues)), EPSILON(0.0_dp))
    2172              : 
    2173              :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%BASIS_MOLOPT_QUANTITIES", &
    2174           28 :                                         extension=".molopt")
    2175              : 
    2176           28 :          IF (unit_nr > 0) THEN
    2177              :             ! please keep this format fixed, needs to be grepable for molopt
    2178              :             ! optimizations
    2179           14 :             WRITE (unit_nr, '(T2,A28,2A25)') "", "Tot. Ener.", "S Cond. Numb."
    2180           14 :             WRITE (unit_nr, '(T2,A28,2E25.17)') "BASIS_MOLOPT_QUANTITIES", energy%total, S_cond_number
    2181              :          END IF
    2182              : 
    2183              :          CALL cp_print_key_finished_output(unit_nr, logger, input, &
    2184           84 :                                            "DFT%PRINT%BASIS_MOLOPT_QUANTITIES")
    2185              : 
    2186              :       END IF
    2187              : 
    2188        12077 :       CALL timestop(handle)
    2189              : 
    2190        24154 :    END SUBROUTINE qs_scf_post_molopt
    2191              : 
    2192              : ! **************************************************************************************************
    2193              : !> \brief Dumps EPR
    2194              : !> \param input ...
    2195              : !> \param logger ...
    2196              : !> \param qs_env the qs_env in which the qs_env lives
    2197              : ! **************************************************************************************************
    2198        12077 :    SUBROUTINE qs_scf_post_epr(input, logger, qs_env)
    2199              :       TYPE(section_vals_type), POINTER                   :: input
    2200              :       TYPE(cp_logger_type), POINTER                      :: logger
    2201              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2202              : 
    2203              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_scf_post_epr'
    2204              : 
    2205              :       INTEGER                                            :: handle
    2206              :       TYPE(section_vals_type), POINTER                   :: print_key
    2207              : 
    2208        12077 :       CALL timeset(routineN, handle)
    2209              : 
    2210              :       print_key => section_vals_get_subs_vals(section_vals=input, &
    2211        12077 :                                               subsection_name="DFT%PRINT%HYPERFINE_COUPLING_TENSOR")
    2212        12077 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
    2213              :                 cp_p_file)) THEN
    2214           30 :          CALL qs_epr_hyp_calc(qs_env=qs_env)
    2215              :       END IF
    2216              : 
    2217        12077 :       CALL timestop(handle)
    2218              : 
    2219        12077 :    END SUBROUTINE qs_scf_post_epr
    2220              : 
    2221              : ! **************************************************************************************************
    2222              : !> \brief Interface routine to trigger writing of results available from normal
    2223              : !>        SCF. Can write MO-dependent and MO free results (needed for call from
    2224              : !>        the linear scaling code)
    2225              : !> \param qs_env the qs_env in which the qs_env lives
    2226              : !> \param scf_env ...
    2227              : ! **************************************************************************************************
    2228        12077 :    SUBROUTINE write_available_results(qs_env, scf_env)
    2229              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2230              :       TYPE(qs_scf_env_type), OPTIONAL, POINTER           :: scf_env
    2231              : 
    2232              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results'
    2233              : 
    2234              :       INTEGER                                            :: handle
    2235              : 
    2236        12077 :       CALL timeset(routineN, handle)
    2237              : 
    2238              :       ! those properties that require MOs (not suitable density matrix based methods)
    2239        12077 :       CALL write_mo_dependent_results(qs_env, scf_env)
    2240              : 
    2241              :       ! those that depend only on the density matrix, they should be linear scaling in their implementation
    2242        12077 :       CALL write_mo_free_results(qs_env)
    2243              : 
    2244        12077 :       CALL timestop(handle)
    2245              : 
    2246        12077 :    END SUBROUTINE write_available_results
    2247              : 
    2248              : ! **************************************************************************************************
    2249              : !> \brief Write QS results available if MO's are present (if switched on through the print_keys)
    2250              : !>        Writes only MO dependent results. Split is necessary as ls_scf does not
    2251              : !>        provide MO's
    2252              : !> \param qs_env the qs_env in which the qs_env lives
    2253              : !> \param scf_env ...
    2254              : ! **************************************************************************************************
    2255        12425 :    SUBROUTINE write_mo_dependent_results(qs_env, scf_env)
    2256              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2257              :       TYPE(qs_scf_env_type), OPTIONAL, POINTER           :: scf_env
    2258              : 
    2259              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_mo_dependent_results'
    2260              : 
    2261              :       INTEGER                                            :: handle, homo, ispin, nlumo_dos, &
    2262              :                                                             nlumo_molden, nlumo_required, nlumos, &
    2263              :                                                             nmo, output_unit
    2264              :       LOGICAL                                            :: all_equal, defer_molden, do_curve, &
    2265              :                                                             do_dos, do_kpoints, do_pdos, &
    2266              :                                                             do_projected_dos, explicit
    2267              :       REAL(KIND=dp)                                      :: maxocc, s_square, s_square_ideal, &
    2268              :                                                             total_abs_spin_dens, total_spin_dens
    2269        12425 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues, occupation_numbers
    2270              :       TYPE(admm_type), POINTER                           :: admm_env
    2271        12425 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2272              :       TYPE(cell_type), POINTER                           :: cell
    2273        12425 :       TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER        :: unoccupied_evals
    2274        12425 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: unoccupied_orbs
    2275              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    2276              :       TYPE(cp_logger_type), POINTER                      :: logger
    2277        12425 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_rmpv, matrix_s
    2278              :       TYPE(dbcsr_type), POINTER                          :: mo_coeff_deriv
    2279              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2280        12425 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    2281        12425 :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
    2282              :       TYPE(particle_list_type), POINTER                  :: particles
    2283        12425 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2284              :       TYPE(pw_env_type), POINTER                         :: pw_env
    2285        12425 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
    2286              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    2287              :       TYPE(pw_r3d_rs_type)                               :: wf_r
    2288        12425 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
    2289              :       TYPE(qs_charges_type), POINTER                     :: qs_charges
    2290              :       TYPE(qs_energy_type), POINTER                      :: energy
    2291        12425 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2292              :       TYPE(qs_rho_type), POINTER                         :: rho
    2293              :       TYPE(qs_subsys_type), POINTER                      :: subsys
    2294              :       TYPE(scf_control_type), POINTER                    :: scf_control
    2295              :       TYPE(section_vals_type), POINTER                   :: casino_section, dft_section, &
    2296              :                                                             dos_section, input, sprint_section, &
    2297              :                                                             trexio_section
    2298              : 
    2299              : ! TYPE(kpoint_type), POINTER                         :: kpoints
    2300              : 
    2301        12425 :       CALL timeset(routineN, handle)
    2302              : 
    2303        12425 :       NULLIFY (cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, mo_coeff, &
    2304        12425 :                mo_coeff_deriv, mo_eigenvalues, mos, atomic_kind_set, qs_kind_set, &
    2305        12425 :                particle_set, rho, ks_rmpv, matrix_s, scf_control, dft_section, &
    2306        12425 :                molecule_set, input, particles, subsys, rho_r, unoccupied_orbs, &
    2307        12425 :                unoccupied_evals, casino_section, dos_section)
    2308              : 
    2309        12425 :       logger => cp_get_default_logger()
    2310        12425 :       output_unit = cp_logger_get_default_io_unit(logger)
    2311              : 
    2312        12425 :       CPASSERT(ASSOCIATED(qs_env))
    2313              :       CALL get_qs_env(qs_env, &
    2314              :                       dft_control=dft_control, &
    2315              :                       molecule_set=molecule_set, &
    2316              :                       atomic_kind_set=atomic_kind_set, &
    2317              :                       particle_set=particle_set, &
    2318              :                       qs_kind_set=qs_kind_set, &
    2319              :                       admm_env=admm_env, &
    2320              :                       scf_control=scf_control, &
    2321              :                       input=input, &
    2322              :                       cell=cell, &
    2323        12425 :                       subsys=subsys)
    2324        12425 :       CALL qs_subsys_get(subsys, particles=particles)
    2325        12425 :       CALL get_qs_env(qs_env, rho=rho)
    2326        12425 :       CALL qs_rho_get(rho, rho_r=rho_r)
    2327              : 
    2328              :       ! k points
    2329        12425 :       CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
    2330              : 
    2331              :       ! Write last MO information to output file if requested
    2332        12425 :       dft_section => section_vals_get_subs_vals(input, "DFT")
    2333        12425 :       IF (.NOT. qs_env%run_rtp) THEN
    2334        12077 :          CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.TRUE.)
    2335        12077 :          trexio_section => section_vals_get_subs_vals(dft_section, "PRINT%TREXIO")
    2336        12077 :          CALL section_vals_get(trexio_section, explicit=explicit)
    2337        12077 :          IF (explicit) THEN
    2338           10 :             CALL write_trexio(qs_env, trexio_section)
    2339              :          END IF
    2340        12077 :          casino_section => section_vals_get_subs_vals(dft_section, "PRINT%CASINO")
    2341        12077 :          CALL section_vals_get(casino_section, explicit=explicit)
    2342        12077 :          IF (explicit) THEN
    2343           10 :             CALL write_casino(qs_env, casino_section)
    2344              :          END IF
    2345        12077 :          sprint_section => section_vals_get_subs_vals(dft_section, "PRINT%MO_MOLDEN")
    2346        12077 :          defer_molden = .FALSE.
    2347        12077 :          IF (.NOT. do_kpoints) THEN
    2348        11409 :             CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv)
    2349        11409 :             CALL write_dm_binary_restart(mos, dft_section, ks_rmpv)
    2350              :             ! Check if molden write should be deferred for OT unoccupied orbitals
    2351        11409 :             CALL section_vals_val_get(sprint_section, "NLUMO", i_val=nlumo_molden)
    2352        11409 :             IF (nlumo_molden /= 0 .AND. PRESENT(scf_env)) THEN
    2353            0 :                IF (scf_env%method == ot_method_nr) defer_molden = .TRUE.
    2354              :             END IF
    2355              :             IF (.NOT. defer_molden) THEN
    2356              :                CALL write_mos_molden(mos, qs_kind_set, particle_set, sprint_section, cell=cell, &
    2357        11409 :                                      qs_env=qs_env, calc_energies=.TRUE.)
    2358              :             END IF
    2359              :             ! Write Chargemol .wfx
    2360        11409 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%CHARGEMOL"), &
    2361              :                       cp_p_file)) THEN
    2362            2 :                CALL write_wfx(qs_env, dft_section)
    2363              :             END IF
    2364              :          ELSE
    2365          668 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, sprint_section, ""), cp_p_file)) THEN
    2366            0 :                CPWARN("Molden format output is not possible for k-point calculations.")
    2367              :             END IF
    2368          668 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%CHARGEMOL"), &
    2369              :                       cp_p_file)) THEN
    2370            0 :                CPWARN("Chargemol .wfx format output is not possible for k-point calculations.")
    2371              :             END IF
    2372              :          END IF
    2373              : 
    2374              :          ! K-point MO wavefunction dump
    2375        12077 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_KP"), &
    2376              :                    cp_p_file)) THEN
    2377            0 :             IF (do_kpoints) THEN
    2378              :                CALL write_kpoint_mo_data(qs_env, &
    2379            0 :                                          section_vals_get_subs_vals(input, "DFT%PRINT%MO_KP"))
    2380              :             ELSE
    2381            0 :                CPWARN("MO_KP is only available for k-point calculations, ignored for Gamma-only")
    2382              :             END IF
    2383              :          END IF
    2384              : 
    2385        12077 :          dos_section => section_vals_get_subs_vals(dft_section, "PRINT%DOS")
    2386        12077 :          do_dos = BTEST(cp_print_key_should_output(logger%iter_info, dos_section), cp_p_file)
    2387        12077 :          CALL get_dos_pdos_flags(dos_section, do_dos, do_projected_dos, do_pdos, do_curve)
    2388              : 
    2389              :          ! For OT calculations, collect the largest request for additional unoccupied
    2390              :          ! orbitals among DOS, PDOS, and Molden, and generate them only once.
    2391        12077 :          nlumo_required = 0
    2392        12077 :          IF (.NOT. do_kpoints .AND. PRESENT(scf_env)) THEN
    2393        11409 :             IF (scf_env%method == ot_method_nr) THEN
    2394         4385 :                IF (do_dos) THEN
    2395            8 :                   CALL section_vals_val_get(dft_section, "PRINT%DOS%NLUMO", i_val=nlumo_dos)
    2396            8 :                   IF (nlumo_dos == -1) THEN
    2397            0 :                      nlumo_required = -1
    2398              :                   ELSE IF (nlumo_required /= -1) THEN
    2399            8 :                      nlumo_required = MAX(nlumo_required, nlumo_dos)
    2400              :                   END IF
    2401              :                END IF
    2402              : 
    2403         4385 :                IF (defer_molden) THEN
    2404            0 :                   IF (nlumo_molden == -1) THEN
    2405            0 :                      nlumo_required = -1
    2406            0 :                   ELSE IF (nlumo_required /= -1) THEN
    2407            0 :                      nlumo_required = MAX(nlumo_required, nlumo_molden)
    2408              :                   END IF
    2409              :                END IF
    2410         4385 :                IF (nlumo_required /= 0) THEN
    2411            8 :                   ALLOCATE (unoccupied_orbs(dft_control%nspins))
    2412            8 :                   ALLOCATE (unoccupied_evals(dft_control%nspins))
    2413              :                   CALL make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, &
    2414            2 :                                      nlumo_required, nlumos)
    2415              :                END IF
    2416              : 
    2417         4385 :                IF (do_dos .OR. do_projected_dos) THEN
    2418           16 :                   DO ispin = 1, dft_control%nspins
    2419              :                      ! With ADMM, temporarily modify the Kohn-Sham matrix while updating the
    2420              :                      ! eigenvalues and rotating the occupied OT subspace.
    2421            8 :                      IF (dft_control%do_admm) THEN
    2422            0 :                         CALL admm_correct_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
    2423              :                      END IF
    2424              :                      CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
    2425            8 :                                      eigenvalues=mo_eigenvalues)
    2426            8 :                      IF (ASSOCIATED(qs_env%mo_derivs)) THEN
    2427            8 :                         mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix
    2428              :                      ELSE
    2429            0 :                         mo_coeff_deriv => NULL()
    2430              :                      END IF
    2431              :                      CALL calculate_subspace_eigenvalues(mo_coeff, ks_rmpv(ispin)%matrix, mo_eigenvalues, &
    2432              :                                                          do_rotation=.TRUE., &
    2433            8 :                                                          co_rotate_dbcsr=mo_coeff_deriv)
    2434            8 :                      CALL set_mo_occupation(mo_set=mos(ispin))
    2435         4393 :                      IF (dft_control%do_admm) THEN
    2436            0 :                         CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix)
    2437              :                      END IF
    2438              :                   END DO
    2439              :                END IF
    2440              :             END IF
    2441              :          END IF
    2442              : 
    2443        12077 :          IF (defer_molden) THEN
    2444            0 :             IF (ASSOCIATED(unoccupied_orbs)) THEN
    2445            0 :                IF (output_unit > 0) THEN
    2446              :                   WRITE (output_unit, '(/,T2,A,I6,A)') &
    2447            0 :                      "MO_MOLDEN| Writing ", nlumos, " unoccupied orbitals to molden file"
    2448              :                END IF
    2449              :                CALL write_mos_molden(mos, qs_kind_set, particle_set, sprint_section, cell=cell, &
    2450              :                                      unoccupied_orbs=unoccupied_orbs, &
    2451              :                                      unoccupied_evals=unoccupied_evals, &
    2452            0 :                                      qs_env=qs_env, calc_energies=.TRUE.)
    2453              :             END IF
    2454              :          END IF
    2455              : 
    2456              :          ! DOS printout after the SCF cycle is completed
    2457        12077 :          IF (do_dos) THEN
    2458           64 :             IF (do_kpoints) THEN
    2459            4 :                CALL calculate_dos_kp(qs_env, dft_section)
    2460            4 :                IF (do_curve) CALL calculate_dos_kp(qs_env, dft_section, write_curve_output=.TRUE.)
    2461              :             ELSE
    2462           60 :                CALL get_qs_env(qs_env, mos=mos)
    2463           60 :                IF (ASSOCIATED(unoccupied_evals)) THEN
    2464              :                   CALL calculate_dos(mos, dft_section, unoccupied_evals=unoccupied_evals, &
    2465            2 :                                      smearing_enabled=dft_control%smear)
    2466            2 :                   IF (do_curve) CALL calculate_dos(mos, dft_section, unoccupied_evals=unoccupied_evals, &
    2467            0 :                                                    smearing_enabled=dft_control%smear, write_curve_output=.TRUE.)
    2468              :                ELSE
    2469           58 :                   CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear)
    2470           58 :                   IF (do_curve) CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear, &
    2471            0 :                                                    write_curve_output=.TRUE.)
    2472              :                END IF
    2473              :             END IF
    2474              :          END IF
    2475              : 
    2476              :          ! Print projected density-of-states outputs
    2477        12077 :          IF (do_projected_dos) THEN
    2478           22 :             IF (do_kpoints) THEN
    2479              :                CALL calculate_projected_dos_kp(qs_env, dft_section, pdos_print_key="PRINT%DOS", &
    2480            2 :                                                write_pdos=do_pdos, write_pdos_curve=do_curve)
    2481              :             ELSE
    2482              :                CALL get_qs_env(qs_env, &
    2483              :                                mos=mos, &
    2484           20 :                                matrix_ks=ks_rmpv)
    2485           40 :                DO ispin = 1, dft_control%nspins
    2486           40 :                   IF (dft_control%nspins == 2) THEN
    2487            0 :                      IF (ASSOCIATED(unoccupied_orbs)) THEN
    2488              :                         CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
    2489              :                                                      qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin, &
    2490              :                                                      unoccupied_orbs=unoccupied_orbs(ispin), &
    2491              :                                                      unoccupied_evals=unoccupied_evals(ispin), &
    2492            0 :                                                      pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
    2493              :                      ELSE
    2494              :                         CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
    2495              :                                                      qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin, &
    2496            0 :                                                      pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
    2497              :                      END IF
    2498              :                   ELSE
    2499           20 :                      IF (ASSOCIATED(unoccupied_orbs)) THEN
    2500              :                         CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
    2501              :                                                      qs_kind_set, particle_set, qs_env, dft_section, &
    2502              :                                                      unoccupied_orbs=unoccupied_orbs(ispin), &
    2503              :                                                      unoccupied_evals=unoccupied_evals(ispin), &
    2504            2 :                                                      pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
    2505              :                      ELSE
    2506              :                         CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
    2507              :                                                      qs_kind_set, particle_set, qs_env, dft_section, &
    2508           18 :                                                      pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
    2509              :                      END IF
    2510              :                   END IF
    2511              :                END DO
    2512              :             END IF
    2513              :          END IF
    2514        12077 :          IF (ASSOCIATED(unoccupied_orbs)) THEN
    2515            4 :             DO ispin = 1, dft_control%nspins
    2516            2 :                DEALLOCATE (unoccupied_evals(ispin)%array)
    2517            4 :                CALL cp_fm_release(unoccupied_orbs(ispin))
    2518              :             END DO
    2519            2 :             DEALLOCATE (unoccupied_evals)
    2520            2 :             DEALLOCATE (unoccupied_orbs)
    2521              :          END IF
    2522              :       END IF
    2523              : 
    2524              :       ! Integrated absolute spin density and spin contamination ***
    2525        12425 :       IF (dft_control%nspins == 2) THEN
    2526         2344 :          total_spin_dens = 0.0_dp
    2527         2344 :          total_abs_spin_dens = 0.0_dp
    2528         2344 :          IF (dft_control%qs_control%gapw) THEN
    2529          368 :             CALL get_qs_env(qs_env, qs_charges=qs_charges)
    2530              :             total_spin_dens = qs_charges%total_rho_hard_spin - &
    2531          368 :                               qs_charges%total_rho_soft_spin
    2532              :             total_abs_spin_dens = qs_charges%total_rho_hard_abs_spin - &
    2533          368 :                                   qs_charges%total_rho_soft_abs_spin
    2534              :          END IF
    2535         2344 :          CALL get_qs_env(qs_env, mos=mos)
    2536         2344 :          CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
    2537              :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
    2538         2344 :                          pw_pools=pw_pools)
    2539         2344 :          CALL auxbas_pw_pool%create_pw(wf_r)
    2540         2344 :          CALL pw_copy(rho_r(1), wf_r)
    2541         2344 :          CALL pw_axpy(rho_r(2), wf_r, alpha=-1._dp)
    2542         2344 :          total_spin_dens = total_spin_dens + pw_integrate_function(wf_r)
    2543         2344 :          IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,F20.10))') &
    2544         1195 :             "Integrated spin density: ", total_spin_dens
    2545         2344 :          total_abs_spin_dens = total_abs_spin_dens + pw_integrate_function(wf_r, oprt="ABS")
    2546         2344 :          IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T3,A,T61,F20.10))') &
    2547         1195 :             "Integrated absolute spin density: ", total_abs_spin_dens
    2548         2344 :          CALL auxbas_pw_pool%give_back_pw(wf_r)
    2549              :          !
    2550              :          ! XXX Fix Me XXX
    2551              :          ! should be extended to the case where added MOs are present
    2552              :          ! should be extended to the k-point case
    2553              :          !
    2554         2344 :          IF (.NOT. do_kpoints) THEN
    2555         2264 :             all_equal = .TRUE.
    2556         6792 :             DO ispin = 1, dft_control%nspins
    2557              :                CALL get_mo_set(mo_set=mos(ispin), &
    2558              :                                occupation_numbers=occupation_numbers, &
    2559              :                                homo=homo, &
    2560              :                                nmo=nmo, &
    2561         4528 :                                maxocc=maxocc)
    2562         6792 :                IF (nmo > 0) THEN
    2563              :                   all_equal = all_equal .AND. &
    2564              :                               (ALL(occupation_numbers(1:homo) == maxocc) .AND. &
    2565        25558 :                                ALL(occupation_numbers(homo + 1:nmo) == 0.0_dp))
    2566              :                END IF
    2567              :             END DO
    2568         2264 :             IF (all_equal) THEN
    2569              :                CALL get_qs_env(qs_env=qs_env, &
    2570              :                                matrix_s=matrix_s, &
    2571         2146 :                                energy=energy)
    2572              :                CALL compute_s_square(mos=mos, matrix_s=matrix_s, s_square=s_square, &
    2573         2146 :                                      s_square_ideal=s_square_ideal)
    2574         2146 :                IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(T3,A,T51,2F15.6)') &
    2575         1096 :                   "Ideal and single determinant S**2 : ", s_square_ideal, s_square
    2576         2146 :                energy%s_square = s_square
    2577              :             END IF
    2578              :          END IF
    2579              :       END IF
    2580              : 
    2581        12425 :       CALL timestop(handle)
    2582              : 
    2583        12425 :    END SUBROUTINE write_mo_dependent_results
    2584              : 
    2585              : ! **************************************************************************************************
    2586              : !> \brief Write QS results always available (if switched on through the print_keys)
    2587              : !>        Can be called from ls_scf
    2588              : !> \param qs_env the qs_env in which the qs_env lives
    2589              : ! **************************************************************************************************
    2590        13441 :    SUBROUTINE write_mo_free_results(qs_env)
    2591              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2592              : 
    2593              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_mo_free_results'
    2594              :       CHARACTER(len=1), DIMENSION(3), PARAMETER          :: cdir = ["x", "y", "z"]
    2595              : 
    2596              :       CHARACTER(LEN=2)                                   :: element_symbol
    2597              :       CHARACTER(LEN=default_path_length)                 :: filename, mpi_filename, my_pos_cube, &
    2598              :                                                             my_pos_voro
    2599              :       CHARACTER(LEN=default_string_length)               :: name, print_density
    2600              :       INTEGER :: after, handle, i, iat, id, ikind, img, iso, ispin, iw, l, n_rep_hf, natom, nd(3), &
    2601              :          ngto, niso, nkind, np, nr, output_unit, print_level, should_print_bqb, should_print_voro, &
    2602              :          unit_nr, unit_nr_voro
    2603              :       LOGICAL :: append_cube, append_voro, do_hfx, do_kpoints, mpi_io, omit_headers, print_it, &
    2604              :          rho_r_valid, voro_print_txt, write_ks, write_xc, xrd_interface
    2605              :       REAL(KIND=dp)                                      :: norm_factor, q_max, rho_hard, rho_soft, &
    2606              :                                                             rho_total, rho_total_rspace, udvol, &
    2607              :                                                             volume
    2608        13441 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: zcharge
    2609        13441 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: bfun
    2610        13441 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: aedens, ccdens, ppdens
    2611              :       REAL(KIND=dp), DIMENSION(3)                        :: checksum_hr, dr
    2612        13441 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: my_Q0
    2613        13441 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2614              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    2615              :       TYPE(cell_type), POINTER                           :: cell
    2616              :       TYPE(cp_logger_type), POINTER                      :: logger
    2617              :       TYPE(cp_section_key)                               :: e_density_section
    2618        13441 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_hr
    2619        13441 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ks_rmpv, matrix_vxc, rho_ao
    2620              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2621              :       TYPE(grid_atom_type), POINTER                      :: grid_atom
    2622              :       TYPE(iao_env_type)                                 :: iao_env
    2623              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2624              :       TYPE(particle_list_type), POINTER                  :: particles
    2625        13441 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2626              :       TYPE(pw_c1d_gs_type)                               :: aux_g, rho_elec_gspace
    2627              :       TYPE(pw_c1d_gs_type), POINTER                      :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
    2628              :       TYPE(pw_env_type), POINTER                         :: pw_env
    2629        13441 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
    2630              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    2631              :       TYPE(pw_r3d_rs_type)                               :: aux_r, rho_elec_rspace, wf_r
    2632        13441 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
    2633              :       TYPE(pw_r3d_rs_type), POINTER                      :: mb_rho, v_hartree_rspace, vee
    2634        13441 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2635              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    2636              :       TYPE(qs_rho_type), POINTER                         :: rho
    2637              :       TYPE(qs_subsys_type), POINTER                      :: subsys
    2638              :       TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
    2639        13441 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set
    2640              :       TYPE(rho_atom_type), POINTER                       :: rho_atom
    2641              :       TYPE(section_vals_type), POINTER                   :: dft_section, hfx_section, input, &
    2642              :                                                             print_key, print_key_bqb, &
    2643              :                                                             print_key_voro, xc_section
    2644              : 
    2645        13441 :       CALL timeset(routineN, handle)
    2646        13441 :       NULLIFY (cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, hfx_section, &
    2647        13441 :                atomic_kind_set, qs_kind_set, particle_set, rho, ks_rmpv, rho_ao, rho_r, &
    2648        13441 :                dft_section, xc_section, input, particles, subsys, matrix_vxc, v_hartree_rspace, &
    2649        13441 :                vee)
    2650              : 
    2651        13441 :       logger => cp_get_default_logger()
    2652        13441 :       output_unit = cp_logger_get_default_io_unit(logger)
    2653              : 
    2654        13441 :       CPASSERT(ASSOCIATED(qs_env))
    2655              :       CALL get_qs_env(qs_env, &
    2656              :                       atomic_kind_set=atomic_kind_set, &
    2657              :                       qs_kind_set=qs_kind_set, &
    2658              :                       nkind=nkind, &
    2659              :                       natom=natom, &
    2660              :                       particle_set=particle_set, &
    2661              :                       cell=cell, &
    2662              :                       para_env=para_env, &
    2663              :                       dft_control=dft_control, &
    2664              :                       input=input, &
    2665              :                       do_kpoints=do_kpoints, &
    2666        13441 :                       subsys=subsys)
    2667        13441 :       dft_section => section_vals_get_subs_vals(input, "DFT")
    2668        13441 :       CALL qs_subsys_get(subsys, particles=particles)
    2669              : 
    2670        13441 :       CALL get_qs_env(qs_env, rho=rho)
    2671        13441 :       CALL qs_rho_get(rho, rho_r=rho_r)
    2672              : 
    2673        13441 :       CALL get_effective_core_charges(qs_env, zcharge)
    2674              : 
    2675              :       ! Print the total density (electronic + core charge)
    2676        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
    2677              :                                            "DFT%PRINT%TOT_DENSITY_CUBE"), cp_p_file)) THEN
    2678           82 :          NULLIFY (rho_core, rho0_s_gs, rhoz_cneo_s_gs)
    2679           82 :          append_cube = section_get_lval(input, "DFT%PRINT%TOT_DENSITY_CUBE%APPEND")
    2680           82 :          my_pos_cube = "REWIND"
    2681           82 :          IF (append_cube) THEN
    2682            0 :             my_pos_cube = "APPEND"
    2683              :          END IF
    2684              : 
    2685              :          CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho_core=rho_core, &
    2686           82 :                          rho0_s_gs=rho0_s_gs, rhoz_cneo_s_gs=rhoz_cneo_s_gs)
    2687              :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
    2688           82 :                          pw_pools=pw_pools)
    2689           82 :          CALL auxbas_pw_pool%create_pw(wf_r)
    2690           82 :          IF (dft_control%qs_control%gapw) THEN
    2691            0 :             IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
    2692            0 :                CALL pw_axpy(rho_core, rho0_s_gs)
    2693            0 :                IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
    2694            0 :                   CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
    2695              :                END IF
    2696            0 :                CALL pw_transfer(rho0_s_gs, wf_r)
    2697            0 :                CALL pw_axpy(rho_core, rho0_s_gs, -1.0_dp)
    2698            0 :                IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
    2699            0 :                   CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
    2700              :                END IF
    2701              :             ELSE
    2702            0 :                IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
    2703            0 :                   CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
    2704              :                END IF
    2705            0 :                CALL pw_transfer(rho0_s_gs, wf_r)
    2706            0 :                IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
    2707            0 :                   CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
    2708              :                END IF
    2709              :             END IF
    2710              :          ELSE
    2711           82 :             CALL pw_transfer(rho_core, wf_r)
    2712              :          END IF
    2713          164 :          DO ispin = 1, dft_control%nspins
    2714          164 :             CALL pw_axpy(rho_r(ispin), wf_r)
    2715              :          END DO
    2716           82 :          filename = "TOTAL_DENSITY"
    2717           82 :          mpi_io = .TRUE.
    2718              :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%TOT_DENSITY_CUBE", &
    2719              :                                         extension=".cube", middle_name=TRIM(filename), file_position=my_pos_cube, &
    2720           82 :                                         log_filename=.FALSE., mpi_io=mpi_io)
    2721              :          CALL cp_pw_to_cube(wf_r, unit_nr, "TOTAL DENSITY", &
    2722              :                             particles=particles, zeff=zcharge, &
    2723              :                             stride=section_get_ivals(dft_section, "PRINT%TOT_DENSITY_CUBE%STRIDE"), &
    2724              :                             max_file_size_mb=section_get_rval(dft_section, "PRINT%TOT_DENSITY_CUBE%MAX_FILE_SIZE_MB"), &
    2725           82 :                             mpi_io=mpi_io)
    2726              :          CALL cp_print_key_finished_output(unit_nr, logger, input, &
    2727           82 :                                            "DFT%PRINT%TOT_DENSITY_CUBE", mpi_io=mpi_io)
    2728           82 :          CALL auxbas_pw_pool%give_back_pw(wf_r)
    2729              :       END IF
    2730              : 
    2731        13441 :       e_density_section = cube_or_openpmd(input, str_e_density_cubes, str_e_density_openpmd, logger)
    2732              : 
    2733              :       ! Write cube file with electron density
    2734        13441 :       IF (e_density_section%do_output) THEN
    2735              :          CALL section_vals_val_get(dft_section, &
    2736              :                                    keyword_name=e_density_section%concat_to_relative("%DENSITY_INCLUDE"), &
    2737          154 :                                    c_val=print_density)
    2738              :          print_density = TRIM(print_density)
    2739              :          ! For openPMD, refer to access modes instead of APPEND key
    2740          154 :          IF (e_density_section%grid_output == grid_output_cubes) THEN
    2741          154 :             append_cube = section_get_lval(input, e_density_section%concat_to_absolute("%APPEND"))
    2742              :          END IF
    2743          154 :          my_pos_cube = "REWIND"
    2744          154 :          IF (append_cube) THEN
    2745            0 :             my_pos_cube = "APPEND"
    2746              :          END IF
    2747              :          ! Write the info on core densities for the interface between cp2k and the XRD code
    2748              :          ! together with the valence density they are used to compute the form factor (Fourier transform)
    2749          154 :          IF (e_density_section%grid_output == grid_output_cubes) THEN
    2750          154 :             xrd_interface = section_get_lval(input, e_density_section%concat_to_absolute("%XRD_INTERFACE"))
    2751              :          ELSE
    2752              :             ! Unimplemented for openPMD, since this does not use the regular routines
    2753              :             xrd_interface = .FALSE.
    2754              :          END IF
    2755              : 
    2756          154 :          IF (xrd_interface) THEN
    2757              :             !cube file only contains soft density (GAPW)
    2758            2 :             IF (dft_control%qs_control%gapw) print_density = "SOFT_DENSITY"
    2759              : 
    2760            2 :             filename = "ELECTRON_DENSITY"
    2761              :             unit_nr = cp_print_key_unit_nr(logger, input, e_density_section%absolute_section_key, &
    2762              :                                            extension=".xrd", middle_name=TRIM(filename), &
    2763            2 :                                            file_position=my_pos_cube, log_filename=.FALSE.)
    2764            2 :             ngto = section_get_ival(input, e_density_section%concat_to_absolute("%NGAUSS"))
    2765            2 :             IF (output_unit > 0) THEN
    2766            1 :                INQUIRE (UNIT=unit_nr, NAME=filename)
    2767              :                WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") &
    2768            1 :                   "The electron density (atomic part) is written to the file:", &
    2769            2 :                   TRIM(filename)
    2770              :             END IF
    2771              : 
    2772            2 :             xc_section => section_vals_get_subs_vals(input, "DFT%XC")
    2773            2 :             nkind = SIZE(atomic_kind_set)
    2774            2 :             IF (unit_nr > 0) THEN
    2775            1 :                WRITE (unit_nr, *) "Atomic (core) densities"
    2776            1 :                WRITE (unit_nr, *) "Unit cell"
    2777            1 :                WRITE (unit_nr, FMT="(3F20.12)") cell%hmat(1, 1), cell%hmat(1, 2), cell%hmat(1, 3)
    2778            1 :                WRITE (unit_nr, FMT="(3F20.12)") cell%hmat(2, 1), cell%hmat(2, 2), cell%hmat(2, 3)
    2779            1 :                WRITE (unit_nr, FMT="(3F20.12)") cell%hmat(3, 1), cell%hmat(3, 2), cell%hmat(3, 3)
    2780            1 :                WRITE (unit_nr, *) "Atomic types"
    2781            1 :                WRITE (unit_nr, *) nkind
    2782              :             END IF
    2783              :             ! calculate atomic density and core density
    2784           16 :             ALLOCATE (ppdens(ngto, 2, nkind), aedens(ngto, 2, nkind), ccdens(ngto, 2, nkind))
    2785            6 :             DO ikind = 1, nkind
    2786            4 :                atomic_kind => atomic_kind_set(ikind)
    2787            4 :                qs_kind => qs_kind_set(ikind)
    2788            4 :                CALL get_atomic_kind(atomic_kind, name=name, element_symbol=element_symbol)
    2789              :                CALL calculate_atomic_density(ppdens(:, :, ikind), atomic_kind, qs_kind, ngto, &
    2790            4 :                                              iunit=output_unit, confine=.TRUE.)
    2791              :                CALL calculate_atomic_density(aedens(:, :, ikind), atomic_kind, qs_kind, ngto, &
    2792            4 :                                              iunit=output_unit, allelectron=.TRUE., confine=.TRUE.)
    2793           52 :                ccdens(:, 1, ikind) = aedens(:, 1, ikind)
    2794           52 :                ccdens(:, 2, ikind) = 0._dp
    2795              :                CALL project_function_a(ccdens(1:ngto, 2, ikind), ccdens(1:ngto, 1, ikind), &
    2796            4 :                                        ppdens(1:ngto, 2, ikind), ppdens(1:ngto, 1, ikind), 0)
    2797           52 :                ccdens(:, 2, ikind) = aedens(:, 2, ikind) - ccdens(:, 2, ikind)
    2798            4 :                IF (unit_nr > 0) THEN
    2799            2 :                   WRITE (unit_nr, FMT="(I6,A10,A20)") ikind, TRIM(element_symbol), TRIM(name)
    2800            2 :                   WRITE (unit_nr, FMT="(I6)") ngto
    2801            2 :                   WRITE (unit_nr, *) "   Total density"
    2802           26 :                   WRITE (unit_nr, FMT="(2G24.12)") (aedens(i, 1, ikind), aedens(i, 2, ikind), i=1, ngto)
    2803            2 :                   WRITE (unit_nr, *) "    Core density"
    2804           26 :                   WRITE (unit_nr, FMT="(2G24.12)") (ccdens(i, 1, ikind), ccdens(i, 2, ikind), i=1, ngto)
    2805              :                END IF
    2806            6 :                NULLIFY (atomic_kind)
    2807              :             END DO
    2808              : 
    2809            2 :             IF (dft_control%qs_control%gapw) THEN
    2810            2 :                CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set)
    2811              : 
    2812            2 :                IF (unit_nr > 0) THEN
    2813            1 :                   WRITE (unit_nr, *) "Coordinates and GAPW density"
    2814              :                END IF
    2815            2 :                np = particles%n_els
    2816            6 :                DO iat = 1, np
    2817            4 :                   CALL get_atomic_kind(particles%els(iat)%atomic_kind, kind_number=ikind)
    2818            4 :                   CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid_atom)
    2819            4 :                   rho_atom => rho_atom_set(iat)
    2820            4 :                   IF (ASSOCIATED(rho_atom%rho_rad_h(1)%r_coef)) THEN
    2821            2 :                      nr = SIZE(rho_atom%rho_rad_h(1)%r_coef, 1)
    2822            2 :                      niso = SIZE(rho_atom%rho_rad_h(1)%r_coef, 2)
    2823              :                   ELSE
    2824            2 :                      nr = 0
    2825            2 :                      niso = 0
    2826              :                   END IF
    2827            4 :                   CALL para_env%sum(nr)
    2828            4 :                   CALL para_env%sum(niso)
    2829              : 
    2830           16 :                   ALLOCATE (bfun(nr, niso))
    2831            4 :                   bfun = 0._dp
    2832            8 :                   DO ispin = 1, dft_control%nspins
    2833            8 :                      IF (ASSOCIATED(rho_atom%rho_rad_h(1)%r_coef)) THEN
    2834          920 :                         bfun(:, :) = bfun + rho_atom%rho_rad_h(ispin)%r_coef - rho_atom%rho_rad_s(ispin)%r_coef
    2835              :                      END IF
    2836              :                   END DO
    2837            4 :                   CALL para_env%sum(bfun)
    2838           52 :                   ccdens(:, 1, ikind) = ppdens(:, 1, ikind)
    2839           52 :                   ccdens(:, 2, ikind) = 0._dp
    2840            4 :                   IF (unit_nr > 0) THEN
    2841            8 :                      WRITE (unit_nr, '(I10,I5,3f12.6)') iat, ikind, particles%els(iat)%r
    2842              :                   END IF
    2843           40 :                   DO iso = 1, niso
    2844           36 :                      l = indso(1, iso)
    2845           36 :                      CALL project_function_b(ccdens(:, 2, ikind), ccdens(:, 1, ikind), bfun(:, iso), grid_atom, l)
    2846           40 :                      IF (unit_nr > 0) THEN
    2847           18 :                         WRITE (unit_nr, FMT="(3I6)") iso, l, ngto
    2848          234 :                         WRITE (unit_nr, FMT="(2G24.12)") (ccdens(i, 1, ikind), ccdens(i, 2, ikind), i=1, ngto)
    2849              :                      END IF
    2850              :                   END DO
    2851           10 :                   DEALLOCATE (bfun)
    2852              :                END DO
    2853              :             ELSE
    2854            0 :                IF (unit_nr > 0) THEN
    2855            0 :                   WRITE (unit_nr, *) "Coordinates"
    2856            0 :                   np = particles%n_els
    2857            0 :                   DO iat = 1, np
    2858            0 :                      CALL get_atomic_kind(particles%els(iat)%atomic_kind, kind_number=ikind)
    2859            0 :                      WRITE (unit_nr, '(I10,I5,3f12.6)') iat, ikind, particles%els(iat)%r
    2860              :                   END DO
    2861              :                END IF
    2862              :             END IF
    2863              : 
    2864            2 :             DEALLOCATE (ppdens, aedens, ccdens)
    2865              : 
    2866              :             CALL cp_print_key_finished_output(unit_nr, logger, input, &
    2867            2 :                                               e_density_section%absolute_section_key)
    2868              : 
    2869              :          END IF
    2870          154 :          IF (dft_control%qs_control%gapw .AND. print_density == "TOTAL_DENSITY") THEN
    2871              :             ! total density in g-space not implemented for k-points
    2872            8 :             CPASSERT(.NOT. do_kpoints)
    2873              :             ! Print total electronic density
    2874              :             CALL get_qs_env(qs_env=qs_env, &
    2875            8 :                             pw_env=pw_env)
    2876              :             CALL pw_env_get(pw_env=pw_env, &
    2877              :                             auxbas_pw_pool=auxbas_pw_pool, &
    2878            8 :                             pw_pools=pw_pools)
    2879            8 :             CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
    2880            8 :             CALL pw_zero(rho_elec_rspace)
    2881            8 :             CALL auxbas_pw_pool%create_pw(pw=rho_elec_gspace)
    2882            8 :             CALL pw_zero(rho_elec_gspace)
    2883              :             CALL get_pw_grid_info(pw_grid=rho_elec_gspace%pw_grid, &
    2884              :                                   dr=dr, &
    2885            8 :                                   vol=volume)
    2886           32 :             q_max = SQRT(SUM((pi/dr(:))**2))
    2887              :             CALL calculate_rhotot_elec_gspace(qs_env=qs_env, &
    2888              :                                               auxbas_pw_pool=auxbas_pw_pool, &
    2889              :                                               rhotot_elec_gspace=rho_elec_gspace, &
    2890              :                                               q_max=q_max, &
    2891              :                                               rho_hard=rho_hard, &
    2892            8 :                                               rho_soft=rho_soft)
    2893            8 :             rho_total = rho_hard + rho_soft
    2894              :             CALL get_pw_grid_info(pw_grid=rho_elec_gspace%pw_grid, &
    2895            8 :                                   vol=volume)
    2896              :             ! rhotot pw coefficients are by default scaled by grid volume
    2897              :             ! need to undo this to get proper charge from printed cube
    2898            8 :             CALL pw_scale(rho_elec_gspace, 1.0_dp/volume)
    2899              : 
    2900            8 :             CALL pw_transfer(rho_elec_gspace, rho_elec_rspace)
    2901            8 :             rho_total_rspace = pw_integrate_function(rho_elec_rspace, isign=-1)
    2902            8 :             filename = "TOTAL_ELECTRON_DENSITY"
    2903            8 :             mpi_io = .TRUE.
    2904              :             unit_nr = e_density_section%print_key_unit_nr( &
    2905              :                       logger, &
    2906              :                       input, &
    2907              :                       e_density_section%absolute_section_key, &
    2908              :                       extension=".cube", &
    2909              :                       middle_name=TRIM(filename), &
    2910              :                       file_position=my_pos_cube, &
    2911              :                       log_filename=.FALSE., &
    2912              :                       mpi_io=mpi_io, &
    2913              :                       fout=mpi_filename, &
    2914              :                       openpmd_basename="dft-total-electron-density", &
    2915              :                       openpmd_unit_dimension=openpmd_unit_dimension_density, &
    2916              :                       openpmd_unit_si=openpmd_unit_si_density, &
    2917            8 :                       sim_time=qs_env%sim_time)
    2918            8 :             IF (output_unit > 0) THEN
    2919            4 :                IF (.NOT. mpi_io) THEN
    2920            0 :                   INQUIRE (UNIT=unit_nr, NAME=filename)
    2921              :                ELSE
    2922            4 :                   filename = mpi_filename
    2923              :                END IF
    2924              :                CALL print_density_output_message(output_unit, "The total electron density", &
    2925            4 :                                                  e_density_section, filename)
    2926              :                WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") &
    2927            4 :                   "q(max) [1/Angstrom]              :", q_max/angstrom, &
    2928            4 :                   "Soft electronic charge (G-space) :", rho_soft, &
    2929            4 :                   "Hard electronic charge (G-space) :", rho_hard, &
    2930            4 :                   "Total electronic charge (G-space):", rho_total, &
    2931            8 :                   "Total electronic charge (R-space):", rho_total_rspace
    2932              :             END IF
    2933              :             CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "TOTAL ELECTRON DENSITY", &
    2934              :                                             particles=particles, zeff=zcharge, &
    2935            8 :                               stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
    2936              :             CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
    2937            8 :                                                              e_density_section%absolute_section_key, mpi_io=mpi_io)
    2938              :             ! Print total spin density for spin-polarized systems
    2939            8 :             IF (dft_control%nspins > 1) THEN
    2940            6 :                CALL pw_zero(rho_elec_gspace)
    2941            6 :                CALL pw_zero(rho_elec_rspace)
    2942              :                CALL calculate_rhotot_elec_gspace(qs_env=qs_env, &
    2943              :                                                  auxbas_pw_pool=auxbas_pw_pool, &
    2944              :                                                  rhotot_elec_gspace=rho_elec_gspace, &
    2945              :                                                  q_max=q_max, &
    2946              :                                                  rho_hard=rho_hard, &
    2947              :                                                  rho_soft=rho_soft, &
    2948            6 :                                                  fsign=-1.0_dp)
    2949            6 :                rho_total = rho_hard + rho_soft
    2950              : 
    2951              :                ! rhotot pw coefficients are by default scaled by grid volume
    2952              :                ! need to undo this to get proper charge from printed cube
    2953            6 :                CALL pw_scale(rho_elec_gspace, 1.0_dp/volume)
    2954              : 
    2955            6 :                CALL pw_transfer(rho_elec_gspace, rho_elec_rspace)
    2956            6 :                rho_total_rspace = pw_integrate_function(rho_elec_rspace, isign=-1)
    2957            6 :                filename = "TOTAL_SPIN_DENSITY"
    2958            6 :                mpi_io = .TRUE.
    2959              :                unit_nr = e_density_section%print_key_unit_nr( &
    2960              :                          logger, &
    2961              :                          input, &
    2962              :                          e_density_section%absolute_section_key, &
    2963              :                          extension=".cube", &
    2964              :                          middle_name=TRIM(filename), &
    2965              :                          file_position=my_pos_cube, &
    2966              :                          log_filename=.FALSE., &
    2967              :                          mpi_io=mpi_io, &
    2968              :                          fout=mpi_filename, &
    2969              :                          openpmd_basename="dft-total-spin-density", &
    2970              :                          openpmd_unit_dimension=openpmd_unit_dimension_density, &
    2971              :                          openpmd_unit_si=openpmd_unit_si_density, &
    2972            6 :                          sim_time=qs_env%sim_time)
    2973            6 :                IF (output_unit > 0) THEN
    2974            3 :                   IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
    2975            0 :                      INQUIRE (UNIT=unit_nr, NAME=filename)
    2976              :                   ELSE
    2977            3 :                      filename = mpi_filename
    2978              :                   END IF
    2979              :                   CALL print_density_output_message(output_unit, "The total spin density", &
    2980            3 :                                                     e_density_section, filename)
    2981              :                   WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") &
    2982            3 :                      "q(max) [1/Angstrom]                    :", q_max/angstrom, &
    2983            3 :                      "Soft part of the spin density (G-space):", rho_soft, &
    2984            3 :                      "Hard part of the spin density (G-space):", rho_hard, &
    2985            3 :                      "Total spin density (G-space)           :", rho_total, &
    2986            6 :                      "Total spin density (R-space)           :", rho_total_rspace
    2987              :                END IF
    2988              :                CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "TOTAL SPIN DENSITY", &
    2989              :                                                particles=particles, zeff=zcharge, &
    2990            6 :                               stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
    2991              :                CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
    2992            6 :                                                                 e_density_section%absolute_section_key, mpi_io=mpi_io)
    2993              :             END IF
    2994            8 :             CALL auxbas_pw_pool%give_back_pw(rho_elec_gspace)
    2995            8 :             CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
    2996              : 
    2997          146 :          ELSE IF (print_density == "SOFT_DENSITY" .OR. .NOT. dft_control%qs_control%gapw) THEN
    2998          142 :             IF (dft_control%nspins > 1) THEN
    2999              :                CALL get_qs_env(qs_env=qs_env, &
    3000           48 :                                pw_env=pw_env)
    3001              :                CALL pw_env_get(pw_env=pw_env, &
    3002              :                                auxbas_pw_pool=auxbas_pw_pool, &
    3003           48 :                                pw_pools=pw_pools)
    3004           48 :                CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
    3005           48 :                CALL pw_copy(rho_r(1), rho_elec_rspace)
    3006           48 :                CALL pw_axpy(rho_r(2), rho_elec_rspace)
    3007           48 :                filename = "ELECTRON_DENSITY"
    3008           48 :                mpi_io = .TRUE.
    3009              :                unit_nr = e_density_section%print_key_unit_nr( &
    3010              :                          logger, &
    3011              :                          input, &
    3012              :                          e_density_section%absolute_section_key, &
    3013              :                          extension=".cube", &
    3014              :                          middle_name=TRIM(filename), &
    3015              :                          file_position=my_pos_cube, &
    3016              :                          log_filename=.FALSE., &
    3017              :                          mpi_io=mpi_io, &
    3018              :                          fout=mpi_filename, &
    3019              :                          openpmd_basename="dft-electron-density", &
    3020              :                          openpmd_unit_dimension=openpmd_unit_dimension_density, &
    3021              :                          openpmd_unit_si=openpmd_unit_si_density, &
    3022           48 :                          sim_time=qs_env%sim_time)
    3023           48 :                IF (output_unit > 0) THEN
    3024           24 :                   IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
    3025            0 :                      INQUIRE (UNIT=unit_nr, NAME=filename)
    3026              :                   ELSE
    3027           24 :                      filename = mpi_filename
    3028              :                   END IF
    3029              :                   CALL print_density_output_message(output_unit, "The sum of alpha and beta density", &
    3030           24 :                                                     e_density_section, filename)
    3031              :                END IF
    3032              :                CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "SUM OF ALPHA AND BETA DENSITY", &
    3033              :         particles=particles, zeff=zcharge, stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), &
    3034           48 :                                                mpi_io=mpi_io)
    3035              :                CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
    3036           48 :                                                                 e_density_section%absolute_section_key, mpi_io=mpi_io)
    3037           48 :                CALL pw_copy(rho_r(1), rho_elec_rspace)
    3038           48 :                CALL pw_axpy(rho_r(2), rho_elec_rspace, alpha=-1.0_dp)
    3039           48 :                filename = "SPIN_DENSITY"
    3040           48 :                mpi_io = .TRUE.
    3041              :                unit_nr = e_density_section%print_key_unit_nr( &
    3042              :                          logger, &
    3043              :                          input, &
    3044              :                          e_density_section%absolute_section_key, &
    3045              :                          extension=".cube", &
    3046              :                          middle_name=TRIM(filename), &
    3047              :                          file_position=my_pos_cube, &
    3048              :                          log_filename=.FALSE., &
    3049              :                          mpi_io=mpi_io, &
    3050              :                          fout=mpi_filename, &
    3051              :                          openpmd_basename="dft-spin-density", &
    3052              :                          openpmd_unit_dimension=openpmd_unit_dimension_density, &
    3053              :                          openpmd_unit_si=openpmd_unit_si_density, &
    3054           48 :                          sim_time=qs_env%sim_time)
    3055           48 :                IF (output_unit > 0) THEN
    3056           24 :                   IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
    3057            0 :                      INQUIRE (UNIT=unit_nr, NAME=filename)
    3058              :                   ELSE
    3059           24 :                      filename = mpi_filename
    3060              :                   END IF
    3061              :                   CALL print_density_output_message(output_unit, "The spin density", &
    3062           24 :                                                     e_density_section, filename)
    3063              :                END IF
    3064              :                CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "SPIN DENSITY", &
    3065              :                                                particles=particles, zeff=zcharge, &
    3066           48 :                               stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
    3067              :                CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
    3068           48 :                                                                 e_density_section%absolute_section_key, mpi_io=mpi_io)
    3069           48 :                CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
    3070              :             ELSE
    3071           94 :                filename = "ELECTRON_DENSITY"
    3072           94 :                mpi_io = .TRUE.
    3073              :                unit_nr = e_density_section%print_key_unit_nr( &
    3074              :                          logger, &
    3075              :                          input, &
    3076              :                          e_density_section%absolute_section_key, &
    3077              :                          extension=".cube", &
    3078              :                          middle_name=TRIM(filename), &
    3079              :                          file_position=my_pos_cube, &
    3080              :                          log_filename=.FALSE., &
    3081              :                          mpi_io=mpi_io, &
    3082              :                          fout=mpi_filename, &
    3083              :                          openpmd_basename="dft-electron-density", &
    3084              :                          openpmd_unit_dimension=openpmd_unit_dimension_density, &
    3085              :                          openpmd_unit_si=openpmd_unit_si_density, &
    3086           94 :                          sim_time=qs_env%sim_time)
    3087           94 :                IF (output_unit > 0) THEN
    3088           47 :                   IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
    3089            0 :                      INQUIRE (UNIT=unit_nr, NAME=filename)
    3090              :                   ELSE
    3091           47 :                      filename = mpi_filename
    3092              :                   END IF
    3093              :                   CALL print_density_output_message(output_unit, "The electron density", &
    3094           47 :                                                     e_density_section, filename)
    3095              :                END IF
    3096              :                CALL e_density_section%write_pw(rho_r(1), unit_nr, "ELECTRON DENSITY", &
    3097              :                                                particles=particles, zeff=zcharge, &
    3098           94 :                               stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
    3099              :                CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
    3100           94 :                                                                 e_density_section%absolute_section_key, mpi_io=mpi_io)
    3101              :             END IF ! nspins
    3102              : 
    3103            4 :          ELSE IF (dft_control%qs_control%gapw .AND. print_density == "TOTAL_HARD_APPROX") THEN
    3104            4 :             CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho0_mpole=rho0_mpole, natom=natom)
    3105            4 :             CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
    3106            4 :             CALL auxbas_pw_pool%create_pw(rho_elec_rspace)
    3107              : 
    3108            4 :             NULLIFY (my_Q0)
    3109           12 :             ALLOCATE (my_Q0(natom))
    3110           16 :             my_Q0 = 0.0_dp
    3111              : 
    3112              :             ! (eta/pi)**3: normalization for 3d gaussian of form exp(-eta*r**2)
    3113            4 :             norm_factor = SQRT((rho0_mpole%zet0_h/pi)**3)
    3114              : 
    3115              :             ! store hard part of electronic density in array
    3116           16 :             DO iat = 1, natom
    3117           34 :                my_Q0(iat) = SUM(rho0_mpole%mp_rho(iat)%Q0(1:dft_control%nspins))*norm_factor
    3118              :             END DO
    3119              :             ! multiply coeff with gaussian and put on realspace grid
    3120              :             ! coeff is the gaussian prefactor, eta the gaussian exponent
    3121            4 :             CALL calculate_rho_resp_all(rho_elec_rspace, coeff=my_Q0, natom=natom, eta=rho0_mpole%zet0_h, qs_env=qs_env)
    3122            4 :             rho_hard = pw_integrate_function(rho_elec_rspace, isign=-1)
    3123              : 
    3124            4 :             rho_soft = 0.0_dp
    3125           10 :             DO ispin = 1, dft_control%nspins
    3126            6 :                CALL pw_axpy(rho_r(ispin), rho_elec_rspace)
    3127           10 :                rho_soft = rho_soft + pw_integrate_function(rho_r(ispin), isign=-1)
    3128              :             END DO
    3129              : 
    3130            4 :             rho_total_rspace = rho_soft + rho_hard
    3131              : 
    3132            4 :             filename = "ELECTRON_DENSITY"
    3133            4 :             mpi_io = .TRUE.
    3134              :             unit_nr = e_density_section%print_key_unit_nr( &
    3135              :                       logger, &
    3136              :                       input, &
    3137              :                       e_density_section%absolute_section_key, &
    3138              :                       extension=".cube", &
    3139              :                       middle_name=TRIM(filename), &
    3140              :                       file_position=my_pos_cube, &
    3141              :                       log_filename=.FALSE., &
    3142              :                       mpi_io=mpi_io, &
    3143              :                       fout=mpi_filename, &
    3144              :                       openpmd_basename="dft-electron-density", &
    3145              :                       openpmd_unit_dimension=openpmd_unit_dimension_density, &
    3146              :                       openpmd_unit_si=openpmd_unit_si_density, &
    3147            4 :                       sim_time=qs_env%sim_time)
    3148            4 :             IF (output_unit > 0) THEN
    3149            2 :                IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
    3150            0 :                   INQUIRE (UNIT=unit_nr, NAME=filename)
    3151              :                ELSE
    3152            2 :                   filename = mpi_filename
    3153              :                END IF
    3154              :                CALL print_density_output_message(output_unit, "The electron density", &
    3155            2 :                                                  e_density_section, filename)
    3156              :                WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") &
    3157            2 :                   "Soft electronic charge (R-space) :", rho_soft, &
    3158            2 :                   "Hard electronic charge (R-space) :", rho_hard, &
    3159            4 :                   "Total electronic charge (R-space):", rho_total_rspace
    3160              :             END IF
    3161              :             CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "ELECTRON DENSITY", &
    3162              :         particles=particles, zeff=zcharge, stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), &
    3163            4 :                                             mpi_io=mpi_io)
    3164              :             CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
    3165            4 :                                                              e_density_section%absolute_section_key, mpi_io=mpi_io)
    3166              : 
    3167              :             !------------
    3168            4 :             IF (dft_control%nspins > 1) THEN
    3169            8 :             DO iat = 1, natom
    3170            8 :                my_Q0(iat) = (rho0_mpole%mp_rho(iat)%Q0(1) - rho0_mpole%mp_rho(iat)%Q0(2))*norm_factor
    3171              :             END DO
    3172            2 :             CALL pw_zero(rho_elec_rspace)
    3173            2 :             CALL calculate_rho_resp_all(rho_elec_rspace, coeff=my_Q0, natom=natom, eta=rho0_mpole%zet0_h, qs_env=qs_env)
    3174            2 :             rho_hard = pw_integrate_function(rho_elec_rspace, isign=-1)
    3175              : 
    3176            2 :             CALL pw_axpy(rho_r(1), rho_elec_rspace)
    3177            2 :             CALL pw_axpy(rho_r(2), rho_elec_rspace, alpha=-1.0_dp)
    3178              :             rho_soft = pw_integrate_function(rho_r(1), isign=-1) &
    3179            2 :                        - pw_integrate_function(rho_r(2), isign=-1)
    3180              : 
    3181            2 :             rho_total_rspace = rho_soft + rho_hard
    3182              : 
    3183            2 :             filename = "SPIN_DENSITY"
    3184            2 :             mpi_io = .TRUE.
    3185              :             unit_nr = e_density_section%print_key_unit_nr( &
    3186              :                       logger, &
    3187              :                       input, &
    3188              :                       e_density_section%absolute_section_key, &
    3189              :                       extension=".cube", &
    3190              :                       middle_name=TRIM(filename), &
    3191              :                       file_position=my_pos_cube, &
    3192              :                       log_filename=.FALSE., &
    3193              :                       mpi_io=mpi_io, &
    3194              :                       fout=mpi_filename, &
    3195              :                       openpmd_basename="dft-spin-density", &
    3196              :                       openpmd_unit_dimension=openpmd_unit_dimension_density, &
    3197              :                       openpmd_unit_si=openpmd_unit_si_density, &
    3198            2 :                       sim_time=qs_env%sim_time)
    3199            2 :             IF (output_unit > 0) THEN
    3200            1 :                IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
    3201            0 :                   INQUIRE (UNIT=unit_nr, NAME=filename)
    3202              :                ELSE
    3203            1 :                   filename = mpi_filename
    3204              :                END IF
    3205              :                CALL print_density_output_message(output_unit, "The spin density", &
    3206            1 :                                                  e_density_section, filename)
    3207              :                WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") &
    3208            1 :                   "Soft part of the spin density          :", rho_soft, &
    3209            1 :                   "Hard part of the spin density          :", rho_hard, &
    3210            2 :                   "Total spin density (R-space)           :", rho_total_rspace
    3211              :             END IF
    3212              :             CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "SPIN DENSITY", &
    3213              :                                             particles=particles, zeff=zcharge, &
    3214            2 :                               stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
    3215              :             CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
    3216            2 :                                                              e_density_section%absolute_section_key, mpi_io=mpi_io)
    3217              :             END IF ! nspins
    3218            4 :             CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
    3219            4 :             DEALLOCATE (my_Q0)
    3220              :          END IF ! print_density
    3221              :       END IF ! print key
    3222              : 
    3223              :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
    3224        13441 :                                            dft_section, "PRINT%ENERGY_WINDOWS"), cp_p_file) .AND. .NOT. do_kpoints) THEN
    3225           90 :          CALL energy_windows(qs_env)
    3226              :       END IF
    3227              : 
    3228              :       ! Print the hartree potential
    3229        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
    3230              :                                            "DFT%PRINT%V_HARTREE_CUBE"), cp_p_file)) THEN
    3231              : 
    3232              :          CALL get_qs_env(qs_env=qs_env, &
    3233              :                          pw_env=pw_env, &
    3234          114 :                          v_hartree_rspace=v_hartree_rspace)
    3235          114 :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
    3236          114 :          CALL auxbas_pw_pool%create_pw(aux_r)
    3237              : 
    3238          114 :          append_cube = section_get_lval(input, "DFT%PRINT%V_HARTREE_CUBE%APPEND")
    3239          114 :          my_pos_cube = "REWIND"
    3240          114 :          IF (append_cube) THEN
    3241            0 :             my_pos_cube = "APPEND"
    3242              :          END IF
    3243          114 :          mpi_io = .TRUE.
    3244          114 :          CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
    3245          114 :          CALL pw_env_get(pw_env)
    3246              :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%V_HARTREE_CUBE", &
    3247          114 :                                         extension=".cube", middle_name="v_hartree", file_position=my_pos_cube, mpi_io=mpi_io)
    3248          114 :          udvol = 1.0_dp/v_hartree_rspace%pw_grid%dvol
    3249              : 
    3250          114 :          CALL pw_copy(v_hartree_rspace, aux_r)
    3251          114 :          CALL pw_scale(aux_r, udvol)
    3252              : 
    3253              :          CALL cp_pw_to_cube(aux_r, unit_nr, "HARTREE POTENTIAL", particles=particles, zeff=zcharge, &
    3254              :                             stride=section_get_ivals(dft_section, "PRINT%V_HARTREE_CUBE%STRIDE"), &
    3255              :                             max_file_size_mb=section_get_rval(dft_section, "PRINT%V_HARTREE_CUBE%MAX_FILE_SIZE_MB"), &
    3256          114 :                             mpi_io=mpi_io)
    3257              :          CALL cp_print_key_finished_output(unit_nr, logger, input, &
    3258          114 :                                            "DFT%PRINT%V_HARTREE_CUBE", mpi_io=mpi_io)
    3259              : 
    3260          114 :          CALL auxbas_pw_pool%give_back_pw(aux_r)
    3261              :       END IF
    3262              : 
    3263              :       ! Print the external potential
    3264        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
    3265              :                                            "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE"), cp_p_file)) THEN
    3266           86 :          IF (dft_control%apply_external_potential) THEN
    3267            4 :             CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, vee=vee)
    3268            4 :             CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
    3269            4 :             CALL auxbas_pw_pool%create_pw(aux_r)
    3270              : 
    3271            4 :             append_cube = section_get_lval(input, "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE%APPEND")
    3272            4 :             my_pos_cube = "REWIND"
    3273            4 :             IF (append_cube) THEN
    3274            0 :                my_pos_cube = "APPEND"
    3275              :             END IF
    3276            4 :             mpi_io = .TRUE.
    3277            4 :             CALL pw_env_get(pw_env)
    3278              :             unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE", &
    3279            4 :                                            extension=".cube", middle_name="ext_pot", file_position=my_pos_cube, mpi_io=mpi_io)
    3280              : 
    3281            4 :             CALL pw_copy(vee, aux_r)
    3282              : 
    3283              :             CALL cp_pw_to_cube(aux_r, unit_nr, "EXTERNAL POTENTIAL", particles=particles, zeff=zcharge, &
    3284              :                                stride=section_get_ivals(dft_section, "PRINT%EXTERNAL_POTENTIAL_CUBE%STRIDE"), &
    3285              :                                max_file_size_mb=section_get_rval(dft_section, "PRINT%EXTERNAL_POTENTIAL_CUBE%MAX_FILE_SIZE_MB"), &
    3286            4 :                                mpi_io=mpi_io)
    3287              :             CALL cp_print_key_finished_output(unit_nr, logger, input, &
    3288            4 :                                               "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE", mpi_io=mpi_io)
    3289              : 
    3290            4 :             CALL auxbas_pw_pool%give_back_pw(aux_r)
    3291              :          END IF
    3292              :       END IF
    3293              : 
    3294              :       ! Print the Electrical Field Components
    3295        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
    3296              :                                            "DFT%PRINT%EFIELD_CUBE"), cp_p_file)) THEN
    3297              : 
    3298           82 :          CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
    3299           82 :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
    3300           82 :          CALL auxbas_pw_pool%create_pw(aux_r)
    3301           82 :          CALL auxbas_pw_pool%create_pw(aux_g)
    3302              : 
    3303           82 :          append_cube = section_get_lval(input, "DFT%PRINT%EFIELD_CUBE%APPEND")
    3304           82 :          my_pos_cube = "REWIND"
    3305           82 :          IF (append_cube) THEN
    3306            0 :             my_pos_cube = "APPEND"
    3307              :          END IF
    3308              :          CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, &
    3309           82 :                          v_hartree_rspace=v_hartree_rspace)
    3310           82 :          CALL pw_env_get(pw_env)
    3311           82 :          udvol = 1.0_dp/v_hartree_rspace%pw_grid%dvol
    3312          328 :          DO id = 1, 3
    3313          246 :             mpi_io = .TRUE.
    3314              :             unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EFIELD_CUBE", &
    3315              :                                            extension=".cube", middle_name="efield_"//cdir(id), file_position=my_pos_cube, &
    3316          246 :                                            mpi_io=mpi_io)
    3317              : 
    3318          246 :             CALL pw_transfer(v_hartree_rspace, aux_g)
    3319          246 :             nd = 0
    3320          246 :             nd(id) = 1
    3321          246 :             CALL pw_derive(aux_g, nd)
    3322          246 :             CALL pw_transfer(aux_g, aux_r)
    3323          246 :             CALL pw_scale(aux_r, udvol)
    3324              : 
    3325              :             CALL cp_pw_to_cube(aux_r, unit_nr, "ELECTRIC FIELD", particles=particles, zeff=zcharge, &
    3326              :                                stride=section_get_ivals(dft_section, "PRINT%EFIELD_CUBE%STRIDE"), &
    3327              :                                max_file_size_mb=section_get_rval(dft_section, "PRINT%EFIELD_CUBE%MAX_FILE_SIZE_MB"), &
    3328          246 :                                mpi_io=mpi_io)
    3329              :             CALL cp_print_key_finished_output(unit_nr, logger, input, &
    3330          328 :                                               "DFT%PRINT%EFIELD_CUBE", mpi_io=mpi_io)
    3331              :          END DO
    3332              : 
    3333           82 :          CALL auxbas_pw_pool%give_back_pw(aux_r)
    3334           82 :          CALL auxbas_pw_pool%give_back_pw(aux_g)
    3335              :       END IF
    3336              : 
    3337              :       ! Write cube files from the local energy
    3338        13441 :       CALL qs_scf_post_local_energy(input, logger, qs_env)
    3339              : 
    3340              :       ! Write cube files from the local stress tensor
    3341        13441 :       CALL qs_scf_post_local_stress(input, logger, qs_env)
    3342              : 
    3343              :       ! Write cube files from the implicit Poisson solver
    3344        13441 :       CALL qs_scf_post_ps_implicit(input, logger, qs_env)
    3345              : 
    3346              :       ! post SCF finite-volume Kubo transport
    3347        13441 :       CALL qs_scf_post_kubo_transport(qs_env)
    3348              : 
    3349              :       ! post SCF Transport
    3350        13441 :       CALL qs_scf_post_transport(qs_env)
    3351              : 
    3352        13441 :       CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
    3353              :       ! Write the density matrices
    3354        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
    3355              :                                            "DFT%PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN
    3356              :          iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/DENSITY", &
    3357            4 :                                    extension=".Log")
    3358            4 :          CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
    3359            4 :          CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
    3360            4 :          after = MIN(MAX(after, 1), 16)
    3361            8 :          DO ispin = 1, dft_control%nspins
    3362           12 :             DO img = 1, dft_control%nimages
    3363              :                CALL cp_dbcsr_write_sparse_matrix(rho_ao(ispin, img)%matrix, 4, after, qs_env, &
    3364            8 :                                                  para_env, output_unit=iw, omit_headers=omit_headers)
    3365              :             END DO
    3366              :          END DO
    3367              :          CALL cp_print_key_finished_output(iw, logger, input, &
    3368            4 :                                            "DFT%PRINT%AO_MATRICES/DENSITY")
    3369              :       END IF
    3370              : 
    3371              :       ! Write the Kohn-Sham matrices
    3372              :       write_ks = BTEST(cp_print_key_should_output(logger%iter_info, input, &
    3373        13441 :                                                   "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)
    3374              :       write_xc = BTEST(cp_print_key_should_output(logger%iter_info, input, &
    3375        13441 :                                                   "DFT%PRINT%AO_MATRICES/MATRIX_VXC"), cp_p_file)
    3376              :       ! we need to update stuff before writing, potentially computing the matrix_vxc
    3377        13441 :       IF (write_ks .OR. write_xc) THEN
    3378            4 :          IF (write_xc) qs_env%requires_matrix_vxc = .TRUE.
    3379            4 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
    3380              :          CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., &
    3381            4 :                                   just_energy=.FALSE.)
    3382            4 :          IF (write_xc) qs_env%requires_matrix_vxc = .FALSE.
    3383              :       END IF
    3384              : 
    3385              :       ! Write the Kohn-Sham matrices
    3386        13441 :       IF (write_ks) THEN
    3387              :          iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", &
    3388            4 :                                    extension=".Log")
    3389            4 :          CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=ks_rmpv)
    3390            4 :          CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
    3391            4 :          after = MIN(MAX(after, 1), 16)
    3392            8 :          DO ispin = 1, dft_control%nspins
    3393           12 :             DO img = 1, dft_control%nimages
    3394              :                CALL cp_dbcsr_write_sparse_matrix(ks_rmpv(ispin, img)%matrix, 4, after, qs_env, &
    3395            8 :                                                  para_env, output_unit=iw, omit_headers=omit_headers)
    3396              :             END DO
    3397              :          END DO
    3398              :          CALL cp_print_key_finished_output(iw, logger, input, &
    3399            4 :                                            "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
    3400              :       END IF
    3401              : 
    3402              :       ! write csr matrices
    3403              :       ! matrices in terms of the PAO basis will be taken care of in pao_post_scf.
    3404        13441 :       IF (.NOT. dft_control%qs_control%pao) THEN
    3405        12929 :          CALL write_ks_matrix_csr(qs_env, input)
    3406        12929 :          CALL write_s_matrix_csr(qs_env, input)
    3407        12929 :          CALL write_hcore_matrix_csr(qs_env, input)
    3408        12929 :          CALL write_p_matrix_csr(qs_env, input)
    3409              :       END IF
    3410              : 
    3411              :       ! write adjacency matrix
    3412        13441 :       CALL write_adjacency_matrix(qs_env, input)
    3413              : 
    3414              :       ! Write the xc matrix
    3415        13441 :       IF (write_xc) THEN
    3416            0 :          CALL get_qs_env(qs_env=qs_env, matrix_vxc_kp=matrix_vxc)
    3417            0 :          CPASSERT(ASSOCIATED(matrix_vxc))
    3418              :          iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/MATRIX_VXC", &
    3419            0 :                                    extension=".Log")
    3420            0 :          CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
    3421            0 :          after = MIN(MAX(after, 1), 16)
    3422            0 :          DO ispin = 1, dft_control%nspins
    3423            0 :             DO img = 1, dft_control%nimages
    3424              :                CALL cp_dbcsr_write_sparse_matrix(matrix_vxc(ispin, img)%matrix, 4, after, qs_env, &
    3425            0 :                                                  para_env, output_unit=iw, omit_headers=omit_headers)
    3426              :             END DO
    3427              :          END DO
    3428              :          CALL cp_print_key_finished_output(iw, logger, input, &
    3429            0 :                                            "DFT%PRINT%AO_MATRICES/MATRIX_VXC")
    3430              :       END IF
    3431              : 
    3432              :       ! Write the [H,r] commutator matrices
    3433        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
    3434              :                                            "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR"), cp_p_file)) THEN
    3435              :          iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR", &
    3436            2 :                                    extension=".Log")
    3437            2 :          CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
    3438            2 :          NULLIFY (matrix_hr)
    3439            2 :          CALL build_com_hr_matrix(qs_env, matrix_hr)
    3440            8 :          DO img = 1, 3
    3441            8 :             checksum_hr(img) = dbcsr_checksum(matrix_hr(img)%matrix)
    3442              :          END DO
    3443            2 :          IF (output_unit > 0) THEN
    3444            1 :             WRITE (output_unit, '(T2,A,E23.16)') 'COMMUTATOR_HR| CheckSum X =', checksum_hr(1)
    3445            1 :             WRITE (output_unit, '(T2,A,E23.16)') 'COMMUTATOR_HR| CheckSum Y =', checksum_hr(2)
    3446            1 :             WRITE (output_unit, '(T2,A,E23.16)') 'COMMUTATOR_HR| CheckSum Z =', checksum_hr(3)
    3447              :          END IF
    3448            2 :          after = MIN(MAX(after, 1), 16)
    3449            8 :          DO img = 1, 3
    3450              :             CALL cp_dbcsr_write_sparse_matrix(matrix_hr(img)%matrix, 4, after, qs_env, &
    3451            8 :                                               para_env, output_unit=iw, omit_headers=omit_headers)
    3452              :          END DO
    3453            2 :          CALL dbcsr_deallocate_matrix_set(matrix_hr)
    3454              :          CALL cp_print_key_finished_output(iw, logger, input, &
    3455            2 :                                            "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR")
    3456              :       END IF
    3457              : 
    3458              :       ! Compute the Mulliken charges
    3459        13441 :       print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MULLIKEN")
    3460        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    3461         5346 :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MULLIKEN", extension=".mulliken", log_filename=.FALSE.)
    3462         5346 :          print_level = 1
    3463         5346 :          CALL section_vals_val_get(print_key, "PRINT_GOP", l_val=print_it)
    3464         5346 :          IF (print_it) print_level = 2
    3465         5346 :          CALL section_vals_val_get(print_key, "PRINT_ALL", l_val=print_it)
    3466         5346 :          IF (print_it) print_level = 3
    3467         5346 :          CALL mulliken_population_analysis(qs_env, unit_nr, print_level)
    3468         5346 :          CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MULLIKEN")
    3469              :       END IF
    3470              : 
    3471              :       ! Compute the Hirshfeld charges
    3472        13441 :       print_key => section_vals_get_subs_vals(input, "DFT%PRINT%HIRSHFELD")
    3473        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    3474              :          ! we check if real space density is available
    3475         5428 :          NULLIFY (rho)
    3476         5428 :          CALL get_qs_env(qs_env=qs_env, rho=rho)
    3477         5428 :          CALL qs_rho_get(rho, rho_r_valid=rho_r_valid)
    3478         5428 :          IF (rho_r_valid) THEN
    3479         5354 :             unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%HIRSHFELD", extension=".hirshfeld", log_filename=.FALSE.)
    3480         5354 :             CALL hirshfeld_charges(qs_env, print_key, unit_nr)
    3481         5354 :             CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%HIRSHFELD")
    3482              :          END IF
    3483              :       END IF
    3484              : 
    3485              :       ! Compute EEQ charges
    3486        13441 :       print_key => section_vals_get_subs_vals(input, "DFT%PRINT%EEQ_CHARGES")
    3487        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    3488           30 :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EEQ_CHARGES", extension=".eeq", log_filename=.FALSE.)
    3489           30 :          print_level = 1
    3490           30 :          CALL eeq_print(qs_env, unit_nr, print_level, ext=.FALSE.)
    3491           30 :          CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MULLIKEN")
    3492              :       END IF
    3493              : 
    3494              :       ! Do a Voronoi Integration or write a compressed BQB File
    3495        13441 :       print_key_voro => section_vals_get_subs_vals(input, "DFT%PRINT%VORONOI")
    3496        13441 :       print_key_bqb => section_vals_get_subs_vals(input, "DFT%PRINT%E_DENSITY_BQB")
    3497        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key_voro), cp_p_file)) THEN
    3498           24 :          should_print_voro = 1
    3499              :       ELSE
    3500        13417 :          should_print_voro = 0
    3501              :       END IF
    3502        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key_bqb), cp_p_file)) THEN
    3503            2 :          should_print_bqb = 1
    3504              :       ELSE
    3505        13439 :          should_print_bqb = 0
    3506              :       END IF
    3507        13441 :       IF ((should_print_voro /= 0) .OR. (should_print_bqb /= 0)) THEN
    3508              : 
    3509              :          ! we check if real space density is available
    3510           26 :          NULLIFY (rho)
    3511           26 :          CALL get_qs_env(qs_env=qs_env, rho=rho)
    3512           26 :          CALL qs_rho_get(rho, rho_r_valid=rho_r_valid)
    3513           26 :          IF (rho_r_valid) THEN
    3514              : 
    3515           26 :             IF (dft_control%nspins > 1) THEN
    3516              :                CALL get_qs_env(qs_env=qs_env, &
    3517            0 :                                pw_env=pw_env)
    3518              :                CALL pw_env_get(pw_env=pw_env, &
    3519              :                                auxbas_pw_pool=auxbas_pw_pool, &
    3520            0 :                                pw_pools=pw_pools)
    3521            0 :                NULLIFY (mb_rho)
    3522            0 :                ALLOCATE (mb_rho)
    3523            0 :                CALL auxbas_pw_pool%create_pw(pw=mb_rho)
    3524            0 :                CALL pw_copy(rho_r(1), mb_rho)
    3525            0 :                CALL pw_axpy(rho_r(2), mb_rho)
    3526              :                !CALL voronoi_analysis(qs_env, rho_elec_rspace, print_key, unit_nr)
    3527              :             ELSE
    3528           26 :                mb_rho => rho_r(1)
    3529              :                !CALL voronoi_analysis( qs_env, rho_r(1), print_key, unit_nr )
    3530              :             END IF ! nspins
    3531              : 
    3532           26 :             IF (should_print_voro /= 0) THEN
    3533           24 :                CALL section_vals_val_get(print_key_voro, "OUTPUT_TEXT", l_val=voro_print_txt)
    3534           24 :                IF (voro_print_txt) THEN
    3535           24 :                   append_voro = section_get_lval(input, "DFT%PRINT%VORONOI%APPEND")
    3536           24 :                   my_pos_voro = "REWIND"
    3537           24 :                   IF (append_voro) THEN
    3538            0 :                      my_pos_voro = "APPEND"
    3539              :                   END IF
    3540              :                   unit_nr_voro = cp_print_key_unit_nr(logger, input, "DFT%PRINT%VORONOI", extension=".voronoi", &
    3541           24 :                                                       file_position=my_pos_voro, log_filename=.FALSE.)
    3542              :                ELSE
    3543            0 :                   unit_nr_voro = 0
    3544              :                END IF
    3545              :             ELSE
    3546            2 :                unit_nr_voro = 0
    3547              :             END IF
    3548              : 
    3549              :             CALL entry_voronoi_or_bqb(should_print_voro, should_print_bqb, print_key_voro, print_key_bqb, &
    3550           26 :                                       unit_nr_voro, qs_env, mb_rho)
    3551              : 
    3552           26 :             IF (dft_control%nspins > 1) THEN
    3553            0 :                CALL auxbas_pw_pool%give_back_pw(mb_rho)
    3554            0 :                DEALLOCATE (mb_rho)
    3555              :             END IF
    3556              : 
    3557           26 :             IF (unit_nr_voro > 0) THEN
    3558           12 :                CALL cp_print_key_finished_output(unit_nr_voro, logger, input, "DFT%PRINT%VORONOI")
    3559              :             END IF
    3560              : 
    3561              :          END IF
    3562              :       END IF
    3563              : 
    3564              :       ! MAO analysis
    3565        13441 :       print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MAO_ANALYSIS")
    3566        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    3567           38 :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MAO_ANALYSIS", extension=".mao", log_filename=.FALSE.)
    3568           38 :          CALL mao_analysis(qs_env, print_key, unit_nr)
    3569           38 :          CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MAO_ANALYSIS")
    3570              :       END IF
    3571              : 
    3572              :       ! MINBAS analysis
    3573        13441 :       print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MINBAS_ANALYSIS")
    3574        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    3575           28 :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MINBAS_ANALYSIS", extension=".mao", log_filename=.FALSE.)
    3576           28 :          CALL minbas_analysis(qs_env, print_key, unit_nr)
    3577           28 :          CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MINBAS_ANALYSIS")
    3578              :       END IF
    3579              : 
    3580              :       ! IAO analysis
    3581        13441 :       print_key => section_vals_get_subs_vals(input, "DFT%PRINT%IAO_ANALYSIS")
    3582        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    3583           34 :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IAO_ANALYSIS", extension=".iao", log_filename=.FALSE.)
    3584           34 :          CALL iao_read_input(iao_env, print_key, cell)
    3585           34 :          IF (particle_set(1)%fragment_index /= 0) iao_env%do_fragments = .TRUE.
    3586           34 :          IF (iao_env%do_iao) THEN
    3587            6 :             CALL iao_wfn_analysis(qs_env, iao_env, unit_nr)
    3588              :          END IF
    3589           34 :          CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%IAO_ANALYSIS")
    3590              :       END IF
    3591              : 
    3592              :       ! Energy Decomposition Analysis
    3593        13441 :       print_key => section_vals_get_subs_vals(input, "DFT%PRINT%ENERGY_DECOMPOSITION_ANALYSIS")
    3594        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
    3595              :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%ENERGY_DECOMPOSITION_ANALYSIS", &
    3596           58 :                                         extension=".mao", log_filename=.FALSE.)
    3597           58 :          CALL edmf_analysis(qs_env, print_key, unit_nr)
    3598           58 :          CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%ENERGY_DECOMPOSITION_ANALYSIS")
    3599              :       END IF
    3600              : 
    3601              :       ! Print the density in the RI-HFX basis
    3602        13441 :       hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
    3603        13441 :       CALL section_vals_get(hfx_section, explicit=do_hfx)
    3604        13441 :       CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
    3605        13441 :       IF (do_hfx) THEN
    3606         5238 :          DO i = 1, n_rep_hf
    3607         5238 :             IF (qs_env%x_data(i, 1)%do_hfx_ri) CALL print_ri_hfx(qs_env%x_data(i, 1)%ri_data, qs_env)
    3608              :          END DO
    3609              :       END IF
    3610              : 
    3611        13441 :       DEALLOCATE (zcharge)
    3612              : 
    3613        13441 :       CALL timestop(handle)
    3614              : 
    3615        26882 :    END SUBROUTINE write_mo_free_results
    3616              : 
    3617              : ! **************************************************************************************************
    3618              : !> \brief Calculates Hirshfeld charges
    3619              : !> \param qs_env the qs_env where to calculate the charges
    3620              : !> \param input_section the input section for Hirshfeld charges
    3621              : !> \param unit_nr the output unit number
    3622              : ! **************************************************************************************************
    3623         5354 :    SUBROUTINE hirshfeld_charges(qs_env, input_section, unit_nr)
    3624              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    3625              :       TYPE(section_vals_type), POINTER                   :: input_section
    3626              :       INTEGER, INTENT(IN)                                :: unit_nr
    3627              : 
    3628              :       INTEGER                                            :: i, iat, ikind, natom, nkind, nspin, &
    3629              :                                                             radius_type, refc, shapef
    3630         5354 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list
    3631              :       LOGICAL                                            :: do_radius, do_sc, paw_atom
    3632              :       REAL(KIND=dp)                                      :: zeff
    3633         5354 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: radii
    3634         5354 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: charges
    3635         5354 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    3636              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    3637         5354 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p, matrix_s
    3638              :       TYPE(dft_control_type), POINTER                    :: dft_control
    3639              :       TYPE(hirshfeld_type), POINTER                      :: hirshfeld_env
    3640              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    3641         5354 :       TYPE(mpole_rho_atom), DIMENSION(:), POINTER        :: mp_rho
    3642         5354 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    3643         5354 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    3644              :       TYPE(qs_rho_type), POINTER                         :: rho
    3645              :       TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
    3646              : 
    3647         5354 :       NULLIFY (hirshfeld_env)
    3648         5354 :       NULLIFY (radii)
    3649         5354 :       CALL create_hirshfeld_type(hirshfeld_env)
    3650              :       !
    3651         5354 :       CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
    3652        16062 :       ALLOCATE (hirshfeld_env%charges(natom))
    3653              :       ! input options
    3654         5354 :       CALL section_vals_val_get(input_section, "SELF_CONSISTENT", l_val=do_sc)
    3655         5354 :       CALL section_vals_val_get(input_section, "USER_RADIUS", l_val=do_radius)
    3656         5354 :       CALL section_vals_val_get(input_section, "SHAPE_FUNCTION", i_val=shapef)
    3657         5354 :       CALL section_vals_val_get(input_section, "REFERENCE_CHARGE", i_val=refc)
    3658         5354 :       IF (do_radius) THEN
    3659            0 :          radius_type = radius_user
    3660            0 :          CALL section_vals_val_get(input_section, "ATOMIC_RADII", r_vals=radii)
    3661            0 :          IF (.NOT. SIZE(radii) == nkind) THEN
    3662              :             CALL cp_abort(__LOCATION__, &
    3663              :                           "Length of keyword HIRSHFELD\ATOMIC_RADII does not "// &
    3664            0 :                           "match number of atomic kinds in the input coordinate file.")
    3665              :          END IF
    3666              :       ELSE
    3667         5354 :          radius_type = radius_covalent
    3668              :       END IF
    3669              :       CALL set_hirshfeld_info(hirshfeld_env, shape_function_type=shapef, &
    3670              :                               iterative=do_sc, ref_charge=refc, &
    3671         5354 :                               radius_type=radius_type)
    3672              :       ! shape function
    3673         5354 :       CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
    3674              :       CALL create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, &
    3675         5354 :                                  radii_list=radii)
    3676              :       ! reference charges
    3677         5354 :       CALL get_qs_env(qs_env, rho=rho)
    3678         5354 :       CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
    3679         5354 :       nspin = SIZE(matrix_p, 1)
    3680        21416 :       ALLOCATE (charges(natom, nspin))
    3681         5342 :       SELECT CASE (refc)
    3682              :       CASE (ref_charge_atomic)
    3683        14546 :          DO ikind = 1, nkind
    3684         9204 :             CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
    3685         9204 :             atomic_kind => atomic_kind_set(ikind)
    3686         9204 :             CALL get_atomic_kind(atomic_kind, atom_list=atom_list)
    3687        46060 :             DO iat = 1, SIZE(atom_list)
    3688        22310 :                i = atom_list(iat)
    3689        31514 :                hirshfeld_env%charges(i) = zeff
    3690              :             END DO
    3691              :          END DO
    3692              :       CASE (ref_charge_mulliken)
    3693           12 :          CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env)
    3694           12 :          CALL mulliken_charges(matrix_p, matrix_s, para_env, charges)
    3695           48 :          DO iat = 1, natom
    3696          108 :             hirshfeld_env%charges(iat) = SUM(charges(iat, :))
    3697              :          END DO
    3698              :       CASE DEFAULT
    3699         5354 :          CPABORT("Unknown type of reference charge for Hirshfeld partitioning.")
    3700              :       END SELECT
    3701              :       !
    3702        36976 :       charges = 0.0_dp
    3703         5354 :       IF (hirshfeld_env%iterative) THEN
    3704              :          ! Hirshfeld-I charges
    3705           22 :          CALL comp_hirshfeld_i_charges(qs_env, hirshfeld_env, charges, unit_nr)
    3706              :       ELSE
    3707              :          ! Hirshfeld charges
    3708         5332 :          CALL comp_hirshfeld_charges(qs_env, hirshfeld_env, charges)
    3709              :       END IF
    3710         5354 :       CALL get_qs_env(qs_env, particle_set=particle_set, dft_control=dft_control)
    3711         5354 :       IF (dft_control%qs_control%gapw) THEN
    3712              :          ! GAPW: add core charges (rho_hard - rho_soft)
    3713          964 :          CALL get_qs_env(qs_env, rho0_mpole=rho0_mpole)
    3714          964 :          CALL get_rho0_mpole(rho0_mpole, mp_rho=mp_rho)
    3715         4012 :          DO iat = 1, natom
    3716         3048 :             atomic_kind => particle_set(iat)%atomic_kind
    3717         3048 :             CALL get_atomic_kind(atomic_kind, kind_number=ikind)
    3718         3048 :             CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
    3719         4012 :             IF (paw_atom) THEN
    3720         5684 :                charges(iat, 1:nspin) = charges(iat, 1:nspin) + mp_rho(iat)%q0(1:nspin)
    3721              :             END IF
    3722              :          END DO
    3723              :       END IF
    3724              :       !
    3725         5354 :       IF (unit_nr > 0) THEN
    3726              :          CALL write_hirshfeld_charges(charges, hirshfeld_env, particle_set, &
    3727         2692 :                                       qs_kind_set, unit_nr)
    3728              :       END IF
    3729              :       ! Save the charges to the results under the tag [HIRSHFELD-CHARGES]
    3730         5354 :       CALL save_hirshfeld_charges(charges, particle_set, qs_kind_set, qs_env)
    3731              :       !
    3732         5354 :       CALL release_hirshfeld_type(hirshfeld_env)
    3733         5354 :       DEALLOCATE (charges)
    3734              : 
    3735        10708 :    END SUBROUTINE hirshfeld_charges
    3736              : 
    3737              : ! **************************************************************************************************
    3738              : !> \brief ...
    3739              : !> \param ca ...
    3740              : !> \param a ...
    3741              : !> \param cb ...
    3742              : !> \param b ...
    3743              : !> \param l ...
    3744              : ! **************************************************************************************************
    3745            4 :    SUBROUTINE project_function_a(ca, a, cb, b, l)
    3746              :       ! project function cb on ca
    3747              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: ca
    3748              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: a, cb, b
    3749              :       INTEGER, INTENT(IN)                                :: l
    3750              : 
    3751              :       INTEGER                                            :: info, n
    3752            4 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: ipiv
    3753            4 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: smat, tmat, v
    3754              : 
    3755            4 :       n = SIZE(ca)
    3756           40 :       ALLOCATE (smat(n, n), tmat(n, n), v(n, 1), ipiv(n))
    3757              : 
    3758            4 :       CALL sg_overlap(smat, l, a, a)
    3759            4 :       CALL sg_overlap(tmat, l, a, b)
    3760         1252 :       v(:, 1) = MATMUL(tmat, cb)
    3761            4 :       CALL dgesv(n, 1, smat, n, ipiv, v, n, info)
    3762            4 :       CPASSERT(info == 0)
    3763           52 :       ca(:) = v(:, 1)
    3764              : 
    3765            4 :       DEALLOCATE (smat, tmat, v, ipiv)
    3766              : 
    3767            4 :    END SUBROUTINE project_function_a
    3768              : 
    3769              : ! **************************************************************************************************
    3770              : !> \brief ...
    3771              : !> \param ca ...
    3772              : !> \param a ...
    3773              : !> \param bfun ...
    3774              : !> \param grid_atom ...
    3775              : !> \param l ...
    3776              : ! **************************************************************************************************
    3777           36 :    SUBROUTINE project_function_b(ca, a, bfun, grid_atom, l)
    3778              :       ! project function f on ca
    3779              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: ca
    3780              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: a, bfun
    3781              :       TYPE(grid_atom_type), POINTER                      :: grid_atom
    3782              :       INTEGER, INTENT(IN)                                :: l
    3783              : 
    3784              :       INTEGER                                            :: i, info, n, nr
    3785           36 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: ipiv
    3786           36 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: afun
    3787           36 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: smat, v
    3788              : 
    3789           36 :       n = SIZE(ca)
    3790           36 :       nr = grid_atom%nr
    3791          360 :       ALLOCATE (smat(n, n), v(n, 1), ipiv(n), afun(nr))
    3792              : 
    3793           36 :       CALL sg_overlap(smat, l, a, a)
    3794          468 :       DO i = 1, n
    3795        22032 :          afun(:) = grid_atom%rad(:)**l*EXP(-a(i)*grid_atom%rad2(:))
    3796        22068 :          v(i, 1) = SUM(afun(:)*bfun(:)*grid_atom%wr(:))
    3797              :       END DO
    3798           36 :       CALL dgesv(n, 1, smat, n, ipiv, v, n, info)
    3799           36 :       CPASSERT(info == 0)
    3800          468 :       ca(:) = v(:, 1)
    3801              : 
    3802           36 :       DEALLOCATE (smat, v, ipiv, afun)
    3803              : 
    3804           36 :    END SUBROUTINE project_function_b
    3805              : 
    3806              : ! **************************************************************************************************
    3807              : !> \brief Performs printing of cube files from local energy
    3808              : !> \param input input
    3809              : !> \param logger the logger
    3810              : !> \param qs_env the qs_env in which the qs_env lives
    3811              : !> \par History
    3812              : !>      07.2019 created
    3813              : !> \author JGH
    3814              : ! **************************************************************************************************
    3815        13441 :    SUBROUTINE qs_scf_post_local_energy(input, logger, qs_env)
    3816              :       TYPE(section_vals_type), POINTER                   :: input
    3817              :       TYPE(cp_logger_type), POINTER                      :: logger
    3818              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    3819              : 
    3820              :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_local_energy'
    3821              : 
    3822              :       CHARACTER(LEN=default_path_length)                 :: filename, my_pos_cube
    3823              :       INTEGER                                            :: handle, io_unit, natom, unit_nr
    3824              :       LOGICAL                                            :: append_cube, gapw, gapw_xc, mpi_io
    3825        13441 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: zcharge
    3826              :       TYPE(dft_control_type), POINTER                    :: dft_control
    3827              :       TYPE(particle_list_type), POINTER                  :: particles
    3828              :       TYPE(pw_env_type), POINTER                         :: pw_env
    3829              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    3830              :       TYPE(pw_r3d_rs_type)                               :: eden
    3831              :       TYPE(qs_subsys_type), POINTER                      :: subsys
    3832              :       TYPE(section_vals_type), POINTER                   :: dft_section
    3833              : 
    3834        13441 :       CALL timeset(routineN, handle)
    3835        13441 :       io_unit = cp_logger_get_default_io_unit(logger)
    3836        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
    3837              :                                            "DFT%PRINT%LOCAL_ENERGY_CUBE"), cp_p_file)) THEN
    3838           34 :          dft_section => section_vals_get_subs_vals(input, "DFT")
    3839           34 :          CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, natom=natom)
    3840           34 :          gapw = dft_control%qs_control%gapw
    3841           34 :          gapw_xc = dft_control%qs_control%gapw_xc
    3842           34 :          CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
    3843           34 :          CALL qs_subsys_get(subsys, particles=particles)
    3844           34 :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
    3845           34 :          CALL auxbas_pw_pool%create_pw(eden)
    3846              :          !
    3847           34 :          CALL qs_local_energy(qs_env, eden)
    3848           34 :          CALL get_effective_core_charges(qs_env, zcharge)
    3849              :          !
    3850           34 :          append_cube = section_get_lval(input, "DFT%PRINT%LOCAL_ENERGY_CUBE%APPEND")
    3851           34 :          IF (append_cube) THEN
    3852            0 :             my_pos_cube = "APPEND"
    3853              :          ELSE
    3854           34 :             my_pos_cube = "REWIND"
    3855              :          END IF
    3856           34 :          mpi_io = .TRUE.
    3857              :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%LOCAL_ENERGY_CUBE", &
    3858              :                                         extension=".cube", middle_name="local_energy", &
    3859           34 :                                         file_position=my_pos_cube, mpi_io=mpi_io)
    3860              :          CALL cp_pw_to_cube(eden, unit_nr, "LOCAL ENERGY", particles=particles, zeff=zcharge, &
    3861              :                             stride=section_get_ivals(dft_section, "PRINT%LOCAL_ENERGY_CUBE%STRIDE"), &
    3862              :                             max_file_size_mb=section_get_rval(dft_section, "PRINT%LOCAL_ENERGY_CUBE%MAX_FILE_SIZE_MB"), &
    3863           34 :                             mpi_io=mpi_io)
    3864           34 :          IF (io_unit > 0) THEN
    3865           17 :             INQUIRE (UNIT=unit_nr, NAME=filename)
    3866           17 :             IF (gapw .OR. gapw_xc) THEN
    3867              :                WRITE (UNIT=io_unit, FMT="(/,T3,A,A)") &
    3868            1 :                   "The soft part of the local energy is written to the file: ", TRIM(ADJUSTL(filename))
    3869              :             ELSE
    3870              :                WRITE (UNIT=io_unit, FMT="(/,T3,A,A)") &
    3871           16 :                   "The local energy is written to the file: ", TRIM(ADJUSTL(filename))
    3872              :             END IF
    3873              :          END IF
    3874              :          CALL cp_print_key_finished_output(unit_nr, logger, input, &
    3875           34 :                                            "DFT%PRINT%LOCAL_ENERGY_CUBE", mpi_io=mpi_io)
    3876              :          !
    3877           34 :          CALL auxbas_pw_pool%give_back_pw(eden)
    3878           34 :          DEALLOCATE (zcharge)
    3879              :       END IF
    3880        13441 :       CALL timestop(handle)
    3881              : 
    3882        26882 :    END SUBROUTINE qs_scf_post_local_energy
    3883              : 
    3884              : ! **************************************************************************************************
    3885              : !> \brief Performs printing of cube files from local energy
    3886              : !> \param input input
    3887              : !> \param logger the logger
    3888              : !> \param qs_env the qs_env in which the qs_env lives
    3889              : !> \par History
    3890              : !>      07.2019 created
    3891              : !> \author JGH
    3892              : ! **************************************************************************************************
    3893        13441 :    SUBROUTINE qs_scf_post_local_stress(input, logger, qs_env)
    3894              :       TYPE(section_vals_type), POINTER                   :: input
    3895              :       TYPE(cp_logger_type), POINTER                      :: logger
    3896              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    3897              : 
    3898              :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_local_stress'
    3899              : 
    3900              :       CHARACTER(LEN=default_path_length)                 :: filename, my_pos_cube
    3901              :       INTEGER                                            :: handle, io_unit, natom, unit_nr
    3902              :       LOGICAL                                            :: append_cube, gapw, gapw_xc, mpi_io
    3903              :       REAL(KIND=dp)                                      :: beta
    3904        13441 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: zcharge
    3905              :       TYPE(dft_control_type), POINTER                    :: dft_control
    3906              :       TYPE(particle_list_type), POINTER                  :: particles
    3907              :       TYPE(pw_env_type), POINTER                         :: pw_env
    3908              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    3909              :       TYPE(pw_r3d_rs_type)                               :: stress
    3910              :       TYPE(qs_subsys_type), POINTER                      :: subsys
    3911              :       TYPE(section_vals_type), POINTER                   :: dft_section
    3912              : 
    3913        13441 :       CALL timeset(routineN, handle)
    3914        13441 :       io_unit = cp_logger_get_default_io_unit(logger)
    3915        13441 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
    3916              :                                            "DFT%PRINT%LOCAL_STRESS_CUBE"), cp_p_file)) THEN
    3917              :          CALL cp_warn(__LOCATION__, &
    3918           30 :                       "LOCAL_STRESS_CUBE uses the existing experimental local stress implementation")
    3919           30 :          dft_section => section_vals_get_subs_vals(input, "DFT")
    3920           30 :          CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, natom=natom)
    3921           30 :          gapw = dft_control%qs_control%gapw
    3922           30 :          gapw_xc = dft_control%qs_control%gapw_xc
    3923           30 :          CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
    3924           30 :          CALL qs_subsys_get(subsys, particles=particles)
    3925           30 :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
    3926           30 :          CALL auxbas_pw_pool%create_pw(stress)
    3927              :          !
    3928              :          ! use beta=0: kinetic energy density in symmetric form
    3929           30 :          beta = 0.0_dp
    3930           30 :          CALL qs_local_stress(qs_env, beta=beta)
    3931           30 :          CALL get_effective_core_charges(qs_env, zcharge)
    3932              :          !
    3933           30 :          append_cube = section_get_lval(input, "DFT%PRINT%LOCAL_STRESS_CUBE%APPEND")
    3934           30 :          IF (append_cube) THEN
    3935            0 :             my_pos_cube = "APPEND"
    3936              :          ELSE
    3937           30 :             my_pos_cube = "REWIND"
    3938              :          END IF
    3939           30 :          mpi_io = .TRUE.
    3940              :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%LOCAL_STRESS_CUBE", &
    3941              :                                         extension=".cube", middle_name="local_stress", &
    3942           30 :                                         file_position=my_pos_cube, mpi_io=mpi_io)
    3943              :          CALL cp_pw_to_cube(stress, unit_nr, "LOCAL STRESS", particles=particles, zeff=zcharge, &
    3944              :                             stride=section_get_ivals(dft_section, "PRINT%LOCAL_STRESS_CUBE%STRIDE"), &
    3945              :                             max_file_size_mb=section_get_rval(dft_section, "PRINT%LOCAL_STRESS_CUBE%MAX_FILE_SIZE_MB"), &
    3946           30 :                             mpi_io=mpi_io)
    3947           30 :          IF (io_unit > 0) THEN
    3948           15 :             INQUIRE (UNIT=unit_nr, NAME=filename)
    3949           15 :             WRITE (UNIT=io_unit, FMT="(/,T3,A)") "Write 1/3*Tr(sigma) to cube file"
    3950           15 :             IF (gapw .OR. gapw_xc) THEN
    3951              :                WRITE (UNIT=io_unit, FMT="(T3,A,A)") &
    3952            0 :                   "The soft part of the local stress is written to the file: ", TRIM(ADJUSTL(filename))
    3953              :             ELSE
    3954              :                WRITE (UNIT=io_unit, FMT="(T3,A,A)") &
    3955           15 :                   "The local stress is written to the file: ", TRIM(ADJUSTL(filename))
    3956              :             END IF
    3957              :          END IF
    3958              :          CALL cp_print_key_finished_output(unit_nr, logger, input, &
    3959           30 :                                            "DFT%PRINT%LOCAL_STRESS_CUBE", mpi_io=mpi_io)
    3960              :          !
    3961           30 :          CALL auxbas_pw_pool%give_back_pw(stress)
    3962           30 :          DEALLOCATE (zcharge)
    3963              :       END IF
    3964              : 
    3965        13441 :       CALL timestop(handle)
    3966              : 
    3967        26882 :    END SUBROUTINE qs_scf_post_local_stress
    3968              : 
    3969              : ! **************************************************************************************************
    3970              : !> \brief Performs printing of cube files related to the implicit Poisson solver
    3971              : !> \param input input
    3972              : !> \param logger the logger
    3973              : !> \param qs_env the qs_env in which the qs_env lives
    3974              : !> \par History
    3975              : !>      03.2016 refactored from write_mo_free_results [Hossein Bani-Hashemian]
    3976              : !> \author Mohammad Hossein Bani-Hashemian
    3977              : ! **************************************************************************************************
    3978        13441 :    SUBROUTINE qs_scf_post_ps_implicit(input, logger, qs_env)
    3979              :       TYPE(section_vals_type), POINTER                   :: input
    3980              :       TYPE(cp_logger_type), POINTER                      :: logger
    3981              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    3982              : 
    3983              :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_ps_implicit'
    3984              : 
    3985              :       CHARACTER(LEN=default_path_length)                 :: filename, my_pos_cube
    3986              :       INTEGER                                            :: boundary_condition, handle, i, j, &
    3987              :                                                             n_cstr, n_tiles, unit_nr
    3988              :       LOGICAL :: append_cube, do_cstr_charge_cube, do_dielectric_cube, do_dirichlet_bc_cube, &
    3989              :          has_dirichlet_bc, has_implicit_ps, mpi_io, tile_cubes
    3990        13441 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: zcharge
    3991              :       TYPE(particle_list_type), POINTER                  :: particles
    3992              :       TYPE(pw_env_type), POINTER                         :: pw_env
    3993              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
    3994              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    3995              :       TYPE(pw_r3d_rs_type)                               :: aux_r
    3996              :       TYPE(pw_r3d_rs_type), POINTER                      :: dirichlet_tile
    3997              :       TYPE(qs_subsys_type), POINTER                      :: subsys
    3998              :       TYPE(section_vals_type), POINTER                   :: dft_section
    3999              : 
    4000        13441 :       CALL timeset(routineN, handle)
    4001              : 
    4002        13441 :       NULLIFY (pw_env, auxbas_pw_pool, dft_section, particles)
    4003              : 
    4004        13441 :       dft_section => section_vals_get_subs_vals(input, "DFT")
    4005        13441 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
    4006        13441 :       CALL qs_subsys_get(subsys, particles=particles)
    4007        13441 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
    4008              : 
    4009        13441 :       has_implicit_ps = .FALSE.
    4010        13441 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
    4011        13441 :       IF (pw_env%poisson_env%parameters%solver == pw_poisson_implicit) has_implicit_ps = .TRUE.
    4012              : 
    4013              :       ! Write the dielectric constant into a cube file
    4014              :       do_dielectric_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, &
    4015        13441 :                                                             "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE"), cp_p_file)
    4016        13441 :       IF (has_implicit_ps .AND. do_dielectric_cube) THEN
    4017            2 :          IF (.NOT. ALLOCATED(zcharge)) CALL get_effective_core_charges(qs_env, zcharge)
    4018            2 :          append_cube = section_get_lval(input, "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%APPEND")
    4019            2 :          my_pos_cube = "REWIND"
    4020            2 :          IF (append_cube) THEN
    4021            0 :             my_pos_cube = "APPEND"
    4022              :          END IF
    4023            2 :          mpi_io = .TRUE.
    4024              :          unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE", &
    4025              :                                         extension=".cube", middle_name="DIELECTRIC_CONSTANT", file_position=my_pos_cube, &
    4026            2 :                                         mpi_io=mpi_io)
    4027            2 :          CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
    4028            2 :          CALL auxbas_pw_pool%create_pw(aux_r)
    4029              : 
    4030            2 :          boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
    4031            2 :          SELECT CASE (boundary_condition)
    4032              :          CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
    4033            2 :             CALL pw_copy(poisson_env%implicit_env%dielectric%eps, aux_r)
    4034              :          CASE (MIXED_BC, NEUMANN_BC)
    4035              :             CALL pw_shrink(pw_env%poisson_env%parameters%ps_implicit_params%neumann_directions, &
    4036              :                            pw_env%poisson_env%implicit_env%dct_env%dests_shrink, &
    4037              :                            pw_env%poisson_env%implicit_env%dct_env%srcs_shrink, &
    4038              :                            pw_env%poisson_env%implicit_env%dct_env%bounds_local_shftd, &
    4039            2 :                            poisson_env%implicit_env%dielectric%eps, aux_r)
    4040              :          END SELECT
    4041              : 
    4042              :          CALL cp_pw_to_cube(aux_r, unit_nr, "DIELECTRIC CONSTANT", particles=particles, zeff=zcharge, &
    4043              :                             stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%STRIDE"), &
    4044              :                         max_file_size_mb=section_get_rval(dft_section, "PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%MAX_FILE_SIZE_MB"), &
    4045            2 :                             mpi_io=mpi_io)
    4046              :          CALL cp_print_key_finished_output(unit_nr, logger, input, &
    4047            2 :                                            "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE", mpi_io=mpi_io)
    4048              : 
    4049            2 :          CALL auxbas_pw_pool%give_back_pw(aux_r)
    4050              :       END IF
    4051              : 
    4052              :       ! Write Dirichlet constraint charges into a cube file
    4053              :       do_cstr_charge_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, &
    4054        13441 :                                                              "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE"), cp_p_file)
    4055              : 
    4056        13441 :       has_dirichlet_bc = .FALSE.
    4057        13441 :       IF (has_implicit_ps) THEN
    4058           86 :          boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
    4059           86 :          IF (boundary_condition == MIXED_PERIODIC_BC .OR. boundary_condition == MIXED_BC) THEN
    4060           60 :             has_dirichlet_bc = .TRUE.
    4061              :          END IF
    4062              :       END IF
    4063              : 
    4064           86 :       IF (has_implicit_ps .AND. do_cstr_charge_cube .AND. has_dirichlet_bc) THEN
    4065            2 :          IF (.NOT. ALLOCATED(zcharge)) CALL get_effective_core_charges(qs_env, zcharge)
    4066              :          append_cube = section_get_lval(input, &
    4067            2 :                                         "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%APPEND")
    4068            2 :          my_pos_cube = "REWIND"
    4069            2 :          IF (append_cube) THEN
    4070            0 :             my_pos_cube = "APPEND"
    4071              :          END IF
    4072            2 :          mpi_io = .TRUE.
    4073              :          unit_nr = cp_print_key_unit_nr(logger, input, &
    4074              :                                         "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE", &
    4075              :                                         extension=".cube", middle_name="dirichlet_cstr_charge", file_position=my_pos_cube, &
    4076            2 :                                         mpi_io=mpi_io)
    4077            2 :          CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
    4078            2 :          CALL auxbas_pw_pool%create_pw(aux_r)
    4079              : 
    4080            2 :          boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
    4081            2 :          SELECT CASE (boundary_condition)
    4082              :          CASE (MIXED_PERIODIC_BC)
    4083            2 :             CALL pw_copy(poisson_env%implicit_env%cstr_charge, aux_r)
    4084              :          CASE (MIXED_BC)
    4085              :             CALL pw_shrink(pw_env%poisson_env%parameters%ps_implicit_params%neumann_directions, &
    4086              :                            pw_env%poisson_env%implicit_env%dct_env%dests_shrink, &
    4087              :                            pw_env%poisson_env%implicit_env%dct_env%srcs_shrink, &
    4088              :                            pw_env%poisson_env%implicit_env%dct_env%bounds_local_shftd, &
    4089            2 :                            poisson_env%implicit_env%cstr_charge, aux_r)
    4090              :          END SELECT
    4091              : 
    4092              :          CALL cp_pw_to_cube(aux_r, unit_nr, "DIRICHLET CONSTRAINT CHARGE", particles=particles, zeff=zcharge, &
    4093              :                             stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%STRIDE"), &
    4094              :              max_file_size_mb=section_get_rval(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%MAX_FILE_SIZE_MB"), &
    4095            2 :                             mpi_io=mpi_io)
    4096              :          CALL cp_print_key_finished_output(unit_nr, logger, input, &
    4097            2 :                                            "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE", mpi_io=mpi_io)
    4098              : 
    4099            2 :          CALL auxbas_pw_pool%give_back_pw(aux_r)
    4100              :       END IF
    4101              : 
    4102              :       ! Write Dirichlet type constranits into cube files
    4103              :       do_dirichlet_bc_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, &
    4104        13441 :                                                               "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE"), cp_p_file)
    4105        13441 :       has_dirichlet_bc = .FALSE.
    4106        13441 :       IF (has_implicit_ps) THEN
    4107           86 :          boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
    4108           86 :          IF (boundary_condition == MIXED_PERIODIC_BC .OR. boundary_condition == MIXED_BC) THEN
    4109              :             has_dirichlet_bc = .TRUE.
    4110              :          END IF
    4111              :       END IF
    4112              : 
    4113           60 :       IF (has_implicit_ps .AND. has_dirichlet_bc .AND. do_dirichlet_bc_cube) THEN
    4114            2 :          IF (.NOT. ALLOCATED(zcharge)) CALL get_effective_core_charges(qs_env, zcharge)
    4115            2 :          append_cube = section_get_lval(input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%APPEND")
    4116            2 :          my_pos_cube = "REWIND"
    4117            2 :          IF (append_cube) THEN
    4118            0 :             my_pos_cube = "APPEND"
    4119              :          END IF
    4120            2 :          tile_cubes = section_get_lval(input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%TILE_CUBES")
    4121              : 
    4122            2 :          CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
    4123            2 :          CALL auxbas_pw_pool%create_pw(aux_r)
    4124            2 :          CALL pw_zero(aux_r)
    4125              : 
    4126            2 :          IF (tile_cubes) THEN
    4127              :             ! one cube file per tile
    4128            0 :             n_cstr = SIZE(poisson_env%implicit_env%contacts)
    4129            0 :             DO j = 1, n_cstr
    4130            0 :                n_tiles = poisson_env%implicit_env%contacts(j)%dirichlet_bc%n_tiles
    4131            0 :                DO i = 1, n_tiles
    4132              :                   filename = "dirichlet_cstr_"//TRIM(ADJUSTL(cp_to_string(j)))// &
    4133            0 :                              "_tile_"//TRIM(ADJUSTL(cp_to_string(i)))
    4134            0 :                   mpi_io = .TRUE.
    4135              :                   unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", &
    4136              :                                                  extension=".cube", middle_name=filename, file_position=my_pos_cube, &
    4137            0 :                                                  mpi_io=mpi_io)
    4138              : 
    4139            0 :                   CALL pw_copy(poisson_env%implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, aux_r)
    4140              : 
    4141              :                   CALL cp_pw_to_cube(aux_r, unit_nr, "DIRICHLET TYPE CONSTRAINT", particles=particles, zeff=zcharge, &
    4142              :                                      stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%STRIDE"), &
    4143              :                       max_file_size_mb=section_get_rval(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%MAX_FILE_SIZE_MB"), &
    4144            0 :                                      mpi_io=mpi_io)
    4145              :                   CALL cp_print_key_finished_output(unit_nr, logger, input, &
    4146            0 :                                                     "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", mpi_io=mpi_io)
    4147              :                END DO
    4148              :             END DO
    4149              :          ELSE
    4150              :             ! a single cube file
    4151            2 :             NULLIFY (dirichlet_tile)
    4152            2 :             ALLOCATE (dirichlet_tile)
    4153            2 :             CALL auxbas_pw_pool%create_pw(dirichlet_tile)
    4154            2 :             CALL pw_zero(dirichlet_tile)
    4155            2 :             mpi_io = .TRUE.
    4156              :             unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", &
    4157              :                                            extension=".cube", middle_name="DIRICHLET_CSTR", file_position=my_pos_cube, &
    4158            2 :                                            mpi_io=mpi_io)
    4159              : 
    4160            2 :             n_cstr = SIZE(poisson_env%implicit_env%contacts)
    4161            6 :             DO j = 1, n_cstr
    4162            4 :                n_tiles = poisson_env%implicit_env%contacts(j)%dirichlet_bc%n_tiles
    4163           10 :                DO i = 1, n_tiles
    4164            4 :                   CALL pw_copy(poisson_env%implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, dirichlet_tile)
    4165            8 :                   CALL pw_axpy(dirichlet_tile, aux_r)
    4166              :                END DO
    4167              :             END DO
    4168              : 
    4169              :             CALL cp_pw_to_cube(aux_r, unit_nr, "DIRICHLET TYPE CONSTRAINT", particles=particles, zeff=zcharge, &
    4170              :                                stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%STRIDE"), &
    4171              :                       max_file_size_mb=section_get_rval(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%MAX_FILE_SIZE_MB"), &
    4172            2 :                                mpi_io=mpi_io)
    4173              :             CALL cp_print_key_finished_output(unit_nr, logger, input, &
    4174            2 :                                               "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", mpi_io=mpi_io)
    4175            2 :             CALL auxbas_pw_pool%give_back_pw(dirichlet_tile)
    4176            2 :             DEALLOCATE (dirichlet_tile)
    4177              :          END IF
    4178              : 
    4179            2 :          CALL auxbas_pw_pool%give_back_pw(aux_r)
    4180              :       END IF
    4181              : 
    4182        13441 :       CALL timestop(handle)
    4183              : 
    4184        26882 :    END SUBROUTINE qs_scf_post_ps_implicit
    4185              : 
    4186              : !**************************************************************************************************
    4187              : !> \brief write an adjacency (interaction) matrix
    4188              : !> \param qs_env qs environment
    4189              : !> \param input the input
    4190              : !> \author Mohammad Hossein Bani-Hashemian
    4191              : ! **************************************************************************************************
    4192        13441 :    SUBROUTINE write_adjacency_matrix(qs_env, input)
    4193              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    4194              :       TYPE(section_vals_type), POINTER                   :: input
    4195              : 
    4196              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_adjacency_matrix'
    4197              : 
    4198              :       INTEGER                                            :: adjm_size, colind, handle, iatom, ikind, &
    4199              :                                                             ind, jatom, jkind, k, natom, nkind, &
    4200              :                                                             output_unit, rowind, unit_nr
    4201        13441 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: interact_adjm
    4202              :       LOGICAL                                            :: do_adjm_write, do_symmetric
    4203              :       TYPE(cp_logger_type), POINTER                      :: logger
    4204        13441 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list_a, basis_set_list_b
    4205              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
    4206              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    4207              :       TYPE(neighbor_list_iterator_p_type), &
    4208        13441 :          DIMENSION(:), POINTER                           :: nl_iterator
    4209              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    4210        13441 :          POINTER                                         :: nl
    4211        13441 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    4212              :       TYPE(section_vals_type), POINTER                   :: dft_section
    4213              : 
    4214        13441 :       CALL timeset(routineN, handle)
    4215              : 
    4216        13441 :       NULLIFY (dft_section)
    4217              : 
    4218        13441 :       logger => cp_get_default_logger()
    4219        13441 :       output_unit = cp_logger_get_default_io_unit(logger)
    4220              : 
    4221        13441 :       dft_section => section_vals_get_subs_vals(input, "DFT")
    4222              :       do_adjm_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
    4223        13441 :                                                        "PRINT%ADJMAT_WRITE"), cp_p_file)
    4224              : 
    4225        13441 :       IF (do_adjm_write) THEN
    4226           28 :          NULLIFY (qs_kind_set, nl_iterator)
    4227           28 :          NULLIFY (basis_set_list_a, basis_set_list_b, basis_set_a, basis_set_b)
    4228              : 
    4229           28 :          CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, sab_orb=nl, natom=natom, para_env=para_env)
    4230              : 
    4231           28 :          nkind = SIZE(qs_kind_set)
    4232           28 :          CPASSERT(SIZE(nl) > 0)
    4233           28 :          CALL get_neighbor_list_set_p(neighbor_list_sets=nl, symmetric=do_symmetric)
    4234           28 :          CPASSERT(do_symmetric)
    4235          216 :          ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
    4236           28 :          CALL basis_set_list_setup(basis_set_list_a, "ORB", qs_kind_set)
    4237           28 :          CALL basis_set_list_setup(basis_set_list_b, "ORB", qs_kind_set)
    4238              : 
    4239           28 :          adjm_size = ((natom + 1)*natom)/2
    4240           84 :          ALLOCATE (interact_adjm(4*adjm_size))
    4241           28 :          interact_adjm = 0
    4242              : 
    4243           28 :          NULLIFY (nl_iterator)
    4244           28 :          CALL neighbor_list_iterator_create(nl_iterator, nl)
    4245         2021 :          DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    4246              :             CALL get_iterator_info(nl_iterator, &
    4247              :                                    ikind=ikind, jkind=jkind, &
    4248         1993 :                                    iatom=iatom, jatom=jatom)
    4249              : 
    4250         1993 :             basis_set_a => basis_set_list_a(ikind)%gto_basis_set
    4251         1993 :             IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
    4252         1993 :             basis_set_b => basis_set_list_b(jkind)%gto_basis_set
    4253         1993 :             IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
    4254              : 
    4255              :             ! move everything to the upper triangular part
    4256         1993 :             IF (iatom <= jatom) THEN
    4257              :                rowind = iatom
    4258              :                colind = jatom
    4259              :             ELSE
    4260          670 :                rowind = jatom
    4261          670 :                colind = iatom
    4262              :                ! swap the kinds too
    4263              :                ikind = ikind + jkind
    4264          670 :                jkind = ikind - jkind
    4265          670 :                ikind = ikind - jkind
    4266              :             END IF
    4267              : 
    4268              :             ! indexing upper triangular matrix
    4269         1993 :             ind = adjm_size - (natom - rowind + 1)*((natom - rowind + 1) + 1)/2 + colind - rowind + 1
    4270              :             ! convert the upper triangular matrix into a adjm_size x 4 matrix
    4271              :             ! columns are: iatom, jatom, ikind, jkind
    4272         1993 :             interact_adjm((ind - 1)*4 + 1) = rowind
    4273         1993 :             interact_adjm((ind - 1)*4 + 2) = colind
    4274         1993 :             interact_adjm((ind - 1)*4 + 3) = ikind
    4275         1993 :             interact_adjm((ind - 1)*4 + 4) = jkind
    4276              :          END DO
    4277              : 
    4278           28 :          CALL para_env%sum(interact_adjm)
    4279              : 
    4280              :          unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%ADJMAT_WRITE", &
    4281              :                                         extension=".adjmat", file_form="FORMATTED", &
    4282           28 :                                         file_status="REPLACE")
    4283           28 :          IF (unit_nr > 0) THEN
    4284           14 :             WRITE (unit_nr, "(1A,2X,1A,5X,1A,4X,A5,3X,A5)") "#", "iatom", "jatom", "ikind", "jkind"
    4285           88 :             DO k = 1, 4*adjm_size, 4
    4286              :                ! print only the interacting atoms
    4287           88 :                IF (interact_adjm(k) > 0 .AND. interact_adjm(k + 1) > 0) THEN
    4288           74 :                   WRITE (unit_nr, "(I8,2X,I8,3X,I6,2X,I6)") interact_adjm(k:k + 3)
    4289              :                END IF
    4290              :             END DO
    4291              :          END IF
    4292              : 
    4293           28 :          CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%ADJMAT_WRITE")
    4294              : 
    4295           28 :          CALL neighbor_list_iterator_release(nl_iterator)
    4296           56 :          DEALLOCATE (basis_set_list_a, basis_set_list_b)
    4297              :       END IF
    4298              : 
    4299        13441 :       CALL timestop(handle)
    4300              : 
    4301        26882 :    END SUBROUTINE write_adjacency_matrix
    4302              : 
    4303              : ! **************************************************************************************************
    4304              : !> \brief Updates Hartree potential with MP2 density. Important for REPEAT charges
    4305              : !> \param rho ...
    4306              : !> \param qs_env ...
    4307              : !> \author Vladimir Rybkin
    4308              : ! **************************************************************************************************
    4309          322 :    SUBROUTINE update_hartree_with_mp2(rho, qs_env)
    4310              :       TYPE(qs_rho_type), POINTER                         :: rho
    4311              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    4312              : 
    4313              :       LOGICAL                                            :: use_virial
    4314              :       TYPE(pw_c1d_gs_type)                               :: rho_tot_gspace, v_hartree_gspace
    4315              :       TYPE(pw_c1d_gs_type), POINTER                      :: rho_core
    4316              :       TYPE(pw_env_type), POINTER                         :: pw_env
    4317              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
    4318              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    4319              :       TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_rspace
    4320              :       TYPE(qs_energy_type), POINTER                      :: energy
    4321              :       TYPE(virial_type), POINTER                         :: virial
    4322              : 
    4323          322 :       NULLIFY (auxbas_pw_pool, pw_env, poisson_env, energy, rho_core, v_hartree_rspace, virial)
    4324              :       CALL get_qs_env(qs_env, pw_env=pw_env, energy=energy, &
    4325              :                       rho_core=rho_core, virial=virial, &
    4326          322 :                       v_hartree_rspace=v_hartree_rspace)
    4327              : 
    4328          322 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
    4329              : 
    4330              :       IF (.NOT. use_virial) THEN
    4331              : 
    4332              :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
    4333          268 :                          poisson_env=poisson_env)
    4334          268 :          CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
    4335          268 :          CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
    4336              : 
    4337          268 :          CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
    4338              :          CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%hartree, &
    4339          268 :                                v_hartree_gspace, rho_core=rho_core)
    4340              : 
    4341          268 :          CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
    4342          268 :          CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
    4343              : 
    4344          268 :          CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
    4345          268 :          CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
    4346              :       END IF
    4347              : 
    4348          322 :    END SUBROUTINE update_hartree_with_mp2
    4349              : 
    4350            0 : END MODULE qs_scf_post_gpw
        

Generated by: LCOV version 2.0-1