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