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