(git:ccc2433)
almo_scf.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 for all ALMO-based SCF methods
10 !> 'RZK-warning' marks unresolved issues
11 !> \par History
12 !> 2011.05 created [Rustam Z Khaliullin]
13 !> \author Rustam Z Khaliullin
14 ! **************************************************************************************************
15 MODULE almo_scf
26  USE almo_scf_qs, ONLY: almo_dm_to_almo_ks,&
38  almo_scf_env_type,&
39  optimizer_options_type,&
41  USE atomic_kind_types, ONLY: atomic_kind_type
42  USE bibliography, ONLY: khaliullin2013,&
43  kolafa2004,&
44  kuhne2007,&
45  scheiber2018,&
46  staub2019,&
47  cite_reference
49  USE cp_control_types, ONLY: dft_control_type
50  USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
52  USE cp_fm_types, ONLY: cp_fm_type
55  cp_logger_type
56  USE dbcsr_api, ONLY: &
57  dbcsr_add, dbcsr_add_on_diag, dbcsr_binary_read, dbcsr_checksum, dbcsr_copy, dbcsr_create, &
58  dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
59  dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_init_random, &
60  dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
61  dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_nblkcols_total, &
62  dbcsr_nblkrows_total, dbcsr_p_type, dbcsr_release, dbcsr_reserve_block2d, dbcsr_scale, &
63  dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_work_create
64  USE domain_submatrix_methods, ONLY: init_submatrices,&
65  release_submatrices
66  USE input_constants, ONLY: &
74  USE input_section_types, ONLY: section_vals_type
75  USE iterate_matrix, ONLY: invert_hotelling,&
77  USE kinds, ONLY: default_path_length,&
78  dp
79  USE mathlib, ONLY: binomial
80  USE message_passing, ONLY: mp_comm_type,&
82  mp_para_env_type
84  molecule_type
86  molecular_scf_guess_env_type
87  USE particle_types, ONLY: particle_type
89  USE qs_environment_types, ONLY: get_qs_env,&
90  qs_environment_type
92  USE qs_kind_types, ONLY: qs_kind_type
93  USE qs_mo_types, ONLY: get_mo_set,&
94  mo_set_type
95  USE qs_rho_types, ONLY: qs_rho_get,&
96  qs_rho_type
98  USE qs_scf_types, ONLY: qs_scf_env_type
99 #include "./base/base_uses.f90"
100 
101  IMPLICIT NONE
102 
103  PRIVATE
104 
105  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'almo_scf'
106 
107  PUBLIC :: almo_entry_scf
108 
109  LOGICAL, PARAMETER :: debug_mode = .false.
110  LOGICAL, PARAMETER :: safe_mode = .false.
111 
112 CONTAINS
113 
114 ! **************************************************************************************************
115 !> \brief The entry point into ALMO SCF routines
116 !> \param qs_env pointer to the QS environment
117 !> \param calc_forces calculate forces?
118 !> \par History
119 !> 2011.05 created [Rustam Z Khaliullin]
120 !> \author Rustam Z Khaliullin
121 ! **************************************************************************************************
122  SUBROUTINE almo_entry_scf(qs_env, calc_forces)
123  TYPE(qs_environment_type), POINTER :: qs_env
124  LOGICAL, INTENT(IN) :: calc_forces
125 
126  CHARACTER(len=*), PARAMETER :: routinen = 'almo_entry_scf'
127 
128  INTEGER :: handle
129  TYPE(almo_scf_env_type), POINTER :: almo_scf_env
130 
131  CALL timeset(routinen, handle)
132 
133  CALL cite_reference(khaliullin2013)
134 
135  ! get a pointer to the almo environment
136  CALL get_qs_env(qs_env, almo_scf_env=almo_scf_env)
137 
138  ! initialize scf
139  CALL almo_scf_init(qs_env, almo_scf_env, calc_forces)
140 
141  ! create the initial guess for ALMOs
142  CALL almo_scf_initial_guess(qs_env, almo_scf_env)
143 
144  ! perform SCF for block diagonal ALMOs
145  CALL almo_scf_main(qs_env, almo_scf_env)
146 
147  ! allow electron delocalization
148  CALL almo_scf_delocalization(qs_env, almo_scf_env)
149 
150  ! construct NLMOs
151  CALL construct_nlmos(qs_env, almo_scf_env)
152 
153  ! electron correlation methods
154  !CALL almo_correlation_main(qs_env,almo_scf_env)
155 
156  ! do post scf processing
157  CALL almo_scf_post(qs_env, almo_scf_env)
158 
159  ! clean up the mess
160  CALL almo_scf_clean_up(almo_scf_env)
161 
162  CALL timestop(handle)
163 
164  END SUBROUTINE almo_entry_scf
165 
166 ! **************************************************************************************************
167 !> \brief Initialization of the almo_scf_env_type.
168 !> \param qs_env ...
169 !> \param almo_scf_env ...
170 !> \param calc_forces ...
171 !> \par History
172 !> 2011.05 created [Rustam Z Khaliullin]
173 !> 2018.09 smearing support [Ruben Staub]
174 !> \author Rustam Z Khaliullin
175 ! **************************************************************************************************
176  SUBROUTINE almo_scf_init(qs_env, almo_scf_env, calc_forces)
177  TYPE(qs_environment_type), POINTER :: qs_env
178  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
179  LOGICAL, INTENT(IN) :: calc_forces
180 
181  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_init'
182 
183  INTEGER :: ao, handle, i, iao, idomain, ispin, &
184  multip, naos, natoms, ndomains, nelec, &
185  nelec_a, nelec_b, nmols, nspins, &
186  unit_nr
187  TYPE(cp_logger_type), POINTER :: logger
188  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
189  TYPE(dft_control_type), POINTER :: dft_control
190  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
191  TYPE(section_vals_type), POINTER :: input
192 
193  CALL timeset(routinen, handle)
194 
195  ! define the output_unit
196  logger => cp_get_default_logger()
197  IF (logger%para_env%is_source()) THEN
198  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
199  ELSE
200  unit_nr = -1
201  END IF
202 
203  ! set optimizers' types
204  almo_scf_env%opt_block_diag_diis%optimizer_type = optimizer_diis
205  almo_scf_env%opt_block_diag_pcg%optimizer_type = optimizer_pcg
206  almo_scf_env%opt_xalmo_diis%optimizer_type = optimizer_diis
207  almo_scf_env%opt_xalmo_pcg%optimizer_type = optimizer_pcg
208  almo_scf_env%opt_xalmo_trustr%optimizer_type = optimizer_trustr
209  almo_scf_env%opt_nlmo_pcg%optimizer_type = optimizer_pcg
210  almo_scf_env%opt_block_diag_trustr%optimizer_type = optimizer_trustr
211  almo_scf_env%opt_xalmo_newton_pcg_solver%optimizer_type = optimizer_lin_eq_pcg
212 
213  ! get info from the qs_env
214  CALL get_qs_env(qs_env, &
215  nelectron_total=almo_scf_env%nelectrons_total, &
216  matrix_s=matrix_s, &
217  dft_control=dft_control, &
218  molecule_set=molecule_set, &
219  input=input, &
220  has_unit_metric=almo_scf_env%orthogonal_basis, &
221  para_env=almo_scf_env%para_env, &
222  blacs_env=almo_scf_env%blacs_env, &
223  nelectron_spin=almo_scf_env%nelectrons_spin)
224  CALL almo_scf_env%para_env%retain()
225  CALL almo_scf_env%blacs_env%retain()
226 
227  ! copy basic quantities
228  almo_scf_env%nspins = dft_control%nspins
229  almo_scf_env%natoms = dbcsr_nblkrows_total(matrix_s(1)%matrix)
230  almo_scf_env%nmolecules = SIZE(molecule_set)
231  CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=naos)
232  almo_scf_env%naos = naos
233 
234  !! retrieve smearing parameters, and check compatibility of methods requested
235  almo_scf_env%smear = dft_control%smear
236  IF (almo_scf_env%smear) THEN
237  CALL cite_reference(staub2019)
238  IF ((almo_scf_env%almo_update_algorithm .NE. almo_scf_diag) .OR. &
239  ((almo_scf_env%deloc_method .NE. almo_deloc_none) .AND. &
240  (almo_scf_env%xalmo_update_algorithm .NE. almo_scf_diag))) THEN
241  cpabort("ALMO smearing is currently implemented for DIAG algorithm only")
242  END IF
243  IF (qs_env%scf_control%smear%method .NE. smear_fermi_dirac) THEN
244  cpabort("Only Fermi-Dirac smearing is currently compatible with ALMO")
245  END IF
246  almo_scf_env%smear_e_temp = qs_env%scf_control%smear%electronic_temperature
247  IF ((almo_scf_env%mat_distr_aos .NE. almo_mat_distr_molecular) .OR. &
248  (almo_scf_env%domain_layout_mos .NE. almo_domain_layout_molecular)) THEN
249  cpabort("ALMO smearing was designed to work with molecular fragments only")
250  END IF
251  END IF
252 
253  ! convenient local varibales
254  nspins = almo_scf_env%nspins
255  nmols = almo_scf_env%nmolecules
256  natoms = almo_scf_env%natoms
257 
258  ! Define groups: either atomic or molecular
259  IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
260  almo_scf_env%ndomains = almo_scf_env%nmolecules
261  ELSE
262  almo_scf_env%ndomains = almo_scf_env%natoms
263  END IF
264 
265  ! allocate domain descriptors
266  ndomains = almo_scf_env%ndomains
267  ALLOCATE (almo_scf_env%domain_index_of_atom(natoms))
268  ALLOCATE (almo_scf_env%domain_index_of_ao(naos))
269  ALLOCATE (almo_scf_env%first_atom_of_domain(ndomains))
270  ALLOCATE (almo_scf_env%last_atom_of_domain(ndomains))
271  ALLOCATE (almo_scf_env%nbasis_of_domain(ndomains))
272  ALLOCATE (almo_scf_env%nocc_of_domain(ndomains, nspins)) !! with smearing, nb of available orbitals for occupation
273  ALLOCATE (almo_scf_env%real_ne_of_domain(ndomains, nspins)) !! with smearing, nb of fully-occupied orbitals
274  ALLOCATE (almo_scf_env%nvirt_full_of_domain(ndomains, nspins))
275  ALLOCATE (almo_scf_env%nvirt_of_domain(ndomains, nspins))
276  ALLOCATE (almo_scf_env%nvirt_disc_of_domain(ndomains, nspins))
277  ALLOCATE (almo_scf_env%mu_of_domain(ndomains, nspins))
278  ALLOCATE (almo_scf_env%cpu_of_domain(ndomains))
279  ALLOCATE (almo_scf_env%charge_of_domain(ndomains))
280  ALLOCATE (almo_scf_env%multiplicity_of_domain(ndomains))
281 
282  ! fill out domain descriptors and group descriptors
283  IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
284  ! get domain info from molecule_set
285  CALL get_molecule_set_info(molecule_set, &
286  atom_to_mol=almo_scf_env%domain_index_of_atom, &
287  mol_to_first_atom=almo_scf_env%first_atom_of_domain, &
288  mol_to_last_atom=almo_scf_env%last_atom_of_domain, &
289  mol_to_nelectrons=almo_scf_env%nocc_of_domain(1:ndomains, 1), &
290  mol_to_nbasis=almo_scf_env%nbasis_of_domain, &
291  mol_to_charge=almo_scf_env%charge_of_domain, &
292  mol_to_multiplicity=almo_scf_env%multiplicity_of_domain)
293  ! calculate number of alpha and beta occupied orbitals from
294  ! the number of electrons and multiplicity of each molecule
295  ! Na + Nb = Ne
296  ! Na - Nb = Mult - 1 (assume Na > Nb as we do not have more info from get_molecule_set_info)
297  DO idomain = 1, ndomains
298  nelec = almo_scf_env%nocc_of_domain(idomain, 1)
299  multip = almo_scf_env%multiplicity_of_domain(idomain)
300  nelec_a = (nelec + multip - 1)/2
301  nelec_b = nelec - nelec_a
302  !! Initializing an occupation-rescaling trick if smearing is on
303  IF (almo_scf_env%smear) THEN
304  IF (multip .GT. 1) THEN
305  cpwarn("BEWARE: Non singlet state detected, treating it as closed-shell")
306  END IF
307  !! Save real number of electrons of each spin, as it is required for Fermi-dirac smearing
308  !! BEWARE : Non singlet states are allowed but treated as closed-shell
309  almo_scf_env%real_ne_of_domain(idomain, :) = real(nelec, kind=dp)/2.0_dp
310  !! Add a number of added_mos equal to the number of atoms in domain
311  !! (since fragments were computed this way with smearing)
312  almo_scf_env%nocc_of_domain(idomain, :) = ceiling(almo_scf_env%real_ne_of_domain(idomain, :)) &
313  + (almo_scf_env%last_atom_of_domain(idomain) &
314  - almo_scf_env%first_atom_of_domain(idomain) + 1)
315  ELSE
316  almo_scf_env%nocc_of_domain(idomain, 1) = nelec_a
317  IF (nelec_a .NE. nelec_b) THEN
318  IF (nspins .EQ. 1) THEN
319  WRITE (*, *) "Domain ", idomain, " out of ", ndomains, ". Electrons = ", nelec
320  cpabort("odd e- -- use unrestricted methods")
321  END IF
322  almo_scf_env%nocc_of_domain(idomain, 2) = nelec_b
323  ! RZK-warning: open-shell procedures have not been tested yet
324  ! Stop the program now
325  cpabort("Unrestricted ALMO methods are NYI")
326  END IF
327  END IF
328  END DO
329  DO ispin = 1, nspins
330  ! take care of the full virtual subspace
331  almo_scf_env%nvirt_full_of_domain(:, ispin) = &
332  almo_scf_env%nbasis_of_domain(:) - &
333  almo_scf_env%nocc_of_domain(:, ispin)
334  ! and the truncated virtual subspace
335  SELECT CASE (almo_scf_env%deloc_truncate_virt)
336  CASE (virt_full)
337  almo_scf_env%nvirt_of_domain(:, ispin) = &
338  almo_scf_env%nvirt_full_of_domain(:, ispin)
339  almo_scf_env%nvirt_disc_of_domain(:, ispin) = 0
340  CASE (virt_number)
341  DO idomain = 1, ndomains
342  almo_scf_env%nvirt_of_domain(idomain, ispin) = &
343  min(almo_scf_env%deloc_virt_per_domain, &
344  almo_scf_env%nvirt_full_of_domain(idomain, ispin))
345  almo_scf_env%nvirt_disc_of_domain(idomain, ispin) = &
346  almo_scf_env%nvirt_full_of_domain(idomain, ispin) - &
347  almo_scf_env%nvirt_of_domain(idomain, ispin)
348  END DO
349  CASE (virt_occ_size)
350  DO idomain = 1, ndomains
351  almo_scf_env%nvirt_of_domain(idomain, ispin) = &
352  min(almo_scf_env%nocc_of_domain(idomain, ispin), &
353  almo_scf_env%nvirt_full_of_domain(idomain, ispin))
354  almo_scf_env%nvirt_disc_of_domain(idomain, ispin) = &
355  almo_scf_env%nvirt_full_of_domain(idomain, ispin) - &
356  almo_scf_env%nvirt_of_domain(idomain, ispin)
357  END DO
358  CASE DEFAULT
359  cpabort("illegal method for virtual space truncation")
360  END SELECT
361  END DO ! spin
362  ELSE ! domains are atomic
363  ! RZK-warning do the same for atomic domains/groups
364  almo_scf_env%domain_index_of_atom(1:natoms) = (/(i, i=1, natoms)/)
365  END IF
366 
367  ao = 1
368  DO idomain = 1, ndomains
369  DO iao = 1, almo_scf_env%nbasis_of_domain(idomain)
370  almo_scf_env%domain_index_of_ao(ao) = idomain
371  ao = ao + 1
372  END DO
373  END DO
374 
375  almo_scf_env%mu_of_domain(:, :) = almo_scf_env%mu
376 
377  ! build domain (i.e. layout) indices for distribution blocks
378  ! ao blocks
379  IF (almo_scf_env%mat_distr_aos == almo_mat_distr_atomic) THEN
380  ALLOCATE (almo_scf_env%domain_index_of_ao_block(natoms))
381  almo_scf_env%domain_index_of_ao_block(:) = &
382  almo_scf_env%domain_index_of_atom(:)
383  ELSE IF (almo_scf_env%mat_distr_aos == almo_mat_distr_molecular) THEN
384  ALLOCATE (almo_scf_env%domain_index_of_ao_block(nmols))
385  ! if distr blocks are molecular then domain layout is also molecular
386  almo_scf_env%domain_index_of_ao_block(:) = (/(i, i=1, nmols)/)
387  END IF
388  ! mo blocks
389  IF (almo_scf_env%mat_distr_mos == almo_mat_distr_atomic) THEN
390  ALLOCATE (almo_scf_env%domain_index_of_mo_block(natoms))
391  almo_scf_env%domain_index_of_mo_block(:) = &
392  almo_scf_env%domain_index_of_atom(:)
393  ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN
394  ALLOCATE (almo_scf_env%domain_index_of_mo_block(nmols))
395  ! if distr blocks are molecular then domain layout is also molecular
396  almo_scf_env%domain_index_of_mo_block(:) = (/(i, i=1, nmols)/)
397  END IF
398 
399  ! set all flags
400  !almo_scf_env%need_previous_ks=.FALSE.
401  !IF (almo_scf_env%deloc_method==almo_deloc_harris) THEN
402  almo_scf_env%need_previous_ks = .true.
403  !ENDIF
404 
405  !almo_scf_env%need_virtuals=.FALSE.
406  !almo_scf_env%need_orbital_energies=.FALSE.
407  !IF (almo_scf_env%almo_update_algorithm==almo_scf_diag) THEN
408  almo_scf_env%need_virtuals = .true.
409  almo_scf_env%need_orbital_energies = .true.
410  !ENDIF
411 
412  almo_scf_env%calc_forces = calc_forces
413  IF (calc_forces) THEN
414  CALL cite_reference(scheiber2018)
415  IF (almo_scf_env%deloc_method .EQ. almo_deloc_x .OR. &
416  almo_scf_env%deloc_method .EQ. almo_deloc_xalmo_x .OR. &
417  almo_scf_env%deloc_method .EQ. almo_deloc_xalmo_1diag) THEN
418  cpabort("Forces for perturbative methods are NYI. Change DELOCALIZE_METHOD")
419  END IF
420  ! switch to ASPC after a certain number of exact steps is done
421  IF (almo_scf_env%almo_history%istore .GT. (almo_scf_env%almo_history%nstore + 1)) THEN
422  IF (almo_scf_env%opt_block_diag_pcg%eps_error_early .GT. 0.0_dp) THEN
423  almo_scf_env%opt_block_diag_pcg%eps_error = almo_scf_env%opt_block_diag_pcg%eps_error_early
424  almo_scf_env%opt_block_diag_pcg%early_stopping_on = .true.
425  IF (unit_nr > 0) WRITE (*, *) "ALMO_OPTIMIZER_PCG: EPS_ERROR_EARLY is on"
426  END IF
427  IF (almo_scf_env%opt_block_diag_diis%eps_error_early .GT. 0.0_dp) THEN
428  almo_scf_env%opt_block_diag_diis%eps_error = almo_scf_env%opt_block_diag_diis%eps_error_early
429  almo_scf_env%opt_block_diag_diis%early_stopping_on = .true.
430  IF (unit_nr > 0) WRITE (*, *) "ALMO_OPTIMIZER_DIIS: EPS_ERROR_EARLY is on"
431  END IF
432  IF (almo_scf_env%opt_block_diag_pcg%max_iter_early .GT. 0) THEN
433  almo_scf_env%opt_block_diag_pcg%max_iter = almo_scf_env%opt_block_diag_pcg%max_iter_early
434  almo_scf_env%opt_block_diag_pcg%early_stopping_on = .true.
435  IF (unit_nr > 0) WRITE (*, *) "ALMO_OPTIMIZER_PCG: MAX_ITER_EARLY is on"
436  END IF
437  IF (almo_scf_env%opt_block_diag_diis%max_iter_early .GT. 0) THEN
438  almo_scf_env%opt_block_diag_diis%max_iter = almo_scf_env%opt_block_diag_diis%max_iter_early
439  almo_scf_env%opt_block_diag_diis%early_stopping_on = .true.
440  IF (unit_nr > 0) WRITE (*, *) "ALMO_OPTIMIZER_DIIS: MAX_ITER_EARLY is on"
441  END IF
442  ELSE
443  almo_scf_env%opt_block_diag_diis%early_stopping_on = .false.
444  almo_scf_env%opt_block_diag_pcg%early_stopping_on = .false.
445  END IF
446  IF (almo_scf_env%xalmo_history%istore .GT. (almo_scf_env%xalmo_history%nstore + 1)) THEN
447  IF (almo_scf_env%opt_xalmo_pcg%eps_error_early .GT. 0.0_dp) THEN
448  almo_scf_env%opt_xalmo_pcg%eps_error = almo_scf_env%opt_xalmo_pcg%eps_error_early
449  almo_scf_env%opt_xalmo_pcg%early_stopping_on = .true.
450  IF (unit_nr > 0) WRITE (*, *) "XALMO_OPTIMIZER_PCG: EPS_ERROR_EARLY is on"
451  END IF
452  IF (almo_scf_env%opt_xalmo_pcg%max_iter_early .GT. 0.0_dp) THEN
453  almo_scf_env%opt_xalmo_pcg%max_iter = almo_scf_env%opt_xalmo_pcg%max_iter_early
454  almo_scf_env%opt_xalmo_pcg%early_stopping_on = .true.
455  IF (unit_nr > 0) WRITE (*, *) "XALMO_OPTIMIZER_PCG: MAX_ITER_EARLY is on"
456  END IF
457  ELSE
458  almo_scf_env%opt_xalmo_pcg%early_stopping_on = .false.
459  END IF
460  END IF
461 
462  ! create all matrices
463  CALL almo_scf_env_create_matrices(almo_scf_env, matrix_s(1)%matrix)
464 
465  ! set up matrix S and all required functions of S
466  almo_scf_env%s_inv_done = .false.
467  almo_scf_env%s_sqrt_done = .false.
468  CALL almo_scf_init_ao_overlap(matrix_s(1)%matrix, almo_scf_env)
469 
470  ! create the quencher (imposes sparsity template)
471  CALL almo_scf_construct_quencher(qs_env, almo_scf_env)
472  CALL distribute_domains(almo_scf_env)
473 
474  ! FINISH setting job parameters here, print out job info
475  CALL almo_scf_print_job_info(almo_scf_env, unit_nr)
476 
477  ! allocate and init the domain preconditioner
478  ALLOCATE (almo_scf_env%domain_preconditioner(ndomains, nspins))
479  CALL init_submatrices(almo_scf_env%domain_preconditioner)
480 
481  ! allocate and init projected KS for domains
482  ALLOCATE (almo_scf_env%domain_ks_xx(ndomains, nspins))
483  CALL init_submatrices(almo_scf_env%domain_ks_xx)
484 
485  ! init ao-overlap subblocks
486  ALLOCATE (almo_scf_env%domain_s_inv(ndomains, nspins))
487  CALL init_submatrices(almo_scf_env%domain_s_inv)
488  ALLOCATE (almo_scf_env%domain_s_sqrt_inv(ndomains, nspins))
489  CALL init_submatrices(almo_scf_env%domain_s_sqrt_inv)
490  ALLOCATE (almo_scf_env%domain_s_sqrt(ndomains, nspins))
491  CALL init_submatrices(almo_scf_env%domain_s_sqrt)
492  ALLOCATE (almo_scf_env%domain_t(ndomains, nspins))
493  CALL init_submatrices(almo_scf_env%domain_t)
494  ALLOCATE (almo_scf_env%domain_err(ndomains, nspins))
495  CALL init_submatrices(almo_scf_env%domain_err)
496  ALLOCATE (almo_scf_env%domain_r_down_up(ndomains, nspins))
497  CALL init_submatrices(almo_scf_env%domain_r_down_up)
498 
499  ! initialization of the KS matrix
500  CALL init_almo_ks_matrix_via_qs(qs_env, &
501  almo_scf_env%matrix_ks, &
502  almo_scf_env%mat_distr_aos, &
503  almo_scf_env%eps_filter)
504  CALL construct_qs_mos(qs_env, almo_scf_env)
505 
506  CALL timestop(handle)
507 
508  END SUBROUTINE almo_scf_init
509 
510 ! **************************************************************************************************
511 !> \brief create the scf initial guess for ALMOs
512 !> \param qs_env ...
513 !> \param almo_scf_env ...
514 !> \par History
515 !> 2016.11 created [Rustam Z Khaliullin]
516 !> 2018.09 smearing support [Ruben Staub]
517 !> \author Rustam Z Khaliullin
518 ! **************************************************************************************************
519  SUBROUTINE almo_scf_initial_guess(qs_env, almo_scf_env)
520  TYPE(qs_environment_type), POINTER :: qs_env
521  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
522 
523  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_initial_guess'
524 
525  CHARACTER(LEN=default_path_length) :: file_name, project_name
526  INTEGER :: handle, iaspc, ispin, istore, naspc, &
527  nspins, unit_nr
528  INTEGER, DIMENSION(2) :: nelectron_spin
529  LOGICAL :: aspc_guess, has_unit_metric
530  REAL(kind=dp) :: alpha, cs_pos, energy, kts_sum
531  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
532  TYPE(cp_logger_type), POINTER :: logger
533  TYPE(dbcsr_distribution_type) :: dist
534  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
535  TYPE(dft_control_type), POINTER :: dft_control
536  TYPE(molecular_scf_guess_env_type), POINTER :: mscfg_env
537  TYPE(mp_para_env_type), POINTER :: para_env
538  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
539  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
540  TYPE(qs_rho_type), POINTER :: rho
541 
542  CALL timeset(routinen, handle)
543 
544  NULLIFY (rho, rho_ao)
545 
546  ! get a useful output_unit
547  logger => cp_get_default_logger()
548  IF (logger%para_env%is_source()) THEN
549  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
550  ELSE
551  unit_nr = -1
552  END IF
553 
554  ! get basic quantities from the qs_env
555  CALL get_qs_env(qs_env, &
556  dft_control=dft_control, &
557  matrix_s=matrix_s, &
558  atomic_kind_set=atomic_kind_set, &
559  qs_kind_set=qs_kind_set, &
560  particle_set=particle_set, &
561  has_unit_metric=has_unit_metric, &
562  para_env=para_env, &
563  nelectron_spin=nelectron_spin, &
564  mscfg_env=mscfg_env, &
565  rho=rho)
566 
567  CALL qs_rho_get(rho, rho_ao=rho_ao)
568  cpassert(ASSOCIATED(mscfg_env))
569 
570  ! initial guess on the first simulation step is determined by almo_scf_env%almo_scf_guess
571  ! the subsequent simulation steps are determined by extrapolation_order
572  ! if extrapolation order is zero then again almo_scf_env%almo_scf_guess is used
573  ! ... the number of stored history points will remain zero if extrapolation order is zero
574  IF (almo_scf_env%almo_history%istore == 0) THEN
575  aspc_guess = .false.
576  ELSE
577  aspc_guess = .true.
578  END IF
579 
580  nspins = almo_scf_env%nspins
581 
582  ! create an initial guess
583  IF (.NOT. aspc_guess) THEN
584 
585  SELECT CASE (almo_scf_env%almo_scf_guess)
586  CASE (molecular_guess)
587 
588  DO ispin = 1, nspins
589 
590  ! the calculations on "isolated" molecules has already been done
591  ! all we need to do is convert the MOs of molecules into
592  ! the ALMO matrix taking into account different distributions
593  CALL get_matrix_from_submatrices(mscfg_env, &
594  almo_scf_env%matrix_t_blk(ispin), ispin)
595  CALL dbcsr_filter(almo_scf_env%matrix_t_blk(ispin), &
596  almo_scf_env%eps_filter)
597 
598  END DO
599 
600  CASE (atomic_guess)
601 
602  IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical .OR. &
603  dft_control%qs_control%xtb) THEN
604  CALL calculate_mopac_dm(rho_ao, &
605  matrix_s(1)%matrix, has_unit_metric, &
606  dft_control, particle_set, atomic_kind_set, qs_kind_set, &
607  nspins, nelectron_spin, &
608  para_env)
609  ELSE
610  CALL calculate_atomic_block_dm(rho_ao, matrix_s(1)%matrix, atomic_kind_set, qs_kind_set, &
611  nspins, nelectron_spin, unit_nr, para_env)
612  END IF
613 
614  DO ispin = 1, nspins
615  ! copy the atomic-block dm into matrix_p_blk
616  CALL matrix_qs_to_almo(rho_ao(ispin)%matrix, &
617  almo_scf_env%matrix_p_blk(ispin), almo_scf_env%mat_distr_aos, &
618  .false.)
619  CALL dbcsr_filter(almo_scf_env%matrix_p_blk(ispin), &
620  almo_scf_env%eps_filter)
621  END DO ! ispin
622 
623  ! obtain orbitals from the density matrix
624  ! (the current version of ALMO SCF needs orbitals)
625  CALL almo_scf_p_blk_to_t_blk(almo_scf_env, ionic=.false.)
626 
627  CASE (restart_guess)
628 
629  project_name = logger%iter_info%project_name
630 
631  DO ispin = 1, nspins
632  WRITE (file_name, '(A,I0,A)') trim(project_name)//"_ALMO_SPIN_", ispin, "_RESTART.mo"
633  CALL dbcsr_get_info(almo_scf_env%matrix_t_blk(ispin), distribution=dist)
634  CALL dbcsr_binary_read(file_name, distribution=dist, matrix_new=almo_scf_env%matrix_t_blk(ispin))
635  cs_pos = dbcsr_checksum(almo_scf_env%matrix_t_blk(ispin), pos=.true.)
636  IF (unit_nr > 0) THEN
637  WRITE (unit_nr, '(T2,A,E20.8)') "Read restart ALMO "//trim(file_name)//" with checksum: ", cs_pos
638  END IF
639  END DO
640  END SELECT
641 
642  ELSE !aspc_guess
643 
644  CALL cite_reference(kolafa2004)
645  CALL cite_reference(kuhne2007)
646 
647  naspc = min(almo_scf_env%almo_history%istore, almo_scf_env%almo_history%nstore)
648  IF (unit_nr > 0) THEN
649  WRITE (unit_nr, fmt="(/,T2,A,/,/,T3,A,I0)") &
650  "Parameters for the always stable predictor-corrector (ASPC) method:", &
651  "ASPC order: ", naspc
652  END IF
653 
654  DO ispin = 1, nspins
655 
656  ! extrapolation
657  DO iaspc = 1, naspc
658  istore = mod(almo_scf_env%almo_history%istore - iaspc, almo_scf_env%almo_history%nstore) + 1
659  alpha = (-1.0_dp)**(iaspc + 1)*real(iaspc, kind=dp)* &
660  binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
661  IF (unit_nr > 0) THEN
662  WRITE (unit_nr, fmt="(T3,A2,I0,A4,F10.6)") &
663  "B(", iaspc, ") = ", alpha
664  END IF
665  IF (iaspc == 1) THEN
666  CALL dbcsr_copy(almo_scf_env%matrix_t_blk(ispin), &
667  almo_scf_env%almo_history%matrix_t(ispin), &
668  keep_sparsity=.true.)
669  CALL dbcsr_scale(almo_scf_env%matrix_t_blk(ispin), alpha)
670  ELSE
671  CALL dbcsr_multiply("N", "N", alpha, &
672  almo_scf_env%almo_history%matrix_p_up_down(ispin, istore), &
673  almo_scf_env%almo_history%matrix_t(ispin), &
674  1.0_dp, almo_scf_env%matrix_t_blk(ispin), &
675  retain_sparsity=.true.)
676  END IF
677  END DO !iaspc
678 
679  END DO !ispin
680 
681  END IF !aspc_guess?
682 
683  DO ispin = 1, nspins
684 
685  CALL orthogonalize_mos(ket=almo_scf_env%matrix_t_blk(ispin), &
686  overlap=almo_scf_env%matrix_sigma_blk(ispin), &
687  metric=almo_scf_env%matrix_s_blk(1), &
688  retain_locality=.true., &
689  only_normalize=.false., &
690  nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
691  eps_filter=almo_scf_env%eps_filter, &
692  order_lanczos=almo_scf_env%order_lanczos, &
693  eps_lanczos=almo_scf_env%eps_lanczos, &
694  max_iter_lanczos=almo_scf_env%max_iter_lanczos)
695 
696  !! Application of an occupation-rescaling trick for smearing, if requested
697  IF (almo_scf_env%smear) THEN
698  CALL almo_scf_t_rescaling(matrix_t=almo_scf_env%matrix_t_blk(ispin), &
699  mo_energies=almo_scf_env%mo_energies(:, ispin), &
700  mu_of_domain=almo_scf_env%mu_of_domain(:, ispin), &
701  real_ne_of_domain=almo_scf_env%real_ne_of_domain(:, ispin), &
702  spin_kts=almo_scf_env%kTS(ispin), &
703  smear_e_temp=almo_scf_env%smear_e_temp, &
704  ndomains=almo_scf_env%ndomains, &
705  nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin))
706  END IF
707 
708  CALL almo_scf_t_to_proj(t=almo_scf_env%matrix_t_blk(ispin), &
709  p=almo_scf_env%matrix_p(ispin), &
710  eps_filter=almo_scf_env%eps_filter, &
711  orthog_orbs=.false., &
712  nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
713  s=almo_scf_env%matrix_s(1), &
714  sigma=almo_scf_env%matrix_sigma(ispin), &
715  sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
716  use_guess=.false., &
717  smear=almo_scf_env%smear, &
718  algorithm=almo_scf_env%sigma_inv_algorithm, &
719  eps_lanczos=almo_scf_env%eps_lanczos, &
720  max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
721  inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
722  para_env=almo_scf_env%para_env, &
723  blacs_env=almo_scf_env%blacs_env)
724 
725  END DO
726 
727  ! compute dm from the projector(s)
728  IF (nspins == 1) THEN
729  CALL dbcsr_scale(almo_scf_env%matrix_p(1), 2.0_dp)
730  !! Rescaling electronic entropy contribution by spin_factor
731  IF (almo_scf_env%smear) THEN
732  almo_scf_env%kTS(1) = almo_scf_env%kTS(1)*2.0_dp
733  END IF
734  END IF
735 
736  IF (almo_scf_env%smear) THEN
737  kts_sum = sum(almo_scf_env%kTS)
738  ELSE
739  kts_sum = 0.0_dp
740  END IF
741 
742  CALL almo_dm_to_almo_ks(qs_env, &
743  almo_scf_env%matrix_p, &
744  almo_scf_env%matrix_ks, &
745  energy, &
746  almo_scf_env%eps_filter, &
747  almo_scf_env%mat_distr_aos, &
748  smear=almo_scf_env%smear, &
749  kts_sum=kts_sum)
750 
751  IF (unit_nr > 0) THEN
752  IF (almo_scf_env%almo_scf_guess .EQ. molecular_guess) THEN
753  WRITE (unit_nr, '(T2,A38,F40.10)') "Single-molecule energy:", &
754  sum(mscfg_env%energy_of_frag)
755  END IF
756  WRITE (unit_nr, '(T2,A38,F40.10)') "Energy of the initial guess:", energy
757  WRITE (unit_nr, '()')
758  END IF
759 
760  CALL timestop(handle)
761 
762  END SUBROUTINE almo_scf_initial_guess
763 
764 ! **************************************************************************************************
765 !> \brief store a history of matrices for later use in almo_scf_initial_guess
766 !> \param almo_scf_env ...
767 !> \par History
768 !> 2016.11 created [Rustam Z Khaliullin]
769 !> \author Rustam Khaliullin
770 ! **************************************************************************************************
771  SUBROUTINE almo_scf_store_extrapolation_data(almo_scf_env)
772  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
773 
774  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_store_extrapolation_data'
775 
776  INTEGER :: handle, ispin, istore, unit_nr
777  LOGICAL :: delocalization_uses_extrapolation
778  TYPE(cp_logger_type), POINTER :: logger
779  TYPE(dbcsr_type) :: matrix_no_tmp1, matrix_no_tmp2, &
780  matrix_no_tmp3, matrix_no_tmp4
781 
782  CALL timeset(routinen, handle)
783 
784  ! get a useful output_unit
785  logger => cp_get_default_logger()
786  IF (logger%para_env%is_source()) THEN
787  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
788  ELSE
789  unit_nr = -1
790  END IF
791 
792  IF (almo_scf_env%almo_history%nstore > 0) THEN
793 
794  almo_scf_env%almo_history%istore = almo_scf_env%almo_history%istore + 1
795 
796  DO ispin = 1, SIZE(almo_scf_env%matrix_t_blk)
797 
798  istore = mod(almo_scf_env%almo_history%istore - 1, almo_scf_env%almo_history%nstore) + 1
799 
800  IF (almo_scf_env%almo_history%istore == 1) THEN
801  CALL dbcsr_create(almo_scf_env%almo_history%matrix_t(ispin), &
802  template=almo_scf_env%matrix_t_blk(ispin), &
803  matrix_type=dbcsr_type_no_symmetry)
804  END IF
805  CALL dbcsr_copy(almo_scf_env%almo_history%matrix_t(ispin), &
806  almo_scf_env%matrix_t_blk(ispin))
807 
808  IF (almo_scf_env%almo_history%istore <= almo_scf_env%almo_history%nstore) THEN
809  CALL dbcsr_create(almo_scf_env%almo_history%matrix_p_up_down(ispin, istore), &
810  template=almo_scf_env%matrix_s(1), &
811  matrix_type=dbcsr_type_no_symmetry)
812  END IF
813 
814  CALL dbcsr_create(matrix_no_tmp1, template=almo_scf_env%matrix_t_blk(ispin), &
815  matrix_type=dbcsr_type_no_symmetry)
816  CALL dbcsr_create(matrix_no_tmp2, template=almo_scf_env%matrix_t_blk(ispin), &
817  matrix_type=dbcsr_type_no_symmetry)
818 
819  ! compute contra-covariant density matrix
820  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s(1), &
821  almo_scf_env%matrix_t_blk(ispin), &
822  0.0_dp, matrix_no_tmp1, &
823  filter_eps=almo_scf_env%eps_filter)
824  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_no_tmp1, &
825  almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
826  0.0_dp, matrix_no_tmp2, &
827  filter_eps=almo_scf_env%eps_filter)
828  CALL dbcsr_multiply("N", "T", 1.0_dp, &
829  almo_scf_env%matrix_t_blk(ispin), &
830  matrix_no_tmp2, &
831  0.0_dp, almo_scf_env%almo_history%matrix_p_up_down(ispin, istore), &
832  filter_eps=almo_scf_env%eps_filter)
833 
834  CALL dbcsr_release(matrix_no_tmp1)
835  CALL dbcsr_release(matrix_no_tmp2)
836 
837  END DO
838 
839  END IF
840 
841  ! exrapolate xalmos?
842  delocalization_uses_extrapolation = &
843  .NOT. ((almo_scf_env%deloc_method .EQ. almo_deloc_none) .OR. &
844  (almo_scf_env%deloc_method .EQ. almo_deloc_xalmo_1diag))
845  IF (almo_scf_env%xalmo_history%nstore > 0 .AND. &
846  delocalization_uses_extrapolation) THEN
847 
848  almo_scf_env%xalmo_history%istore = almo_scf_env%xalmo_history%istore + 1
849 
850  DO ispin = 1, SIZE(almo_scf_env%matrix_t)
851 
852  istore = mod(almo_scf_env%xalmo_history%istore - 1, almo_scf_env%xalmo_history%nstore) + 1
853 
854  IF (almo_scf_env%xalmo_history%istore == 1) THEN
855  CALL dbcsr_create(almo_scf_env%xalmo_history%matrix_t(ispin), &
856  template=almo_scf_env%matrix_t(ispin), &
857  matrix_type=dbcsr_type_no_symmetry)
858  END IF
859  CALL dbcsr_copy(almo_scf_env%xalmo_history%matrix_t(ispin), &
860  almo_scf_env%matrix_t(ispin))
861 
862  IF (almo_scf_env%xalmo_history%istore <= almo_scf_env%xalmo_history%nstore) THEN
863  !CALL dbcsr_init(almo_scf_env%xalmo_history%matrix_x(ispin, istore))
864  !CALL dbcsr_create(almo_scf_env%xalmo_history%matrix_x(ispin, istore), &
865  ! template=almo_scf_env%matrix_t(ispin), &
866  ! matrix_type=dbcsr_type_no_symmetry)
867  CALL dbcsr_create(almo_scf_env%xalmo_history%matrix_p_up_down(ispin, istore), &
868  template=almo_scf_env%matrix_s(1), &
869  matrix_type=dbcsr_type_no_symmetry)
870  END IF
871 
872  CALL dbcsr_create(matrix_no_tmp3, template=almo_scf_env%matrix_t(ispin), &
873  matrix_type=dbcsr_type_no_symmetry)
874  CALL dbcsr_create(matrix_no_tmp4, template=almo_scf_env%matrix_t(ispin), &
875  matrix_type=dbcsr_type_no_symmetry)
876 
877  ! compute contra-covariant density matrix
878  CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s(1), &
879  almo_scf_env%matrix_t(ispin), &
880  0.0_dp, matrix_no_tmp3, &
881  filter_eps=almo_scf_env%eps_filter)
882  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_no_tmp3, &
883  almo_scf_env%matrix_sigma_inv(ispin), &
884  0.0_dp, matrix_no_tmp4, &
885  filter_eps=almo_scf_env%eps_filter)
886  CALL dbcsr_multiply("N", "T", 1.0_dp, &
887  almo_scf_env%matrix_t(ispin), &
888  matrix_no_tmp4, &
889  0.0_dp, almo_scf_env%xalmo_history%matrix_p_up_down(ispin, istore), &
890  filter_eps=almo_scf_env%eps_filter)
891 
892  ! store the difference between t and t0
893  !CALL dbcsr_copy(almo_scf_env%xalmo_history%matrix_x(ispin, istore),&
894  ! almo_scf_env%matrix_t(ispin))
895  !CALL dbcsr_add(almo_scf_env%xalmo_history%matrix_x(ispin, istore),&
896  ! almo_scf_env%matrix_t_blk(ispin),1.0_dp,-1.0_dp)
897 
898  CALL dbcsr_release(matrix_no_tmp3)
899  CALL dbcsr_release(matrix_no_tmp4)
900 
901  END DO
902 
903  END IF
904 
905  CALL timestop(handle)
906 
907  END SUBROUTINE almo_scf_store_extrapolation_data
908 
909 ! **************************************************************************************************
910 !> \brief Prints out a short summary about the ALMO SCF job
911 !> \param almo_scf_env ...
912 !> \param unit_nr ...
913 !> \par History
914 !> 2011.10 created [Rustam Z Khaliullin]
915 !> \author Rustam Z Khaliullin
916 ! **************************************************************************************************
917  SUBROUTINE almo_scf_print_job_info(almo_scf_env, unit_nr)
918 
919  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
920  INTEGER, INTENT(IN) :: unit_nr
921 
922  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_print_job_info'
923 
924  CHARACTER(len=13) :: neig_string
925  CHARACTER(len=33) :: deloc_method_string
926  INTEGER :: handle, idomain, index1_prev, sum_temp
927  INTEGER, ALLOCATABLE, DIMENSION(:) :: nneighbors
928 
929  CALL timeset(routinen, handle)
930 
931  IF (unit_nr > 0) THEN
932  WRITE (unit_nr, '()')
933  WRITE (unit_nr, '(T2,A,A,A)') repeat("-", 32), " ALMO SETTINGS ", repeat("-", 32)
934 
935  WRITE (unit_nr, '(T2,A,T48,E33.3)') "eps_filter:", almo_scf_env%eps_filter
936 
937  IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_skip) THEN
938  WRITE (unit_nr, '(T2,A)') "skip optimization of block-diagonal ALMOs"
939  ELSE
940  WRITE (unit_nr, '(T2,A)') "optimization of block-diagonal ALMOs:"
941  SELECT CASE (almo_scf_env%almo_update_algorithm)
942  CASE (almo_scf_diag)
943  ! the DIIS algorith is the only choice for the diagonlaization-based algorithm
944  CALL print_optimizer_options(almo_scf_env%opt_block_diag_diis, unit_nr)
945  CASE (almo_scf_pcg)
946  ! print out PCG options
947  CALL print_optimizer_options(almo_scf_env%opt_block_diag_pcg, unit_nr)
948  CASE (almo_scf_trustr)
949  ! print out TRUST REGION options
950  CALL print_optimizer_options(almo_scf_env%opt_block_diag_trustr, unit_nr)
951  END SELECT
952  END IF
953 
954  SELECT CASE (almo_scf_env%deloc_method)
955  CASE (almo_deloc_none)
956  deloc_method_string = "NONE"
957  CASE (almo_deloc_x)
958  deloc_method_string = "FULL_X"
959  CASE (almo_deloc_scf)
960  deloc_method_string = "FULL_SCF"
961  CASE (almo_deloc_x_then_scf)
962  deloc_method_string = "FULL_X_THEN_SCF"
964  deloc_method_string = "XALMO_1DIAG"
965  CASE (almo_deloc_xalmo_x)
966  deloc_method_string = "XALMO_X"
967  CASE (almo_deloc_xalmo_scf)
968  deloc_method_string = "XALMO_SCF"
969  END SELECT
970  WRITE (unit_nr, '(T2,A,T48,A33)') "delocalization:", trim(deloc_method_string)
971 
972  IF (almo_scf_env%deloc_method .NE. almo_deloc_none) THEN
973 
974  SELECT CASE (almo_scf_env%deloc_method)
976  WRITE (unit_nr, '(T2,A,T48,A33)') "delocalization cutoff radius:", &
977  "infinite"
978  deloc_method_string = "FULL_X_THEN_SCF"
980  WRITE (unit_nr, '(T2,A,T48,F33.5)') "XALMO cutoff radius:", &
981  almo_scf_env%quencher_r0_factor
982  END SELECT
983 
984  IF (almo_scf_env%deloc_method .EQ. almo_deloc_xalmo_1diag) THEN
985  ! print nothing because no actual optimization is done
986  ELSE
987  WRITE (unit_nr, '(T2,A)') "optimization of extended orbitals:"
988  SELECT CASE (almo_scf_env%xalmo_update_algorithm)
989  CASE (almo_scf_diag)
990  CALL print_optimizer_options(almo_scf_env%opt_xalmo_diis, unit_nr)
991  CASE (almo_scf_trustr)
992  CALL print_optimizer_options(almo_scf_env%opt_xalmo_trustr, unit_nr)
993  CASE (almo_scf_pcg)
994  CALL print_optimizer_options(almo_scf_env%opt_xalmo_pcg, unit_nr)
995  END SELECT
996  END IF
997 
998  END IF
999 
1000  !SELECT CASE(almo_scf_env%domain_layout_mos)
1001  !CASE(almo_domain_layout_orbital)
1002  ! WRITE(unit_nr,'(T2,A,T48,A33)') "Delocalization constraints","ORBITAL"
1003  !CASE(almo_domain_layout_atomic)
1004  ! WRITE(unit_nr,'(T2,A,T48,A33)') "Delocalization constraints","ATOMIC"
1005  !CASE(almo_domain_layout_molecular)
1006  ! WRITE(unit_nr,'(T2,A,T48,A33)') "Delocalization constraints","MOLECULAR"
1007  !END SELECT
1008 
1009  !SELECT CASE(almo_scf_env%domain_layout_aos)
1010  !CASE(almo_domain_layout_atomic)
1011  ! WRITE(unit_nr,'(T2,A,T48,A33)') "Basis function domains","ATOMIC"
1012  !CASE(almo_domain_layout_molecular)
1013  ! WRITE(unit_nr,'(T2,A,T48,A33)') "Basis function domains","MOLECULAR"
1014  !END SELECT
1015 
1016  !SELECT CASE(almo_scf_env%mat_distr_aos)
1017  !CASE(almo_mat_distr_atomic)
1018  ! WRITE(unit_nr,'(T2,A,T48,A33)') "Parallel distribution for AOs","ATOMIC"
1019  !CASE(almo_mat_distr_molecular)
1020  ! WRITE(unit_nr,'(T2,A,T48,A33)') "Parallel distribution for AOs","MOLECULAR"
1021  !END SELECT
1022 
1023  !SELECT CASE(almo_scf_env%mat_distr_mos)
1024  !CASE(almo_mat_distr_atomic)
1025  ! WRITE(unit_nr,'(T2,A,T48,A33)') "Parallel distribution for MOs","ATOMIC"
1026  !CASE(almo_mat_distr_molecular)
1027  ! WRITE(unit_nr,'(T2,A,T48,A33)') "Parallel distribution for MOs","MOLECULAR"
1028  !END SELECT
1029 
1030  ! print fragment's statistics
1031  WRITE (unit_nr, '(T2,A)') repeat("-", 79)
1032  WRITE (unit_nr, '(T2,A,T48,I33)') "Total fragments:", &
1033  almo_scf_env%ndomains
1034 
1035  sum_temp = sum(almo_scf_env%nbasis_of_domain(:))
1036  WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1037  "Basis set size per fragment (min, av, max, total):", &
1038  minval(almo_scf_env%nbasis_of_domain(:)), &
1039  (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1040  maxval(almo_scf_env%nbasis_of_domain(:)), &
1041  sum_temp
1042  !WRITE (unit_nr, '(T2,I13,F13.3,I13,I13)') &
1043  ! MINVAL(almo_scf_env%nbasis_of_domain(:)), &
1044  ! (1.0_dp*sum_temp) / almo_scf_env%ndomains, &
1045  ! MAXVAL(almo_scf_env%nbasis_of_domain(:)), &
1046  ! sum_temp
1047 
1048  sum_temp = sum(almo_scf_env%nocc_of_domain(:, :))
1049  WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1050  "Occupied MOs per fragment (min, av, max, total):", &
1051  minval(sum(almo_scf_env%nocc_of_domain, dim=2)), &
1052  (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1053  maxval(sum(almo_scf_env%nocc_of_domain, dim=2)), &
1054  sum_temp
1055  !WRITE (unit_nr, '(T2,I13,F13.3,I13,I13)') &
1056  ! MINVAL( SUM(almo_scf_env%nocc_of_domain, DIM=2) ), &
1057  ! (1.0_dp*sum_temp) / almo_scf_env%ndomains, &
1058  ! MAXVAL( SUM(almo_scf_env%nocc_of_domain, DIM=2) ), &
1059  ! sum_temp
1060 
1061  sum_temp = sum(almo_scf_env%nvirt_of_domain(:, :))
1062  WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1063  "Virtual MOs per fragment (min, av, max, total):", &
1064  minval(sum(almo_scf_env%nvirt_of_domain, dim=2)), &
1065  (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1066  maxval(sum(almo_scf_env%nvirt_of_domain, dim=2)), &
1067  sum_temp
1068  !WRITE (unit_nr, '(T2,I13,F13.3,I13,I13)') &
1069  ! MINVAL( SUM(almo_scf_env%nvirt_of_domain, DIM=2) ), &
1070  ! (1.0_dp*sum_temp) / almo_scf_env%ndomains, &
1071  ! MAXVAL( SUM(almo_scf_env%nvirt_of_domain, DIM=2) ), &
1072  ! sum_temp
1073 
1074  sum_temp = sum(almo_scf_env%charge_of_domain(:))
1075  WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1076  "Charges per fragment (min, av, max, total):", &
1077  minval(almo_scf_env%charge_of_domain(:)), &
1078  (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1079  maxval(almo_scf_env%charge_of_domain(:)), &
1080  sum_temp
1081  !WRITE (unit_nr, '(T2,I13,F13.3,I13,I13)') &
1082  ! MINVAL(almo_scf_env%charge_of_domain(:)), &
1083  ! (1.0_dp*sum_temp) / almo_scf_env%ndomains, &
1084  ! MAXVAL(almo_scf_env%charge_of_domain(:)), &
1085  ! sum_temp
1086 
1087  ! compute the number of neighbors of each fragment
1088  ALLOCATE (nneighbors(almo_scf_env%ndomains))
1089 
1090  DO idomain = 1, almo_scf_env%ndomains
1091 
1092  IF (idomain .EQ. 1) THEN
1093  index1_prev = 1
1094  ELSE
1095  index1_prev = almo_scf_env%domain_map(1)%index1(idomain - 1)
1096  END IF
1097 
1098  SELECT CASE (almo_scf_env%deloc_method)
1099  CASE (almo_deloc_none)
1100  nneighbors(idomain) = 0
1102  nneighbors(idomain) = almo_scf_env%ndomains - 1 ! minus self
1104  nneighbors(idomain) = almo_scf_env%domain_map(1)%index1(idomain) - index1_prev - 1 ! minus self
1105  CASE DEFAULT
1106  nneighbors(idomain) = -1
1107  END SELECT
1108 
1109  END DO ! cycle over domains
1110 
1111  sum_temp = sum(nneighbors(:))
1112  WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1113  "Deloc. neighbors of fragment (min, av, max, total):", &
1114  minval(nneighbors(:)), &
1115  (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1116  maxval(nneighbors(:)), &
1117  sum_temp
1118 
1119  WRITE (unit_nr, '(T2,A)') repeat("-", 79)
1120  WRITE (unit_nr, '()')
1121 
1122  IF (almo_scf_env%ndomains .LE. 64) THEN
1123 
1124  ! print fragment info
1125  WRITE (unit_nr, '(T2,A10,A13,A13,A13,A13,A13)') &
1126  "Fragment", "Basis Set", "Occupied", "Virtual", "Charge", "Deloc Neig" !,"Discarded Virt"
1127  WRITE (unit_nr, '(T2,A)') repeat("-", 79)
1128  DO idomain = 1, almo_scf_env%ndomains
1129 
1130  SELECT CASE (almo_scf_env%deloc_method)
1131  CASE (almo_deloc_none)
1132  neig_string = "NONE"
1134  neig_string = "ALL"
1136  WRITE (neig_string, '(I13)') nneighbors(idomain)
1137  CASE DEFAULT
1138  neig_string = "N/A"
1139  END SELECT
1140 
1141  WRITE (unit_nr, '(T2,I10,I13,I13,I13,I13,A13)') &
1142  idomain, almo_scf_env%nbasis_of_domain(idomain), &
1143  sum(almo_scf_env%nocc_of_domain(idomain, :)), &
1144  sum(almo_scf_env%nvirt_of_domain(idomain, :)), &
1145  !SUM(almo_scf_env%nvirt_disc_of_domain(idomain,:)),&
1146  almo_scf_env%charge_of_domain(idomain), &
1147  adjustr(trim(neig_string))
1148 
1149  END DO ! cycle over domains
1150 
1151  SELECT CASE (almo_scf_env%deloc_method)
1153 
1154  WRITE (unit_nr, '(T2,A)') repeat("-", 79)
1155 
1156  ! print fragment neighbors
1157  WRITE (unit_nr, '(T2,A78)') &
1158  "Neighbor lists (including self)"
1159  WRITE (unit_nr, '(T2,A)') repeat("-", 79)
1160  DO idomain = 1, almo_scf_env%ndomains
1161 
1162  IF (idomain .EQ. 1) THEN
1163  index1_prev = 1
1164  ELSE
1165  index1_prev = almo_scf_env%domain_map(1)%index1(idomain - 1)
1166  END IF
1167 
1168  WRITE (unit_nr, '(T2,I10,":")') idomain
1169  WRITE (unit_nr, '(T12,11I6)') &
1170  almo_scf_env%domain_map(1)%pairs &
1171  (index1_prev:almo_scf_env%domain_map(1)%index1(idomain) - 1, 1) ! includes self
1172 
1173  END DO ! cycle over domains
1174 
1175  END SELECT
1176 
1177  ELSE ! too big to print details for each fragment
1178 
1179  WRITE (unit_nr, '(T2,A)') "The system is too big to print details for each fragment."
1180 
1181  END IF ! how many fragments?
1182 
1183  WRITE (unit_nr, '(T2,A)') repeat("-", 79)
1184 
1185  WRITE (unit_nr, '()')
1186 
1187  DEALLOCATE (nneighbors)
1188 
1189  END IF ! unit_nr > 0
1190 
1191  CALL timestop(handle)
1192 
1193  END SUBROUTINE almo_scf_print_job_info
1194 
1195 ! **************************************************************************************************
1196 !> \brief Initializes the ALMO SCF copy of the AO overlap matrix
1197 !> and all necessary functions (sqrt, inverse...)
1198 !> \param matrix_s ...
1199 !> \param almo_scf_env ...
1200 !> \par History
1201 !> 2011.06 created [Rustam Z Khaliullin]
1202 !> \author Rustam Z Khaliullin
1203 ! **************************************************************************************************
1204  SUBROUTINE almo_scf_init_ao_overlap(matrix_s, almo_scf_env)
1205  TYPE(dbcsr_type), INTENT(IN) :: matrix_s
1206  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1207 
1208  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_init_ao_overlap'
1209 
1210  INTEGER :: handle, unit_nr
1211  TYPE(cp_logger_type), POINTER :: logger
1212 
1213  CALL timeset(routinen, handle)
1214 
1215  ! get a useful output_unit
1216  logger => cp_get_default_logger()
1217  IF (logger%para_env%is_source()) THEN
1218  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1219  ELSE
1220  unit_nr = -1
1221  END IF
1222 
1223  ! make almo copy of S
1224  ! also copy S to S_blk (i.e. to S with the domain structure imposed)
1225  IF (almo_scf_env%orthogonal_basis) THEN
1226  CALL dbcsr_set(almo_scf_env%matrix_s(1), 0.0_dp)
1227  CALL dbcsr_add_on_diag(almo_scf_env%matrix_s(1), 1.0_dp)
1228  CALL dbcsr_set(almo_scf_env%matrix_s_blk(1), 0.0_dp)
1229  CALL dbcsr_add_on_diag(almo_scf_env%matrix_s_blk(1), 1.0_dp)
1230  ELSE
1231  CALL matrix_qs_to_almo(matrix_s, almo_scf_env%matrix_s(1), &
1232  almo_scf_env%mat_distr_aos, .false.)
1233  CALL dbcsr_copy(almo_scf_env%matrix_s_blk(1), &
1234  almo_scf_env%matrix_s(1), keep_sparsity=.true.)
1235  END IF
1236 
1237  CALL dbcsr_filter(almo_scf_env%matrix_s(1), almo_scf_env%eps_filter)
1238  CALL dbcsr_filter(almo_scf_env%matrix_s_blk(1), almo_scf_env%eps_filter)
1239 
1240  IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_diag) THEN
1241  CALL matrix_sqrt_newton_schulz(almo_scf_env%matrix_s_blk_sqrt(1), &
1242  almo_scf_env%matrix_s_blk_sqrt_inv(1), &
1243  almo_scf_env%matrix_s_blk(1), &
1244  threshold=almo_scf_env%eps_filter, &
1245  order=almo_scf_env%order_lanczos, &
1246  !order=0, &
1247  eps_lanczos=almo_scf_env%eps_lanczos, &
1248  max_iter_lanczos=almo_scf_env%max_iter_lanczos)
1249  ELSE IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_dm_sign) THEN
1250  CALL invert_hotelling(almo_scf_env%matrix_s_blk_inv(1), &
1251  almo_scf_env%matrix_s_blk(1), &
1252  threshold=almo_scf_env%eps_filter, &
1253  filter_eps=almo_scf_env%eps_filter)
1254  END IF
1255 
1256  CALL timestop(handle)
1257 
1258  END SUBROUTINE almo_scf_init_ao_overlap
1259 
1260 ! **************************************************************************************************
1261 !> \brief Selects the subroutine for the optimization of block-daigonal ALMOs.
1262 !> Keep it short and clean.
1263 !> \param qs_env ...
1264 !> \param almo_scf_env ...
1265 !> \par History
1266 !> 2011.11 created [Rustam Z Khaliullin]
1267 !> \author Rustam Z Khaliullin
1268 ! **************************************************************************************************
1269  SUBROUTINE almo_scf_main(qs_env, almo_scf_env)
1270  TYPE(qs_environment_type), POINTER :: qs_env
1271  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1272 
1273  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_main'
1274 
1275  INTEGER :: handle, ispin, unit_nr
1276  TYPE(cp_logger_type), POINTER :: logger
1277 
1278  CALL timeset(routinen, handle)
1279 
1280  ! get a useful output_unit
1281  logger => cp_get_default_logger()
1282  IF (logger%para_env%is_source()) THEN
1283  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1284  ELSE
1285  unit_nr = -1
1286  END IF
1287 
1288  SELECT CASE (almo_scf_env%almo_update_algorithm)
1290 
1291  SELECT CASE (almo_scf_env%almo_update_algorithm)
1292  CASE (almo_scf_pcg)
1293 
1294  ! ALMO PCG optimizer as a special case of XALMO PCG
1295  CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1296  almo_scf_env=almo_scf_env, &
1297  optimizer=almo_scf_env%opt_block_diag_pcg, &
1298  quench_t=almo_scf_env%quench_t_blk, &
1299  matrix_t_in=almo_scf_env%matrix_t_blk, &
1300  matrix_t_out=almo_scf_env%matrix_t_blk, &
1301  assume_t0_q0x=.false., &
1302  perturbation_only=.false., &
1303  special_case=xalmo_case_block_diag)
1304 
1305  CASE (almo_scf_trustr)
1306 
1307  CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
1308  almo_scf_env=almo_scf_env, &
1309  optimizer=almo_scf_env%opt_block_diag_trustr, &
1310  quench_t=almo_scf_env%quench_t_blk, &
1311  matrix_t_in=almo_scf_env%matrix_t_blk, &
1312  matrix_t_out=almo_scf_env%matrix_t_blk, &
1313  perturbation_only=.false., &
1314  special_case=xalmo_case_block_diag)
1315 
1316  END SELECT
1317 
1318  DO ispin = 1, almo_scf_env%nspins
1319  CALL orthogonalize_mos(ket=almo_scf_env%matrix_t_blk(ispin), &
1320  overlap=almo_scf_env%matrix_sigma_blk(ispin), &
1321  metric=almo_scf_env%matrix_s_blk(1), &
1322  retain_locality=.true., &
1323  only_normalize=.false., &
1324  nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1325  eps_filter=almo_scf_env%eps_filter, &
1326  order_lanczos=almo_scf_env%order_lanczos, &
1327  eps_lanczos=almo_scf_env%eps_lanczos, &
1328  max_iter_lanczos=almo_scf_env%max_iter_lanczos)
1329  END DO
1330 
1331  CASE (almo_scf_diag)
1332 
1333  ! mixing/DIIS optimizer
1334  CALL almo_scf_block_diagonal(qs_env, almo_scf_env, &
1335  almo_scf_env%opt_block_diag_diis)
1336 
1337  END SELECT
1338 
1339  ! we might need a copy of the converged KS and sigma_inv
1340  DO ispin = 1, almo_scf_env%nspins
1341  CALL dbcsr_copy(almo_scf_env%matrix_ks_0deloc(ispin), &
1342  almo_scf_env%matrix_ks(ispin))
1343  CALL dbcsr_copy(almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
1344  almo_scf_env%matrix_sigma_inv(ispin))
1345  END DO
1346 
1347  CALL timestop(handle)
1348 
1349  END SUBROUTINE almo_scf_main
1350 
1351 ! **************************************************************************************************
1352 !> \brief selects various post scf routines
1353 !> \param qs_env ...
1354 !> \param almo_scf_env ...
1355 !> \par History
1356 !> 2011.06 created [Rustam Z Khaliullin]
1357 !> \author Rustam Z Khaliullin
1358 ! **************************************************************************************************
1359  SUBROUTINE almo_scf_delocalization(qs_env, almo_scf_env)
1360 
1361  TYPE(qs_environment_type), POINTER :: qs_env
1362  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1363 
1364  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_delocalization'
1365 
1366  INTEGER :: col, handle, hold, iblock_col, &
1367  iblock_row, ispin, mynode, &
1368  nblkcols_tot, nblkrows_tot, row, &
1369  unit_nr
1370  LOGICAL :: almo_experimental, tr
1371  REAL(kind=dp), DIMENSION(:, :), POINTER :: p_new_block
1372  TYPE(cp_logger_type), POINTER :: logger
1373  TYPE(dbcsr_distribution_type) :: dist
1374  TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: no_quench
1375  TYPE(optimizer_options_type) :: arbitrary_optimizer
1376 
1377  CALL timeset(routinen, handle)
1378 
1379  ! get a useful output_unit
1380  logger => cp_get_default_logger()
1381  IF (logger%para_env%is_source()) THEN
1382  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1383  ELSE
1384  unit_nr = -1
1385  END IF
1386 
1387  ! create a local optimizer that handles XALMO DIIS
1388  ! the options of this optimizer are arbitrary because
1389  ! XALMO DIIS SCF does not converge for yet unknown reasons and
1390  ! currently used in the code to get perturbative estimates only
1391  arbitrary_optimizer%optimizer_type = optimizer_diis
1392  arbitrary_optimizer%max_iter = 3
1393  arbitrary_optimizer%eps_error = 1.0e-6_dp
1394  arbitrary_optimizer%ndiis = 2
1395 
1396  SELECT CASE (almo_scf_env%deloc_method)
1398 
1399  ! RZK-warning hack into the quenched routine:
1400  ! create a quench matrix with all-all-all blocks 1.0
1401  ! it is a waste of memory but since matrices are distributed
1402  ! we can tolerate it for now
1403  ALLOCATE (no_quench(almo_scf_env%nspins))
1404  CALL dbcsr_create(no_quench(1), &
1405  template=almo_scf_env%matrix_t(1), &
1406  matrix_type=dbcsr_type_no_symmetry)
1407  CALL dbcsr_get_info(no_quench(1), distribution=dist)
1408  CALL dbcsr_distribution_get(dist, mynode=mynode)
1409  CALL dbcsr_work_create(no_quench(1), &
1410  work_mutable=.true.)
1411  nblkrows_tot = dbcsr_nblkrows_total(no_quench(1))
1412  nblkcols_tot = dbcsr_nblkcols_total(no_quench(1))
1413  ! RZK-warning: is it a quadratic-scaling routine?
1414  ! As a matter of fact it is! But this block treats
1415  ! fully delocalized MOs. So it is unavoidable.
1416  ! C'est la vie
1417  DO row = 1, nblkrows_tot
1418  DO col = 1, nblkcols_tot
1419  tr = .false.
1420  iblock_row = row
1421  iblock_col = col
1422  CALL dbcsr_get_stored_coordinates(no_quench(1), &
1423  iblock_row, iblock_col, hold)
1424  IF (hold .EQ. mynode) THEN
1425  NULLIFY (p_new_block)
1426  CALL dbcsr_reserve_block2d(no_quench(1), &
1427  iblock_row, iblock_col, p_new_block)
1428  cpassert(ASSOCIATED(p_new_block))
1429  p_new_block(:, :) = 1.0_dp
1430  END IF
1431  END DO
1432  END DO
1433  CALL dbcsr_finalize(no_quench(1))
1434  IF (almo_scf_env%nspins .GT. 1) THEN
1435  DO ispin = 2, almo_scf_env%nspins
1436  CALL dbcsr_create(no_quench(ispin), &
1437  template=almo_scf_env%matrix_t(1), &
1438  matrix_type=dbcsr_type_no_symmetry)
1439  CALL dbcsr_copy(no_quench(ispin), no_quench(1))
1440  END DO
1441  END IF
1442 
1443  END SELECT
1444 
1445  SELECT CASE (almo_scf_env%deloc_method)
1447 
1448  DO ispin = 1, almo_scf_env%nspins
1449  CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
1450  almo_scf_env%matrix_t_blk(ispin))
1451  END DO
1452 
1454 
1455  !!!! RZK-warning a whole class of delocalization methods
1456  !!!! are commented out at the moment because some of their
1457  !!!! routines have not been thoroughly tested.
1458 
1459  !!!! if we have virtuals pre-optimize and truncate them
1460  !!!IF (almo_scf_env%need_virtuals) THEN
1461  !!! SELECT CASE (almo_scf_env%deloc_truncate_virt)
1462  !!! CASE (virt_full)
1463  !!! ! simply copy virtual orbitals from matrix_v_full_blk to matrix_v_blk
1464  !!! DO ispin=1,almo_scf_env%nspins
1465  !!! CALL dbcsr_copy(almo_scf_env%matrix_v_blk(ispin),&
1466  !!! almo_scf_env%matrix_v_full_blk(ispin))
1467  !!! ENDDO
1468  !!! CASE (virt_number,virt_occ_size)
1469  !!! CALL split_v_blk(almo_scf_env)
1470  !!! !CALL truncate_subspace_v_blk(qs_env,almo_scf_env)
1471  !!! CASE DEFAULT
1472  !!! CPErrorMessage(cp_failure_level,routineP,"illegal method for virtual space truncation")
1473  !!! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
1474  !!! END SELECT
1475  !!!ENDIF
1476  !!!CALL harris_foulkes_correction(qs_env,almo_scf_env)
1477 
1478  IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_pcg) THEN
1479 
1480  CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1481  almo_scf_env=almo_scf_env, &
1482  optimizer=almo_scf_env%opt_xalmo_pcg, &
1483  quench_t=no_quench, &
1484  matrix_t_in=almo_scf_env%matrix_t_blk, &
1485  matrix_t_out=almo_scf_env%matrix_t, &
1486  assume_t0_q0x=(almo_scf_env%xalmo_trial_wf .EQ. xalmo_trial_r0_out), &
1487  perturbation_only=.true., &
1488  special_case=xalmo_case_fully_deloc)
1489 
1490  ELSE IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_trustr) THEN
1491 
1492  CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
1493  almo_scf_env=almo_scf_env, &
1494  optimizer=almo_scf_env%opt_xalmo_trustr, &
1495  quench_t=no_quench, &
1496  matrix_t_in=almo_scf_env%matrix_t_blk, &
1497  matrix_t_out=almo_scf_env%matrix_t, &
1498  perturbation_only=.true., &
1499  special_case=xalmo_case_fully_deloc)
1500 
1501  ELSE
1502 
1503  cpabort("Other algorithms do not exist")
1504 
1505  END IF
1506 
1507  CASE (almo_deloc_xalmo_1diag)
1508 
1509  IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_diag) THEN
1510 
1511  almo_scf_env%perturbative_delocalization = .true.
1512  DO ispin = 1, almo_scf_env%nspins
1513  CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
1514  almo_scf_env%matrix_t_blk(ispin))
1515  END DO
1516  CALL almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, &
1517  arbitrary_optimizer)
1518 
1519  ELSE
1520 
1521  cpabort("Other algorithms do not exist")
1522 
1523  END IF
1524 
1525  CASE (almo_deloc_xalmo_x)
1526 
1527  IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_pcg) THEN
1528 
1529  CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1530  almo_scf_env=almo_scf_env, &
1531  optimizer=almo_scf_env%opt_xalmo_pcg, &
1532  quench_t=almo_scf_env%quench_t, &
1533  matrix_t_in=almo_scf_env%matrix_t_blk, &
1534  matrix_t_out=almo_scf_env%matrix_t, &
1535  assume_t0_q0x=(almo_scf_env%xalmo_trial_wf .EQ. xalmo_trial_r0_out), &
1536  perturbation_only=.true., &
1537  special_case=xalmo_case_normal)
1538 
1539  ELSE IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_trustr) THEN
1540 
1541  CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
1542  almo_scf_env=almo_scf_env, &
1543  optimizer=almo_scf_env%opt_xalmo_trustr, &
1544  quench_t=almo_scf_env%quench_t, &
1545  matrix_t_in=almo_scf_env%matrix_t_blk, &
1546  matrix_t_out=almo_scf_env%matrix_t, &
1547  perturbation_only=.true., &
1548  special_case=xalmo_case_normal)
1549 
1550  ELSE
1551 
1552  cpabort("Other algorithms do not exist")
1553 
1554  END IF
1555 
1556  CASE (almo_deloc_xalmo_scf)
1557 
1558  IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_diag) THEN
1559 
1560  cpabort("Should not be here: convergence will fail!")
1561 
1562  almo_scf_env%perturbative_delocalization = .false.
1563  DO ispin = 1, almo_scf_env%nspins
1564  CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
1565  almo_scf_env%matrix_t_blk(ispin))
1566  END DO
1567  CALL almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, &
1568  arbitrary_optimizer)
1569 
1570  ELSE IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_pcg) THEN
1571 
1572  CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1573  almo_scf_env=almo_scf_env, &
1574  optimizer=almo_scf_env%opt_xalmo_pcg, &
1575  quench_t=almo_scf_env%quench_t, &
1576  matrix_t_in=almo_scf_env%matrix_t_blk, &
1577  matrix_t_out=almo_scf_env%matrix_t, &
1578  assume_t0_q0x=(almo_scf_env%xalmo_trial_wf .EQ. xalmo_trial_r0_out), &
1579  perturbation_only=.false., &
1580  special_case=xalmo_case_normal)
1581 
1582  ! RZK-warning THIS IS A HACK TO GET ORBITAL ENERGIES
1583  almo_experimental = .false.
1584  IF (almo_experimental) THEN
1585  almo_scf_env%perturbative_delocalization = .true.
1586  !DO ispin=1,almo_scf_env%nspins
1587  ! CALL dbcsr_copy(almo_scf_env%matrix_t(ispin),&
1588  ! almo_scf_env%matrix_t_blk(ispin))
1589  !ENDDO
1590  CALL almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, &
1591  arbitrary_optimizer)
1592  END IF ! experimental
1593 
1594  ELSE IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_trustr) THEN
1595 
1596  CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
1597  almo_scf_env=almo_scf_env, &
1598  optimizer=almo_scf_env%opt_xalmo_trustr, &
1599  quench_t=almo_scf_env%quench_t, &
1600  matrix_t_in=almo_scf_env%matrix_t_blk, &
1601  matrix_t_out=almo_scf_env%matrix_t, &
1602  perturbation_only=.false., &
1603  special_case=xalmo_case_normal)
1604 
1605  ELSE
1606 
1607  cpabort("Other algorithms do not exist")
1608 
1609  END IF
1610 
1611  CASE DEFAULT
1612 
1613  cpabort("Illegal delocalization method")
1614 
1615  END SELECT
1616 
1617  SELECT CASE (almo_scf_env%deloc_method)
1619 
1620  IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
1621  cpabort("full scf is NYI for truncated virtual space")
1622  END IF
1623 
1624  IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_pcg) THEN
1625 
1626  CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1627  almo_scf_env=almo_scf_env, &
1628  optimizer=almo_scf_env%opt_xalmo_pcg, &
1629  quench_t=no_quench, &
1630  matrix_t_in=almo_scf_env%matrix_t, &
1631  matrix_t_out=almo_scf_env%matrix_t, &
1632  assume_t0_q0x=.false., &
1633  perturbation_only=.false., &
1634  special_case=xalmo_case_fully_deloc)
1635 
1636  ELSE IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_trustr) THEN
1637 
1638  CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
1639  almo_scf_env=almo_scf_env, &
1640  optimizer=almo_scf_env%opt_xalmo_trustr, &
1641  quench_t=no_quench, &
1642  matrix_t_in=almo_scf_env%matrix_t, &
1643  matrix_t_out=almo_scf_env%matrix_t, &
1644  perturbation_only=.false., &
1645  special_case=xalmo_case_fully_deloc)
1646 
1647  ELSE
1648 
1649  cpabort("Other algorithms do not exist")
1650 
1651  END IF
1652 
1653  END SELECT
1654 
1655  ! clean up
1656  SELECT CASE (almo_scf_env%deloc_method)
1658  DO ispin = 1, almo_scf_env%nspins
1659  CALL dbcsr_release(no_quench(ispin))
1660  END DO
1661  DEALLOCATE (no_quench)
1662  END SELECT
1663 
1664  CALL timestop(handle)
1665 
1666  END SUBROUTINE almo_scf_delocalization
1667 
1668 ! **************************************************************************************************
1669 !> \brief orbital localization
1670 !> \param qs_env ...
1671 !> \param almo_scf_env ...
1672 !> \par History
1673 !> 2018.09 created [Ziling Luo]
1674 !> \author Ziling Luo
1675 ! **************************************************************************************************
1676  SUBROUTINE construct_nlmos(qs_env, almo_scf_env)
1677 
1678  TYPE(qs_environment_type), POINTER :: qs_env
1679  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1680 
1681  INTEGER :: ispin
1682 
1683  IF (almo_scf_env%construct_nlmos) THEN
1684 
1685  DO ispin = 1, almo_scf_env%nspins
1686 
1687  CALL orthogonalize_mos(ket=almo_scf_env%matrix_t(ispin), &
1688  overlap=almo_scf_env%matrix_sigma(ispin), &
1689  metric=almo_scf_env%matrix_s(1), &
1690  retain_locality=.false., &
1691  only_normalize=.false., &
1692  nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1693  eps_filter=almo_scf_env%eps_filter, &
1694  order_lanczos=almo_scf_env%order_lanczos, &
1695  eps_lanczos=almo_scf_env%eps_lanczos, &
1696  max_iter_lanczos=almo_scf_env%max_iter_lanczos)
1697  END DO
1698 
1699  CALL construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals=.false.)
1700 
1701  IF (almo_scf_env%opt_nlmo_pcg%opt_penalty%virtual_nlmos) THEN
1702  CALL construct_virtuals(almo_scf_env)
1703  CALL construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals=.true.)
1704  END IF
1705 
1706  IF (almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start .GT. 0.0_dp) THEN
1707  CALL nlmo_compactification(qs_env, almo_scf_env, almo_scf_env%matrix_t)
1708  END IF
1709 
1710  END IF
1711 
1712  END SUBROUTINE construct_nlmos
1713 
1714 ! **************************************************************************************************
1715 !> \brief Calls NLMO optimization
1716 !> \param qs_env ...
1717 !> \param almo_scf_env ...
1718 !> \param virtuals ...
1719 !> \par History
1720 !> 2019.10 created [Ziling Luo]
1721 !> \author Ziling Luo
1722 ! **************************************************************************************************
1723  SUBROUTINE construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals)
1724 
1725  TYPE(qs_environment_type), POINTER :: qs_env
1726  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1727  LOGICAL, INTENT(IN) :: virtuals
1728 
1729  REAL(kind=dp) :: det_diff, prev_determinant
1730 
1731  almo_scf_env%overlap_determinant = 1.0
1732  ! KEEP: initial_vol_coeff = almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength
1733  almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength = &
1734  -1.0_dp*almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength !NEW1
1735 
1736  ! loop over the strength of the orthogonalization penalty
1737  prev_determinant = 10.0_dp
1738  DO WHILE (almo_scf_env%overlap_determinant .GT. almo_scf_env%opt_nlmo_pcg%opt_penalty%final_determinant)
1739 
1740  IF (.NOT. virtuals) THEN
1741  CALL almo_scf_construct_nlmos(qs_env=qs_env, &
1742  optimizer=almo_scf_env%opt_nlmo_pcg, &
1743  matrix_s=almo_scf_env%matrix_s(1), &
1744  matrix_mo_in=almo_scf_env%matrix_t, &
1745  matrix_mo_out=almo_scf_env%matrix_t, &
1746  template_matrix_sigma=almo_scf_env%matrix_sigma_inv, &
1747  overlap_determinant=almo_scf_env%overlap_determinant, &
1748  mat_distr_aos=almo_scf_env%mat_distr_aos, &
1749  virtuals=virtuals, &
1750  eps_filter=almo_scf_env%eps_filter)
1751  ELSE
1752  CALL almo_scf_construct_nlmos(qs_env=qs_env, &
1753  optimizer=almo_scf_env%opt_nlmo_pcg, &
1754  matrix_s=almo_scf_env%matrix_s(1), &
1755  matrix_mo_in=almo_scf_env%matrix_v, &
1756  matrix_mo_out=almo_scf_env%matrix_v, &
1757  template_matrix_sigma=almo_scf_env%matrix_sigma_vv, &
1758  overlap_determinant=almo_scf_env%overlap_determinant, &
1759  mat_distr_aos=almo_scf_env%mat_distr_aos, &
1760  virtuals=virtuals, &
1761  eps_filter=almo_scf_env%eps_filter)
1762 
1763  END IF
1764 
1765  det_diff = prev_determinant - almo_scf_env%overlap_determinant
1766  almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength = &
1767  almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength/ &
1768  abs(almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength_dec_factor)
1769 
1770  IF (det_diff < almo_scf_env%opt_nlmo_pcg%opt_penalty%determinant_tolerance) THEN
1771  EXIT
1772  END IF
1773  prev_determinant = almo_scf_env%overlap_determinant
1774 
1775  END DO
1776 
1777  END SUBROUTINE construct_nlmos_wrapper
1778 
1779 ! **************************************************************************************************
1780 !> \brief Construct virtual orbitals
1781 !> \param almo_scf_env ...
1782 !> \par History
1783 !> 2019.10 created [Ziling Luo]
1784 !> \author Ziling Luo
1785 ! **************************************************************************************************
1786  SUBROUTINE construct_virtuals(almo_scf_env)
1787 
1788  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1789 
1790  INTEGER :: ispin, n
1791  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1792  TYPE(dbcsr_type) :: tempnv1, tempvocc1, tempvocc2, tempvv1, &
1793  tempvv2
1794 
1795  DO ispin = 1, almo_scf_env%nspins
1796 
1797  CALL dbcsr_create(tempnv1, &
1798  template=almo_scf_env%matrix_v(ispin), &
1799  matrix_type=dbcsr_type_no_symmetry)
1800  CALL dbcsr_create(tempvocc1, &
1801  template=almo_scf_env%matrix_vo(ispin), &
1802  matrix_type=dbcsr_type_no_symmetry)
1803  CALL dbcsr_create(tempvocc2, &
1804  template=almo_scf_env%matrix_vo(ispin), &
1805  matrix_type=dbcsr_type_no_symmetry)
1806  CALL dbcsr_create(tempvv1, &
1807  template=almo_scf_env%matrix_sigma_vv(ispin), &
1808  matrix_type=dbcsr_type_no_symmetry)
1809  CALL dbcsr_create(tempvv2, &
1810  template=almo_scf_env%matrix_sigma_vv(ispin), &
1811  matrix_type=dbcsr_type_no_symmetry)
1812 
1813  ! Generate random virtual matrix
1814  CALL dbcsr_init_random(almo_scf_env%matrix_v(ispin), &
1815  keep_sparsity=.false.)
1816 
1817  ! Project the orbital subspace out
1818  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1819  almo_scf_env%matrix_s(1), &
1820  almo_scf_env%matrix_v(ispin), &
1821  0.0_dp, tempnv1, &
1822  filter_eps=almo_scf_env%eps_filter)
1823 
1824  CALL dbcsr_multiply("T", "N", 1.0_dp, &
1825  tempnv1, &
1826  almo_scf_env%matrix_t(ispin), &
1827  0.0_dp, tempvocc1, &
1828  filter_eps=almo_scf_env%eps_filter)
1829 
1830  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1831  tempvocc1, &
1832  almo_scf_env%matrix_sigma_inv(ispin), &
1833  0.0_dp, tempvocc2, &
1834  filter_eps=almo_scf_env%eps_filter)
1835 
1836  CALL dbcsr_multiply("N", "T", 1.0_dp, &
1837  almo_scf_env%matrix_t(ispin), &
1838  tempvocc2, &
1839  0.0_dp, tempnv1, &
1840  filter_eps=almo_scf_env%eps_filter)
1841 
1842  CALL dbcsr_add(almo_scf_env%matrix_v(ispin), tempnv1, 1.0_dp, -1.0_dp)
1843 
1844  ! compute VxV overlap
1845  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1846  almo_scf_env%matrix_s(1), &
1847  almo_scf_env%matrix_v(ispin), &
1848  0.0_dp, tempnv1, &
1849  filter_eps=almo_scf_env%eps_filter)
1850 
1851  CALL dbcsr_multiply("T", "N", 1.0_dp, &
1852  almo_scf_env%matrix_v(ispin), &
1853  tempnv1, &
1854  0.0_dp, tempvv1, &
1855  filter_eps=almo_scf_env%eps_filter)
1856 
1857  CALL orthogonalize_mos(ket=almo_scf_env%matrix_v(ispin), &
1858  overlap=tempvv1, &
1859  metric=almo_scf_env%matrix_s(1), &
1860  retain_locality=.false., &
1861  only_normalize=.false., &
1862  nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1863  eps_filter=almo_scf_env%eps_filter, &
1864  order_lanczos=almo_scf_env%order_lanczos, &
1865  eps_lanczos=almo_scf_env%eps_lanczos, &
1866  max_iter_lanczos=almo_scf_env%max_iter_lanczos)
1867 
1868  ! compute VxV block of the KS matrix
1869  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1870  almo_scf_env%matrix_ks(ispin), &
1871  almo_scf_env%matrix_v(ispin), &
1872  0.0_dp, tempnv1, &
1873  filter_eps=almo_scf_env%eps_filter)
1874 
1875  CALL dbcsr_multiply("T", "N", 1.0_dp, &
1876  almo_scf_env%matrix_v(ispin), &
1877  tempnv1, &
1878  0.0_dp, tempvv1, &
1879  filter_eps=almo_scf_env%eps_filter)
1880 
1881  CALL dbcsr_get_info(tempvv1, nfullrows_total=n)
1882  ALLOCATE (eigenvalues(n))
1883  CALL cp_dbcsr_syevd(tempvv1, tempvv2, &
1884  eigenvalues, &
1885  para_env=almo_scf_env%para_env, &
1886  blacs_env=almo_scf_env%blacs_env)
1887  DEALLOCATE (eigenvalues)
1888 
1889  CALL dbcsr_multiply("N", "N", 1.0_dp, &
1890  almo_scf_env%matrix_v(ispin), &
1891  tempvv2, &
1892  0.0_dp, tempnv1, &
1893  filter_eps=almo_scf_env%eps_filter)
1894 
1895  CALL dbcsr_copy(almo_scf_env%matrix_v(ispin), tempnv1)
1896 
1897  CALL dbcsr_release(tempnv1)
1898  CALL dbcsr_release(tempvocc1)
1899  CALL dbcsr_release(tempvocc2)
1900  CALL dbcsr_release(tempvv1)
1901  CALL dbcsr_release(tempvv2)
1902 
1903  END DO
1904 
1905  END SUBROUTINE construct_virtuals
1906 
1907 ! **************************************************************************************************
1908 !> \brief Compactify (set small blocks to zero) orbitals
1909 !> \param qs_env ...
1910 !> \param almo_scf_env ...
1911 !> \param matrix ...
1912 !> \par History
1913 !> 2019.10 created [Ziling Luo]
1914 !> \author Ziling Luo
1915 ! **************************************************************************************************
1916  SUBROUTINE nlmo_compactification(qs_env, almo_scf_env, matrix)
1917 
1918  TYPE(qs_environment_type), POINTER :: qs_env
1919  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1920  TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:), &
1921  INTENT(IN) :: matrix
1922 
1923  INTEGER :: iblock_col, iblock_col_size, iblock_row, iblock_row_size, icol, irow, ispin, &
1924  ncols, nrows, nspins, para_group_handle, unit_nr
1925  LOGICAL :: element_by_element
1926  REAL(kind=dp) :: energy, eps_local, eps_start, &
1927  max_element, spin_factor
1928  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: occ, retained
1929  REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p
1930  TYPE(cp_logger_type), POINTER :: logger
1931  TYPE(dbcsr_iterator_type) :: iter
1932  TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_p_tmp, matrix_t_tmp
1933  TYPE(mp_comm_type) :: para_group
1934 
1935  ! define the output_unit
1936  logger => cp_get_default_logger()
1937  IF (logger%para_env%is_source()) THEN
1938  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1939  ELSE
1940  unit_nr = -1
1941  END IF
1942 
1943  nspins = SIZE(matrix)
1944  element_by_element = .false.
1945 
1946  IF (nspins .EQ. 1) THEN
1947  spin_factor = 2.0_dp
1948  ELSE
1949  spin_factor = 1.0_dp
1950  END IF
1951 
1952  ALLOCATE (matrix_t_tmp(nspins))
1953  ALLOCATE (matrix_p_tmp(nspins))
1954  ALLOCATE (retained(nspins))
1955  ALLOCATE (occ(2))
1956 
1957  DO ispin = 1, nspins
1958 
1959  ! init temporary storage
1960  CALL dbcsr_create(matrix_t_tmp(ispin), &
1961  template=matrix(ispin), &
1962  matrix_type=dbcsr_type_no_symmetry)
1963  CALL dbcsr_copy(matrix_t_tmp(ispin), matrix(ispin))
1964 
1965  CALL dbcsr_create(matrix_p_tmp(ispin), &
1966  template=almo_scf_env%matrix_p(ispin), &
1967  matrix_type=dbcsr_type_no_symmetry)
1968  CALL dbcsr_copy(matrix_p_tmp(ispin), almo_scf_env%matrix_p(ispin))
1969 
1970  END DO
1971 
1972  IF (unit_nr > 0) THEN
1973  WRITE (unit_nr, *)
1974  WRITE (unit_nr, '(T2,A)') &
1975  "Energy dependence on the (block-by-block) filtering of the NLMO coefficients"
1976  IF (unit_nr > 0) WRITE (unit_nr, '(T2,A13,A20,A20,A25)') &
1977  "EPS filter", "Occupation Alpha", "Occupation Beta", "Energy"
1978  END IF
1979 
1980  eps_start = almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start
1981  eps_local = max(eps_start, 10e-14_dp)
1982 
1983  DO
1984 
1985  IF (eps_local > 0.11_dp) EXIT
1986 
1987  DO ispin = 1, nspins
1988 
1989  retained(ispin) = 0
1990  CALL dbcsr_work_create(matrix_t_tmp(ispin), work_mutable=.true.)
1991  CALL dbcsr_iterator_start(iter, matrix_t_tmp(ispin))
1992  DO WHILE (dbcsr_iterator_blocks_left(iter))
1993  CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, &
1994  row_size=iblock_row_size, col_size=iblock_col_size)
1995  DO icol = 1, iblock_col_size
1996 
1997  IF (element_by_element) THEN
1998 
1999  DO irow = 1, iblock_row_size
2000  IF (abs(data_p(irow, icol)) .LT. eps_local) THEN
2001  data_p(irow, icol) = 0.0_dp
2002  ELSE
2003  retained(ispin) = retained(ispin) + 1
2004  END IF
2005  END DO
2006 
2007  ELSE ! rows are blocked
2008 
2009  max_element = 0.0_dp
2010  DO irow = 1, iblock_row_size
2011  IF (abs(data_p(irow, icol)) .GT. max_element) THEN
2012  max_element = abs(data_p(irow, icol))
2013  END IF
2014  END DO
2015  IF (max_element .LT. eps_local) THEN
2016  DO irow = 1, iblock_row_size
2017  data_p(irow, icol) = 0.0_dp
2018  END DO
2019  ELSE
2020  retained(ispin) = retained(ispin) + iblock_row_size
2021  END IF
2022 
2023  END IF ! block rows?
2024  END DO ! icol
2025 
2026  END DO ! iterator
2027  CALL dbcsr_iterator_stop(iter)
2028  CALL dbcsr_finalize(matrix_t_tmp(ispin))
2029  CALL dbcsr_filter(matrix_t_tmp(ispin), eps_local)
2030 
2031  CALL dbcsr_get_info(matrix_t_tmp(ispin), group=para_group_handle, &
2032  nfullrows_total=nrows, &
2033  nfullcols_total=ncols)
2034  CALL para_group%set_handle(para_group_handle)
2035  CALL para_group%sum(retained(ispin))
2036 
2037  !devide by the total no. elements
2038  occ(ispin) = retained(ispin)/nrows/ncols
2039 
2040  ! compute the global projectors (for the density matrix)
2041  CALL almo_scf_t_to_proj( &
2042  t=matrix_t_tmp(ispin), &
2043  p=matrix_p_tmp(ispin), &
2044  eps_filter=almo_scf_env%eps_filter, &
2045  orthog_orbs=.false., &
2046  nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
2047  s=almo_scf_env%matrix_s(1), &
2048  sigma=almo_scf_env%matrix_sigma(ispin), &
2049  sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
2050  use_guess=.false., &
2051  algorithm=almo_scf_env%sigma_inv_algorithm, &
2052  inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
2053  inverse_accelerator=almo_scf_env%order_lanczos, &
2054  eps_lanczos=almo_scf_env%eps_lanczos, &
2055  max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
2056  para_env=almo_scf_env%para_env, &
2057  blacs_env=almo_scf_env%blacs_env)
2058 
2059  ! compute dm from the projector(s)
2060  CALL dbcsr_scale(matrix_p_tmp(ispin), spin_factor)
2061 
2062  END DO
2063 
2064  ! the KS matrix is updated outside the spin loop
2065  CALL almo_dm_to_almo_ks(qs_env, &
2066  matrix_p_tmp, &
2067  almo_scf_env%matrix_ks, &
2068  energy, &
2069  almo_scf_env%eps_filter, &
2070  almo_scf_env%mat_distr_aos)
2071 
2072  IF (nspins .LT. 2) occ(2) = occ(1)
2073  IF (unit_nr > 0) WRITE (unit_nr, '(T2,E13.3,F20.10,F20.10,F25.15)') &
2074  eps_local, occ(1), occ(2), energy
2075 
2076  eps_local = 2.0_dp*eps_local
2077 
2078  END DO
2079 
2080  DO ispin = 1, nspins
2081 
2082  CALL dbcsr_release(matrix_t_tmp(ispin))
2083  CALL dbcsr_release(matrix_p_tmp(ispin))
2084 
2085  END DO
2086 
2087  DEALLOCATE (matrix_t_tmp)
2088  DEALLOCATE (matrix_p_tmp)
2089  DEALLOCATE (occ)
2090  DEALLOCATE (retained)
2091 
2092  END SUBROUTINE nlmo_compactification
2093 
2094 ! *****************************************************************************
2095 !> \brief after SCF we have the final density and KS matrices compute various
2096 !> post-scf quantities
2097 !> \param qs_env ...
2098 !> \param almo_scf_env ...
2099 !> \par History
2100 !> 2015.03 created [Rustam Z. Khaliullin]
2101 !> \author Rustam Z. Khaliullin
2102 ! **************************************************************************************************
2103  SUBROUTINE almo_scf_post(qs_env, almo_scf_env)
2104  TYPE(qs_environment_type), POINTER :: qs_env
2105  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
2106 
2107  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_post'
2108 
2109  INTEGER :: handle, ispin
2110  TYPE(cp_fm_type), POINTER :: mo_coeff
2111  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w
2112  TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_t_processed
2113  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2114  TYPE(qs_scf_env_type), POINTER :: scf_env
2115 
2116  CALL timeset(routinen, handle)
2117 
2118  ! store matrices to speed up the next scf run
2119  CALL almo_scf_store_extrapolation_data(almo_scf_env)
2120 
2121  ! orthogonalize orbitals before returning them to QS
2122  ALLOCATE (matrix_t_processed(almo_scf_env%nspins))
2123  !ALLOCATE (matrix_v_processed(almo_scf_env%nspins))
2124 
2125  DO ispin = 1, almo_scf_env%nspins
2126 
2127  CALL dbcsr_create(matrix_t_processed(ispin), &
2128  template=almo_scf_env%matrix_t(ispin), &
2129  matrix_type=dbcsr_type_no_symmetry)
2130 
2131  CALL dbcsr_copy(matrix_t_processed(ispin), &
2132  almo_scf_env%matrix_t(ispin))
2133 
2134  !CALL dbcsr_create(matrix_v_processed(ispin), &
2135  ! template=almo_scf_env%matrix_v(ispin), &
2136  ! matrix_type=dbcsr_type_no_symmetry)
2137 
2138  !CALL dbcsr_copy(matrix_v_processed(ispin), &
2139  ! almo_scf_env%matrix_v(ispin))
2140 
2141  IF (almo_scf_env%return_orthogonalized_mos) THEN
2142 
2143  CALL orthogonalize_mos(ket=matrix_t_processed(ispin), &
2144  overlap=almo_scf_env%matrix_sigma(ispin), &
2145  metric=almo_scf_env%matrix_s(1), &
2146  retain_locality=.false., &
2147  only_normalize=.false., &
2148  nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
2149  eps_filter=almo_scf_env%eps_filter, &
2150  order_lanczos=almo_scf_env%order_lanczos, &
2151  eps_lanczos=almo_scf_env%eps_lanczos, &
2152  max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
2153  smear=almo_scf_env%smear)
2154  END IF
2155 
2156  END DO
2157 
2158  !! RS-WARNING: If smearing ALMO is requested, rescaled fully-occupied orbitals are returned to QS
2159  !! RS-WARNING: Beware that QS will not be informed about electronic entropy.
2160  !! If you want a quick and dirty transfer to QS energy, uncomment the following hack:
2161  !! IF (almo_scf_env%smear) THEN
2162  !! qs_env%energy%kTS = 0.0_dp
2163  !! DO ispin = 1, almo_scf_env%nspins
2164  !! qs_env%energy%kTS = qs_env%energy%kTS + almo_scf_env%kTS(ispin)
2165  !! END DO
2166  !! END IF
2167 
2168  ! return orbitals to QS
2169  NULLIFY (mos, mo_coeff, scf_env)
2170 
2171  CALL get_qs_env(qs_env, mos=mos, scf_env=scf_env)
2172 
2173  DO ispin = 1, almo_scf_env%nspins
2174  ! Currently only fm version of mo_set is usable.
2175  ! First transform the matrix_t to fm version
2176  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
2177  CALL copy_dbcsr_to_fm(matrix_t_processed(ispin), mo_coeff)
2178  CALL dbcsr_release(matrix_t_processed(ispin))
2179  END DO
2180  DO ispin = 1, almo_scf_env%nspins
2181  CALL dbcsr_release(matrix_t_processed(ispin))
2182  END DO
2183  DEALLOCATE (matrix_t_processed)
2184 
2185  ! calculate post scf properties
2186  ! CALL almo_post_scf_compute_properties(qs_env, almo_scf_env)
2187  CALL almo_post_scf_compute_properties(qs_env)
2188 
2189  ! compute the W matrix if associated
2190  IF (almo_scf_env%calc_forces) THEN
2191  CALL get_qs_env(qs_env, matrix_w=matrix_w)
2192  IF (ASSOCIATED(matrix_w)) THEN
2193  CALL calculate_w_matrix_almo(matrix_w, almo_scf_env)
2194  ELSE
2195  cpabort("Matrix W is needed but not associated")
2196  END IF
2197  END IF
2198 
2199  CALL timestop(handle)
2200 
2201  END SUBROUTINE almo_scf_post
2202 
2203 ! **************************************************************************************************
2204 !> \brief create various matrices
2205 !> \param almo_scf_env ...
2206 !> \param matrix_s0 ...
2207 !> \par History
2208 !> 2011.07 created [Rustam Z Khaliullin]
2209 !> \author Rustam Z Khaliullin
2210 ! **************************************************************************************************
2211  SUBROUTINE almo_scf_env_create_matrices(almo_scf_env, matrix_s0)
2212 
2213  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
2214  TYPE(dbcsr_type), INTENT(IN) :: matrix_s0
2215 
2216  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_env_create_matrices'
2217 
2218  INTEGER :: handle, ispin, nspins
2219 
2220  CALL timeset(routinen, handle)
2221 
2222  nspins = almo_scf_env%nspins
2223 
2224  ! AO overlap matrix and its various functions
2225  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s(1), &
2226  matrix_qs=matrix_s0, &
2227  almo_scf_env=almo_scf_env, &
2228  name_new="S", &
2229  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2230  symmetry_new=dbcsr_type_symmetric, &
2231  spin_key=0, &
2232  init_domains=.false.)
2233  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk(1), &
2234  matrix_qs=matrix_s0, &
2235  almo_scf_env=almo_scf_env, &
2236  name_new="S_BLK", &
2237  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2238  symmetry_new=dbcsr_type_symmetric, &
2239  spin_key=0, &
2240  init_domains=.true.)
2241  IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_diag) THEN
2242  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk_sqrt_inv(1), &
2243  matrix_qs=matrix_s0, &
2244  almo_scf_env=almo_scf_env, &
2245  name_new="S_BLK_SQRT_INV", &
2246  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2247  symmetry_new=dbcsr_type_symmetric, &
2248  spin_key=0, &
2249  init_domains=.true.)
2250  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk_sqrt(1), &
2251  matrix_qs=matrix_s0, &
2252  almo_scf_env=almo_scf_env, &
2253  name_new="S_BLK_SQRT", &
2254  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2255  symmetry_new=dbcsr_type_symmetric, &
2256  spin_key=0, &
2257  init_domains=.true.)
2258  ELSE IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_dm_sign) THEN
2259  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk_inv(1), &
2260  matrix_qs=matrix_s0, &
2261  almo_scf_env=almo_scf_env, &
2262  name_new="S_BLK_INV", &
2263  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2264  symmetry_new=dbcsr_type_symmetric, &
2265  spin_key=0, &
2266  init_domains=.true.)
2267  END IF
2268 
2269  ! MO coeff matrices and their derivatives
2270  ALLOCATE (almo_scf_env%matrix_t_blk(nspins))
2271  ALLOCATE (almo_scf_env%quench_t_blk(nspins))
2272  ALLOCATE (almo_scf_env%matrix_err_blk(nspins))
2273  ALLOCATE (almo_scf_env%matrix_err_xx(nspins))
2274  ALLOCATE (almo_scf_env%matrix_sigma(nspins))
2275  ALLOCATE (almo_scf_env%matrix_sigma_inv(nspins))
2276  ALLOCATE (almo_scf_env%matrix_sigma_sqrt(nspins))
2277  ALLOCATE (almo_scf_env%matrix_sigma_sqrt_inv(nspins))
2278  ALLOCATE (almo_scf_env%matrix_sigma_blk(nspins))
2279  ALLOCATE (almo_scf_env%matrix_sigma_inv_0deloc(nspins))
2280  ALLOCATE (almo_scf_env%matrix_t(nspins))
2281  ALLOCATE (almo_scf_env%matrix_t_tr(nspins))
2282  DO ispin = 1, nspins
2283  ! create the blocked quencher
2284  CALL matrix_almo_create(matrix_new=almo_scf_env%quench_t_blk(ispin), &
2285  matrix_qs=matrix_s0, &
2286  almo_scf_env=almo_scf_env, &
2287  name_new="Q_BLK", &
2288  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_occ/), &
2289  symmetry_new=dbcsr_type_no_symmetry, &
2290  spin_key=ispin, &
2291  init_domains=.true.)
2292  ! create ALMO coefficient matrix
2293  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_t_blk(ispin), &
2294  matrix_qs=matrix_s0, &
2295  almo_scf_env=almo_scf_env, &
2296  name_new="T_BLK", &
2297  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_occ/), &
2298  symmetry_new=dbcsr_type_no_symmetry, &
2299  spin_key=ispin, &
2300  init_domains=.true.)
2301  ! create the error matrix
2302  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_err_blk(ispin), &
2303  matrix_qs=matrix_s0, &
2304  almo_scf_env=almo_scf_env, &
2305  name_new="ERR_BLK", &
2306  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2307  symmetry_new=dbcsr_type_no_symmetry, &
2308  spin_key=ispin, &
2309  init_domains=.true.)
2310  ! create the error matrix for the quenched ALMOs
2311  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_err_xx(ispin), &
2312  matrix_qs=matrix_s0, &
2313  almo_scf_env=almo_scf_env, &
2314  name_new="ERR_XX", &
2315  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_occ/), &
2316  symmetry_new=dbcsr_type_no_symmetry, &
2317  spin_key=ispin, &
2318  init_domains=.false.)
2319  ! create a matrix with dimensions of a transposed mo coefficient matrix
2320  ! it might be necessary to perform the correction step using cayley
2321  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_t_tr(ispin), &
2322  matrix_qs=matrix_s0, &
2323  almo_scf_env=almo_scf_env, &
2324  name_new="T_TR", &
2325  size_keys=(/almo_mat_dim_occ, almo_mat_dim_aobasis/), &
2326  symmetry_new=dbcsr_type_no_symmetry, &
2327  spin_key=ispin, &
2328  init_domains=.false.)
2329  ! create mo overlap matrix
2330  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma(ispin), &
2331  matrix_qs=matrix_s0, &
2332  almo_scf_env=almo_scf_env, &
2333  name_new="SIG", &
2334  size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
2335  symmetry_new=dbcsr_type_symmetric, &
2336  spin_key=ispin, &
2337  init_domains=.false.)
2338  ! create blocked mo overlap matrix
2339  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_blk(ispin), &
2340  matrix_qs=matrix_s0, &
2341  almo_scf_env=almo_scf_env, &
2342  name_new="SIG_BLK", &
2343  size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
2344  symmetry_new=dbcsr_type_symmetric, &
2345  spin_key=ispin, &
2346  init_domains=.true.)
2347  ! create blocked inverse mo overlap matrix
2348  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
2349  matrix_qs=matrix_s0, &
2350  almo_scf_env=almo_scf_env, &
2351  name_new="SIGINV_BLK", &
2352  size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
2353  symmetry_new=dbcsr_type_symmetric, &
2354  spin_key=ispin, &
2355  init_domains=.true.)
2356  ! create inverse mo overlap matrix
2357  CALL matrix_almo_create( &
2358  matrix_new=almo_scf_env%matrix_sigma_inv(ispin), &
2359  matrix_qs=matrix_s0, &
2360  almo_scf_env=almo_scf_env, &
2361  name_new="SIGINV", &
2362  size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
2363  symmetry_new=dbcsr_type_symmetric, &
2364  spin_key=ispin, &
2365  init_domains=.false.)
2366  ! create various templates that will be necessary later
2367  CALL matrix_almo_create( &
2368  matrix_new=almo_scf_env%matrix_t(ispin), &
2369  matrix_qs=matrix_s0, &
2370  almo_scf_env=almo_scf_env, &
2371  name_new="T", &
2372  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_occ/), &
2373  symmetry_new=dbcsr_type_no_symmetry, &
2374  spin_key=ispin, &
2375  init_domains=.false.)
2376  CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt(ispin), &
2377  template=almo_scf_env%matrix_sigma(ispin), &
2378  matrix_type=dbcsr_type_no_symmetry)
2379  CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
2380  template=almo_scf_env%matrix_sigma(ispin), &
2381  matrix_type=dbcsr_type_no_symmetry)
2382  END DO
2383 
2384  ! create virtual orbitals if necessary
2385  IF (almo_scf_env%need_virtuals) THEN
2386  ALLOCATE (almo_scf_env%matrix_v_blk(nspins))
2387  ALLOCATE (almo_scf_env%matrix_v_full_blk(nspins))
2388  ALLOCATE (almo_scf_env%matrix_v(nspins))
2389  ALLOCATE (almo_scf_env%matrix_vo(nspins))
2390  ALLOCATE (almo_scf_env%matrix_x(nspins))
2391  ALLOCATE (almo_scf_env%matrix_ov(nspins))
2392  ALLOCATE (almo_scf_env%matrix_ov_full(nspins))
2393  ALLOCATE (almo_scf_env%matrix_sigma_vv(nspins))
2394  ALLOCATE (almo_scf_env%matrix_sigma_vv_blk(nspins))
2395  ALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt(nspins))
2396  ALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt_inv(nspins))
2397  ALLOCATE (almo_scf_env%matrix_vv_full_blk(nspins))
2398 
2399  IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
2400  ALLOCATE (almo_scf_env%matrix_k_blk(nspins))
2401  ALLOCATE (almo_scf_env%matrix_k_blk_ones(nspins))
2402  ALLOCATE (almo_scf_env%matrix_k_tr(nspins))
2403  ALLOCATE (almo_scf_env%matrix_v_disc(nspins))
2404  ALLOCATE (almo_scf_env%matrix_v_disc_blk(nspins))
2405  ALLOCATE (almo_scf_env%matrix_ov_disc(nspins))
2406  ALLOCATE (almo_scf_env%matrix_vv_disc_blk(nspins))
2407  ALLOCATE (almo_scf_env%matrix_vv_disc(nspins))
2408  ALLOCATE (almo_scf_env%opt_k_t_dd(nspins))
2409  ALLOCATE (almo_scf_env%opt_k_t_rr(nspins))
2410  ALLOCATE (almo_scf_env%opt_k_denom(nspins))
2411  END IF
2412 
2413  DO ispin = 1, nspins
2414  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_full_blk(ispin), &
2415  matrix_qs=matrix_s0, &
2416  almo_scf_env=almo_scf_env, &
2417  name_new="V_FULL_BLK", &
2419  symmetry_new=dbcsr_type_no_symmetry, &
2420  spin_key=ispin, &
2421  init_domains=.false.)
2422  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_blk(ispin), &
2423  matrix_qs=matrix_s0, &
2424  almo_scf_env=almo_scf_env, &
2425  name_new="V_BLK", &
2426  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_virt/), &
2427  symmetry_new=dbcsr_type_no_symmetry, &
2428  spin_key=ispin, &
2429  init_domains=.false.)
2430  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v(ispin), &
2431  matrix_qs=matrix_s0, &
2432  almo_scf_env=almo_scf_env, &
2433  name_new="V", &
2434  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_virt/), &
2435  symmetry_new=dbcsr_type_no_symmetry, &
2436  spin_key=ispin, &
2437  init_domains=.false.)
2438  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ov_full(ispin), &
2439  matrix_qs=matrix_s0, &
2440  almo_scf_env=almo_scf_env, &
2441  name_new="OV_FULL", &
2442  size_keys=(/almo_mat_dim_occ, almo_mat_dim_virt_full/), &
2443  symmetry_new=dbcsr_type_no_symmetry, &
2444  spin_key=ispin, &
2445  init_domains=.false.)
2446  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ov(ispin), &
2447  matrix_qs=matrix_s0, &
2448  almo_scf_env=almo_scf_env, &
2449  name_new="OV", &
2450  size_keys=(/almo_mat_dim_occ, almo_mat_dim_virt/), &
2451  symmetry_new=dbcsr_type_no_symmetry, &
2452  spin_key=ispin, &
2453  init_domains=.false.)
2454  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vo(ispin), &
2455  matrix_qs=matrix_s0, &
2456  almo_scf_env=almo_scf_env, &
2457  name_new="VO", &
2458  size_keys=(/almo_mat_dim_virt, almo_mat_dim_occ/), &
2459  symmetry_new=dbcsr_type_no_symmetry, &
2460  spin_key=ispin, &
2461  init_domains=.false.)
2462  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_x(ispin), &
2463  matrix_qs=matrix_s0, &
2464  almo_scf_env=almo_scf_env, &
2465  name_new="VO", &
2466  size_keys=(/almo_mat_dim_virt, almo_mat_dim_occ/), &
2467  symmetry_new=dbcsr_type_no_symmetry, &
2468  spin_key=ispin, &
2469  init_domains=.false.)
2470  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_vv(ispin), &
2471  matrix_qs=matrix_s0, &
2472  almo_scf_env=almo_scf_env, &
2473  name_new="SIG_VV", &
2474  size_keys=(/almo_mat_dim_virt, almo_mat_dim_virt/), &
2475  symmetry_new=dbcsr_type_symmetric, &
2476  spin_key=ispin, &
2477  init_domains=.false.)
2478  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vv_full_blk(ispin), &
2479  matrix_qs=matrix_s0, &
2480  almo_scf_env=almo_scf_env, &
2481  name_new="VV_FULL_BLK", &
2483  symmetry_new=dbcsr_type_no_symmetry, &
2484  spin_key=ispin, &
2485  init_domains=.true.)
2486  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_vv_blk(ispin), &
2487  matrix_qs=matrix_s0, &
2488  almo_scf_env=almo_scf_env, &
2489  name_new="SIG_VV_BLK", &
2490  size_keys=(/almo_mat_dim_virt, almo_mat_dim_virt/), &
2491  symmetry_new=dbcsr_type_symmetric, &
2492  spin_key=ispin, &
2493  init_domains=.true.)
2494  CALL dbcsr_create(almo_scf_env%matrix_sigma_vv_sqrt(ispin), &
2495  template=almo_scf_env%matrix_sigma_vv(ispin), &
2496  matrix_type=dbcsr_type_no_symmetry)
2497  CALL dbcsr_create(almo_scf_env%matrix_sigma_vv_sqrt_inv(ispin), &
2498  template=almo_scf_env%matrix_sigma_vv(ispin), &
2499  matrix_type=dbcsr_type_no_symmetry)
2500 
2501  IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
2502  CALL matrix_almo_create(matrix_new=almo_scf_env%opt_k_t_rr(ispin), &
2503  matrix_qs=matrix_s0, &
2504  almo_scf_env=almo_scf_env, &
2505  name_new="OPT_K_U_RR", &
2506  size_keys=(/almo_mat_dim_virt, almo_mat_dim_virt/), &
2507  symmetry_new=dbcsr_type_no_symmetry, &
2508  spin_key=ispin, &
2509  init_domains=.false.)
2510  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vv_disc(ispin), &
2511  matrix_qs=matrix_s0, &
2512  almo_scf_env=almo_scf_env, &
2513  name_new="VV_DISC", &
2515  symmetry_new=dbcsr_type_symmetric, &
2516  spin_key=ispin, &
2517  init_domains=.false.)
2518  CALL matrix_almo_create(matrix_new=almo_scf_env%opt_k_t_dd(ispin), &
2519  matrix_qs=matrix_s0, &
2520  almo_scf_env=almo_scf_env, &
2521  name_new="OPT_K_U_DD", &
2523  symmetry_new=dbcsr_type_no_symmetry, &
2524  spin_key=ispin, &
2525  init_domains=.false.)
2526  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vv_disc_blk(ispin), &
2527  matrix_qs=matrix_s0, &
2528  almo_scf_env=almo_scf_env, &
2529  name_new="VV_DISC_BLK", &
2531  symmetry_new=dbcsr_type_symmetric, &
2532  spin_key=ispin, &
2533  init_domains=.true.)
2534  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_k_blk(ispin), &
2535  matrix_qs=matrix_s0, &
2536  almo_scf_env=almo_scf_env, &
2537  name_new="K_BLK", &
2538  size_keys=(/almo_mat_dim_virt_disc, almo_mat_dim_virt/), &
2539  symmetry_new=dbcsr_type_no_symmetry, &
2540  spin_key=ispin, &
2541  init_domains=.true.)
2542  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_k_blk_ones(ispin), &
2543  matrix_qs=matrix_s0, &
2544  almo_scf_env=almo_scf_env, &
2545  name_new="K_BLK_1", &
2546  size_keys=(/almo_mat_dim_virt_disc, almo_mat_dim_virt/), &
2547  symmetry_new=dbcsr_type_no_symmetry, &
2548  spin_key=ispin, &
2549  init_domains=.true.)
2550  CALL matrix_almo_create(matrix_new=almo_scf_env%opt_k_denom(ispin), &
2551  matrix_qs=matrix_s0, &
2552  almo_scf_env=almo_scf_env, &
2553  name_new="OPT_K_DENOM", &
2554  size_keys=(/almo_mat_dim_virt_disc, almo_mat_dim_virt/), &
2555  symmetry_new=dbcsr_type_no_symmetry, &
2556  spin_key=ispin, &
2557  init_domains=.false.)
2558  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_k_tr(ispin), &
2559  matrix_qs=matrix_s0, &
2560  almo_scf_env=almo_scf_env, &
2561  name_new="K_TR", &
2562  size_keys=(/almo_mat_dim_virt, almo_mat_dim_virt_disc/), &
2563  symmetry_new=dbcsr_type_no_symmetry, &
2564  spin_key=ispin, &
2565  init_domains=.false.)
2566  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_disc_blk(ispin), &
2567  matrix_qs=matrix_s0, &
2568  almo_scf_env=almo_scf_env, &
2569  name_new="V_DISC_BLK", &
2571  symmetry_new=dbcsr_type_no_symmetry, &
2572  spin_key=ispin, &
2573  init_domains=.false.)
2574  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_disc(ispin), &
2575  matrix_qs=matrix_s0, &
2576  almo_scf_env=almo_scf_env, &
2577  name_new="V_DISC", &
2579  symmetry_new=dbcsr_type_no_symmetry, &
2580  spin_key=ispin, &
2581  init_domains=.false.)
2582  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ov_disc(ispin), &
2583  matrix_qs=matrix_s0, &
2584  almo_scf_env=almo_scf_env, &
2585  name_new="OV_DISC", &
2586  size_keys=(/almo_mat_dim_occ, almo_mat_dim_virt_disc/), &
2587  symmetry_new=dbcsr_type_no_symmetry, &
2588  spin_key=ispin, &
2589  init_domains=.false.)
2590 
2591  END IF ! end need_discarded_virtuals
2592 
2593  END DO ! spin
2594  END IF
2595 
2596  ! create matrices of orbital energies if necessary
2597  IF (almo_scf_env%need_orbital_energies) THEN
2598  ALLOCATE (almo_scf_env%matrix_eoo(nspins))
2599  ALLOCATE (almo_scf_env%matrix_evv_full(nspins))
2600  DO ispin = 1, nspins
2601  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_eoo(ispin), &
2602  matrix_qs=matrix_s0, &
2603  almo_scf_env=almo_scf_env, &
2604  name_new="E_OCC", &
2605  size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
2606  symmetry_new=dbcsr_type_no_symmetry, &
2607  spin_key=ispin, &
2608  init_domains=.false.)
2609  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_evv_full(ispin), &
2610  matrix_qs=matrix_s0, &
2611  almo_scf_env=almo_scf_env, &
2612  name_new="E_VIRT", &
2614  symmetry_new=dbcsr_type_no_symmetry, &
2615  spin_key=ispin, &
2616  init_domains=.false.)
2617  END DO
2618  END IF
2619 
2620  ! Density and KS matrices
2621  ALLOCATE (almo_scf_env%matrix_p(nspins))
2622  ALLOCATE (almo_scf_env%matrix_p_blk(nspins))
2623  ALLOCATE (almo_scf_env%matrix_ks(nspins))
2624  ALLOCATE (almo_scf_env%matrix_ks_blk(nspins))
2625  IF (almo_scf_env%need_previous_ks) &
2626  ALLOCATE (almo_scf_env%matrix_ks_0deloc(nspins))
2627  DO ispin = 1, nspins
2628  ! RZK-warning copy with symmery but remember that this might cause problems
2629  CALL dbcsr_create(almo_scf_env%matrix_p(ispin), &
2630  template=almo_scf_env%matrix_s(1), &
2631  matrix_type=dbcsr_type_symmetric)
2632  CALL dbcsr_create(almo_scf_env%matrix_ks(ispin), &
2633  template=almo_scf_env%matrix_s(1), &
2634  matrix_type=dbcsr_type_symmetric)
2635  IF (almo_scf_env%need_previous_ks) THEN
2636  CALL dbcsr_create(almo_scf_env%matrix_ks_0deloc(ispin), &
2637  template=almo_scf_env%matrix_s(1), &
2638  matrix_type=dbcsr_type_symmetric)
2639  END IF
2640  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_p_blk(ispin), &
2641  matrix_qs=matrix_s0, &
2642  almo_scf_env=almo_scf_env, &
2643  name_new="P_BLK", &
2644  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2645  symmetry_new=dbcsr_type_symmetric, &
2646  spin_key=ispin, &
2647  init_domains=.true.)
2648  CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ks_blk(ispin), &
2649  matrix_qs=matrix_s0, &
2650  almo_scf_env=almo_scf_env, &
2651  name_new="KS_BLK", &
2652  size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2653  symmetry_new=dbcsr_type_symmetric, &
2654  spin_key=ispin, &
2655  init_domains=.true.)
2656  END DO
2657 
2658  CALL timestop(handle)
2659 
2660  END SUBROUTINE almo_scf_env_create_matrices
2661 
2662 ! **************************************************************************************************
2663 !> \brief clean up procedures for almo scf
2664 !> \param almo_scf_env ...
2665 !> \par History
2666 !> 2011.06 created [Rustam Z Khaliullin]
2667 !> 2018.09 smearing support [Ruben Staub]
2668 !> \author Rustam Z Khaliullin
2669 ! **************************************************************************************************
2670  SUBROUTINE almo_scf_clean_up(almo_scf_env)
2671 
2672  TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
2673 
2674  CHARACTER(len=*), PARAMETER :: routinen = 'almo_scf_clean_up'
2675 
2676  INTEGER :: handle, ispin, unit_nr
2677  TYPE(cp_logger_type), POINTER :: logger
2678 
2679  CALL timeset(routinen, handle)
2680 
2681  ! get a useful output_unit
2682  logger => cp_get_default_logger()
2683  IF (logger%para_env%is_source()) THEN
2684  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
2685  ELSE
2686  unit_nr = -1
2687  END IF
2688 
2689  ! release matrices
2690  CALL dbcsr_release(almo_scf_env%matrix_s(1))
2691  CALL dbcsr_release(almo_scf_env%matrix_s_blk(1))
2692  IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_diag) THEN
2693  CALL dbcsr_release(almo_scf_env%matrix_s_blk_sqrt_inv(1))
2694  CALL dbcsr_release(almo_scf_env%matrix_s_blk_sqrt(1))
2695  ELSE IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_dm_sign) THEN
2696  CALL dbcsr_release(almo_scf_env%matrix_s_blk_inv(1))
2697  END IF
2698  DO ispin = 1, almo_scf_env%nspins
2699  CALL dbcsr_release(almo_scf_env%quench_t(ispin))
2700  CALL dbcsr_release(almo_scf_env%quench_t_blk(ispin))
2701  CALL dbcsr_release(almo_scf_env%matrix_t_blk(ispin))
2702  CALL dbcsr_release(almo_scf_env%matrix_err_blk(ispin))
2703  CALL dbcsr_release(almo_scf_env%matrix_err_xx(ispin))
2704  CALL dbcsr_release(almo_scf_env%matrix_t_tr(ispin))
2705  CALL dbcsr_release(almo_scf_env%matrix_sigma(ispin))
2706  CALL dbcsr_release(almo_scf_env%matrix_sigma_blk(ispin))
2707  CALL dbcsr_release(almo_scf_env%matrix_sigma_inv_0deloc(ispin))
2708  CALL dbcsr_release(almo_scf_env%matrix_sigma_inv(ispin))
2709  CALL dbcsr_release(almo_scf_env%matrix_t(ispin))
2710  CALL dbcsr_release(almo_scf_env%matrix_sigma_sqrt(ispin))
2711  CALL dbcsr_release(almo_scf_env%matrix_sigma_sqrt_inv(ispin))
2712  CALL dbcsr_release(almo_scf_env%matrix_p(ispin))
2713  CALL dbcsr_release(almo_scf_env%matrix_ks(ispin))
2714  CALL dbcsr_release(almo_scf_env%matrix_p_blk(ispin))
2715  CALL dbcsr_release(almo_scf_env%matrix_ks_blk(ispin))
2716  IF (almo_scf_env%need_previous_ks) THEN
2717  CALL dbcsr_release(almo_scf_env%matrix_ks_0deloc(ispin))
2718  END IF
2719  IF (almo_scf_env%need_virtuals) THEN
2720  CALL dbcsr_release(almo_scf_env%matrix_v_blk(ispin))
2721  CALL dbcsr_release(almo_scf_env%matrix_v_full_blk(ispin))
2722  CALL dbcsr_release(almo_scf_env%matrix_v(ispin))
2723  CALL dbcsr_release(almo_scf_env%matrix_vo(ispin))
2724  CALL dbcsr_release(almo_scf_env%matrix_x(ispin))
2725  CALL dbcsr_release(almo_scf_env%matrix_ov(ispin))
2726  CALL dbcsr_release(almo_scf_env%matrix_ov_full(ispin))
2727  CALL dbcsr_release(almo_scf_env%matrix_sigma_vv(ispin))
2728  CALL dbcsr_release(almo_scf_env%matrix_sigma_vv_blk(ispin))
2729  CALL dbcsr_release(almo_scf_env%matrix_sigma_vv_sqrt(ispin))
2730  CALL dbcsr_release(almo_scf_env%matrix_sigma_vv_sqrt_inv(ispin))
2731  CALL dbcsr_release(almo_scf_env%matrix_vv_full_blk(ispin))
2732  IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
2733  CALL dbcsr_release(almo_scf_env%matrix_k_tr(ispin))
2734  CALL dbcsr_release(almo_scf_env%matrix_k_blk(ispin))
2735  CALL dbcsr_release(almo_scf_env%matrix_k_blk_ones(ispin))
2736  CALL dbcsr_release(almo_scf_env%matrix_v_disc(ispin))
2737  CALL dbcsr_release(almo_scf_env%matrix_v_disc_blk(ispin))
2738  CALL dbcsr_release(almo_scf_env%matrix_ov_disc(ispin))
2739  CALL dbcsr_release(almo_scf_env%matrix_vv_disc_blk(ispin))
2740  CALL dbcsr_release(almo_scf_env%matrix_vv_disc(ispin))
2741  CALL dbcsr_release(almo_scf_env%opt_k_t_dd(ispin))
2742  CALL dbcsr_release(almo_scf_env%opt_k_t_rr(ispin))
2743  CALL dbcsr_release(almo_scf_env%opt_k_denom(ispin))
2744  END IF
2745  END IF
2746  IF (almo_scf_env%need_orbital_energies) THEN
2747  CALL dbcsr_release(almo_scf_env%matrix_eoo(ispin))
2748  CALL dbcsr_release(almo_scf_env%matrix_evv_full(ispin))
2749  END IF
2750  END DO
2751 
2752  ! deallocate matrices
2753  DEALLOCATE (almo_scf_env%matrix_p)
2754  DEALLOCATE (almo_scf_env%matrix_p_blk)
2755  DEALLOCATE (almo_scf_env%matrix_ks)
2756  DEALLOCATE (almo_scf_env%matrix_ks_blk)
2757  DEALLOCATE (almo_scf_env%matrix_t_blk)
2758  DEALLOCATE (almo_scf_env%matrix_err_blk)
2759  DEALLOCATE (almo_scf_env%matrix_err_xx)
2760  DEALLOCATE (almo_scf_env%matrix_t)
2761  DEALLOCATE (almo_scf_env%matrix_t_tr)
2762  DEALLOCATE (almo_scf_env%matrix_sigma)
2763  DEALLOCATE (almo_scf_env%matrix_sigma_blk)
2764  DEALLOCATE (almo_scf_env%matrix_sigma_inv_0deloc)
2765  DEALLOCATE (almo_scf_env%matrix_sigma_sqrt)
2766  DEALLOCATE (almo_scf_env%matrix_sigma_sqrt_inv)
2767  DEALLOCATE (almo_scf_env%matrix_sigma_inv)
2768  DEALLOCATE (almo_scf_env%quench_t)
2769  DEALLOCATE (almo_scf_env%quench_t_blk)
2770  IF (almo_scf_env%need_virtuals) THEN
2771  DEALLOCATE (almo_scf_env%matrix_v_blk)
2772  DEALLOCATE (almo_scf_env%matrix_v_full_blk)
2773  DEALLOCATE (almo_scf_env%matrix_v)
2774  DEALLOCATE (almo_scf_env%matrix_vo)
2775  DEALLOCATE (almo_scf_env%matrix_x)
2776  DEALLOCATE (almo_scf_env%matrix_ov)
2777  DEALLOCATE (almo_scf_env%matrix_ov_full)
2778  DEALLOCATE (almo_scf_env%matrix_sigma_vv)
2779  DEALLOCATE (almo_scf_env%matrix_sigma_vv_blk)
2780  DEALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt)
2781  DEALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt_inv)
2782  DEALLOCATE (almo_scf_env%matrix_vv_full_blk)
2783  IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
2784  DEALLOCATE (almo_scf_env%matrix_k_tr)
2785  DEALLOCATE (almo_scf_env%matrix_k_blk)
2786  DEALLOCATE (almo_scf_env%matrix_v_disc)
2787  DEALLOCATE (almo_scf_env%matrix_v_disc_blk)
2788  DEALLOCATE (almo_scf_env%matrix_ov_disc)
2789  DEALLOCATE (almo_scf_env%matrix_vv_disc_blk)
2790  DEALLOCATE (almo_scf_env%matrix_vv_disc)
2791  DEALLOCATE (almo_scf_env%matrix_k_blk_ones)
2792  DEALLOCATE (almo_scf_env%opt_k_t_dd)
2793  DEALLOCATE (almo_scf_env%opt_k_t_rr)
2794  DEALLOCATE (almo_scf_env%opt_k_denom)
2795  END IF
2796  END IF
2797  IF (almo_scf_env%need_previous_ks) THEN
2798  DEALLOCATE (almo_scf_env%matrix_ks_0deloc)
2799  END IF
2800  IF (almo_scf_env%need_orbital_energies) THEN
2801  DEALLOCATE (almo_scf_env%matrix_eoo)
2802  DEALLOCATE (almo_scf_env%matrix_evv_full)
2803  END IF
2804 
2805  ! clean up other variables
2806  DO ispin = 1, almo_scf_env%nspins
2807  CALL release_submatrices( &
2808  almo_scf_env%domain_preconditioner(:, ispin))
2809  CALL release_submatrices(almo_scf_env%domain_s_inv(:, ispin))
2810  CALL release_submatrices(almo_scf_env%domain_s_sqrt_inv(:, ispin))
2811  CALL release_submatrices(almo_scf_env%domain_s_sqrt(:, ispin))
2812  CALL release_submatrices(almo_scf_env%domain_ks_xx(:, ispin))
2813  CALL release_submatrices(almo_scf_env%domain_t(:, ispin))
2814  CALL release_submatrices(almo_scf_env%domain_err(:, ispin))
2815  CALL release_submatrices(almo_scf_env%domain_r_down_up(:, ispin))
2816  END DO
2817  DEALLOCATE (almo_scf_env%domain_preconditioner)
2818  DEALLOCATE (almo_scf_env%domain_s_inv)
2819  DEALLOCATE (almo_scf_env%domain_s_sqrt_inv)
2820  DEALLOCATE (almo_scf_env%domain_s_sqrt)
2821  DEALLOCATE (almo_scf_env%domain_ks_xx)
2822  DEALLOCATE (almo_scf_env%domain_t)
2823  DEALLOCATE (almo_scf_env%domain_err)
2824  DEALLOCATE (almo_scf_env%domain_r_down_up)
2825  DO ispin = 1, almo_scf_env%nspins
2826  DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs)
2827  DEALLOCATE (almo_scf_env%domain_map(ispin)%index1)
2828  END DO
2829  DEALLOCATE (almo_scf_env%domain_map)
2830  DEALLOCATE (almo_scf_env%domain_index_of_ao)
2831  DEALLOCATE (almo_scf_env%domain_index_of_atom)
2832  DEALLOCATE (almo_scf_env%first_atom_of_domain)
2833  DEALLOCATE (almo_scf_env%last_atom_of_domain)
2834  DEALLOCATE (almo_scf_env%nbasis_of_domain)
2835  DEALLOCATE (almo_scf_env%nocc_of_domain)
2836  DEALLOCATE (almo_scf_env%real_ne_of_domain)
2837  DEALLOCATE (almo_scf_env%nvirt_full_of_domain)
2838  DEALLOCATE (almo_scf_env%nvirt_of_domain)
2839  DEALLOCATE (almo_scf_env%nvirt_disc_of_domain)
2840  DEALLOCATE (almo_scf_env%mu_of_domain)
2841  DEALLOCATE (almo_scf_env%cpu_of_domain)
2842  DEALLOCATE (almo_scf_env%charge_of_domain)
2843  DEALLOCATE (almo_scf_env%multiplicity_of_domain)
2844  IF (almo_scf_env%smear) THEN
2845  DEALLOCATE (almo_scf_env%mo_energies)
2846  DEALLOCATE (almo_scf_env%kTS)
2847  END IF
2848 
2849  DEALLOCATE (almo_scf_env%domain_index_of_ao_block)
2850  DEALLOCATE (almo_scf_env%domain_index_of_mo_block)
2851 
2852  CALL mp_para_env_release(almo_scf_env%para_env)
2853  CALL cp_blacs_env_release(almo_scf_env%blacs_env)
2854 
2855  CALL timestop(handle)
2856 
2857  END SUBROUTINE almo_scf_clean_up
2858 
2859 ! **************************************************************************************************
2860 !> \brief Do post scf calculations with ALMO
2861 !> WARNING: ALMO post scf calculation may not work for certain quantities,
2862 !> like forces, since ALMO wave function is only 'partially' optimized
2863 !> \param qs_env ...
2864 !> \par History
2865 !> 2016.12 created [Yifei Shi]
2866 !> \author Yifei Shi
2867 ! **************************************************************************************************
2868  SUBROUTINE almo_post_scf_compute_properties(qs_env)
2869  TYPE(qs_environment_type), POINTER :: qs_env
2870 
2871  CALL qs_scf_compute_properties(qs_env)
2872 
2873  END SUBROUTINE almo_post_scf_compute_properties
2874 
2875 END MODULE almo_scf
2876 
Subroutines for ALMO SCF.
subroutine, public distribute_domains(almo_scf_env)
Load balancing of the submatrix computations.
subroutine, public almo_scf_p_blk_to_t_blk(almo_scf_env, ionic)
computes occupied ALMOs from the superimposed atomic density blocks
subroutine, public almo_scf_t_rescaling(matrix_t, mo_energies, mu_of_domain, real_ne_of_domain, spin_kTS, smear_e_temp, ndomains, nocc_of_domain)
Apply an occupation-rescaling trick to ALMOs for smearing. Partially occupied orbitals are considered...
subroutine, public almo_scf_t_to_proj(t, p, eps_filter, orthog_orbs, nocc_of_domain, s, sigma, sigma_inv, use_guess, smear, algorithm, para_env, blacs_env, eps_lanczos, max_iter_lanczos, inverse_accelerator, inv_eps_factor)
computes the idempotent density matrix from MOs MOs can be either orthogonal or non-orthogonal
subroutine, public orthogonalize_mos(ket, overlap, metric, retain_locality, only_normalize, nocc_of_domain, eps_filter, order_lanczos, eps_lanczos, max_iter_lanczos, overlap_sqrti, smear)
orthogonalize MOs
Optimization routines for all ALMO-based SCF methods.
subroutine, public almo_scf_xalmo_trustr(qs_env, almo_scf_env, optimizer, quench_t, matrix_t_in, matrix_t_out, perturbation_only, special_case)
Optimization of ALMOs using trust region minimizers.
subroutine, public almo_scf_xalmo_pcg(qs_env, almo_scf_env, optimizer, quench_t, matrix_t_in, matrix_t_out, assume_t0_q0x, perturbation_only, special_case)
Optimization of ALMOs using PCG-like minimizers.
subroutine, public almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, optimizer)
An eigensolver-based SCF to optimize extended ALMOs (i.e. ALMOs on overlapping domains)
subroutine, public almo_scf_construct_nlmos(qs_env, optimizer, matrix_s, matrix_mo_in, matrix_mo_out, template_matrix_sigma, overlap_determinant, mat_distr_aos, virtuals, eps_filter)
Optimization of NLMOs using PCG minimizers.
subroutine, public almo_scf_block_diagonal(qs_env, almo_scf_env, optimizer)
An SCF procedure that optimizes block-diagonal ALMOs using DIIS.
Interface between ALMO SCF and QS.
Definition: almo_scf_qs.F:14
subroutine, public construct_qs_mos(qs_env, almo_scf_env)
Create MOs in the QS env to be able to return ALMOs to QS.
Definition: almo_scf_qs.F:574
subroutine, public matrix_qs_to_almo(matrix_qs, matrix_almo, mat_distr_aos, keep_sparsity)
convert between two types of matrices: QS style to ALMO style
Definition: almo_scf_qs.F:423
subroutine, public calculate_w_matrix_almo(matrix_w, almo_scf_env)
Compute matrix W (energy-weighted density matrix) that is needed for the evaluation of forces.
Definition: almo_scf_qs.F:1544
subroutine, public matrix_almo_create(matrix_new, matrix_qs, almo_scf_env, name_new, size_keys, symmetry_new, spin_key, init_domains)
create the ALMO matrix templates
Definition: almo_scf_qs.F:117
subroutine, public init_almo_ks_matrix_via_qs(qs_env, matrix_ks, mat_distr_aos, eps_filter)
Initialization of the QS and ALMO KS matrix.
Definition: almo_scf_qs.F:512
subroutine, public almo_dm_to_almo_ks(qs_env, matrix_p, matrix_ks, energy_total, eps_filter, mat_distr_aos, smear, kTS_sum)
uses the ALMO density matrix to compute ALMO KS matrix and the new energy
Definition: almo_scf_qs.F:752
subroutine, public almo_scf_construct_quencher(qs_env, almo_scf_env)
Creates the matrix that imposes absolute locality on MOs.
Definition: almo_scf_qs.F:842
Types for all ALMO-based methods.
integer, parameter, public almo_mat_dim_occ
integer, parameter, public almo_mat_dim_virt_full
integer, parameter, public almo_mat_dim_aobasis
subroutine, public print_optimizer_options(optimizer, unit_nr)
Prints out the options of an optimizer.
integer, parameter, public almo_mat_dim_virt
integer, parameter, public almo_mat_dim_virt_disc
Routines for all ALMO-based SCF methods 'RZK-warning' marks unresolved issues.
Definition: almo_scf.F:15
subroutine, public almo_entry_scf(qs_env, calc_forces)
The entry point into ALMO SCF routines.
Definition: almo_scf.F:123
Define the atomic kind types and their sub types.
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public scheiber2018
Definition: bibliography.F:43
integer, save, public kuhne2007
Definition: bibliography.F:43
integer, save, public staub2019
Definition: bibliography.F:43
integer, save, public khaliullin2013
Definition: bibliography.F:43
integer, save, public kolafa2004
Definition: bibliography.F:43
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition: cp_blacs_env.F:282
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Interface to (sca)lapack for the Cholesky based procedures.
Definition: cp_dbcsr_diag.F:17
subroutine, public cp_dbcsr_syevd(matrix, eigenvectors, eigenvalues, para_env, blacs_env)
...
Definition: cp_dbcsr_diag.F:67
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Subroutines to handle submatrices.
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public smear_fermi_dirac
integer, parameter, public xalmo_case_normal
integer, parameter, public xalmo_trial_r0_out
integer, parameter, public molecular_guess
integer, parameter, public almo_deloc_scf
integer, parameter, public xalmo_case_fully_deloc
integer, parameter, public virt_number
integer, parameter, public xalmo_case_block_diag
integer, parameter, public almo_mat_distr_molecular
integer, parameter, public almo_domain_layout_molecular
integer, parameter, public atomic_guess
integer, parameter, public almo_scf_diag
integer, parameter, public almo_deloc_xk
integer, parameter, public optimizer_diis
integer, parameter, public almo_scf_skip
integer, parameter, public optimizer_lin_eq_pcg
integer, parameter, public almo_deloc_xalmo_x
integer, parameter, public almo_scf_trustr
integer, parameter, public almo_deloc_x
integer, parameter, public almo_scf_dm_sign
integer, parameter, public almo_scf_pcg
integer, parameter, public virt_full
integer, parameter, public almo_mat_distr_atomic
integer, parameter, public almo_deloc_xalmo_scf
integer, parameter, public almo_deloc_xalmo_1diag
integer, parameter, public virt_occ_size
integer, parameter, public almo_deloc_none
integer, parameter, public optimizer_trustr
integer, parameter, public restart_guess
integer, parameter, public almo_deloc_x_then_scf
integer, parameter, public optimizer_pcg
objects that represent the structure of input sections and the data contained in an input section
Routines useful for iterative matrix calculations.
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_path_length
Definition: kinds.F:58
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
Definition: mathlib.F:205
Interface to the message passing library MPI.
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Define the data structure for the molecule information.
subroutine, public get_molecule_set_info(molecule_set, atom_to_mol, mol_to_first_atom, mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity)
returns information about molecules in the set.
Types used to generate the molecular SCF guess.
Definition: mscfg_types.F:14
subroutine, public get_matrix_from_submatrices(mscfg_env, matrix_out, iset)
Creates a distributed matrix from MOs on fragments.
Definition: mscfg_types.F:125
Define the data structure for the particle information.
Routine to return block diagonal density matrix. Blocks correspond to the atomic densities.
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.
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.
Routines to somehow generate an initial guess.
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.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
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
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
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
Utility routines for qs_scf.
subroutine, public qs_scf_compute_properties(qs_env, wf_type, do_mp2)
computes properties for a given hamilonian using the current wfn
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14