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