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