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 24684 : 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 24684 : IF (LEN(TRIM(extend_by)) > 0 .AND. extend_by(1:1) == "%") THEN
305 24684 : 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 24684 : 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 530 : FUNCTION cp_section_key_do_cubes(self) RESULT(res)
317 : CLASS(cp_section_key) :: self
318 : LOGICAL :: res
319 :
320 530 : res = self%do_output .AND. self%grid_output == grid_output_cubes
321 530 : 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 530 : FUNCTION cp_section_key_do_openpmd(self) RESULT(res)
329 : CLASS(cp_section_key) :: self
330 : LOGICAL :: res
331 :
332 530 : res = self%do_output .AND. self%grid_output == grid_output_openpmd
333 530 : 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 36227 : 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 36227 : "DFT%"//TRIM(ADJUSTL(str_cubes))), cp_p_file)
538 : do_openpmd = BTEST(cp_print_key_should_output( &
539 : logger%iter_info, input, &
540 36227 : "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 36227 : CPASSERT(.NOT. (do_cubes .AND. do_openpmd))
545 36227 : res%do_output = do_cubes .OR. do_openpmd
546 36227 : 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 36227 : res%grid_output = grid_output_cubes
552 36227 : res%relative_section_key = TRIM(ADJUSTL(str_cubes))
553 36227 : res%format_name = "Cube"
554 : END IF
555 36227 : res%absolute_section_key = "DFT%"//TRIM(ADJUSTL(res%relative_section_key))
556 36227 : 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 11629 : 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 11629 : 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 11629 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
642 : TYPE(admm_type), POINTER :: admm_env
643 11629 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
644 11629 : TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: mixed_evals, occupied_evals, &
645 11629 : unoccupied_evals, unoccupied_evals_stm
646 11629 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mixed_orbs, occupied_orbs
647 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
648 11629 : TARGET :: homo_localized, lumo_localized, &
649 11629 : mixed_localized
650 11629 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: lumo_ptr, mo_loc_history, &
651 11629 : 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 11629 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, matrix_p_mp2, matrix_s, &
656 11629 : mo_derivs
657 11629 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: kinetic_m, rho_ao
658 : TYPE(dft_control_type), POINTER :: dft_control
659 11629 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
660 11629 : TYPE(molecule_type), POINTER :: molecule_set(:)
661 : TYPE(mp_para_env_type), POINTER :: para_env
662 : TYPE(particle_list_type), POINTER :: particles
663 11629 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
664 : TYPE(pw_c1d_gs_type) :: wf_g
665 : TYPE(pw_env_type), POINTER :: pw_env
666 11629 : 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 11629 : 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 11629 : CALL timeset(routineN, handle)
682 :
683 11629 : logger => cp_get_default_logger()
684 11629 : 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 11629 : my_do_mp2 = .FALSE.
688 11629 : IF (PRESENT(do_mp2)) my_do_mp2 = do_mp2
689 11629 : 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 11629 : CALL get_qs_env(qs_env, scf_env=scf_env)
699 :
700 11629 : my_localized_wfn = .FALSE.
701 11629 : NULLIFY (admm_env, dft_control, pw_env, auxbas_pw_pool, pw_pools, mos, rho, &
702 11629 : mo_coeff, ks_rmpv, matrix_s, qs_loc_env_homo, qs_loc_env_lumo, scf_control, &
703 11629 : unoccupied_orbs, mo_eigenvalues, unoccupied_evals, &
704 11629 : unoccupied_evals_stm, molecule_set, mo_derivs, &
705 11629 : subsys, particles, input, print_key, kinetic_m, marked_states, &
706 11629 : mixed_evals, qs_loc_env_mixed)
707 11629 : NULLIFY (lumo_ptr, rho_ao)
708 :
709 11629 : has_homo = .FALSE.
710 11629 : has_lumo = .FALSE.
711 11629 : p_loc = .FALSE.
712 11629 : p_loc_homo = .FALSE.
713 11629 : p_loc_lumo = .FALSE.
714 11629 : p_loc_mixed = .FALSE.
715 :
716 11629 : CPASSERT(ASSOCIATED(scf_env))
717 11629 : 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 11629 : qs_kind_set=qs_kind_set)
731 11629 : rtp_control => dft_control%rtp_control
732 11629 : CALL qs_subsys_get(subsys, particles=particles)
733 :
734 11629 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
735 :
736 11629 : 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 11629 : CALL write_available_results(qs_env, scf_env)
749 :
750 : ! **** the kinetic energy
751 11629 : 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 11629 : CALL qs_scf_post_charges(input, logger, qs_env)
768 :
769 : ! Moments of charge distribution
770 11629 : 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 11629 : dft_section => section_vals_get_subs_vals(input, "DFT")
774 11629 : localize_section => section_vals_get_subs_vals(dft_section, "LOCALIZE")
775 11629 : loc_print_section => section_vals_get_subs_vals(localize_section, "PRINT")
776 11629 : CALL section_vals_get(localize_section, explicit=loc_explicit)
777 11629 : CALL section_vals_get(loc_print_section, explicit=loc_print_explicit)
778 :
779 : ! Print_keys controlled by localization
780 11629 : 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 11629 : 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 11529 : p_loc_homo = .FALSE.
811 11529 : p_loc_lumo = .FALSE.
812 11529 : p_loc_mixed = .FALSE.
813 11529 : n_rep = 0
814 : END IF
815 :
816 11629 : 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 11629 : stm_section => section_vals_get_subs_vals(input, "DFT%PRINT%STM")
824 11629 : CALL section_vals_get(stm_section, explicit=do_stm)
825 11629 : nlumo_stm = 0
826 11629 : IF (do_stm) nlumo_stm = section_get_ival(stm_section, "NLUMO")
827 :
828 : ! check for CUBES or openPMD (MOs and WANNIERS)
829 11629 : mo_section = cube_or_openpmd(input, str_mo_cubes, str_mo_openpmd, logger)
830 :
831 11629 : 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 11629 : nlumo = section_get_ival(dft_section, mo_section%concat_to_relative("%NLUMO"))
838 11629 : 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 11629 : 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 11629 : IF (dft_control%restricted) THEN
849 : !For ROKS useful only first term
850 74 : nspins = 1
851 : ELSE
852 11555 : nspins = dft_control%nspins
853 : END IF
854 : !Some info about ROKS
855 11629 : 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 11629 : IF (do_kpoints) THEN
861 530 : CPWARN_IF(mo_section%do_cubes(), warning_cube_kpoint)
862 530 : CPWARN_IF(mo_section%do_openpmd(), warning_openpmd_kpoint)
863 : ELSE
864 : CALL get_qs_env(qs_env, &
865 : mos=mos, &
866 11099 : matrix_ks=ks_rmpv)
867 11099 : 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 11099 : 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 11629 : IF (p_loc_homo) THEN
904 94 : IF (do_kpoints) THEN
905 0 : CPWARN("Localization not implemented for k-point calculations!")
906 : ELSE IF (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 11629 : IF (do_kpoints) THEN
955 530 : 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 11099 : compute_lumos = mo_section%do_output .AND. nlumo /= 0
961 11099 : compute_lumos = compute_lumos .OR. p_loc_lumo
962 :
963 24256 : DO ispin = 1, dft_control%nspins
964 13157 : CALL get_mo_set(mo_set=mos(ispin), homo=homo, nmo=nmo)
965 37365 : compute_lumos = compute_lumos .AND. homo == nmo
966 : END DO
967 :
968 11099 : 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 11067 : 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 11629 : IF (p_loc_mixed) THEN
1054 2 : IF (do_kpoints) THEN
1055 0 : CPWARN("Localization not implemented for k-point calculations!")
1056 2 : ELSE IF (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 11629 : 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 11629 : IF (.NOT. do_kpoints) THEN
1114 11099 : 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 11099 : 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 11099 : 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 11629 : IF (do_kpoints) THEN
1130 : ! nothing at the moment, not implemented
1131 : ELSE
1132 11099 : 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 11099 : matrix_s=matrix_s, marked_states=marked_states)
1136 :
1137 11099 : IF (p_loc_lumo) CALL cp_fm_release(lumo_localized)
1138 : END IF
1139 11629 : 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 11629 : IF (.NOT. do_kpoints) THEN
1145 11099 : 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 11629 : 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 11629 : CALL qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit)
1183 :
1184 : ! Calculation of Electric Field Gradients
1185 11629 : CALL qs_scf_post_efg(input, logger, qs_env)
1186 :
1187 : ! Calculation of ET
1188 11629 : CALL qs_scf_post_et(input, qs_env, dft_control)
1189 :
1190 : ! Calculation of EPR Hyperfine Coupling Tensors
1191 11629 : CALL qs_scf_post_epr(input, logger, qs_env)
1192 :
1193 : ! Calculation of properties needed for BASIS_MOLOPT optimizations
1194 11629 : CALL qs_scf_post_molopt(input, logger, qs_env)
1195 :
1196 : ! Calculate ELF
1197 11629 : CALL qs_scf_post_elf(input, logger, qs_env)
1198 :
1199 : ! Use Wannier90 interface
1200 11629 : CALL wannier90_interface(input, logger, qs_env)
1201 :
1202 11629 : 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 11629 : CALL cp_openpmd_close_iterations()
1212 :
1213 11629 : CALL timestop(handle)
1214 :
1215 23258 : 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 11629 : 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 11629 : CALL timeset(routineN, handle)
1344 :
1345 11629 : 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 11629 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%LOWDIN")
1351 11629 : 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 11629 : CALL resp_fit(qs_env)
1365 :
1366 : ! Compute the Density Derived Atomic Point charges with the Bloechl scheme
1367 11629 : print_key => section_vals_get_subs_vals(input, "PROPERTIES%FIT_CHARGE")
1368 11629 : 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 11629 : CALL timestop(handle)
1377 :
1378 11629 : 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 12909 : 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 12909 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
1637 : TYPE(section_vals_type), POINTER :: print_key
1638 :
1639 12909 : CALL timeset(routineN, handle)
1640 :
1641 : print_key => section_vals_get_subs_vals(section_vals=input, &
1642 12909 : subsection_name="DFT%PRINT%MOMENTS")
1643 :
1644 12909 : 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 1466 : keyword_name="DFT%PRINT%MOMENTS%MAX_MOMENT")
1648 : periodic = section_get_lval(section_vals=input, &
1649 1466 : keyword_name="DFT%PRINT%MOMENTS%PERIODIC")
1650 : reference = section_get_ival(section_vals=input, &
1651 1466 : keyword_name="DFT%PRINT%MOMENTS%REFERENCE")
1652 : magnetic = section_get_lval(section_vals=input, &
1653 1466 : keyword_name="DFT%PRINT%MOMENTS%MAGNETIC")
1654 : vel_reprs = section_get_lval(section_vals=input, &
1655 1466 : keyword_name="DFT%PRINT%MOMENTS%VEL_REPRS")
1656 : com_nl = section_get_lval(section_vals=input, &
1657 1466 : keyword_name="DFT%PRINT%MOMENTS%COM_NL")
1658 : second_ref_point = section_get_lval(section_vals=input, &
1659 1466 : keyword_name="DFT%PRINT%MOMENTS%SECOND_REFERENCE_POINT")
1660 : do_kg = section_get_lval(section_vals=input, &
1661 1466 : keyword_name="DFT%PRINT%MOMENTS%KG")
1662 : max_nmo = section_get_ival(section_vals=input, &
1663 1466 : keyword_name="DFT%PRINT%MOMENTS%MAX_NMO")
1664 :
1665 1466 : NULLIFY (ref_point)
1666 1466 : 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 1466 : middle_name="moments", log_filename=.FALSE.)
1670 :
1671 1466 : IF (output_unit > 0) THEN
1672 743 : 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 704 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS"
1679 : END IF
1680 : END IF
1681 :
1682 1466 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
1683 :
1684 1466 : IF (do_kpoints) THEN
1685 10 : CALL qs_moment_kpoints(qs_env, maxmom, reference, ref_point, max_nmo, unit_nr)
1686 : ELSE
1687 1456 : IF (periodic) THEN
1688 472 : CALL qs_moment_berry_phase(qs_env, magnetic, maxmom, reference, ref_point, unit_nr)
1689 : ELSE
1690 984 : CALL qs_moment_locop(qs_env, magnetic, maxmom, reference, ref_point, unit_nr, vel_reprs, com_nl)
1691 : END IF
1692 1456 : 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 1466 : basis_section=input, print_key_path="DFT%PRINT%MOMENTS")
1699 :
1700 1466 : 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 12909 : CALL timestop(handle)
1736 :
1737 12909 : 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 11629 : 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 11629 : CALL timeset(routineN, handle)
1762 :
1763 : print_key => section_vals_get_subs_vals(section_vals=input, &
1764 11629 : subsection_name="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM")
1765 :
1766 11629 : 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 11629 : CALL timestop(handle)
1795 :
1796 11629 : 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 11629 : 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 11629 : CALL timeset(routineN, handle)
1815 :
1816 : print_key => section_vals_get_subs_vals(section_vals=input, &
1817 11629 : subsection_name="DFT%PRINT%ELECTRIC_FIELD_GRADIENT")
1818 11629 : 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 11629 : CALL timestop(handle)
1824 :
1825 11629 : 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 23258 : 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 11629 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: my_mos
1843 : TYPE(section_vals_type), POINTER :: et_section
1844 :
1845 11629 : CALL timeset(routineN, handle)
1846 :
1847 : do_et = .FALSE.
1848 11629 : et_section => section_vals_get_subs_vals(input, "PROPERTIES%ET_COUPLING")
1849 11629 : CALL section_vals_get(et_section, explicit=do_et)
1850 11629 : 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 11629 : CALL timestop(handle)
1868 :
1869 11629 : 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 11629 : 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 11629 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1897 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1898 11629 : 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 11629 : CALL timeset(routineN, handle)
1903 11629 : output_unit = cp_logger_get_default_io_unit(logger)
1904 :
1905 11629 : elf_section_key = cube_or_openpmd(input, str_elf_cubes, str_elf_openpmd, logger)
1906 :
1907 11629 : elf_section => section_vals_get_subs_vals(input, elf_section_key%absolute_section_key)
1908 11629 : 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 11629 : CALL timestop(handle)
1995 :
1996 23258 : 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 11629 : 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 11629 : 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 11629 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2022 11629 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2023 : TYPE(qs_energy_type), POINTER :: energy
2024 : TYPE(section_vals_type), POINTER :: print_key
2025 :
2026 11629 : CALL timeset(routineN, handle)
2027 :
2028 : print_key => section_vals_get_subs_vals(section_vals=input, &
2029 11629 : subsection_name="DFT%PRINT%BASIS_MOLOPT_QUANTITIES")
2030 11629 : 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 11629 : CALL timestop(handle)
2071 :
2072 23258 : 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 11629 : 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 11629 : CALL timeset(routineN, handle)
2091 :
2092 : print_key => section_vals_get_subs_vals(section_vals=input, &
2093 11629 : subsection_name="DFT%PRINT%HYPERFINE_COUPLING_TENSOR")
2094 11629 : 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 11629 : CALL timestop(handle)
2100 :
2101 11629 : 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 11629 : 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 11629 : CALL timeset(routineN, handle)
2119 :
2120 : ! those properties that require MOs (not suitable density matrix based methods)
2121 11629 : 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 11629 : CALL write_mo_free_results(qs_env)
2125 :
2126 11629 : CALL timestop(handle)
2127 :
2128 11629 : 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 11953 : 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_curve, &
2147 : do_dos, do_kpoints, do_pdos, &
2148 : do_projected_dos, explicit
2149 : REAL(KIND=dp) :: maxocc, s_square, s_square_ideal, &
2150 : total_abs_spin_dens, total_spin_dens
2151 11953 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, occupation_numbers
2152 : TYPE(admm_type), POINTER :: admm_env
2153 11953 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2154 : TYPE(cell_type), POINTER :: cell
2155 11953 : TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: unoccupied_evals
2156 11953 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: unoccupied_orbs
2157 : TYPE(cp_fm_type), POINTER :: mo_coeff
2158 : TYPE(cp_logger_type), POINTER :: logger
2159 11953 : 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 11953 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2163 11953 : TYPE(molecule_type), POINTER :: molecule_set(:)
2164 : TYPE(particle_list_type), POINTER :: particles
2165 11953 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2166 : TYPE(pw_env_type), POINTER :: pw_env
2167 11953 : 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 11953 : 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 11953 : 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 11953 : CALL timeset(routineN, handle)
2184 :
2185 11953 : NULLIFY (cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, mo_coeff, &
2186 11953 : mo_coeff_deriv, mo_eigenvalues, mos, atomic_kind_set, qs_kind_set, &
2187 11953 : particle_set, rho, ks_rmpv, matrix_s, scf_control, dft_section, &
2188 11953 : molecule_set, input, particles, subsys, rho_r, unoccupied_orbs, &
2189 11953 : unoccupied_evals, casino_section, dos_section)
2190 :
2191 11953 : logger => cp_get_default_logger()
2192 11953 : output_unit = cp_logger_get_default_io_unit(logger)
2193 :
2194 11953 : 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 11953 : subsys=subsys)
2206 11953 : CALL qs_subsys_get(subsys, particles=particles)
2207 11953 : CALL get_qs_env(qs_env, rho=rho)
2208 11953 : CALL qs_rho_get(rho, rho_r=rho_r)
2209 :
2210 : ! k points
2211 11953 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
2212 :
2213 : ! Write last MO information to output file if requested
2214 11953 : dft_section => section_vals_get_subs_vals(input, "DFT")
2215 11953 : IF (.NOT. qs_env%run_rtp) THEN
2216 11629 : CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.TRUE.)
2217 11629 : trexio_section => section_vals_get_subs_vals(dft_section, "PRINT%TREXIO")
2218 11629 : CALL section_vals_get(trexio_section, explicit=explicit)
2219 11629 : IF (explicit) THEN
2220 10 : CALL write_trexio(qs_env, trexio_section)
2221 : END IF
2222 11629 : casino_section => section_vals_get_subs_vals(dft_section, "PRINT%CASINO")
2223 11629 : CALL section_vals_get(casino_section, explicit=explicit)
2224 11629 : IF (explicit) THEN
2225 10 : CALL write_casino(qs_env, casino_section)
2226 : END IF
2227 11629 : sprint_section => section_vals_get_subs_vals(dft_section, "PRINT%MO_MOLDEN")
2228 11629 : defer_molden = .FALSE.
2229 11629 : IF (.NOT. do_kpoints) THEN
2230 11099 : CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv)
2231 11099 : CALL write_dm_binary_restart(mos, dft_section, ks_rmpv)
2232 : ! Check if molden write should be deferred for OT unoccupied orbitals
2233 11099 : CALL section_vals_val_get(sprint_section, "NLUMO", i_val=nlumo_molden)
2234 11099 : 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 11099 : qs_env=qs_env, calc_energies=.TRUE.)
2240 : END IF
2241 : ! Write Chargemol .wfx
2242 11099 : 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 530 : 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 530 : 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 11629 : 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 11629 : dos_section => section_vals_get_subs_vals(dft_section, "PRINT%DOS")
2268 11629 : do_dos = BTEST(cp_print_key_should_output(logger%iter_info, dos_section), cp_p_file)
2269 11629 : CALL get_dos_pdos_flags(dos_section, do_dos, do_projected_dos, do_pdos, do_curve)
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 11629 : nlumo_required = 0
2274 11629 : IF (.NOT. do_kpoints .AND. PRESENT(scf_env)) THEN
2275 11099 : IF (scf_env%method == ot_method_nr) THEN
2276 4287 : 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 : ELSE IF (nlumo_required /= -1) THEN
2281 8 : nlumo_required = MAX(nlumo_required, nlumo_dos)
2282 : END IF
2283 : END IF
2284 :
2285 4287 : IF (defer_molden) THEN
2286 0 : IF (nlumo_molden == -1) THEN
2287 0 : nlumo_required = -1
2288 0 : ELSE IF (nlumo_required /= -1) THEN
2289 0 : nlumo_required = MAX(nlumo_required, nlumo_molden)
2290 : END IF
2291 : END IF
2292 4287 : 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 4287 : 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 4295 : 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 11629 : 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 11629 : IF (do_dos) THEN
2340 64 : IF (do_kpoints) THEN
2341 4 : CALL calculate_dos_kp(qs_env, dft_section)
2342 4 : IF (do_curve) CALL calculate_dos_kp(qs_env, dft_section, write_curve_output=.TRUE.)
2343 : ELSE
2344 60 : CALL get_qs_env(qs_env, mos=mos)
2345 60 : IF (ASSOCIATED(unoccupied_evals)) THEN
2346 : CALL calculate_dos(mos, dft_section, unoccupied_evals=unoccupied_evals, &
2347 2 : smearing_enabled=dft_control%smear)
2348 2 : IF (do_curve) CALL calculate_dos(mos, dft_section, unoccupied_evals=unoccupied_evals, &
2349 0 : smearing_enabled=dft_control%smear, write_curve_output=.TRUE.)
2350 : ELSE
2351 58 : CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear)
2352 58 : IF (do_curve) CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear, &
2353 0 : write_curve_output=.TRUE.)
2354 : END IF
2355 : END IF
2356 : END IF
2357 :
2358 : ! Print projected density-of-states outputs
2359 11629 : IF (do_projected_dos) THEN
2360 22 : IF (do_kpoints) THEN
2361 : CALL calculate_projected_dos_kp(qs_env, dft_section, pdos_print_key="PRINT%DOS", &
2362 2 : write_pdos=do_pdos, write_pdos_curve=do_curve)
2363 : ELSE
2364 : CALL get_qs_env(qs_env, &
2365 : mos=mos, &
2366 20 : matrix_ks=ks_rmpv)
2367 40 : DO ispin = 1, dft_control%nspins
2368 40 : IF (dft_control%nspins == 2) THEN
2369 0 : IF (ASSOCIATED(unoccupied_orbs)) THEN
2370 : CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
2371 : qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin, &
2372 : unoccupied_orbs=unoccupied_orbs(ispin), &
2373 : unoccupied_evals=unoccupied_evals(ispin), &
2374 0 : pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
2375 : ELSE
2376 : CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
2377 : qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin, &
2378 0 : pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
2379 : END IF
2380 : ELSE
2381 20 : IF (ASSOCIATED(unoccupied_orbs)) THEN
2382 : CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
2383 : qs_kind_set, particle_set, qs_env, dft_section, &
2384 : unoccupied_orbs=unoccupied_orbs(ispin), &
2385 : unoccupied_evals=unoccupied_evals(ispin), &
2386 2 : pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
2387 : ELSE
2388 : CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
2389 : qs_kind_set, particle_set, qs_env, dft_section, &
2390 18 : pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
2391 : END IF
2392 : END IF
2393 : END DO
2394 : END IF
2395 : END IF
2396 11629 : IF (ASSOCIATED(unoccupied_orbs)) THEN
2397 4 : DO ispin = 1, dft_control%nspins
2398 2 : DEALLOCATE (unoccupied_evals(ispin)%array)
2399 4 : CALL cp_fm_release(unoccupied_orbs(ispin))
2400 : END DO
2401 2 : DEALLOCATE (unoccupied_evals)
2402 2 : DEALLOCATE (unoccupied_orbs)
2403 : END IF
2404 : END IF
2405 :
2406 : ! Integrated absolute spin density and spin contamination ***
2407 11953 : IF (dft_control%nspins == 2) THEN
2408 2164 : total_spin_dens = 0.0_dp
2409 2164 : total_abs_spin_dens = 0.0_dp
2410 2164 : IF (dft_control%qs_control%gapw) THEN
2411 336 : CALL get_qs_env(qs_env, qs_charges=qs_charges)
2412 : total_spin_dens = qs_charges%total_rho_hard_spin - &
2413 336 : qs_charges%total_rho_soft_spin
2414 : total_abs_spin_dens = qs_charges%total_rho_hard_abs_spin - &
2415 336 : qs_charges%total_rho_soft_abs_spin
2416 : END IF
2417 2164 : CALL get_qs_env(qs_env, mos=mos)
2418 2164 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
2419 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
2420 2164 : pw_pools=pw_pools)
2421 2164 : CALL auxbas_pw_pool%create_pw(wf_r)
2422 2164 : CALL pw_copy(rho_r(1), wf_r)
2423 2164 : CALL pw_axpy(rho_r(2), wf_r, alpha=-1._dp)
2424 2164 : total_spin_dens = total_spin_dens + pw_integrate_function(wf_r)
2425 2164 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,F20.10))') &
2426 1105 : "Integrated spin density: ", total_spin_dens
2427 2164 : total_abs_spin_dens = total_abs_spin_dens + pw_integrate_function(wf_r, oprt="ABS")
2428 2164 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T3,A,T61,F20.10))') &
2429 1105 : "Integrated absolute spin density: ", total_abs_spin_dens
2430 2164 : CALL auxbas_pw_pool%give_back_pw(wf_r)
2431 : !
2432 : ! XXX Fix Me XXX
2433 : ! should be extended to the case where added MOs are present
2434 : ! should be extended to the k-point case
2435 : !
2436 2164 : IF (.NOT. do_kpoints) THEN
2437 2124 : all_equal = .TRUE.
2438 6372 : DO ispin = 1, dft_control%nspins
2439 : CALL get_mo_set(mo_set=mos(ispin), &
2440 : occupation_numbers=occupation_numbers, &
2441 : homo=homo, &
2442 : nmo=nmo, &
2443 4248 : maxocc=maxocc)
2444 6372 : IF (nmo > 0) THEN
2445 : all_equal = all_equal .AND. &
2446 : (ALL(occupation_numbers(1:homo) == maxocc) .AND. &
2447 24482 : ALL(occupation_numbers(homo + 1:nmo) == 0.0_dp))
2448 : END IF
2449 : END DO
2450 2124 : IF (all_equal) THEN
2451 : CALL get_qs_env(qs_env=qs_env, &
2452 : matrix_s=matrix_s, &
2453 2006 : energy=energy)
2454 : CALL compute_s_square(mos=mos, matrix_s=matrix_s, s_square=s_square, &
2455 2006 : s_square_ideal=s_square_ideal)
2456 2006 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(T3,A,T51,2F15.6)') &
2457 1026 : "Ideal and single determinant S**2 : ", s_square_ideal, s_square
2458 2006 : energy%s_square = s_square
2459 : END IF
2460 : END IF
2461 : END IF
2462 :
2463 11953 : CALL timestop(handle)
2464 :
2465 11953 : END SUBROUTINE write_mo_dependent_results
2466 :
2467 : ! **************************************************************************************************
2468 : !> \brief Write QS results always available (if switched on through the print_keys)
2469 : !> Can be called from ls_scf
2470 : !> \param qs_env the qs_env in which the qs_env lives
2471 : ! **************************************************************************************************
2472 12969 : SUBROUTINE write_mo_free_results(qs_env)
2473 : TYPE(qs_environment_type), POINTER :: qs_env
2474 :
2475 : CHARACTER(len=*), PARAMETER :: routineN = 'write_mo_free_results'
2476 : CHARACTER(len=1), DIMENSION(3), PARAMETER :: cdir = ["x", "y", "z"]
2477 :
2478 : CHARACTER(LEN=2) :: element_symbol
2479 : CHARACTER(LEN=default_path_length) :: filename, mpi_filename, my_pos_cube, &
2480 : my_pos_voro
2481 : CHARACTER(LEN=default_string_length) :: name, print_density
2482 : INTEGER :: after, handle, i, iat, iatom, id, ikind, img, iso, ispin, iw, l, n_rep_hf, nat, &
2483 : natom, nd(3), ngto, niso, nkind, np, nr, output_unit, print_level, should_print_bqb, &
2484 : should_print_voro, unit_nr, unit_nr_voro
2485 : LOGICAL :: append_cube, append_voro, do_hfx, do_kpoints, mpi_io, omit_headers, print_it, &
2486 : rho_r_valid, voro_print_txt, write_ks, write_xc, xrd_interface
2487 : REAL(KIND=dp) :: norm_factor, q_max, rho_hard, rho_soft, &
2488 : rho_total, rho_total_rspace, udvol, &
2489 : volume, zeff
2490 12969 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zcharge
2491 12969 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: bfun
2492 12969 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: aedens, ccdens, ppdens
2493 : REAL(KIND=dp), DIMENSION(3) :: dr
2494 12969 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_Q0
2495 12969 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2496 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2497 : TYPE(cell_type), POINTER :: cell
2498 : TYPE(cp_logger_type), POINTER :: logger
2499 : TYPE(cp_section_key) :: e_density_section
2500 12969 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hr
2501 12969 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_rmpv, matrix_vxc, rho_ao
2502 : TYPE(dft_control_type), POINTER :: dft_control
2503 : TYPE(grid_atom_type), POINTER :: grid_atom
2504 : TYPE(iao_env_type) :: iao_env
2505 : TYPE(mp_para_env_type), POINTER :: para_env
2506 : TYPE(particle_list_type), POINTER :: particles
2507 12969 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2508 : TYPE(pw_c1d_gs_type) :: aux_g, rho_elec_gspace
2509 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
2510 : TYPE(pw_env_type), POINTER :: pw_env
2511 12969 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
2512 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2513 : TYPE(pw_r3d_rs_type) :: aux_r, rho_elec_rspace, wf_r
2514 12969 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
2515 : TYPE(pw_r3d_rs_type), POINTER :: mb_rho, v_hartree_rspace, vee
2516 12969 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2517 : TYPE(qs_kind_type), POINTER :: qs_kind
2518 : TYPE(qs_rho_type), POINTER :: rho
2519 : TYPE(qs_subsys_type), POINTER :: subsys
2520 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
2521 12969 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
2522 : TYPE(rho_atom_type), POINTER :: rho_atom
2523 : TYPE(section_vals_type), POINTER :: dft_section, hfx_section, input, &
2524 : print_key, print_key_bqb, &
2525 : print_key_voro, xc_section
2526 :
2527 12969 : CALL timeset(routineN, handle)
2528 12969 : NULLIFY (cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, hfx_section, &
2529 12969 : atomic_kind_set, qs_kind_set, particle_set, rho, ks_rmpv, rho_ao, rho_r, &
2530 12969 : dft_section, xc_section, input, particles, subsys, matrix_vxc, v_hartree_rspace, &
2531 12969 : vee)
2532 :
2533 12969 : logger => cp_get_default_logger()
2534 12969 : output_unit = cp_logger_get_default_io_unit(logger)
2535 :
2536 12969 : CPASSERT(ASSOCIATED(qs_env))
2537 : CALL get_qs_env(qs_env, &
2538 : atomic_kind_set=atomic_kind_set, &
2539 : qs_kind_set=qs_kind_set, &
2540 : particle_set=particle_set, &
2541 : cell=cell, &
2542 : para_env=para_env, &
2543 : dft_control=dft_control, &
2544 : input=input, &
2545 : do_kpoints=do_kpoints, &
2546 12969 : subsys=subsys)
2547 12969 : dft_section => section_vals_get_subs_vals(input, "DFT")
2548 12969 : CALL qs_subsys_get(subsys, particles=particles)
2549 :
2550 12969 : CALL get_qs_env(qs_env, rho=rho)
2551 12969 : CALL qs_rho_get(rho, rho_r=rho_r)
2552 :
2553 12969 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
2554 38907 : ALLOCATE (zcharge(natom))
2555 36085 : DO ikind = 1, nkind
2556 23116 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
2557 23116 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
2558 85830 : DO iatom = 1, nat
2559 49745 : iat = atomic_kind_set(ikind)%atom_list(iatom)
2560 72861 : zcharge(iat) = zeff
2561 : END DO
2562 : END DO
2563 :
2564 : ! Print the total density (electronic + core charge)
2565 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
2566 : "DFT%PRINT%TOT_DENSITY_CUBE"), cp_p_file)) THEN
2567 82 : NULLIFY (rho_core, rho0_s_gs, rhoz_cneo_s_gs)
2568 82 : append_cube = section_get_lval(input, "DFT%PRINT%TOT_DENSITY_CUBE%APPEND")
2569 82 : my_pos_cube = "REWIND"
2570 82 : IF (append_cube) THEN
2571 0 : my_pos_cube = "APPEND"
2572 : END IF
2573 :
2574 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho_core=rho_core, &
2575 82 : rho0_s_gs=rho0_s_gs, rhoz_cneo_s_gs=rhoz_cneo_s_gs)
2576 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
2577 82 : pw_pools=pw_pools)
2578 82 : CALL auxbas_pw_pool%create_pw(wf_r)
2579 82 : IF (dft_control%qs_control%gapw) THEN
2580 0 : IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
2581 0 : CALL pw_axpy(rho_core, rho0_s_gs)
2582 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
2583 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
2584 : END IF
2585 0 : CALL pw_transfer(rho0_s_gs, wf_r)
2586 0 : CALL pw_axpy(rho_core, rho0_s_gs, -1.0_dp)
2587 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
2588 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
2589 : END IF
2590 : ELSE
2591 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
2592 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
2593 : END IF
2594 0 : CALL pw_transfer(rho0_s_gs, wf_r)
2595 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
2596 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
2597 : END IF
2598 : END IF
2599 : ELSE
2600 82 : CALL pw_transfer(rho_core, wf_r)
2601 : END IF
2602 164 : DO ispin = 1, dft_control%nspins
2603 164 : CALL pw_axpy(rho_r(ispin), wf_r)
2604 : END DO
2605 82 : filename = "TOTAL_DENSITY"
2606 82 : mpi_io = .TRUE.
2607 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%TOT_DENSITY_CUBE", &
2608 : extension=".cube", middle_name=TRIM(filename), file_position=my_pos_cube, &
2609 82 : log_filename=.FALSE., mpi_io=mpi_io)
2610 : CALL cp_pw_to_cube(wf_r, unit_nr, "TOTAL DENSITY", &
2611 : particles=particles, zeff=zcharge, &
2612 : stride=section_get_ivals(dft_section, "PRINT%TOT_DENSITY_CUBE%STRIDE"), &
2613 : max_file_size_mb=section_get_rval(dft_section, "PRINT%TOT_DENSITY_CUBE%MAX_FILE_SIZE_MB"), &
2614 82 : mpi_io=mpi_io)
2615 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
2616 82 : "DFT%PRINT%TOT_DENSITY_CUBE", mpi_io=mpi_io)
2617 82 : CALL auxbas_pw_pool%give_back_pw(wf_r)
2618 : END IF
2619 :
2620 12969 : e_density_section = cube_or_openpmd(input, str_e_density_cubes, str_e_density_openpmd, logger)
2621 :
2622 : ! Write cube file with electron density
2623 12969 : IF (e_density_section%do_output) THEN
2624 : CALL section_vals_val_get(dft_section, &
2625 : keyword_name=e_density_section%concat_to_relative("%DENSITY_INCLUDE"), &
2626 150 : c_val=print_density)
2627 : print_density = TRIM(print_density)
2628 : ! For openPMD, refer to access modes instead of APPEND key
2629 150 : IF (e_density_section%grid_output == grid_output_cubes) THEN
2630 150 : append_cube = section_get_lval(input, e_density_section%concat_to_absolute("%APPEND"))
2631 : END IF
2632 150 : my_pos_cube = "REWIND"
2633 150 : IF (append_cube) THEN
2634 0 : my_pos_cube = "APPEND"
2635 : END IF
2636 : ! Write the info on core densities for the interface between cp2k and the XRD code
2637 : ! together with the valence density they are used to compute the form factor (Fourier transform)
2638 150 : IF (e_density_section%grid_output == grid_output_cubes) THEN
2639 150 : xrd_interface = section_get_lval(input, e_density_section%concat_to_absolute("%XRD_INTERFACE"))
2640 : ELSE
2641 : ! Unimplemented for openPMD, since this does not use the regular routines
2642 : xrd_interface = .FALSE.
2643 : END IF
2644 :
2645 150 : IF (xrd_interface) THEN
2646 : !cube file only contains soft density (GAPW)
2647 2 : IF (dft_control%qs_control%gapw) print_density = "SOFT_DENSITY"
2648 :
2649 2 : filename = "ELECTRON_DENSITY"
2650 : unit_nr = cp_print_key_unit_nr(logger, input, e_density_section%absolute_section_key, &
2651 : extension=".xrd", middle_name=TRIM(filename), &
2652 2 : file_position=my_pos_cube, log_filename=.FALSE.)
2653 2 : ngto = section_get_ival(input, e_density_section%concat_to_absolute("%NGAUSS"))
2654 2 : IF (output_unit > 0) THEN
2655 1 : INQUIRE (UNIT=unit_nr, NAME=filename)
2656 : WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") &
2657 1 : "The electron density (atomic part) is written to the file:", &
2658 2 : TRIM(filename)
2659 : END IF
2660 :
2661 2 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
2662 2 : nkind = SIZE(atomic_kind_set)
2663 2 : IF (unit_nr > 0) THEN
2664 1 : WRITE (unit_nr, *) "Atomic (core) densities"
2665 1 : WRITE (unit_nr, *) "Unit cell"
2666 1 : WRITE (unit_nr, FMT="(3F20.12)") cell%hmat(1, 1), cell%hmat(1, 2), cell%hmat(1, 3)
2667 1 : WRITE (unit_nr, FMT="(3F20.12)") cell%hmat(2, 1), cell%hmat(2, 2), cell%hmat(2, 3)
2668 1 : WRITE (unit_nr, FMT="(3F20.12)") cell%hmat(3, 1), cell%hmat(3, 2), cell%hmat(3, 3)
2669 1 : WRITE (unit_nr, *) "Atomic types"
2670 1 : WRITE (unit_nr, *) nkind
2671 : END IF
2672 : ! calculate atomic density and core density
2673 16 : ALLOCATE (ppdens(ngto, 2, nkind), aedens(ngto, 2, nkind), ccdens(ngto, 2, nkind))
2674 6 : DO ikind = 1, nkind
2675 4 : atomic_kind => atomic_kind_set(ikind)
2676 4 : qs_kind => qs_kind_set(ikind)
2677 4 : CALL get_atomic_kind(atomic_kind, name=name, element_symbol=element_symbol)
2678 : CALL calculate_atomic_density(ppdens(:, :, ikind), atomic_kind, qs_kind, ngto, &
2679 4 : iunit=output_unit, confine=.TRUE.)
2680 : CALL calculate_atomic_density(aedens(:, :, ikind), atomic_kind, qs_kind, ngto, &
2681 4 : iunit=output_unit, allelectron=.TRUE., confine=.TRUE.)
2682 52 : ccdens(:, 1, ikind) = aedens(:, 1, ikind)
2683 52 : ccdens(:, 2, ikind) = 0._dp
2684 : CALL project_function_a(ccdens(1:ngto, 2, ikind), ccdens(1:ngto, 1, ikind), &
2685 4 : ppdens(1:ngto, 2, ikind), ppdens(1:ngto, 1, ikind), 0)
2686 52 : ccdens(:, 2, ikind) = aedens(:, 2, ikind) - ccdens(:, 2, ikind)
2687 4 : IF (unit_nr > 0) THEN
2688 2 : WRITE (unit_nr, FMT="(I6,A10,A20)") ikind, TRIM(element_symbol), TRIM(name)
2689 2 : WRITE (unit_nr, FMT="(I6)") ngto
2690 2 : WRITE (unit_nr, *) " Total density"
2691 26 : WRITE (unit_nr, FMT="(2G24.12)") (aedens(i, 1, ikind), aedens(i, 2, ikind), i=1, ngto)
2692 2 : WRITE (unit_nr, *) " Core density"
2693 26 : WRITE (unit_nr, FMT="(2G24.12)") (ccdens(i, 1, ikind), ccdens(i, 2, ikind), i=1, ngto)
2694 : END IF
2695 6 : NULLIFY (atomic_kind)
2696 : END DO
2697 :
2698 2 : IF (dft_control%qs_control%gapw) THEN
2699 2 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set)
2700 :
2701 2 : IF (unit_nr > 0) THEN
2702 1 : WRITE (unit_nr, *) "Coordinates and GAPW density"
2703 : END IF
2704 2 : np = particles%n_els
2705 6 : DO iat = 1, np
2706 4 : CALL get_atomic_kind(particles%els(iat)%atomic_kind, kind_number=ikind)
2707 4 : CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid_atom)
2708 4 : rho_atom => rho_atom_set(iat)
2709 4 : IF (ASSOCIATED(rho_atom%rho_rad_h(1)%r_coef)) THEN
2710 2 : nr = SIZE(rho_atom%rho_rad_h(1)%r_coef, 1)
2711 2 : niso = SIZE(rho_atom%rho_rad_h(1)%r_coef, 2)
2712 : ELSE
2713 2 : nr = 0
2714 2 : niso = 0
2715 : END IF
2716 4 : CALL para_env%sum(nr)
2717 4 : CALL para_env%sum(niso)
2718 :
2719 16 : ALLOCATE (bfun(nr, niso))
2720 4 : bfun = 0._dp
2721 8 : DO ispin = 1, dft_control%nspins
2722 8 : IF (ASSOCIATED(rho_atom%rho_rad_h(1)%r_coef)) THEN
2723 920 : bfun(:, :) = bfun + rho_atom%rho_rad_h(ispin)%r_coef - rho_atom%rho_rad_s(ispin)%r_coef
2724 : END IF
2725 : END DO
2726 4 : CALL para_env%sum(bfun)
2727 52 : ccdens(:, 1, ikind) = ppdens(:, 1, ikind)
2728 52 : ccdens(:, 2, ikind) = 0._dp
2729 4 : IF (unit_nr > 0) THEN
2730 8 : WRITE (unit_nr, '(I10,I5,3f12.6)') iat, ikind, particles%els(iat)%r
2731 : END IF
2732 40 : DO iso = 1, niso
2733 36 : l = indso(1, iso)
2734 36 : CALL project_function_b(ccdens(:, 2, ikind), ccdens(:, 1, ikind), bfun(:, iso), grid_atom, l)
2735 40 : IF (unit_nr > 0) THEN
2736 18 : WRITE (unit_nr, FMT="(3I6)") iso, l, ngto
2737 234 : WRITE (unit_nr, FMT="(2G24.12)") (ccdens(i, 1, ikind), ccdens(i, 2, ikind), i=1, ngto)
2738 : END IF
2739 : END DO
2740 10 : DEALLOCATE (bfun)
2741 : END DO
2742 : ELSE
2743 0 : IF (unit_nr > 0) THEN
2744 0 : WRITE (unit_nr, *) "Coordinates"
2745 0 : np = particles%n_els
2746 0 : DO iat = 1, np
2747 0 : CALL get_atomic_kind(particles%els(iat)%atomic_kind, kind_number=ikind)
2748 0 : WRITE (unit_nr, '(I10,I5,3f12.6)') iat, ikind, particles%els(iat)%r
2749 : END DO
2750 : END IF
2751 : END IF
2752 :
2753 2 : DEALLOCATE (ppdens, aedens, ccdens)
2754 :
2755 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
2756 2 : e_density_section%absolute_section_key)
2757 :
2758 : END IF
2759 150 : IF (dft_control%qs_control%gapw .AND. print_density == "TOTAL_DENSITY") THEN
2760 : ! total density in g-space not implemented for k-points
2761 4 : CPASSERT(.NOT. do_kpoints)
2762 : ! Print total electronic density
2763 : CALL get_qs_env(qs_env=qs_env, &
2764 4 : pw_env=pw_env)
2765 : CALL pw_env_get(pw_env=pw_env, &
2766 : auxbas_pw_pool=auxbas_pw_pool, &
2767 4 : pw_pools=pw_pools)
2768 4 : CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
2769 4 : CALL pw_zero(rho_elec_rspace)
2770 4 : CALL auxbas_pw_pool%create_pw(pw=rho_elec_gspace)
2771 4 : CALL pw_zero(rho_elec_gspace)
2772 : CALL get_pw_grid_info(pw_grid=rho_elec_gspace%pw_grid, &
2773 : dr=dr, &
2774 4 : vol=volume)
2775 16 : q_max = SQRT(SUM((pi/dr(:))**2))
2776 : CALL calculate_rhotot_elec_gspace(qs_env=qs_env, &
2777 : auxbas_pw_pool=auxbas_pw_pool, &
2778 : rhotot_elec_gspace=rho_elec_gspace, &
2779 : q_max=q_max, &
2780 : rho_hard=rho_hard, &
2781 4 : rho_soft=rho_soft)
2782 4 : rho_total = rho_hard + rho_soft
2783 : CALL get_pw_grid_info(pw_grid=rho_elec_gspace%pw_grid, &
2784 4 : vol=volume)
2785 : ! rhotot pw coefficients are by default scaled by grid volume
2786 : ! need to undo this to get proper charge from printed cube
2787 4 : CALL pw_scale(rho_elec_gspace, 1.0_dp/volume)
2788 :
2789 4 : CALL pw_transfer(rho_elec_gspace, rho_elec_rspace)
2790 4 : rho_total_rspace = pw_integrate_function(rho_elec_rspace, isign=-1)
2791 4 : filename = "TOTAL_ELECTRON_DENSITY"
2792 4 : mpi_io = .TRUE.
2793 : unit_nr = e_density_section%print_key_unit_nr( &
2794 : logger, &
2795 : input, &
2796 : e_density_section%absolute_section_key, &
2797 : extension=".cube", &
2798 : middle_name=TRIM(filename), &
2799 : file_position=my_pos_cube, &
2800 : log_filename=.FALSE., &
2801 : mpi_io=mpi_io, &
2802 : fout=mpi_filename, &
2803 : openpmd_basename="dft-total-electron-density", &
2804 : openpmd_unit_dimension=openpmd_unit_dimension_density, &
2805 : openpmd_unit_si=openpmd_unit_si_density, &
2806 4 : sim_time=qs_env%sim_time)
2807 4 : IF (output_unit > 0) THEN
2808 2 : IF (.NOT. mpi_io) THEN
2809 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
2810 : ELSE
2811 2 : filename = mpi_filename
2812 : END IF
2813 : CALL print_density_output_message(output_unit, "The total electron density", &
2814 2 : e_density_section, filename)
2815 : WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") &
2816 2 : "q(max) [1/Angstrom] :", q_max/angstrom, &
2817 2 : "Soft electronic charge (G-space) :", rho_soft, &
2818 2 : "Hard electronic charge (G-space) :", rho_hard, &
2819 2 : "Total electronic charge (G-space):", rho_total, &
2820 4 : "Total electronic charge (R-space):", rho_total_rspace
2821 : END IF
2822 : CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "TOTAL ELECTRON DENSITY", &
2823 : particles=particles, zeff=zcharge, &
2824 4 : stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
2825 : CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
2826 4 : e_density_section%absolute_section_key, mpi_io=mpi_io)
2827 : ! Print total spin density for spin-polarized systems
2828 4 : IF (dft_control%nspins > 1) THEN
2829 2 : CALL pw_zero(rho_elec_gspace)
2830 2 : CALL pw_zero(rho_elec_rspace)
2831 : CALL calculate_rhotot_elec_gspace(qs_env=qs_env, &
2832 : auxbas_pw_pool=auxbas_pw_pool, &
2833 : rhotot_elec_gspace=rho_elec_gspace, &
2834 : q_max=q_max, &
2835 : rho_hard=rho_hard, &
2836 : rho_soft=rho_soft, &
2837 2 : fsign=-1.0_dp)
2838 2 : rho_total = rho_hard + rho_soft
2839 :
2840 : ! rhotot pw coefficients are by default scaled by grid volume
2841 : ! need to undo this to get proper charge from printed cube
2842 2 : CALL pw_scale(rho_elec_gspace, 1.0_dp/volume)
2843 :
2844 2 : CALL pw_transfer(rho_elec_gspace, rho_elec_rspace)
2845 2 : rho_total_rspace = pw_integrate_function(rho_elec_rspace, isign=-1)
2846 2 : filename = "TOTAL_SPIN_DENSITY"
2847 2 : mpi_io = .TRUE.
2848 : unit_nr = e_density_section%print_key_unit_nr( &
2849 : logger, &
2850 : input, &
2851 : e_density_section%absolute_section_key, &
2852 : extension=".cube", &
2853 : middle_name=TRIM(filename), &
2854 : file_position=my_pos_cube, &
2855 : log_filename=.FALSE., &
2856 : mpi_io=mpi_io, &
2857 : fout=mpi_filename, &
2858 : openpmd_basename="dft-total-spin-density", &
2859 : openpmd_unit_dimension=openpmd_unit_dimension_density, &
2860 : openpmd_unit_si=openpmd_unit_si_density, &
2861 2 : sim_time=qs_env%sim_time)
2862 2 : IF (output_unit > 0) THEN
2863 1 : IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
2864 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
2865 : ELSE
2866 1 : filename = mpi_filename
2867 : END IF
2868 : CALL print_density_output_message(output_unit, "The total spin density", &
2869 1 : e_density_section, filename)
2870 : WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") &
2871 1 : "q(max) [1/Angstrom] :", q_max/angstrom, &
2872 1 : "Soft part of the spin density (G-space):", rho_soft, &
2873 1 : "Hard part of the spin density (G-space):", rho_hard, &
2874 1 : "Total spin density (G-space) :", rho_total, &
2875 2 : "Total spin density (R-space) :", rho_total_rspace
2876 : END IF
2877 : CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "TOTAL SPIN DENSITY", &
2878 : particles=particles, zeff=zcharge, &
2879 2 : stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
2880 : CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
2881 2 : e_density_section%absolute_section_key, mpi_io=mpi_io)
2882 : END IF
2883 4 : CALL auxbas_pw_pool%give_back_pw(rho_elec_gspace)
2884 4 : CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
2885 :
2886 146 : ELSE IF (print_density == "SOFT_DENSITY" .OR. .NOT. dft_control%qs_control%gapw) THEN
2887 142 : IF (dft_control%nspins > 1) THEN
2888 : CALL get_qs_env(qs_env=qs_env, &
2889 48 : pw_env=pw_env)
2890 : CALL pw_env_get(pw_env=pw_env, &
2891 : auxbas_pw_pool=auxbas_pw_pool, &
2892 48 : pw_pools=pw_pools)
2893 48 : CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
2894 48 : CALL pw_copy(rho_r(1), rho_elec_rspace)
2895 48 : CALL pw_axpy(rho_r(2), rho_elec_rspace)
2896 48 : filename = "ELECTRON_DENSITY"
2897 48 : mpi_io = .TRUE.
2898 : unit_nr = e_density_section%print_key_unit_nr( &
2899 : logger, &
2900 : input, &
2901 : e_density_section%absolute_section_key, &
2902 : extension=".cube", &
2903 : middle_name=TRIM(filename), &
2904 : file_position=my_pos_cube, &
2905 : log_filename=.FALSE., &
2906 : mpi_io=mpi_io, &
2907 : fout=mpi_filename, &
2908 : openpmd_basename="dft-electron-density", &
2909 : openpmd_unit_dimension=openpmd_unit_dimension_density, &
2910 : openpmd_unit_si=openpmd_unit_si_density, &
2911 48 : sim_time=qs_env%sim_time)
2912 48 : IF (output_unit > 0) THEN
2913 24 : IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
2914 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
2915 : ELSE
2916 24 : filename = mpi_filename
2917 : END IF
2918 : CALL print_density_output_message(output_unit, "The sum of alpha and beta density", &
2919 24 : e_density_section, filename)
2920 : END IF
2921 : CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "SUM OF ALPHA AND BETA DENSITY", &
2922 : particles=particles, zeff=zcharge, stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), &
2923 48 : mpi_io=mpi_io)
2924 : CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
2925 48 : e_density_section%absolute_section_key, mpi_io=mpi_io)
2926 48 : CALL pw_copy(rho_r(1), rho_elec_rspace)
2927 48 : CALL pw_axpy(rho_r(2), rho_elec_rspace, alpha=-1.0_dp)
2928 48 : filename = "SPIN_DENSITY"
2929 48 : mpi_io = .TRUE.
2930 : unit_nr = e_density_section%print_key_unit_nr( &
2931 : logger, &
2932 : input, &
2933 : e_density_section%absolute_section_key, &
2934 : extension=".cube", &
2935 : middle_name=TRIM(filename), &
2936 : file_position=my_pos_cube, &
2937 : log_filename=.FALSE., &
2938 : mpi_io=mpi_io, &
2939 : fout=mpi_filename, &
2940 : openpmd_basename="dft-spin-density", &
2941 : openpmd_unit_dimension=openpmd_unit_dimension_density, &
2942 : openpmd_unit_si=openpmd_unit_si_density, &
2943 48 : sim_time=qs_env%sim_time)
2944 48 : IF (output_unit > 0) THEN
2945 24 : IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
2946 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
2947 : ELSE
2948 24 : filename = mpi_filename
2949 : END IF
2950 : CALL print_density_output_message(output_unit, "The spin density", &
2951 24 : e_density_section, filename)
2952 : END IF
2953 : CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "SPIN DENSITY", &
2954 : particles=particles, zeff=zcharge, &
2955 48 : stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
2956 : CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
2957 48 : e_density_section%absolute_section_key, mpi_io=mpi_io)
2958 48 : CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
2959 : ELSE
2960 94 : filename = "ELECTRON_DENSITY"
2961 94 : mpi_io = .TRUE.
2962 : unit_nr = e_density_section%print_key_unit_nr( &
2963 : logger, &
2964 : input, &
2965 : e_density_section%absolute_section_key, &
2966 : extension=".cube", &
2967 : middle_name=TRIM(filename), &
2968 : file_position=my_pos_cube, &
2969 : log_filename=.FALSE., &
2970 : mpi_io=mpi_io, &
2971 : fout=mpi_filename, &
2972 : openpmd_basename="dft-electron-density", &
2973 : openpmd_unit_dimension=openpmd_unit_dimension_density, &
2974 : openpmd_unit_si=openpmd_unit_si_density, &
2975 94 : sim_time=qs_env%sim_time)
2976 94 : IF (output_unit > 0) THEN
2977 47 : IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
2978 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
2979 : ELSE
2980 47 : filename = mpi_filename
2981 : END IF
2982 : CALL print_density_output_message(output_unit, "The electron density", &
2983 47 : e_density_section, filename)
2984 : END IF
2985 : CALL e_density_section%write_pw(rho_r(1), unit_nr, "ELECTRON DENSITY", &
2986 : particles=particles, zeff=zcharge, &
2987 94 : stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
2988 : CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
2989 94 : e_density_section%absolute_section_key, mpi_io=mpi_io)
2990 : END IF ! nspins
2991 :
2992 4 : ELSE IF (dft_control%qs_control%gapw .AND. print_density == "TOTAL_HARD_APPROX") THEN
2993 4 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho0_mpole=rho0_mpole, natom=natom)
2994 4 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
2995 4 : CALL auxbas_pw_pool%create_pw(rho_elec_rspace)
2996 :
2997 4 : NULLIFY (my_Q0)
2998 12 : ALLOCATE (my_Q0(natom))
2999 16 : my_Q0 = 0.0_dp
3000 :
3001 : ! (eta/pi)**3: normalization for 3d gaussian of form exp(-eta*r**2)
3002 4 : norm_factor = SQRT((rho0_mpole%zet0_h/pi)**3)
3003 :
3004 : ! store hard part of electronic density in array
3005 16 : DO iat = 1, natom
3006 34 : my_Q0(iat) = SUM(rho0_mpole%mp_rho(iat)%Q0(1:dft_control%nspins))*norm_factor
3007 : END DO
3008 : ! multiply coeff with gaussian and put on realspace grid
3009 : ! coeff is the gaussian prefactor, eta the gaussian exponent
3010 4 : CALL calculate_rho_resp_all(rho_elec_rspace, coeff=my_Q0, natom=natom, eta=rho0_mpole%zet0_h, qs_env=qs_env)
3011 4 : rho_hard = pw_integrate_function(rho_elec_rspace, isign=-1)
3012 :
3013 4 : rho_soft = 0.0_dp
3014 10 : DO ispin = 1, dft_control%nspins
3015 6 : CALL pw_axpy(rho_r(ispin), rho_elec_rspace)
3016 10 : rho_soft = rho_soft + pw_integrate_function(rho_r(ispin), isign=-1)
3017 : END DO
3018 :
3019 4 : rho_total_rspace = rho_soft + rho_hard
3020 :
3021 4 : filename = "ELECTRON_DENSITY"
3022 4 : mpi_io = .TRUE.
3023 : unit_nr = e_density_section%print_key_unit_nr( &
3024 : logger, &
3025 : input, &
3026 : e_density_section%absolute_section_key, &
3027 : extension=".cube", &
3028 : middle_name=TRIM(filename), &
3029 : file_position=my_pos_cube, &
3030 : log_filename=.FALSE., &
3031 : mpi_io=mpi_io, &
3032 : fout=mpi_filename, &
3033 : openpmd_basename="dft-electron-density", &
3034 : openpmd_unit_dimension=openpmd_unit_dimension_density, &
3035 : openpmd_unit_si=openpmd_unit_si_density, &
3036 4 : sim_time=qs_env%sim_time)
3037 4 : IF (output_unit > 0) THEN
3038 2 : IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
3039 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
3040 : ELSE
3041 2 : filename = mpi_filename
3042 : END IF
3043 : CALL print_density_output_message(output_unit, "The electron density", &
3044 2 : e_density_section, filename)
3045 : WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") &
3046 2 : "Soft electronic charge (R-space) :", rho_soft, &
3047 2 : "Hard electronic charge (R-space) :", rho_hard, &
3048 4 : "Total electronic charge (R-space):", rho_total_rspace
3049 : END IF
3050 : CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "ELECTRON DENSITY", &
3051 : particles=particles, zeff=zcharge, stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), &
3052 4 : mpi_io=mpi_io)
3053 : CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
3054 4 : e_density_section%absolute_section_key, mpi_io=mpi_io)
3055 :
3056 : !------------
3057 4 : IF (dft_control%nspins > 1) THEN
3058 8 : DO iat = 1, natom
3059 8 : my_Q0(iat) = (rho0_mpole%mp_rho(iat)%Q0(1) - rho0_mpole%mp_rho(iat)%Q0(2))*norm_factor
3060 : END DO
3061 2 : CALL pw_zero(rho_elec_rspace)
3062 2 : CALL calculate_rho_resp_all(rho_elec_rspace, coeff=my_Q0, natom=natom, eta=rho0_mpole%zet0_h, qs_env=qs_env)
3063 2 : rho_hard = pw_integrate_function(rho_elec_rspace, isign=-1)
3064 :
3065 2 : CALL pw_axpy(rho_r(1), rho_elec_rspace)
3066 2 : CALL pw_axpy(rho_r(2), rho_elec_rspace, alpha=-1.0_dp)
3067 : rho_soft = pw_integrate_function(rho_r(1), isign=-1) &
3068 2 : - pw_integrate_function(rho_r(2), isign=-1)
3069 :
3070 2 : rho_total_rspace = rho_soft + rho_hard
3071 :
3072 2 : filename = "SPIN_DENSITY"
3073 2 : mpi_io = .TRUE.
3074 : unit_nr = e_density_section%print_key_unit_nr( &
3075 : logger, &
3076 : input, &
3077 : e_density_section%absolute_section_key, &
3078 : extension=".cube", &
3079 : middle_name=TRIM(filename), &
3080 : file_position=my_pos_cube, &
3081 : log_filename=.FALSE., &
3082 : mpi_io=mpi_io, &
3083 : fout=mpi_filename, &
3084 : openpmd_basename="dft-spin-density", &
3085 : openpmd_unit_dimension=openpmd_unit_dimension_density, &
3086 : openpmd_unit_si=openpmd_unit_si_density, &
3087 2 : sim_time=qs_env%sim_time)
3088 2 : IF (output_unit > 0) THEN
3089 1 : IF (.NOT. mpi_io .AND. e_density_section%grid_output == grid_output_cubes) THEN
3090 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
3091 : ELSE
3092 1 : filename = mpi_filename
3093 : END IF
3094 : CALL print_density_output_message(output_unit, "The spin density", &
3095 1 : e_density_section, filename)
3096 : WRITE (UNIT=output_unit, FMT="(/,(T2,A,F20.10))") &
3097 1 : "Soft part of the spin density :", rho_soft, &
3098 1 : "Hard part of the spin density :", rho_hard, &
3099 2 : "Total spin density (R-space) :", rho_total_rspace
3100 : END IF
3101 : CALL e_density_section%write_pw(rho_elec_rspace, unit_nr, "SPIN DENSITY", &
3102 : particles=particles, zeff=zcharge, &
3103 2 : stride=section_get_ivals(dft_section, e_density_section%concat_to_relative("%STRIDE")), mpi_io=mpi_io)
3104 : CALL e_density_section%print_key_finished_output(unit_nr, logger, input, &
3105 2 : e_density_section%absolute_section_key, mpi_io=mpi_io)
3106 : END IF ! nspins
3107 4 : CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
3108 4 : DEALLOCATE (my_Q0)
3109 : END IF ! print_density
3110 : END IF ! print key
3111 :
3112 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
3113 12969 : dft_section, "PRINT%ENERGY_WINDOWS"), cp_p_file) .AND. .NOT. do_kpoints) THEN
3114 90 : CALL energy_windows(qs_env)
3115 : END IF
3116 :
3117 : ! Print the hartree potential
3118 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
3119 : "DFT%PRINT%V_HARTREE_CUBE"), cp_p_file)) THEN
3120 :
3121 : CALL get_qs_env(qs_env=qs_env, &
3122 : pw_env=pw_env, &
3123 114 : v_hartree_rspace=v_hartree_rspace)
3124 114 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3125 114 : CALL auxbas_pw_pool%create_pw(aux_r)
3126 :
3127 114 : append_cube = section_get_lval(input, "DFT%PRINT%V_HARTREE_CUBE%APPEND")
3128 114 : my_pos_cube = "REWIND"
3129 114 : IF (append_cube) THEN
3130 0 : my_pos_cube = "APPEND"
3131 : END IF
3132 114 : mpi_io = .TRUE.
3133 114 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
3134 114 : CALL pw_env_get(pw_env)
3135 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%V_HARTREE_CUBE", &
3136 114 : extension=".cube", middle_name="v_hartree", file_position=my_pos_cube, mpi_io=mpi_io)
3137 114 : udvol = 1.0_dp/v_hartree_rspace%pw_grid%dvol
3138 :
3139 114 : CALL pw_copy(v_hartree_rspace, aux_r)
3140 114 : CALL pw_scale(aux_r, udvol)
3141 :
3142 : CALL cp_pw_to_cube(aux_r, unit_nr, "HARTREE POTENTIAL", particles=particles, zeff=zcharge, &
3143 : stride=section_get_ivals(dft_section, "PRINT%V_HARTREE_CUBE%STRIDE"), &
3144 : max_file_size_mb=section_get_rval(dft_section, "PRINT%V_HARTREE_CUBE%MAX_FILE_SIZE_MB"), &
3145 114 : mpi_io=mpi_io)
3146 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
3147 114 : "DFT%PRINT%V_HARTREE_CUBE", mpi_io=mpi_io)
3148 :
3149 114 : CALL auxbas_pw_pool%give_back_pw(aux_r)
3150 : END IF
3151 :
3152 : ! Print the external potential
3153 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
3154 : "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE"), cp_p_file)) THEN
3155 86 : IF (dft_control%apply_external_potential) THEN
3156 4 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, vee=vee)
3157 4 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3158 4 : CALL auxbas_pw_pool%create_pw(aux_r)
3159 :
3160 4 : append_cube = section_get_lval(input, "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE%APPEND")
3161 4 : my_pos_cube = "REWIND"
3162 4 : IF (append_cube) THEN
3163 0 : my_pos_cube = "APPEND"
3164 : END IF
3165 4 : mpi_io = .TRUE.
3166 4 : CALL pw_env_get(pw_env)
3167 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE", &
3168 4 : extension=".cube", middle_name="ext_pot", file_position=my_pos_cube, mpi_io=mpi_io)
3169 :
3170 4 : CALL pw_copy(vee, aux_r)
3171 :
3172 : CALL cp_pw_to_cube(aux_r, unit_nr, "EXTERNAL POTENTIAL", particles=particles, zeff=zcharge, &
3173 : stride=section_get_ivals(dft_section, "PRINT%EXTERNAL_POTENTIAL_CUBE%STRIDE"), &
3174 : max_file_size_mb=section_get_rval(dft_section, "PRINT%EXTERNAL_POTENTIAL_CUBE%MAX_FILE_SIZE_MB"), &
3175 4 : mpi_io=mpi_io)
3176 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
3177 4 : "DFT%PRINT%EXTERNAL_POTENTIAL_CUBE", mpi_io=mpi_io)
3178 :
3179 4 : CALL auxbas_pw_pool%give_back_pw(aux_r)
3180 : END IF
3181 : END IF
3182 :
3183 : ! Print the Electrical Field Components
3184 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
3185 : "DFT%PRINT%EFIELD_CUBE"), cp_p_file)) THEN
3186 :
3187 82 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
3188 82 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3189 82 : CALL auxbas_pw_pool%create_pw(aux_r)
3190 82 : CALL auxbas_pw_pool%create_pw(aux_g)
3191 :
3192 82 : append_cube = section_get_lval(input, "DFT%PRINT%EFIELD_CUBE%APPEND")
3193 82 : my_pos_cube = "REWIND"
3194 82 : IF (append_cube) THEN
3195 0 : my_pos_cube = "APPEND"
3196 : END IF
3197 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, &
3198 82 : v_hartree_rspace=v_hartree_rspace)
3199 82 : CALL pw_env_get(pw_env)
3200 82 : udvol = 1.0_dp/v_hartree_rspace%pw_grid%dvol
3201 328 : DO id = 1, 3
3202 246 : mpi_io = .TRUE.
3203 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EFIELD_CUBE", &
3204 : extension=".cube", middle_name="efield_"//cdir(id), file_position=my_pos_cube, &
3205 246 : mpi_io=mpi_io)
3206 :
3207 246 : CALL pw_transfer(v_hartree_rspace, aux_g)
3208 246 : nd = 0
3209 246 : nd(id) = 1
3210 246 : CALL pw_derive(aux_g, nd)
3211 246 : CALL pw_transfer(aux_g, aux_r)
3212 246 : CALL pw_scale(aux_r, udvol)
3213 :
3214 : CALL cp_pw_to_cube(aux_r, unit_nr, "ELECTRIC FIELD", particles=particles, zeff=zcharge, &
3215 : stride=section_get_ivals(dft_section, "PRINT%EFIELD_CUBE%STRIDE"), &
3216 : max_file_size_mb=section_get_rval(dft_section, "PRINT%EFIELD_CUBE%MAX_FILE_SIZE_MB"), &
3217 246 : mpi_io=mpi_io)
3218 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
3219 328 : "DFT%PRINT%EFIELD_CUBE", mpi_io=mpi_io)
3220 : END DO
3221 :
3222 82 : CALL auxbas_pw_pool%give_back_pw(aux_r)
3223 82 : CALL auxbas_pw_pool%give_back_pw(aux_g)
3224 : END IF
3225 :
3226 : ! Write cube files from the local energy
3227 12969 : CALL qs_scf_post_local_energy(input, logger, qs_env)
3228 :
3229 : ! Write cube files from the local stress tensor
3230 12969 : CALL qs_scf_post_local_stress(input, logger, qs_env)
3231 :
3232 : ! Write cube files from the implicit Poisson solver
3233 12969 : CALL qs_scf_post_ps_implicit(input, logger, qs_env)
3234 :
3235 : ! post SCF finite-volume Kubo transport
3236 12969 : CALL qs_scf_post_kubo_transport(qs_env)
3237 :
3238 : ! post SCF Transport
3239 12969 : CALL qs_scf_post_transport(qs_env)
3240 :
3241 12969 : CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
3242 : ! Write the density matrices
3243 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
3244 : "DFT%PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN
3245 : iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/DENSITY", &
3246 4 : extension=".Log")
3247 4 : CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
3248 4 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
3249 4 : after = MIN(MAX(after, 1), 16)
3250 8 : DO ispin = 1, dft_control%nspins
3251 12 : DO img = 1, dft_control%nimages
3252 : CALL cp_dbcsr_write_sparse_matrix(rho_ao(ispin, img)%matrix, 4, after, qs_env, &
3253 8 : para_env, output_unit=iw, omit_headers=omit_headers)
3254 : END DO
3255 : END DO
3256 : CALL cp_print_key_finished_output(iw, logger, input, &
3257 4 : "DFT%PRINT%AO_MATRICES/DENSITY")
3258 : END IF
3259 :
3260 : ! Write the Kohn-Sham matrices
3261 : write_ks = BTEST(cp_print_key_should_output(logger%iter_info, input, &
3262 12969 : "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)
3263 : write_xc = BTEST(cp_print_key_should_output(logger%iter_info, input, &
3264 12969 : "DFT%PRINT%AO_MATRICES/MATRIX_VXC"), cp_p_file)
3265 : ! we need to update stuff before writing, potentially computing the matrix_vxc
3266 12969 : IF (write_ks .OR. write_xc) THEN
3267 4 : IF (write_xc) qs_env%requires_matrix_vxc = .TRUE.
3268 4 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
3269 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., &
3270 4 : just_energy=.FALSE.)
3271 4 : IF (write_xc) qs_env%requires_matrix_vxc = .FALSE.
3272 : END IF
3273 :
3274 : ! Write the Kohn-Sham matrices
3275 12969 : IF (write_ks) THEN
3276 : iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", &
3277 4 : extension=".Log")
3278 4 : CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=ks_rmpv)
3279 4 : CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
3280 4 : after = MIN(MAX(after, 1), 16)
3281 8 : DO ispin = 1, dft_control%nspins
3282 12 : DO img = 1, dft_control%nimages
3283 : CALL cp_dbcsr_write_sparse_matrix(ks_rmpv(ispin, img)%matrix, 4, after, qs_env, &
3284 8 : para_env, output_unit=iw, omit_headers=omit_headers)
3285 : END DO
3286 : END DO
3287 : CALL cp_print_key_finished_output(iw, logger, input, &
3288 4 : "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
3289 : END IF
3290 :
3291 : ! write csr matrices
3292 : ! matrices in terms of the PAO basis will be taken care of in pao_post_scf.
3293 12969 : IF (.NOT. dft_control%qs_control%pao) THEN
3294 12457 : CALL write_ks_matrix_csr(qs_env, input)
3295 12457 : CALL write_s_matrix_csr(qs_env, input)
3296 12457 : CALL write_hcore_matrix_csr(qs_env, input)
3297 12457 : CALL write_p_matrix_csr(qs_env, input)
3298 : END IF
3299 :
3300 : ! write adjacency matrix
3301 12969 : CALL write_adjacency_matrix(qs_env, input)
3302 :
3303 : ! Write the xc matrix
3304 12969 : IF (write_xc) THEN
3305 0 : CALL get_qs_env(qs_env=qs_env, matrix_vxc_kp=matrix_vxc)
3306 0 : CPASSERT(ASSOCIATED(matrix_vxc))
3307 : iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/MATRIX_VXC", &
3308 0 : extension=".Log")
3309 0 : CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
3310 0 : after = MIN(MAX(after, 1), 16)
3311 0 : DO ispin = 1, dft_control%nspins
3312 0 : DO img = 1, dft_control%nimages
3313 : CALL cp_dbcsr_write_sparse_matrix(matrix_vxc(ispin, img)%matrix, 4, after, qs_env, &
3314 0 : para_env, output_unit=iw, omit_headers=omit_headers)
3315 : END DO
3316 : END DO
3317 : CALL cp_print_key_finished_output(iw, logger, input, &
3318 0 : "DFT%PRINT%AO_MATRICES/MATRIX_VXC")
3319 : END IF
3320 :
3321 : ! Write the [H,r] commutator matrices
3322 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
3323 : "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR"), cp_p_file)) THEN
3324 : iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR", &
3325 0 : extension=".Log")
3326 0 : CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
3327 0 : NULLIFY (matrix_hr)
3328 0 : CALL build_com_hr_matrix(qs_env, matrix_hr)
3329 0 : after = MIN(MAX(after, 1), 16)
3330 0 : DO img = 1, 3
3331 : CALL cp_dbcsr_write_sparse_matrix(matrix_hr(img)%matrix, 4, after, qs_env, &
3332 0 : para_env, output_unit=iw, omit_headers=omit_headers)
3333 : END DO
3334 0 : CALL dbcsr_deallocate_matrix_set(matrix_hr)
3335 : CALL cp_print_key_finished_output(iw, logger, input, &
3336 0 : "DFT%PRINT%AO_MATRICES/COMMUTATOR_HR")
3337 : END IF
3338 :
3339 : ! Compute the Mulliken charges
3340 12969 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MULLIKEN")
3341 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
3342 5302 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MULLIKEN", extension=".mulliken", log_filename=.FALSE.)
3343 5302 : print_level = 1
3344 5302 : CALL section_vals_val_get(print_key, "PRINT_GOP", l_val=print_it)
3345 5302 : IF (print_it) print_level = 2
3346 5302 : CALL section_vals_val_get(print_key, "PRINT_ALL", l_val=print_it)
3347 5302 : IF (print_it) print_level = 3
3348 5302 : CALL mulliken_population_analysis(qs_env, unit_nr, print_level)
3349 5302 : CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MULLIKEN")
3350 : END IF
3351 :
3352 : ! Compute the Hirshfeld charges
3353 12969 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%HIRSHFELD")
3354 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
3355 : ! we check if real space density is available
3356 5378 : NULLIFY (rho)
3357 5378 : CALL get_qs_env(qs_env=qs_env, rho=rho)
3358 5378 : CALL qs_rho_get(rho, rho_r_valid=rho_r_valid)
3359 5378 : IF (rho_r_valid) THEN
3360 5304 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%HIRSHFELD", extension=".hirshfeld", log_filename=.FALSE.)
3361 5304 : CALL hirshfeld_charges(qs_env, print_key, unit_nr)
3362 5304 : CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%HIRSHFELD")
3363 : END IF
3364 : END IF
3365 :
3366 : ! Compute EEQ charges
3367 12969 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%EEQ_CHARGES")
3368 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
3369 30 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EEQ_CHARGES", extension=".eeq", log_filename=.FALSE.)
3370 30 : print_level = 1
3371 30 : CALL eeq_print(qs_env, unit_nr, print_level, ext=.FALSE.)
3372 30 : CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MULLIKEN")
3373 : END IF
3374 :
3375 : ! Do a Voronoi Integration or write a compressed BQB File
3376 12969 : print_key_voro => section_vals_get_subs_vals(input, "DFT%PRINT%VORONOI")
3377 12969 : print_key_bqb => section_vals_get_subs_vals(input, "DFT%PRINT%E_DENSITY_BQB")
3378 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key_voro), cp_p_file)) THEN
3379 24 : should_print_voro = 1
3380 : ELSE
3381 12945 : should_print_voro = 0
3382 : END IF
3383 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key_bqb), cp_p_file)) THEN
3384 2 : should_print_bqb = 1
3385 : ELSE
3386 12967 : should_print_bqb = 0
3387 : END IF
3388 12969 : IF ((should_print_voro /= 0) .OR. (should_print_bqb /= 0)) THEN
3389 :
3390 : ! we check if real space density is available
3391 26 : NULLIFY (rho)
3392 26 : CALL get_qs_env(qs_env=qs_env, rho=rho)
3393 26 : CALL qs_rho_get(rho, rho_r_valid=rho_r_valid)
3394 26 : IF (rho_r_valid) THEN
3395 :
3396 26 : IF (dft_control%nspins > 1) THEN
3397 : CALL get_qs_env(qs_env=qs_env, &
3398 0 : pw_env=pw_env)
3399 : CALL pw_env_get(pw_env=pw_env, &
3400 : auxbas_pw_pool=auxbas_pw_pool, &
3401 0 : pw_pools=pw_pools)
3402 0 : NULLIFY (mb_rho)
3403 0 : ALLOCATE (mb_rho)
3404 0 : CALL auxbas_pw_pool%create_pw(pw=mb_rho)
3405 0 : CALL pw_copy(rho_r(1), mb_rho)
3406 0 : CALL pw_axpy(rho_r(2), mb_rho)
3407 : !CALL voronoi_analysis(qs_env, rho_elec_rspace, print_key, unit_nr)
3408 : ELSE
3409 26 : mb_rho => rho_r(1)
3410 : !CALL voronoi_analysis( qs_env, rho_r(1), print_key, unit_nr )
3411 : END IF ! nspins
3412 :
3413 26 : IF (should_print_voro /= 0) THEN
3414 24 : CALL section_vals_val_get(print_key_voro, "OUTPUT_TEXT", l_val=voro_print_txt)
3415 24 : IF (voro_print_txt) THEN
3416 24 : append_voro = section_get_lval(input, "DFT%PRINT%VORONOI%APPEND")
3417 24 : my_pos_voro = "REWIND"
3418 24 : IF (append_voro) THEN
3419 0 : my_pos_voro = "APPEND"
3420 : END IF
3421 : unit_nr_voro = cp_print_key_unit_nr(logger, input, "DFT%PRINT%VORONOI", extension=".voronoi", &
3422 24 : file_position=my_pos_voro, log_filename=.FALSE.)
3423 : ELSE
3424 0 : unit_nr_voro = 0
3425 : END IF
3426 : ELSE
3427 2 : unit_nr_voro = 0
3428 : END IF
3429 :
3430 : CALL entry_voronoi_or_bqb(should_print_voro, should_print_bqb, print_key_voro, print_key_bqb, &
3431 26 : unit_nr_voro, qs_env, mb_rho)
3432 :
3433 26 : IF (dft_control%nspins > 1) THEN
3434 0 : CALL auxbas_pw_pool%give_back_pw(mb_rho)
3435 0 : DEALLOCATE (mb_rho)
3436 : END IF
3437 :
3438 26 : IF (unit_nr_voro > 0) THEN
3439 12 : CALL cp_print_key_finished_output(unit_nr_voro, logger, input, "DFT%PRINT%VORONOI")
3440 : END IF
3441 :
3442 : END IF
3443 : END IF
3444 :
3445 : ! MAO analysis
3446 12969 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MAO_ANALYSIS")
3447 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
3448 38 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MAO_ANALYSIS", extension=".mao", log_filename=.FALSE.)
3449 38 : CALL mao_analysis(qs_env, print_key, unit_nr)
3450 38 : CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MAO_ANALYSIS")
3451 : END IF
3452 :
3453 : ! MINBAS analysis
3454 12969 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MINBAS_ANALYSIS")
3455 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
3456 28 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MINBAS_ANALYSIS", extension=".mao", log_filename=.FALSE.)
3457 28 : CALL minbas_analysis(qs_env, print_key, unit_nr)
3458 28 : CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MINBAS_ANALYSIS")
3459 : END IF
3460 :
3461 : ! IAO analysis
3462 12969 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%IAO_ANALYSIS")
3463 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
3464 34 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IAO_ANALYSIS", extension=".iao", log_filename=.FALSE.)
3465 34 : CALL iao_read_input(iao_env, print_key, cell)
3466 34 : IF (particle_set(1)%fragment_index /= 0) iao_env%do_fragments = .TRUE.
3467 34 : IF (iao_env%do_iao) THEN
3468 6 : CALL iao_wfn_analysis(qs_env, iao_env, unit_nr)
3469 : END IF
3470 34 : CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%IAO_ANALYSIS")
3471 : END IF
3472 :
3473 : ! Energy Decomposition Analysis
3474 12969 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%ENERGY_DECOMPOSITION_ANALYSIS")
3475 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
3476 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%ENERGY_DECOMPOSITION_ANALYSIS", &
3477 58 : extension=".mao", log_filename=.FALSE.)
3478 58 : CALL edmf_analysis(qs_env, print_key, unit_nr)
3479 58 : CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%ENERGY_DECOMPOSITION_ANALYSIS")
3480 : END IF
3481 :
3482 : ! Print the density in the RI-HFX basis
3483 12969 : hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
3484 12969 : CALL section_vals_get(hfx_section, explicit=do_hfx)
3485 12969 : CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
3486 12969 : IF (do_hfx) THEN
3487 5042 : DO i = 1, n_rep_hf
3488 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)
3489 : END DO
3490 : END IF
3491 :
3492 12969 : DEALLOCATE (zcharge)
3493 :
3494 12969 : CALL timestop(handle)
3495 :
3496 51876 : END SUBROUTINE write_mo_free_results
3497 :
3498 : ! **************************************************************************************************
3499 : !> \brief Calculates Hirshfeld charges
3500 : !> \param qs_env the qs_env where to calculate the charges
3501 : !> \param input_section the input section for Hirshfeld charges
3502 : !> \param unit_nr the output unit number
3503 : ! **************************************************************************************************
3504 5304 : SUBROUTINE hirshfeld_charges(qs_env, input_section, unit_nr)
3505 : TYPE(qs_environment_type), POINTER :: qs_env
3506 : TYPE(section_vals_type), POINTER :: input_section
3507 : INTEGER, INTENT(IN) :: unit_nr
3508 :
3509 : INTEGER :: i, iat, ikind, natom, nkind, nspin, &
3510 : radius_type, refc, shapef
3511 5304 : INTEGER, DIMENSION(:), POINTER :: atom_list
3512 : LOGICAL :: do_radius, do_sc, paw_atom
3513 : REAL(KIND=dp) :: zeff
3514 5304 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
3515 5304 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges
3516 5304 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3517 : TYPE(atomic_kind_type), POINTER :: atomic_kind
3518 5304 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
3519 : TYPE(dft_control_type), POINTER :: dft_control
3520 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
3521 : TYPE(mp_para_env_type), POINTER :: para_env
3522 5304 : TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho
3523 5304 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3524 5304 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3525 : TYPE(qs_rho_type), POINTER :: rho
3526 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
3527 :
3528 5304 : NULLIFY (hirshfeld_env)
3529 5304 : NULLIFY (radii)
3530 5304 : CALL create_hirshfeld_type(hirshfeld_env)
3531 : !
3532 5304 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
3533 15912 : ALLOCATE (hirshfeld_env%charges(natom))
3534 : ! input options
3535 5304 : CALL section_vals_val_get(input_section, "SELF_CONSISTENT", l_val=do_sc)
3536 5304 : CALL section_vals_val_get(input_section, "USER_RADIUS", l_val=do_radius)
3537 5304 : CALL section_vals_val_get(input_section, "SHAPE_FUNCTION", i_val=shapef)
3538 5304 : CALL section_vals_val_get(input_section, "REFERENCE_CHARGE", i_val=refc)
3539 5304 : IF (do_radius) THEN
3540 0 : radius_type = radius_user
3541 0 : CALL section_vals_val_get(input_section, "ATOMIC_RADII", r_vals=radii)
3542 0 : IF (.NOT. SIZE(radii) == nkind) THEN
3543 : CALL cp_abort(__LOCATION__, &
3544 : "Length of keyword HIRSHFELD\ATOMIC_RADII does not "// &
3545 0 : "match number of atomic kinds in the input coordinate file.")
3546 : END IF
3547 : ELSE
3548 5304 : radius_type = radius_covalent
3549 : END IF
3550 : CALL set_hirshfeld_info(hirshfeld_env, shape_function_type=shapef, &
3551 : iterative=do_sc, ref_charge=refc, &
3552 5304 : radius_type=radius_type)
3553 : ! shape function
3554 5304 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
3555 : CALL create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, &
3556 5304 : radii_list=radii)
3557 : ! reference charges
3558 5304 : CALL get_qs_env(qs_env, rho=rho)
3559 5304 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
3560 5304 : nspin = SIZE(matrix_p, 1)
3561 21216 : ALLOCATE (charges(natom, nspin))
3562 5292 : SELECT CASE (refc)
3563 : CASE (ref_charge_atomic)
3564 14528 : DO ikind = 1, nkind
3565 9236 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
3566 9236 : atomic_kind => atomic_kind_set(ikind)
3567 9236 : CALL get_atomic_kind(atomic_kind, atom_list=atom_list)
3568 46076 : DO iat = 1, SIZE(atom_list)
3569 22312 : i = atom_list(iat)
3570 31548 : hirshfeld_env%charges(i) = zeff
3571 : END DO
3572 : END DO
3573 : CASE (ref_charge_mulliken)
3574 12 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env)
3575 12 : CALL mulliken_charges(matrix_p, matrix_s, para_env, charges)
3576 48 : DO iat = 1, natom
3577 108 : hirshfeld_env%charges(iat) = SUM(charges(iat, :))
3578 : END DO
3579 : CASE DEFAULT
3580 5304 : CPABORT("Unknown type of reference charge for Hirshfeld partitioning.")
3581 : END SELECT
3582 : !
3583 36736 : charges = 0.0_dp
3584 5304 : IF (hirshfeld_env%iterative) THEN
3585 : ! Hirshfeld-I charges
3586 22 : CALL comp_hirshfeld_i_charges(qs_env, hirshfeld_env, charges, unit_nr)
3587 : ELSE
3588 : ! Hirshfeld charges
3589 5282 : CALL comp_hirshfeld_charges(qs_env, hirshfeld_env, charges)
3590 : END IF
3591 5304 : CALL get_qs_env(qs_env, particle_set=particle_set, dft_control=dft_control)
3592 5304 : IF (dft_control%qs_control%gapw) THEN
3593 : ! GAPW: add core charges (rho_hard - rho_soft)
3594 866 : CALL get_qs_env(qs_env, rho0_mpole=rho0_mpole)
3595 866 : CALL get_rho0_mpole(rho0_mpole, mp_rho=mp_rho)
3596 3724 : DO iat = 1, natom
3597 2858 : atomic_kind => particle_set(iat)%atomic_kind
3598 2858 : CALL get_atomic_kind(atomic_kind, kind_number=ikind)
3599 2858 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
3600 3724 : IF (paw_atom) THEN
3601 5514 : charges(iat, 1:nspin) = charges(iat, 1:nspin) + mp_rho(iat)%q0(1:nspin)
3602 : END IF
3603 : END DO
3604 : END IF
3605 : !
3606 5304 : IF (unit_nr > 0) THEN
3607 : CALL write_hirshfeld_charges(charges, hirshfeld_env, particle_set, &
3608 2667 : qs_kind_set, unit_nr)
3609 : END IF
3610 : ! Save the charges to the results under the tag [HIRSHFELD-CHARGES]
3611 5304 : CALL save_hirshfeld_charges(charges, particle_set, qs_kind_set, qs_env)
3612 : !
3613 5304 : CALL release_hirshfeld_type(hirshfeld_env)
3614 5304 : DEALLOCATE (charges)
3615 :
3616 10608 : END SUBROUTINE hirshfeld_charges
3617 :
3618 : ! **************************************************************************************************
3619 : !> \brief ...
3620 : !> \param ca ...
3621 : !> \param a ...
3622 : !> \param cb ...
3623 : !> \param b ...
3624 : !> \param l ...
3625 : ! **************************************************************************************************
3626 4 : SUBROUTINE project_function_a(ca, a, cb, b, l)
3627 : ! project function cb on ca
3628 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ca
3629 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, cb, b
3630 : INTEGER, INTENT(IN) :: l
3631 :
3632 : INTEGER :: info, n
3633 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
3634 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smat, tmat, v
3635 :
3636 4 : n = SIZE(ca)
3637 40 : ALLOCATE (smat(n, n), tmat(n, n), v(n, 1), ipiv(n))
3638 :
3639 4 : CALL sg_overlap(smat, l, a, a)
3640 4 : CALL sg_overlap(tmat, l, a, b)
3641 1252 : v(:, 1) = MATMUL(tmat, cb)
3642 4 : CALL dgesv(n, 1, smat, n, ipiv, v, n, info)
3643 4 : CPASSERT(info == 0)
3644 52 : ca(:) = v(:, 1)
3645 :
3646 4 : DEALLOCATE (smat, tmat, v, ipiv)
3647 :
3648 4 : END SUBROUTINE project_function_a
3649 :
3650 : ! **************************************************************************************************
3651 : !> \brief ...
3652 : !> \param ca ...
3653 : !> \param a ...
3654 : !> \param bfun ...
3655 : !> \param grid_atom ...
3656 : !> \param l ...
3657 : ! **************************************************************************************************
3658 36 : SUBROUTINE project_function_b(ca, a, bfun, grid_atom, l)
3659 : ! project function f on ca
3660 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ca
3661 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, bfun
3662 : TYPE(grid_atom_type), POINTER :: grid_atom
3663 : INTEGER, INTENT(IN) :: l
3664 :
3665 : INTEGER :: i, info, n, nr
3666 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
3667 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: afun
3668 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smat, v
3669 :
3670 36 : n = SIZE(ca)
3671 36 : nr = grid_atom%nr
3672 360 : ALLOCATE (smat(n, n), v(n, 1), ipiv(n), afun(nr))
3673 :
3674 36 : CALL sg_overlap(smat, l, a, a)
3675 468 : DO i = 1, n
3676 22032 : afun(:) = grid_atom%rad(:)**l*EXP(-a(i)*grid_atom%rad2(:))
3677 22068 : v(i, 1) = SUM(afun(:)*bfun(:)*grid_atom%wr(:))
3678 : END DO
3679 36 : CALL dgesv(n, 1, smat, n, ipiv, v, n, info)
3680 36 : CPASSERT(info == 0)
3681 468 : ca(:) = v(:, 1)
3682 :
3683 36 : DEALLOCATE (smat, v, ipiv, afun)
3684 :
3685 36 : END SUBROUTINE project_function_b
3686 :
3687 : ! **************************************************************************************************
3688 : !> \brief Performs printing of cube files from local energy
3689 : !> \param input input
3690 : !> \param logger the logger
3691 : !> \param qs_env the qs_env in which the qs_env lives
3692 : !> \par History
3693 : !> 07.2019 created
3694 : !> \author JGH
3695 : ! **************************************************************************************************
3696 12969 : SUBROUTINE qs_scf_post_local_energy(input, logger, qs_env)
3697 : TYPE(section_vals_type), POINTER :: input
3698 : TYPE(cp_logger_type), POINTER :: logger
3699 : TYPE(qs_environment_type), POINTER :: qs_env
3700 :
3701 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_local_energy'
3702 :
3703 : CHARACTER(LEN=default_path_length) :: filename, my_pos_cube
3704 : INTEGER :: handle, io_unit, natom, unit_nr
3705 : LOGICAL :: append_cube, gapw, gapw_xc, mpi_io
3706 : TYPE(dft_control_type), POINTER :: dft_control
3707 : TYPE(particle_list_type), POINTER :: particles
3708 : TYPE(pw_env_type), POINTER :: pw_env
3709 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3710 : TYPE(pw_r3d_rs_type) :: eden
3711 : TYPE(qs_subsys_type), POINTER :: subsys
3712 : TYPE(section_vals_type), POINTER :: dft_section
3713 :
3714 12969 : CALL timeset(routineN, handle)
3715 12969 : io_unit = cp_logger_get_default_io_unit(logger)
3716 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
3717 : "DFT%PRINT%LOCAL_ENERGY_CUBE"), cp_p_file)) THEN
3718 34 : dft_section => section_vals_get_subs_vals(input, "DFT")
3719 34 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, natom=natom)
3720 34 : gapw = dft_control%qs_control%gapw
3721 34 : gapw_xc = dft_control%qs_control%gapw_xc
3722 34 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
3723 34 : CALL qs_subsys_get(subsys, particles=particles)
3724 34 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3725 34 : CALL auxbas_pw_pool%create_pw(eden)
3726 : !
3727 34 : CALL qs_local_energy(qs_env, eden)
3728 : !
3729 34 : append_cube = section_get_lval(input, "DFT%PRINT%LOCAL_ENERGY_CUBE%APPEND")
3730 34 : IF (append_cube) THEN
3731 0 : my_pos_cube = "APPEND"
3732 : ELSE
3733 34 : my_pos_cube = "REWIND"
3734 : END IF
3735 34 : mpi_io = .TRUE.
3736 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%LOCAL_ENERGY_CUBE", &
3737 : extension=".cube", middle_name="local_energy", &
3738 34 : file_position=my_pos_cube, mpi_io=mpi_io)
3739 : CALL cp_pw_to_cube(eden, unit_nr, "LOCAL ENERGY", particles=particles, &
3740 : stride=section_get_ivals(dft_section, "PRINT%LOCAL_ENERGY_CUBE%STRIDE"), &
3741 : max_file_size_mb=section_get_rval(dft_section, "PRINT%LOCAL_ENERGY_CUBE%MAX_FILE_SIZE_MB"), &
3742 34 : mpi_io=mpi_io)
3743 34 : IF (io_unit > 0) THEN
3744 17 : INQUIRE (UNIT=unit_nr, NAME=filename)
3745 17 : IF (gapw .OR. gapw_xc) THEN
3746 : WRITE (UNIT=io_unit, FMT="(/,T3,A,A)") &
3747 1 : "The soft part of the local energy is written to the file: ", TRIM(ADJUSTL(filename))
3748 : ELSE
3749 : WRITE (UNIT=io_unit, FMT="(/,T3,A,A)") &
3750 16 : "The local energy is written to the file: ", TRIM(ADJUSTL(filename))
3751 : END IF
3752 : END IF
3753 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
3754 34 : "DFT%PRINT%LOCAL_ENERGY_CUBE", mpi_io=mpi_io)
3755 : !
3756 34 : CALL auxbas_pw_pool%give_back_pw(eden)
3757 : END IF
3758 12969 : CALL timestop(handle)
3759 :
3760 12969 : END SUBROUTINE qs_scf_post_local_energy
3761 :
3762 : ! **************************************************************************************************
3763 : !> \brief Performs printing of cube files from local energy
3764 : !> \param input input
3765 : !> \param logger the logger
3766 : !> \param qs_env the qs_env in which the qs_env lives
3767 : !> \par History
3768 : !> 07.2019 created
3769 : !> \author JGH
3770 : ! **************************************************************************************************
3771 12969 : SUBROUTINE qs_scf_post_local_stress(input, logger, qs_env)
3772 : TYPE(section_vals_type), POINTER :: input
3773 : TYPE(cp_logger_type), POINTER :: logger
3774 : TYPE(qs_environment_type), POINTER :: qs_env
3775 :
3776 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_local_stress'
3777 :
3778 : CHARACTER(LEN=default_path_length) :: filename, my_pos_cube
3779 : INTEGER :: handle, io_unit, natom, unit_nr
3780 : LOGICAL :: append_cube, gapw, gapw_xc, mpi_io
3781 : REAL(KIND=dp) :: beta
3782 : TYPE(dft_control_type), POINTER :: dft_control
3783 : TYPE(particle_list_type), POINTER :: particles
3784 : TYPE(pw_env_type), POINTER :: pw_env
3785 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3786 : TYPE(pw_r3d_rs_type) :: stress
3787 : TYPE(qs_subsys_type), POINTER :: subsys
3788 : TYPE(section_vals_type), POINTER :: dft_section
3789 :
3790 12969 : CALL timeset(routineN, handle)
3791 12969 : io_unit = cp_logger_get_default_io_unit(logger)
3792 12969 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
3793 : "DFT%PRINT%LOCAL_STRESS_CUBE"), cp_p_file)) THEN
3794 : CALL cp_warn(__LOCATION__, &
3795 30 : "LOCAL_STRESS_CUBE uses the existing experimental local stress implementation")
3796 30 : dft_section => section_vals_get_subs_vals(input, "DFT")
3797 30 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, natom=natom)
3798 30 : gapw = dft_control%qs_control%gapw
3799 30 : gapw_xc = dft_control%qs_control%gapw_xc
3800 30 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
3801 30 : CALL qs_subsys_get(subsys, particles=particles)
3802 30 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3803 30 : CALL auxbas_pw_pool%create_pw(stress)
3804 : !
3805 : ! use beta=0: kinetic energy density in symmetric form
3806 30 : beta = 0.0_dp
3807 30 : CALL qs_local_stress(qs_env, beta=beta)
3808 : !
3809 30 : append_cube = section_get_lval(input, "DFT%PRINT%LOCAL_STRESS_CUBE%APPEND")
3810 30 : IF (append_cube) THEN
3811 0 : my_pos_cube = "APPEND"
3812 : ELSE
3813 30 : my_pos_cube = "REWIND"
3814 : END IF
3815 30 : mpi_io = .TRUE.
3816 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%LOCAL_STRESS_CUBE", &
3817 : extension=".cube", middle_name="local_stress", &
3818 30 : file_position=my_pos_cube, mpi_io=mpi_io)
3819 : CALL cp_pw_to_cube(stress, unit_nr, "LOCAL STRESS", particles=particles, &
3820 : stride=section_get_ivals(dft_section, "PRINT%LOCAL_STRESS_CUBE%STRIDE"), &
3821 : max_file_size_mb=section_get_rval(dft_section, "PRINT%LOCAL_STRESS_CUBE%MAX_FILE_SIZE_MB"), &
3822 30 : mpi_io=mpi_io)
3823 30 : IF (io_unit > 0) THEN
3824 15 : INQUIRE (UNIT=unit_nr, NAME=filename)
3825 15 : WRITE (UNIT=io_unit, FMT="(/,T3,A)") "Write 1/3*Tr(sigma) to cube file"
3826 15 : IF (gapw .OR. gapw_xc) THEN
3827 : WRITE (UNIT=io_unit, FMT="(T3,A,A)") &
3828 0 : "The soft part of the local stress is written to the file: ", TRIM(ADJUSTL(filename))
3829 : ELSE
3830 : WRITE (UNIT=io_unit, FMT="(T3,A,A)") &
3831 15 : "The local stress is written to the file: ", TRIM(ADJUSTL(filename))
3832 : END IF
3833 : END IF
3834 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
3835 30 : "DFT%PRINT%LOCAL_STRESS_CUBE", mpi_io=mpi_io)
3836 : !
3837 30 : CALL auxbas_pw_pool%give_back_pw(stress)
3838 : END IF
3839 :
3840 12969 : CALL timestop(handle)
3841 :
3842 12969 : END SUBROUTINE qs_scf_post_local_stress
3843 :
3844 : ! **************************************************************************************************
3845 : !> \brief Performs printing of cube files related to the implicit Poisson solver
3846 : !> \param input input
3847 : !> \param logger the logger
3848 : !> \param qs_env the qs_env in which the qs_env lives
3849 : !> \par History
3850 : !> 03.2016 refactored from write_mo_free_results [Hossein Bani-Hashemian]
3851 : !> \author Mohammad Hossein Bani-Hashemian
3852 : ! **************************************************************************************************
3853 12969 : SUBROUTINE qs_scf_post_ps_implicit(input, logger, qs_env)
3854 : TYPE(section_vals_type), POINTER :: input
3855 : TYPE(cp_logger_type), POINTER :: logger
3856 : TYPE(qs_environment_type), POINTER :: qs_env
3857 :
3858 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_ps_implicit'
3859 :
3860 : CHARACTER(LEN=default_path_length) :: filename, my_pos_cube
3861 : INTEGER :: boundary_condition, handle, i, j, &
3862 : n_cstr, n_tiles, unit_nr
3863 : LOGICAL :: append_cube, do_cstr_charge_cube, do_dielectric_cube, do_dirichlet_bc_cube, &
3864 : has_dirichlet_bc, has_implicit_ps, mpi_io, tile_cubes
3865 : TYPE(particle_list_type), POINTER :: particles
3866 : TYPE(pw_env_type), POINTER :: pw_env
3867 : TYPE(pw_poisson_type), POINTER :: poisson_env
3868 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3869 : TYPE(pw_r3d_rs_type) :: aux_r
3870 : TYPE(pw_r3d_rs_type), POINTER :: dirichlet_tile
3871 : TYPE(qs_subsys_type), POINTER :: subsys
3872 : TYPE(section_vals_type), POINTER :: dft_section
3873 :
3874 12969 : CALL timeset(routineN, handle)
3875 :
3876 12969 : NULLIFY (pw_env, auxbas_pw_pool, dft_section, particles)
3877 :
3878 12969 : dft_section => section_vals_get_subs_vals(input, "DFT")
3879 12969 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
3880 12969 : CALL qs_subsys_get(subsys, particles=particles)
3881 12969 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3882 :
3883 12969 : has_implicit_ps = .FALSE.
3884 12969 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
3885 12969 : IF (pw_env%poisson_env%parameters%solver == pw_poisson_implicit) has_implicit_ps = .TRUE.
3886 :
3887 : ! Write the dielectric constant into a cube file
3888 : do_dielectric_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, &
3889 12969 : "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE"), cp_p_file)
3890 12969 : IF (has_implicit_ps .AND. do_dielectric_cube) THEN
3891 0 : append_cube = section_get_lval(input, "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%APPEND")
3892 0 : my_pos_cube = "REWIND"
3893 0 : IF (append_cube) THEN
3894 0 : my_pos_cube = "APPEND"
3895 : END IF
3896 0 : mpi_io = .TRUE.
3897 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE", &
3898 : extension=".cube", middle_name="DIELECTRIC_CONSTANT", file_position=my_pos_cube, &
3899 0 : mpi_io=mpi_io)
3900 0 : CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
3901 0 : CALL auxbas_pw_pool%create_pw(aux_r)
3902 :
3903 0 : boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
3904 0 : SELECT CASE (boundary_condition)
3905 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
3906 0 : CALL pw_copy(poisson_env%implicit_env%dielectric%eps, aux_r)
3907 : CASE (MIXED_BC, NEUMANN_BC)
3908 : CALL pw_shrink(pw_env%poisson_env%parameters%ps_implicit_params%neumann_directions, &
3909 : pw_env%poisson_env%implicit_env%dct_env%dests_shrink, &
3910 : pw_env%poisson_env%implicit_env%dct_env%srcs_shrink, &
3911 : pw_env%poisson_env%implicit_env%dct_env%bounds_local_shftd, &
3912 0 : poisson_env%implicit_env%dielectric%eps, aux_r)
3913 : END SELECT
3914 :
3915 : CALL cp_pw_to_cube(aux_r, unit_nr, "DIELECTRIC CONSTANT", particles=particles, &
3916 : stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%STRIDE"), &
3917 : max_file_size_mb=section_get_rval(dft_section, "PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE%MAX_FILE_SIZE_MB"), &
3918 0 : mpi_io=mpi_io)
3919 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
3920 0 : "DFT%PRINT%IMPLICIT_PSOLVER%DIELECTRIC_CUBE", mpi_io=mpi_io)
3921 :
3922 0 : CALL auxbas_pw_pool%give_back_pw(aux_r)
3923 : END IF
3924 :
3925 : ! Write Dirichlet constraint charges into a cube file
3926 : do_cstr_charge_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, &
3927 12969 : "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE"), cp_p_file)
3928 :
3929 12969 : has_dirichlet_bc = .FALSE.
3930 12969 : IF (has_implicit_ps) THEN
3931 86 : boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
3932 86 : IF (boundary_condition == MIXED_PERIODIC_BC .OR. boundary_condition == MIXED_BC) THEN
3933 60 : has_dirichlet_bc = .TRUE.
3934 : END IF
3935 : END IF
3936 :
3937 86 : IF (has_implicit_ps .AND. do_cstr_charge_cube .AND. has_dirichlet_bc) THEN
3938 : append_cube = section_get_lval(input, &
3939 0 : "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%APPEND")
3940 0 : my_pos_cube = "REWIND"
3941 0 : IF (append_cube) THEN
3942 0 : my_pos_cube = "APPEND"
3943 : END IF
3944 0 : mpi_io = .TRUE.
3945 : unit_nr = cp_print_key_unit_nr(logger, input, &
3946 : "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE", &
3947 : extension=".cube", middle_name="dirichlet_cstr_charge", file_position=my_pos_cube, &
3948 0 : mpi_io=mpi_io)
3949 0 : CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
3950 0 : CALL auxbas_pw_pool%create_pw(aux_r)
3951 :
3952 0 : boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
3953 0 : SELECT CASE (boundary_condition)
3954 : CASE (MIXED_PERIODIC_BC)
3955 0 : CALL pw_copy(poisson_env%implicit_env%cstr_charge, aux_r)
3956 : CASE (MIXED_BC)
3957 : CALL pw_shrink(pw_env%poisson_env%parameters%ps_implicit_params%neumann_directions, &
3958 : pw_env%poisson_env%implicit_env%dct_env%dests_shrink, &
3959 : pw_env%poisson_env%implicit_env%dct_env%srcs_shrink, &
3960 : pw_env%poisson_env%implicit_env%dct_env%bounds_local_shftd, &
3961 0 : poisson_env%implicit_env%cstr_charge, aux_r)
3962 : END SELECT
3963 :
3964 : CALL cp_pw_to_cube(aux_r, unit_nr, "DIRICHLET CONSTRAINT CHARGE", particles=particles, &
3965 : stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%STRIDE"), &
3966 : max_file_size_mb=section_get_rval(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE%MAX_FILE_SIZE_MB"), &
3967 0 : mpi_io=mpi_io)
3968 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
3969 0 : "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_CSTR_CHARGE_CUBE", mpi_io=mpi_io)
3970 :
3971 0 : CALL auxbas_pw_pool%give_back_pw(aux_r)
3972 : END IF
3973 :
3974 : ! Write Dirichlet type constranits into cube files
3975 : do_dirichlet_bc_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, &
3976 12969 : "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE"), cp_p_file)
3977 12969 : has_dirichlet_bc = .FALSE.
3978 12969 : IF (has_implicit_ps) THEN
3979 86 : boundary_condition = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
3980 86 : IF (boundary_condition == MIXED_PERIODIC_BC .OR. boundary_condition == MIXED_BC) THEN
3981 : has_dirichlet_bc = .TRUE.
3982 : END IF
3983 : END IF
3984 :
3985 60 : IF (has_implicit_ps .AND. has_dirichlet_bc .AND. do_dirichlet_bc_cube) THEN
3986 0 : append_cube = section_get_lval(input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%APPEND")
3987 0 : my_pos_cube = "REWIND"
3988 0 : IF (append_cube) THEN
3989 0 : my_pos_cube = "APPEND"
3990 : END IF
3991 0 : tile_cubes = section_get_lval(input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%TILE_CUBES")
3992 :
3993 0 : CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
3994 0 : CALL auxbas_pw_pool%create_pw(aux_r)
3995 0 : CALL pw_zero(aux_r)
3996 :
3997 0 : IF (tile_cubes) THEN
3998 : ! one cube file per tile
3999 0 : n_cstr = SIZE(poisson_env%implicit_env%contacts)
4000 0 : DO j = 1, n_cstr
4001 0 : n_tiles = poisson_env%implicit_env%contacts(j)%dirichlet_bc%n_tiles
4002 0 : DO i = 1, n_tiles
4003 : filename = "dirichlet_cstr_"//TRIM(ADJUSTL(cp_to_string(j)))// &
4004 0 : "_tile_"//TRIM(ADJUSTL(cp_to_string(i)))
4005 0 : mpi_io = .TRUE.
4006 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", &
4007 : extension=".cube", middle_name=filename, file_position=my_pos_cube, &
4008 0 : mpi_io=mpi_io)
4009 :
4010 0 : CALL pw_copy(poisson_env%implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, aux_r)
4011 :
4012 : CALL cp_pw_to_cube(aux_r, unit_nr, "DIRICHLET TYPE CONSTRAINT", particles=particles, &
4013 : stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%STRIDE"), &
4014 : max_file_size_mb=section_get_rval(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%MAX_FILE_SIZE_MB"), &
4015 0 : mpi_io=mpi_io)
4016 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
4017 0 : "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", mpi_io=mpi_io)
4018 : END DO
4019 : END DO
4020 : ELSE
4021 : ! a single cube file
4022 0 : NULLIFY (dirichlet_tile)
4023 0 : ALLOCATE (dirichlet_tile)
4024 0 : CALL auxbas_pw_pool%create_pw(dirichlet_tile)
4025 0 : CALL pw_zero(dirichlet_tile)
4026 0 : mpi_io = .TRUE.
4027 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", &
4028 : extension=".cube", middle_name="DIRICHLET_CSTR", file_position=my_pos_cube, &
4029 0 : mpi_io=mpi_io)
4030 :
4031 0 : n_cstr = SIZE(poisson_env%implicit_env%contacts)
4032 0 : DO j = 1, n_cstr
4033 0 : n_tiles = poisson_env%implicit_env%contacts(j)%dirichlet_bc%n_tiles
4034 0 : DO i = 1, n_tiles
4035 0 : CALL pw_copy(poisson_env%implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw, dirichlet_tile)
4036 0 : CALL pw_axpy(dirichlet_tile, aux_r)
4037 : END DO
4038 : END DO
4039 :
4040 : CALL cp_pw_to_cube(aux_r, unit_nr, "DIRICHLET TYPE CONSTRAINT", particles=particles, &
4041 : stride=section_get_ivals(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%STRIDE"), &
4042 : max_file_size_mb=section_get_rval(dft_section, "PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE%MAX_FILE_SIZE_MB"), &
4043 0 : mpi_io=mpi_io)
4044 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
4045 0 : "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE", mpi_io=mpi_io)
4046 0 : CALL auxbas_pw_pool%give_back_pw(dirichlet_tile)
4047 0 : DEALLOCATE (dirichlet_tile)
4048 : END IF
4049 :
4050 0 : CALL auxbas_pw_pool%give_back_pw(aux_r)
4051 : END IF
4052 :
4053 12969 : CALL timestop(handle)
4054 :
4055 12969 : END SUBROUTINE qs_scf_post_ps_implicit
4056 :
4057 : !**************************************************************************************************
4058 : !> \brief write an adjacency (interaction) matrix
4059 : !> \param qs_env qs environment
4060 : !> \param input the input
4061 : !> \author Mohammad Hossein Bani-Hashemian
4062 : ! **************************************************************************************************
4063 12969 : SUBROUTINE write_adjacency_matrix(qs_env, input)
4064 : TYPE(qs_environment_type), POINTER :: qs_env
4065 : TYPE(section_vals_type), POINTER :: input
4066 :
4067 : CHARACTER(len=*), PARAMETER :: routineN = 'write_adjacency_matrix'
4068 :
4069 : INTEGER :: adjm_size, colind, handle, iatom, ikind, &
4070 : ind, jatom, jkind, k, natom, nkind, &
4071 : output_unit, rowind, unit_nr
4072 12969 : INTEGER, ALLOCATABLE, DIMENSION(:) :: interact_adjm
4073 : LOGICAL :: do_adjm_write, do_symmetric
4074 : TYPE(cp_logger_type), POINTER :: logger
4075 12969 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
4076 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
4077 : TYPE(mp_para_env_type), POINTER :: para_env
4078 : TYPE(neighbor_list_iterator_p_type), &
4079 12969 : DIMENSION(:), POINTER :: nl_iterator
4080 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
4081 12969 : POINTER :: nl
4082 12969 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
4083 : TYPE(section_vals_type), POINTER :: dft_section
4084 :
4085 12969 : CALL timeset(routineN, handle)
4086 :
4087 12969 : NULLIFY (dft_section)
4088 :
4089 12969 : logger => cp_get_default_logger()
4090 12969 : output_unit = cp_logger_get_default_io_unit(logger)
4091 :
4092 12969 : dft_section => section_vals_get_subs_vals(input, "DFT")
4093 : do_adjm_write = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
4094 12969 : "PRINT%ADJMAT_WRITE"), cp_p_file)
4095 :
4096 12969 : IF (do_adjm_write) THEN
4097 28 : NULLIFY (qs_kind_set, nl_iterator)
4098 28 : NULLIFY (basis_set_list_a, basis_set_list_b, basis_set_a, basis_set_b)
4099 :
4100 28 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, sab_orb=nl, natom=natom, para_env=para_env)
4101 :
4102 28 : nkind = SIZE(qs_kind_set)
4103 28 : CPASSERT(SIZE(nl) > 0)
4104 28 : CALL get_neighbor_list_set_p(neighbor_list_sets=nl, symmetric=do_symmetric)
4105 28 : CPASSERT(do_symmetric)
4106 216 : ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
4107 28 : CALL basis_set_list_setup(basis_set_list_a, "ORB", qs_kind_set)
4108 28 : CALL basis_set_list_setup(basis_set_list_b, "ORB", qs_kind_set)
4109 :
4110 28 : adjm_size = ((natom + 1)*natom)/2
4111 84 : ALLOCATE (interact_adjm(4*adjm_size))
4112 28 : interact_adjm = 0
4113 :
4114 28 : NULLIFY (nl_iterator)
4115 28 : CALL neighbor_list_iterator_create(nl_iterator, nl)
4116 2021 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
4117 : CALL get_iterator_info(nl_iterator, &
4118 : ikind=ikind, jkind=jkind, &
4119 1993 : iatom=iatom, jatom=jatom)
4120 :
4121 1993 : basis_set_a => basis_set_list_a(ikind)%gto_basis_set
4122 1993 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
4123 1993 : basis_set_b => basis_set_list_b(jkind)%gto_basis_set
4124 1993 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
4125 :
4126 : ! move everything to the upper triangular part
4127 1993 : IF (iatom <= jatom) THEN
4128 : rowind = iatom
4129 : colind = jatom
4130 : ELSE
4131 670 : rowind = jatom
4132 670 : colind = iatom
4133 : ! swap the kinds too
4134 : ikind = ikind + jkind
4135 670 : jkind = ikind - jkind
4136 670 : ikind = ikind - jkind
4137 : END IF
4138 :
4139 : ! indexing upper triangular matrix
4140 1993 : ind = adjm_size - (natom - rowind + 1)*((natom - rowind + 1) + 1)/2 + colind - rowind + 1
4141 : ! convert the upper triangular matrix into a adjm_size x 4 matrix
4142 : ! columns are: iatom, jatom, ikind, jkind
4143 1993 : interact_adjm((ind - 1)*4 + 1) = rowind
4144 1993 : interact_adjm((ind - 1)*4 + 2) = colind
4145 1993 : interact_adjm((ind - 1)*4 + 3) = ikind
4146 1993 : interact_adjm((ind - 1)*4 + 4) = jkind
4147 : END DO
4148 :
4149 28 : CALL para_env%sum(interact_adjm)
4150 :
4151 : unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%ADJMAT_WRITE", &
4152 : extension=".adjmat", file_form="FORMATTED", &
4153 28 : file_status="REPLACE")
4154 28 : IF (unit_nr > 0) THEN
4155 14 : WRITE (unit_nr, "(1A,2X,1A,5X,1A,4X,A5,3X,A5)") "#", "iatom", "jatom", "ikind", "jkind"
4156 88 : DO k = 1, 4*adjm_size, 4
4157 : ! print only the interacting atoms
4158 88 : IF (interact_adjm(k) > 0 .AND. interact_adjm(k + 1) > 0) THEN
4159 74 : WRITE (unit_nr, "(I8,2X,I8,3X,I6,2X,I6)") interact_adjm(k:k + 3)
4160 : END IF
4161 : END DO
4162 : END IF
4163 :
4164 28 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%ADJMAT_WRITE")
4165 :
4166 28 : CALL neighbor_list_iterator_release(nl_iterator)
4167 56 : DEALLOCATE (basis_set_list_a, basis_set_list_b)
4168 : END IF
4169 :
4170 12969 : CALL timestop(handle)
4171 :
4172 25938 : END SUBROUTINE write_adjacency_matrix
4173 :
4174 : ! **************************************************************************************************
4175 : !> \brief Updates Hartree potential with MP2 density. Important for REPEAT charges
4176 : !> \param rho ...
4177 : !> \param qs_env ...
4178 : !> \author Vladimir Rybkin
4179 : ! **************************************************************************************************
4180 322 : SUBROUTINE update_hartree_with_mp2(rho, qs_env)
4181 : TYPE(qs_rho_type), POINTER :: rho
4182 : TYPE(qs_environment_type), POINTER :: qs_env
4183 :
4184 : LOGICAL :: use_virial
4185 : TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
4186 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
4187 : TYPE(pw_env_type), POINTER :: pw_env
4188 : TYPE(pw_poisson_type), POINTER :: poisson_env
4189 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
4190 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
4191 : TYPE(qs_energy_type), POINTER :: energy
4192 : TYPE(virial_type), POINTER :: virial
4193 :
4194 322 : NULLIFY (auxbas_pw_pool, pw_env, poisson_env, energy, rho_core, v_hartree_rspace, virial)
4195 : CALL get_qs_env(qs_env, pw_env=pw_env, energy=energy, &
4196 : rho_core=rho_core, virial=virial, &
4197 322 : v_hartree_rspace=v_hartree_rspace)
4198 :
4199 322 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
4200 :
4201 : IF (.NOT. use_virial) THEN
4202 :
4203 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
4204 268 : poisson_env=poisson_env)
4205 268 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
4206 268 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
4207 :
4208 268 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
4209 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%hartree, &
4210 268 : v_hartree_gspace, rho_core=rho_core)
4211 :
4212 268 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
4213 268 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
4214 :
4215 268 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
4216 268 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
4217 : END IF
4218 :
4219 322 : END SUBROUTINE update_hartree_with_mp2
4220 :
4221 0 : END MODULE qs_scf_post_gpw
|
