(git:b279b6b)
qs_scf_post_tb.F
Go to the documentation of this file.
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  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  REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
183  REAL(kind=dp), DIMENSION(:, :), POINTER :: charges
184  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
185  TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: unoccupied_evals_stm
186  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, mo_derivs
190  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  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
195  TYPE(mp_para_env_type), POINTER :: para_env
196  TYPE(particle_list_type), POINTER :: particles
197  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
198  TYPE(qs_dftb_atom_type), POINTER :: dftb_kind
199  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  CALL timeset(routinen, handle)
210 
211  logger => cp_get_default_logger()
212 
213  cpassert(ASSOCIATED(qs_env))
214  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  particle_set=particle_set, do_kpoints=do_kpoints, matrix_s_kp=matrix_s)
218  nspins = dft_control%nspins
219  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
220  ! Mulliken charges
221  ALLOCATE (charges(natom, nspins), mcharge(natom))
222  !
223  CALL mulliken_charges(matrix_p, matrix_s, para_env, charges)
224  !
225  nkind = SIZE(atomic_kind_set)
226  DO ikind = 1, nkind
227  CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
228  SELECT CASE (trim(tb_type))
229  CASE ("DFTB")
230  CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
231  CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
232  CASE ("xTB")
233  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
234  CALL get_xtb_atom_param(xtb_kind, zeff=zeff)
235  CASE DEFAULT
236  cpabort("unknown TB type")
237  END SELECT
238  DO iatom = 1, nat
239  iat = atomic_kind_set(ikind)%atom_list(iatom)
240  mcharge(iat) = zeff - sum(charges(iat, 1:nspins))
241  END DO
242  END DO
243 
244  dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
245  print_section => section_vals_get_subs_vals(dft_section, "PRINT")
246 
247  ! Mulliken
248  print_key => section_vals_get_subs_vals(print_section, "MULLIKEN")
249  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  extension=".mulliken", log_filename=.false.)
252  IF (unit_nr > 0) THEN
253  WRITE (unit=unit_nr, fmt="(/,/,T2,A)") "MULLIKEN POPULATION ANALYSIS"
254  IF (nspins == 1) THEN
255  WRITE (unit=unit_nr, fmt="(/,T2,A,T70,A)") &
256  " # Atom Element Kind Atomic population", " Net charge"
257  DO ikind = 1, nkind
258  CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
259  aname = ""
260  SELECT CASE (tb_type)
261  CASE ("DFTB")
262  CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
263  CALL get_dftb_atom_param(dftb_kind, name=aname)
264  CASE ("xTB")
265  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
266  CALL get_xtb_atom_param(xtb_kind, symbol=aname)
267  CASE DEFAULT
268  cpabort("unknown TB type")
269  END SELECT
270  ana = adjustr(trim(adjustl(aname)))
271  DO iatom = 1, nat
272  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  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  "# Total charge", sum(charges(:, 1)), sum(mcharge(:))
281  ELSE
282  WRITE (unit=unit_nr, fmt="(/,T2,A)") &
283  "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
284  DO ikind = 1, nkind
285  CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
286  aname = ""
287  SELECT CASE (tb_type)
288  CASE ("DFTB")
289  CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
290  CALL get_dftb_atom_param(dftb_kind, name=aname)
291  CASE ("xTB")
292  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
293  CALL get_xtb_atom_param(xtb_kind, symbol=aname)
294  CASE DEFAULT
295  cpabort("unknown TB type")
296  END SELECT
297  ana = adjustr(trim(adjustl(aname)))
298  DO iatom = 1, nat
299  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  iat, adjustl(ana), ikind, charges(iat, 1:2), mcharge(iat), &
303  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  "# Total charge and spin", sum(charges(:, 1)), sum(charges(:, 2)), sum(mcharge(:))
309  END IF
310  CALL m_flush(unit_nr)
311  END IF
312  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
313  END IF
314 
315  ! Compute the Lowdin charges
316  print_key => section_vals_get_subs_vals(print_section, "LOWDIN")
317  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
318  SELECT CASE (tb_type)
319  CASE ("DFTB")
320  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  log_filename=.false.)
324  print_level = 1
325  CALL section_vals_val_get(print_key, "PRINT_GOP", l_val=print_it)
326  IF (print_it) print_level = 2
327  CALL section_vals_val_get(print_key, "PRINT_ALL", l_val=print_it)
328  IF (print_it) print_level = 3
329  IF (do_kpoints) THEN
330  cpwarn("Lowdin charges not implemented for k-point calculations!")
331  ELSE
332  CALL lowdin_population_analysis(qs_env, unit_nr, print_level)
333  END IF
334  CALL cp_print_key_finished_output(unit_nr, logger, print_section, "LOWDIN")
335  CASE DEFAULT
336  cpabort("unknown TB type")
337  END SELECT
338  END IF
339 
340  ! Hirshfeld
341  print_key => section_vals_get_subs_vals(print_section, "HIRSHFELD")
342  CALL section_vals_get(print_key, explicit=explicit)
343  IF (explicit) THEN
344  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
345  cpwarn("Hirshfeld charges not available for TB methods.")
346  END IF
347  END IF
348 
349  ! MAO
350  print_key => section_vals_get_subs_vals(print_section, "MAO_ANALYSIS")
351  CALL section_vals_get(print_key, explicit=explicit)
352  IF (explicit) THEN
353  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
354  cpwarn("MAO analysis not available for TB methods.")
355  END IF
356  END IF
357 
358  ! ED
359  print_key => section_vals_get_subs_vals(print_section, "ENERGY_DECOMPOSITION_ANALYSIS")
360  CALL section_vals_get(print_key, explicit=explicit)
361  IF (explicit) THEN
362  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
363  cpwarn("ED analysis not available for TB methods.")
364  END IF
365  END IF
366 
367  ! Dipole Moments
368  print_key => section_vals_get_subs_vals(print_section, "MOMENTS")
369  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  extension=".data", middle_name="tb_dipole", log_filename=.false.)
372  moments_section => section_vals_get_subs_vals(print_section, "MOMENTS")
373  CALL tb_dipole(qs_env, moments_section, unit_nr, mcharge)
374  CALL cp_print_key_finished_output(unit_nr, logger, print_key)
375  END IF
376 
377  DEALLOCATE (charges, mcharge)
378 
379  ! MO
380  IF (.NOT. no_mos) THEN
381  print_key => section_vals_get_subs_vals(print_section, "MO")
382  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
383  CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.true.)
384  IF (.NOT. do_kpoints) THEN
385  SELECT CASE (tb_type)
386  CASE ("DFTB")
387  CASE ("xTB")
388  sprint_section => section_vals_get_subs_vals(dft_section, "PRINT%MO_MOLDEN")
389  CALL get_qs_env(qs_env, mos=mos)
390  CALL write_mos_molden(mos, qs_kind_set, particle_set, sprint_section)
391  CASE DEFAULT
392  cpabort("Unknown TB type")
393  END SELECT
394  END IF
395  END IF
396  END IF
397 
398  ! Wavefunction mixing
399  IF (.NOT. no_mos) THEN
400  wfn_mix_section => section_vals_get_subs_vals(dft_section, "PRINT%WFN_MIX")
401  CALL section_vals_get(wfn_mix_section, explicit=explicit)
402  IF (explicit .AND. .NOT. qs_env%run_rtp) CALL wfn_mix_tb(qs_env, dft_section, scf_env)
403  END IF
404 
405  IF (.NOT. no_mos) THEN
406  print_key => section_vals_get_subs_vals(print_section, "DOS")
407  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
408  IF (do_kpoints) THEN
409  CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
410  CALL calculate_dos_kp(kpoints, qs_env, dft_section)
411  ELSE
412  CALL get_qs_env(qs_env, mos=mos)
413  CALL calculate_dos(mos, dft_section)
414  END IF
415  END IF
416  END IF
417 
418  ! PDOS
419  IF (.NOT. no_mos) THEN
420  print_key => section_vals_get_subs_vals(print_section, "PDOS")
421  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
422  IF (do_kpoints) THEN
423  cpwarn("Projected density of states not implemented for k-points.")
424  ELSE
425  CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv)
426  DO ispin = 1, dft_control%nspins
427  IF (scf_env%method == ot_method_nr) THEN
428  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
429  eigenvalues=mo_eigenvalues)
430  IF (ASSOCIATED(qs_env%mo_derivs)) THEN
431  mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix
432  ELSE
433  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  co_rotate_dbcsr=mo_coeff_deriv)
438  CALL set_mo_occupation(mo_set=mos(ispin))
439  END IF
440  IF (dft_control%nspins == 2) THEN
441  CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
442  qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin)
443  ELSE
444  CALL calculate_projected_dos(mos(ispin), atomic_kind_set, &
445  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  rebuild = .false.
457  CASE ("xTB")
458  do_cube = .true.
459  rebuild = .true.
460  CASE DEFAULT
461  cpabort("unknown TB type")
462  END SELECT
463 
464  ! Energy Windows for LS code
465  print_key => section_vals_get_subs_vals(print_section, "ENERGY_WINDOWS")
466  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
467  IF (do_cube) THEN
468  IF (do_kpoints) THEN
469  cpwarn("Energy Windows not implemented for k-points.")
470  ELSE
471  IF (rebuild) THEN
472  CALL rebuild_pw_env(qs_env)
473  rebuild = .false.
474  END IF
475  CALL energy_windows(qs_env)
476  END IF
477  ELSE
478  cpwarn("Energy Windows not implemented for TB methods.")
479  END IF
480  END IF
481 
482  ! DENSITY CUBE FILE
483  print_key => section_vals_get_subs_vals(print_section, "E_DENSITY_CUBE")
484  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
485  IF (do_cube) THEN
486  IF (rebuild) THEN
487  CALL rebuild_pw_env(qs_env)
488  rebuild = .false.
489  END IF
490  CALL print_e_density(qs_env, print_key)
491  ELSE
492  cpwarn("Electronic density cube file not implemented for TB methods.")
493  END IF
494  END IF
495 
496  ! TOTAL DENSITY CUBE FILE
497  print_key => section_vals_get_subs_vals(print_section, "TOT_DENSITY_CUBE")
498  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
499  IF (do_cube) THEN
500  IF (rebuild) THEN
501  CALL rebuild_pw_env(qs_env)
502  rebuild = .false.
503  END IF
504  CALL print_density_cubes(qs_env, print_key, total_density=.true.)
505  ELSE
506  cpwarn("Total density cube file not implemented for TB methods.")
507  END IF
508  END IF
509 
510  ! V_Hartree CUBE FILE
511  print_key => section_vals_get_subs_vals(print_section, "V_HARTREE_CUBE")
512  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
513  IF (do_cube) THEN
514  IF (rebuild) THEN
515  CALL rebuild_pw_env(qs_env)
516  rebuild = .false.
517  END IF
518  CALL print_density_cubes(qs_env, print_key, v_hartree=.true.)
519  ELSE
520  cpwarn("Hartree potential cube file not implemented for TB methods.")
521  END IF
522  END IF
523 
524  ! EFIELD CUBE FILE
525  print_key => section_vals_get_subs_vals(print_section, "EFIELD_CUBE")
526  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
527  IF (do_cube) THEN
528  IF (rebuild) THEN
529  CALL rebuild_pw_env(qs_env)
530  rebuild = .false.
531  END IF
532  CALL print_density_cubes(qs_env, print_key, efield=.true.)
533  ELSE
534  cpwarn("Efield cube file not implemented for TB methods.")
535  END IF
536  END IF
537 
538  ! ELF
539  print_key => section_vals_get_subs_vals(print_section, "ELF_CUBE")
540  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
541  IF (do_cube) THEN
542  IF (rebuild) THEN
543  CALL rebuild_pw_env(qs_env)
544  rebuild = .false.
545  END IF
546  CALL print_elf(qs_env, print_key)
547  ELSE
548  cpwarn("ELF not implemented for TB methods.")
549  END IF
550  END IF
551 
552  ! MO CUBES
553  IF (.NOT. no_mos) THEN
554  print_key => section_vals_get_subs_vals(print_section, "MO_CUBES")
555  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
556  IF (do_cube) THEN
557  IF (rebuild) THEN
558  CALL rebuild_pw_env(qs_env)
559  rebuild = .false.
560  END IF
561  CALL print_mo_cubes(qs_env, print_key)
562  ELSE
563  cpwarn("Printing of MO cube files not implemented for TB methods.")
564  END IF
565  END IF
566  END IF
567 
568  ! STM
569  IF (.NOT. no_mos) THEN
570  print_key => section_vals_get_subs_vals(print_section, "STM")
571  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
572  IF (do_cube) THEN
573  IF (rebuild) THEN
574  CALL rebuild_pw_env(qs_env)
575  rebuild = .false.
576  END IF
577  IF (do_kpoints) THEN
578  cpwarn("STM not implemented for k-point calculations!")
579  ELSE
580  nlumo_stm = section_get_ival(print_key, "NLUMO")
581  cpassert(.NOT. dft_control%restricted)
582  CALL get_qs_env(qs_env, mos=mos, mo_derivs=mo_derivs, &
583  scf_control=scf_control, matrix_ks=ks_rmpv)
584  CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs)
585  DO ispin = 1, dft_control%nspins
586  CALL get_mo_set(mo_set=mos(ispin), eigenvalues=mo_eigenvalues, homo=homo)
587  homo_lumo(ispin, 1) = mo_eigenvalues(homo)
588  END DO
589  has_homo = .true.
590  NULLIFY (unoccupied_orbs_stm, unoccupied_evals_stm)
591  IF (nlumo_stm > 0) THEN
592  ALLOCATE (unoccupied_orbs_stm(dft_control%nspins))
593  ALLOCATE (unoccupied_evals_stm(dft_control%nspins))
594  CALL make_lumo_tb(qs_env, scf_env, unoccupied_orbs_stm, unoccupied_evals_stm, &
595  nlumo_stm, nlumos)
596  END IF
597 
598  CALL get_qs_env(qs_env, subsys=subsys)
599  CALL qs_subsys_get(subsys, particles=particles)
600  CALL th_stm_image(qs_env, print_key, particles, unoccupied_orbs_stm, &
601  unoccupied_evals_stm)
602 
603  IF (nlumo_stm > 0) THEN
604  DO ispin = 1, dft_control%nspins
605  DEALLOCATE (unoccupied_evals_stm(ispin)%array)
606  END DO
607  DEALLOCATE (unoccupied_evals_stm)
608  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  CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks)
617  CALL section_vals_val_get(print_section, "AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
618  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  extension=".Log")
622  CALL section_vals_val_get(print_section, "AO_MATRICES%NDIGITS", i_val=after)
623  after = min(max(after, 1), 16)
624  DO ispin = 1, dft_control%nspins
625  DO img = 1, SIZE(matrix_p, 2)
626  CALL cp_dbcsr_write_sparse_matrix(matrix_p(ispin, img)%matrix, 4, after, qs_env, &
627  para_env, output_unit=iw, omit_headers=omit_headers)
628  END DO
629  END DO
630  CALL cp_print_key_finished_output(iw, logger, print_section, "AO_MATRICES/DENSITY")
631  END IF
632 
633  ! The xTB matrix itself
634  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  extension=".Log")
638  CALL section_vals_val_get(print_section, "AO_MATRICES%NDIGITS", i_val=after)
639  after = min(max(after, 1), 16)
640  DO ispin = 1, dft_control%nspins
641  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  output_unit=iw, omit_headers=omit_headers)
644  END DO
645  END DO
646  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  print_key => section_vals_get_subs_vals(print_section, "V_XC_CUBE")
653  CALL section_vals_get(print_key, explicit=explicit)
654  IF (explicit) THEN
655  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
656  cpwarn("XC potential cube file not available for TB methods.")
657  END IF
658  END IF
659 
660  ! Electric field gradients
661  print_key => section_vals_get_subs_vals(print_section, "ELECTRIC_FIELD_GRADIENT")
662  CALL section_vals_get(print_key, explicit=explicit)
663  IF (explicit) THEN
664  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
665  cpwarn("Electric field gradient not implemented for TB methods.")
666  END IF
667  END IF
668 
669  ! KINETIC ENERGY
670  print_key => section_vals_get_subs_vals(print_section, "KINETIC_ENERGY")
671  CALL section_vals_get(print_key, explicit=explicit)
672  IF (explicit) THEN
673  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
674  cpwarn("Kinetic energy not available for TB methods.")
675  END IF
676  END IF
677 
678  ! Xray diffraction spectrum
679  print_key => section_vals_get_subs_vals(print_section, "XRAY_DIFFRACTION_SPECTRUM")
680  CALL section_vals_get(print_key, explicit=explicit)
681  IF (explicit) THEN
682  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
683  cpwarn("Xray diffraction spectrum not implemented for TB methods.")
684  END IF
685  END IF
686 
687  ! EPR Hyperfine Coupling
688  print_key => section_vals_get_subs_vals(print_section, "HYPERFINE_COUPLING_TENSOR")
689  CALL section_vals_get(print_key, explicit=explicit)
690  IF (explicit) THEN
691  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
692  cpwarn("Hyperfine Coupling not implemented for TB methods.")
693  END IF
694  END IF
695 
696  ! PLUS_U
697  print_key => section_vals_get_subs_vals(print_section, "PLUS_U")
698  CALL section_vals_get(print_key, explicit=explicit)
699  IF (explicit) THEN
700  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
701  cpwarn("DFT+U method not implemented for TB methods.")
702  END IF
703  END IF
704 
705  CALL write_ks_matrix_csr(qs_env, qs_env%input)
706  CALL write_s_matrix_csr(qs_env, qs_env%input)
707 
708  CALL timestop(handle)
709 
710  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
734  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
735  TYPE(cell_type), POINTER :: cell
736  TYPE(cp_result_type), POINTER :: results
737  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
738 
739  NULLIFY (atomic_kind_set, cell, results)
740  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, &
741  particle_set=particle_set, cell=cell, results=results)
742 
743  ! Reference point
744  reference = section_get_ival(input, keyword_name="REFERENCE")
745  NULLIFY (ref_point)
746  description = '[DIPOLE]'
747  CALL section_vals_val_get(input, "REF_POINT", r_vals=ref_point)
748  CALL section_vals_val_get(input, "PERIODIC", l_val=do_berry)
749 
750  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  dipole_deriv = 0.0_dp
754  dipole = 0.0_dp
755  IF (do_berry) THEN
756  dipole_type = "periodic (Berry phase)"
757  rcc = pbc(rcc, cell)
758  charge_tot = 0._dp
759  charge_tot = sum(charges)
760  ria = twopi*matmul(cell%h_inv, rcc)
761  zphase = cmplx(cos(ria), sin(ria), dp)**charge_tot
762 
763  dria = twopi*matmul(cell%h_inv, drcc)
764  dzphase = charge_tot*cmplx(-sin(ria), cos(ria), dp)**(charge_tot - 1.0_dp)*dria
765 
766  ggamma = cmplx(1.0_dp, 0.0_dp, kind=dp)
767  dggamma = cmplx(0.0_dp, 0.0_dp, kind=dp)
768  DO ikind = 1, SIZE(atomic_kind_set)
769  CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
770  DO i = 1, nat
771  iat = atomic_kind_set(ikind)%atom_list(i)
772  ria = particle_set(iat)%r(:)
773  ria = pbc(ria, cell)
774  via = particle_set(iat)%v(:)
775  q = charges(iat)
776  DO j = 1, 3
777  gvec = twopi*cell%h_inv(j, :)
778  theta = sum(ria(:)*gvec(:))
779  dtheta = sum(via(:)*gvec(:))
780  zeta = cmplx(cos(theta), sin(theta), kind=dp)**(-q)
781  dzeta = -q*cmplx(-sin(theta), cos(theta), kind=dp)**(-q - 1.0_dp)*dtheta
782  dggamma(j) = dggamma(j)*zeta + ggamma(j)*dzeta
783  ggamma(j) = ggamma(j)*zeta
784  END DO
785  END DO
786  END DO
787  dggamma = dggamma*zphase + ggamma*dzphase
788  ggamma = ggamma*zphase
789  IF (all(real(ggamma, kind=dp) /= 0.0_dp)) THEN
790  tmp = aimag(ggamma)/real(ggamma, kind=dp)
791  ci = -atan(tmp)
792  dci = -(1.0_dp/(1.0_dp + tmp**2))* &
793  (aimag(dggamma)*real(ggamma, kind=dp) - aimag(ggamma)*real(dggamma, kind=dp))/(real(ggamma, kind=dp))**2
794  dipole = matmul(cell%hmat, ci)/twopi
795  dipole_deriv = matmul(cell%hmat, dci)/twopi
796  END IF
797  ELSE
798  dipole_type = "non-periodic"
799  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  ria = particle_set(i)%r(:)
802  q = charges(i)
803  dipole = dipole + q*(ria - rcc)
804  dipole_deriv(:) = dipole_deriv(:) + q*(particle_set(i)%v(:) - drcc)
805  END DO
806  END IF
807  CALL cp_results_erase(results=results, description=description)
808  CALL put_results(results=results, description=description, &
809  values=dipole(1:3))
810  IF (unit_nr > 0) THEN
811  WRITE (unit_nr, '(/,T2,A,T31,A50)') &
812  'TB_DIPOLE| Dipole type', adjustr(trim(dipole_type))
813  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  'TB_DIPOLE| Moment [a.u.]', dipole(1:3)
816  WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
817  'TB_DIPOLE| Moment [Debye]', dipole(1:3)*debye
818  WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
819  'TB_DIPOLE| Derivative [a.u.]', dipole_deriv(1:3)
820  END IF
821 
822  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  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  REAL(dp), DIMENSION(:), POINTER :: mo_eigenvalues
841  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  TYPE(cp_fm_type), DIMENSION(:), POINTER :: lumos
845  TYPE(cp_fm_type), POINTER :: mo_coeff
846  TYPE(cp_logger_type), POINTER :: logger
847  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
848  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
849  TYPE(mp_para_env_type), POINTER :: para_env
850  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
851  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
852  TYPE(section_vals_type), POINTER :: wfn_mix_section
853 
854  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  qs_kind_set=qs_kind_set, mos=mos, para_env=para_env)
858 
859  wfn_mix_section => section_vals_get_subs_vals(dft_section, "PRINT%WFN_MIX")
860 
861  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  template_fmstruct=mo_coeff%matrix_struct)
865  CALL cp_fm_create(s_tmp, matrix_struct=ao_ao_fmstruct)
866  CALL cp_fm_create(ks_tmp, matrix_struct=ao_ao_fmstruct)
867  CALL cp_fm_create(mo_tmp, matrix_struct=ao_ao_fmstruct)
868  CALL cp_fm_create(work, matrix_struct=ao_ao_fmstruct)
869  ALLOCATE (lumos(SIZE(mos)))
870 
871  CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, s_tmp)
872  CALL cp_fm_cholesky_decompose(s_tmp)
873 
874  DO ispin = 1, SIZE(mos)
875  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  template_fmstruct=mo_coeff%matrix_struct)
878 
879  CALL cp_fm_create(lumos(ispin), matrix_struct=ao_lumo_struct)
880  CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, ks_tmp)
881  CALL cp_fm_cholesky_reduce(ks_tmp, s_tmp)
882  CALL choose_eigv_solver(ks_tmp, work, mo_eigenvalues)
883  CALL cp_fm_cholesky_restore(work, nao, s_tmp, mo_tmp, "SOLVE")
884  CALL cp_fm_to_fm_submat(mo_tmp, mo_coeff, nao, nmo, 1, 1, 1, 1)
885  CALL cp_fm_to_fm_submat(mo_tmp, lumos(ispin), nao, nao - nmo, 1, nmo + 1, 1, 1)
886 
887  CALL cp_fm_struct_release(ao_lumo_struct)
888  END DO
889 
890  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  unoccupied_orbs=lumos, scf_env=scf_env, matrix_s=matrix_s)
893 
894  CALL cp_fm_release(lumos)
895  CALL cp_fm_release(s_tmp)
896  CALL cp_fm_release(mo_tmp)
897  CALL cp_fm_release(ks_tmp)
898  CALL cp_fm_release(work)
899  CALL cp_fm_struct_release(ao_ao_fmstruct)
900 
901  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, matrix_s
927  TYPE(dft_control_type), POINTER :: dft_control
928  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  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  blacs_env=blacs_env)
942 
943  logger => cp_get_default_logger()
944  iounit = cp_logger_get_default_io_unit(logger)
945 
946  DO ispin = 1, dft_control%nspins
947  NULLIFY (unoccupied_evals(ispin)%array)
948  ! Always write eigenvalues
949  IF (iounit > 0) WRITE (iounit, *) " "
950  IF (iounit > 0) WRITE (iounit, *) " Lowest Eigenvalues of the unoccupied subspace spin ", ispin
951  IF (iounit > 0) WRITE (iounit, fmt='(1X,A)') "-----------------------------------------------------"
952  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=homo, nao=nao, nmo=nmo)
953  CALL cp_fm_get_info(mo_coeff, nrow_global=n)
954  nlumos = max(1, min(nlumo, nao - nmo))
955  IF (nlumo == -1) nlumos = nao - nmo
956  ALLOCATE (unoccupied_evals(ispin)%array(nlumos))
957  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
958  nrow_global=n, ncol_global=nlumos)
959  CALL cp_fm_create(unoccupied_orbs(ispin), fm_struct_tmp, name="lumos")
960  CALL cp_fm_struct_release(fm_struct_tmp)
961  CALL cp_fm_init_random(unoccupied_orbs(ispin), nlumos)
962 
963  ! the full_all preconditioner makes not much sense for lumos search
964  NULLIFY (local_preconditioner)
965  IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
966  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  IF (local_preconditioner%in_use == ot_precond_full_all) THEN
969  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  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  ionode=iounit > 0)
984 
985  END DO
986 
987  END SUBROUTINE make_lumo_tb
988 
989 ! **************************************************************************************************
990 !> \brief ...
991 !> \param qs_env ...
992 ! **************************************************************************************************
993  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  CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=new_pw_env)
1006  IF (.NOT. ASSOCIATED(new_pw_env)) THEN
1007  CALL pw_env_create(new_pw_env)
1008  CALL set_ks_env(ks_env, pw_env=new_pw_env)
1009  CALL pw_env_release(new_pw_env)
1010  END IF
1011  CALL get_qs_env(qs_env, pw_env=new_pw_env, dft_control=dft_control, cell=cell)
1012 
1013  new_pw_env%cell_hmat = cell%hmat
1014  CALL pw_env_rebuild(new_pw_env, qs_env=qs_env)
1015 
1016  NULLIFY (task_list)
1017  CALL get_ks_env(ks_env, task_list=task_list)
1018  IF (.NOT. ASSOCIATED(task_list)) THEN
1019  CALL allocate_task_list(task_list)
1020  CALL set_ks_env(ks_env, task_list=task_list)
1021  END IF
1022  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  skip_load_balance_distributed=skip_load_balance_distributed)
1026  CALL get_qs_env(qs_env, rho=rho)
1027  CALL qs_rho_rebuild(rho, qs_env=qs_env, rebuild_ao=.false., rebuild_grids=.true.)
1028 
1029  END SUBROUTINE rebuild_pw_env
1030 
1031 ! **************************************************************************************************
1032 !> \brief ...
1033 !> \param qs_env ...
1034 !> \param cube_section ...
1035 ! **************************************************************************************************
1036  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  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
1045  TYPE(cp_logger_type), POINTER :: logger
1046  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
1047  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  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1051  TYPE(pw_env_type), POINTER :: pw_env
1052  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1053  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1054  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  CALL get_qs_env(qs_env, dft_control=dft_control)
1060 
1061  append_cube = section_get_lval(cube_section, "APPEND")
1062  my_pos_cube = "REWIND"
1063  IF (append_cube) my_pos_cube = "APPEND"
1064 
1065  logger => cp_get_default_logger()
1066  iounit = cp_logger_get_default_io_unit(logger)
1067 
1068  ! we need to construct the density on a realspace grid
1069  CALL get_qs_env(qs_env, ks_env=ks_env, rho=rho)
1070  NULLIFY (rho_r, rho_g, tot_rho_r)
1071  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, &
1072  rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
1073  DO ispin = 1, dft_control%nspins
1074  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  ks_env=ks_env)
1080  END DO
1081  CALL qs_rho_set(rho, rho_r_valid=.true., rho_g_valid=.true.)
1082 
1083  CALL get_qs_env(qs_env, subsys=subsys)
1084  CALL qs_subsys_get(subsys, particles=particles)
1085 
1086  IF (dft_control%nspins > 1) THEN
1087  IF (iounit > 0) THEN
1088  WRITE (unit=iounit, fmt="(/,T2,A,T51,2F15.6)") &
1089  "Integrated alpha and beta electronic density:", tot_rho_r(1:2)
1090  END IF
1091  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1092  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  CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
1096  CALL pw_copy(rho_r(1), rho_elec_rspace)
1097  CALL pw_axpy(rho_r(2), rho_elec_rspace)
1098  filename = "ELECTRON_DENSITY"
1099  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  fout=mpi_filename)
1104  IF (iounit > 0) THEN
1105  IF (.NOT. mpi_io) THEN
1106  INQUIRE (unit=unit_nr, name=filename)
1107  ELSE
1108  filename = mpi_filename
1109  END IF
1110  WRITE (unit=iounit, fmt="(T2,A,/,T2,A79)") &
1111  "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  mpi_io=mpi_io)
1116  CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1117  CALL pw_copy(rho_r(1), rho_elec_rspace)
1118  CALL pw_axpy(rho_r(2), rho_elec_rspace, alpha=-1.0_dp)
1119  filename = "SPIN_DENSITY"
1120  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  fout=mpi_filename)
1125  IF (iounit > 0) THEN
1126  IF (.NOT. mpi_io) THEN
1127  INQUIRE (unit=unit_nr, name=filename)
1128  ELSE
1129  filename = mpi_filename
1130  END IF
1131  WRITE (unit=iounit, fmt="(T2,A,/,T2,A79)") &
1132  "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  stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1137  CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1138  CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
1139  END block
1140  ELSE
1141  IF (iounit > 0) THEN
1142  WRITE (unit=iounit, fmt="(/,T2,A,T66,F15.6)") &
1143  "Integrated electronic density:", tot_rho_r(1)
1144  END IF
1145  filename = "ELECTRON_DENSITY"
1146  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  fout=mpi_filename)
1151  IF (iounit > 0) THEN
1152  IF (.NOT. mpi_io) THEN
1153  INQUIRE (unit=unit_nr, name=filename)
1154  ELSE
1155  filename = mpi_filename
1156  END IF
1157  WRITE (unit=iounit, fmt="(T2,A,/,T2,A79)") &
1158  "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  stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1163  CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1164  END IF ! nspins
1165 
1166  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
1189  TYPE(cell_type), POINTER :: cell
1190  TYPE(cp_logger_type), POINTER :: logger
1191  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
1192  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  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  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  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  CALL get_qs_env(qs_env, cell=cell, dft_control=dft_control)
1208 
1209  append_cube = section_get_lval(cube_section, "APPEND")
1210  my_pos_cube = "REWIND"
1211  IF (append_cube) my_pos_cube = "APPEND"
1212 
1213  IF (PRESENT(total_density)) THEN
1214  my_total_density = total_density
1215  ELSE
1216  my_total_density = .false.
1217  END IF
1218  IF (PRESENT(v_hartree)) THEN
1219  my_v_hartree = v_hartree
1220  ELSE
1221  my_v_hartree = .false.
1222  END IF
1223  IF (PRESENT(efield)) THEN
1224  my_efield = efield
1225  ELSE
1226  my_efield = .false.
1227  END IF
1228 
1229  logger => cp_get_default_logger()
1230  iounit = cp_logger_get_default_io_unit(logger)
1231 
1232  ! we need to construct the density on a realspace grid
1233  CALL get_qs_env(qs_env, ks_env=ks_env, rho=rho)
1234  NULLIFY (rho_r, rho_g, tot_rho_r)
1235  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, &
1236  rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
1237  DO ispin = 1, dft_control%nspins
1238  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  ks_env=ks_env)
1244  END DO
1245  CALL qs_rho_set(rho, rho_r_valid=.true., rho_g_valid=.true.)
1246 
1247  CALL get_qs_env(qs_env, subsys=subsys)
1248  CALL qs_subsys_get(subsys, particles=particles)
1249 
1250  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1251  CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1252  CALL auxbas_pw_pool%create_pw(pw=rho_core)
1253  CALL calculate_rho_core(rho_core, total_rho_core_rspace, qs_env)
1254 
1255  IF (iounit > 0) THEN
1256  WRITE (unit=iounit, fmt="(/,T2,A,T66,F15.6)") &
1257  "Integrated electronic density:", sum(tot_rho_r(:))
1258  WRITE (unit=iounit, fmt="(T2,A,T66,F15.6)") &
1259  "Integrated core density:", total_rho_core_rspace
1260  END IF
1261 
1262  CALL auxbas_pw_pool%create_pw(pw=rho_tot_rspace)
1263  CALL pw_transfer(rho_core, rho_tot_rspace)
1264  DO ispin = 1, dft_control%nspins
1265  CALL pw_axpy(rho_r(ispin), rho_tot_rspace)
1266  END DO
1267 
1268  IF (my_total_density) THEN
1269  filename = "TOTAL_DENSITY"
1270  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  log_filename=.false., mpi_io=mpi_io, fout=mpi_filename)
1274  IF (iounit > 0) THEN
1275  IF (.NOT. mpi_io) THEN
1276  INQUIRE (unit=unit_nr, name=filename)
1277  ELSE
1278  filename = mpi_filename
1279  END IF
1280  WRITE (unit=iounit, fmt="(T2,A,/,T2,A79)") &
1281  "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  stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1286  CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1287  END IF
1288  IF (my_v_hartree .OR. my_efield) THEN
1289  block
1290  TYPE(pw_c1d_gs_type) :: rho_tot_gspace
1291  CALL auxbas_pw_pool%create_pw(pw=rho_tot_gspace)
1292  CALL pw_transfer(rho_tot_rspace, rho_tot_gspace)
1293  poisson_params%solver = pw_poisson_analytic
1294  poisson_params%periodic = cell%perd
1295  poisson_params%ewald_type = do_ewald_none
1296  block
1297  TYPE(greens_fn_type) :: green_fft
1298  TYPE(pw_grid_type), POINTER :: pwdummy
1299  NULLIFY (pwdummy)
1300  CALL pw_green_create(green_fft, poisson_params, cell%hmat, auxbas_pw_pool, pwdummy, pwdummy)
1301  rho_tot_gspace%array(:) = rho_tot_gspace%array(:)*green_fft%influence_fn%array(:)
1302  CALL pw_green_release(green_fft, auxbas_pw_pool)
1303  END block
1304  IF (my_v_hartree) THEN
1305  block
1306  TYPE(pw_r3d_rs_type) :: vhartree
1307  CALL auxbas_pw_pool%create_pw(pw=vhartree)
1308  CALL pw_transfer(rho_tot_gspace, vhartree)
1309  filename = "V_HARTREE"
1310  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  log_filename=.false., mpi_io=mpi_io, fout=mpi_filename)
1314  IF (iounit > 0) THEN
1315  IF (.NOT. mpi_io) THEN
1316  INQUIRE (unit=unit_nr, name=filename)
1317  ELSE
1318  filename = mpi_filename
1319  END IF
1320  WRITE (unit=iounit, fmt="(T2,A,/,T2,A79)") &
1321  "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  stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1326  CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1327  CALL auxbas_pw_pool%give_back_pw(vhartree)
1328  END block
1329  END IF
1330  IF (my_efield) THEN
1331  block
1332  TYPE(pw_c1d_gs_type) :: vhartree
1333  CALL auxbas_pw_pool%create_pw(pw=vhartree)
1334  udvol = 1.0_dp/rho_tot_rspace%pw_grid%dvol
1335  DO id = 1, 3
1336  CALL pw_transfer(rho_tot_gspace, vhartree)
1337  nd = 0
1338  nd(id) = 1
1339  CALL pw_derive(vhartree, nd)
1340  CALL pw_transfer(vhartree, rho_tot_rspace)
1341  CALL pw_scale(rho_tot_rspace, udvol)
1342 
1343  filename = "EFIELD_"//cdir(id)
1344  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  log_filename=.false., mpi_io=mpi_io, fout=mpi_filename)
1348  IF (iounit > 0) THEN
1349  IF (.NOT. mpi_io) THEN
1350  INQUIRE (unit=unit_nr, name=filename)
1351  ELSE
1352  filename = mpi_filename
1353  END IF
1354  WRITE (unit=iounit, fmt="(T2,A,/,T2,A79)") &
1355  "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  stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1360  CALL cp_print_key_finished_output(unit_nr, logger, cube_section, '', mpi_io=mpi_io)
1361  END DO
1362  CALL auxbas_pw_pool%give_back_pw(vhartree)
1363  END block
1364  END IF
1365  CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1366  END block
1367  END IF
1368 
1369  CALL auxbas_pw_pool%give_back_pw(rho_tot_rspace)
1370  CALL auxbas_pw_pool%give_back_pw(rho_core)
1371 
1372  END SUBROUTINE print_density_cubes
1373 
1374 ! **************************************************************************************************
1375 !> \brief ...
1376 !> \param qs_env ...
1377 !> \param elf_section ...
1378 ! **************************************************************************************************
1379  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  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
1390  TYPE(cp_logger_type), POINTER :: logger
1391  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
1392  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  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1396  TYPE(pw_env_type), POINTER :: pw_env
1397  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1398  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1399  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: elf_r
1400  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  logger => cp_get_default_logger()
1406  iounit = cp_logger_get_default_io_unit(logger)
1407 
1408  ! we need to construct the density on a realspace grid
1409  CALL get_qs_env(qs_env, dft_control=dft_control, ks_env=ks_env, rho=rho)
1410  NULLIFY (rho_r, rho_g, tot_rho_r)
1411  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, &
1412  rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
1413  DO ispin = 1, dft_control%nspins
1414  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  ks_env=ks_env)
1420  END DO
1421  CALL qs_rho_set(rho, rho_r_valid=.true., rho_g_valid=.true.)
1422 
1423  CALL get_qs_env(qs_env, subsys=subsys)
1424  CALL qs_subsys_get(subsys, particles=particles)
1425 
1426  ALLOCATE (elf_r(dft_control%nspins))
1427  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1428  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1429  DO ispin = 1, dft_control%nspins
1430  CALL auxbas_pw_pool%create_pw(elf_r(ispin))
1431  CALL pw_zero(elf_r(ispin))
1432  END DO
1433 
1434  IF (iounit > 0) THEN
1435  WRITE (unit=iounit, fmt="(/,T2,A)") &
1436  "ELF is computed on the real space grid -----"
1437  END IF
1438  rho_cutoff = section_get_rval(elf_section, "density_cutoff")
1439  CALL qs_elf_calc(qs_env, elf_r, rho_cutoff)
1440 
1441  ! write ELF into cube file
1442  append_cube = section_get_lval(elf_section, "APPEND")
1443  my_pos_cube = "REWIND"
1444  IF (append_cube) my_pos_cube = "APPEND"
1445  DO ispin = 1, dft_control%nspins
1446  WRITE (filename, '(a5,I1.1)') "ELF_S", ispin
1447  WRITE (title, *) "ELF spin ", ispin
1448  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  log_filename=.false., mpi_io=mpi_io, fout=mpi_filename)
1452  IF (iounit > 0) THEN
1453  IF (.NOT. mpi_io) THEN
1454  INQUIRE (unit=unit_nr, name=filename)
1455  ELSE
1456  filename = mpi_filename
1457  END IF
1458  WRITE (unit=iounit, fmt="(T2,A,/,T2,A79)") &
1459  "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  stride=section_get_ivals(elf_section, "STRIDE"), mpi_io=mpi_io)
1464  CALL cp_print_key_finished_output(unit_nr, logger, elf_section, '', mpi_io=mpi_io)
1465 
1466  CALL auxbas_pw_pool%give_back_pw(elf_r(ispin))
1467  END DO
1468 
1469  DEALLOCATE (elf_r)
1470 
1471  END SUBROUTINE print_elf
1472 ! **************************************************************************************************
1473 !> \brief ...
1474 !> \param qs_env ...
1475 !> \param cube_section ...
1476 ! **************************************************************************************************
1477  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  INTEGER, DIMENSION(:), POINTER :: list, list_index
1487  LOGICAL :: append_cube, mpi_io, write_cube
1488  REAL(kind=dp) :: homo_lumo(2, 2)
1489  REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
1490  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, mo_derivs
1495  TYPE(dft_control_type), POINTER :: dft_control
1496  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1497  TYPE(particle_list_type), POINTER :: particles
1498  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1499  TYPE(pw_c1d_gs_type) :: wf_g
1500  TYPE(pw_env_type), POINTER :: pw_env
1501  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  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  logger => cp_get_default_logger()
1509  iounit = cp_logger_get_default_io_unit(logger)
1510 
1511  CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv, scf_control=scf_control)
1512  CALL get_qs_env(qs_env, dft_control=dft_control, mo_derivs=mo_derivs)
1513  CALL make_mo_eig(mos, dft_control%nspins, ks_rmpv, scf_control, mo_derivs)
1514  NULLIFY (mo_eigenvalues)
1515  homo = 0
1516  DO ispin = 1, dft_control%nspins
1517  CALL get_mo_set(mo_set=mos(ispin), eigenvalues=mo_eigenvalues, homo=shomo)
1518  homo_lumo(ispin, 1) = mo_eigenvalues(shomo)
1519  homo = max(homo, shomo)
1520  END DO
1521  write_cube = section_get_lval(cube_section, "WRITE_CUBE")
1522  nlumo = section_get_ival(cube_section, "NLUMO")
1523  nhomo = section_get_ival(cube_section, "NHOMO")
1524  NULLIFY (list_index)
1525  CALL section_vals_val_get(cube_section, "HOMO_LIST", n_rep_val=n_rep)
1526  IF (n_rep > 0) THEN
1527  nlist = 0
1528  DO ir = 1, n_rep
1529  NULLIFY (list)
1530  CALL section_vals_val_get(cube_section, "HOMO_LIST", i_rep_val=ir, i_vals=list)
1531  IF (ASSOCIATED(list)) THEN
1532  CALL reallocate(list_index, 1, nlist + SIZE(list))
1533  DO i = 1, SIZE(list)
1534  list_index(i + nlist) = list(i)
1535  END DO
1536  nlist = nlist + SIZE(list)
1537  END IF
1538  END DO
1539  nhomo = maxval(list_index)
1540  ELSE
1541  IF (nhomo == -1) nhomo = homo
1542  nlist = homo - max(1, homo - nhomo + 1) + 1
1543  ALLOCATE (list_index(nlist))
1544  DO i = 1, nlist
1545  list_index(i) = max(1, homo - nhomo + 1) + i - 1
1546  END DO
1547  END IF
1548 
1549  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1550  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1551  CALL auxbas_pw_pool%create_pw(wf_r)
1552  CALL auxbas_pw_pool%create_pw(wf_g)
1553 
1554  CALL get_qs_env(qs_env, subsys=subsys)
1555  CALL qs_subsys_get(subsys, particles=particles)
1556 
1557  append_cube = section_get_lval(cube_section, "APPEND")
1558  my_pos_cube = "REWIND"
1559  IF (append_cube) THEN
1560  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  particle_set=particle_set)
1568 
1569  IF (nhomo >= 0) THEN
1570  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  eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
1574  IF (write_cube) THEN
1575  DO i = 1, nlist
1576  ivector = list_index(i)
1577  IF (ivector > homo) cycle
1578  CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1579  cell, dft_control, particle_set, pw_env)
1580  WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", ivector, "_", ispin
1581  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  log_filename=.false., mpi_io=mpi_io)
1585  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  stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1588  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  IF (nlumo /= 0) THEN
1595  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  eigenvalues=mo_eigenvalues, homo=homo, nmo=nmo)
1599  IF (write_cube) THEN
1600  ifirst = homo + 1
1601  IF (nlumo == -1) THEN
1602  ilast = nmo
1603  ELSE
1604  ilast = ifirst + nlumo - 1
1605  ilast = min(nmo, ilast)
1606  END IF
1607  DO ivector = ifirst, ilast
1608  CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, &
1609  qs_kind_set, cell, dft_control, particle_set, pw_env)
1610  WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", ivector, "_", ispin
1611  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  log_filename=.false., mpi_io=mpi_io)
1615  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  stride=section_get_ivals(cube_section, "STRIDE"), mpi_io=mpi_io)
1618  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  CALL auxbas_pw_pool%give_back_pw(wf_g)
1625  CALL auxbas_pw_pool%give_back_pw(wf_r)
1626  IF (ASSOCIATED(list_index)) DEALLOCATE (list_index)
1627 
1628  END SUBROUTINE print_mo_cubes
1629 
1630 ! **************************************************************************************************
1631 
1632 END MODULE qs_scf_post_tb
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition: cp_array_utils.F:18
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_dbcsr_output
DBCSR output in CP2K.
Definition: cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
Definition: cp_dbcsr_output.F:163
cp_fm_cholesky
various cholesky decomposition related routines
Definition: cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition: cp_fm_cholesky.F:50
cp_fm_cholesky::cp_fm_cholesky_reduce
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
Definition: cp_fm_cholesky.F:192
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa)
...
Definition: cp_fm_cholesky.F:254
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition: cp_fm_diag.F:208
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition: cp_fm_types.F:1473
cp_fm_types::cp_fm_init_random
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
Definition: cp_fm_types.F:452
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition: cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition: cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
Definition: cp_realspace_grid_cube.F:45
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition: cp_result_methods.F:17
cp_result_methods::cp_results_erase
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
Definition: cp_result_methods.F:362
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition: cp_result_types.F:18
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition: input_constants.F:513
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_get_rval
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
Definition: input_section_types.F:941
input_section_types::section_get_ivals
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
Definition: input_section_types.F:989
input_section_types::section_get_ival
integer function, public section_get_ival(section_vals, keyword_name)
...
Definition: input_section_types.F:973
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition: input_section_types.F:1005
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
molden_utils
Functions handling the MOLDEN format. Split from mode_selective.
Definition: molden_utils.F:12
molden_utils::write_mos_molden
subroutine, public write_mos_molden(mos, qs_kind_set, particle_set, print_section)
Write out the MOs in molden format for visualisation.
Definition: molden_utils.F:65
moments_utils
Calculates the moment integrals <a|r^m|b>
Definition: moments_utils.F:15
moments_utils::get_reference_point
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition: moments_utils.F:61
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition: mulliken.F:13
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::debye
real(kind=dp), parameter, public debye
Definition: physcon.F:201
population_analyses
Provide various population analyses and print the requested output information.
Definition: population_analyses.F:16
population_analyses::lowdin_population_analysis
subroutine, public lowdin_population_analysis(qs_env, output_unit, print_level)
Perform a Lowdin population analysis based on a symmetric orthogonalisation of the density matrix usi...
Definition: population_analyses.F:83
preconditioner_types
types of preconditioners
Definition: preconditioner_types.F:14
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition: preconditioner.F:15
pw_env_methods
methods of pw_env that have dependence on qs_env
Definition: pw_env_methods.F:17
pw_env_methods::pw_env_rebuild
subroutine, public pw_env_rebuild(pw_env, qs_env, external_para_env)
rebuilds the pw_env data (necessary if cell or cutoffs change)
Definition: pw_env_methods.F:160
pw_env_methods::pw_env_create
subroutine, public pw_env_create(pw_env)
creates a pw_env, if qs_env is given calls pw_env_rebuild
Definition: pw_env_methods.F:129
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_release
subroutine, public pw_env_release(pw_env, para_env)
releases the given pw_env (see doc/ReferenceCounting.html)
Definition: pw_env_types.F:176
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_grid_types
Definition: pw_grid_types.F:13
pw_methods
Definition: pw_methods.F:25
pw_methods::pw_derive
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition: pw_methods.F:10106
pw_poisson_types
functions related to the poisson solver on regular grids
Definition: pw_poisson_types.F:22
pw_poisson_types::pw_green_create
subroutine, public pw_green_create(green, poisson_params, cell_hmat, pw_pool, mt_super_ref_pw_grid, dct_pw_grid)
Allocates and sets up the green functions for the fft based poisson solvers.
Definition: pw_poisson_types.F:182
pw_poisson_types::pw_green_release
subroutine, public pw_green_release(gftype, pw_pool)
destroys the type (deallocates data)
Definition: pw_poisson_types.F:418
pw_poisson_types::do_ewald_none
integer, parameter, public do_ewald_none
Definition: pw_poisson_types.F:75
pw_poisson_types::pw_poisson_analytic
integer, parameter, public pw_poisson_analytic
Definition: pw_poisson_types.F:66
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition: qs_collocate_density.F:1710
qs_collocate_density::calculate_wavefunction
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type, external_vector)
maps a given wavefunction on the grid
Definition: qs_collocate_density.F:2860
qs_dftb_types
Definition of the DFTB parameter types.
Definition: qs_dftb_types.F:12
qs_dftb_utils
Working with the DFTB parameter types.
Definition: qs_dftb_utils.F:12
qs_dftb_utils::get_dftb_atom_param
subroutine, public get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq)
...
Definition: qs_dftb_utils.F:123
qs_dos
Calculation and writing of density of states.
Definition: qs_dos.F:14
qs_dos::calculate_dos_kp
subroutine, public calculate_dos_kp(kpoints, qs_env, dft_section)
Compute and write density of states (kpoints)
Definition: qs_dos.F:205
qs_dos::calculate_dos
subroutine, public calculate_dos(mos, dft_section)
Compute and write density of states.
Definition: qs_dos.F:57
qs_elf_methods
Does all kind of post scf calculations for GPW/GAPW.
Definition: qs_elf_methods.F:14
qs_elf_methods::qs_elf_calc
subroutine, public qs_elf_calc(qs_env, elf_r, rho_cutoff)
...
Definition: qs_elf_methods.F:54
qs_energy_window
Does all kind of post scf calculations for GPW/GAPW.
Definition: qs_energy_window.F:15
qs_energy_window::energy_windows
subroutine, public energy_windows(qs_env)
...
Definition: qs_energy_window.F:80
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_ks_types
Definition: qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition: qs_ks_types.F:330
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition: qs_ks_types.F:567
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition: qs_mo_methods.F:14
qs_mo_methods::make_mo_eig
subroutine, public make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env)
Calculate KS eigenvalues starting from OF MOS.
Definition: qs_mo_methods.F:817
qs_mo_occupation
Set occupation of molecular orbitals.
Definition: qs_mo_occupation.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_ot_eigensolver
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
Definition: qs_ot_eigensolver.F:13
qs_ot_eigensolver::ot_eigensolver
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
Definition: qs_ot_eigensolver.F:73
qs_pdos
Calculation and writing of projected density of states The DOS is computed per angular momentum and p...
Definition: qs_pdos.F:15
qs_pdos::calculate_projected_dos
subroutine, public calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, dft_section, ispin, xas_mittle, external_matrix_shalf)
Compute and write projected density of states.
Definition: qs_pdos.F:139
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition: qs_rho_methods.F:15
qs_rho_methods::qs_rho_rebuild
subroutine, public qs_rho_rebuild(rho, qs_env, rebuild_ao, rebuild_grids, admm, pw_env_external)
rebuilds rho (if necessary allocating and initializing it)
Definition: qs_rho_methods.F:96
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_set
subroutine, public qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
...
Definition: qs_rho_types.F:308
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_scf_csr_write
Functions to print the KS and S matrix in the CSR format to file.
Definition: qs_scf_csr_write.F:14
qs_scf_csr_write::write_s_matrix_csr
subroutine, public write_s_matrix_csr(qs_env, input)
writing the overlap matrix in csr format into a file
Definition: qs_scf_csr_write.F:122
qs_scf_csr_write::write_ks_matrix_csr
subroutine, public write_ks_matrix_csr(qs_env, input)
writing the KS matrix in csr format into a file
Definition: qs_scf_csr_write.F:72
qs_scf_output
Definition: qs_scf_output.F:8
qs_scf_output::qs_scf_write_mos
subroutine, public qs_scf_write_mos(qs_env, scf_env, final_mos)
Write the MO eigenvector, eigenvalues, and occupation numbers to the output unit.
Definition: qs_scf_output.F:178
qs_scf_post_tb
Does all kind of post scf calculations for DFTB.
Definition: qs_scf_post_tb.F:15
qs_scf_post_tb::scf_post_calculation_tb
subroutine, public scf_post_calculation_tb(qs_env, tb_type, no_mos)
collects possible post - scf calculations and prints info / computes properties.
Definition: qs_scf_post_tb.F:164
qs_scf_post_tb::make_lumo_tb
subroutine, public make_lumo_tb(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, nlumo, nlumos)
Gets the lumos, and eigenvalues for the lumos.
Definition: qs_scf_post_tb.F:913
qs_scf_types
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
qs_scf_types::ot_method_nr
integer, parameter, public ot_method_nr
Definition: qs_scf_types.F:51
qs_scf_wfn_mix
Does all kind of post scf calculations for GPW/GAPW.
Definition: qs_scf_wfn_mix.F:14
qs_scf_wfn_mix::wfn_mix
subroutine, public wfn_mix(mos, particle_set, dft_section, qs_kind_set, para_env, output_unit, unoccupied_orbs, scf_env, matrix_s, marked_states, for_rtp)
writes a new 'mixed' set of mos to restart file, without touching the current MOs
Definition: qs_scf_wfn_mix.F:74
qs_subsys_types
types that represent a quickstep subsys
Definition: qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition: qs_subsys_types.F:134
scf_control_types
parameters that control an scf iteration
Definition: scf_control_types.F:17
stm_images
Calculation of STM image as post processing of an electronic structure calculation,...
Definition: stm_images.F:15
stm_images::th_stm_image
subroutine, public th_stm_image(qs_env, stm_section, particles, unoccupied_orbs, unoccupied_evals)
Driver for the calculation of STM image, as post processing of a ground-state electronic structure ca...
Definition: stm_images.F:90
task_list_methods
generate the tasks lists used by collocate and integrate routines
Definition: task_list_methods.F:14
task_list_methods::generate_qs_task_list
subroutine, public generate_qs_task_list(ks_env, task_list, reorder_rs_grid_ranks, skip_load_balance_distributed, soft_valid, basis_type, pw_env_external, sab_orb_external)
...
Definition: task_list_methods.F:119
task_list_types
types for task lists
Definition: task_list_types.F:14
task_list_types::allocate_task_list
subroutine, public allocate_task_list(task_list)
allocates and initialised the components of the task_list_type
Definition: task_list_types.F:97
xtb_types
Definition of the xTB parameter types.
Definition: xtb_types.F:20
xtb_types::get_xtb_atom_param
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, occupation, electronegativity, chmax)
...
Definition: xtb_types.F:175