1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Routines for all ALMO-based SCF methods
8!>        RZK-warning marks unresolved issues
9!> \par History
10!>       2011.05 created [Rustam Z Khaliullin]
11!> \author Rustam Z Khaliullin
12! **************************************************************************************************
13MODULE almo_scf
14   USE almo_scf_methods,                ONLY: almo_scf_p_blk_to_t_blk,&
15                                              almo_scf_t_rescaling,&
16                                              almo_scf_t_to_proj,&
17                                              distribute_domains,&
18                                              orthogonalize_mos
19   USE almo_scf_optimizer,              ONLY: almo_scf_block_diagonal,&
20                                              almo_scf_xalmo_eigensolver,&
21                                              almo_scf_xalmo_pcg
22   USE almo_scf_qs,                     ONLY: almo_dm_to_almo_ks,&
23                                              almo_scf_construct_quencher,&
24                                              calculate_w_matrix_almo,&
25                                              construct_qs_mos,&
26                                              init_almo_ks_matrix_via_qs,&
27                                              matrix_almo_create,&
28                                              matrix_qs_to_almo
29   USE almo_scf_types,                  ONLY: almo_mat_dim_aobasis,&
30                                              almo_mat_dim_occ,&
31                                              almo_mat_dim_virt,&
32                                              almo_mat_dim_virt_disc,&
33                                              almo_mat_dim_virt_full,&
34                                              almo_scf_env_type,&
35                                              optimizer_options_type,&
36                                              print_optimizer_options
37   USE atomic_kind_types,               ONLY: atomic_kind_type
38   USE bibliography,                    ONLY: Khaliullin2013,&
39                                              Kolafa2004,&
40                                              Scheiber2018,&
41                                              Staub2019,&
42                                              cite_reference
43   USE cp_blacs_env,                    ONLY: cp_blacs_env_release,&
44                                              cp_blacs_env_retain
45   USE cp_control_types,                ONLY: dft_control_type
46   USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm
47   USE cp_fm_types,                     ONLY: cp_fm_type
48   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
49                                              cp_logger_get_default_unit_nr,&
50                                              cp_logger_type
51   USE cp_para_env,                     ONLY: cp_para_env_release,&
52                                              cp_para_env_retain
53   USE cp_para_types,                   ONLY: cp_para_env_type
54   USE dbcsr_api,                       ONLY: &
55        dbcsr_add_on_diag, dbcsr_binary_read, dbcsr_checksum, dbcsr_copy, dbcsr_create, &
56        dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
57        dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_multiply, dbcsr_nblkcols_total, &
58        dbcsr_nblkrows_total, dbcsr_p_type, dbcsr_release, dbcsr_reserve_block2d, dbcsr_scale, &
59        dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_work_create
60   USE domain_submatrix_methods,        ONLY: init_submatrices,&
61                                              release_submatrices
62   USE input_constants,                 ONLY: &
63        almo_deloc_none, almo_deloc_scf, almo_deloc_x, almo_deloc_x_then_scf, &
64        almo_deloc_xalmo_1diag, almo_deloc_xalmo_scf, almo_deloc_xalmo_x, almo_deloc_xk, &
65        almo_domain_layout_molecular, almo_mat_distr_atomic, almo_mat_distr_molecular, &
66        almo_scf_diag, almo_scf_dm_sign, almo_scf_pcg, almo_scf_skip, atomic_guess, &
67        molecular_guess, optimizer_diis, optimizer_lin_eq_pcg, optimizer_pcg, restart_guess, &
68        smear_fermi_dirac, virt_full, virt_number, virt_occ_size, xalmo_case_block_diag, &
69        xalmo_case_fully_deloc, xalmo_case_normal, xalmo_trial_r0_out
70   USE input_section_types,             ONLY: section_vals_type
71   USE iterate_matrix,                  ONLY: invert_Hotelling,&
72                                              matrix_sqrt_Newton_Schulz
73   USE kinds,                           ONLY: default_path_length,&
74                                              dp
75   USE mathlib,                         ONLY: binomial
76   USE molecule_types,                  ONLY: get_molecule_set_info,&
77                                              molecule_type
78   USE mscfg_types,                     ONLY: get_matrix_from_submatrices,&
79                                              molecular_scf_guess_env_type
80   USE particle_types,                  ONLY: particle_type
81   USE qs_environment_types,            ONLY: get_qs_env,&
82                                              qs_environment_type
83   USE qs_initial_guess,                ONLY: calculate_atomic_block_dm,&
84                                              calculate_mopac_dm
85   USE qs_kind_types,                   ONLY: qs_kind_type
86   USE qs_mo_types,                     ONLY: get_mo_set,&
87                                              mo_set_p_type
88   USE qs_rho_types,                    ONLY: qs_rho_get,&
89                                              qs_rho_type
90   USE qs_scf_post_scf,                 ONLY: qs_scf_compute_properties
91   USE qs_scf_types,                    ONLY: qs_scf_env_type,&
92                                              scf_env_release
93#include "./base/base_uses.f90"
94
95   IMPLICIT NONE
96
97   PRIVATE
98
99   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'almo_scf'
100
101   PUBLIC :: almo_entry_scf
102
103   LOGICAL, PARAMETER :: debug_mode = .FALSE.
104   LOGICAL, PARAMETER :: safe_mode = .FALSE.
105
106CONTAINS
107
108! **************************************************************************************************
109!> \brief The entry point into ALMO SCF routines
110!> \param qs_env   pointer to the QS environment
111!> \param calc_forces   calculate forces?
112!> \par History
113!>       2011.05 created [Rustam Z Khaliullin]
114!> \author Rustam Z Khaliullin
115! **************************************************************************************************
116   SUBROUTINE almo_entry_scf(qs_env, calc_forces)
117      TYPE(qs_environment_type), POINTER                 :: qs_env
118      LOGICAL                                            :: calc_forces
119
120      CHARACTER(len=*), PARAMETER :: routineN = 'almo_entry_scf', routineP = moduleN//':'//routineN
121
122      INTEGER                                            :: handle
123      TYPE(almo_scf_env_type), POINTER                   :: almo_scf_env
124
125      CALL timeset(routineN, handle)
126
127      CALL cite_reference(Khaliullin2013)
128
129      ! get a pointer to the almo environment
130      CALL get_qs_env(qs_env, almo_scf_env=almo_scf_env)
131
132      ! initialize scf
133      CALL almo_scf_init(qs_env, almo_scf_env, calc_forces)
134
135      ! create the initial guess for ALMOs
136      CALL almo_scf_initial_guess(qs_env, almo_scf_env)
137
138      ! perform SCF for block diagonal ALMOs
139      CALL almo_scf_main(qs_env, almo_scf_env)
140
141      ! allow electron delocalization
142      CALL almo_scf_delocalization(qs_env, almo_scf_env)
143
144      ! electron correlation methods
145      !CALL almo_correlation_main(qs_env,almo_scf_env)
146
147      ! do post scf processing
148      CALL almo_scf_post(qs_env, almo_scf_env)
149
150      ! clean up the mess
151      CALL almo_scf_clean_up(almo_scf_env)
152
153      CALL timestop(handle)
154
155   END SUBROUTINE almo_entry_scf
156
157! **************************************************************************************************
158!> \brief Initialization of the almo_scf_env_type.
159!> \param qs_env ...
160!> \param almo_scf_env ...
161!> \param calc_forces ...
162!> \par History
163!>       2011.05 created [Rustam Z Khaliullin]
164!>       2018.09 smearing support [Ruben Staub]
165!> \author Rustam Z Khaliullin
166! **************************************************************************************************
167   SUBROUTINE almo_scf_init(qs_env, almo_scf_env, calc_forces)
168      TYPE(qs_environment_type), POINTER                 :: qs_env
169      TYPE(almo_scf_env_type)                            :: almo_scf_env
170      LOGICAL                                            :: calc_forces
171
172      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_init', routineP = moduleN//':'//routineN
173
174      INTEGER                                            :: ao, handle, i, iao, idomain, ispin, &
175                                                            multip, naos, natoms, ndomains, nelec, &
176                                                            nelec_a, nelec_b, nmols, nspins, &
177                                                            unit_nr
178      TYPE(cp_logger_type), POINTER                      :: logger
179      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
180      TYPE(dft_control_type), POINTER                    :: dft_control
181      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
182      TYPE(section_vals_type), POINTER                   :: input
183
184      CALL timeset(routineN, handle)
185
186      ! define the output_unit
187      logger => cp_get_default_logger()
188      IF (logger%para_env%ionode) THEN
189         unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
190      ELSE
191         unit_nr = -1
192      ENDIF
193
194      ! set optimizers' types
195      almo_scf_env%opt_block_diag_diis%optimizer_type = optimizer_diis
196      almo_scf_env%opt_block_diag_pcg%optimizer_type = optimizer_pcg
197      almo_scf_env%opt_xalmo_diis%optimizer_type = optimizer_diis
198      almo_scf_env%opt_xalmo_pcg%optimizer_type = optimizer_pcg
199      almo_scf_env%opt_xalmo_newton_pcg_solver%optimizer_type = optimizer_lin_eq_pcg
200
201      ! get info from the qs_env
202      CALL get_qs_env(qs_env, &
203                      nelectron_total=almo_scf_env%nelectrons_total, &
204                      matrix_s=matrix_s, &
205                      dft_control=dft_control, &
206                      molecule_set=molecule_set, &
207                      input=input, &
208                      has_unit_metric=almo_scf_env%orthogonal_basis, &
209                      para_env=almo_scf_env%para_env, &
210                      blacs_env=almo_scf_env%blacs_env, &
211                      nelectron_spin=almo_scf_env%nelectrons_spin)
212      CALL cp_para_env_retain(almo_scf_env%para_env)
213      CALL cp_blacs_env_retain(almo_scf_env%blacs_env)
214
215      ! copy basic quantities
216      almo_scf_env%nspins = dft_control%nspins
217      almo_scf_env%natoms = dbcsr_nblkrows_total(matrix_s(1)%matrix)
218      almo_scf_env%nmolecules = SIZE(molecule_set)
219      CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=naos)
220      almo_scf_env%naos = naos
221
222      !! retrieve smearing parameters, and check compatibility of methods requested
223      almo_scf_env%smear = dft_control%smear
224      IF (almo_scf_env%smear) THEN
225         CALL cite_reference(Staub2019)
226         IF ((almo_scf_env%almo_update_algorithm .NE. almo_scf_diag) .OR. &
227             ((almo_scf_env%deloc_method .NE. almo_deloc_none) .AND. &
228              (almo_scf_env%xalmo_update_algorithm .NE. almo_scf_diag))) THEN
229            CPABORT("ALMO smearing is currently implemented for DIAG algorithm only")
230         END IF
231         IF (qs_env%scf_control%smear%method .NE. smear_fermi_dirac) THEN
232            CPABORT("Only Fermi-Dirac smearing is currently compatible with ALMO")
233         END IF
234         almo_scf_env%smear_e_temp = qs_env%scf_control%smear%electronic_temperature
235         IF ((almo_scf_env%mat_distr_aos .NE. almo_mat_distr_molecular) .OR. &
236             (almo_scf_env%domain_layout_mos .NE. almo_domain_layout_molecular)) THEN
237            CPABORT("ALMO smearing was designed to work with molecular fragments only")
238         END IF
239      END IF
240
241      ! convenient local varibales
242      nspins = almo_scf_env%nspins
243      nmols = almo_scf_env%nmolecules
244      natoms = almo_scf_env%natoms
245
246      ! Define groups: either atomic or molecular
247      IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
248         almo_scf_env%ndomains = almo_scf_env%nmolecules
249      ELSE
250         almo_scf_env%ndomains = almo_scf_env%natoms
251      ENDIF
252
253      ! allocate domain descriptors
254      ndomains = almo_scf_env%ndomains
255      ALLOCATE (almo_scf_env%domain_index_of_atom(natoms))
256      ALLOCATE (almo_scf_env%domain_index_of_ao(naos))
257      ALLOCATE (almo_scf_env%first_atom_of_domain(ndomains))
258      ALLOCATE (almo_scf_env%last_atom_of_domain(ndomains))
259      ALLOCATE (almo_scf_env%nbasis_of_domain(ndomains))
260      ALLOCATE (almo_scf_env%nocc_of_domain(ndomains, nspins)) !! with smearing, nb of available orbitals for occupation
261      ALLOCATE (almo_scf_env%real_ne_of_domain(ndomains, nspins)) !! with smearing, nb of fully-occupied orbitals
262      ALLOCATE (almo_scf_env%nvirt_full_of_domain(ndomains, nspins))
263      ALLOCATE (almo_scf_env%nvirt_of_domain(ndomains, nspins))
264      ALLOCATE (almo_scf_env%nvirt_disc_of_domain(ndomains, nspins))
265      ALLOCATE (almo_scf_env%mu_of_domain(ndomains, nspins))
266      ALLOCATE (almo_scf_env%cpu_of_domain(ndomains))
267      ALLOCATE (almo_scf_env%charge_of_domain(ndomains))
268      ALLOCATE (almo_scf_env%multiplicity_of_domain(ndomains))
269
270      ! fill out domain descriptors and group descriptors
271      IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
272         ! get domain info from molecule_set
273         CALL get_molecule_set_info(molecule_set, &
274                                    atom_to_mol=almo_scf_env%domain_index_of_atom, &
275                                    mol_to_first_atom=almo_scf_env%first_atom_of_domain, &
276                                    mol_to_last_atom=almo_scf_env%last_atom_of_domain, &
277                                    mol_to_nelectrons=almo_scf_env%nocc_of_domain(1:ndomains, 1), &
278                                    mol_to_nbasis=almo_scf_env%nbasis_of_domain, &
279                                    mol_to_charge=almo_scf_env%charge_of_domain, &
280                                    mol_to_multiplicity=almo_scf_env%multiplicity_of_domain)
281         ! calculate number of alpha and beta occupied orbitals from
282         ! the number of electrons and multiplicity of each molecule
283         ! Na + Nb = Ne
284         ! Na - Nb = Mult - 1 (assume Na > Nb as we do not have more info from get_molecule_set_info)
285         DO idomain = 1, ndomains
286            nelec = almo_scf_env%nocc_of_domain(idomain, 1)
287            multip = almo_scf_env%multiplicity_of_domain(idomain)
288            nelec_a = (nelec + multip - 1)/2
289            nelec_b = nelec - nelec_a
290            !! Initializing an occupation-rescaling trick if smearing is on
291            IF (almo_scf_env%smear) THEN
292               IF (multip .GT. 1) THEN
293                  CPWARN("BEWARE: Non singlet state detected, treating it as closed-shell")
294               ENDIF
295               !! Save real number of electrons of each spin, as it is required for Fermi-dirac smearing
296               !! BEWARE : Non singlet states are allowed but treated as closed-shell
297               almo_scf_env%real_ne_of_domain(idomain, :) = REAL(nelec, KIND=dp)/2.0_dp
298               !! Add a number of added_mos equal to the number of atoms in domain
299               !! (since fragments were computed this way with smearing)
300               almo_scf_env%nocc_of_domain(idomain, :) = CEILING(almo_scf_env%real_ne_of_domain(idomain, :)) &
301                                                         + (almo_scf_env%last_atom_of_domain(idomain) &
302                                                            - almo_scf_env%first_atom_of_domain(idomain) + 1)
303            ELSE
304               almo_scf_env%nocc_of_domain(idomain, 1) = nelec_a
305               IF (nelec_a .NE. nelec_b) THEN
306                  IF (nspins .EQ. 1) THEN
307                     WRITE (*, *) "Domain ", idomain, " out of ", ndomains, ". Electrons = ", nelec
308                     CPABORT("odd e- -- use unrestricted methods")
309                  ENDIF
310                  almo_scf_env%nocc_of_domain(idomain, 2) = nelec_b
311                  ! RZK-warning: open-shell procedures have not been tested yet
312                  ! Stop the program now
313                  CPABORT("Unrestricted ALMO methods are NYI")
314               ENDIF
315            END IF
316         ENDDO
317         DO ispin = 1, nspins
318            ! take care of the full virtual subspace
319            almo_scf_env%nvirt_full_of_domain(:, ispin) = &
320               almo_scf_env%nbasis_of_domain(:) - &
321               almo_scf_env%nocc_of_domain(:, ispin)
322            ! and the truncated virtual subspace
323            SELECT CASE (almo_scf_env%deloc_truncate_virt)
324            CASE (virt_full)
325               almo_scf_env%nvirt_of_domain(:, ispin) = &
326                  almo_scf_env%nvirt_full_of_domain(:, ispin)
327               almo_scf_env%nvirt_disc_of_domain(:, ispin) = 0
328            CASE (virt_number)
329               DO idomain = 1, ndomains
330                  almo_scf_env%nvirt_of_domain(idomain, ispin) = &
331                     MIN(almo_scf_env%deloc_virt_per_domain, &
332                         almo_scf_env%nvirt_full_of_domain(idomain, ispin))
333                  almo_scf_env%nvirt_disc_of_domain(idomain, ispin) = &
334                     almo_scf_env%nvirt_full_of_domain(idomain, ispin) - &
335                     almo_scf_env%nvirt_of_domain(idomain, ispin)
336               ENDDO
337            CASE (virt_occ_size)
338               DO idomain = 1, ndomains
339                  almo_scf_env%nvirt_of_domain(idomain, ispin) = &
340                     MIN(almo_scf_env%nocc_of_domain(idomain, ispin), &
341                         almo_scf_env%nvirt_full_of_domain(idomain, ispin))
342                  almo_scf_env%nvirt_disc_of_domain(idomain, ispin) = &
343                     almo_scf_env%nvirt_full_of_domain(idomain, ispin) - &
344                     almo_scf_env%nvirt_of_domain(idomain, ispin)
345               ENDDO
346            CASE DEFAULT
347               CPABORT("illegal method for virtual space truncation")
348            END SELECT
349         ENDDO ! spin
350      ELSE ! domains are atomic
351         ! RZK-warning do the same for atomic domains/groups
352         almo_scf_env%domain_index_of_atom(1:natoms) = (/(i, i=1, natoms)/)
353      ENDIF
354
355      ao = 1
356      DO idomain = 1, ndomains
357         DO iao = 1, almo_scf_env%nbasis_of_domain(idomain)
358            almo_scf_env%domain_index_of_ao(ao) = idomain
359            ao = ao + 1
360         ENDDO
361      ENDDO
362
363      almo_scf_env%mu_of_domain(:, :) = almo_scf_env%mu
364
365      ! build domain (i.e. layout) indices for distribution blocks
366      ! ao blocks
367      IF (almo_scf_env%mat_distr_aos == almo_mat_distr_atomic) THEN
368         ALLOCATE (almo_scf_env%domain_index_of_ao_block(natoms))
369         almo_scf_env%domain_index_of_ao_block(:) = &
370            almo_scf_env%domain_index_of_atom(:)
371      ELSE IF (almo_scf_env%mat_distr_aos == almo_mat_distr_molecular) THEN
372         ALLOCATE (almo_scf_env%domain_index_of_ao_block(nmols))
373         ! if distr blocks are molecular then domain layout is also molecular
374         almo_scf_env%domain_index_of_ao_block(:) = (/(i, i=1, nmols)/)
375      ENDIF
376      ! mo blocks
377      IF (almo_scf_env%mat_distr_mos == almo_mat_distr_atomic) THEN
378         ALLOCATE (almo_scf_env%domain_index_of_mo_block(natoms))
379         almo_scf_env%domain_index_of_mo_block(:) = &
380            almo_scf_env%domain_index_of_atom(:)
381      ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN
382         ALLOCATE (almo_scf_env%domain_index_of_mo_block(nmols))
383         ! if distr blocks are molecular then domain layout is also molecular
384         almo_scf_env%domain_index_of_mo_block(:) = (/(i, i=1, nmols)/)
385      ENDIF
386
387      ! set all flags
388      !almo_scf_env%need_previous_ks=.FALSE.
389      !IF (almo_scf_env%deloc_method==almo_deloc_harris) THEN
390      almo_scf_env%need_previous_ks = .TRUE.
391      !ENDIF
392
393      !almo_scf_env%need_virtuals=.FALSE.
394      !almo_scf_env%need_orbital_energies=.FALSE.
395      !IF (almo_scf_env%almo_update_algorithm==almo_scf_diag) THEN
396      almo_scf_env%need_virtuals = .TRUE.
397      almo_scf_env%need_orbital_energies = .TRUE.
398      !ENDIF
399
400      almo_scf_env%calc_forces = calc_forces
401      IF (calc_forces) THEN
402         CALL cite_reference(Scheiber2018)
403         IF (almo_scf_env%deloc_method .EQ. almo_deloc_x .OR. &
404             almo_scf_env%deloc_method .EQ. almo_deloc_xalmo_x .OR. &
405             almo_scf_env%deloc_method .EQ. almo_deloc_xalmo_1diag) THEN
406            CPABORT("Forces for perturbative methods are NYI. Change DELOCALIZE_METHOD")
407         ENDIF
408         ! switch to ASPC after a certain number of exact steps is done
409         IF (almo_scf_env%almo_history%istore .GT. (almo_scf_env%almo_history%nstore + 1)) THEN
410            IF (almo_scf_env%opt_block_diag_pcg%eps_error_early .GT. 0.0_dp) THEN
411               almo_scf_env%opt_block_diag_pcg%eps_error = almo_scf_env%opt_block_diag_pcg%eps_error_early
412               almo_scf_env%opt_block_diag_pcg%early_stopping_on = .TRUE.
413               IF (unit_nr > 0) WRITE (*, *) "ALMO_OPTIMIZER_PCG: EPS_ERROR_EARLY is on"
414            ENDIF
415            IF (almo_scf_env%opt_block_diag_diis%eps_error_early .GT. 0.0_dp) THEN
416               almo_scf_env%opt_block_diag_diis%eps_error = almo_scf_env%opt_block_diag_diis%eps_error_early
417               almo_scf_env%opt_block_diag_diis%early_stopping_on = .TRUE.
418               IF (unit_nr > 0) WRITE (*, *) "ALMO_OPTIMIZER_DIIS: EPS_ERROR_EARLY is on"
419            ENDIF
420            IF (almo_scf_env%opt_block_diag_pcg%max_iter_early .GT. 0) THEN
421               almo_scf_env%opt_block_diag_pcg%max_iter = almo_scf_env%opt_block_diag_pcg%max_iter_early
422               almo_scf_env%opt_block_diag_pcg%early_stopping_on = .TRUE.
423               IF (unit_nr > 0) WRITE (*, *) "ALMO_OPTIMIZER_PCG: MAX_ITER_EARLY is on"
424            ENDIF
425            IF (almo_scf_env%opt_block_diag_diis%max_iter_early .GT. 0) THEN
426               almo_scf_env%opt_block_diag_diis%max_iter = almo_scf_env%opt_block_diag_diis%max_iter_early
427               almo_scf_env%opt_block_diag_diis%early_stopping_on = .TRUE.
428               IF (unit_nr > 0) WRITE (*, *) "ALMO_OPTIMIZER_DIIS: MAX_ITER_EARLY is on"
429            ENDIF
430         ELSE
431            almo_scf_env%opt_block_diag_diis%early_stopping_on = .FALSE.
432            almo_scf_env%opt_block_diag_pcg%early_stopping_on = .FALSE.
433         ENDIF
434         IF (almo_scf_env%xalmo_history%istore .GT. (almo_scf_env%xalmo_history%nstore + 1)) THEN
435            IF (almo_scf_env%opt_xalmo_pcg%eps_error_early .GT. 0.0_dp) THEN
436               almo_scf_env%opt_xalmo_pcg%eps_error = almo_scf_env%opt_xalmo_pcg%eps_error_early
437               almo_scf_env%opt_xalmo_pcg%early_stopping_on = .TRUE.
438               IF (unit_nr > 0) WRITE (*, *) "XALMO_OPTIMIZER_PCG: EPS_ERROR_EARLY is on"
439            ENDIF
440            IF (almo_scf_env%opt_xalmo_pcg%max_iter_early .GT. 0.0_dp) THEN
441               almo_scf_env%opt_xalmo_pcg%max_iter = almo_scf_env%opt_xalmo_pcg%max_iter_early
442               almo_scf_env%opt_xalmo_pcg%early_stopping_on = .TRUE.
443               IF (unit_nr > 0) WRITE (*, *) "XALMO_OPTIMIZER_PCG: MAX_ITER_EARLY is on"
444            ENDIF
445         ELSE
446            almo_scf_env%opt_xalmo_pcg%early_stopping_on = .FALSE.
447         ENDIF
448      ENDIF
449
450      ! create all matrices
451      CALL almo_scf_env_create_matrices(almo_scf_env, matrix_s(1)%matrix)
452
453      ! set up matrix S and all required functions of S
454      almo_scf_env%s_inv_done = .FALSE.
455      almo_scf_env%s_sqrt_done = .FALSE.
456      CALL almo_scf_init_ao_overlap(matrix_s(1)%matrix, almo_scf_env)
457
458      ! create the quencher (imposes sparsity template)
459      CALL almo_scf_construct_quencher(qs_env, almo_scf_env)
460      CALL distribute_domains(almo_scf_env)
461
462      ! FINISH setting job parameters here, print out job info
463      CALL almo_scf_print_job_info(almo_scf_env, unit_nr)
464
465      ! allocate and init the domain preconditioner
466      ALLOCATE (almo_scf_env%domain_preconditioner(ndomains, nspins))
467      CALL init_submatrices(almo_scf_env%domain_preconditioner)
468
469      ! allocate and init projected KS for domains
470      ALLOCATE (almo_scf_env%domain_ks_xx(ndomains, nspins))
471      CALL init_submatrices(almo_scf_env%domain_ks_xx)
472
473      ! init ao-overlap subblocks
474      ALLOCATE (almo_scf_env%domain_s_inv(ndomains, nspins))
475      CALL init_submatrices(almo_scf_env%domain_s_inv)
476      ALLOCATE (almo_scf_env%domain_s_sqrt_inv(ndomains, nspins))
477      CALL init_submatrices(almo_scf_env%domain_s_sqrt_inv)
478      ALLOCATE (almo_scf_env%domain_s_sqrt(ndomains, nspins))
479      CALL init_submatrices(almo_scf_env%domain_s_sqrt)
480      ALLOCATE (almo_scf_env%domain_t(ndomains, nspins))
481      CALL init_submatrices(almo_scf_env%domain_t)
482      ALLOCATE (almo_scf_env%domain_err(ndomains, nspins))
483      CALL init_submatrices(almo_scf_env%domain_err)
484      ALLOCATE (almo_scf_env%domain_r_down_up(ndomains, nspins))
485      CALL init_submatrices(almo_scf_env%domain_r_down_up)
486
487      ! initialization of the KS matrix
488      CALL init_almo_ks_matrix_via_qs(qs_env, &
489                                      almo_scf_env%matrix_ks, &
490                                      almo_scf_env%mat_distr_aos, &
491                                      almo_scf_env%eps_filter)
492      CALL construct_qs_mos(qs_env, almo_scf_env)
493
494      CALL timestop(handle)
495
496   END SUBROUTINE almo_scf_init
497
498! **************************************************************************************************
499!> \brief create the scf initial guess for ALMOs
500!> \param qs_env ...
501!> \param almo_scf_env ...
502!> \par History
503!>       2016.11 created [Rustam Z Khaliullin]
504!>       2018.09 smearing support [Ruben Staub]
505!> \author Rustam Z Khaliullin
506! **************************************************************************************************
507   SUBROUTINE almo_scf_initial_guess(qs_env, almo_scf_env)
508      TYPE(qs_environment_type), POINTER                 :: qs_env
509      TYPE(almo_scf_env_type)                            :: almo_scf_env
510
511      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_initial_guess', &
512         routineP = moduleN//':'//routineN
513
514      CHARACTER(LEN=default_path_length)                 :: file_name, project_name
515      INTEGER                                            :: handle, iaspc, ispin, istore, naspc, &
516                                                            nspins, unit_nr
517      INTEGER, DIMENSION(2)                              :: nelectron_spin
518      LOGICAL                                            :: aspc_guess, has_unit_metric
519      REAL(KIND=dp)                                      :: alpha, cs_pos, energy, kTS_sum
520      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
521      TYPE(cp_logger_type), POINTER                      :: logger
522      TYPE(cp_para_env_type), POINTER                    :: para_env
523      TYPE(dbcsr_distribution_type)                      :: dist
524      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, rho_ao
525      TYPE(dft_control_type), POINTER                    :: dft_control
526      TYPE(molecular_scf_guess_env_type), POINTER        :: mscfg_env
527      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
528      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
529      TYPE(qs_rho_type), POINTER                         :: rho
530
531      CALL timeset(routineN, handle)
532
533      NULLIFY (rho, rho_ao)
534
535      ! get a useful output_unit
536      logger => cp_get_default_logger()
537      IF (logger%para_env%ionode) THEN
538         unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
539      ELSE
540         unit_nr = -1
541      ENDIF
542
543      ! get basic quantities from the qs_env
544      CALL get_qs_env(qs_env, &
545                      dft_control=dft_control, &
546                      matrix_s=matrix_s, &
547                      atomic_kind_set=atomic_kind_set, &
548                      qs_kind_set=qs_kind_set, &
549                      particle_set=particle_set, &
550                      has_unit_metric=has_unit_metric, &
551                      para_env=para_env, &
552                      nelectron_spin=nelectron_spin, &
553                      mscfg_env=mscfg_env, &
554                      rho=rho)
555
556      CALL qs_rho_get(rho, rho_ao=rho_ao)
557      CPASSERT(ASSOCIATED(mscfg_env))
558
559      ! initial guess on the first simulation step is determined by almo_scf_env%almo_scf_guess
560      ! the subsequent simulation steps are determined by extrapolation_order
561      ! if extrapolation order is zero then again almo_scf_env%almo_scf_guess is used
562      ! ... the number of stored history points will remain zero if extrapolation order is zero
563      IF (almo_scf_env%almo_history%istore == 0) THEN
564         aspc_guess = .FALSE.
565      ELSE
566         aspc_guess = .TRUE.
567      ENDIF
568
569      nspins = almo_scf_env%nspins
570
571      ! create an initial guess
572      IF (.NOT. aspc_guess) THEN
573
574         SELECT CASE (almo_scf_env%almo_scf_guess)
575         CASE (molecular_guess)
576
577            DO ispin = 1, nspins
578
579               ! the calculations on "isolated" molecules has already been done
580               ! all we need to do is convert the MOs of molecules into
581               ! the ALMO matrix taking into account different distributions
582               CALL get_matrix_from_submatrices(mscfg_env, &
583                                                almo_scf_env%matrix_t_blk(ispin), ispin)
584               CALL dbcsr_filter(almo_scf_env%matrix_t_blk(ispin), &
585                                 almo_scf_env%eps_filter)
586
587            ENDDO
588
589         CASE (atomic_guess)
590
591            IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical .OR. &
592                dft_control%qs_control%xtb) THEN
593               CALL calculate_mopac_dm(rho_ao, &
594                                       matrix_s(1)%matrix, has_unit_metric, &
595                                       dft_control, particle_set, atomic_kind_set, qs_kind_set, &
596                                       nspins, nelectron_spin, &
597                                       para_env)
598            ELSE
599               CALL calculate_atomic_block_dm(rho_ao, &
600                                              matrix_s(1)%matrix, particle_set, atomic_kind_set, qs_kind_set, &
601                                              nspins, nelectron_spin, unit_nr, para_env)
602            ENDIF
603
604            DO ispin = 1, nspins
605               ! copy the atomic-block dm into matrix_p_blk
606               CALL matrix_qs_to_almo(rho_ao(ispin)%matrix, &
607                                      almo_scf_env%matrix_p_blk(ispin), almo_scf_env%mat_distr_aos, &
608                                      .FALSE.)
609               CALL dbcsr_filter(almo_scf_env%matrix_p_blk(ispin), &
610                                 almo_scf_env%eps_filter)
611            ENDDO ! ispin
612
613            ! obtain orbitals from the density matrix
614            ! (the current version of ALMO SCF needs orbitals)
615            CALL almo_scf_p_blk_to_t_blk(almo_scf_env, ionic=.FALSE.)
616
617         CASE (restart_guess)
618
619            project_name = logger%iter_info%project_name
620
621            DO ispin = 1, nspins
622               WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_ALMO_SPIN_", ispin, "_RESTART.mo"
623               CALL dbcsr_get_info(almo_scf_env%matrix_t_blk(ispin), distribution=dist)
624               CALL dbcsr_binary_read(file_name, distribution=dist, matrix_new=almo_scf_env%matrix_t_blk(ispin))
625               cs_pos = dbcsr_checksum(almo_scf_env%matrix_t_blk(ispin), pos=.TRUE.)
626               IF (unit_nr > 0) THEN
627                  WRITE (unit_nr, '(T2,A,E20.8)') "Read restart ALMO "//TRIM(file_name)//" with checksum: ", cs_pos
628               ENDIF
629            ENDDO
630         END SELECT
631
632      ELSE !aspc_guess
633
634         CALL cite_reference(Kolafa2004)
635
636         naspc = MIN(almo_scf_env%almo_history%istore, almo_scf_env%almo_history%nstore)
637         IF (unit_nr > 0) THEN
638            WRITE (unit_nr, FMT="(/,T2,A,/,/,T3,A,I0)") &
639               "Parameters for the always stable predictor-corrector (ASPC) method:", &
640               "ASPC order: ", naspc
641         END IF
642
643         DO ispin = 1, nspins
644
645            ! extrapolation
646            DO iaspc = 1, naspc
647               istore = MOD(almo_scf_env%almo_history%istore - iaspc, almo_scf_env%almo_history%nstore) + 1
648               alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
649                       binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
650               IF (unit_nr > 0) THEN
651                  WRITE (unit_nr, FMT="(T3,A2,I0,A4,F10.6)") &
652                     "B(", iaspc, ") = ", alpha
653               END IF
654               IF (iaspc == 1) THEN
655                  CALL dbcsr_copy(almo_scf_env%matrix_t_blk(ispin), &
656                                  almo_scf_env%almo_history%matrix_t(ispin), &
657                                  keep_sparsity=.TRUE.)
658                  CALL dbcsr_scale(almo_scf_env%matrix_t_blk(ispin), alpha)
659               ELSE
660                  CALL dbcsr_multiply("N", "N", alpha, &
661                                      almo_scf_env%almo_history%matrix_p_up_down(ispin, istore), &
662                                      almo_scf_env%almo_history%matrix_t(ispin), &
663                                      1.0_dp, almo_scf_env%matrix_t_blk(ispin), &
664                                      retain_sparsity=.TRUE.)
665               ENDIF
666            ENDDO !iaspc
667
668         ENDDO !ispin
669
670      ENDIF !aspc_guess?
671
672      DO ispin = 1, nspins
673
674         CALL orthogonalize_mos(ket=almo_scf_env%matrix_t_blk(ispin), &
675                                overlap=almo_scf_env%matrix_sigma_blk(ispin), &
676                                metric=almo_scf_env%matrix_s_blk(1), &
677                                retain_locality=.TRUE., &
678                                only_normalize=.FALSE., &
679                                nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
680                                eps_filter=almo_scf_env%eps_filter, &
681                                order_lanczos=almo_scf_env%order_lanczos, &
682                                eps_lanczos=almo_scf_env%eps_lanczos, &
683                                max_iter_lanczos=almo_scf_env%max_iter_lanczos)
684
685         !! Application of an occupation-rescaling trick for smearing, if requested
686         IF (almo_scf_env%smear) THEN
687            CALL almo_scf_t_rescaling(matrix_t=almo_scf_env%matrix_t_blk(ispin), &
688                                      mo_energies=almo_scf_env%mo_energies(:, ispin), &
689                                      mu_of_domain=almo_scf_env%mu_of_domain(:, ispin), &
690                                      real_ne_of_domain=almo_scf_env%real_ne_of_domain(:, ispin), &
691                                      spin_kTS=almo_scf_env%kTS(ispin), &
692                                      smear_e_temp=almo_scf_env%smear_e_temp, &
693                                      ndomains=almo_scf_env%ndomains, &
694                                      nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin))
695         END IF
696
697         CALL almo_scf_t_to_proj(t=almo_scf_env%matrix_t_blk(ispin), &
698                                 p=almo_scf_env%matrix_p(ispin), &
699                                 eps_filter=almo_scf_env%eps_filter, &
700                                 orthog_orbs=.FALSE., &
701                                 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
702                                 s=almo_scf_env%matrix_s(1), &
703                                 sigma=almo_scf_env%matrix_sigma(ispin), &
704                                 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
705                                 use_guess=.FALSE., &
706                                 smear=almo_scf_env%smear, &
707                                 algorithm=almo_scf_env%sigma_inv_algorithm, &
708                                 eps_lanczos=almo_scf_env%eps_lanczos, &
709                                 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
710                                 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
711                                 para_env=almo_scf_env%para_env, &
712                                 blacs_env=almo_scf_env%blacs_env)
713
714      ENDDO
715
716      ! compute dm from the projector(s)
717      IF (nspins == 1) THEN
718         CALL dbcsr_scale(almo_scf_env%matrix_p(1), 2.0_dp)
719         !! Rescaling electronic entropy contribution by spin_factor
720         IF (almo_scf_env%smear) THEN
721            almo_scf_env%kTS(1) = almo_scf_env%kTS(1)*2.0_dp
722         END IF
723      ENDIF
724
725      IF (almo_scf_env%smear) THEN
726         kTS_sum = SUM(almo_scf_env%kTS)
727      ELSE
728         kTS_sum = 0.0_dp
729      ENDIF
730
731      CALL almo_dm_to_almo_ks(qs_env, &
732                              almo_scf_env%matrix_p, &
733                              almo_scf_env%matrix_ks, &
734                              energy, &
735                              almo_scf_env%eps_filter, &
736                              almo_scf_env%mat_distr_aos, &
737                              smear=almo_scf_env%smear, &
738                              kTS_sum=kTS_sum)
739
740      IF (unit_nr > 0) THEN
741         IF (almo_scf_env%almo_scf_guess .EQ. molecular_guess) THEN
742            WRITE (unit_nr, '(T2,A38,F40.10)') "Single-molecule energy:", &
743               SUM(mscfg_env%energy_of_frag)
744         ENDIF
745         WRITE (unit_nr, '(T2,A38,F40.10)') "Energy of the initial guess:", energy
746         WRITE (unit_nr, '()')
747      ENDIF
748
749      CALL timestop(handle)
750
751   END SUBROUTINE almo_scf_initial_guess
752
753! **************************************************************************************************
754!> \brief store a history of matrices for later use in almo_scf_initial_guess
755!> \param almo_scf_env ...
756!> \par History
757!>       2016.11 created [Rustam Z Khaliullin]
758!> \author Rustam Khaliullin
759! **************************************************************************************************
760   SUBROUTINE almo_scf_store_extrapolation_data(almo_scf_env)
761      TYPE(almo_scf_env_type)                            :: almo_scf_env
762
763      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_store_extrapolation_data', &
764         routineP = moduleN//':'//routineN
765
766      INTEGER                                            :: handle, ispin, istore, unit_nr
767      LOGICAL :: delocalization_uses_extrapolation
768      TYPE(cp_logger_type), POINTER                      :: logger
769      TYPE(dbcsr_type)                                   :: matrix_no_tmp1, matrix_no_tmp2, &
770                                                            matrix_no_tmp3, matrix_no_tmp4
771
772      CALL timeset(routineN, handle)
773
774      ! get a useful output_unit
775      logger => cp_get_default_logger()
776      IF (logger%para_env%ionode) THEN
777         unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
778      ELSE
779         unit_nr = -1
780      ENDIF
781
782      IF (almo_scf_env%almo_history%nstore > 0) THEN
783
784         almo_scf_env%almo_history%istore = almo_scf_env%almo_history%istore + 1
785
786         DO ispin = 1, SIZE(almo_scf_env%matrix_t_blk)
787
788            istore = MOD(almo_scf_env%almo_history%istore - 1, almo_scf_env%almo_history%nstore) + 1
789
790            IF (almo_scf_env%almo_history%istore == 1) THEN
791               CALL dbcsr_create(almo_scf_env%almo_history%matrix_t(ispin), &
792                                 template=almo_scf_env%matrix_t_blk(ispin), &
793                                 matrix_type=dbcsr_type_no_symmetry)
794            ENDIF
795            CALL dbcsr_copy(almo_scf_env%almo_history%matrix_t(ispin), &
796                            almo_scf_env%matrix_t_blk(ispin))
797
798            IF (almo_scf_env%almo_history%istore <= almo_scf_env%almo_history%nstore) THEN
799               CALL dbcsr_create(almo_scf_env%almo_history%matrix_p_up_down(ispin, istore), &
800                                 template=almo_scf_env%matrix_s(1), &
801                                 matrix_type=dbcsr_type_no_symmetry)
802            ENDIF
803
804            CALL dbcsr_create(matrix_no_tmp1, template=almo_scf_env%matrix_t_blk(ispin), &
805                              matrix_type=dbcsr_type_no_symmetry)
806            CALL dbcsr_create(matrix_no_tmp2, template=almo_scf_env%matrix_t_blk(ispin), &
807                              matrix_type=dbcsr_type_no_symmetry)
808
809            ! compute contra-covariant density matrix
810            CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s(1), &
811                                almo_scf_env%matrix_t_blk(ispin), &
812                                0.0_dp, matrix_no_tmp1, &
813                                filter_eps=almo_scf_env%eps_filter)
814            CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_no_tmp1, &
815                                almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
816                                0.0_dp, matrix_no_tmp2, &
817                                filter_eps=almo_scf_env%eps_filter)
818            CALL dbcsr_multiply("N", "T", 1.0_dp, &
819                                almo_scf_env%matrix_t_blk(ispin), &
820                                matrix_no_tmp2, &
821                                0.0_dp, almo_scf_env%almo_history%matrix_p_up_down(ispin, istore), &
822                                filter_eps=almo_scf_env%eps_filter)
823
824            CALL dbcsr_release(matrix_no_tmp1)
825            CALL dbcsr_release(matrix_no_tmp2)
826
827         ENDDO
828
829      ENDIF
830
831      ! exrapolate xalmos?
832      delocalization_uses_extrapolation = &
833         .NOT. ((almo_scf_env%deloc_method .EQ. almo_deloc_none) .OR. &
834                (almo_scf_env%deloc_method .EQ. almo_deloc_xalmo_1diag))
835      IF (almo_scf_env%xalmo_history%nstore > 0 .AND. &
836          delocalization_uses_extrapolation) THEN
837
838         almo_scf_env%xalmo_history%istore = almo_scf_env%xalmo_history%istore + 1
839
840         DO ispin = 1, SIZE(almo_scf_env%matrix_t)
841
842            istore = MOD(almo_scf_env%xalmo_history%istore - 1, almo_scf_env%xalmo_history%nstore) + 1
843
844            IF (almo_scf_env%xalmo_history%istore == 1) THEN
845               CALL dbcsr_create(almo_scf_env%xalmo_history%matrix_t(ispin), &
846                                 template=almo_scf_env%matrix_t(ispin), &
847                                 matrix_type=dbcsr_type_no_symmetry)
848            ENDIF
849            CALL dbcsr_copy(almo_scf_env%xalmo_history%matrix_t(ispin), &
850                            almo_scf_env%matrix_t(ispin))
851
852            IF (almo_scf_env%xalmo_history%istore <= almo_scf_env%xalmo_history%nstore) THEN
853               !CALL dbcsr_init(almo_scf_env%xalmo_history%matrix_x(ispin, istore))
854               !CALL dbcsr_create(almo_scf_env%xalmo_history%matrix_x(ispin, istore), &
855               !        template=almo_scf_env%matrix_t(ispin), &
856               !        matrix_type=dbcsr_type_no_symmetry)
857               CALL dbcsr_create(almo_scf_env%xalmo_history%matrix_p_up_down(ispin, istore), &
858                                 template=almo_scf_env%matrix_s(1), &
859                                 matrix_type=dbcsr_type_no_symmetry)
860            ENDIF
861
862            CALL dbcsr_create(matrix_no_tmp3, template=almo_scf_env%matrix_t(ispin), &
863                              matrix_type=dbcsr_type_no_symmetry)
864            CALL dbcsr_create(matrix_no_tmp4, template=almo_scf_env%matrix_t(ispin), &
865                              matrix_type=dbcsr_type_no_symmetry)
866
867            ! compute contra-covariant density matrix
868            CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s(1), &
869                                almo_scf_env%matrix_t(ispin), &
870                                0.0_dp, matrix_no_tmp3, &
871                                filter_eps=almo_scf_env%eps_filter)
872            CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_no_tmp3, &
873                                almo_scf_env%matrix_sigma_inv(ispin), &
874                                0.0_dp, matrix_no_tmp4, &
875                                filter_eps=almo_scf_env%eps_filter)
876            CALL dbcsr_multiply("N", "T", 1.0_dp, &
877                                almo_scf_env%matrix_t(ispin), &
878                                matrix_no_tmp4, &
879                                0.0_dp, almo_scf_env%xalmo_history%matrix_p_up_down(ispin, istore), &
880                                filter_eps=almo_scf_env%eps_filter)
881
882            ! store the difference between t and t0
883            !CALL dbcsr_copy(almo_scf_env%xalmo_history%matrix_x(ispin, istore),&
884            !        almo_scf_env%matrix_t(ispin))
885            !CALL dbcsr_add(almo_scf_env%xalmo_history%matrix_x(ispin, istore),&
886            !        almo_scf_env%matrix_t_blk(ispin),1.0_dp,-1.0_dp)
887
888            CALL dbcsr_release(matrix_no_tmp3)
889            CALL dbcsr_release(matrix_no_tmp4)
890
891         ENDDO
892
893      ENDIF
894
895      CALL timestop(handle)
896
897   END SUBROUTINE almo_scf_store_extrapolation_data
898
899! **************************************************************************************************
900!> \brief Prints out a short summary about the ALMO SCF job
901!> \param almo_scf_env ...
902!> \param unit_nr ...
903!> \par History
904!>       2011.10 created [Rustam Z Khaliullin]
905!> \author Rustam Z Khaliullin
906! **************************************************************************************************
907   SUBROUTINE almo_scf_print_job_info(almo_scf_env, unit_nr)
908
909      TYPE(almo_scf_env_type), INTENT(INOUT)             :: almo_scf_env
910      INTEGER, INTENT(IN)                                :: unit_nr
911
912      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_print_job_info', &
913         routineP = moduleN//':'//routineN
914
915      CHARACTER(len=13)                                  :: neig_string
916      CHARACTER(len=33)                                  :: deloc_method_string
917      INTEGER                                            :: handle, idomain, index1_prev, sum_temp
918      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: nneighbors
919
920      CALL timeset(routineN, handle)
921
922      IF (unit_nr > 0) THEN
923         WRITE (unit_nr, '()')
924         WRITE (unit_nr, '(T2,A,A,A)') REPEAT("-", 32), " ALMO SETTINGS ", REPEAT("-", 32)
925
926         WRITE (unit_nr, '(T2,A,T48,E33.3)') "eps_filter:", almo_scf_env%eps_filter
927
928         IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_skip) THEN
929            WRITE (unit_nr, '(T2,A)') "skip optimization of block-diagonal ALMOs"
930         ELSE
931            WRITE (unit_nr, '(T2,A)') "optimization of block-diagonal ALMOs:"
932            SELECT CASE (almo_scf_env%almo_update_algorithm)
933            CASE (almo_scf_diag)
934               ! the DIIS algorith is the only choice for the diagonlaization-based algorithm
935               CALL print_optimizer_options(almo_scf_env%opt_block_diag_diis, unit_nr)
936            CASE (almo_scf_pcg)
937               ! print out PCG options
938               CALL print_optimizer_options(almo_scf_env%opt_block_diag_pcg, unit_nr)
939            END SELECT
940         ENDIF
941
942         SELECT CASE (almo_scf_env%deloc_method)
943         CASE (almo_deloc_none)
944            deloc_method_string = "NONE"
945         CASE (almo_deloc_x)
946            deloc_method_string = "FULL_X"
947         CASE (almo_deloc_scf)
948            deloc_method_string = "FULL_SCF"
949         CASE (almo_deloc_x_then_scf)
950            deloc_method_string = "FULL_X_THEN_SCF"
951         CASE (almo_deloc_xalmo_1diag)
952            deloc_method_string = "XALMO_1DIAG"
953         CASE (almo_deloc_xalmo_x)
954            deloc_method_string = "XALMO_X"
955         CASE (almo_deloc_xalmo_scf)
956            deloc_method_string = "XALMO_SCF"
957         END SELECT
958         WRITE (unit_nr, '(T2,A,T48,A33)') "delocalization:", TRIM(deloc_method_string)
959
960         IF (almo_scf_env%deloc_method .NE. almo_deloc_none) THEN
961
962            SELECT CASE (almo_scf_env%deloc_method)
963            CASE (almo_deloc_x, almo_deloc_scf, almo_deloc_x_then_scf)
964               WRITE (unit_nr, '(T2,A,T48,A33)') "delocalization cutoff radius:", &
965                  "infinite"
966               deloc_method_string = "FULL_X_THEN_SCF"
967            CASE (almo_deloc_xalmo_1diag, almo_deloc_xalmo_x, almo_deloc_xalmo_scf)
968               WRITE (unit_nr, '(T2,A,T48,F33.5)') "XALMO cutoff radius:", &
969                  almo_scf_env%quencher_r0_factor
970            END SELECT
971
972            IF (almo_scf_env%deloc_method .EQ. almo_deloc_xalmo_1diag) THEN
973               ! print nothing because no actual optimization is done
974            ELSE
975               WRITE (unit_nr, '(T2,A)') "optimization of extended orbitals:"
976               SELECT CASE (almo_scf_env%xalmo_update_algorithm)
977               CASE (almo_scf_diag)
978                  CALL print_optimizer_options(almo_scf_env%opt_xalmo_diis, unit_nr)
979               CASE (almo_scf_pcg)
980                  CALL print_optimizer_options(almo_scf_env%opt_xalmo_pcg, unit_nr)
981               END SELECT
982            ENDIF
983
984         ENDIF
985
986         !SELECT CASE(almo_scf_env%domain_layout_mos)
987         !CASE(almo_domain_layout_orbital)
988         !    WRITE(unit_nr,'(T2,A,T48,A33)') "Delocalization constraints","ORBITAL"
989         !CASE(almo_domain_layout_atomic)
990         !    WRITE(unit_nr,'(T2,A,T48,A33)') "Delocalization constraints","ATOMIC"
991         !CASE(almo_domain_layout_molecular)
992         !    WRITE(unit_nr,'(T2,A,T48,A33)') "Delocalization constraints","MOLECULAR"
993         !END SELECT
994
995         !SELECT CASE(almo_scf_env%domain_layout_aos)
996         !CASE(almo_domain_layout_atomic)
997         !    WRITE(unit_nr,'(T2,A,T48,A33)') "Basis function domains","ATOMIC"
998         !CASE(almo_domain_layout_molecular)
999         !    WRITE(unit_nr,'(T2,A,T48,A33)') "Basis function domains","MOLECULAR"
1000         !END SELECT
1001
1002         !SELECT CASE(almo_scf_env%mat_distr_aos)
1003         !CASE(almo_mat_distr_atomic)
1004         !    WRITE(unit_nr,'(T2,A,T48,A33)') "Parallel distribution for AOs","ATOMIC"
1005         !CASE(almo_mat_distr_molecular)
1006         !    WRITE(unit_nr,'(T2,A,T48,A33)') "Parallel distribution for AOs","MOLECULAR"
1007         !END SELECT
1008
1009         !SELECT CASE(almo_scf_env%mat_distr_mos)
1010         !CASE(almo_mat_distr_atomic)
1011         !    WRITE(unit_nr,'(T2,A,T48,A33)') "Parallel distribution for MOs","ATOMIC"
1012         !CASE(almo_mat_distr_molecular)
1013         !    WRITE(unit_nr,'(T2,A,T48,A33)') "Parallel distribution for MOs","MOLECULAR"
1014         !END SELECT
1015
1016         ! print fragment's statistics
1017         WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
1018         WRITE (unit_nr, '(T2,A,T48,I33)') "Total fragments:", &
1019            almo_scf_env%ndomains
1020
1021         sum_temp = SUM(almo_scf_env%nbasis_of_domain(:))
1022         WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1023            "Basis set size per fragment (min, av, max, total):", &
1024            MINVAL(almo_scf_env%nbasis_of_domain(:)), &
1025            (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1026            MAXVAL(almo_scf_env%nbasis_of_domain(:)), &
1027            sum_temp
1028         !WRITE (unit_nr, '(T2,I13,F13.3,I13,I13)') &
1029         !         MINVAL(almo_scf_env%nbasis_of_domain(:)), &
1030         !         (1.0_dp*sum_temp) / almo_scf_env%ndomains, &
1031         !         MAXVAL(almo_scf_env%nbasis_of_domain(:)), &
1032         !         sum_temp
1033
1034         sum_temp = SUM(almo_scf_env%nocc_of_domain(:, :))
1035         WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1036            "Occupied MOs per fragment (min, av, max, total):", &
1037            MINVAL(SUM(almo_scf_env%nocc_of_domain, DIM=2)), &
1038            (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1039            MAXVAL(SUM(almo_scf_env%nocc_of_domain, DIM=2)), &
1040            sum_temp
1041         !WRITE (unit_nr, '(T2,I13,F13.3,I13,I13)') &
1042         !         MINVAL( SUM(almo_scf_env%nocc_of_domain, DIM=2) ), &
1043         !         (1.0_dp*sum_temp) / almo_scf_env%ndomains, &
1044         !         MAXVAL( SUM(almo_scf_env%nocc_of_domain, DIM=2) ), &
1045         !         sum_temp
1046
1047         sum_temp = SUM(almo_scf_env%nvirt_of_domain(:, :))
1048         WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1049            "Virtual MOs per fragment (min, av, max, total):", &
1050            MINVAL(SUM(almo_scf_env%nvirt_of_domain, DIM=2)), &
1051            (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1052            MAXVAL(SUM(almo_scf_env%nvirt_of_domain, DIM=2)), &
1053            sum_temp
1054         !WRITE (unit_nr, '(T2,I13,F13.3,I13,I13)') &
1055         !         MINVAL( SUM(almo_scf_env%nvirt_of_domain, DIM=2) ), &
1056         !         (1.0_dp*sum_temp) / almo_scf_env%ndomains, &
1057         !         MAXVAL( SUM(almo_scf_env%nvirt_of_domain, DIM=2) ), &
1058         !         sum_temp
1059
1060         sum_temp = SUM(almo_scf_env%charge_of_domain(:))
1061         WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1062            "Charges per fragment (min, av, max, total):", &
1063            MINVAL(almo_scf_env%charge_of_domain(:)), &
1064            (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1065            MAXVAL(almo_scf_env%charge_of_domain(:)), &
1066            sum_temp
1067         !WRITE (unit_nr, '(T2,I13,F13.3,I13,I13)') &
1068         !         MINVAL(almo_scf_env%charge_of_domain(:)), &
1069         !         (1.0_dp*sum_temp) / almo_scf_env%ndomains, &
1070         !         MAXVAL(almo_scf_env%charge_of_domain(:)), &
1071         !         sum_temp
1072
1073         ! compute the number of neighbors of each fragment
1074         ALLOCATE (nneighbors(almo_scf_env%ndomains))
1075
1076         DO idomain = 1, almo_scf_env%ndomains
1077
1078            IF (idomain .EQ. 1) THEN
1079               index1_prev = 1
1080            ELSE
1081               index1_prev = almo_scf_env%domain_map(1)%index1(idomain - 1)
1082            ENDIF
1083
1084            SELECT CASE (almo_scf_env%deloc_method)
1085            CASE (almo_deloc_none)
1086               nneighbors(idomain) = 0
1087            CASE (almo_deloc_x, almo_deloc_scf, almo_deloc_x_then_scf)
1088               nneighbors(idomain) = almo_scf_env%ndomains - 1 ! minus self
1089            CASE (almo_deloc_xalmo_1diag, almo_deloc_xalmo_x, almo_deloc_xalmo_scf)
1090               nneighbors(idomain) = almo_scf_env%domain_map(1)%index1(idomain) - index1_prev - 1 ! minus self
1091            CASE DEFAULT
1092               nneighbors(idomain) = -1
1093            END SELECT
1094
1095         ENDDO ! cycle over domains
1096
1097         sum_temp = SUM(nneighbors(:))
1098         WRITE (unit_nr, '(T2,A,T53,I5,F9.2,I5,I9)') &
1099            "Deloc. neighbors of fragment (min, av, max, total):", &
1100            MINVAL(nneighbors(:)), &
1101            (1.0_dp*sum_temp)/almo_scf_env%ndomains, &
1102            MAXVAL(nneighbors(:)), &
1103            sum_temp
1104
1105         WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
1106         WRITE (unit_nr, '()')
1107
1108         IF (almo_scf_env%ndomains .LE. 64) THEN
1109
1110            ! print fragment info
1111            WRITE (unit_nr, '(T2,A10,A13,A13,A13,A13,A13)') &
1112               "Fragment", "Basis Set", "Occupied", "Virtual", "Charge", "Deloc Neig" !,"Discarded Virt"
1113            WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
1114            DO idomain = 1, almo_scf_env%ndomains
1115
1116               SELECT CASE (almo_scf_env%deloc_method)
1117               CASE (almo_deloc_none)
1118                  neig_string = "NONE"
1119               CASE (almo_deloc_x, almo_deloc_scf, almo_deloc_x_then_scf)
1120                  neig_string = "ALL"
1121               CASE (almo_deloc_xalmo_1diag, almo_deloc_xalmo_x, almo_deloc_xalmo_scf)
1122                  WRITE (neig_string, '(I13)') nneighbors(idomain)
1123               CASE DEFAULT
1124                  neig_string = "N/A"
1125               END SELECT
1126
1127               WRITE (unit_nr, '(T2,I10,I13,I13,I13,I13,A13)') &
1128                  idomain, almo_scf_env%nbasis_of_domain(idomain), &
1129                  SUM(almo_scf_env%nocc_of_domain(idomain, :)), &
1130                  SUM(almo_scf_env%nvirt_of_domain(idomain, :)), &
1131                  !SUM(almo_scf_env%nvirt_disc_of_domain(idomain,:)),&
1132                  almo_scf_env%charge_of_domain(idomain), &
1133                  ADJUSTR(TRIM(neig_string))
1134
1135            ENDDO ! cycle over domains
1136
1137            SELECT CASE (almo_scf_env%deloc_method)
1138            CASE (almo_deloc_xalmo_1diag, almo_deloc_xalmo_x, almo_deloc_xalmo_scf)
1139
1140               WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
1141
1142               ! print fragment neigbors
1143               WRITE (unit_nr, '(T2,A78)') &
1144                  "Neighbor lists (including self)"
1145               WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
1146               DO idomain = 1, almo_scf_env%ndomains
1147
1148                  IF (idomain .EQ. 1) THEN
1149                     index1_prev = 1
1150                  ELSE
1151                     index1_prev = almo_scf_env%domain_map(1)%index1(idomain - 1)
1152                  ENDIF
1153
1154                  WRITE (unit_nr, '(T2,I10,":")') idomain
1155                  WRITE (unit_nr, '(T12,11I6)') &
1156                     almo_scf_env%domain_map(1)%pairs &
1157                     (index1_prev:almo_scf_env%domain_map(1)%index1(idomain) - 1, 1) ! includes self
1158
1159               ENDDO ! cycle over domains
1160
1161            END SELECT
1162
1163         ELSE ! too big to print details for each fragment
1164
1165            WRITE (unit_nr, '(T2,A)') "The system is too big to print details for each fragment."
1166
1167         ENDIF ! how many fragments?
1168
1169         WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
1170
1171         WRITE (unit_nr, '()')
1172
1173         DEALLOCATE (nneighbors)
1174
1175      ENDIF ! unit_nr > 0
1176
1177      CALL timestop(handle)
1178
1179   END SUBROUTINE almo_scf_print_job_info
1180
1181! **************************************************************************************************
1182!> \brief Initializes the ALMO SCF copy of the AO overlap matrix
1183!>        and all necessary functions (sqrt, inverse...)
1184!> \param matrix_s ...
1185!> \param almo_scf_env ...
1186!> \par History
1187!>       2011.06 created [Rustam Z Khaliullin]
1188!> \author Rustam Z Khaliullin
1189! **************************************************************************************************
1190   SUBROUTINE almo_scf_init_ao_overlap(matrix_s, almo_scf_env)
1191      TYPE(dbcsr_type), INTENT(IN)                       :: matrix_s
1192      TYPE(almo_scf_env_type), INTENT(INOUT)             :: almo_scf_env
1193
1194      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_init_ao_overlap', &
1195         routineP = moduleN//':'//routineN
1196
1197      INTEGER                                            :: handle, unit_nr
1198      TYPE(cp_logger_type), POINTER                      :: logger
1199
1200      CALL timeset(routineN, handle)
1201
1202      ! get a useful output_unit
1203      logger => cp_get_default_logger()
1204      IF (logger%para_env%ionode) THEN
1205         unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1206      ELSE
1207         unit_nr = -1
1208      ENDIF
1209
1210      ! make almo copy of S
1211      ! also copy S to S_blk (i.e. to S with the domain structure imposed)
1212      IF (almo_scf_env%orthogonal_basis) THEN
1213         CALL dbcsr_set(almo_scf_env%matrix_s(1), 0.0_dp)
1214         CALL dbcsr_add_on_diag(almo_scf_env%matrix_s(1), 1.0_dp)
1215         CALL dbcsr_set(almo_scf_env%matrix_s_blk(1), 0.0_dp)
1216         CALL dbcsr_add_on_diag(almo_scf_env%matrix_s_blk(1), 1.0_dp)
1217      ELSE
1218         CALL matrix_qs_to_almo(matrix_s, almo_scf_env%matrix_s(1), &
1219                                almo_scf_env%mat_distr_aos, .FALSE.)
1220         CALL dbcsr_copy(almo_scf_env%matrix_s_blk(1), &
1221                         almo_scf_env%matrix_s(1), keep_sparsity=.TRUE.)
1222      ENDIF
1223
1224      CALL dbcsr_filter(almo_scf_env%matrix_s(1), almo_scf_env%eps_filter)
1225      CALL dbcsr_filter(almo_scf_env%matrix_s_blk(1), almo_scf_env%eps_filter)
1226
1227      IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_diag) THEN
1228         CALL matrix_sqrt_Newton_Schulz(almo_scf_env%matrix_s_blk_sqrt(1), &
1229                                        almo_scf_env%matrix_s_blk_sqrt_inv(1), &
1230                                        almo_scf_env%matrix_s_blk(1), &
1231                                        threshold=almo_scf_env%eps_filter, &
1232                                        order=almo_scf_env%order_lanczos, &
1233                                        !order=0, &
1234                                        eps_lanczos=almo_scf_env%eps_lanczos, &
1235                                        max_iter_lanczos=almo_scf_env%max_iter_lanczos)
1236      ELSE IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_dm_sign) THEN
1237         CALL invert_Hotelling(almo_scf_env%matrix_s_blk_inv(1), &
1238                               almo_scf_env%matrix_s_blk(1), &
1239                               threshold=almo_scf_env%eps_filter, &
1240                               filter_eps=almo_scf_env%eps_filter)
1241      ENDIF
1242
1243      CALL timestop(handle)
1244
1245   END SUBROUTINE almo_scf_init_ao_overlap
1246
1247! **************************************************************************************************
1248!> \brief Selects the subroutine for the optimization of block-daigonal ALMOs.
1249!>        Keep it short and clean.
1250!> \param qs_env ...
1251!> \param almo_scf_env ...
1252!> \par History
1253!>       2011.11 created [Rustam Z Khaliullin]
1254!> \author Rustam Z Khaliullin
1255! **************************************************************************************************
1256   SUBROUTINE almo_scf_main(qs_env, almo_scf_env)
1257      TYPE(qs_environment_type), POINTER                 :: qs_env
1258      TYPE(almo_scf_env_type)                            :: almo_scf_env
1259
1260      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_main', routineP = moduleN//':'//routineN
1261
1262      INTEGER                                            :: handle, ispin, unit_nr
1263      TYPE(cp_logger_type), POINTER                      :: logger
1264
1265      CALL timeset(routineN, handle)
1266
1267      ! get a useful output_unit
1268      logger => cp_get_default_logger()
1269      IF (logger%para_env%ionode) THEN
1270         unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1271      ELSE
1272         unit_nr = -1
1273      ENDIF
1274
1275      SELECT CASE (almo_scf_env%almo_update_algorithm)
1276      CASE (almo_scf_pcg)
1277
1278         ! ALMO PCG optimizer as a special case of XALMO PCG
1279         CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1280                                 almo_scf_env=almo_scf_env, &
1281                                 optimizer=almo_scf_env%opt_block_diag_pcg, &
1282                                 quench_t=almo_scf_env%quench_t_blk, &
1283                                 matrix_t_in=almo_scf_env%matrix_t_blk, &
1284                                 matrix_t_out=almo_scf_env%matrix_t_blk, &
1285                                 assume_t0_q0x=.FALSE., &
1286                                 perturbation_only=.FALSE., &
1287                                 special_case=xalmo_case_block_diag)
1288
1289         DO ispin = 1, almo_scf_env%nspins
1290            CALL orthogonalize_mos(ket=almo_scf_env%matrix_t_blk(ispin), &
1291                                   overlap=almo_scf_env%matrix_sigma_blk(ispin), &
1292                                   metric=almo_scf_env%matrix_s_blk(1), &
1293                                   retain_locality=.TRUE., &
1294                                   only_normalize=.FALSE., &
1295                                   nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1296                                   eps_filter=almo_scf_env%eps_filter, &
1297                                   order_lanczos=almo_scf_env%order_lanczos, &
1298                                   eps_lanczos=almo_scf_env%eps_lanczos, &
1299                                   max_iter_lanczos=almo_scf_env%max_iter_lanczos)
1300         ENDDO
1301         !CALL almo_scf_t_blk_to_t_blk_orthonormal(almo_scf_env)
1302
1303      CASE (almo_scf_diag)
1304
1305         ! mixing/DIIS optimizer
1306         CALL almo_scf_block_diagonal(qs_env, almo_scf_env, &
1307                                      almo_scf_env%opt_block_diag_diis)
1308
1309      CASE (almo_scf_skip)
1310
1311         DO ispin = 1, almo_scf_env%nspins
1312            CALL orthogonalize_mos(ket=almo_scf_env%matrix_t_blk(ispin), &
1313                                   overlap=almo_scf_env%matrix_sigma_blk(ispin), &
1314                                   metric=almo_scf_env%matrix_s_blk(1), &
1315                                   retain_locality=.TRUE., &
1316                                   only_normalize=.FALSE., &
1317                                   nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1318                                   eps_filter=almo_scf_env%eps_filter, &
1319                                   order_lanczos=almo_scf_env%order_lanczos, &
1320                                   eps_lanczos=almo_scf_env%eps_lanczos, &
1321                                   max_iter_lanczos=almo_scf_env%max_iter_lanczos)
1322         ENDDO
1323         !CALL almo_scf_t_blk_to_t_blk_orthonormal(almo_scf_env)
1324
1325      END SELECT
1326
1327      ! we might need a copy of the converged KS and sigma_inv
1328      DO ispin = 1, almo_scf_env%nspins
1329         CALL dbcsr_copy(almo_scf_env%matrix_ks_0deloc(ispin), &
1330                         almo_scf_env%matrix_ks(ispin))
1331         CALL dbcsr_copy(almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
1332                         almo_scf_env%matrix_sigma_inv(ispin))
1333      ENDDO
1334
1335      CALL timestop(handle)
1336
1337   END SUBROUTINE almo_scf_main
1338
1339! **************************************************************************************************
1340!> \brief selects various post scf routines
1341!> \param qs_env ...
1342!> \param almo_scf_env ...
1343!> \par History
1344!>       2011.06 created [Rustam Z Khaliullin]
1345!> \author Rustam Z Khaliullin
1346! **************************************************************************************************
1347   SUBROUTINE almo_scf_delocalization(qs_env, almo_scf_env)
1348
1349      TYPE(qs_environment_type), POINTER                 :: qs_env
1350      TYPE(almo_scf_env_type)                            :: almo_scf_env
1351
1352      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_delocalization', &
1353         routineP = moduleN//':'//routineN
1354
1355      INTEGER                                            :: col, handle, hold, iblock_col, &
1356                                                            iblock_row, ispin, mynode, &
1357                                                            nblkcols_tot, nblkrows_tot, row, &
1358                                                            unit_nr
1359      LOGICAL                                            :: almo_experimental, tr
1360      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_new_block
1361      TYPE(cp_logger_type), POINTER                      :: logger
1362      TYPE(dbcsr_distribution_type)                      :: dist
1363      TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:)        :: no_quench
1364      TYPE(optimizer_options_type)                       :: arbitrary_optimizer
1365
1366      CALL timeset(routineN, handle)
1367
1368      ! get a useful output_unit
1369      logger => cp_get_default_logger()
1370      IF (logger%para_env%ionode) THEN
1371         unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1372      ELSE
1373         unit_nr = -1
1374      ENDIF
1375
1376      ! create a local optimizer that handles XALMO DIIS
1377      ! the options of this optimizer are arbitrary because
1378      ! XALMO DIIS SCF does not converge for yet unknown reasons and
1379      ! currently used in the code to get perturbative estimates only
1380      arbitrary_optimizer%optimizer_type = optimizer_diis
1381      arbitrary_optimizer%max_iter = 3
1382      arbitrary_optimizer%eps_error = 1.0E-6_dp
1383      arbitrary_optimizer%ndiis = 2
1384
1385      SELECT CASE (almo_scf_env%deloc_method)
1386      CASE (almo_deloc_x, almo_deloc_scf, almo_deloc_x_then_scf)
1387
1388         ! RZK-warning hack into the quenched routine:
1389         ! create a quench matrix with all-all-all blocks 1.0
1390         ! it is a waste of memory but since matrices are distributed
1391         ! we can tolerate it for now
1392         ALLOCATE (no_quench(almo_scf_env%nspins))
1393         CALL dbcsr_create(no_quench(1), &
1394                           template=almo_scf_env%matrix_t(1), &
1395                           matrix_type=dbcsr_type_no_symmetry)
1396         CALL dbcsr_get_info(no_quench(1), distribution=dist)
1397         CALL dbcsr_distribution_get(dist, mynode=mynode)
1398         CALL dbcsr_work_create(no_quench(1), &
1399                                work_mutable=.TRUE.)
1400         nblkrows_tot = dbcsr_nblkrows_total(no_quench(1))
1401         nblkcols_tot = dbcsr_nblkcols_total(no_quench(1))
1402         ! RZK-warning: is it a quadratic-scaling routine?
1403         ! As a matter of fact it is! But this block treats
1404         ! fully delocalized MOs. So it is unavoidable.
1405         ! C'est la vie
1406         DO row = 1, nblkrows_tot
1407            DO col = 1, nblkcols_tot
1408               tr = .FALSE.
1409               iblock_row = row
1410               iblock_col = col
1411               CALL dbcsr_get_stored_coordinates(no_quench(1), &
1412                                                 iblock_row, iblock_col, hold)
1413               IF (hold .EQ. mynode) THEN
1414                  NULLIFY (p_new_block)
1415                  CALL dbcsr_reserve_block2d(no_quench(1), &
1416                                             iblock_row, iblock_col, p_new_block)
1417                  CPASSERT(ASSOCIATED(p_new_block))
1418                  p_new_block(:, :) = 1.0_dp
1419               ENDIF
1420            ENDDO
1421         ENDDO
1422         CALL dbcsr_finalize(no_quench(1))
1423         IF (almo_scf_env%nspins .GT. 1) THEN
1424            DO ispin = 2, almo_scf_env%nspins
1425               CALL dbcsr_create(no_quench(ispin), &
1426                                 template=almo_scf_env%matrix_t(1), &
1427                                 matrix_type=dbcsr_type_no_symmetry)
1428               CALL dbcsr_copy(no_quench(ispin), no_quench(1))
1429            ENDDO
1430         ENDIF
1431
1432      END SELECT
1433
1434      SELECT CASE (almo_scf_env%deloc_method)
1435      CASE (almo_deloc_none, almo_deloc_scf)
1436
1437         DO ispin = 1, almo_scf_env%nspins
1438            CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
1439                            almo_scf_env%matrix_t_blk(ispin))
1440         ENDDO
1441
1442      CASE (almo_deloc_x, almo_deloc_xk, almo_deloc_x_then_scf)
1443
1444         !!!! RZK-warning a whole class of delocalization methods
1445         !!!! are commented out at the moment because some of their
1446         !!!! routines have not been thoroughly tested.
1447
1448         !!!! if we have virtuals pre-optimize and truncate them
1449         !!!IF (almo_scf_env%need_virtuals) THEN
1450         !!!   SELECT CASE (almo_scf_env%deloc_truncate_virt)
1451         !!!   CASE (virt_full)
1452         !!!      ! simply copy virtual orbitals from matrix_v_full_blk to matrix_v_blk
1453         !!!      DO ispin=1,almo_scf_env%nspins
1454         !!!         CALL dbcsr_copy(almo_scf_env%matrix_v_blk(ispin),&
1455         !!!                 almo_scf_env%matrix_v_full_blk(ispin))
1456         !!!      ENDDO
1457         !!!   CASE (virt_number,virt_occ_size)
1458         !!!      CALL split_v_blk(almo_scf_env)
1459         !!!      !CALL truncate_subspace_v_blk(qs_env,almo_scf_env)
1460         !!!   CASE DEFAULT
1461         !!!      CPErrorMessage(cp_failure_level,routineP,"illegal method for virtual space truncation")
1462         !!!      CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
1463         !!!   END SELECT
1464         !!!ENDIF
1465         !!!CALL harris_foulkes_correction(qs_env,almo_scf_env)
1466
1467         CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1468                                 almo_scf_env=almo_scf_env, &
1469                                 optimizer=almo_scf_env%opt_xalmo_pcg, &
1470                                 quench_t=no_quench, &
1471                                 matrix_t_in=almo_scf_env%matrix_t_blk, &
1472                                 matrix_t_out=almo_scf_env%matrix_t, &
1473                                 !assume_t0_q0x=.TRUE., &
1474                                 assume_t0_q0x=(almo_scf_env%xalmo_trial_wf .EQ. xalmo_trial_r0_out), &
1475                                 perturbation_only=.TRUE., &
1476                                 special_case=xalmo_case_fully_deloc)
1477
1478      CASE (almo_deloc_xalmo_1diag)
1479
1480         IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_diag) THEN
1481
1482            almo_scf_env%perturbative_delocalization = .TRUE.
1483            DO ispin = 1, almo_scf_env%nspins
1484               CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
1485                               almo_scf_env%matrix_t_blk(ispin))
1486            ENDDO
1487            CALL almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, &
1488                                            arbitrary_optimizer)
1489
1490         ELSE
1491
1492            CPABORT("Other algorithms do not exist")
1493
1494         ENDIF
1495
1496      CASE (almo_deloc_xalmo_x)
1497
1498         IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_pcg) THEN
1499
1500            CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1501                                    almo_scf_env=almo_scf_env, &
1502                                    optimizer=almo_scf_env%opt_xalmo_pcg, &
1503                                    quench_t=almo_scf_env%quench_t, &
1504                                    matrix_t_in=almo_scf_env%matrix_t_blk, &
1505                                    matrix_t_out=almo_scf_env%matrix_t, &
1506                                    !assume_t0_q0x=.TRUE., &
1507                                    assume_t0_q0x=(almo_scf_env%xalmo_trial_wf .EQ. xalmo_trial_r0_out), &
1508                                    perturbation_only=.TRUE., &
1509                                    special_case=xalmo_case_normal)
1510
1511         ELSE
1512
1513            CPABORT("Other algorithms do not exist")
1514
1515         ENDIF
1516
1517      CASE (almo_deloc_xalmo_scf)
1518
1519         IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_diag) THEN
1520
1521            CPABORT("Should not be here: convergence will fail!")
1522
1523            almo_scf_env%perturbative_delocalization = .FALSE.
1524            DO ispin = 1, almo_scf_env%nspins
1525               CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
1526                               almo_scf_env%matrix_t_blk(ispin))
1527            ENDDO
1528            CALL almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, &
1529                                            arbitrary_optimizer)
1530
1531         ELSE IF (almo_scf_env%xalmo_update_algorithm .EQ. almo_scf_pcg) THEN
1532
1533            CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1534                                    almo_scf_env=almo_scf_env, &
1535                                    optimizer=almo_scf_env%opt_xalmo_pcg, &
1536                                    quench_t=almo_scf_env%quench_t, &
1537                                    matrix_t_in=almo_scf_env%matrix_t_blk, &
1538                                    matrix_t_out=almo_scf_env%matrix_t, &
1539                                    !assume_t0_q0x=.TRUE., &
1540                                    assume_t0_q0x=(almo_scf_env%xalmo_trial_wf .EQ. xalmo_trial_r0_out), &
1541                                    perturbation_only=.FALSE., &
1542                                    special_case=xalmo_case_normal)
1543
1544            ! RZK-warning THIS IS A HACK TO GET ORBITAL ENERGIES
1545            almo_experimental = .FALSE.
1546            IF (almo_experimental) THEN
1547               almo_scf_env%perturbative_delocalization = .TRUE.
1548               !DO ispin=1,almo_scf_env%nspins
1549               !   CALL dbcsr_copy(almo_scf_env%matrix_t(ispin),&
1550               !           almo_scf_env%matrix_t_blk(ispin))
1551               !ENDDO
1552               CALL almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, &
1553                                               arbitrary_optimizer)
1554            ENDIF ! experimental
1555
1556         ENDIF
1557
1558      CASE DEFAULT
1559
1560         CPABORT("Illegal delocalization method")
1561
1562      END SELECT
1563
1564      SELECT CASE (almo_scf_env%deloc_method)
1565      CASE (almo_deloc_scf, almo_deloc_x_then_scf)
1566
1567         IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
1568            CPABORT("full scf is NYI for truncated virtual space")
1569         ENDIF
1570
1571         CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1572                                 almo_scf_env=almo_scf_env, &
1573                                 optimizer=almo_scf_env%opt_xalmo_pcg, &
1574                                 quench_t=no_quench, &
1575                                 matrix_t_in=almo_scf_env%matrix_t, &
1576                                 matrix_t_out=almo_scf_env%matrix_t, &
1577                                 assume_t0_q0x=.FALSE., &
1578                                 perturbation_only=.FALSE., &
1579                                 special_case=xalmo_case_fully_deloc)
1580
1581      END SELECT
1582
1583      ! clean up
1584      SELECT CASE (almo_scf_env%deloc_method)
1585      CASE (almo_deloc_x, almo_deloc_scf, almo_deloc_x_then_scf)
1586         DO ispin = 1, almo_scf_env%nspins
1587            CALL dbcsr_release(no_quench(ispin))
1588         ENDDO
1589         DEALLOCATE (no_quench)
1590      END SELECT
1591
1592      CALL timestop(handle)
1593
1594   END SUBROUTINE almo_scf_delocalization
1595
1596! *****************************************************************************
1597!> \brief after SCF we have the final density and KS matrices compute various
1598!>        post-scf quantities
1599!> \param qs_env ...
1600!> \param almo_scf_env ...
1601!> \par History
1602!>       2015.03 created  [Rustam Z. Khaliullin]
1603!> \author Rustam Z. Khaliullin
1604! **************************************************************************************************
1605   SUBROUTINE almo_scf_post(qs_env, almo_scf_env)
1606      TYPE(qs_environment_type), POINTER                 :: qs_env
1607      TYPE(almo_scf_env_type)                            :: almo_scf_env
1608
1609      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_post', routineP = moduleN//':'//routineN
1610
1611      INTEGER                                            :: handle, ispin
1612      TYPE(cp_fm_type), POINTER                          :: mo_coeff
1613      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_w
1614      TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:)        :: matrix_t_processed
1615      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos
1616      TYPE(qs_scf_env_type), POINTER                     :: scf_env
1617
1618      CALL timeset(routineN, handle)
1619
1620      ! store matrices to speed up the next scf run
1621      CALL almo_scf_store_extrapolation_data(almo_scf_env)
1622
1623      ! orthogonalize orbitals before returning them to QS
1624      ALLOCATE (matrix_t_processed(almo_scf_env%nspins))
1625
1626      DO ispin = 1, almo_scf_env%nspins
1627
1628         CALL dbcsr_create(matrix_t_processed(ispin), &
1629                           template=almo_scf_env%matrix_t(ispin), &
1630                           matrix_type=dbcsr_type_no_symmetry)
1631
1632         CALL dbcsr_copy(matrix_t_processed(ispin), &
1633                         almo_scf_env%matrix_t(ispin))
1634
1635         IF (almo_scf_env%return_orthogonalized_mos) THEN
1636
1637            CALL orthogonalize_mos(ket=matrix_t_processed(ispin), &
1638                                   overlap=almo_scf_env%matrix_sigma(ispin), &
1639                                   metric=almo_scf_env%matrix_s(1), &
1640                                   retain_locality=.FALSE., &
1641                                   only_normalize=.FALSE., &
1642                                   nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1643                                   eps_filter=almo_scf_env%eps_filter, &
1644                                   order_lanczos=almo_scf_env%order_lanczos, &
1645                                   eps_lanczos=almo_scf_env%eps_lanczos, &
1646                                   max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
1647                                   smear=almo_scf_env%smear)
1648         ENDIF
1649
1650      ENDDO
1651
1652      !! RS-WARNING: If smearing ALMO is requested, rescaled fully-occupied orbitals are returned to QS
1653      !! RS-WARNING: Beware that QS will not be informed about electronic entropy.
1654      !!             If you want a quick and dirty transfer to QS energy, uncomment the following hack:
1655      !! IF (almo_scf_env%smear) THEN
1656      !!    qs_env%energy%kTS = 0.0_dp
1657      !!    DO ispin = 1, almo_scf_env%nspins
1658      !!       qs_env%energy%kTS = qs_env%energy%kTS + almo_scf_env%kTS(ispin)
1659      !!    END DO
1660      !! END IF
1661
1662      ! return orbitals to QS
1663      NULLIFY (mos, mo_coeff, scf_env)
1664
1665      CALL get_qs_env(qs_env, mos=mos, scf_env=scf_env)
1666
1667      DO ispin = 1, almo_scf_env%nspins
1668         ! Currently only fm version of mo_set is usable.
1669         ! First transform the matrix_t to fm version
1670         CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff)
1671         CALL copy_dbcsr_to_fm(matrix_t_processed(ispin), mo_coeff)
1672         CALL dbcsr_release(matrix_t_processed(ispin))
1673      ENDDO
1674      DEALLOCATE (matrix_t_processed)
1675
1676      ! calculate post scf properties
1677      ! CALL almo_post_scf_compute_properties(qs_env, almo_scf_env)
1678      CALL almo_post_scf_compute_properties(qs_env)
1679
1680      ! In normal scf procedure this is done by scf,
1681      ! but it's not called by ALMO
1682      ! best to release it somewhere else
1683      IF (ASSOCIATED(scf_env)) THEN
1684         CALL scf_env_release(scf_env)
1685      ENDIF
1686
1687      ! compute the W matrix if associated
1688      IF (almo_scf_env%calc_forces) THEN
1689         CALL get_qs_env(qs_env, matrix_w=matrix_w)
1690         IF (ASSOCIATED(matrix_w)) THEN
1691            CALL calculate_w_matrix_almo(matrix_w, almo_scf_env)
1692         ELSE
1693            CPABORT("Matrix W is needed but not associated")
1694         ENDIF
1695      ENDIF
1696
1697      CALL timestop(handle)
1698
1699   END SUBROUTINE almo_scf_post
1700
1701! **************************************************************************************************
1702!> \brief create various matrices
1703!> \param almo_scf_env ...
1704!> \param matrix_s0 ...
1705!> \par History
1706!>       2011.07 created [Rustam Z Khaliullin]
1707!> \author Rustam Z Khaliullin
1708! **************************************************************************************************
1709   SUBROUTINE almo_scf_env_create_matrices(almo_scf_env, matrix_s0)
1710
1711      TYPE(almo_scf_env_type), INTENT(INOUT)             :: almo_scf_env
1712      TYPE(dbcsr_type), INTENT(IN)                       :: matrix_s0
1713
1714      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_env_create_matrices', &
1715         routineP = moduleN//':'//routineN
1716
1717      INTEGER                                            :: handle, ispin, nspins
1718
1719      CALL timeset(routineN, handle)
1720
1721      nspins = almo_scf_env%nspins
1722
1723      ! AO overlap matrix and its various functions
1724      CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s(1), &
1725                              matrix_qs=matrix_s0, &
1726                              almo_scf_env=almo_scf_env, &
1727                              name_new="S", &
1728                              size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
1729                              symmetry_new=dbcsr_type_symmetric, &
1730                              spin_key=0, &
1731                              init_domains=.FALSE.)
1732      CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk(1), &
1733                              matrix_qs=matrix_s0, &
1734                              almo_scf_env=almo_scf_env, &
1735                              name_new="S_BLK", &
1736                              size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
1737                              symmetry_new=dbcsr_type_symmetric, &
1738                              spin_key=0, &
1739                              init_domains=.TRUE.)
1740      IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_diag) THEN
1741         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk_sqrt_inv(1), &
1742                                 matrix_qs=matrix_s0, &
1743                                 almo_scf_env=almo_scf_env, &
1744                                 name_new="S_BLK_SQRT_INV", &
1745                                 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
1746                                 symmetry_new=dbcsr_type_symmetric, &
1747                                 spin_key=0, &
1748                                 init_domains=.TRUE.)
1749         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk_sqrt(1), &
1750                                 matrix_qs=matrix_s0, &
1751                                 almo_scf_env=almo_scf_env, &
1752                                 name_new="S_BLK_SQRT", &
1753                                 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
1754                                 symmetry_new=dbcsr_type_symmetric, &
1755                                 spin_key=0, &
1756                                 init_domains=.TRUE.)
1757      ELSE IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_dm_sign) THEN
1758         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk_inv(1), &
1759                                 matrix_qs=matrix_s0, &
1760                                 almo_scf_env=almo_scf_env, &
1761                                 name_new="S_BLK_INV", &
1762                                 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
1763                                 symmetry_new=dbcsr_type_symmetric, &
1764                                 spin_key=0, &
1765                                 init_domains=.TRUE.)
1766      ENDIF
1767
1768      ! MO coeff matrices and their derivatives
1769      ALLOCATE (almo_scf_env%matrix_t_blk(nspins))
1770      ALLOCATE (almo_scf_env%quench_t_blk(nspins))
1771      ALLOCATE (almo_scf_env%matrix_err_blk(nspins))
1772      ALLOCATE (almo_scf_env%matrix_err_xx(nspins))
1773      ALLOCATE (almo_scf_env%matrix_sigma(nspins))
1774      ALLOCATE (almo_scf_env%matrix_sigma_inv(nspins))
1775      ALLOCATE (almo_scf_env%matrix_sigma_sqrt(nspins))
1776      ALLOCATE (almo_scf_env%matrix_sigma_sqrt_inv(nspins))
1777      ALLOCATE (almo_scf_env%matrix_sigma_blk(nspins))
1778      ALLOCATE (almo_scf_env%matrix_sigma_inv_0deloc(nspins))
1779      ALLOCATE (almo_scf_env%matrix_t(nspins))
1780      ALLOCATE (almo_scf_env%matrix_t_tr(nspins))
1781      DO ispin = 1, nspins
1782         ! create the blocked quencher
1783         CALL matrix_almo_create(matrix_new=almo_scf_env%quench_t_blk(ispin), &
1784                                 matrix_qs=matrix_s0, &
1785                                 almo_scf_env=almo_scf_env, &
1786                                 name_new="Q_BLK", &
1787                                 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_occ/), &
1788                                 symmetry_new=dbcsr_type_no_symmetry, &
1789                                 spin_key=ispin, &
1790                                 init_domains=.TRUE.)
1791         ! create ALMO coefficient matrix
1792         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_t_blk(ispin), &
1793                                 matrix_qs=matrix_s0, &
1794                                 almo_scf_env=almo_scf_env, &
1795                                 name_new="T_BLK", &
1796                                 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_occ/), &
1797                                 symmetry_new=dbcsr_type_no_symmetry, &
1798                                 spin_key=ispin, &
1799                                 init_domains=.TRUE.)
1800         ! create the error matrix
1801         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_err_blk(ispin), &
1802                                 matrix_qs=matrix_s0, &
1803                                 almo_scf_env=almo_scf_env, &
1804                                 name_new="ERR_BLK", &
1805                                 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
1806                                 symmetry_new=dbcsr_type_no_symmetry, &
1807                                 spin_key=ispin, &
1808                                 init_domains=.TRUE.)
1809         ! create the error matrix for the quenched ALMOs
1810         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_err_xx(ispin), &
1811                                 matrix_qs=matrix_s0, &
1812                                 almo_scf_env=almo_scf_env, &
1813                                 name_new="ERR_XX", &
1814                                 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_occ/), &
1815                                 symmetry_new=dbcsr_type_no_symmetry, &
1816                                 spin_key=ispin, &
1817                                 init_domains=.FALSE.)
1818         ! create a matrix with dimensions of a transposed mo coefficient matrix
1819         ! it might be necessary to perform the correction step using cayley
1820         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_t_tr(ispin), &
1821                                 matrix_qs=matrix_s0, &
1822                                 almo_scf_env=almo_scf_env, &
1823                                 name_new="T_TR", &
1824                                 size_keys=(/almo_mat_dim_occ, almo_mat_dim_aobasis/), &
1825                                 symmetry_new=dbcsr_type_no_symmetry, &
1826                                 spin_key=ispin, &
1827                                 init_domains=.FALSE.)
1828         ! create mo overlap matrix
1829         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma(ispin), &
1830                                 matrix_qs=matrix_s0, &
1831                                 almo_scf_env=almo_scf_env, &
1832                                 name_new="SIG", &
1833                                 size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
1834                                 symmetry_new=dbcsr_type_symmetric, &
1835                                 spin_key=ispin, &
1836                                 init_domains=.FALSE.)
1837         ! create blocked mo overlap matrix
1838         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_blk(ispin), &
1839                                 matrix_qs=matrix_s0, &
1840                                 almo_scf_env=almo_scf_env, &
1841                                 name_new="SIG_BLK", &
1842                                 size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
1843                                 symmetry_new=dbcsr_type_symmetric, &
1844                                 spin_key=ispin, &
1845                                 init_domains=.TRUE.)
1846         ! create blocked inverse mo overlap matrix
1847         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
1848                                 matrix_qs=matrix_s0, &
1849                                 almo_scf_env=almo_scf_env, &
1850                                 name_new="SIGINV_BLK", &
1851                                 size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
1852                                 symmetry_new=dbcsr_type_symmetric, &
1853                                 spin_key=ispin, &
1854                                 init_domains=.TRUE.)
1855         ! create inverse mo overlap matrix
1856         CALL matrix_almo_create( &
1857            matrix_new=almo_scf_env%matrix_sigma_inv(ispin), &
1858            matrix_qs=matrix_s0, &
1859            almo_scf_env=almo_scf_env, &
1860            name_new="SIGINV", &
1861            size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
1862            symmetry_new=dbcsr_type_symmetric, &
1863            spin_key=ispin, &
1864            init_domains=.FALSE.)
1865         ! create various templates that will be necessary later
1866         CALL matrix_almo_create( &
1867            matrix_new=almo_scf_env%matrix_t(ispin), &
1868            matrix_qs=matrix_s0, &
1869            almo_scf_env=almo_scf_env, &
1870            name_new="T", &
1871            size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_occ/), &
1872            symmetry_new=dbcsr_type_no_symmetry, &
1873            spin_key=ispin, &
1874            init_domains=.FALSE.)
1875         CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt(ispin), &
1876                           template=almo_scf_env%matrix_sigma(ispin), &
1877                           matrix_type=dbcsr_type_no_symmetry)
1878         CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
1879                           template=almo_scf_env%matrix_sigma(ispin), &
1880                           matrix_type=dbcsr_type_no_symmetry)
1881      ENDDO
1882
1883      ! create virtual orbitals if necessary
1884      IF (almo_scf_env%need_virtuals) THEN
1885         ALLOCATE (almo_scf_env%matrix_v_blk(nspins))
1886         ALLOCATE (almo_scf_env%matrix_v_full_blk(nspins))
1887         ALLOCATE (almo_scf_env%matrix_v(nspins))
1888         ALLOCATE (almo_scf_env%matrix_vo(nspins))
1889         ALLOCATE (almo_scf_env%matrix_x(nspins))
1890         ALLOCATE (almo_scf_env%matrix_ov(nspins))
1891         ALLOCATE (almo_scf_env%matrix_ov_full(nspins))
1892         ALLOCATE (almo_scf_env%matrix_sigma_vv(nspins))
1893         ALLOCATE (almo_scf_env%matrix_sigma_vv_blk(nspins))
1894         ALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt(nspins))
1895         ALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt_inv(nspins))
1896         ALLOCATE (almo_scf_env%matrix_vv_full_blk(nspins))
1897
1898         IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
1899            ALLOCATE (almo_scf_env%matrix_k_blk(nspins))
1900            ALLOCATE (almo_scf_env%matrix_k_blk_ones(nspins))
1901            ALLOCATE (almo_scf_env%matrix_k_tr(nspins))
1902            ALLOCATE (almo_scf_env%matrix_v_disc(nspins))
1903            ALLOCATE (almo_scf_env%matrix_v_disc_blk(nspins))
1904            ALLOCATE (almo_scf_env%matrix_ov_disc(nspins))
1905            ALLOCATE (almo_scf_env%matrix_vv_disc_blk(nspins))
1906            ALLOCATE (almo_scf_env%matrix_vv_disc(nspins))
1907            ALLOCATE (almo_scf_env%opt_k_t_dd(nspins))
1908            ALLOCATE (almo_scf_env%opt_k_t_rr(nspins))
1909            ALLOCATE (almo_scf_env%opt_k_denom(nspins))
1910         ENDIF
1911
1912         DO ispin = 1, nspins
1913            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_full_blk(ispin), &
1914                                    matrix_qs=matrix_s0, &
1915                                    almo_scf_env=almo_scf_env, &
1916                                    name_new="V_FULL_BLK", &
1917                                    size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_virt_full/), &
1918                                    symmetry_new=dbcsr_type_no_symmetry, &
1919                                    spin_key=ispin, &
1920                                    init_domains=.FALSE.)
1921            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_blk(ispin), &
1922                                    matrix_qs=matrix_s0, &
1923                                    almo_scf_env=almo_scf_env, &
1924                                    name_new="V_BLK", &
1925                                    size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_virt/), &
1926                                    symmetry_new=dbcsr_type_no_symmetry, &
1927                                    spin_key=ispin, &
1928                                    init_domains=.FALSE.)
1929            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v(ispin), &
1930                                    matrix_qs=matrix_s0, &
1931                                    almo_scf_env=almo_scf_env, &
1932                                    name_new="V", &
1933                                    size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_virt/), &
1934                                    symmetry_new=dbcsr_type_no_symmetry, &
1935                                    spin_key=ispin, &
1936                                    init_domains=.FALSE.)
1937            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ov_full(ispin), &
1938                                    matrix_qs=matrix_s0, &
1939                                    almo_scf_env=almo_scf_env, &
1940                                    name_new="OV_FULL", &
1941                                    size_keys=(/almo_mat_dim_occ, almo_mat_dim_virt_full/), &
1942                                    symmetry_new=dbcsr_type_no_symmetry, &
1943                                    spin_key=ispin, &
1944                                    init_domains=.FALSE.)
1945            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ov(ispin), &
1946                                    matrix_qs=matrix_s0, &
1947                                    almo_scf_env=almo_scf_env, &
1948                                    name_new="OV", &
1949                                    size_keys=(/almo_mat_dim_occ, almo_mat_dim_virt/), &
1950                                    symmetry_new=dbcsr_type_no_symmetry, &
1951                                    spin_key=ispin, &
1952                                    init_domains=.FALSE.)
1953            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vo(ispin), &
1954                                    matrix_qs=matrix_s0, &
1955                                    almo_scf_env=almo_scf_env, &
1956                                    name_new="VO", &
1957                                    size_keys=(/almo_mat_dim_virt, almo_mat_dim_occ/), &
1958                                    symmetry_new=dbcsr_type_no_symmetry, &
1959                                    spin_key=ispin, &
1960                                    init_domains=.FALSE.)
1961            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_x(ispin), &
1962                                    matrix_qs=matrix_s0, &
1963                                    almo_scf_env=almo_scf_env, &
1964                                    name_new="VO", &
1965                                    size_keys=(/almo_mat_dim_virt, almo_mat_dim_occ/), &
1966                                    symmetry_new=dbcsr_type_no_symmetry, &
1967                                    spin_key=ispin, &
1968                                    init_domains=.FALSE.)
1969            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_vv(ispin), &
1970                                    matrix_qs=matrix_s0, &
1971                                    almo_scf_env=almo_scf_env, &
1972                                    name_new="SIG_VV", &
1973                                    size_keys=(/almo_mat_dim_virt, almo_mat_dim_virt/), &
1974                                    symmetry_new=dbcsr_type_symmetric, &
1975                                    spin_key=ispin, &
1976                                    init_domains=.FALSE.)
1977            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vv_full_blk(ispin), &
1978                                    matrix_qs=matrix_s0, &
1979                                    almo_scf_env=almo_scf_env, &
1980                                    name_new="VV_FULL_BLK", &
1981                                    size_keys=(/almo_mat_dim_virt_full, almo_mat_dim_virt_full/), &
1982                                    symmetry_new=dbcsr_type_no_symmetry, &
1983                                    spin_key=ispin, &
1984                                    init_domains=.TRUE.)
1985            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_vv_blk(ispin), &
1986                                    matrix_qs=matrix_s0, &
1987                                    almo_scf_env=almo_scf_env, &
1988                                    name_new="SIG_VV_BLK", &
1989                                    size_keys=(/almo_mat_dim_virt, almo_mat_dim_virt/), &
1990                                    symmetry_new=dbcsr_type_symmetric, &
1991                                    spin_key=ispin, &
1992                                    init_domains=.TRUE.)
1993            CALL dbcsr_create(almo_scf_env%matrix_sigma_vv_sqrt(ispin), &
1994                              template=almo_scf_env%matrix_sigma_vv(ispin), &
1995                              matrix_type=dbcsr_type_no_symmetry)
1996            CALL dbcsr_create(almo_scf_env%matrix_sigma_vv_sqrt_inv(ispin), &
1997                              template=almo_scf_env%matrix_sigma_vv(ispin), &
1998                              matrix_type=dbcsr_type_no_symmetry)
1999
2000            IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
2001               CALL matrix_almo_create(matrix_new=almo_scf_env%opt_k_t_rr(ispin), &
2002                                       matrix_qs=matrix_s0, &
2003                                       almo_scf_env=almo_scf_env, &
2004                                       name_new="OPT_K_U_RR", &
2005                                       size_keys=(/almo_mat_dim_virt, almo_mat_dim_virt/), &
2006                                       symmetry_new=dbcsr_type_no_symmetry, &
2007                                       spin_key=ispin, &
2008                                       init_domains=.FALSE.)
2009               CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vv_disc(ispin), &
2010                                       matrix_qs=matrix_s0, &
2011                                       almo_scf_env=almo_scf_env, &
2012                                       name_new="VV_DISC", &
2013                                       size_keys=(/almo_mat_dim_virt_disc, almo_mat_dim_virt_disc/), &
2014                                       symmetry_new=dbcsr_type_symmetric, &
2015                                       spin_key=ispin, &
2016                                       init_domains=.FALSE.)
2017               CALL matrix_almo_create(matrix_new=almo_scf_env%opt_k_t_dd(ispin), &
2018                                       matrix_qs=matrix_s0, &
2019                                       almo_scf_env=almo_scf_env, &
2020                                       name_new="OPT_K_U_DD", &
2021                                       size_keys=(/almo_mat_dim_virt_disc, almo_mat_dim_virt_disc/), &
2022                                       symmetry_new=dbcsr_type_no_symmetry, &
2023                                       spin_key=ispin, &
2024                                       init_domains=.FALSE.)
2025               CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vv_disc_blk(ispin), &
2026                                       matrix_qs=matrix_s0, &
2027                                       almo_scf_env=almo_scf_env, &
2028                                       name_new="VV_DISC_BLK", &
2029                                       size_keys=(/almo_mat_dim_virt_disc, almo_mat_dim_virt_disc/), &
2030                                       symmetry_new=dbcsr_type_symmetric, &
2031                                       spin_key=ispin, &
2032                                       init_domains=.TRUE.)
2033               CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_k_blk(ispin), &
2034                                       matrix_qs=matrix_s0, &
2035                                       almo_scf_env=almo_scf_env, &
2036                                       name_new="K_BLK", &
2037                                       size_keys=(/almo_mat_dim_virt_disc, almo_mat_dim_virt/), &
2038                                       symmetry_new=dbcsr_type_no_symmetry, &
2039                                       spin_key=ispin, &
2040                                       init_domains=.TRUE.)
2041               CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_k_blk_ones(ispin), &
2042                                       matrix_qs=matrix_s0, &
2043                                       almo_scf_env=almo_scf_env, &
2044                                       name_new="K_BLK_1", &
2045                                       size_keys=(/almo_mat_dim_virt_disc, almo_mat_dim_virt/), &
2046                                       symmetry_new=dbcsr_type_no_symmetry, &
2047                                       spin_key=ispin, &
2048                                       init_domains=.TRUE.)
2049               CALL matrix_almo_create(matrix_new=almo_scf_env%opt_k_denom(ispin), &
2050                                       matrix_qs=matrix_s0, &
2051                                       almo_scf_env=almo_scf_env, &
2052                                       name_new="OPT_K_DENOM", &
2053                                       size_keys=(/almo_mat_dim_virt_disc, almo_mat_dim_virt/), &
2054                                       symmetry_new=dbcsr_type_no_symmetry, &
2055                                       spin_key=ispin, &
2056                                       init_domains=.FALSE.)
2057               CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_k_tr(ispin), &
2058                                       matrix_qs=matrix_s0, &
2059                                       almo_scf_env=almo_scf_env, &
2060                                       name_new="K_TR", &
2061                                       size_keys=(/almo_mat_dim_virt, almo_mat_dim_virt_disc/), &
2062                                       symmetry_new=dbcsr_type_no_symmetry, &
2063                                       spin_key=ispin, &
2064                                       init_domains=.FALSE.)
2065               CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_disc_blk(ispin), &
2066                                       matrix_qs=matrix_s0, &
2067                                       almo_scf_env=almo_scf_env, &
2068                                       name_new="V_DISC_BLK", &
2069                                       size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_virt_disc/), &
2070                                       symmetry_new=dbcsr_type_no_symmetry, &
2071                                       spin_key=ispin, &
2072                                       init_domains=.FALSE.)
2073               CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_disc(ispin), &
2074                                       matrix_qs=matrix_s0, &
2075                                       almo_scf_env=almo_scf_env, &
2076                                       name_new="V_DISC", &
2077                                       size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_virt_disc/), &
2078                                       symmetry_new=dbcsr_type_no_symmetry, &
2079                                       spin_key=ispin, &
2080                                       init_domains=.FALSE.)
2081               CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ov_disc(ispin), &
2082                                       matrix_qs=matrix_s0, &
2083                                       almo_scf_env=almo_scf_env, &
2084                                       name_new="OV_DISC", &
2085                                       size_keys=(/almo_mat_dim_occ, almo_mat_dim_virt_disc/), &
2086                                       symmetry_new=dbcsr_type_no_symmetry, &
2087                                       spin_key=ispin, &
2088                                       init_domains=.FALSE.)
2089
2090            ENDIF ! end need_discarded_virtuals
2091
2092         ENDDO ! spin
2093      ENDIF
2094
2095      ! create matrices of orbital energies if necessary
2096      IF (almo_scf_env%need_orbital_energies) THEN
2097         ALLOCATE (almo_scf_env%matrix_eoo(nspins))
2098         ALLOCATE (almo_scf_env%matrix_evv_full(nspins))
2099         DO ispin = 1, nspins
2100            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_eoo(ispin), &
2101                                    matrix_qs=matrix_s0, &
2102                                    almo_scf_env=almo_scf_env, &
2103                                    name_new="E_OCC", &
2104                                    size_keys=(/almo_mat_dim_occ, almo_mat_dim_occ/), &
2105                                    symmetry_new=dbcsr_type_no_symmetry, &
2106                                    spin_key=ispin, &
2107                                    init_domains=.FALSE.)
2108            CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_evv_full(ispin), &
2109                                    matrix_qs=matrix_s0, &
2110                                    almo_scf_env=almo_scf_env, &
2111                                    name_new="E_VIRT", &
2112                                    size_keys=(/almo_mat_dim_virt_full, almo_mat_dim_virt_full/), &
2113                                    symmetry_new=dbcsr_type_no_symmetry, &
2114                                    spin_key=ispin, &
2115                                    init_domains=.FALSE.)
2116         ENDDO
2117      ENDIF
2118
2119      ! Density and KS matrices
2120      ALLOCATE (almo_scf_env%matrix_p(nspins))
2121      ALLOCATE (almo_scf_env%matrix_p_blk(nspins))
2122      ALLOCATE (almo_scf_env%matrix_ks(nspins))
2123      ALLOCATE (almo_scf_env%matrix_ks_blk(nspins))
2124      IF (almo_scf_env%need_previous_ks) &
2125         ALLOCATE (almo_scf_env%matrix_ks_0deloc(nspins))
2126      DO ispin = 1, nspins
2127         ! RZK-warning copy with symmery but remember that this might cause problems
2128         CALL dbcsr_create(almo_scf_env%matrix_p(ispin), &
2129                           template=almo_scf_env%matrix_s(1), &
2130                           matrix_type=dbcsr_type_symmetric)
2131         CALL dbcsr_create(almo_scf_env%matrix_ks(ispin), &
2132                           template=almo_scf_env%matrix_s(1), &
2133                           matrix_type=dbcsr_type_symmetric)
2134         IF (almo_scf_env%need_previous_ks) THEN
2135            CALL dbcsr_create(almo_scf_env%matrix_ks_0deloc(ispin), &
2136                              template=almo_scf_env%matrix_s(1), &
2137                              matrix_type=dbcsr_type_symmetric)
2138         ENDIF
2139         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_p_blk(ispin), &
2140                                 matrix_qs=matrix_s0, &
2141                                 almo_scf_env=almo_scf_env, &
2142                                 name_new="P_BLK", &
2143                                 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2144                                 symmetry_new=dbcsr_type_symmetric, &
2145                                 spin_key=ispin, &
2146                                 init_domains=.TRUE.)
2147         CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ks_blk(ispin), &
2148                                 matrix_qs=matrix_s0, &
2149                                 almo_scf_env=almo_scf_env, &
2150                                 name_new="KS_BLK", &
2151                                 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_aobasis/), &
2152                                 symmetry_new=dbcsr_type_symmetric, &
2153                                 spin_key=ispin, &
2154                                 init_domains=.TRUE.)
2155      ENDDO
2156
2157      CALL timestop(handle)
2158
2159   END SUBROUTINE almo_scf_env_create_matrices
2160
2161! **************************************************************************************************
2162!> \brief clean up procedures for almo scf
2163!> \param almo_scf_env ...
2164!> \par History
2165!>       2011.06 created [Rustam Z Khaliullin]
2166!>       2018.09 smearing support [Ruben Staub]
2167!> \author Rustam Z Khaliullin
2168! **************************************************************************************************
2169   SUBROUTINE almo_scf_clean_up(almo_scf_env)
2170
2171      TYPE(almo_scf_env_type)                            :: almo_scf_env
2172
2173      CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_clean_up', &
2174         routineP = moduleN//':'//routineN
2175
2176      INTEGER                                            :: handle, ispin, unit_nr
2177      TYPE(cp_logger_type), POINTER                      :: logger
2178
2179      CALL timeset(routineN, handle)
2180
2181      ! get a useful output_unit
2182      logger => cp_get_default_logger()
2183      IF (logger%para_env%ionode) THEN
2184         unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
2185      ELSE
2186         unit_nr = -1
2187      ENDIF
2188
2189      ! release matrices
2190      CALL dbcsr_release(almo_scf_env%matrix_s(1))
2191      CALL dbcsr_release(almo_scf_env%matrix_s_blk(1))
2192      IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_diag) THEN
2193         CALL dbcsr_release(almo_scf_env%matrix_s_blk_sqrt_inv(1))
2194         CALL dbcsr_release(almo_scf_env%matrix_s_blk_sqrt(1))
2195      ELSE IF (almo_scf_env%almo_update_algorithm .EQ. almo_scf_dm_sign) THEN
2196         CALL dbcsr_release(almo_scf_env%matrix_s_blk_inv(1))
2197      ENDIF
2198      DO ispin = 1, almo_scf_env%nspins
2199         CALL dbcsr_release(almo_scf_env%quench_t(ispin))
2200         CALL dbcsr_release(almo_scf_env%quench_t_blk(ispin))
2201         CALL dbcsr_release(almo_scf_env%matrix_t_blk(ispin))
2202         CALL dbcsr_release(almo_scf_env%matrix_err_blk(ispin))
2203         CALL dbcsr_release(almo_scf_env%matrix_err_xx(ispin))
2204         CALL dbcsr_release(almo_scf_env%matrix_t_tr(ispin))
2205         CALL dbcsr_release(almo_scf_env%matrix_sigma(ispin))
2206         CALL dbcsr_release(almo_scf_env%matrix_sigma_blk(ispin))
2207         CALL dbcsr_release(almo_scf_env%matrix_sigma_inv_0deloc(ispin))
2208         CALL dbcsr_release(almo_scf_env%matrix_sigma_inv(ispin))
2209         CALL dbcsr_release(almo_scf_env%matrix_t(ispin))
2210         CALL dbcsr_release(almo_scf_env%matrix_sigma_sqrt(ispin))
2211         CALL dbcsr_release(almo_scf_env%matrix_sigma_sqrt_inv(ispin))
2212         CALL dbcsr_release(almo_scf_env%matrix_p(ispin))
2213         CALL dbcsr_release(almo_scf_env%matrix_ks(ispin))
2214         CALL dbcsr_release(almo_scf_env%matrix_p_blk(ispin))
2215         CALL dbcsr_release(almo_scf_env%matrix_ks_blk(ispin))
2216         IF (almo_scf_env%need_previous_ks) THEN
2217            CALL dbcsr_release(almo_scf_env%matrix_ks_0deloc(ispin))
2218         ENDIF
2219         IF (almo_scf_env%need_virtuals) THEN
2220            CALL dbcsr_release(almo_scf_env%matrix_v_blk(ispin))
2221            CALL dbcsr_release(almo_scf_env%matrix_v_full_blk(ispin))
2222            CALL dbcsr_release(almo_scf_env%matrix_v(ispin))
2223            CALL dbcsr_release(almo_scf_env%matrix_vo(ispin))
2224            CALL dbcsr_release(almo_scf_env%matrix_x(ispin))
2225            CALL dbcsr_release(almo_scf_env%matrix_ov(ispin))
2226            CALL dbcsr_release(almo_scf_env%matrix_ov_full(ispin))
2227            CALL dbcsr_release(almo_scf_env%matrix_sigma_vv(ispin))
2228            CALL dbcsr_release(almo_scf_env%matrix_sigma_vv_blk(ispin))
2229            CALL dbcsr_release(almo_scf_env%matrix_sigma_vv_sqrt(ispin))
2230            CALL dbcsr_release(almo_scf_env%matrix_sigma_vv_sqrt_inv(ispin))
2231            CALL dbcsr_release(almo_scf_env%matrix_vv_full_blk(ispin))
2232            IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
2233               CALL dbcsr_release(almo_scf_env%matrix_k_tr(ispin))
2234               CALL dbcsr_release(almo_scf_env%matrix_k_blk(ispin))
2235               CALL dbcsr_release(almo_scf_env%matrix_k_blk_ones(ispin))
2236               CALL dbcsr_release(almo_scf_env%matrix_v_disc(ispin))
2237               CALL dbcsr_release(almo_scf_env%matrix_v_disc_blk(ispin))
2238               CALL dbcsr_release(almo_scf_env%matrix_ov_disc(ispin))
2239               CALL dbcsr_release(almo_scf_env%matrix_vv_disc_blk(ispin))
2240               CALL dbcsr_release(almo_scf_env%matrix_vv_disc(ispin))
2241               CALL dbcsr_release(almo_scf_env%opt_k_t_dd(ispin))
2242               CALL dbcsr_release(almo_scf_env%opt_k_t_rr(ispin))
2243               CALL dbcsr_release(almo_scf_env%opt_k_denom(ispin))
2244            ENDIF
2245         ENDIF
2246         IF (almo_scf_env%need_orbital_energies) THEN
2247            CALL dbcsr_release(almo_scf_env%matrix_eoo(ispin))
2248            CALL dbcsr_release(almo_scf_env%matrix_evv_full(ispin))
2249         ENDIF
2250      ENDDO
2251
2252      ! deallocate matrices
2253      DEALLOCATE (almo_scf_env%matrix_p)
2254      DEALLOCATE (almo_scf_env%matrix_p_blk)
2255      DEALLOCATE (almo_scf_env%matrix_ks)
2256      DEALLOCATE (almo_scf_env%matrix_ks_blk)
2257      DEALLOCATE (almo_scf_env%matrix_t_blk)
2258      DEALLOCATE (almo_scf_env%matrix_err_blk)
2259      DEALLOCATE (almo_scf_env%matrix_err_xx)
2260      DEALLOCATE (almo_scf_env%matrix_t)
2261      DEALLOCATE (almo_scf_env%matrix_t_tr)
2262      DEALLOCATE (almo_scf_env%matrix_sigma)
2263      DEALLOCATE (almo_scf_env%matrix_sigma_blk)
2264      DEALLOCATE (almo_scf_env%matrix_sigma_inv_0deloc)
2265      DEALLOCATE (almo_scf_env%matrix_sigma_sqrt)
2266      DEALLOCATE (almo_scf_env%matrix_sigma_sqrt_inv)
2267      DEALLOCATE (almo_scf_env%matrix_sigma_inv)
2268      DEALLOCATE (almo_scf_env%quench_t)
2269      DEALLOCATE (almo_scf_env%quench_t_blk)
2270      IF (almo_scf_env%need_virtuals) THEN
2271         DEALLOCATE (almo_scf_env%matrix_v_blk)
2272         DEALLOCATE (almo_scf_env%matrix_v_full_blk)
2273         DEALLOCATE (almo_scf_env%matrix_v)
2274         DEALLOCATE (almo_scf_env%matrix_vo)
2275         DEALLOCATE (almo_scf_env%matrix_x)
2276         DEALLOCATE (almo_scf_env%matrix_ov)
2277         DEALLOCATE (almo_scf_env%matrix_ov_full)
2278         DEALLOCATE (almo_scf_env%matrix_sigma_vv)
2279         DEALLOCATE (almo_scf_env%matrix_sigma_vv_blk)
2280         DEALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt)
2281         DEALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt_inv)
2282         DEALLOCATE (almo_scf_env%matrix_vv_full_blk)
2283         IF (almo_scf_env%deloc_truncate_virt .NE. virt_full) THEN
2284            DEALLOCATE (almo_scf_env%matrix_k_tr)
2285            DEALLOCATE (almo_scf_env%matrix_k_blk)
2286            DEALLOCATE (almo_scf_env%matrix_v_disc)
2287            DEALLOCATE (almo_scf_env%matrix_v_disc_blk)
2288            DEALLOCATE (almo_scf_env%matrix_ov_disc)
2289            DEALLOCATE (almo_scf_env%matrix_vv_disc_blk)
2290            DEALLOCATE (almo_scf_env%matrix_vv_disc)
2291            DEALLOCATE (almo_scf_env%matrix_k_blk_ones)
2292            DEALLOCATE (almo_scf_env%opt_k_t_dd)
2293            DEALLOCATE (almo_scf_env%opt_k_t_rr)
2294            DEALLOCATE (almo_scf_env%opt_k_denom)
2295         ENDIF
2296      ENDIF
2297      IF (almo_scf_env%need_previous_ks) THEN
2298         DEALLOCATE (almo_scf_env%matrix_ks_0deloc)
2299      ENDIF
2300      IF (almo_scf_env%need_orbital_energies) THEN
2301         DEALLOCATE (almo_scf_env%matrix_eoo)
2302         DEALLOCATE (almo_scf_env%matrix_evv_full)
2303      ENDIF
2304
2305      ! clean up other variables
2306      DO ispin = 1, almo_scf_env%nspins
2307         CALL release_submatrices( &
2308            almo_scf_env%domain_preconditioner(:, ispin))
2309         CALL release_submatrices(almo_scf_env%domain_s_inv(:, ispin))
2310         CALL release_submatrices(almo_scf_env%domain_s_sqrt_inv(:, ispin))
2311         CALL release_submatrices(almo_scf_env%domain_s_sqrt(:, ispin))
2312         CALL release_submatrices(almo_scf_env%domain_ks_xx(:, ispin))
2313         CALL release_submatrices(almo_scf_env%domain_t(:, ispin))
2314         CALL release_submatrices(almo_scf_env%domain_err(:, ispin))
2315         CALL release_submatrices(almo_scf_env%domain_r_down_up(:, ispin))
2316      ENDDO
2317      DEALLOCATE (almo_scf_env%domain_preconditioner)
2318      DEALLOCATE (almo_scf_env%domain_s_inv)
2319      DEALLOCATE (almo_scf_env%domain_s_sqrt_inv)
2320      DEALLOCATE (almo_scf_env%domain_s_sqrt)
2321      DEALLOCATE (almo_scf_env%domain_ks_xx)
2322      DEALLOCATE (almo_scf_env%domain_t)
2323      DEALLOCATE (almo_scf_env%domain_err)
2324      DEALLOCATE (almo_scf_env%domain_r_down_up)
2325      DO ispin = 1, almo_scf_env%nspins
2326         DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs)
2327         DEALLOCATE (almo_scf_env%domain_map(ispin)%index1)
2328      ENDDO
2329      DEALLOCATE (almo_scf_env%domain_map)
2330      DEALLOCATE (almo_scf_env%domain_index_of_ao)
2331      DEALLOCATE (almo_scf_env%domain_index_of_atom)
2332      DEALLOCATE (almo_scf_env%first_atom_of_domain)
2333      DEALLOCATE (almo_scf_env%last_atom_of_domain)
2334      DEALLOCATE (almo_scf_env%nbasis_of_domain)
2335      DEALLOCATE (almo_scf_env%nocc_of_domain)
2336      DEALLOCATE (almo_scf_env%real_ne_of_domain)
2337      DEALLOCATE (almo_scf_env%nvirt_full_of_domain)
2338      DEALLOCATE (almo_scf_env%nvirt_of_domain)
2339      DEALLOCATE (almo_scf_env%nvirt_disc_of_domain)
2340      DEALLOCATE (almo_scf_env%mu_of_domain)
2341      DEALLOCATE (almo_scf_env%cpu_of_domain)
2342      DEALLOCATE (almo_scf_env%charge_of_domain)
2343      DEALLOCATE (almo_scf_env%multiplicity_of_domain)
2344      IF (almo_scf_env%smear) THEN
2345         DEALLOCATE (almo_scf_env%mo_energies)
2346         DEALLOCATE (almo_scf_env%kTS)
2347      END IF
2348
2349      DEALLOCATE (almo_scf_env%domain_index_of_ao_block)
2350      DEALLOCATE (almo_scf_env%domain_index_of_mo_block)
2351
2352      CALL cp_para_env_release(almo_scf_env%para_env)
2353      CALL cp_blacs_env_release(almo_scf_env%blacs_env)
2354
2355      CALL timestop(handle)
2356
2357   END SUBROUTINE almo_scf_clean_up
2358
2359! **************************************************************************************************
2360!> \brief Do post scf calculations with ALMO
2361!>        WARNING: ALMO post scf calculation may not work for certain quantities,
2362!>        like forces, since ALMO wave function is only 'partially' optimized
2363!> \param qs_env ...
2364!> \par History
2365!>       2016.12 created [Yifei Shi]
2366!> \author Yifei Shi
2367! **************************************************************************************************
2368   SUBROUTINE almo_post_scf_compute_properties(qs_env)
2369      TYPE(qs_environment_type), POINTER                 :: qs_env
2370
2371      CHARACTER(len=*), PARAMETER :: routineN = 'almo_post_scf_compute_properties', &
2372         routineP = moduleN//':'//routineN
2373
2374      CALL qs_scf_compute_properties(qs_env)
2375
2376   END SUBROUTINE almo_post_scf_compute_properties
2377
2378END MODULE almo_scf
2379
2380