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