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