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