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