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