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