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