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