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