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