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