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