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