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