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