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