(git:b279b6b)
qs_initial_guess.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 Routines to somehow generate an initial guess
10 !> \par History
11 !> 2006.03 Moved here from qs_scf.F [Joost VandeVondele]
12 ! **************************************************************************************************
13 MODULE qs_initial_guess
14  USE atom_kind_orbitals, ONLY: calculate_atomic_orbitals
15  USE atomic_kind_types, ONLY: atomic_kind_type,&
16  get_atomic_kind,&
17  get_atomic_kind_set
18  USE basis_set_types, ONLY: get_gto_basis_set,&
19  gto_basis_set_type
20  USE cp_control_types, ONLY: dft_control_type
21  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
22  copy_fm_to_dbcsr,&
23  cp_dbcsr_sm_fm_multiply,&
24  cp_fm_to_dbcsr_row_template
25  USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
26  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
27  cp_fm_struct_get,&
28  cp_fm_struct_release,&
29  cp_fm_struct_type
30  USE cp_fm_types, ONLY: &
31  cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_init_random, cp_fm_release, &
32  cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type
33  USE cp_log_handling, ONLY: cp_get_default_logger,&
34  cp_logger_get_default_io_unit,&
35  cp_logger_type,&
36  cp_to_string
37  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
38  cp_print_key_unit_nr
39  USE dbcsr_api, ONLY: &
40  dbcsr_checksum, dbcsr_copy, dbcsr_dot, dbcsr_filter, dbcsr_get_diag, dbcsr_get_num_blocks, &
41  dbcsr_get_occupation, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
42  dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
43  dbcsr_nfullrows_total, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, &
44  dbcsr_set_diag, dbcsr_type, dbcsr_verify_matrix
45  USE external_potential_types, ONLY: all_potential_type,&
46  gth_potential_type,&
47  sgp_potential_type
48  USE hfx_types, ONLY: hfx_type
49  USE input_constants, ONLY: atomic_guess,&
50  core_guess,&
51  history_guess,&
52  mopac_guess,&
53  no_guess,&
54  random_guess,&
55  restart_guess,&
56  sparse_guess
57  USE input_cp2k_hfx, ONLY: ri_mo
58  USE input_section_types, ONLY: section_vals_get_subs_vals,&
59  section_vals_type,&
60  section_vals_val_get
61  USE kinds, ONLY: default_path_length,&
62  dp
63  USE kpoint_io, ONLY: read_kpoints_restart
64  USE kpoint_types, ONLY: kpoint_type
65  USE message_passing, ONLY: mp_para_env_type
66  USE particle_methods, ONLY: get_particle_set
67  USE particle_types, ONLY: particle_type
68  USE qs_atomic_block, ONLY: calculate_atomic_block_dm
69  USE qs_density_matrices, ONLY: calculate_density_matrix
70  USE qs_dftb_utils, ONLY: get_dftb_atom_param
71  USE qs_environment_types, ONLY: get_qs_env,&
72  qs_environment_type
73  USE qs_kind_types, ONLY: get_qs_kind,&
74  get_qs_kind_set,&
75  qs_kind_type
76  USE qs_mo_io, ONLY: read_mo_set_from_restart,&
77  wfn_restart_file_name
78  USE qs_mo_methods, ONLY: make_basis_lowdin,&
79  make_basis_simple,&
80  make_basis_sm
81  USE qs_mo_occupation, ONLY: set_mo_occupation
82  USE qs_mo_types, ONLY: get_mo_set,&
83  mo_set_restrict,&
84  mo_set_type,&
85  reassign_allocated_mos
86  USE qs_mom_methods, ONLY: do_mom_guess
87  USE qs_rho_methods, ONLY: qs_rho_update_rho
88  USE qs_rho_types, ONLY: qs_rho_get,&
89  qs_rho_type
90  USE qs_scf_methods, ONLY: eigensolver,&
91  eigensolver_simple
92  USE qs_scf_types, ONLY: block_davidson_diag_method_nr,&
93  block_krylov_diag_method_nr,&
94  general_diag_method_nr,&
95  ot_diag_method_nr,&
96  qs_scf_env_type
97  USE qs_wf_history_methods, ONLY: wfi_update
98  USE scf_control_types, ONLY: scf_control_type
99  USE util, ONLY: sort
100  USE xtb_types, ONLY: get_xtb_atom_param,&
101  xtb_atom_type
102 #include "./base/base_uses.f90"
103 
104  IMPLICIT NONE
105 
106  PRIVATE
107 
108  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_initial_guess'
109 
110  PUBLIC :: calculate_first_density_matrix, calculate_mopac_dm
111  PUBLIC :: calculate_atomic_fock_matrix
112 
113  TYPE atom_matrix_type
114  REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mat
115  END TYPE atom_matrix_type
116 
117 CONTAINS
118 
119 ! **************************************************************************************************
120 !> \brief can use a variety of methods to come up with an initial
121 !> density matrix and optionally an initial wavefunction
122 !> \param scf_env SCF environment information
123 !> \param qs_env QS environment
124 !> \par History
125 !> 03.2006 moved here from qs_scf [Joost VandeVondele]
126 !> 06.2007 allow to skip the initial guess [jgh]
127 !> 08.2014 kpoints [JGH]
128 !> 10.2019 tot_corr_zeff, switch_surf_dip [SGh]
129 !> \note
130 !> badly needs to be split in subroutines each doing one of the possible
131 !> schemes
132 ! **************************************************************************************************
133  SUBROUTINE calculate_first_density_matrix(scf_env, qs_env)
134 
135  TYPE(qs_scf_env_type), POINTER :: scf_env
136  TYPE(qs_environment_type), POINTER :: qs_env
137 
138  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_first_density_matrix'
139 
140  CHARACTER(LEN=default_path_length) :: file_name, filename
141  INTEGER :: atom_a, blk, density_guess, handle, homo, i, iatom, ic, icol, id_nr, ikind, irow, &
142  iseed(4), ispin, istart_col, istart_row, j, last_read, n, n_cols, n_rows, nao, natom, &
143  natoms, natoms_tmp, nblocks, nelectron, nmo, nmo_tmp, not_read, nsgf, nspin, nvec, ounit, &
144  safe_density_guess, size_atomic_kind_set, z
145  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_of, last_sgf
146  INTEGER, DIMENSION(2) :: nelectron_spin
147  INTEGER, DIMENSION(:), POINTER :: atom_list, elec_conf, nelec_kind, &
148  sort_kind
149  LOGICAL :: did_guess, do_hfx_ri_mo, do_kpoints, do_std_diag, exist, has_unit_metric, &
150  natom_mismatch, need_mos, need_wm, ofgpw, owns_ortho, print_history_log, print_log
151  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: buff, buff2
152  REAL(dp), DIMENSION(:, :), POINTER :: pdata
153  REAL(kind=dp) :: checksum, eps, length, maxocc, occ, &
154  rscale, tot_corr_zeff, trps1, zeff
155  REAL(kind=dp), DIMENSION(0:3) :: edftb
156  TYPE(atom_matrix_type), DIMENSION(:), POINTER :: pmat
157  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
158  TYPE(atomic_kind_type), POINTER :: atomic_kind
159  TYPE(cp_fm_struct_type), POINTER :: ao_ao_struct, ao_mo_struct
160  TYPE(cp_fm_type) :: sv
161  TYPE(cp_fm_type), DIMENSION(:), POINTER :: work1
162  TYPE(cp_fm_type), POINTER :: mo_coeff, moa, mob, ortho, work2
163  TYPE(cp_logger_type), POINTER :: logger
164  TYPE(dbcsr_iterator_type) :: iter
165  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: h_core_sparse, matrix_ks, p_rmpv, &
166  s_sparse
167  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h_kp, matrix_ks_kp, matrix_s_kp, &
168  rho_ao_kp
169  TYPE(dbcsr_type) :: mo_dbcsr, mo_tmp_dbcsr
170  TYPE(dft_control_type), POINTER :: dft_control
171  TYPE(gto_basis_set_type), POINTER :: orb_basis_set
172  TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
173  TYPE(kpoint_type), POINTER :: kpoints
174  TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array, mos_last_converged
175  TYPE(mp_para_env_type), POINTER :: para_env
176  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
177  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
178  TYPE(qs_kind_type), POINTER :: qs_kind
179  TYPE(qs_rho_type), POINTER :: rho
180  TYPE(scf_control_type), POINTER :: scf_control
181  TYPE(section_vals_type), POINTER :: dft_section, input, subsys_section
182 
183  logger => cp_get_default_logger()
184  NULLIFY (atomic_kind, qs_kind, mo_coeff, orb_basis_set, atomic_kind_set, &
185  qs_kind_set, particle_set, ortho, work2, work1, mo_array, s_sparse, &
186  scf_control, dft_control, p_rmpv, para_env, h_core_sparse, matrix_ks, rho, &
187  mos_last_converged)
188  NULLIFY (dft_section, input, subsys_section)
189  NULLIFY (matrix_s_kp, matrix_h_kp, matrix_ks_kp, rho_ao_kp)
190  NULLIFY (moa, mob)
191  NULLIFY (atom_list, elec_conf, kpoints)
192  edftb = 0.0_dp
193  tot_corr_zeff = 0.0_dp
194 
195  CALL timeset(routinen, handle)
196 
197  CALL get_qs_env(qs_env, &
198  atomic_kind_set=atomic_kind_set, &
199  qs_kind_set=qs_kind_set, &
200  particle_set=particle_set, &
201  mos=mo_array, &
202  matrix_s_kp=matrix_s_kp, &
203  matrix_h_kp=matrix_h_kp, &
204  matrix_ks_kp=matrix_ks_kp, &
205  input=input, &
206  scf_control=scf_control, &
207  dft_control=dft_control, &
208  has_unit_metric=has_unit_metric, &
209  do_kpoints=do_kpoints, &
210  kpoints=kpoints, &
211  rho=rho, &
212  nelectron_spin=nelectron_spin, &
213  para_env=para_env, &
214  x_data=x_data)
215 
216  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
217 
218  IF (dft_control%switch_surf_dip) THEN
219  CALL get_qs_env(qs_env, mos_last_converged=mos_last_converged)
220  END IF
221 
222  ! just initialize the first image, the other density are set to zero
223  DO ispin = 1, dft_control%nspins
224  DO ic = 1, SIZE(rho_ao_kp, 2)
225  CALL dbcsr_set(rho_ao_kp(ispin, ic)%matrix, 0.0_dp)
226  END DO
227  END DO
228  s_sparse => matrix_s_kp(:, 1)
229  h_core_sparse => matrix_h_kp(:, 1)
230  matrix_ks => matrix_ks_kp(:, 1)
231  p_rmpv => rho_ao_kp(:, 1)
232 
233  work1 => scf_env%scf_work1
234  work2 => scf_env%scf_work2
235  ortho => scf_env%ortho
236 
237  dft_section => section_vals_get_subs_vals(input, "DFT")
238 
239  nspin = dft_control%nspins
240  ofgpw = dft_control%qs_control%ofgpw
241  density_guess = scf_control%density_guess
242  do_std_diag = .false.
243 
244  do_hfx_ri_mo = .false.
245  IF (ASSOCIATED(x_data)) THEN
246  IF (x_data(1, 1)%do_hfx_ri) THEN
247  IF (x_data(1, 1)%ri_data%flavor == ri_mo) do_hfx_ri_mo = .true.
248  END IF
249  END IF
250 
251  IF (ASSOCIATED(scf_env%krylov_space)) do_std_diag = (scf_env%krylov_space%eps_std_diag > 0.0_dp)
252 
253  need_mos = scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr .OR. &
254  (scf_env%method == block_krylov_diag_method_nr .AND. .NOT. do_std_diag) &
255  .OR. dft_control%do_admm .OR. scf_env%method == block_davidson_diag_method_nr &
256  .OR. do_hfx_ri_mo
257 
258  safe_density_guess = atomic_guess
259  IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb) THEN
260  IF (density_guess == atomic_guess) density_guess = mopac_guess
261  safe_density_guess = mopac_guess
262  END IF
263  IF (dft_control%qs_control%xtb) THEN
264  IF (do_kpoints) THEN
265  IF (density_guess == atomic_guess) density_guess = mopac_guess
266  safe_density_guess = mopac_guess
267  ELSE
268  IF (density_guess == atomic_guess) density_guess = core_guess
269  safe_density_guess = core_guess
270  END IF
271  END IF
272 
273  IF (scf_control%use_ot .AND. &
274  (.NOT. ((density_guess == random_guess) .OR. &
275  (density_guess == atomic_guess) .OR. &
276  (density_guess == core_guess) .OR. &
277  (density_guess == mopac_guess) .OR. &
278  (density_guess == sparse_guess) .OR. &
279  (((density_guess == restart_guess) .OR. &
280  (density_guess == history_guess)) .AND. &
281  (scf_control%level_shift == 0.0_dp))))) THEN
282  CALL cp_abort(__location__, &
283  "OT needs GUESS ATOMIC / CORE / RANDOM / SPARSE / RESTART / HISTORY RESTART: other options NYI")
284  END IF
285 
286  ! if a restart was requested, check that the file exists,
287  ! if not we fall back to an atomic guess. No kidding, the file name should remain
288  ! in sync with read_mo_set_from_restart
289  id_nr = 0
290  IF (density_guess == restart_guess) THEN
291  ! only check existence on I/O node, otherwise if file exists there but
292  ! not on compute nodes, everything goes crazy even though only I/O
293  ! node actually reads the file
294  IF (do_kpoints) THEN
295  IF (para_env%is_source()) THEN
296  CALL wfn_restart_file_name(file_name, exist, dft_section, logger, kp=.true.)
297  END IF
298  ELSE
299  IF (para_env%is_source()) THEN
300  CALL wfn_restart_file_name(file_name, exist, dft_section, logger)
301  END IF
302  END IF
303  CALL para_env%bcast(exist)
304  CALL para_env%bcast(file_name)
305  IF (.NOT. exist) THEN
306  CALL cp_warn(__location__, &
307  "User requested to restart the wavefunction from the file named: "// &
308  trim(file_name)//". This file does not exist. Please check the existence of"// &
309  " the file or change properly the value of the keyword WFN_RESTART_FILE_NAME."// &
310  " Calculation continues using ATOMIC GUESS. ")
311  density_guess = safe_density_guess
312  END IF
313  ELSE IF (density_guess == history_guess) THEN
314  IF (do_kpoints) THEN
315  cpabort("calculate_first_density_matrix: history_guess not implemented for k-points")
316  END IF
317  IF (para_env%is_source()) THEN
318  CALL wfn_restart_file_name(file_name, exist, dft_section, logger)
319  END IF
320  CALL para_env%bcast(exist)
321  CALL para_env%bcast(file_name)
322  nvec = qs_env%wf_history%memory_depth
323  not_read = nvec + 1
324  ! At this level we read the saved backup RESTART files..
325  DO i = 1, nvec
326  j = i - 1
327  filename = trim(file_name)
328  IF (j /= 0) THEN
329  filename = trim(file_name)//".bak-"//adjustl(cp_to_string(j))
330  END IF
331  IF (para_env%is_source()) &
332  INQUIRE (file=filename, exist=exist)
333  CALL para_env%bcast(exist)
334  IF ((.NOT. exist) .AND. (i < not_read)) THEN
335  not_read = i
336  END IF
337  END DO
338  IF (not_read == 1) THEN
339  density_guess = restart_guess
340  filename = trim(file_name)
341  IF (para_env%is_source()) INQUIRE (file=filename, exist=exist)
342  CALL para_env%bcast(exist)
343  IF (.NOT. exist) THEN
344  CALL cp_warn(__location__, &
345  "User requested to restart the wavefunction from a series of restart files named: "// &
346  trim(file_name)//" with extensions (.bak-n). These files do not exist."// &
347  " Even trying to switch to a plain restart wave-function failes because the"// &
348  " file named: "//trim(file_name)//" does not exist. Please check the existence of"// &
349  " the file or change properly the value of the keyword WFN_RESTART_FILE_NAME."// &
350  " Calculation continues using ATOMIC GUESS. ")
351  density_guess = safe_density_guess
352  END IF
353  END IF
354  last_read = not_read - 1
355  END IF
356 
357  did_guess = .false.
358 
359  IF (dft_control%correct_el_density_dip) THEN
360  tot_corr_zeff = qs_env%total_zeff_corr
361  !WRITE(*,*) "tot_corr_zeff = ", tot_corr_zeff
362  IF ((abs(tot_corr_zeff) > 0.0_dp) .AND. (density_guess /= restart_guess)) THEN
363  CALL cp_warn(__location__, &
364  "Use SCF_GUESS RESTART in conjunction with "// &
365  "CORE_CORRECTION /= 0.0 and SURFACE_DIPOLE_CORRECTION TRUE. "// &
366  "It is always advisable to perform SURFACE_DIPOLE_CORRECTION "// &
367  "after a simulation without the surface dipole correction "// &
368  "and using the ensuing wavefunction restart file. ")
369  END IF
370  END IF
371 
372  ounit = -1
373  print_log = .false.
374  print_history_log = .false.
375  IF (para_env%is_source()) THEN
376  CALL section_vals_val_get(dft_section, &
377  "SCF%PRINT%RESTART%LOG_PRINT_KEY", &
378  l_val=print_log)
379  CALL section_vals_val_get(dft_section, &
380  "SCF%PRINT%RESTART_HISTORY%LOG_PRINT_KEY", &
381  l_val=print_history_log)
382  IF (print_log .OR. print_history_log) THEN
383  ounit = cp_logger_get_default_io_unit(logger)
384  END IF
385  END IF
386 
387  IF (density_guess == restart_guess) THEN
388  IF (ounit > 0) THEN
389  WRITE (unit=ounit, fmt="(/,T2,A)") &
390  "WFN_RESTART| Reading restart file"
391  END IF
392  IF (do_kpoints) THEN
393  natoms = SIZE(particle_set)
394  CALL read_kpoints_restart(rho_ao_kp, kpoints, work1, &
395  natoms, para_env, id_nr, dft_section, natom_mismatch)
396  IF (natom_mismatch) density_guess = safe_density_guess
397  ELSE
398  CALL read_mo_set_from_restart(mo_array, atomic_kind_set, qs_kind_set, particle_set, para_env, &
399  id_nr=id_nr, multiplicity=dft_control%multiplicity, dft_section=dft_section, &
400  natom_mismatch=natom_mismatch, out_unit=ounit)
401 
402  IF (natom_mismatch) THEN
403  density_guess = safe_density_guess
404  ELSE
405  DO ispin = 1, nspin
406  IF (scf_control%level_shift /= 0.0_dp) THEN
407  CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
408  CALL cp_fm_to_fm(mo_coeff, ortho)
409  END IF
410 
411  ! make all nmo vectors present orthonormal
412  CALL get_mo_set(mo_set=mo_array(ispin), &
413  mo_coeff=mo_coeff, nmo=nmo, homo=homo)
414 
415  IF (has_unit_metric) THEN
416  CALL make_basis_simple(mo_coeff, nmo)
417  ELSEIF (dft_control%smear) THEN
418  CALL make_basis_lowdin(vmatrix=mo_coeff, ncol=nmo, &
419  matrix_s=s_sparse(1)%matrix)
420  ELSE
421  ! ortho so that one can restart for different positions (basis sets?)
422  CALL make_basis_sm(mo_coeff, homo, s_sparse(1)%matrix)
423  END IF
424  ! only alpha spin is kept for restricted
425  IF (dft_control%restricted) EXIT
426  END DO
427  IF (dft_control%restricted) CALL mo_set_restrict(mo_array)
428 
429  IF (.NOT. scf_control%diagonalization%mom) THEN
430  IF (dft_control%correct_surf_dip) THEN
431  IF (abs(tot_corr_zeff) > 0.0_dp) THEN
432  CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear, &
433  tot_zeff_corr=tot_corr_zeff)
434  ELSE
435  CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear)
436  END IF
437  ELSE
438  CALL set_mo_occupation(mo_array, smear=qs_env%scf_control%smear)
439  END IF
440  END IF
441 
442  DO ispin = 1, nspin
443 
444  IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr) THEN !fm->dbcsr
445  CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
446  mo_array(ispin)%mo_coeff_b) !fm->dbcsr
447  END IF !fm->dbcsr
448 
449  CALL calculate_density_matrix(mo_array(ispin), &
450  p_rmpv(ispin)%matrix)
451  END DO
452  END IF ! natom_mismatch
453 
454  END IF
455 
456  ! Maximum Overlap Method
457  IF (scf_control%diagonalization%mom) THEN
458  CALL do_mom_guess(nspin, mo_array, scf_control, p_rmpv)
459  END IF
460 
461  did_guess = .true.
462  END IF
463 
464  IF (density_guess == history_guess) THEN
465  IF (not_read > 1) THEN
466  IF (ounit > 0) THEN
467  WRITE (unit=ounit, fmt="(/,T2,A)") &
468  "WFN_RESTART| Reading restart file history"
469  END IF
470  DO i = 1, last_read
471  j = last_read - i
472  CALL read_mo_set_from_restart(mo_array, atomic_kind_set, qs_kind_set, particle_set, para_env, &
473  id_nr=j, multiplicity=dft_control%multiplicity, &
474  dft_section=dft_section, out_unit=ounit)
475 
476  DO ispin = 1, nspin
477  IF (scf_control%level_shift /= 0.0_dp) THEN
478  CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff)
479  CALL cp_fm_to_fm(mo_coeff, ortho)
480  END IF
481 
482  ! make all nmo vectors present orthonormal
483  CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, nmo=nmo, homo=homo)
484 
485  IF (has_unit_metric) THEN
486  CALL make_basis_simple(mo_coeff, nmo)
487  ELSE
488  ! ortho so that one can restart for different positions (basis sets?)
489  CALL make_basis_sm(mo_coeff, homo, s_sparse(1)%matrix)
490  END IF
491  ! only alpha spin is kept for restricted
492  IF (dft_control%restricted) EXIT
493  END DO
494  IF (dft_control%restricted) CALL mo_set_restrict(mo_array)
495 
496  DO ispin = 1, nspin
497  CALL set_mo_occupation(mo_set=mo_array(ispin), &
498  smear=qs_env%scf_control%smear)
499  END DO
500 
501  DO ispin = 1, nspin
502  IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr) THEN !fm->dbcsr
503  CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
504  mo_array(ispin)%mo_coeff_b) !fm->dbcsr
505  END IF !fm->dbcsr
506  CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
507  END DO
508 
509  ! Write to extrapolation pipeline
510  CALL wfi_update(wf_history=qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
511  END DO
512  END IF
513 
514  did_guess = .true.
515  END IF
516 
517  IF (density_guess == random_guess) THEN
518 
519  DO ispin = 1, nspin
520  CALL get_mo_set(mo_set=mo_array(ispin), &
521  mo_coeff=mo_coeff, nmo=nmo)
522  CALL cp_fm_init_random(mo_coeff, nmo)
523  IF (has_unit_metric) THEN
524  CALL make_basis_simple(mo_coeff, nmo)
525  ELSE
526  CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
527  END IF
528  ! only alpha spin is kept for restricted
529  IF (dft_control%restricted) EXIT
530  END DO
531  IF (dft_control%restricted) CALL mo_set_restrict(mo_array)
532 
533  DO ispin = 1, nspin
534  CALL set_mo_occupation(mo_set=mo_array(ispin), &
535  smear=qs_env%scf_control%smear)
536  END DO
537 
538  DO ispin = 1, nspin
539 
540  IF (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr) THEN !fm->dbcsr
541  CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
542  mo_array(ispin)%mo_coeff_b) !fm->dbcsr
543  END IF !fm->dbcsr
544 
545  CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
546  END DO
547 
548  did_guess = .true.
549  END IF
550 
551  IF (density_guess == core_guess) THEN
552 
553  IF (do_kpoints) THEN
554  cpabort("calculate_first_density_matrix: core_guess not implemented for k-points")
555  END IF
556 
557  owns_ortho = .false.
558  IF (.NOT. ASSOCIATED(work1)) THEN
559  need_wm = .true.
560  cpassert(.NOT. ASSOCIATED(work2))
561  cpassert(.NOT. ASSOCIATED(ortho))
562  ELSE
563  need_wm = .false.
564  cpassert(ASSOCIATED(work2))
565  IF (.NOT. ASSOCIATED(ortho)) THEN
566  ALLOCATE (ortho)
567  owns_ortho = .true.
568  END IF
569  END IF
570 
571  IF (need_wm) THEN
572  CALL get_mo_set(mo_set=mo_array(1), mo_coeff=moa)
573  CALL cp_fm_get_info(moa, matrix_struct=ao_mo_struct)
574  CALL cp_fm_struct_get(ao_mo_struct, nrow_global=nao, nrow_block=nblocks)
575  CALL cp_fm_struct_create(fmstruct=ao_ao_struct, &
576  nrow_block=nblocks, &
577  ncol_block=nblocks, &
578  nrow_global=nao, &
579  ncol_global=nao, &
580  template_fmstruct=ao_mo_struct)
581  ALLOCATE (work1(1))
582  ALLOCATE (work2, ortho)
583  CALL cp_fm_create(work1(1), ao_ao_struct)
584  CALL cp_fm_create(work2, ao_ao_struct)
585  CALL cp_fm_create(ortho, ao_ao_struct)
586  CALL copy_dbcsr_to_fm(matrix_s_kp(1, 1)%matrix, ortho)
587  CALL cp_fm_cholesky_decompose(ortho)
588  CALL cp_fm_struct_release(ao_ao_struct)
589  END IF
590 
591  ispin = 1
592  ! Load core Hamiltonian into work matrix
593  CALL copy_dbcsr_to_fm(h_core_sparse(1)%matrix, work1(ispin))
594 
595  ! Diagonalize the core Hamiltonian matrix and retrieve a first set of
596  ! molecular orbitals (MOs)
597  IF (has_unit_metric) THEN
598  CALL eigensolver_simple(matrix_ks=work1(ispin), &
599  mo_set=mo_array(ispin), &
600  work=work2, &
601  do_level_shift=.false., &
602  level_shift=0.0_dp, &
603  use_jacobi=.false., jacobi_threshold=0._dp)
604  ELSE
605  CALL eigensolver(matrix_ks_fm=work1(ispin), &
606  mo_set=mo_array(ispin), &
607  ortho=ortho, &
608  work=work2, &
609  cholesky_method=scf_env%cholesky_method, &
610  do_level_shift=.false., &
611  level_shift=0.0_dp, &
612  use_jacobi=.false.)
613  END IF
614 
615  ! Open shell case: copy alpha MOs to beta MOs
616  IF (nspin == 2) THEN
617  CALL get_mo_set(mo_set=mo_array(1), mo_coeff=moa)
618  CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mob, nmo=nmo)
619  CALL cp_fm_to_fm(moa, mob, nmo)
620  END IF
621 
622  ! Build an initial density matrix (for each spin in the case of
623  ! an open shell calculation) from the first MOs set
624  DO ispin = 1, nspin
625  CALL set_mo_occupation(mo_set=mo_array(ispin), smear=scf_control%smear)
626  CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
627  END DO
628 
629  ! release intermediate matrices
630  IF (need_wm) THEN
631  CALL cp_fm_release(ortho)
632  CALL cp_fm_release(work2)
633  CALL cp_fm_release(work1(1))
634  DEALLOCATE (ortho, work2)
635  DEALLOCATE (work1)
636  NULLIFY (work1, work2, ortho)
637  ELSE IF (owns_ortho) THEN
638  DEALLOCATE (ortho)
639  END IF
640 
641  did_guess = .true.
642  END IF
643 
644  IF (density_guess == atomic_guess) THEN
645 
646  subsys_section => section_vals_get_subs_vals(input, "SUBSYS")
647  ounit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", extension=".Log")
648  IF (ounit > 0) THEN
649  WRITE (unit=ounit, fmt="(/,(T2,A))") &
650  "Atomic guess: The first density matrix is obtained in terms of atomic orbitals", &
651  " and electronic configurations assigned to each atomic kind"
652  END IF
653 
654  CALL calculate_atomic_block_dm(p_rmpv, s_sparse(1)%matrix, atomic_kind_set, qs_kind_set, &
655  nspin, nelectron_spin, ounit, para_env)
656 
657  DO ispin = 1, nspin
658 
659  ! The orbital transformation method (OT) requires not only an
660  ! initial density matrix, but also an initial wavefunction (MO set)
661  IF (ofgpw .AND. (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr)) THEN
662  ! get orbitals later
663  ELSE
664  IF (need_mos) THEN
665 
666  IF (dft_control%restricted .AND. (ispin == 2)) THEN
667  CALL mo_set_restrict(mo_array)
668  ELSE
669  CALL get_mo_set(mo_set=mo_array(ispin), &
670  mo_coeff=mo_coeff, &
671  nmo=nmo, nao=nao, homo=homo)
672 
673  CALL cp_fm_set_all(mo_coeff, 0.0_dp)
674  CALL cp_fm_init_random(mo_coeff, nmo)
675 
676  CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
677  ! multiply times PS
678  IF (has_unit_metric) THEN
679  CALL cp_fm_to_fm(mo_coeff, sv)
680  ELSE
681  ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
682  CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo)
683  END IF
684  CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
685 
686  CALL cp_fm_release(sv)
687  ! and ortho the result
688  IF (has_unit_metric) THEN
689  CALL make_basis_simple(mo_coeff, nmo)
690  ELSE
691  CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
692  END IF
693  END IF
694 
695  CALL set_mo_occupation(mo_set=mo_array(ispin), &
696  smear=qs_env%scf_control%smear)
697 
698  CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
699  mo_array(ispin)%mo_coeff_b) !fm->dbcsr
700 
701  CALL calculate_density_matrix(mo_array(ispin), &
702  p_rmpv(ispin)%matrix)
703  END IF
704  ! adjust el_density in case surface_dipole_correction is switched
705  ! on and CORE_CORRECTION is non-zero
706  IF (scf_env%method == general_diag_method_nr) THEN
707  IF (dft_control%correct_surf_dip) THEN
708  IF (abs(tot_corr_zeff) > 0.0_dp) THEN
709  CALL get_mo_set(mo_set=mo_array(ispin), &
710  mo_coeff=mo_coeff, &
711  nmo=nmo, nao=nao, homo=homo)
712 
713  CALL cp_fm_set_all(mo_coeff, 0.0_dp)
714  CALL cp_fm_init_random(mo_coeff, nmo)
715 
716  CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
717  ! multiply times PS
718  IF (has_unit_metric) THEN
719  CALL cp_fm_to_fm(mo_coeff, sv)
720  ELSE
721  ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
722  CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo)
723  END IF
724  CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
725 
726  CALL cp_fm_release(sv)
727  ! and ortho the result
728  IF (has_unit_metric) THEN
729  CALL make_basis_simple(mo_coeff, nmo)
730  ELSE
731  CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
732  END IF
733 
734  CALL set_mo_occupation(mo_set=mo_array(ispin), smear=qs_env%scf_control%smear, &
735  tot_zeff_corr=tot_corr_zeff)
736 
737  CALL calculate_density_matrix(mo_array(ispin), &
738  p_rmpv(ispin)%matrix)
739  END IF
740  END IF
741  END IF
742 
743  END IF
744 
745  END DO
746 
747  IF (ofgpw .AND. (scf_control%use_ot .OR. scf_env%method == ot_diag_method_nr)) THEN
748  ! We fit a function to the square root of the density
749  CALL qs_rho_update_rho(rho, qs_env)
750  cpassert(1 == 0)
751 ! CALL cp_fm_create(sv,mo_coeff%matrix_struct,"SV")
752 ! DO ispin=1,nspin
753 ! CALL integrate_ppl_rspace(qs%rho%rho_r(ispin),qs_env)
754 ! CALL cp_cfm_solve(overlap,mos)
755 ! CALL get_mo_set(mo_set=mo_array(ispin),&
756 ! mo_coeff=mo_coeff, nmo=nmo, nao=nao)
757 ! CALL cp_fm_init_random(mo_coeff,nmo)
758 ! END DO
759 ! CALL cp_fm_release(sv)
760  END IF
761 
762  IF (scf_control%diagonalization%mom) THEN
763  CALL do_mom_guess(nspin, mo_array, scf_control, p_rmpv)
764  END IF
765 
766  CALL cp_print_key_finished_output(ounit, logger, subsys_section, &
767  "PRINT%KINDS")
768 
769  did_guess = .true.
770  END IF
771 
772  IF (density_guess == sparse_guess) THEN
773 
774  IF (ofgpw) cpabort("SPARSE_GUESS not implemented for OFGPW")
775  IF (.NOT. scf_control%use_ot) cpabort("OT needed!")
776  IF (do_kpoints) THEN
777  cpabort("calculate_first_density_matrix: sparse_guess not implemented for k-points")
778  END IF
779 
780  eps = 1.0e-5_dp
781 
782  ounit = cp_logger_get_default_io_unit(logger)
783  natoms = SIZE(particle_set)
784  ALLOCATE (kind_of(natoms))
785  ALLOCATE (first_sgf(natoms), last_sgf(natoms))
786 
787  checksum = dbcsr_checksum(s_sparse(1)%matrix)
788  i = dbcsr_get_num_blocks(s_sparse(1)%matrix); CALL para_env%sum(i)
789  IF (ounit > 0) WRITE (ounit, *) 'S nblks', i, ' checksum', checksum
790  CALL dbcsr_filter(s_sparse(1)%matrix, eps)
791  checksum = dbcsr_checksum(s_sparse(1)%matrix)
792  i = dbcsr_get_num_blocks(s_sparse(1)%matrix); CALL para_env%sum(i)
793  IF (ounit > 0) WRITE (ounit, *) 'S nblks', i, ' checksum', checksum
794 
795  CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, &
796  last_sgf=last_sgf)
797  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
798 
799  ALLOCATE (pmat(SIZE(atomic_kind_set)))
800 
801  rscale = 1._dp
802  IF (nspin == 2) rscale = 0.5_dp
803  DO ikind = 1, SIZE(atomic_kind_set)
804  atomic_kind => atomic_kind_set(ikind)
805  qs_kind => qs_kind_set(ikind)
806  NULLIFY (pmat(ikind)%mat)
807  CALL calculate_atomic_orbitals(atomic_kind, qs_kind, pmat=pmat(ikind)%mat)
808  NULLIFY (atomic_kind)
809  END DO
810 
811  DO ispin = 1, nspin
812  CALL get_mo_set(mo_set=mo_array(ispin), &
813  maxocc=maxocc, &
814  nelectron=nelectron)
815  !
816  CALL dbcsr_iterator_start(iter, p_rmpv(ispin)%matrix)
817  DO WHILE (dbcsr_iterator_blocks_left(iter))
818  CALL dbcsr_iterator_next_block(iter, irow, icol, pdata, blk)
819  ikind = kind_of(irow)
820  IF (icol .EQ. irow) THEN
821  IF (ispin == 1) THEN
822  pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + &
823  pmat(ikind)%mat(:, :, 2)*rscale
824  ELSE
825  pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - &
826  pmat(ikind)%mat(:, :, 2)*rscale
827  END IF
828  END IF
829  END DO
830  CALL dbcsr_iterator_stop(iter)
831 
832  !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix)
833  checksum = dbcsr_checksum(p_rmpv(ispin)%matrix)
834  occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix)
835  IF (ounit > 0) WRITE (ounit, *) 'P_init occ', occ, ' checksum', checksum
836  ! so far p needs to have the same sparsity as S
837  !CALL dbcsr_filter(p_rmpv(ispin)%matrix, eps)
838  !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix)
839  checksum = dbcsr_checksum(p_rmpv(ispin)%matrix)
840  occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix)
841  IF (ounit > 0) WRITE (ounit, *) 'P_init occ', occ, ' checksum', checksum
842 
843  CALL dbcsr_dot(p_rmpv(ispin)%matrix, s_sparse(1)%matrix, trps1)
844  rscale = real(nelectron, dp)/trps1
845  CALL dbcsr_scale(p_rmpv(ispin)%matrix, rscale)
846 
847  !CALL dbcsr_verify_matrix(p_rmpv(ispin)%matrix)
848  checksum = dbcsr_checksum(p_rmpv(ispin)%matrix)
849  occ = dbcsr_get_occupation(p_rmpv(ispin)%matrix)
850  IF (ounit > 0) WRITE (ounit, *) 'P occ', occ, ' checksum', checksum
851  !
852  ! The orbital transformation method (OT) requires not only an
853  ! initial density matrix, but also an initial wavefunction (MO set)
854  IF (dft_control%restricted .AND. (ispin == 2)) THEN
855  CALL mo_set_restrict(mo_array)
856  ELSE
857  CALL get_mo_set(mo_set=mo_array(ispin), &
858  mo_coeff=mo_coeff, &
859  nmo=nmo, nao=nao, homo=homo)
860  CALL cp_fm_set_all(mo_coeff, 0.0_dp)
861 
862  n = maxval(last_sgf - first_sgf) + 1
863  size_atomic_kind_set = SIZE(atomic_kind_set)
864 
865  ALLOCATE (buff(n, n), sort_kind(size_atomic_kind_set), &
866  nelec_kind(size_atomic_kind_set))
867  !
868  ! sort kind vs nbr electron
869  DO ikind = 1, size_atomic_kind_set
870  atomic_kind => atomic_kind_set(ikind)
871  qs_kind => qs_kind_set(ikind)
872  CALL get_atomic_kind(atomic_kind=atomic_kind, &
873  natom=natom, &
874  atom_list=atom_list, &
875  z=z)
876  CALL get_qs_kind(qs_kind, nsgf=nsgf, elec_conf=elec_conf, &
877  basis_set=orb_basis_set, zeff=zeff)
878  nelec_kind(ikind) = sum(elec_conf)
879  END DO
880  CALL sort(nelec_kind, size_atomic_kind_set, sort_kind)
881  !
882  ! a -very- naive sparse guess
883  nmo_tmp = nmo
884  natoms_tmp = natoms
885  istart_col = 1
886  iseed(1) = 4; iseed(2) = 3; iseed(3) = 2; iseed(4) = 1 ! set the seed for dlarnv
887  DO i = 1, size_atomic_kind_set
888  ikind = sort_kind(i)
889  atomic_kind => atomic_kind_set(ikind)
890  CALL get_atomic_kind(atomic_kind=atomic_kind, &
891  natom=natom, atom_list=atom_list)
892  DO iatom = 1, natom
893  !
894  atom_a = atom_list(iatom)
895  istart_row = first_sgf(atom_a)
896  n_rows = last_sgf(atom_a) - first_sgf(atom_a) + 1
897  !
898  ! compute the "potential" nbr of states for this atom
899  n_cols = max(int(real(nmo_tmp, dp)/real(natoms_tmp, dp)), 1)
900  IF (n_cols .GT. n_rows) n_cols = n_rows
901  !
902  nmo_tmp = nmo_tmp - n_cols
903  natoms_tmp = natoms_tmp - 1
904  IF (nmo_tmp .LT. 0 .OR. natoms_tmp .LT. 0) THEN
905  cpabort("Wrong1!")
906  END IF
907  DO j = 1, n_cols
908  CALL dlarnv(1, iseed, n_rows, buff(1, j))
909  END DO
910  CALL cp_fm_set_submatrix(mo_coeff, buff, istart_row, istart_col, &
911  n_rows, n_cols)
912  istart_col = istart_col + n_cols
913  END DO
914  END DO
915 
916  IF (istart_col .LE. nmo) THEN
917  cpabort("Wrong2!")
918  END IF
919 
920  DEALLOCATE (buff, nelec_kind, sort_kind)
921 
922  IF (.false.) THEN
923  ALLOCATE (buff(nao, 1), buff2(nao, 1))
924  DO i = 1, nmo
925  CALL cp_fm_get_submatrix(mo_coeff, buff, 1, i, nao, 1)
926  IF (sum(buff**2) .LT. 1e-10_dp) THEN
927  WRITE (*, *) 'wrong', i, sum(buff**2)
928  END IF
929  length = sqrt(dot_product(buff(:, 1), buff(:, 1)))
930  buff(:, :) = buff(:, :)/length
931  DO j = i + 1, nmo
932  CALL cp_fm_get_submatrix(mo_coeff, buff2, 1, j, nao, 1)
933  length = sqrt(dot_product(buff2(:, 1), buff2(:, 1)))
934  buff2(:, :) = buff2(:, :)/length
935  IF (abs(dot_product(buff(:, 1), buff2(:, 1)) - 1.0_dp) .LT. 1e-10_dp) THEN
936  WRITE (*, *) 'wrong2', i, j, dot_product(buff(:, 1), buff2(:, 1))
937  DO ikind = 1, nao
938  IF (abs(mo_coeff%local_data(ikind, i)) .GT. 1e-10_dp) THEN
939  WRITE (*, *) 'c1', ikind, mo_coeff%local_data(ikind, i)
940  END IF
941  IF (abs(mo_coeff%local_data(ikind, j)) .GT. 1e-10_dp) THEN
942  WRITE (*, *) 'c2', ikind, mo_coeff%local_data(ikind, j)
943  END IF
944  END DO
945  cpabort("")
946  END IF
947  END DO
948  END DO
949  DEALLOCATE (buff, buff2)
950 
951  END IF
952  !
953  CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, mo_coeff, s_sparse(1)%matrix)
954  !CALL dbcsr_verify_matrix(mo_dbcsr)
955  checksum = dbcsr_checksum(mo_dbcsr)
956 
957  occ = dbcsr_get_occupation(mo_dbcsr)
958  IF (ounit > 0) WRITE (ounit, *) 'C occ', occ, ' checksum', checksum
959  CALL dbcsr_filter(mo_dbcsr, eps)
960  !CALL dbcsr_verify_matrix(mo_dbcsr)
961  occ = dbcsr_get_occupation(mo_dbcsr)
962  checksum = dbcsr_checksum(mo_dbcsr)
963  IF (ounit > 0) WRITE (ounit, *) 'C occ', occ, ' checksum', checksum
964  !
965  ! multiply times PS
966  IF (has_unit_metric) THEN
967  cpabort("has_unit_metric will be removed soon")
968  END IF
969  !
970  ! S*C
971  CALL dbcsr_copy(mo_tmp_dbcsr, mo_dbcsr, name="mo_tmp")
972  CALL dbcsr_multiply("N", "N", 1.0_dp, s_sparse(1)%matrix, mo_dbcsr, &
973  0.0_dp, mo_tmp_dbcsr, &
974  retain_sparsity=.true.)
975  !CALL dbcsr_verify_matrix(mo_tmp_dbcsr)
976  checksum = dbcsr_checksum(mo_tmp_dbcsr)
977  occ = dbcsr_get_occupation(mo_tmp_dbcsr)
978  IF (ounit > 0) WRITE (ounit, *) 'S*C occ', occ, ' checksum', checksum
979  CALL dbcsr_filter(mo_tmp_dbcsr, eps)
980  !CALL dbcsr_verify_matrix(mo_tmp_dbcsr)
981  checksum = dbcsr_checksum(mo_tmp_dbcsr)
982  occ = dbcsr_get_occupation(mo_tmp_dbcsr)
983  IF (ounit > 0) WRITE (ounit, *) 'S*C occ', occ, ' checksum', checksum
984  !
985  ! P*SC
986  ! the destroy is needed for the moment to avoid memory leaks !
987  ! This one is not needed because _destroy takes care of zeroing.
988  CALL dbcsr_multiply("N", "N", 1.0_dp, p_rmpv(ispin)%matrix, &
989  mo_tmp_dbcsr, 0.0_dp, mo_dbcsr)
990  IF (.false.) CALL dbcsr_verify_matrix(mo_dbcsr)
991  checksum = dbcsr_checksum(mo_dbcsr)
992  occ = dbcsr_get_occupation(mo_dbcsr)
993  IF (ounit > 0) WRITE (ounit, *) 'P*SC occ', occ, ' checksum', checksum
994  CALL dbcsr_filter(mo_dbcsr, eps)
995  !CALL dbcsr_verify_matrix(mo_dbcsr)
996  checksum = dbcsr_checksum(mo_dbcsr)
997  occ = dbcsr_get_occupation(mo_dbcsr)
998  IF (ounit > 0) WRITE (ounit, *) 'P*SC occ', occ, ' checksum', checksum
999  !
1000  CALL copy_dbcsr_to_fm(mo_dbcsr, mo_coeff)
1001 
1002  CALL dbcsr_release(mo_dbcsr)
1003  CALL dbcsr_release(mo_tmp_dbcsr)
1004 
1005  ! and ortho the result
1006  CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
1007  END IF
1008 
1009  CALL set_mo_occupation(mo_set=mo_array(ispin), &
1010  smear=qs_env%scf_control%smear)
1011 
1012  CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
1013  mo_array(ispin)%mo_coeff_b) !fm->dbcsr
1014 
1015  CALL calculate_density_matrix(mo_array(ispin), &
1016  p_rmpv(ispin)%matrix)
1017  DO ikind = 1, SIZE(atomic_kind_set)
1018  IF (ASSOCIATED(pmat(ikind)%mat)) THEN
1019  DEALLOCATE (pmat(ikind)%mat)
1020  END IF
1021  END DO
1022  END DO
1023 
1024  DEALLOCATE (pmat)
1025 
1026  DEALLOCATE (kind_of)
1027 
1028  DEALLOCATE (first_sgf, last_sgf)
1029 
1030  did_guess = .true.
1031  END IF
1032  IF (density_guess == mopac_guess) THEN
1033 
1034  CALL calculate_mopac_dm(p_rmpv, s_sparse(1)%matrix, has_unit_metric, dft_control, &
1035  particle_set, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, para_env)
1036 
1037  DO ispin = 1, nspin
1038  ! The orbital transformation method (OT) requires not only an
1039  ! initial density matrix, but also an initial wavefunction (MO set)
1040  IF (need_mos) THEN
1041  IF (dft_control%restricted .AND. (ispin == 2)) THEN
1042  CALL mo_set_restrict(mo_array)
1043  ELSE
1044  CALL get_mo_set(mo_set=mo_array(ispin), &
1045  mo_coeff=mo_coeff, &
1046  nmo=nmo, homo=homo)
1047  CALL cp_fm_init_random(mo_coeff, nmo)
1048  CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
1049  ! multiply times PS
1050  IF (has_unit_metric) THEN
1051  CALL cp_fm_to_fm(mo_coeff, sv)
1052  ELSE
1053  CALL cp_dbcsr_sm_fm_multiply(s_sparse(1)%matrix, mo_coeff, sv, nmo)
1054  END IF
1055  ! here we could easily multiply with the diag that we actually have replicated already
1056  CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
1057  CALL cp_fm_release(sv)
1058  ! and ortho the result
1059  IF (has_unit_metric) THEN
1060  CALL make_basis_simple(mo_coeff, nmo)
1061  ELSE
1062  CALL make_basis_sm(mo_coeff, nmo, s_sparse(1)%matrix)
1063  END IF
1064  END IF
1065 
1066  CALL set_mo_occupation(mo_set=mo_array(ispin), &
1067  smear=qs_env%scf_control%smear)
1068  CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, &
1069  mo_array(ispin)%mo_coeff_b)
1070 
1071  CALL calculate_density_matrix(mo_array(ispin), &
1072  p_rmpv(ispin)%matrix)
1073  END IF
1074  END DO
1075 
1076  did_guess = .true.
1077  END IF
1078 ! switch_surf_dip [SGh]
1079  IF (dft_control%switch_surf_dip) THEN
1080  DO ispin = 1, nspin
1081  CALL reassign_allocated_mos(mos_last_converged(ispin), &
1082  mo_array(ispin))
1083  END DO
1084  END IF
1085 
1086  IF (density_guess == no_guess) THEN
1087  did_guess = .true.
1088  END IF
1089 
1090  IF (.NOT. did_guess) THEN
1091  cpabort("An invalid keyword for the initial density guess was specified")
1092  END IF
1093 
1094  CALL timestop(handle)
1095 
1096  END SUBROUTINE calculate_first_density_matrix
1097 
1098 ! **************************************************************************************************
1099 !> \brief returns a block diagonal fock matrix.
1100 !> \param matrix_f ...
1101 !> \param atomic_kind_set ...
1102 !> \param qs_kind_set ...
1103 !> \param ounit ...
1104 ! **************************************************************************************************
1105  SUBROUTINE calculate_atomic_fock_matrix(matrix_f, atomic_kind_set, qs_kind_set, ounit)
1106  TYPE(dbcsr_type), INTENT(INOUT) :: matrix_f
1107  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1108  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1109  INTEGER, INTENT(IN) :: ounit
1110 
1111  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_atomic_fock_matrix'
1112 
1113  INTEGER :: handle, icol, ikind, irow
1114  INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
1115  REAL(dp), DIMENSION(:, :), POINTER :: block
1116  TYPE(atom_matrix_type), ALLOCATABLE, DIMENSION(:) :: fmat
1117  TYPE(atomic_kind_type), POINTER :: atomic_kind
1118  TYPE(dbcsr_iterator_type) :: iter
1119  TYPE(qs_kind_type), POINTER :: qs_kind
1120 
1121  CALL timeset(routinen, handle)
1122 
1123  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
1124  ALLOCATE (fmat(SIZE(atomic_kind_set)))
1125 
1126  ! precompute the atomic blocks for each atomic-kind
1127  DO ikind = 1, SIZE(atomic_kind_set)
1128  atomic_kind => atomic_kind_set(ikind)
1129  qs_kind => qs_kind_set(ikind)
1130  NULLIFY (fmat(ikind)%mat)
1131  IF (ounit > 0) WRITE (unit=ounit, fmt="(/,T2,A)") &
1132  "Calculating atomic Fock matrix for atomic kind: "//trim(atomic_kind%name)
1133 
1134  !Currently only ispin=1 is supported
1135  CALL calculate_atomic_orbitals(atomic_kind, qs_kind, iunit=ounit, &
1136  fmat=fmat(ikind)%mat)
1137  END DO
1138 
1139  ! zero result matrix
1140  CALL dbcsr_set(matrix_f, 0.0_dp)
1141 
1142  ! copy precomputed blocks onto diagonal of result matrix
1143  CALL dbcsr_iterator_start(iter, matrix_f)
1144  DO WHILE (dbcsr_iterator_blocks_left(iter))
1145  CALL dbcsr_iterator_next_block(iter, irow, icol, block)
1146  ikind = kind_of(irow)
1147  IF (icol .EQ. irow) block(:, :) = fmat(ikind)%mat(:, :, 1)
1148  END DO
1149  CALL dbcsr_iterator_stop(iter)
1150 
1151  ! cleanup
1152  DO ikind = 1, SIZE(atomic_kind_set)
1153  DEALLOCATE (fmat(ikind)%mat)
1154  END DO
1155  DEALLOCATE (fmat)
1156 
1157  CALL timestop(handle)
1158 
1159  END SUBROUTINE calculate_atomic_fock_matrix
1160 
1161 ! **************************************************************************************************
1162 !> \brief returns a block diagonal density matrix. Blocks correspond to the mopac initial guess.
1163 !> \param pmat ...
1164 !> \param matrix_s ...
1165 !> \param has_unit_metric ...
1166 !> \param dft_control ...
1167 !> \param particle_set ...
1168 !> \param atomic_kind_set ...
1169 !> \param qs_kind_set ...
1170 !> \param nspin ...
1171 !> \param nelectron_spin ...
1172 !> \param para_env ...
1173 ! **************************************************************************************************
1174  SUBROUTINE calculate_mopac_dm(pmat, matrix_s, has_unit_metric, &
1175  dft_control, particle_set, atomic_kind_set, qs_kind_set, &
1176  nspin, nelectron_spin, para_env)
1177  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: pmat
1178  TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s
1179  LOGICAL :: has_unit_metric
1180  TYPE(dft_control_type), POINTER :: dft_control
1181  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1182  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1183  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1184  INTEGER, INTENT(IN) :: nspin
1185  INTEGER, DIMENSION(:), INTENT(IN) :: nelectron_spin
1186  TYPE(mp_para_env_type) :: para_env
1187 
1188  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_mopac_dm'
1189 
1190  INTEGER :: atom_a, handle, iatom, ikind, iset, &
1191  isgf, isgfa, ishell, ispin, la, maxl, &
1192  maxll, na, nao, natom, ncount, nset, &
1193  nsgf, z
1194  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf
1195  INTEGER, DIMENSION(25) :: laox, naox
1196  INTEGER, DIMENSION(5) :: occupation
1197  INTEGER, DIMENSION(:), POINTER :: atom_list, elec_conf, nshell
1198  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, l, last_sgfa
1199  LOGICAL :: has_pot
1200  REAL(kind=dp) :: maxocc, my_sum, nelec, occ, paa, rscale, &
1201  trps1, trps2, yy, zeff
1202  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: econf, pdiag, sdiag
1203  REAL(kind=dp), DIMENSION(0:3) :: edftb
1204  TYPE(all_potential_type), POINTER :: all_potential
1205  TYPE(dbcsr_type), POINTER :: matrix_p
1206  TYPE(gth_potential_type), POINTER :: gth_potential
1207  TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1208  TYPE(sgp_potential_type), POINTER :: sgp_potential
1209  TYPE(xtb_atom_type), POINTER :: xtb_kind
1210 
1211  CALL timeset(routinen, handle)
1212 
1213  DO ispin = 1, nspin
1214  matrix_p => pmat(ispin)%matrix
1215  CALL dbcsr_set(matrix_p, 0.0_dp)
1216  END DO
1217 
1218  natom = SIZE(particle_set)
1219  nao = dbcsr_nfullrows_total(pmat(1)%matrix)
1220  IF (nspin == 1) THEN
1221  maxocc = 2.0_dp
1222  ELSE
1223  maxocc = 1.0_dp
1224  END IF
1225 
1226  ALLOCATE (first_sgf(natom))
1227 
1228  CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf)
1229  CALL get_qs_kind_set(qs_kind_set, maxlgto=maxl)
1230 
1231  ALLOCATE (econf(0:maxl))
1232 
1233  ALLOCATE (pdiag(nao))
1234  pdiag(:) = 0.0_dp
1235 
1236  ALLOCATE (sdiag(nao))
1237 
1238  sdiag(:) = 0.0_dp
1239  IF (has_unit_metric) THEN
1240  sdiag(:) = 1.0_dp
1241  ELSE
1242  CALL dbcsr_get_diag(matrix_s, sdiag)
1243  CALL para_env%sum(sdiag)
1244  END IF
1245 
1246  ncount = 0
1247  trps1 = 0.0_dp
1248  trps2 = 0.0_dp
1249  pdiag(:) = 0.0_dp
1250 
1251  IF (sum(nelectron_spin) /= 0) THEN
1252  DO ikind = 1, SIZE(atomic_kind_set)
1253 
1254  CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
1255  CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
1256  all_potential=all_potential, &
1257  gth_potential=gth_potential, &
1258  sgp_potential=sgp_potential)
1259  has_pot = ASSOCIATED(all_potential) .OR. ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)
1260 
1261  IF (dft_control%qs_control%dftb) THEN
1262  CALL get_dftb_atom_param(qs_kind_set(ikind)%dftb_parameter, &
1263  lmax=maxll, occupation=edftb)
1264  maxll = min(maxll, maxl)
1265  econf(0:maxl) = edftb(0:maxl)
1266  ELSEIF (dft_control%qs_control%xtb) THEN
1267  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1268  CALL get_xtb_atom_param(xtb_kind, z=z, natorb=nsgf, nao=naox, lao=laox, occupation=occupation)
1269  ELSEIF (has_pot) THEN
1270  CALL get_atomic_kind(atomic_kind_set(ikind), z=z)
1271  CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf, elec_conf=elec_conf, zeff=zeff)
1272  maxll = min(SIZE(elec_conf) - 1, maxl)
1273  econf(:) = 0.0_dp
1274  econf(0:maxll) = 0.5_dp*maxocc*real(elec_conf(0:maxll), dp)
1275  ELSE
1276  cycle
1277  END IF
1278 
1279  ! MOPAC TYPE GUESS
1280  IF (dft_control%qs_control%dftb) THEN
1281  DO iatom = 1, natom
1282  atom_a = atom_list(iatom)
1283  isgfa = first_sgf(atom_a)
1284  DO la = 0, maxll
1285  SELECT CASE (la)
1286  CASE (0)
1287  pdiag(isgfa) = econf(0)
1288  CASE (1)
1289  pdiag(isgfa + 1) = econf(1)/3._dp
1290  pdiag(isgfa + 2) = econf(1)/3._dp
1291  pdiag(isgfa + 3) = econf(1)/3._dp
1292  CASE (2)
1293  pdiag(isgfa + 4) = econf(2)/5._dp
1294  pdiag(isgfa + 5) = econf(2)/5._dp
1295  pdiag(isgfa + 6) = econf(2)/5._dp
1296  pdiag(isgfa + 7) = econf(2)/5._dp
1297  pdiag(isgfa + 8) = econf(2)/5._dp
1298  CASE (3)
1299  pdiag(isgfa + 9) = econf(3)/7._dp
1300  pdiag(isgfa + 10) = econf(3)/7._dp
1301  pdiag(isgfa + 11) = econf(3)/7._dp
1302  pdiag(isgfa + 12) = econf(3)/7._dp
1303  pdiag(isgfa + 13) = econf(3)/7._dp
1304  pdiag(isgfa + 14) = econf(3)/7._dp
1305  pdiag(isgfa + 15) = econf(3)/7._dp
1306  CASE DEFAULT
1307  cpabort("")
1308  END SELECT
1309  END DO
1310  END DO
1311  ELSEIF (dft_control%qs_control%xtb) THEN
1312  DO iatom = 1, natom
1313  atom_a = atom_list(iatom)
1314  isgfa = first_sgf(atom_a)
1315  IF (z == 1 .AND. nsgf == 2) THEN
1316  ! Hydrogen 2s basis
1317  pdiag(isgfa) = 1.0_dp
1318  pdiag(isgfa + 1) = 0.0_dp
1319  ELSE
1320  DO isgf = 1, nsgf
1321  na = naox(isgf)
1322  la = laox(isgf)
1323  occ = real(occupation(la + 1), dp)/real(2*la + 1, dp)
1324  pdiag(isgfa + isgf - 1) = occ
1325  END DO
1326  END IF
1327  END DO
1328  ELSEIF (dft_control%qs_control%semi_empirical) THEN
1329  yy = real(dft_control%charge, kind=dp)/real(nao, kind=dp)
1330  DO iatom = 1, natom
1331  atom_a = atom_list(iatom)
1332  isgfa = first_sgf(atom_a)
1333  SELECT CASE (nsgf)
1334  CASE (1) ! s-basis
1335  pdiag(isgfa) = (zeff - yy)*0.5_dp*maxocc
1336  CASE (4) ! sp-basis
1337  IF (z == 1) THEN
1338  ! special case: hydrogen with sp basis
1339  pdiag(isgfa) = (zeff - yy)*0.5_dp*maxocc
1340  pdiag(isgfa + 1) = 0._dp
1341  pdiag(isgfa + 2) = 0._dp
1342  pdiag(isgfa + 3) = 0._dp
1343  ELSE
1344  pdiag(isgfa) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1345  pdiag(isgfa + 1) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1346  pdiag(isgfa + 2) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1347  pdiag(isgfa + 3) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1348  END IF
1349  CASE (9) ! spd-basis
1350  IF (z < 21 .OR. z > 30 .AND. z < 39 .OR. z > 48 .AND. z < 57) THEN
1351  ! Main Group Element: The "d" shell is formally empty.
1352  pdiag(isgfa) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1353  pdiag(isgfa + 1) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1354  pdiag(isgfa + 2) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1355  pdiag(isgfa + 3) = (zeff*0.25_dp - yy)*0.5_dp*maxocc
1356  pdiag(isgfa + 4) = (-yy)*0.5_dp*maxocc
1357  pdiag(isgfa + 5) = (-yy)*0.5_dp*maxocc
1358  pdiag(isgfa + 6) = (-yy)*0.5_dp*maxocc
1359  pdiag(isgfa + 7) = (-yy)*0.5_dp*maxocc
1360  pdiag(isgfa + 8) = (-yy)*0.5_dp*maxocc
1361  ELSE IF (z < 99) THEN
1362  my_sum = zeff - 9.0_dp*yy
1363  ! First, put 2 electrons in the 's' shell
1364  pdiag(isgfa) = (max(0.0_dp, min(my_sum, 2.0_dp)))*0.5_dp*maxocc
1365  my_sum = my_sum - 2.0_dp
1366  IF (my_sum > 0.0_dp) THEN
1367  ! Now put as many electrons as possible into the 'd' shell
1368  pdiag(isgfa + 4) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1369  pdiag(isgfa + 5) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1370  pdiag(isgfa + 6) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1371  pdiag(isgfa + 7) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1372  pdiag(isgfa + 8) = (max(0.0_dp, min(my_sum*0.2_dp, 2.0_dp)))*0.5_dp*maxocc
1373  my_sum = max(0.0_dp, my_sum - 10.0_dp)
1374  ! Put the remaining electrons in the 'p' shell
1375  pdiag(isgfa + 1) = (my_sum/3.0_dp)*0.5_dp*maxocc
1376  pdiag(isgfa + 2) = (my_sum/3.0_dp)*0.5_dp*maxocc
1377  pdiag(isgfa + 3) = (my_sum/3.0_dp)*0.5_dp*maxocc
1378  END IF
1379  END IF
1380  CASE DEFAULT
1381  cpabort("")
1382  END SELECT
1383  END DO
1384  ELSE
1385  CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1386  nset=nset, &
1387  nshell=nshell, &
1388  l=l, &
1389  first_sgf=first_sgfa, &
1390  last_sgf=last_sgfa)
1391 
1392  DO iset = 1, nset
1393  DO ishell = 1, nshell(iset)
1394  la = l(ishell, iset)
1395  nelec = maxocc*real(2*la + 1, dp)
1396  IF (econf(la) > 0.0_dp) THEN
1397  IF (econf(la) >= nelec) THEN
1398  paa = maxocc
1399  econf(la) = econf(la) - nelec
1400  ELSE
1401  paa = maxocc*econf(la)/nelec
1402  econf(la) = 0.0_dp
1403  ncount = ncount + nint(nelec/maxocc)
1404  END IF
1405  DO isgfa = first_sgfa(ishell, iset), last_sgfa(ishell, iset)
1406  DO iatom = 1, natom
1407  atom_a = atom_list(iatom)
1408  isgf = first_sgf(atom_a) + isgfa - 1
1409  pdiag(isgf) = paa
1410  IF (paa == maxocc) THEN
1411  trps1 = trps1 + paa*sdiag(isgf)
1412  ELSE
1413  trps2 = trps2 + paa*sdiag(isgf)
1414  END IF
1415  END DO
1416  END DO
1417  END IF
1418  END DO ! ishell
1419  END DO ! iset
1420  END IF
1421  END DO ! ikind
1422 
1423  IF (trps2 == 0.0_dp) THEN
1424  DO isgf = 1, nao
1425  IF (sdiag(isgf) > 0.0_dp) pdiag(isgf) = pdiag(isgf)/sdiag(isgf)
1426  END DO
1427  DO ispin = 1, nspin
1428  IF (nelectron_spin(ispin) /= 0) THEN
1429  matrix_p => pmat(ispin)%matrix
1430  CALL dbcsr_set_diag(matrix_p, pdiag)
1431  END IF
1432  END DO
1433  ELSE
1434  DO ispin = 1, nspin
1435  IF (nelectron_spin(ispin) /= 0) THEN
1436  rscale = (real(nelectron_spin(ispin), dp) - trps1)/trps2
1437  DO isgf = 1, nao
1438  IF (pdiag(isgf) < maxocc) pdiag(isgf) = rscale*pdiag(isgf)
1439  END DO
1440  matrix_p => pmat(ispin)%matrix
1441  CALL dbcsr_set_diag(matrix_p, pdiag)
1442  DO isgf = 1, nao
1443  IF (pdiag(isgf) < maxocc) pdiag(isgf) = pdiag(isgf)/rscale
1444  END DO
1445  END IF
1446  END DO
1447  END IF
1448  END IF
1449 
1450  DEALLOCATE (econf)
1451 
1452  DEALLOCATE (first_sgf)
1453 
1454  DEALLOCATE (pdiag)
1455 
1456  DEALLOCATE (sdiag)
1457 
1458  CALL timestop(handle)
1459 
1460  END SUBROUTINE calculate_mopac_dm
1461 
1462 END MODULE qs_initial_guess
atom_kind_orbitals
calculate the orbitals for a given atomic kind type
Definition: atom_kind_orbitals.F:11
atom_kind_orbitals::calculate_atomic_orbitals
subroutine, public calculate_atomic_orbitals(atomic_kind, qs_kind, agrid, iunit, pmat, fmat, density, wavefunction, wfninfo, confine, xc_section, nocc)
...
Definition: atom_kind_orbitals.F:89
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
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
basis_set_types
Definition: basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Definition: basis_set_types.F:615
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::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
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_operations::cp_fm_to_dbcsr_row_template
subroutine, public cp_fm_to_dbcsr_row_template(matrix, fm_in, template)
Utility function to copy a specially shaped fm to dbcsr_matrix The result matrix will be the matrix i...
Definition: cp_dbcsr_operations.F:1070
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition: cp_dbcsr_operations.F:118
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_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_get
subroutine, public cp_fm_struct_get(fmstruct, para_env, context, descriptor, ncol_block, nrow_block, nrow_global, ncol_global, first_p_pos, row_indices, col_indices, nrow_local, ncol_local, nrow_locals, ncol_locals, local_leading_dimension)
returns the values of various attributes of the matrix structure
Definition: cp_fm_struct.F:409
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_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition: cp_fm_types.F:768
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition: cp_fm_types.F:901
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
external_potential_types
Definition of the atomic potential types.
Definition: external_potential_types.F:14
hfx_types
Types and set/get functions for HFX.
Definition: hfx_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::core_guess
integer, parameter, public core_guess
Definition: input_constants.F:201
input_constants::mopac_guess
integer, parameter, public mopac_guess
Definition: input_constants.F:201
input_constants::no_guess
integer, parameter, public no_guess
Definition: input_constants.F:201
input_constants::atomic_guess
integer, parameter, public atomic_guess
Definition: input_constants.F:201
input_constants::history_guess
integer, parameter, public history_guess
Definition: input_constants.F:201
input_constants::random_guess
integer, parameter, public random_guess
Definition: input_constants.F:201
input_constants::sparse_guess
integer, parameter, public sparse_guess
Definition: input_constants.F:201
input_constants::restart_guess
integer, parameter, public restart_guess
Definition: input_constants.F:201
input_cp2k_hfx
function that builds the hartree fock exchange section of the input
Definition: input_cp2k_hfx.F:14
input_cp2k_hfx::ri_mo
integer, parameter, public ri_mo
Definition: input_cp2k_hfx.F:44
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_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_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
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
kpoint_io
Restart file for k point calculations.
Definition: kpoint_io.F:13
kpoint_io::read_kpoints_restart
subroutine, public read_kpoints_restart(denmat, kpoints, fmwork, natom, para_env, id_nr, dft_section, natom_mismatch)
...
Definition: kpoint_io.F:297
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
particle_methods
Define methods related to particle_type.
Definition: particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition: particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_atomic_block
Routine to return block diagonal density matrix. Blocks correspond to the atomic densities.
Definition: qs_atomic_block.F:14
qs_atomic_block::calculate_atomic_block_dm
subroutine, public calculate_atomic_block_dm(pmatrix, matrix_s, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, ounit, para_env)
returns a block diagonal density matrix. Blocks correspond to the atomic densities.
Definition: qs_atomic_block.F:54
qs_density_matrices
collects routines that calculate density matrices
Definition: qs_density_matrices.F:14
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_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_initial_guess
Routines to somehow generate an initial guess.
Definition: qs_initial_guess.F:13
qs_initial_guess::calculate_first_density_matrix
subroutine, public calculate_first_density_matrix(scf_env, qs_env)
can use a variety of methods to come up with an initial density matrix and optionally an initial wave...
Definition: qs_initial_guess.F:134
qs_initial_guess::calculate_mopac_dm
subroutine, public calculate_mopac_dm(pmat, matrix_s, has_unit_metric, dft_control, particle_set, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, para_env)
returns a block diagonal density matrix. Blocks correspond to the mopac initial guess.
Definition: qs_initial_guess.F:1177
qs_initial_guess::calculate_atomic_fock_matrix
subroutine, public calculate_atomic_fock_matrix(matrix_f, atomic_kind_set, qs_kind_set, ounit)
returns a block diagonal fock matrix.
Definition: qs_initial_guess.F:1106
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_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_mo_io
Definition and initialisation of the mo data type.
Definition: qs_mo_io.F:21
qs_mo_io::wfn_restart_file_name
subroutine, public wfn_restart_file_name(filename, exist, section, logger, kp, xas, rtp)
...
Definition: qs_mo_io.F:428
qs_mo_io::read_mo_set_from_restart
subroutine, public read_mo_set_from_restart(mo_array, atomic_kind_set, qs_kind_set, particle_set, para_env, id_nr, multiplicity, dft_section, natom_mismatch, cdft, out_unit)
...
Definition: qs_mo_io.F:495
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition: qs_mo_methods.F:14
qs_mo_methods::make_basis_simple
subroutine, public make_basis_simple(vmatrix, ncol)
given a set of vectors, return an orthogonal (C^T C == 1) set spanning the same space (notice,...
Definition: qs_mo_methods.F:352
qs_mo_methods::make_basis_lowdin
subroutine, public make_basis_lowdin(vmatrix, ncol, matrix_s)
return a set of S orthonormal vectors (C^T S C == 1) where a Loedwin transformation is applied to kee...
Definition: qs_mo_methods.F:298
qs_mo_methods::make_basis_sm
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition: qs_mo_methods.F:87
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::mo_set_restrict
subroutine, public mo_set_restrict(mo_array, convert_dbcsr)
make the beta orbitals explicitly equal to the alpha orbitals effectively copying the orbital data
Definition: qs_mo_types.F:315
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_mo_types::reassign_allocated_mos
subroutine, public reassign_allocated_mos(mo_set_new, mo_set_old)
reassign an already allocated mo_set
Definition: qs_mo_types.F:109
qs_mom_methods
methods for deltaSCF calculations
Definition: qs_mom_methods.F:11
qs_mom_methods::do_mom_guess
subroutine, public do_mom_guess(nspins, mos, scf_control, p_rmpv)
initial guess for the maximum overlap method
Definition: qs_mom_methods.F:240
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition: qs_rho_methods.F:15
qs_rho_methods::qs_rho_update_rho
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
Definition: qs_rho_methods.F:375
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_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_methods
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
Definition: qs_scf_methods.F:19
qs_scf_methods::eigensolver_simple
subroutine, public eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, level_shift, use_jacobi, jacobi_threshold)
...
Definition: qs_scf_methods.F:405
qs_scf_methods::eigensolver
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
Definition: qs_scf_methods.F:168
qs_scf_types
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
qs_scf_types::ot_diag_method_nr
integer, parameter, public ot_diag_method_nr
Definition: qs_scf_types.F:51
qs_scf_types::block_davidson_diag_method_nr
integer, parameter, public block_davidson_diag_method_nr
Definition: qs_scf_types.F:51
qs_scf_types::block_krylov_diag_method_nr
integer, parameter, public block_krylov_diag_method_nr
Definition: qs_scf_types.F:51
qs_scf_types::general_diag_method_nr
integer, parameter, public general_diag_method_nr
Definition: qs_scf_types.F:51
qs_wf_history_methods
Storage of past states of the qs_env. Methods to interpolate (or actually normally extrapolate) the n...
Definition: qs_wf_history_methods.F:21
qs_wf_history_methods::wfi_update
subroutine, public wfi_update(wf_history, qs_env, dt)
updates the snapshot buffer, taking a new snapshot
Definition: qs_wf_history_methods.F:968
scf_control_types
parameters that control an scf iteration
Definition: scf_control_types.F:17
util
All kind of helpful little routines.
Definition: util.F:14
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