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