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