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