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