Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Does all kind of post scf calculations for DFTB
10 : !> \par History
11 : !> Started as a copy from the GPW file
12 : !> - Revise MO information printout (10.05.2021, MK)
13 : !> \author JHU (03.2013)
14 : ! **************************************************************************************************
15 : MODULE qs_scf_post_tb
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE cell_types, ONLY: cell_type,&
19 : pbc
20 : USE cp_array_utils, ONLY: cp_1d_r_p_type
21 : USE cp_blacs_env, ONLY: cp_blacs_env_type
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
24 : dbcsr_type
25 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
26 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
27 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
28 : cp_fm_cholesky_reduce,&
29 : cp_fm_cholesky_restore
30 : USE cp_fm_diag, ONLY: choose_eigv_solver
31 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
32 : cp_fm_struct_release,&
33 : cp_fm_struct_type
34 : USE cp_fm_types, ONLY: cp_fm_create,&
35 : cp_fm_get_info,&
36 : cp_fm_init_random,&
37 : cp_fm_release,&
38 : cp_fm_to_fm_submat,&
39 : cp_fm_type
40 : USE cp_log_handling, ONLY: cp_get_default_logger,&
41 : cp_logger_get_default_io_unit,&
42 : cp_logger_type
43 : USE cp_output_handling, ONLY: cp_p_file,&
44 : cp_print_key_finished_output,&
45 : cp_print_key_should_output,&
46 : cp_print_key_unit_nr
47 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
48 : USE cp_result_methods, ONLY: cp_results_erase,&
49 : put_results
50 : USE cp_result_types, ONLY: cp_result_type
51 : USE eeq_method, ONLY: eeq_print
52 : USE input_constants, ONLY: ot_precond_full_all,&
53 : ot_precond_full_all_covariant
54 : USE input_section_types, ONLY: section_get_ival,&
55 : section_get_ivals,&
56 : section_get_lval,&
57 : section_get_rval,&
58 : section_vals_get,&
59 : section_vals_get_subs_vals,&
60 : section_vals_type,&
61 : section_vals_val_get
62 : USE kinds, ONLY: default_path_length,&
63 : default_string_length,&
64 : dp
65 : USE machine, ONLY: m_flush
66 : USE mathconstants, ONLY: twopi,&
67 : z_one,&
68 : z_zero
69 : USE memory_utilities, ONLY: reallocate
70 : USE message_passing, ONLY: mp_para_env_type
71 : USE molden_utils, ONLY: write_mos_molden
72 : USE moments_utils, ONLY: get_reference_point
73 : USE mulliken, ONLY: mulliken_charges
74 : USE particle_list_types, ONLY: particle_list_type
75 : USE particle_types, ONLY: particle_type
76 : USE physcon, ONLY: debye
77 : USE population_analyses, ONLY: lowdin_population_analysis
78 : USE preconditioner_types, ONLY: preconditioner_type
79 : USE pw_env_methods, ONLY: pw_env_create,&
80 : pw_env_rebuild
81 : USE pw_env_types, ONLY: pw_env_get,&
82 : pw_env_release,&
83 : pw_env_type
84 : USE pw_grid_types, ONLY: pw_grid_type
85 : USE pw_methods, ONLY: pw_axpy,&
86 : pw_copy,&
87 : pw_derive,&
88 : pw_scale,&
89 : pw_transfer,&
90 : pw_zero
91 : USE pw_poisson_types, ONLY: do_ewald_none,&
92 : greens_fn_type,&
93 : pw_green_create,&
94 : pw_green_release,&
95 : pw_poisson_analytic,&
96 : pw_poisson_parameter_type
97 : USE pw_pool_types, ONLY: pw_pool_p_type,&
98 : pw_pool_type
99 : USE pw_types, ONLY: pw_c1d_gs_type,&
100 : pw_r3d_rs_type
101 : USE qs_collocate_density, ONLY: calculate_rho_core,&
102 : calculate_rho_elec,&
103 : calculate_wavefunction
104 : USE qs_dftb_types, ONLY: qs_dftb_atom_type
105 : USE qs_dftb_utils, ONLY: get_dftb_atom_param
106 : USE qs_dos, ONLY: calculate_dos,&
107 : calculate_dos_kp
108 : USE qs_dos_utils, ONLY: get_dos_pdos_flags
109 : USE qs_elf_methods, ONLY: qs_elf_calc
110 : USE qs_energy_window, ONLY: energy_windows
111 : USE qs_environment_types, ONLY: get_qs_env,&
112 : qs_environment_type
113 : USE qs_kind_types, ONLY: get_qs_kind,&
114 : qs_kind_type
115 : USE qs_ks_types, ONLY: get_ks_env,&
116 : qs_ks_env_type,&
117 : set_ks_env
118 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues,&
119 : make_mo_eig
120 : USE qs_mo_occupation, ONLY: set_mo_occupation
121 : USE qs_mo_types, ONLY: get_mo_set,&
122 : mo_set_type
123 : USE qs_ot_eigensolver, ONLY: ot_eigensolver
124 : USE qs_pdos, ONLY: calculate_projected_dos,&
125 : calculate_projected_dos_kp
126 : USE qs_rho_methods, ONLY: qs_rho_rebuild
127 : USE qs_rho_types, ONLY: qs_rho_get,&
128 : qs_rho_set,&
129 : qs_rho_type
130 : USE qs_scf_csr_write, ONLY: write_hcore_matrix_csr,&
131 : write_ks_matrix_csr,&
132 : write_p_matrix_csr,&
133 : write_s_matrix_csr
134 : USE qs_scf_output, ONLY: qs_scf_write_mos
135 : USE qs_scf_types, ONLY: ot_method_nr,&
136 : qs_scf_env_type
137 : USE qs_scf_wfn_mix, ONLY: wfn_mix
138 : USE qs_subsys_types, ONLY: qs_subsys_get,&
139 : qs_subsys_type
140 : USE scf_control_types, ONLY: scf_control_type
141 : USE stm_images, ONLY: th_stm_image
142 : USE task_list_methods, ONLY: generate_qs_task_list
143 : USE task_list_types, ONLY: allocate_task_list,&
144 : task_list_type
145 : USE xtb_qresp, ONLY: build_xtb_qresp
146 : USE xtb_types, ONLY: get_xtb_atom_param,&
147 : xtb_atom_type
148 : #include "./base/base_uses.f90"
149 :
150 : IMPLICIT NONE
151 : PRIVATE
152 :
153 : ! Global parameters
154 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_post_tb'
155 : PUBLIC :: scf_post_calculation_tb, make_lumo_tb, rebuild_pw_env
156 :
157 : ! **************************************************************************************************
158 :
159 : CONTAINS
160 :
161 : ! **************************************************************************************************
162 : !> \brief collects possible post - scf calculations and prints info / computes properties.
163 : !> \param qs_env ...
164 : !> \param tb_type ...
165 : !> \param no_mos ...
166 : !> \par History
167 : !> 03.2013 copy of scf_post_gpw
168 : !> \author JHU
169 : !> \note
170 : ! **************************************************************************************************
171 11678 : SUBROUTINE scf_post_calculation_tb(qs_env, tb_type, no_mos)
172 :
173 : TYPE(qs_environment_type), POINTER :: qs_env
174 : CHARACTER(LEN=*) :: tb_type
175 : LOGICAL, INTENT(IN) :: no_mos
176 :
177 : CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_tb'
178 :
179 : CHARACTER(LEN=6) :: ana
180 : CHARACTER(LEN=default_string_length) :: aname
181 : INTEGER :: after, gfn_type, handle, homo, iat, iatom, ikind, img, ispin, iw, nat, natom, &
182 : nkind, nlumo_stm, nlumos, nspins, print_level, unit_nr
183 : LOGICAL :: do_cube, do_curve, do_dos, do_kpoints, do_pdos, do_projected_dos, explicit, gfn0, &
184 : has_homo, omit_headers, print_it, rebuild, vdip
185 : REAL(KIND=dp) :: zeff
186 11678 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge, zcharge
187 : REAL(KIND=dp), DIMENSION(2, 2) :: homo_lumo
188 11678 : REAL(KIND=dp), DIMENSION(:), POINTER :: echarge, mo_eigenvalues
189 11678 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges
190 11678 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
191 : TYPE(cell_type), POINTER :: cell
192 11678 : TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: unoccupied_evals_stm
193 11678 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: unoccupied_orbs_stm
194 : TYPE(cp_fm_type), POINTER :: mo_coeff
195 : TYPE(cp_logger_type), POINTER :: logger
196 11678 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, mo_derivs
197 11678 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_p, matrix_s
198 : TYPE(dbcsr_type), POINTER :: mo_coeff_deriv
199 : TYPE(dft_control_type), POINTER :: dft_control
200 11678 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
201 : TYPE(mp_para_env_type), POINTER :: para_env
202 : TYPE(particle_list_type), POINTER :: particles
203 11678 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
204 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
205 11678 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
206 : TYPE(qs_rho_type), POINTER :: rho
207 : TYPE(qs_scf_env_type), POINTER :: scf_env
208 : TYPE(qs_subsys_type), POINTER :: subsys
209 : TYPE(scf_control_type), POINTER :: scf_control
210 : TYPE(section_vals_type), POINTER :: dft_section, moments_section, print_key, &
211 : print_section, sprint_section, &
212 : wfn_mix_section
213 : TYPE(xtb_atom_type), POINTER :: xtb_kind
214 :
215 11678 : CALL timeset(routineN, handle)
216 :
217 11678 : logger => cp_get_default_logger()
218 :
219 11678 : gfn0 = .FALSE.
220 11678 : vdip = .FALSE.
221 11678 : CALL get_qs_env(qs_env, dft_control=dft_control)
222 19936 : SELECT CASE (TRIM(tb_type))
223 : CASE ("DFTB")
224 : CASE ("xTB")
225 8258 : gfn_type = dft_control%qs_control%xtb_control%gfn_type
226 8258 : gfn0 = (gfn_type == 0)
227 8258 : vdip = dft_control%qs_control%xtb_control%var_dipole
228 : CASE DEFAULT
229 11678 : CPABORT("unknown TB type")
230 : END SELECT
231 :
232 11678 : CPASSERT(ASSOCIATED(qs_env))
233 11678 : NULLIFY (rho, para_env, matrix_s, matrix_p)
234 : CALL get_qs_env(qs_env, scf_env=scf_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
235 : rho=rho, natom=natom, para_env=para_env, &
236 11678 : particle_set=particle_set, do_kpoints=do_kpoints, matrix_s_kp=matrix_s)
237 11678 : nspins = dft_control%nspins
238 11678 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
239 : ! Mulliken charges
240 70068 : ALLOCATE (charges(natom, nspins), mcharge(natom))
241 : !
242 11678 : CALL mulliken_charges(matrix_p, matrix_s, para_env, charges)
243 : !
244 35034 : ALLOCATE (zcharge(natom))
245 11678 : nkind = SIZE(atomic_kind_set)
246 36848 : DO ikind = 1, nkind
247 25170 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
248 32332 : SELECT CASE (TRIM(tb_type))
249 : CASE ("DFTB")
250 7162 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
251 25170 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
252 : CASE ("xTB")
253 18008 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
254 18008 : CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
255 : CASE DEFAULT
256 50340 : CPABORT("unknown TB type")
257 : END SELECT
258 144090 : DO iatom = 1, nat
259 82072 : iat = atomic_kind_set(ikind)%atom_list(iatom)
260 168774 : mcharge(iat) = zeff - SUM(charges(iat, 1:nspins))
261 107242 : zcharge(iat) = zeff
262 : END DO
263 : END DO
264 :
265 11678 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
266 11678 : print_section => section_vals_get_subs_vals(dft_section, "PRINT")
267 :
268 : ! Mulliken
269 11678 : print_key => section_vals_get_subs_vals(print_section, "MULLIKEN")
270 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
271 : unit_nr = cp_print_key_unit_nr(logger, print_section, "MULLIKEN", &
272 3844 : extension=".mulliken", log_filename=.FALSE.)
273 3844 : IF (unit_nr > 0) THEN
274 1933 : WRITE (UNIT=unit_nr, FMT="(/,/,T2,A)") "MULLIKEN POPULATION ANALYSIS"
275 1933 : IF (nspins == 1) THEN
276 : WRITE (UNIT=unit_nr, FMT="(/,T2,A,T70,A)") &
277 1687 : " # Atom Element Kind Atomic population", " Net charge"
278 4550 : DO ikind = 1, nkind
279 2863 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
280 2863 : aname = ""
281 567 : SELECT CASE (tb_type)
282 : CASE ("DFTB")
283 567 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
284 567 : CALL get_dftb_atom_param(dftb_kind, name=aname)
285 : CASE ("xTB")
286 2296 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
287 2296 : CALL get_xtb_atom_param(xtb_kind, symbol=aname(1:2))
288 : CASE DEFAULT
289 2863 : CPABORT("unknown TB type")
290 : END SELECT
291 2863 : ana = ADJUSTR(TRIM(ADJUSTL(aname)))
292 15408 : DO iatom = 1, nat
293 7995 : iat = atomic_kind_set(ikind)%atom_list(iatom)
294 : WRITE (UNIT=unit_nr, &
295 : FMT="(T2,I7,5X,A6,I6,T39,F12.6,T69,F12.6)") &
296 10858 : iat, ADJUSTL(ana), ikind, charges(iat, 1), mcharge(iat)
297 : END DO
298 : END DO
299 : WRITE (UNIT=unit_nr, &
300 : FMT="(T2,A,T39,F12.6,T69,F12.6,/)") &
301 17677 : "# Total charge", SUM(charges(:, 1)), SUM(mcharge(:))
302 : ELSE
303 : WRITE (UNIT=unit_nr, FMT="(/,T2,A)") &
304 246 : "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
305 677 : DO ikind = 1, nkind
306 431 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
307 431 : aname = ""
308 3 : SELECT CASE (tb_type)
309 : CASE ("DFTB")
310 3 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
311 3 : CALL get_dftb_atom_param(dftb_kind, name=aname)
312 : CASE ("xTB")
313 428 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
314 428 : CALL get_xtb_atom_param(xtb_kind, symbol=aname(1:2))
315 : CASE DEFAULT
316 431 : CPABORT("unknown TB type")
317 : END SELECT
318 431 : ana = ADJUSTR(TRIM(ADJUSTL(aname)))
319 1791 : DO iatom = 1, nat
320 683 : iat = atomic_kind_set(ikind)%atom_list(iatom)
321 : WRITE (UNIT=unit_nr, &
322 : FMT="(T2,I6,3X,A6,I6,T29,4(1X,F12.6))") &
323 2049 : iat, ADJUSTL(ana), ikind, charges(iat, 1:2), mcharge(iat), &
324 1797 : charges(iat, 1) - charges(iat, 2)
325 : END DO
326 : END DO
327 : WRITE (UNIT=unit_nr, &
328 : FMT="(T2,A,T29,4(1X,F12.6),/)") &
329 2295 : "# Total charge and spin", SUM(charges(:, 1)), SUM(charges(:, 2)), SUM(mcharge(:))
330 : END IF
331 1933 : CALL m_flush(unit_nr)
332 : END IF
333 3844 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
334 : END IF
335 :
336 : ! Compute the Lowdin charges
337 11678 : print_key => section_vals_get_subs_vals(print_section, "LOWDIN")
338 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
339 48 : SELECT CASE (tb_type)
340 : CASE ("DFTB")
341 48 : CPWARN("Lowdin population analysis not implemented for DFTB method.")
342 : CASE ("xTB")
343 : unit_nr = cp_print_key_unit_nr(logger, print_section, "LOWDIN", extension=".lowdin", &
344 26 : log_filename=.FALSE.)
345 26 : print_level = 1
346 26 : CALL section_vals_val_get(print_key, "PRINT_GOP", l_val=print_it)
347 26 : IF (print_it) print_level = 2
348 26 : CALL section_vals_val_get(print_key, "PRINT_ALL", l_val=print_it)
349 26 : IF (print_it) print_level = 3
350 26 : IF (do_kpoints) THEN
351 2 : CPWARN("Lowdin charges not implemented for k-point calculations!")
352 : ELSE
353 24 : CALL lowdin_population_analysis(qs_env, unit_nr, print_level)
354 : END IF
355 26 : CALL cp_print_key_finished_output(unit_nr, logger, print_section, "LOWDIN")
356 : CASE DEFAULT
357 126 : CPABORT("unknown TB type")
358 : END SELECT
359 : END IF
360 :
361 : ! EEQ Charges
362 11678 : print_key => section_vals_get_subs_vals(print_section, "EEQ_CHARGES")
363 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
364 : unit_nr = cp_print_key_unit_nr(logger, print_section, "EEQ_CHARGES", &
365 4 : extension=".eeq", log_filename=.FALSE.)
366 4 : CALL eeq_print(qs_env, unit_nr, print_level, ext=gfn0)
367 4 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
368 : END IF
369 :
370 : ! Hirshfeld
371 11678 : print_key => section_vals_get_subs_vals(print_section, "HIRSHFELD")
372 11678 : CALL section_vals_get(print_key, explicit=explicit)
373 11678 : IF (explicit) THEN
374 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
375 0 : CPWARN("Hirshfeld charges not available for TB methods.")
376 : END IF
377 : END IF
378 :
379 : ! MAO
380 11678 : print_key => section_vals_get_subs_vals(print_section, "MAO_ANALYSIS")
381 11678 : CALL section_vals_get(print_key, explicit=explicit)
382 11678 : IF (explicit) THEN
383 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
384 0 : CPWARN("MAO analysis not available for TB methods.")
385 : END IF
386 : END IF
387 :
388 : ! ED
389 11678 : print_key => section_vals_get_subs_vals(print_section, "ENERGY_DECOMPOSITION_ANALYSIS")
390 11678 : CALL section_vals_get(print_key, explicit=explicit)
391 11678 : IF (explicit) THEN
392 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
393 0 : CPWARN("ED analysis not available for TB methods.")
394 : END IF
395 : END IF
396 :
397 : ! Dipole Moments
398 11678 : print_key => section_vals_get_subs_vals(print_section, "MOMENTS")
399 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
400 : unit_nr = cp_print_key_unit_nr(logger, print_section, "MOMENTS", &
401 1038 : extension=".data", middle_name="tb_dipole", log_filename=.FALSE.)
402 1038 : moments_section => section_vals_get_subs_vals(print_section, "MOMENTS")
403 1038 : IF (gfn0) THEN
404 158 : NULLIFY (echarge)
405 158 : CALL get_qs_env(qs_env, eeq=echarge)
406 158 : CPASSERT(ASSOCIATED(echarge))
407 158 : IF (vdip) THEN
408 58 : CALL build_xtb_qresp(qs_env, mcharge)
409 290 : mcharge(1:natom) = echarge(1:natom) - mcharge(1:natom)
410 : END IF
411 158 : CALL tb_dipole(qs_env, moments_section, unit_nr, mcharge)
412 : ELSE
413 880 : CALL tb_dipole(qs_env, moments_section, unit_nr, mcharge)
414 : END IF
415 1038 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
416 : END IF
417 :
418 11678 : DEALLOCATE (charges, mcharge)
419 :
420 : ! MO
421 11678 : IF (.NOT. no_mos) THEN
422 11530 : print_key => section_vals_get_subs_vals(print_section, "MO")
423 11530 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
424 154 : CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.TRUE.)
425 154 : IF (.NOT. do_kpoints) THEN
426 102 : SELECT CASE (tb_type)
427 : CASE ("DFTB")
428 : CASE ("xTB")
429 102 : sprint_section => section_vals_get_subs_vals(dft_section, "PRINT%MO_MOLDEN")
430 102 : CALL get_qs_env(qs_env, mos=mos, cell=cell)
431 : CALL write_mos_molden(mos, qs_kind_set, particle_set, sprint_section, cell=cell, &
432 102 : qs_env=qs_env, calc_energies=.TRUE.)
433 : CASE DEFAULT
434 144 : CPABORT("Unknown TB type")
435 : END SELECT
436 : END IF
437 : END IF
438 : END IF
439 :
440 : ! Wavefunction mixing
441 11678 : IF (.NOT. no_mos) THEN
442 11530 : wfn_mix_section => section_vals_get_subs_vals(dft_section, "PRINT%WFN_MIX")
443 11530 : CALL section_vals_get(wfn_mix_section, explicit=explicit)
444 11530 : IF (explicit .AND. .NOT. qs_env%run_rtp) CALL wfn_mix_tb(qs_env, dft_section, scf_env)
445 : END IF
446 :
447 11678 : IF (.NOT. no_mos) THEN
448 11530 : print_key => section_vals_get_subs_vals(print_section, "DOS")
449 11530 : do_dos = BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
450 11530 : CALL get_dos_pdos_flags(print_key, do_dos, do_projected_dos, do_pdos, do_curve)
451 11530 : IF (do_dos) THEN
452 22 : IF (do_kpoints) THEN
453 16 : CALL calculate_dos_kp(qs_env, dft_section)
454 16 : IF (do_curve) CALL calculate_dos_kp(qs_env, dft_section, write_curve_output=.TRUE.)
455 : ELSE
456 6 : CALL get_qs_env(qs_env, mos=mos)
457 6 : CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear)
458 6 : IF (do_curve) CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear, &
459 0 : write_curve_output=.TRUE.)
460 : END IF
461 : END IF
462 :
463 : ! Projected density-of-states outputs
464 11530 : IF (do_projected_dos) THEN
465 18 : IF (do_kpoints) THEN
466 : CALL calculate_projected_dos_kp(qs_env, dft_section, pdos_print_key="PRINT%DOS", &
467 14 : write_pdos=do_pdos, write_pdos_curve=do_curve)
468 : ELSE
469 4 : CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv)
470 8 : DO ispin = 1, dft_control%nspins
471 4 : IF (scf_env%method == ot_method_nr) THEN
472 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
473 0 : eigenvalues=mo_eigenvalues)
474 0 : IF (ASSOCIATED(qs_env%mo_derivs)) THEN
475 0 : mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix
476 : ELSE
477 0 : mo_coeff_deriv => NULL()
478 : END IF
479 : CALL calculate_subspace_eigenvalues(mo_coeff, ks_rmpv(ispin)%matrix, mo_eigenvalues, &
480 : do_rotation=.TRUE., &
481 0 : co_rotate_dbcsr=mo_coeff_deriv)
482 0 : CALL set_mo_occupation(mo_set=mos(ispin))
483 : END IF
484 8 : IF (dft_control%nspins == 2) THEN
485 : CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
486 : qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin, &
487 0 : pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
488 : ELSE
489 : CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
490 : qs_kind_set, particle_set, qs_env, dft_section, &
491 4 : pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve)
492 : END IF
493 : END DO
494 : END IF
495 : END IF
496 : END IF
497 :
498 : ! can we do CUBE files?
499 : SELECT CASE (tb_type)
500 : CASE ("DFTB")
501 : do_cube = .FALSE.
502 8258 : rebuild = .FALSE.
503 : CASE ("xTB")
504 8258 : do_cube = .TRUE.
505 8258 : rebuild = .TRUE.
506 : CASE DEFAULT
507 11678 : CPABORT("unknown TB type")
508 : END SELECT
509 :
510 : ! Energy Windows for LS code
511 11678 : print_key => section_vals_get_subs_vals(print_section, "ENERGY_WINDOWS")
512 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
513 74 : IF (do_cube) THEN
514 26 : IF (do_kpoints) THEN
515 2 : CPWARN("Energy Windows not implemented for k-points.")
516 : ELSE
517 : IF (rebuild) THEN
518 24 : CALL rebuild_pw_env(qs_env)
519 : rebuild = .FALSE.
520 : END IF
521 24 : CALL energy_windows(qs_env)
522 : END IF
523 : ELSE
524 48 : CPWARN("Energy Windows not implemented for TB methods.")
525 : END IF
526 : END IF
527 :
528 : ! DENSITY CUBE FILE
529 11678 : print_key => section_vals_get_subs_vals(print_section, "E_DENSITY_CUBE")
530 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
531 72 : IF (do_cube) THEN
532 24 : IF (rebuild) THEN
533 2 : CALL rebuild_pw_env(qs_env)
534 2 : rebuild = .FALSE.
535 : END IF
536 24 : CALL print_e_density(qs_env, zcharge, print_key)
537 : ELSE
538 48 : CPWARN("Electronic density cube file not implemented for TB methods.")
539 : END IF
540 : END IF
541 :
542 : ! TOTAL DENSITY CUBE FILE
543 11678 : print_key => section_vals_get_subs_vals(print_section, "TOT_DENSITY_CUBE")
544 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
545 74 : IF (do_cube) THEN
546 26 : IF (rebuild) THEN
547 2 : CALL rebuild_pw_env(qs_env)
548 2 : rebuild = .FALSE.
549 : END IF
550 26 : CALL print_density_cubes(qs_env, zcharge, print_key, total_density=.TRUE.)
551 : ELSE
552 48 : CPWARN("Total density cube file not implemented for TB methods.")
553 : END IF
554 : END IF
555 :
556 : ! V_Hartree CUBE FILE
557 11678 : print_key => section_vals_get_subs_vals(print_section, "V_HARTREE_CUBE")
558 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
559 72 : IF (do_cube) THEN
560 24 : IF (rebuild) THEN
561 0 : CALL rebuild_pw_env(qs_env)
562 0 : rebuild = .FALSE.
563 : END IF
564 24 : CALL print_density_cubes(qs_env, zcharge, print_key, v_hartree=.TRUE.)
565 : ELSE
566 48 : CPWARN("Hartree potential cube file not implemented for TB methods.")
567 : END IF
568 : END IF
569 :
570 : ! EFIELD CUBE FILE
571 11678 : print_key => section_vals_get_subs_vals(print_section, "EFIELD_CUBE")
572 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
573 72 : IF (do_cube) THEN
574 24 : IF (rebuild) THEN
575 0 : CALL rebuild_pw_env(qs_env)
576 0 : rebuild = .FALSE.
577 : END IF
578 24 : CALL print_density_cubes(qs_env, zcharge, print_key, efield=.TRUE.)
579 : ELSE
580 48 : CPWARN("Efield cube file not implemented for TB methods.")
581 : END IF
582 : END IF
583 :
584 : ! ELF
585 11678 : print_key => section_vals_get_subs_vals(print_section, "ELF_CUBE")
586 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
587 72 : IF (do_cube) THEN
588 24 : IF (rebuild) THEN
589 0 : CALL rebuild_pw_env(qs_env)
590 0 : rebuild = .FALSE.
591 : END IF
592 24 : CALL print_elf(qs_env, zcharge, print_key)
593 : ELSE
594 48 : CPWARN("ELF not implemented for TB methods.")
595 : END IF
596 : END IF
597 :
598 : ! MO CUBES
599 11678 : IF (.NOT. no_mos) THEN
600 11530 : print_key => section_vals_get_subs_vals(print_section, "MO_CUBES")
601 11530 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
602 74 : IF (do_cube) THEN
603 26 : IF (do_kpoints) THEN
604 2 : CPWARN("Print MO cubes not implemented for k-point calculations")
605 : ELSE
606 24 : IF (rebuild) THEN
607 2 : CALL rebuild_pw_env(qs_env)
608 2 : rebuild = .FALSE.
609 : END IF
610 24 : CALL print_mo_cubes(qs_env, zcharge, print_key)
611 : END IF
612 : ELSE
613 48 : CPWARN("Printing of MO cube files not implemented for TB methods.")
614 : END IF
615 : END IF
616 : END IF
617 :
618 : ! STM
619 11678 : IF (.NOT. no_mos) THEN
620 11530 : print_key => section_vals_get_subs_vals(print_section, "STM")
621 11530 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
622 4 : IF (do_cube) THEN
623 4 : IF (rebuild) THEN
624 2 : CALL rebuild_pw_env(qs_env)
625 2 : rebuild = .FALSE.
626 : END IF
627 4 : IF (do_kpoints) THEN
628 0 : CPWARN("STM not implemented for k-point calculations!")
629 : ELSE
630 4 : nlumo_stm = section_get_ival(print_key, "NLUMO")
631 4 : CPASSERT(.NOT. dft_control%restricted)
632 : CALL get_qs_env(qs_env, mos=mos, mo_derivs=mo_derivs, &
633 4 : scf_control=scf_control, matrix_ks=ks_rmpv)
634 4 : CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs)
635 8 : DO ispin = 1, dft_control%nspins
636 4 : CALL get_mo_set(mo_set=mos(ispin), eigenvalues=mo_eigenvalues, homo=homo)
637 8 : homo_lumo(ispin, 1) = mo_eigenvalues(homo)
638 : END DO
639 4 : has_homo = .TRUE.
640 4 : NULLIFY (unoccupied_orbs_stm, unoccupied_evals_stm)
641 4 : IF (nlumo_stm > 0) THEN
642 8 : ALLOCATE (unoccupied_orbs_stm(dft_control%nspins))
643 8 : ALLOCATE (unoccupied_evals_stm(dft_control%nspins))
644 : CALL make_lumo_tb(qs_env, scf_env, unoccupied_orbs_stm, unoccupied_evals_stm, &
645 2 : nlumo_stm, nlumos)
646 : END IF
647 :
648 4 : CALL get_qs_env(qs_env, subsys=subsys)
649 4 : CALL qs_subsys_get(subsys, particles=particles)
650 : CALL th_stm_image(qs_env, print_key, particles, unoccupied_orbs_stm, &
651 4 : unoccupied_evals_stm)
652 :
653 4 : IF (nlumo_stm > 0) THEN
654 4 : DO ispin = 1, dft_control%nspins
655 4 : DEALLOCATE (unoccupied_evals_stm(ispin)%array)
656 : END DO
657 2 : DEALLOCATE (unoccupied_evals_stm)
658 2 : CALL cp_fm_release(unoccupied_orbs_stm)
659 : END IF
660 : END IF
661 : END IF
662 : END IF
663 : END IF
664 :
665 : ! Write the density matrix
666 11678 : CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks)
667 11678 : CALL section_vals_val_get(print_section, "AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
668 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, &
669 : "AO_MATRICES/DENSITY"), cp_p_file)) THEN
670 : iw = cp_print_key_unit_nr(logger, print_section, "AO_MATRICES/DENSITY", &
671 50 : extension=".Log")
672 50 : CALL section_vals_val_get(print_section, "AO_MATRICES%NDIGITS", i_val=after)
673 50 : after = MIN(MAX(after, 1), 16)
674 100 : DO ispin = 1, dft_control%nspins
675 150 : DO img = 1, SIZE(matrix_p, 2)
676 : CALL cp_dbcsr_write_sparse_matrix(matrix_p(ispin, img)%matrix, 4, after, qs_env, &
677 100 : para_env, output_unit=iw, omit_headers=omit_headers)
678 : END DO
679 : END DO
680 50 : CALL cp_print_key_finished_output(iw, logger, print_section, "AO_MATRICES/DENSITY")
681 : END IF
682 :
683 : ! The xTB matrix itself
684 11678 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, &
685 : "AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)) THEN
686 : iw = cp_print_key_unit_nr(logger, print_section, "AO_MATRICES/KOHN_SHAM_MATRIX", &
687 50 : extension=".Log")
688 50 : CALL section_vals_val_get(print_section, "AO_MATRICES%NDIGITS", i_val=after)
689 50 : after = MIN(MAX(after, 1), 16)
690 100 : DO ispin = 1, dft_control%nspins
691 150 : DO img = 1, SIZE(matrix_ks, 2)
692 : CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, img)%matrix, 4, after, qs_env, para_env, &
693 100 : output_unit=iw, omit_headers=omit_headers)
694 : END DO
695 : END DO
696 50 : CALL cp_print_key_finished_output(iw, logger, print_section, "AO_MATRICES/KOHN_SHAM_MATRIX")
697 : END IF
698 :
699 : ! these print keys are not supported in TB
700 :
701 : ! V_XC CUBE FILE
702 11678 : print_key => section_vals_get_subs_vals(print_section, "V_XC_CUBE")
703 11678 : CALL section_vals_get(print_key, explicit=explicit)
704 11678 : IF (explicit) THEN
705 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
706 0 : CPWARN("XC potential cube file not available for TB methods.")
707 : END IF
708 : END IF
709 :
710 : ! Electric field gradients
711 11678 : print_key => section_vals_get_subs_vals(print_section, "ELECTRIC_FIELD_GRADIENT")
712 11678 : CALL section_vals_get(print_key, explicit=explicit)
713 11678 : IF (explicit) THEN
714 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
715 0 : CPWARN("Electric field gradient not implemented for TB methods.")
716 : END IF
717 : END IF
718 :
719 : ! KINETIC ENERGY
720 11678 : print_key => section_vals_get_subs_vals(print_section, "KINETIC_ENERGY")
721 11678 : CALL section_vals_get(print_key, explicit=explicit)
722 11678 : IF (explicit) THEN
723 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
724 0 : CPWARN("Kinetic energy not available for TB methods.")
725 : END IF
726 : END IF
727 :
728 : ! Xray diffraction spectrum
729 11678 : print_key => section_vals_get_subs_vals(print_section, "XRAY_DIFFRACTION_SPECTRUM")
730 11678 : CALL section_vals_get(print_key, explicit=explicit)
731 11678 : IF (explicit) THEN
732 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
733 0 : CPWARN("Xray diffraction spectrum not implemented for TB methods.")
734 : END IF
735 : END IF
736 :
737 : ! EPR Hyperfine Coupling
738 11678 : print_key => section_vals_get_subs_vals(print_section, "HYPERFINE_COUPLING_TENSOR")
739 11678 : CALL section_vals_get(print_key, explicit=explicit)
740 11678 : IF (explicit) THEN
741 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
742 0 : CPWARN("Hyperfine Coupling not implemented for TB methods.")
743 : END IF
744 : END IF
745 :
746 : ! PLUS_U
747 11678 : print_key => section_vals_get_subs_vals(print_section, "PLUS_U")
748 11678 : CALL section_vals_get(print_key, explicit=explicit)
749 11678 : IF (explicit) THEN
750 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
751 0 : CPWARN("DFT+U method not implemented for TB methods.")
752 : END IF
753 : END IF
754 :
755 11678 : CALL write_ks_matrix_csr(qs_env, qs_env%input)
756 11678 : CALL write_s_matrix_csr(qs_env, qs_env%input)
757 11678 : CALL write_hcore_matrix_csr(qs_env, qs_env%input)
758 11678 : CALL write_p_matrix_csr(qs_env, qs_env%input)
759 :
760 11678 : DEALLOCATE (zcharge)
761 :
762 11678 : CALL timestop(handle)
763 :
764 140136 : END SUBROUTINE scf_post_calculation_tb
765 :
766 : ! **************************************************************************************************
767 : !> \brief ...
768 : !> \param qs_env ...
769 : !> \param input ...
770 : !> \param unit_nr ...
771 : !> \param charges ...
772 : ! **************************************************************************************************
773 1038 : SUBROUTINE tb_dipole(qs_env, input, unit_nr, charges)
774 :
775 : TYPE(qs_environment_type), POINTER :: qs_env
776 : TYPE(section_vals_type), POINTER :: input
777 : INTEGER, INTENT(in) :: unit_nr
778 : REAL(KIND=dp), DIMENSION(:), INTENT(in) :: charges
779 :
780 : CHARACTER(LEN=default_string_length) :: description, dipole_type
781 : COMPLEX(KIND=dp) :: dzeta, dzphase(3), zeta, zphase(3)
782 : COMPLEX(KIND=dp), DIMENSION(3) :: dggamma, ggamma
783 : INTEGER :: i, iat, ikind, j, nat, reference
784 : LOGICAL :: do_berry
785 : REAL(KIND=dp) :: charge_tot, ci(3), dci(3), dipole(3), dipole_deriv(3), drcc(3), dria(3), &
786 : dtheta, gvec(3), q, rcc(3), ria(3), theta, tmp(3), via(3)
787 1038 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
788 1038 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
789 : TYPE(cell_type), POINTER :: cell
790 : TYPE(cp_result_type), POINTER :: results
791 1038 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
792 :
793 1038 : NULLIFY (atomic_kind_set, cell, results)
794 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, &
795 1038 : particle_set=particle_set, cell=cell, results=results)
796 :
797 : ! Reference point
798 1038 : reference = section_get_ival(input, keyword_name="REFERENCE")
799 1038 : NULLIFY (ref_point)
800 1038 : description = '[DIPOLE]'
801 1038 : CALL section_vals_val_get(input, "REF_POINT", r_vals=ref_point)
802 1038 : CALL section_vals_val_get(input, "PERIODIC", l_val=do_berry)
803 :
804 1038 : CALL get_reference_point(rcc, drcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
805 :
806 : ! Dipole deriv will be the derivative of the Dipole(dM/dt=\sum e_j v_j)
807 1038 : dipole_deriv = 0.0_dp
808 1038 : dipole = 0.0_dp
809 1038 : IF (do_berry) THEN
810 640 : dipole_type = "periodic (Berry phase)"
811 2560 : rcc = pbc(rcc, cell)
812 640 : charge_tot = 0._dp
813 4074 : charge_tot = SUM(charges)
814 10240 : ria = twopi*MATMUL(cell%h_inv, rcc)
815 2560 : zphase = CMPLX(COS(ria), SIN(ria), dp)**charge_tot
816 :
817 10240 : dria = twopi*MATMUL(cell%h_inv, drcc)
818 2560 : dzphase = charge_tot*CMPLX(-SIN(ria), COS(ria), dp)**(charge_tot - 1.0_dp)*dria
819 :
820 2560 : ggamma = z_one
821 640 : dggamma = z_zero
822 2152 : DO ikind = 1, SIZE(atomic_kind_set)
823 1512 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
824 5586 : DO i = 1, nat
825 3434 : iat = atomic_kind_set(ikind)%atom_list(i)
826 13736 : ria = particle_set(iat)%r(:)
827 13736 : ria = pbc(ria, cell)
828 13736 : via = particle_set(iat)%v(:)
829 3434 : q = charges(iat)
830 15248 : DO j = 1, 3
831 41208 : gvec = twopi*cell%h_inv(j, :)
832 41208 : theta = SUM(ria(:)*gvec(:))
833 41208 : dtheta = SUM(via(:)*gvec(:))
834 10302 : zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(-q)
835 10302 : dzeta = -q*CMPLX(-SIN(theta), COS(theta), KIND=dp)**(-q - 1.0_dp)*dtheta
836 10302 : dggamma(j) = dggamma(j)*zeta + ggamma(j)*dzeta
837 13736 : ggamma(j) = ggamma(j)*zeta
838 : END DO
839 : END DO
840 : END DO
841 2560 : dggamma = dggamma*zphase + ggamma*dzphase
842 2560 : ggamma = ggamma*zphase
843 2560 : IF (ALL(REAL(ggamma, KIND=dp) /= 0.0_dp)) THEN
844 2560 : tmp = AIMAG(ggamma)/REAL(ggamma, KIND=dp)
845 2560 : ci = -ATAN(tmp)
846 : dci = -(1.0_dp/(1.0_dp + tmp**2))* &
847 2560 : (AIMAG(dggamma)*REAL(ggamma, KIND=dp) - AIMAG(ggamma)*REAL(dggamma, KIND=dp))/(REAL(ggamma, KIND=dp))**2
848 10240 : dipole = MATMUL(cell%hmat, ci)/twopi
849 10240 : dipole_deriv = MATMUL(cell%hmat, dci)/twopi
850 : END IF
851 : ELSE
852 398 : dipole_type = "non-periodic"
853 1864 : DO i = 1, SIZE(particle_set)
854 : ! no pbc(particle_set(i)%r(:),cell) so that the total dipole is the sum of the molecular dipoles
855 5864 : ria = particle_set(i)%r(:)
856 1466 : q = charges(i)
857 5864 : dipole = dipole + q*(ria - rcc)
858 6262 : dipole_deriv(:) = dipole_deriv(:) + q*(particle_set(i)%v(:) - drcc)
859 : END DO
860 : END IF
861 1038 : CALL cp_results_erase(results=results, description=description)
862 : CALL put_results(results=results, description=description, &
863 1038 : values=dipole(1:3))
864 1038 : IF (unit_nr > 0) THEN
865 : WRITE (unit_nr, '(/,T2,A,T31,A50)') &
866 559 : 'TB_DIPOLE| Dipole type', ADJUSTR(TRIM(dipole_type))
867 559 : WRITE (unit_nr, "(T2,A,T30,3(1X,F16.8))") "TB_DIPOLE| Ref. Point [Bohr]", rcc
868 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
869 559 : 'TB_DIPOLE| Moment [a.u.]', dipole(1:3)
870 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
871 2236 : 'TB_DIPOLE| Moment [Debye]', dipole(1:3)*debye
872 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
873 559 : 'TB_DIPOLE| Derivative [a.u.]', dipole_deriv(1:3)
874 : END IF
875 :
876 1038 : END SUBROUTINE tb_dipole
877 :
878 : ! **************************************************************************************************
879 : !> \brief computes the MOs and calls the wavefunction mixing routine.
880 : !> \param qs_env ...
881 : !> \param dft_section ...
882 : !> \param scf_env ...
883 : !> \author Florian Schiffmann
884 : !> \note
885 : ! **************************************************************************************************
886 :
887 2 : SUBROUTINE wfn_mix_tb(qs_env, dft_section, scf_env)
888 :
889 : TYPE(qs_environment_type), POINTER :: qs_env
890 : TYPE(section_vals_type), POINTER :: dft_section
891 : TYPE(qs_scf_env_type), POINTER :: scf_env
892 :
893 : INTEGER :: ispin, nao, nmo, output_unit
894 2 : REAL(dp), DIMENSION(:), POINTER :: mo_eigenvalues
895 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
896 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct, ao_lumo_struct
897 : TYPE(cp_fm_type) :: KS_tmp, MO_tmp, S_tmp, work
898 2 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: lumos
899 : TYPE(cp_fm_type), POINTER :: mo_coeff
900 : TYPE(cp_logger_type), POINTER :: logger
901 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
902 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
903 : TYPE(mp_para_env_type), POINTER :: para_env
904 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
905 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
906 : TYPE(section_vals_type), POINTER :: wfn_mix_section
907 :
908 4 : logger => cp_get_default_logger()
909 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks, &
910 : particle_set=particle_set, atomic_kind_set=atomic_kind_set, &
911 2 : qs_kind_set=qs_kind_set, mos=mos, para_env=para_env)
912 :
913 2 : wfn_mix_section => section_vals_get_subs_vals(dft_section, "PRINT%WFN_MIX")
914 :
915 2 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff, nao=nao)
916 :
917 : CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, nrow_global=nao, ncol_global=nao, &
918 2 : template_fmstruct=mo_coeff%matrix_struct)
919 2 : CALL cp_fm_create(S_tmp, matrix_struct=ao_ao_fmstruct)
920 2 : CALL cp_fm_create(KS_tmp, matrix_struct=ao_ao_fmstruct)
921 2 : CALL cp_fm_create(MO_tmp, matrix_struct=ao_ao_fmstruct)
922 2 : CALL cp_fm_create(work, matrix_struct=ao_ao_fmstruct)
923 10 : ALLOCATE (lumos(SIZE(mos)))
924 :
925 2 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_tmp)
926 2 : CALL cp_fm_cholesky_decompose(S_tmp)
927 :
928 6 : DO ispin = 1, SIZE(mos)
929 4 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues, nmo=nmo)
930 : CALL cp_fm_struct_create(fmstruct=ao_lumo_struct, nrow_global=nao, ncol_global=nao - nmo, &
931 4 : template_fmstruct=mo_coeff%matrix_struct)
932 :
933 4 : CALL cp_fm_create(lumos(ispin), matrix_struct=ao_lumo_struct)
934 4 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, KS_tmp)
935 4 : CALL cp_fm_cholesky_reduce(KS_tmp, S_tmp)
936 4 : CALL choose_eigv_solver(KS_tmp, work, mo_eigenvalues)
937 4 : CALL cp_fm_cholesky_restore(work, nao, S_tmp, MO_tmp, "SOLVE")
938 4 : CALL cp_fm_to_fm_submat(MO_tmp, mo_coeff, nao, nmo, 1, 1, 1, 1)
939 4 : CALL cp_fm_to_fm_submat(MO_tmp, lumos(ispin), nao, nao - nmo, 1, nmo + 1, 1, 1)
940 :
941 10 : CALL cp_fm_struct_release(ao_lumo_struct)
942 : END DO
943 :
944 2 : output_unit = cp_logger_get_default_io_unit(logger)
945 : CALL wfn_mix(mos, particle_set, dft_section, qs_kind_set, para_env, output_unit, &
946 2 : unoccupied_orbs=lumos, scf_env=scf_env, matrix_s=matrix_s)
947 :
948 2 : CALL cp_fm_release(lumos)
949 2 : CALL cp_fm_release(S_tmp)
950 2 : CALL cp_fm_release(MO_tmp)
951 2 : CALL cp_fm_release(KS_tmp)
952 2 : CALL cp_fm_release(work)
953 2 : CALL cp_fm_struct_release(ao_ao_fmstruct)
954 :
955 6 : END SUBROUTINE wfn_mix_tb
956 :
957 : ! **************************************************************************************************
958 : !> \brief Gets the lumos, and eigenvalues for the lumos
959 : !> \param qs_env ...
960 : !> \param scf_env ...
961 : !> \param unoccupied_orbs ...
962 : !> \param unoccupied_evals ...
963 : !> \param nlumo ...
964 : !> \param nlumos ...
965 : ! **************************************************************************************************
966 2 : SUBROUTINE make_lumo_tb(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, nlumo, nlumos)
967 :
968 : TYPE(qs_environment_type), POINTER :: qs_env
969 : TYPE(qs_scf_env_type), POINTER :: scf_env
970 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: unoccupied_orbs
971 : TYPE(cp_1d_r_p_type), DIMENSION(:), INTENT(INOUT) :: unoccupied_evals
972 : INTEGER :: nlumo
973 : INTEGER, INTENT(OUT) :: nlumos
974 :
975 : INTEGER :: homo, iounit, ispin, n, nao, nmo
976 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
977 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
978 : TYPE(cp_fm_type), POINTER :: mo_coeff
979 : TYPE(cp_logger_type), POINTER :: logger
980 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, matrix_s
981 : TYPE(dft_control_type), POINTER :: dft_control
982 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
983 : TYPE(mp_para_env_type), POINTER :: para_env
984 : TYPE(preconditioner_type), POINTER :: local_preconditioner
985 : TYPE(scf_control_type), POINTER :: scf_control
986 :
987 2 : NULLIFY (mos, ks_rmpv, scf_control, dft_control, para_env, blacs_env)
988 : CALL get_qs_env(qs_env, &
989 : mos=mos, &
990 : matrix_ks=ks_rmpv, &
991 : scf_control=scf_control, &
992 : dft_control=dft_control, &
993 : matrix_s=matrix_s, &
994 : para_env=para_env, &
995 2 : blacs_env=blacs_env)
996 :
997 2 : logger => cp_get_default_logger()
998 2 : iounit = cp_logger_get_default_io_unit(logger)
999 :
1000 4 : DO ispin = 1, dft_control%nspins
1001 2 : NULLIFY (unoccupied_evals(ispin)%array)
1002 : ! Always write eigenvalues
1003 2 : IF (iounit > 0) WRITE (iounit, *) " "
1004 2 : IF (iounit > 0) WRITE (iounit, *) " Lowest Eigenvalues of the unoccupied subspace spin ", ispin
1005 2 : IF (iounit > 0) WRITE (iounit, FMT='(1X,A)') "-----------------------------------------------------"
1006 2 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=homo, nao=nao, nmo=nmo)
1007 2 : CALL cp_fm_get_info(mo_coeff, nrow_global=n)
1008 2 : nlumos = MAX(1, MIN(nlumo, nao - nmo))
1009 2 : IF (nlumo == -1) nlumos = nao - nmo
1010 6 : ALLOCATE (unoccupied_evals(ispin)%array(nlumos))
1011 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
1012 2 : nrow_global=n, ncol_global=nlumos)
1013 2 : CALL cp_fm_create(unoccupied_orbs(ispin), fm_struct_tmp, name="lumos")
1014 2 : CALL cp_fm_struct_release(fm_struct_tmp)
1015 2 : CALL cp_fm_init_random(unoccupied_orbs(ispin), nlumos)
1016 :
1017 : ! FULL_ALL has column-dependent occupied-state weights and cannot be reused here.
1018 2 : NULLIFY (local_preconditioner)
1019 2 : IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
1020 2 : local_preconditioner => scf_env%ot_preconditioner(1)%preconditioner
1021 2 : IF (local_preconditioner%in_use == ot_precond_full_all .OR. &
1022 : local_preconditioner%in_use == ot_precond_full_all_covariant) THEN
1023 2 : NULLIFY (local_preconditioner)
1024 : END IF
1025 : END IF
1026 :
1027 : CALL ot_eigensolver(matrix_h=ks_rmpv(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
1028 : matrix_c_fm=unoccupied_orbs(ispin), &
1029 : matrix_orthogonal_space_fm=mo_coeff, &
1030 : eps_gradient=scf_control%eps_lumos, &
1031 : preconditioner=local_preconditioner, &
1032 : iter_max=scf_control%max_iter_lumos, &
1033 2 : size_ortho_space=nmo)
1034 :
1035 : CALL calculate_subspace_eigenvalues(unoccupied_orbs(ispin), ks_rmpv(ispin)%matrix, &
1036 : unoccupied_evals(ispin)%array, scr=iounit, &
1037 6 : ionode=iounit > 0)
1038 :
1039 : END DO
1040 :
1041 2 : END SUBROUTINE make_lumo_tb
1042 :
1043 : ! **************************************************************************************************
1044 : !> \brief ...
1045 : !> \param qs_env ...
1046 : ! **************************************************************************************************
1047 34 : SUBROUTINE rebuild_pw_env(qs_env)
1048 :
1049 : TYPE(qs_environment_type), POINTER :: qs_env
1050 :
1051 : LOGICAL :: skip_load_balance_distributed
1052 : TYPE(cell_type), POINTER :: cell
1053 : TYPE(dft_control_type), POINTER :: dft_control
1054 : TYPE(pw_env_type), POINTER :: new_pw_env
1055 : TYPE(qs_ks_env_type), POINTER :: ks_env
1056 : TYPE(qs_rho_type), POINTER :: rho
1057 : TYPE(task_list_type), POINTER :: task_list
1058 :
1059 34 : CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=new_pw_env)
1060 34 : IF (.NOT. ASSOCIATED(new_pw_env)) THEN
1061 0 : CALL pw_env_create(new_pw_env)
1062 0 : CALL set_ks_env(ks_env, pw_env=new_pw_env)
1063 0 : CALL pw_env_release(new_pw_env)
1064 : END IF
1065 34 : CALL get_qs_env(qs_env, pw_env=new_pw_env, dft_control=dft_control, cell=cell)
1066 :
1067 884 : new_pw_env%cell_hmat = cell%hmat
1068 34 : CALL pw_env_rebuild(new_pw_env, qs_env=qs_env)
1069 :
1070 34 : NULLIFY (task_list)
1071 34 : CALL get_ks_env(ks_env, task_list=task_list)
1072 34 : IF (.NOT. ASSOCIATED(task_list)) THEN
1073 34 : CALL allocate_task_list(task_list)
1074 34 : CALL set_ks_env(ks_env, task_list=task_list)
1075 : END IF
1076 34 : skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
1077 : CALL generate_qs_task_list(ks_env, task_list, basis_type="ORB", &
1078 : reorder_rs_grid_ranks=.TRUE., &
1079 34 : skip_load_balance_distributed=skip_load_balance_distributed)
1080 34 : CALL get_qs_env(qs_env, rho=rho)
1081 34 : CALL qs_rho_rebuild(rho, qs_env=qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.)
1082 :
1083 34 : END SUBROUTINE rebuild_pw_env
1084 :
1085 : ! **************************************************************************************************
1086 : !> \brief ...
1087 : !> \param qs_env ...
1088 : !> \param zcharge ...
1089 : !> \param cube_section ...
1090 : ! **************************************************************************************************
1091 24 : SUBROUTINE print_e_density(qs_env, zcharge, cube_section)
1092 :
1093 : TYPE(qs_environment_type), POINTER :: qs_env
1094 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zcharge
1095 : TYPE(section_vals_type), POINTER :: cube_section
1096 :
1097 : CHARACTER(LEN=default_path_length) :: filename, mpi_filename, my_pos_cube
1098 : INTEGER :: iounit, ispin, unit_nr
1099 : LOGICAL :: append_cube, mpi_io
1100 24 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
1101 : TYPE(cp_logger_type), POINTER :: logger
1102 24 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
1103 24 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1104 : TYPE(dft_control_type), POINTER :: dft_control
1105 : TYPE(particle_list_type), POINTER :: particles
1106 24 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1107 : TYPE(pw_env_type), POINTER :: pw_env
1108 24 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1109 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1110 24 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1111 : TYPE(qs_ks_env_type), POINTER :: ks_env
1112 : TYPE(qs_rho_type), POINTER :: rho
1113 : TYPE(qs_subsys_type), POINTER :: subsys
1114 :
1115 24 : CALL get_qs_env(qs_env, dft_control=dft_control)
1116 :
1117 24 : append_cube = section_get_lval(cube_section, "APPEND")
1118 24 : my_pos_cube = "REWIND"
1119 24 : IF (append_cube) my_pos_cube = "APPEND"
1120 :
1121 24 : logger => cp_get_default_logger()
1122 24 : iounit = cp_logger_get_default_io_unit(logger)
1123 :
1124 : ! we need to construct the density on a realspace grid
1125 24 : CALL get_qs_env(qs_env, ks_env=ks_env, rho=rho)
1126 24 : NULLIFY (rho_r, rho_g, tot_rho_r)
1127 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, &
1128 24 : rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
1129 50 : DO ispin = 1, dft_control%nspins
1130 26 : rho_ao => rho_ao_kp(ispin, :)
1131 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
1132 : rho=rho_r(ispin), &
1133 : rho_gspace=rho_g(ispin), &
1134 : total_rho=tot_rho_r(ispin), &
1135 50 : ks_env=ks_env)
1136 : END DO
1137 24 : CALL qs_rho_set(rho, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
1138 :
1139 24 : CALL get_qs_env(qs_env, subsys=subsys)
1140 24 : CALL qs_subsys_get(subsys, particles=particles)
1141 :
1142 24 : IF (dft_control%nspins > 1) THEN
1143 2 : IF (iounit > 0) THEN
1144 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,2F15.6)") &
1145 3 : "Integrated alpha and beta electronic density:", tot_rho_r(1:2)
1146 : END IF
1147 2 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1148 2 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1149 : BLOCK
1150 : TYPE(pw_r3d_rs_type) :: rho_elec_rspace
1151 2 : CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
1152 2 : CALL pw_copy(rho_r(1), rho_elec_rspace)
1153 2 : CALL pw_axpy(rho_r(2), rho_elec_rspace)
1154 2 : filename = "ELECTRON_DENSITY"
1155 2 : mpi_io = .TRUE.
1156 : unit_nr = cp_print_key_unit_nr(logger, cube_section, '', &
1157 : extension=".cube", middle_name=TRIM(filename), &
1158 : file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, &
1159 2 : fout=mpi_filename)
1160 2 : IF (iounit > 0) THEN
1161 1 : IF (.NOT. mpi_io) THEN
1162 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
1163 : ELSE
1164 1 : filename = mpi_filename
1165 : END IF
1166 : WRITE (UNIT=iounit, FMT="(T2,A,/,T2,A79)") &
1167 1 : "The sum of alpha and beta density is written in cube file format to the file:", ADJUSTR(TRIM(filename))
1168 : END IF
1169 : CALL cp_pw_to_cube(rho_elec_rspace, unit_nr, "SUM OF ALPHA AND BETA DENSITY", &
1170 : particles=particles, zeff=zcharge, stride=section_get_ivals(cube_section, "STRIDE"), &
1171 2 : mpi_io=mpi_io)
1172 2 : CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1173 2 : CALL pw_copy(rho_r(1), rho_elec_rspace)
1174 2 : CALL pw_axpy(rho_r(2), rho_elec_rspace, alpha=-1.0_dp)
1175 2 : filename = "SPIN_DENSITY"
1176 2 : mpi_io = .TRUE.
1177 : unit_nr = cp_print_key_unit_nr(logger, cube_section, '', &
1178 : extension=".cube", middle_name=TRIM(filename), &
1179 : file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, &
1180 2 : fout=mpi_filename)
1181 2 : IF (iounit > 0) THEN
1182 1 : IF (.NOT. mpi_io) THEN
1183 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
1184 : ELSE
1185 1 : filename = mpi_filename
1186 : END IF
1187 : WRITE (UNIT=iounit, FMT="(T2,A,/,T2,A79)") &
1188 1 : "The spin density is written in cube file format to the file:", ADJUSTR(TRIM(filename))
1189 : END IF
1190 : CALL cp_pw_to_cube(rho_elec_rspace, unit_nr, "SPIN DENSITY", &
1191 : particles=particles, zeff=zcharge, &
1192 2 : stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1193 2 : CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1194 2 : CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
1195 : END BLOCK
1196 : ELSE
1197 22 : IF (iounit > 0) THEN
1198 : WRITE (UNIT=iounit, FMT="(/,T2,A,T66,F15.6)") &
1199 11 : "Integrated electronic density:", tot_rho_r(1)
1200 : END IF
1201 22 : filename = "ELECTRON_DENSITY"
1202 22 : mpi_io = .TRUE.
1203 : unit_nr = cp_print_key_unit_nr(logger, cube_section, '', &
1204 : extension=".cube", middle_name=TRIM(filename), &
1205 : file_position=my_pos_cube, log_filename=.FALSE., mpi_io=mpi_io, &
1206 22 : fout=mpi_filename)
1207 22 : IF (iounit > 0) THEN
1208 11 : IF (.NOT. mpi_io) THEN
1209 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
1210 : ELSE
1211 11 : filename = mpi_filename
1212 : END IF
1213 : WRITE (UNIT=iounit, FMT="(T2,A,/,T2,A79)") &
1214 11 : "The electron density is written in cube file format to the file:", ADJUSTR(TRIM(filename))
1215 : END IF
1216 : CALL cp_pw_to_cube(rho_r(1), unit_nr, "ELECTRON DENSITY", &
1217 : particles=particles, zeff=zcharge, &
1218 22 : stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1219 22 : CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1220 : END IF ! nspins
1221 :
1222 24 : END SUBROUTINE print_e_density
1223 : ! **************************************************************************************************
1224 : !> \brief ...
1225 : !> \param qs_env ...
1226 : !> \param zcharge ...
1227 : !> \param cube_section ...
1228 : !> \param total_density ...
1229 : !> \param v_hartree ...
1230 : !> \param efield ...
1231 : ! **************************************************************************************************
1232 74 : SUBROUTINE print_density_cubes(qs_env, zcharge, cube_section, total_density, v_hartree, efield)
1233 :
1234 : TYPE(qs_environment_type), POINTER :: qs_env
1235 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zcharge
1236 : TYPE(section_vals_type), POINTER :: cube_section
1237 : LOGICAL, INTENT(IN), OPTIONAL :: total_density, v_hartree, efield
1238 :
1239 : CHARACTER(len=1), DIMENSION(3), PARAMETER :: cdir = ["x", "y", "z"]
1240 :
1241 : CHARACTER(LEN=default_path_length) :: filename, mpi_filename, my_pos_cube
1242 : INTEGER :: id, iounit, ispin, nd(3), unit_nr
1243 : LOGICAL :: append_cube, mpi_io, my_efield, &
1244 : my_total_density, my_v_hartree
1245 : REAL(KIND=dp) :: total_rho_core_rspace, udvol
1246 74 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
1247 : TYPE(cell_type), POINTER :: cell
1248 : TYPE(cp_logger_type), POINTER :: logger
1249 74 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
1250 74 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1251 : TYPE(dft_control_type), POINTER :: dft_control
1252 : TYPE(particle_list_type), POINTER :: particles
1253 : TYPE(pw_c1d_gs_type) :: rho_core
1254 74 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1255 : TYPE(pw_env_type), POINTER :: pw_env
1256 : TYPE(pw_poisson_parameter_type) :: poisson_params
1257 74 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1258 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1259 : TYPE(pw_r3d_rs_type) :: rho_tot_rspace
1260 74 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1261 : TYPE(qs_ks_env_type), POINTER :: ks_env
1262 : TYPE(qs_rho_type), POINTER :: rho
1263 : TYPE(qs_subsys_type), POINTER :: subsys
1264 :
1265 74 : CALL get_qs_env(qs_env, cell=cell, dft_control=dft_control)
1266 :
1267 74 : append_cube = section_get_lval(cube_section, "APPEND")
1268 74 : my_pos_cube = "REWIND"
1269 74 : IF (append_cube) my_pos_cube = "APPEND"
1270 :
1271 74 : IF (PRESENT(total_density)) THEN
1272 26 : my_total_density = total_density
1273 : ELSE
1274 : my_total_density = .FALSE.
1275 : END IF
1276 74 : IF (PRESENT(v_hartree)) THEN
1277 24 : my_v_hartree = v_hartree
1278 : ELSE
1279 : my_v_hartree = .FALSE.
1280 : END IF
1281 74 : IF (PRESENT(efield)) THEN
1282 24 : my_efield = efield
1283 : ELSE
1284 : my_efield = .FALSE.
1285 : END IF
1286 :
1287 74 : logger => cp_get_default_logger()
1288 74 : iounit = cp_logger_get_default_io_unit(logger)
1289 :
1290 : ! we need to construct the density on a realspace grid
1291 74 : CALL get_qs_env(qs_env, ks_env=ks_env, rho=rho)
1292 74 : NULLIFY (rho_r, rho_g, tot_rho_r)
1293 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, &
1294 74 : rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
1295 150 : DO ispin = 1, dft_control%nspins
1296 76 : rho_ao => rho_ao_kp(ispin, :)
1297 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
1298 : rho=rho_r(ispin), &
1299 : rho_gspace=rho_g(ispin), &
1300 : total_rho=tot_rho_r(ispin), &
1301 150 : ks_env=ks_env)
1302 : END DO
1303 74 : CALL qs_rho_set(rho, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
1304 :
1305 74 : CALL get_qs_env(qs_env, subsys=subsys)
1306 74 : CALL qs_subsys_get(subsys, particles=particles)
1307 :
1308 74 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1309 74 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1310 74 : CALL auxbas_pw_pool%create_pw(pw=rho_core)
1311 74 : CALL calculate_rho_core(rho_core, total_rho_core_rspace, qs_env)
1312 :
1313 74 : IF (iounit > 0) THEN
1314 : WRITE (UNIT=iounit, FMT="(/,T2,A,T66,F15.6)") &
1315 75 : "Integrated electronic density:", SUM(tot_rho_r(:))
1316 : WRITE (UNIT=iounit, FMT="(T2,A,T66,F15.6)") &
1317 37 : "Integrated core density:", total_rho_core_rspace
1318 : END IF
1319 :
1320 74 : CALL auxbas_pw_pool%create_pw(pw=rho_tot_rspace)
1321 74 : CALL pw_transfer(rho_core, rho_tot_rspace)
1322 150 : DO ispin = 1, dft_control%nspins
1323 150 : CALL pw_axpy(rho_r(ispin), rho_tot_rspace)
1324 : END DO
1325 :
1326 74 : IF (my_total_density) THEN
1327 26 : filename = "TOTAL_DENSITY"
1328 26 : mpi_io = .TRUE.
1329 : unit_nr = cp_print_key_unit_nr(logger, cube_section, '', &
1330 : extension=".cube", middle_name=TRIM(filename), file_position=my_pos_cube, &
1331 26 : log_filename=.FALSE., mpi_io=mpi_io, fout=mpi_filename)
1332 26 : IF (iounit > 0) THEN
1333 13 : IF (.NOT. mpi_io) THEN
1334 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
1335 : ELSE
1336 13 : filename = mpi_filename
1337 : END IF
1338 : WRITE (UNIT=iounit, FMT="(T2,A,/,T2,A79)") &
1339 13 : "The total density is written in cube file format to the file:", ADJUSTR(TRIM(filename))
1340 : END IF
1341 : CALL cp_pw_to_cube(rho_tot_rspace, unit_nr, "TOTAL DENSITY", &
1342 : particles=particles, zeff=zcharge, &
1343 26 : stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1344 26 : CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1345 : END IF
1346 74 : IF (my_v_hartree .OR. my_efield) THEN
1347 : BLOCK
1348 : TYPE(pw_c1d_gs_type) :: rho_tot_gspace
1349 48 : CALL auxbas_pw_pool%create_pw(pw=rho_tot_gspace)
1350 48 : CALL pw_transfer(rho_tot_rspace, rho_tot_gspace)
1351 48 : poisson_params%solver = pw_poisson_analytic
1352 192 : poisson_params%periodic = cell%perd
1353 48 : poisson_params%ewald_type = do_ewald_none
1354 96 : BLOCK
1355 48 : TYPE(greens_fn_type) :: green_fft
1356 : TYPE(pw_grid_type), POINTER :: pwdummy
1357 48 : NULLIFY (pwdummy)
1358 48 : CALL pw_green_create(green_fft, poisson_params, cell%hmat, auxbas_pw_pool, pwdummy, pwdummy)
1359 825006 : rho_tot_gspace%array(:) = rho_tot_gspace%array(:)*green_fft%influence_fn%array(:)
1360 96 : CALL pw_green_release(green_fft, auxbas_pw_pool)
1361 : END BLOCK
1362 48 : IF (my_v_hartree) THEN
1363 : BLOCK
1364 : TYPE(pw_r3d_rs_type) :: vhartree
1365 24 : CALL auxbas_pw_pool%create_pw(pw=vhartree)
1366 24 : CALL pw_transfer(rho_tot_gspace, vhartree)
1367 24 : filename = "V_HARTREE"
1368 24 : mpi_io = .TRUE.
1369 : unit_nr = cp_print_key_unit_nr(logger, cube_section, '', &
1370 : extension=".cube", middle_name=TRIM(filename), file_position=my_pos_cube, &
1371 24 : log_filename=.FALSE., mpi_io=mpi_io, fout=mpi_filename)
1372 24 : IF (iounit > 0) THEN
1373 12 : IF (.NOT. mpi_io) THEN
1374 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
1375 : ELSE
1376 12 : filename = mpi_filename
1377 : END IF
1378 : WRITE (UNIT=iounit, FMT="(T2,A,/,T2,A79)") &
1379 12 : "The Hartree potential is written in cube file format to the file:", ADJUSTR(TRIM(filename))
1380 : END IF
1381 : CALL cp_pw_to_cube(vhartree, unit_nr, "Hartree Potential", &
1382 : particles=particles, zeff=zcharge, &
1383 24 : stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1384 24 : CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1385 24 : CALL auxbas_pw_pool%give_back_pw(vhartree)
1386 : END BLOCK
1387 : END IF
1388 48 : IF (my_efield) THEN
1389 : BLOCK
1390 : TYPE(pw_c1d_gs_type) :: vhartree
1391 24 : CALL auxbas_pw_pool%create_pw(pw=vhartree)
1392 24 : udvol = 1.0_dp/rho_tot_rspace%pw_grid%dvol
1393 96 : DO id = 1, 3
1394 72 : CALL pw_transfer(rho_tot_gspace, vhartree)
1395 72 : nd = 0
1396 72 : nd(id) = 1
1397 72 : CALL pw_derive(vhartree, nd)
1398 72 : CALL pw_transfer(vhartree, rho_tot_rspace)
1399 72 : CALL pw_scale(rho_tot_rspace, udvol)
1400 :
1401 72 : filename = "EFIELD_"//cdir(id)
1402 72 : mpi_io = .TRUE.
1403 : unit_nr = cp_print_key_unit_nr(logger, cube_section, '', &
1404 : extension=".cube", middle_name=TRIM(filename), file_position=my_pos_cube, &
1405 72 : log_filename=.FALSE., mpi_io=mpi_io, fout=mpi_filename)
1406 72 : IF (iounit > 0) THEN
1407 36 : IF (.NOT. mpi_io) THEN
1408 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
1409 : ELSE
1410 36 : filename = mpi_filename
1411 : END IF
1412 : WRITE (UNIT=iounit, FMT="(T2,A,/,T2,A79)") &
1413 36 : "The Efield is written in cube file format to the file:", ADJUSTR(TRIM(filename))
1414 : END IF
1415 : CALL cp_pw_to_cube(rho_tot_rspace, unit_nr, "EFIELD "//cdir(id), &
1416 : particles=particles, zeff=zcharge, &
1417 72 : stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1418 96 : CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1419 : END DO
1420 24 : CALL auxbas_pw_pool%give_back_pw(vhartree)
1421 : END BLOCK
1422 : END IF
1423 48 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1424 : END BLOCK
1425 : END IF
1426 :
1427 74 : CALL auxbas_pw_pool%give_back_pw(rho_tot_rspace)
1428 74 : CALL auxbas_pw_pool%give_back_pw(rho_core)
1429 :
1430 296 : END SUBROUTINE print_density_cubes
1431 :
1432 : ! **************************************************************************************************
1433 : !> \brief ...
1434 : !> \param qs_env ...
1435 : !> \param zcharge ...
1436 : !> \param elf_section ...
1437 : ! **************************************************************************************************
1438 24 : SUBROUTINE print_elf(qs_env, zcharge, elf_section)
1439 :
1440 : TYPE(qs_environment_type), POINTER :: qs_env
1441 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zcharge
1442 : TYPE(section_vals_type), POINTER :: elf_section
1443 :
1444 : CHARACTER(LEN=default_path_length) :: filename, mpi_filename, my_pos_cube, &
1445 : title
1446 : INTEGER :: iounit, ispin, unit_nr
1447 : LOGICAL :: append_cube, mpi_io
1448 : REAL(KIND=dp) :: rho_cutoff
1449 24 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
1450 : TYPE(cp_logger_type), POINTER :: logger
1451 24 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
1452 24 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
1453 : TYPE(dft_control_type), POINTER :: dft_control
1454 : TYPE(particle_list_type), POINTER :: particles
1455 24 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1456 : TYPE(pw_env_type), POINTER :: pw_env
1457 24 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1458 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1459 24 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: elf_r
1460 24 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1461 : TYPE(qs_ks_env_type), POINTER :: ks_env
1462 : TYPE(qs_rho_type), POINTER :: rho
1463 : TYPE(qs_subsys_type), POINTER :: subsys
1464 :
1465 48 : logger => cp_get_default_logger()
1466 24 : iounit = cp_logger_get_default_io_unit(logger)
1467 :
1468 : ! we need to construct the density on a realspace grid
1469 24 : CALL get_qs_env(qs_env, dft_control=dft_control, ks_env=ks_env, rho=rho)
1470 24 : NULLIFY (rho_r, rho_g, tot_rho_r)
1471 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, &
1472 24 : rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
1473 50 : DO ispin = 1, dft_control%nspins
1474 26 : rho_ao => rho_ao_kp(ispin, :)
1475 : CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
1476 : rho=rho_r(ispin), &
1477 : rho_gspace=rho_g(ispin), &
1478 : total_rho=tot_rho_r(ispin), &
1479 50 : ks_env=ks_env)
1480 : END DO
1481 24 : CALL qs_rho_set(rho, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
1482 :
1483 24 : CALL get_qs_env(qs_env, subsys=subsys)
1484 24 : CALL qs_subsys_get(subsys, particles=particles)
1485 :
1486 98 : ALLOCATE (elf_r(dft_control%nspins))
1487 24 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1488 24 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1489 50 : DO ispin = 1, dft_control%nspins
1490 26 : CALL auxbas_pw_pool%create_pw(elf_r(ispin))
1491 50 : CALL pw_zero(elf_r(ispin))
1492 : END DO
1493 :
1494 24 : IF (iounit > 0) THEN
1495 : WRITE (UNIT=iounit, FMT="(/,T2,A)") &
1496 12 : "ELF is computed on the real space grid -----"
1497 : END IF
1498 24 : rho_cutoff = section_get_rval(elf_section, "density_cutoff")
1499 24 : CALL qs_elf_calc(qs_env, elf_r, rho_cutoff)
1500 :
1501 : ! write ELF into cube file
1502 24 : append_cube = section_get_lval(elf_section, "APPEND")
1503 24 : my_pos_cube = "REWIND"
1504 24 : IF (append_cube) my_pos_cube = "APPEND"
1505 50 : DO ispin = 1, dft_control%nspins
1506 26 : WRITE (filename, '(a5,I1.1)') "ELF_S", ispin
1507 26 : WRITE (title, *) "ELF spin ", ispin
1508 26 : mpi_io = .TRUE.
1509 : unit_nr = cp_print_key_unit_nr(logger, elf_section, '', extension=".cube", &
1510 : middle_name=TRIM(filename), file_position=my_pos_cube, &
1511 26 : log_filename=.FALSE., mpi_io=mpi_io, fout=mpi_filename)
1512 26 : IF (iounit > 0) THEN
1513 13 : IF (.NOT. mpi_io) THEN
1514 0 : INQUIRE (UNIT=unit_nr, NAME=filename)
1515 : ELSE
1516 13 : filename = mpi_filename
1517 : END IF
1518 : WRITE (UNIT=iounit, FMT="(T2,A,/,T2,A79)") &
1519 13 : "ELF is written in cube file format to the file:", ADJUSTR(TRIM(filename))
1520 : END IF
1521 :
1522 : CALL cp_pw_to_cube(elf_r(ispin), unit_nr, title, particles=particles, zeff=zcharge, &
1523 26 : stride=section_get_ivals(elf_section, "STRIDE"), mpi_io=mpi_io)
1524 26 : CALL cp_print_key_finished_output(unit_nr, logger, elf_section, '', mpi_io=mpi_io)
1525 :
1526 50 : CALL auxbas_pw_pool%give_back_pw(elf_r(ispin))
1527 : END DO
1528 :
1529 24 : DEALLOCATE (elf_r)
1530 :
1531 24 : END SUBROUTINE print_elf
1532 : ! **************************************************************************************************
1533 : !> \brief ...
1534 : !> \param qs_env ...
1535 : !> \param zcharge ...
1536 : !> \param cube_section ...
1537 : ! **************************************************************************************************
1538 24 : SUBROUTINE print_mo_cubes(qs_env, zcharge, cube_section)
1539 :
1540 : TYPE(qs_environment_type), POINTER :: qs_env
1541 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zcharge
1542 : TYPE(section_vals_type), POINTER :: cube_section
1543 :
1544 : CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
1545 : INTEGER :: homo, i, ifirst, ilast, iounit, ir, &
1546 : ispin, ivector, n_rep, nhomo, nlist, &
1547 : nlumo, nmo, shomo, unit_nr
1548 24 : INTEGER, DIMENSION(:), POINTER :: list, list_index
1549 : LOGICAL :: append_cube, mpi_io, write_cube
1550 : REAL(KIND=dp) :: homo_lumo(2, 2)
1551 24 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
1552 24 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1553 : TYPE(cell_type), POINTER :: cell
1554 : TYPE(cp_fm_type), POINTER :: mo_coeff
1555 : TYPE(cp_logger_type), POINTER :: logger
1556 24 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, mo_derivs
1557 : TYPE(dft_control_type), POINTER :: dft_control
1558 24 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1559 : TYPE(particle_list_type), POINTER :: particles
1560 24 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1561 : TYPE(pw_c1d_gs_type) :: wf_g
1562 : TYPE(pw_env_type), POINTER :: pw_env
1563 24 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1564 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1565 : TYPE(pw_r3d_rs_type) :: wf_r
1566 24 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1567 : TYPE(qs_subsys_type), POINTER :: subsys
1568 : TYPE(scf_control_type), POINTER :: scf_control
1569 :
1570 48 : logger => cp_get_default_logger()
1571 24 : iounit = cp_logger_get_default_io_unit(logger)
1572 :
1573 24 : CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv, scf_control=scf_control)
1574 24 : CALL get_qs_env(qs_env, dft_control=dft_control, mo_derivs=mo_derivs)
1575 24 : CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs)
1576 24 : NULLIFY (mo_eigenvalues)
1577 24 : homo = 0
1578 50 : DO ispin = 1, dft_control%nspins
1579 26 : CALL get_mo_set(mo_set=mos(ispin), eigenvalues=mo_eigenvalues, homo=shomo)
1580 26 : homo_lumo(ispin, 1) = mo_eigenvalues(shomo)
1581 50 : homo = MAX(homo, shomo)
1582 : END DO
1583 24 : write_cube = section_get_lval(cube_section, "WRITE_CUBE")
1584 24 : nlumo = section_get_ival(cube_section, "NLUMO")
1585 24 : nhomo = section_get_ival(cube_section, "NHOMO")
1586 24 : NULLIFY (list_index)
1587 24 : CALL section_vals_val_get(cube_section, "HOMO_LIST", n_rep_val=n_rep)
1588 24 : IF (n_rep > 0) THEN
1589 2 : nlist = 0
1590 4 : DO ir = 1, n_rep
1591 2 : NULLIFY (list)
1592 2 : CALL section_vals_val_get(cube_section, "HOMO_LIST", i_rep_val=ir, i_vals=list)
1593 4 : IF (ASSOCIATED(list)) THEN
1594 2 : CALL reallocate(list_index, 1, nlist + SIZE(list))
1595 14 : DO i = 1, SIZE(list)
1596 14 : list_index(i + nlist) = list(i)
1597 : END DO
1598 2 : nlist = nlist + SIZE(list)
1599 : END IF
1600 : END DO
1601 14 : nhomo = MAXVAL(list_index)
1602 : ELSE
1603 22 : IF (nhomo == -1) nhomo = homo
1604 22 : nlist = homo - MAX(1, homo - nhomo + 1) + 1
1605 66 : ALLOCATE (list_index(nlist))
1606 44 : DO i = 1, nlist
1607 44 : list_index(i) = MAX(1, homo - nhomo + 1) + i - 1
1608 : END DO
1609 : END IF
1610 :
1611 24 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1612 24 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1613 24 : CALL auxbas_pw_pool%create_pw(wf_r)
1614 24 : CALL auxbas_pw_pool%create_pw(wf_g)
1615 :
1616 24 : CALL get_qs_env(qs_env, subsys=subsys)
1617 24 : CALL qs_subsys_get(subsys, particles=particles)
1618 :
1619 24 : append_cube = section_get_lval(cube_section, "APPEND")
1620 24 : my_pos_cube = "REWIND"
1621 24 : IF (append_cube) THEN
1622 0 : my_pos_cube = "APPEND"
1623 : END IF
1624 :
1625 : CALL get_qs_env(qs_env=qs_env, &
1626 : atomic_kind_set=atomic_kind_set, &
1627 : qs_kind_set=qs_kind_set, &
1628 : cell=cell, &
1629 24 : particle_set=particle_set)
1630 :
1631 24 : IF (nhomo >= 0) THEN
1632 50 : DO ispin = 1, dft_control%nspins
1633 : ! Prints the cube files of OCCUPIED ORBITALS
1634 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
1635 26 : eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
1636 50 : IF (write_cube) THEN
1637 72 : DO i = 1, nlist
1638 46 : ivector = list_index(i)
1639 46 : IF (ivector > homo) CYCLE
1640 : CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1641 46 : cell, dft_control, particle_set, pw_env)
1642 46 : WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", ivector, "_", ispin
1643 46 : mpi_io = .TRUE.
1644 : unit_nr = cp_print_key_unit_nr(logger, cube_section, '', extension=".cube", &
1645 : middle_name=TRIM(filename), file_position=my_pos_cube, &
1646 46 : log_filename=.FALSE., mpi_io=mpi_io)
1647 46 : WRITE (title, *) "WAVEFUNCTION ", ivector, " spin ", ispin, " i.e. HOMO - ", ivector - homo
1648 : CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, zeff=zcharge, &
1649 46 : stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1650 72 : CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1651 : END DO
1652 : END IF
1653 : END DO
1654 : END IF
1655 :
1656 24 : IF (nlumo /= 0) THEN
1657 6 : DO ispin = 1, dft_control%nspins
1658 : ! Prints the cube files of UNOCCUPIED ORBITALS
1659 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
1660 4 : eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
1661 6 : IF (write_cube) THEN
1662 4 : ifirst = homo + 1
1663 4 : IF (nlumo == -1) THEN
1664 0 : ilast = nmo
1665 : ELSE
1666 4 : ilast = ifirst + nlumo - 1
1667 4 : ilast = MIN(nmo, ilast)
1668 : END IF
1669 12 : DO ivector = ifirst, ilast
1670 : CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, &
1671 8 : qs_kind_set, cell, dft_control, particle_set, pw_env)
1672 8 : WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", ivector, "_", ispin
1673 8 : mpi_io = .TRUE.
1674 : unit_nr = cp_print_key_unit_nr(logger, cube_section, '', extension=".cube", &
1675 : middle_name=TRIM(filename), file_position=my_pos_cube, &
1676 8 : log_filename=.FALSE., mpi_io=mpi_io)
1677 8 : WRITE (title, *) "WAVEFUNCTION ", ivector, " spin ", ispin, " i.e. LUMO + ", ivector - ifirst
1678 : CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, zeff=zcharge, &
1679 8 : stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1680 12 : CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1681 : END DO
1682 : END IF
1683 : END DO
1684 : END IF
1685 :
1686 24 : CALL auxbas_pw_pool%give_back_pw(wf_g)
1687 24 : CALL auxbas_pw_pool%give_back_pw(wf_r)
1688 24 : IF (ASSOCIATED(list_index)) DEALLOCATE (list_index)
1689 :
1690 24 : END SUBROUTINE print_mo_cubes
1691 :
1692 : ! **************************************************************************************************
1693 :
1694 : END MODULE qs_scf_post_tb
|
