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