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