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