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 (unit_nr, "(/,T2,A)") "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 (unit_nr, "(/,T2,A)") "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 (unit_nr, "(/,T2,A)") "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 (unit_nr, "(/,T2,A)") "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 (unit_nr, "(/,T2,A)") "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 (unit_nr, "(/,T2,A)") "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 : TYPE(cp_logger_type), POINTER :: logger
1404 122 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: no_quench
1405 : TYPE(optimizer_options_type) :: arbitrary_optimizer
1406 :
1407 122 : CALL timeset(routineN, handle)
1408 :
1409 : ! get a useful output_unit
1410 122 : logger => cp_get_default_logger()
1411 122 : IF (logger%para_env%is_source()) THEN
1412 61 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1413 : ELSE
1414 : unit_nr = -1
1415 : END IF
1416 :
1417 : ! create a local optimizer that handles XALMO DIIS
1418 : ! the options of this optimizer are arbitrary because
1419 : ! XALMO DIIS SCF does not converge for yet unknown reasons and
1420 : ! currently used in the code to get perturbative estimates only
1421 122 : arbitrary_optimizer%optimizer_type = optimizer_diis
1422 122 : arbitrary_optimizer%max_iter = 3
1423 122 : arbitrary_optimizer%eps_error = 1.0E-6_dp
1424 122 : arbitrary_optimizer%ndiis = 2
1425 :
1426 152 : SELECT CASE (almo_scf_env%deloc_method)
1427 : CASE (almo_deloc_x, almo_deloc_scf, almo_deloc_x_then_scf)
1428 :
1429 : ! RZK-warning hack into the quenched routine:
1430 : ! create a quench matrix with all-all-all blocks 1.0
1431 : ! it is a waste of memory but since matrices are distributed
1432 : ! we can tolerate it for now
1433 120 : ALLOCATE (no_quench(almo_scf_env%nspins))
1434 : CALL dbcsr_create(no_quench(1), &
1435 : template=almo_scf_env%matrix_t(1), &
1436 30 : matrix_type=dbcsr_type_no_symmetry)
1437 30 : CALL dbcsr_reserve_all_blocks(no_quench(1))
1438 30 : CALL dbcsr_set(no_quench(1), 1.0_dp)
1439 152 : IF (almo_scf_env%nspins > 1) THEN
1440 0 : DO ispin = 2, almo_scf_env%nspins
1441 : CALL dbcsr_create(no_quench(ispin), &
1442 : template=almo_scf_env%matrix_t(1), &
1443 0 : matrix_type=dbcsr_type_no_symmetry)
1444 0 : CALL dbcsr_copy(no_quench(ispin), no_quench(1))
1445 : END DO
1446 : END IF
1447 :
1448 : END SELECT
1449 :
1450 170 : SELECT CASE (almo_scf_env%deloc_method)
1451 : CASE (almo_deloc_none, almo_deloc_scf)
1452 :
1453 102 : DO ispin = 1, almo_scf_env%nspins
1454 : CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
1455 102 : almo_scf_env%matrix_t_blk(ispin))
1456 : END DO
1457 :
1458 : CASE (almo_deloc_x, almo_deloc_xk, almo_deloc_x_then_scf)
1459 :
1460 : !!!! RZK-warning a whole class of delocalization methods
1461 : !!!! are commented out at the moment because some of their
1462 : !!!! routines have not been thoroughly tested.
1463 :
1464 18 : IF (almo_scf_env%xalmo_update_algorithm == almo_scf_pcg) THEN
1465 :
1466 : CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1467 : almo_scf_env=almo_scf_env, &
1468 : optimizer=almo_scf_env%opt_xalmo_pcg, &
1469 : quench_t=no_quench, &
1470 : matrix_t_in=almo_scf_env%matrix_t_blk, &
1471 : matrix_t_out=almo_scf_env%matrix_t, &
1472 : assume_t0_q0x=(almo_scf_env%xalmo_trial_wf == xalmo_trial_r0_out), &
1473 : perturbation_only=.TRUE., &
1474 18 : special_case=xalmo_case_fully_deloc)
1475 :
1476 0 : ELSE IF (almo_scf_env%xalmo_update_algorithm == almo_scf_trustr) THEN
1477 :
1478 : CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
1479 : almo_scf_env=almo_scf_env, &
1480 : optimizer=almo_scf_env%opt_xalmo_trustr, &
1481 : quench_t=no_quench, &
1482 : matrix_t_in=almo_scf_env%matrix_t_blk, &
1483 : matrix_t_out=almo_scf_env%matrix_t, &
1484 : perturbation_only=.TRUE., &
1485 0 : special_case=xalmo_case_fully_deloc)
1486 :
1487 : ELSE
1488 :
1489 0 : CPABORT("Other algorithms do not exist")
1490 :
1491 : END IF
1492 :
1493 : CASE (almo_deloc_xalmo_1diag)
1494 :
1495 2 : IF (almo_scf_env%xalmo_update_algorithm == almo_scf_diag) THEN
1496 :
1497 2 : almo_scf_env%perturbative_delocalization = .TRUE.
1498 4 : DO ispin = 1, almo_scf_env%nspins
1499 : CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
1500 4 : almo_scf_env%matrix_t_blk(ispin))
1501 : END DO
1502 : CALL almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, &
1503 2 : arbitrary_optimizer)
1504 :
1505 : ELSE
1506 :
1507 0 : CPABORT("Other algorithms do not exist")
1508 :
1509 : END IF
1510 :
1511 : CASE (almo_deloc_xalmo_x)
1512 :
1513 6 : IF (almo_scf_env%xalmo_update_algorithm == almo_scf_pcg) THEN
1514 :
1515 : CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1516 : almo_scf_env=almo_scf_env, &
1517 : optimizer=almo_scf_env%opt_xalmo_pcg, &
1518 : quench_t=almo_scf_env%quench_t, &
1519 : matrix_t_in=almo_scf_env%matrix_t_blk, &
1520 : matrix_t_out=almo_scf_env%matrix_t, &
1521 : assume_t0_q0x=(almo_scf_env%xalmo_trial_wf == xalmo_trial_r0_out), &
1522 : perturbation_only=.TRUE., &
1523 6 : special_case=xalmo_case_normal)
1524 :
1525 0 : ELSE IF (almo_scf_env%xalmo_update_algorithm == almo_scf_trustr) THEN
1526 :
1527 : CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
1528 : almo_scf_env=almo_scf_env, &
1529 : optimizer=almo_scf_env%opt_xalmo_trustr, &
1530 : quench_t=almo_scf_env%quench_t, &
1531 : matrix_t_in=almo_scf_env%matrix_t_blk, &
1532 : matrix_t_out=almo_scf_env%matrix_t, &
1533 : perturbation_only=.TRUE., &
1534 0 : special_case=xalmo_case_normal)
1535 :
1536 : ELSE
1537 :
1538 0 : CPABORT("Other algorithms do not exist")
1539 :
1540 : END IF
1541 :
1542 : CASE (almo_deloc_xalmo_scf)
1543 :
1544 48 : IF (almo_scf_env%xalmo_update_algorithm == almo_scf_diag) THEN
1545 :
1546 0 : CPABORT("Should not be here: convergence will fail!")
1547 :
1548 0 : almo_scf_env%perturbative_delocalization = .FALSE.
1549 0 : DO ispin = 1, almo_scf_env%nspins
1550 : CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
1551 0 : almo_scf_env%matrix_t_blk(ispin))
1552 : END DO
1553 : CALL almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, &
1554 0 : arbitrary_optimizer)
1555 :
1556 48 : ELSE IF (almo_scf_env%xalmo_update_algorithm == almo_scf_pcg) THEN
1557 :
1558 : CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1559 : almo_scf_env=almo_scf_env, &
1560 : optimizer=almo_scf_env%opt_xalmo_pcg, &
1561 : quench_t=almo_scf_env%quench_t, &
1562 : matrix_t_in=almo_scf_env%matrix_t_blk, &
1563 : matrix_t_out=almo_scf_env%matrix_t, &
1564 : assume_t0_q0x=(almo_scf_env%xalmo_trial_wf == xalmo_trial_r0_out), &
1565 : perturbation_only=.FALSE., &
1566 32 : special_case=xalmo_case_normal)
1567 :
1568 16 : ELSE IF (almo_scf_env%xalmo_update_algorithm == almo_scf_trustr) THEN
1569 :
1570 : CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
1571 : almo_scf_env=almo_scf_env, &
1572 : optimizer=almo_scf_env%opt_xalmo_trustr, &
1573 : quench_t=almo_scf_env%quench_t, &
1574 : matrix_t_in=almo_scf_env%matrix_t_blk, &
1575 : matrix_t_out=almo_scf_env%matrix_t, &
1576 : perturbation_only=.FALSE., &
1577 16 : special_case=xalmo_case_normal)
1578 :
1579 : ELSE
1580 :
1581 0 : CPABORT("Other algorithms do not exist")
1582 :
1583 : END IF
1584 :
1585 : CASE DEFAULT
1586 :
1587 122 : CPABORT("Illegal delocalization method")
1588 :
1589 : END SELECT
1590 :
1591 148 : SELECT CASE (almo_scf_env%deloc_method)
1592 : CASE (almo_deloc_scf, almo_deloc_x_then_scf)
1593 :
1594 26 : IF (almo_scf_env%deloc_truncate_virt /= virt_full) THEN
1595 0 : CPABORT("full scf is NYI for truncated virtual space")
1596 : END IF
1597 :
1598 148 : IF (almo_scf_env%xalmo_update_algorithm == almo_scf_pcg) THEN
1599 :
1600 : CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
1601 : almo_scf_env=almo_scf_env, &
1602 : optimizer=almo_scf_env%opt_xalmo_pcg, &
1603 : quench_t=no_quench, &
1604 : matrix_t_in=almo_scf_env%matrix_t, &
1605 : matrix_t_out=almo_scf_env%matrix_t, &
1606 : assume_t0_q0x=.FALSE., &
1607 : perturbation_only=.FALSE., &
1608 26 : special_case=xalmo_case_fully_deloc)
1609 :
1610 0 : ELSE IF (almo_scf_env%xalmo_update_algorithm == almo_scf_trustr) THEN
1611 :
1612 : CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
1613 : almo_scf_env=almo_scf_env, &
1614 : optimizer=almo_scf_env%opt_xalmo_trustr, &
1615 : quench_t=no_quench, &
1616 : matrix_t_in=almo_scf_env%matrix_t, &
1617 : matrix_t_out=almo_scf_env%matrix_t, &
1618 : perturbation_only=.FALSE., &
1619 0 : special_case=xalmo_case_fully_deloc)
1620 :
1621 : ELSE
1622 :
1623 0 : CPABORT("Other algorithms do not exist")
1624 :
1625 : END IF
1626 :
1627 : END SELECT
1628 :
1629 : ! clean up
1630 152 : SELECT CASE (almo_scf_env%deloc_method)
1631 : CASE (almo_deloc_x, almo_deloc_scf, almo_deloc_x_then_scf)
1632 60 : DO ispin = 1, almo_scf_env%nspins
1633 60 : CALL dbcsr_release(no_quench(ispin))
1634 : END DO
1635 152 : DEALLOCATE (no_quench)
1636 : END SELECT
1637 :
1638 122 : CALL timestop(handle)
1639 :
1640 244 : END SUBROUTINE almo_scf_delocalization
1641 :
1642 : ! **************************************************************************************************
1643 : !> \brief orbital localization
1644 : !> \param qs_env ...
1645 : !> \param almo_scf_env ...
1646 : !> \par History
1647 : !> 2018.09 created [Ziling Luo]
1648 : !> \author Ziling Luo
1649 : ! **************************************************************************************************
1650 122 : SUBROUTINE construct_nlmos(qs_env, almo_scf_env)
1651 :
1652 : TYPE(qs_environment_type), POINTER :: qs_env
1653 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1654 :
1655 : INTEGER :: ispin
1656 :
1657 122 : IF (almo_scf_env%construct_nlmos) THEN
1658 :
1659 8 : DO ispin = 1, almo_scf_env%nspins
1660 :
1661 : CALL orthogonalize_mos(ket=almo_scf_env%matrix_t(ispin), &
1662 : overlap=almo_scf_env%matrix_sigma(ispin), &
1663 : metric=almo_scf_env%matrix_s(1), &
1664 : retain_locality=.FALSE., &
1665 : only_normalize=.FALSE., &
1666 : nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1667 : eps_filter=almo_scf_env%eps_filter, &
1668 : order_lanczos=almo_scf_env%order_lanczos, &
1669 : eps_lanczos=almo_scf_env%eps_lanczos, &
1670 8 : max_iter_lanczos=almo_scf_env%max_iter_lanczos)
1671 : END DO
1672 :
1673 4 : CALL construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals=.FALSE.)
1674 :
1675 4 : IF (almo_scf_env%opt_nlmo_pcg%opt_penalty%virtual_nlmos) THEN
1676 0 : CALL construct_virtuals(almo_scf_env)
1677 0 : CALL construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals=.TRUE.)
1678 : END IF
1679 :
1680 4 : IF (almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start > 0.0_dp) THEN
1681 2 : CALL nlmo_compactification(qs_env, almo_scf_env, almo_scf_env%matrix_t)
1682 : END IF
1683 :
1684 : END IF
1685 :
1686 122 : END SUBROUTINE construct_nlmos
1687 :
1688 : ! **************************************************************************************************
1689 : !> \brief Calls NLMO optimization
1690 : !> \param qs_env ...
1691 : !> \param almo_scf_env ...
1692 : !> \param virtuals ...
1693 : !> \par History
1694 : !> 2019.10 created [Ziling Luo]
1695 : !> \author Ziling Luo
1696 : ! **************************************************************************************************
1697 4 : SUBROUTINE construct_nlmos_wrapper(qs_env, almo_scf_env, virtuals)
1698 :
1699 : TYPE(qs_environment_type), POINTER :: qs_env
1700 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1701 : LOGICAL, INTENT(IN) :: virtuals
1702 :
1703 : REAL(KIND=dp) :: det_diff, prev_determinant
1704 :
1705 4 : almo_scf_env%overlap_determinant = 1.0
1706 : ! KEEP: initial_vol_coeff = almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength
1707 : almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength = &
1708 4 : -1.0_dp*almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength !NEW1
1709 :
1710 : ! loop over the strength of the orthogonalization penalty
1711 4 : prev_determinant = 10.0_dp
1712 10 : DO WHILE (almo_scf_env%overlap_determinant > almo_scf_env%opt_nlmo_pcg%opt_penalty%final_determinant)
1713 :
1714 8 : IF (.NOT. virtuals) THEN
1715 : CALL almo_scf_construct_nlmos(qs_env=qs_env, &
1716 : optimizer=almo_scf_env%opt_nlmo_pcg, &
1717 : matrix_s=almo_scf_env%matrix_s(1), &
1718 : matrix_mo_in=almo_scf_env%matrix_t, &
1719 : matrix_mo_out=almo_scf_env%matrix_t, &
1720 : template_matrix_sigma=almo_scf_env%matrix_sigma_inv, &
1721 : overlap_determinant=almo_scf_env%overlap_determinant, &
1722 : mat_distr_aos=almo_scf_env%mat_distr_aos, &
1723 : virtuals=virtuals, &
1724 8 : eps_filter=almo_scf_env%eps_filter)
1725 : ELSE
1726 : CALL almo_scf_construct_nlmos(qs_env=qs_env, &
1727 : optimizer=almo_scf_env%opt_nlmo_pcg, &
1728 : matrix_s=almo_scf_env%matrix_s(1), &
1729 : matrix_mo_in=almo_scf_env%matrix_v, &
1730 : matrix_mo_out=almo_scf_env%matrix_v, &
1731 : template_matrix_sigma=almo_scf_env%matrix_sigma_vv, &
1732 : overlap_determinant=almo_scf_env%overlap_determinant, &
1733 : mat_distr_aos=almo_scf_env%mat_distr_aos, &
1734 : virtuals=virtuals, &
1735 0 : eps_filter=almo_scf_env%eps_filter)
1736 :
1737 : END IF
1738 :
1739 8 : det_diff = prev_determinant - almo_scf_env%overlap_determinant
1740 : almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength = &
1741 : almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength/ &
1742 8 : ABS(almo_scf_env%opt_nlmo_pcg%opt_penalty%penalty_strength_dec_factor)
1743 :
1744 8 : IF (det_diff < almo_scf_env%opt_nlmo_pcg%opt_penalty%determinant_tolerance) THEN
1745 : EXIT
1746 : END IF
1747 4 : prev_determinant = almo_scf_env%overlap_determinant
1748 :
1749 : END DO
1750 :
1751 4 : END SUBROUTINE construct_nlmos_wrapper
1752 :
1753 : ! **************************************************************************************************
1754 : !> \brief Construct virtual orbitals
1755 : !> \param almo_scf_env ...
1756 : !> \par History
1757 : !> 2019.10 created [Ziling Luo]
1758 : !> \author Ziling Luo
1759 : ! **************************************************************************************************
1760 0 : SUBROUTINE construct_virtuals(almo_scf_env)
1761 :
1762 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1763 :
1764 : INTEGER :: ispin, n
1765 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1766 : TYPE(dbcsr_type) :: tempNV1, tempVOcc1, tempVOcc2, tempVV1, &
1767 : tempVV2
1768 :
1769 0 : DO ispin = 1, almo_scf_env%nspins
1770 :
1771 : CALL dbcsr_create(tempNV1, &
1772 : template=almo_scf_env%matrix_v(ispin), &
1773 0 : matrix_type=dbcsr_type_no_symmetry)
1774 : CALL dbcsr_create(tempVOcc1, &
1775 : template=almo_scf_env%matrix_vo(ispin), &
1776 0 : matrix_type=dbcsr_type_no_symmetry)
1777 : CALL dbcsr_create(tempVOcc2, &
1778 : template=almo_scf_env%matrix_vo(ispin), &
1779 0 : matrix_type=dbcsr_type_no_symmetry)
1780 : CALL dbcsr_create(tempVV1, &
1781 : template=almo_scf_env%matrix_sigma_vv(ispin), &
1782 0 : matrix_type=dbcsr_type_no_symmetry)
1783 : CALL dbcsr_create(tempVV2, &
1784 : template=almo_scf_env%matrix_sigma_vv(ispin), &
1785 0 : matrix_type=dbcsr_type_no_symmetry)
1786 :
1787 : ! Generate random virtual matrix
1788 : CALL dbcsr_init_random(almo_scf_env%matrix_v(ispin), &
1789 0 : keep_sparsity=.FALSE.)
1790 :
1791 : ! Project the orbital subspace out
1792 : CALL dbcsr_multiply("N", "N", 1.0_dp, &
1793 : almo_scf_env%matrix_s(1), &
1794 : almo_scf_env%matrix_v(ispin), &
1795 : 0.0_dp, tempNV1, &
1796 0 : filter_eps=almo_scf_env%eps_filter)
1797 :
1798 : CALL dbcsr_multiply("T", "N", 1.0_dp, &
1799 : tempNV1, &
1800 : almo_scf_env%matrix_t(ispin), &
1801 : 0.0_dp, tempVOcc1, &
1802 0 : filter_eps=almo_scf_env%eps_filter)
1803 :
1804 : CALL dbcsr_multiply("N", "N", 1.0_dp, &
1805 : tempVOcc1, &
1806 : almo_scf_env%matrix_sigma_inv(ispin), &
1807 : 0.0_dp, tempVOcc2, &
1808 0 : filter_eps=almo_scf_env%eps_filter)
1809 :
1810 : CALL dbcsr_multiply("N", "T", 1.0_dp, &
1811 : almo_scf_env%matrix_t(ispin), &
1812 : tempVOcc2, &
1813 : 0.0_dp, tempNV1, &
1814 0 : filter_eps=almo_scf_env%eps_filter)
1815 :
1816 0 : CALL dbcsr_add(almo_scf_env%matrix_v(ispin), tempNV1, 1.0_dp, -1.0_dp)
1817 :
1818 : ! compute VxV overlap
1819 : CALL dbcsr_multiply("N", "N", 1.0_dp, &
1820 : almo_scf_env%matrix_s(1), &
1821 : almo_scf_env%matrix_v(ispin), &
1822 : 0.0_dp, tempNV1, &
1823 0 : filter_eps=almo_scf_env%eps_filter)
1824 :
1825 : CALL dbcsr_multiply("T", "N", 1.0_dp, &
1826 : almo_scf_env%matrix_v(ispin), &
1827 : tempNV1, &
1828 : 0.0_dp, tempVV1, &
1829 0 : filter_eps=almo_scf_env%eps_filter)
1830 :
1831 : CALL orthogonalize_mos(ket=almo_scf_env%matrix_v(ispin), &
1832 : overlap=tempVV1, &
1833 : metric=almo_scf_env%matrix_s(1), &
1834 : retain_locality=.FALSE., &
1835 : only_normalize=.FALSE., &
1836 : nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1837 : eps_filter=almo_scf_env%eps_filter, &
1838 : order_lanczos=almo_scf_env%order_lanczos, &
1839 : eps_lanczos=almo_scf_env%eps_lanczos, &
1840 0 : max_iter_lanczos=almo_scf_env%max_iter_lanczos)
1841 :
1842 : ! compute VxV block of the KS matrix
1843 : CALL dbcsr_multiply("N", "N", 1.0_dp, &
1844 : almo_scf_env%matrix_ks(ispin), &
1845 : almo_scf_env%matrix_v(ispin), &
1846 : 0.0_dp, tempNV1, &
1847 0 : filter_eps=almo_scf_env%eps_filter)
1848 :
1849 : CALL dbcsr_multiply("T", "N", 1.0_dp, &
1850 : almo_scf_env%matrix_v(ispin), &
1851 : tempNV1, &
1852 : 0.0_dp, tempVV1, &
1853 0 : filter_eps=almo_scf_env%eps_filter)
1854 :
1855 0 : CALL dbcsr_get_info(tempVV1, nfullrows_total=n)
1856 0 : ALLOCATE (eigenvalues(n))
1857 : CALL cp_dbcsr_syevd(tempVV1, tempVV2, &
1858 : eigenvalues, &
1859 : para_env=almo_scf_env%para_env, &
1860 0 : blacs_env=almo_scf_env%blacs_env)
1861 0 : DEALLOCATE (eigenvalues)
1862 :
1863 : CALL dbcsr_multiply("N", "N", 1.0_dp, &
1864 : almo_scf_env%matrix_v(ispin), &
1865 : tempVV2, &
1866 : 0.0_dp, tempNV1, &
1867 0 : filter_eps=almo_scf_env%eps_filter)
1868 :
1869 0 : CALL dbcsr_copy(almo_scf_env%matrix_v(ispin), tempNV1)
1870 :
1871 0 : CALL dbcsr_release(tempNV1)
1872 0 : CALL dbcsr_release(tempVOcc1)
1873 0 : CALL dbcsr_release(tempVOcc2)
1874 0 : CALL dbcsr_release(tempVV1)
1875 0 : CALL dbcsr_release(tempVV2)
1876 :
1877 : END DO
1878 :
1879 0 : END SUBROUTINE construct_virtuals
1880 :
1881 : ! **************************************************************************************************
1882 : !> \brief Compactify (set small blocks to zero) orbitals
1883 : !> \param qs_env ...
1884 : !> \param almo_scf_env ...
1885 : !> \param matrix ...
1886 : !> \par History
1887 : !> 2019.10 created [Ziling Luo]
1888 : !> \author Ziling Luo
1889 : ! **************************************************************************************************
1890 2 : SUBROUTINE nlmo_compactification(qs_env, almo_scf_env, matrix)
1891 :
1892 : TYPE(qs_environment_type), POINTER :: qs_env
1893 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
1894 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:), &
1895 : INTENT(IN) :: matrix
1896 :
1897 : INTEGER :: iblock_col, iblock_col_size, iblock_row, &
1898 : iblock_row_size, icol, irow, ispin, &
1899 : Ncols, Nrows, nspins, unit_nr
1900 : LOGICAL :: element_by_element
1901 : REAL(KIND=dp) :: energy, eps_local, eps_start, &
1902 : max_element, spin_factor
1903 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occ, retained
1904 2 : REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p
1905 : TYPE(cp_logger_type), POINTER :: logger
1906 : TYPE(dbcsr_iterator_type) :: iter
1907 2 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_p_tmp, matrix_t_tmp
1908 : TYPE(mp_comm_type) :: group
1909 :
1910 : ! define the output_unit
1911 4 : logger => cp_get_default_logger()
1912 2 : IF (logger%para_env%is_source()) THEN
1913 1 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1914 : ELSE
1915 : unit_nr = -1
1916 : END IF
1917 :
1918 2 : nspins = SIZE(matrix)
1919 2 : element_by_element = .FALSE.
1920 :
1921 2 : IF (nspins == 1) THEN
1922 2 : spin_factor = 2.0_dp
1923 : ELSE
1924 0 : spin_factor = 1.0_dp
1925 : END IF
1926 :
1927 8 : ALLOCATE (matrix_t_tmp(nspins))
1928 6 : ALLOCATE (matrix_p_tmp(nspins))
1929 6 : ALLOCATE (retained(nspins))
1930 2 : ALLOCATE (occ(2))
1931 :
1932 4 : DO ispin = 1, nspins
1933 :
1934 : ! init temporary storage
1935 : CALL dbcsr_create(matrix_t_tmp(ispin), &
1936 : template=matrix(ispin), &
1937 2 : matrix_type=dbcsr_type_no_symmetry)
1938 2 : CALL dbcsr_copy(matrix_t_tmp(ispin), matrix(ispin))
1939 :
1940 : CALL dbcsr_create(matrix_p_tmp(ispin), &
1941 : template=almo_scf_env%matrix_p(ispin), &
1942 2 : matrix_type=dbcsr_type_no_symmetry)
1943 4 : CALL dbcsr_copy(matrix_p_tmp(ispin), almo_scf_env%matrix_p(ispin))
1944 :
1945 : END DO
1946 :
1947 2 : IF (unit_nr > 0) THEN
1948 1 : WRITE (unit_nr, *)
1949 : WRITE (unit_nr, '(T2,A)') &
1950 1 : "Energy dependence on the (block-by-block) filtering of the NLMO coefficients"
1951 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A13,A20,A20,A25)') &
1952 1 : "EPS filter", "Occupation Alpha", "Occupation Beta", "Energy"
1953 : END IF
1954 :
1955 2 : eps_start = almo_scf_env%opt_nlmo_pcg%opt_penalty%compactification_filter_start
1956 2 : eps_local = MAX(eps_start, 10E-14_dp)
1957 :
1958 8 : DO
1959 :
1960 10 : IF (eps_local > 0.11_dp) EXIT
1961 :
1962 16 : DO ispin = 1, nspins
1963 :
1964 8 : retained(ispin) = 0
1965 8 : CALL dbcsr_work_create(matrix_t_tmp(ispin), work_mutable=.TRUE.)
1966 8 : CALL dbcsr_iterator_start(iter, matrix_t_tmp(ispin))
1967 264 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1968 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, &
1969 256 : row_size=iblock_row_size, col_size=iblock_col_size)
1970 776 : DO icol = 1, iblock_col_size
1971 :
1972 256 : IF (element_by_element) THEN
1973 :
1974 : DO irow = 1, iblock_row_size
1975 : IF (ABS(data_p(irow, icol)) < eps_local) THEN
1976 : data_p(irow, icol) = 0.0_dp
1977 : ELSE
1978 : retained(ispin) = retained(ispin) + 1
1979 : END IF
1980 : END DO
1981 :
1982 : ELSE ! rows are blocked
1983 :
1984 512 : max_element = 0.0_dp
1985 2560 : DO irow = 1, iblock_row_size
1986 2560 : IF (ABS(data_p(irow, icol)) > max_element) THEN
1987 : max_element = ABS(data_p(irow, icol))
1988 : END IF
1989 : END DO
1990 512 : IF (max_element < eps_local) THEN
1991 155 : DO irow = 1, iblock_row_size
1992 155 : data_p(irow, icol) = 0.0_dp
1993 : END DO
1994 : ELSE
1995 481 : retained(ispin) = retained(ispin) + iblock_row_size
1996 : END IF
1997 :
1998 : END IF ! block rows?
1999 : END DO ! icol
2000 :
2001 : END DO ! iterator
2002 8 : CALL dbcsr_iterator_stop(iter)
2003 8 : CALL dbcsr_finalize(matrix_t_tmp(ispin))
2004 8 : CALL dbcsr_filter(matrix_t_tmp(ispin), eps_local)
2005 :
2006 : CALL dbcsr_get_info(matrix_t_tmp(ispin), group=group, &
2007 : nfullrows_total=Nrows, &
2008 8 : nfullcols_total=Ncols)
2009 8 : CALL group%sum(retained(ispin))
2010 :
2011 : !devide by the total no. elements
2012 8 : occ(ispin) = retained(ispin)/Nrows/Ncols
2013 :
2014 : ! compute the global projectors (for the density matrix)
2015 : CALL almo_scf_t_to_proj( &
2016 : t=matrix_t_tmp(ispin), &
2017 : p=matrix_p_tmp(ispin), &
2018 : eps_filter=almo_scf_env%eps_filter, &
2019 : orthog_orbs=.FALSE., &
2020 : nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
2021 : s=almo_scf_env%matrix_s(1), &
2022 : sigma=almo_scf_env%matrix_sigma(ispin), &
2023 : sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
2024 : use_guess=.FALSE., &
2025 : algorithm=almo_scf_env%sigma_inv_algorithm, &
2026 : inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
2027 : inverse_accelerator=almo_scf_env%order_lanczos, &
2028 : eps_lanczos=almo_scf_env%eps_lanczos, &
2029 : max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
2030 : para_env=almo_scf_env%para_env, &
2031 8 : blacs_env=almo_scf_env%blacs_env)
2032 :
2033 : ! compute dm from the projector(s)
2034 32 : CALL dbcsr_scale(matrix_p_tmp(ispin), spin_factor)
2035 :
2036 : END DO
2037 :
2038 : ! the KS matrix is updated outside the spin loop
2039 : CALL almo_dm_to_almo_ks(qs_env, &
2040 : matrix_p_tmp, &
2041 : almo_scf_env%matrix_ks, &
2042 : energy, &
2043 : almo_scf_env%eps_filter, &
2044 8 : almo_scf_env%mat_distr_aos)
2045 :
2046 8 : IF (nspins < 2) occ(2) = occ(1)
2047 8 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,E13.3,F20.10,F20.10,F25.15)') &
2048 4 : eps_local, occ(1), occ(2), energy
2049 :
2050 8 : eps_local = 2.0_dp*eps_local
2051 :
2052 : END DO
2053 :
2054 4 : DO ispin = 1, nspins
2055 :
2056 2 : CALL dbcsr_release(matrix_t_tmp(ispin))
2057 4 : CALL dbcsr_release(matrix_p_tmp(ispin))
2058 :
2059 : END DO
2060 :
2061 2 : DEALLOCATE (matrix_t_tmp)
2062 2 : DEALLOCATE (matrix_p_tmp)
2063 2 : DEALLOCATE (occ)
2064 2 : DEALLOCATE (retained)
2065 :
2066 2 : END SUBROUTINE nlmo_compactification
2067 :
2068 : ! *****************************************************************************
2069 : !> \brief after SCF we have the final density and KS matrices compute various
2070 : !> post-scf quantities
2071 : !> \param qs_env ...
2072 : !> \param almo_scf_env ...
2073 : !> \par History
2074 : !> 2015.03 created [Rustam Z. Khaliullin]
2075 : !> \author Rustam Z. Khaliullin
2076 : ! **************************************************************************************************
2077 122 : SUBROUTINE almo_scf_post(qs_env, almo_scf_env)
2078 : TYPE(qs_environment_type), POINTER :: qs_env
2079 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
2080 :
2081 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_post'
2082 :
2083 : INTEGER :: handle, ispin
2084 : TYPE(cp_fm_type), POINTER :: mo_coeff
2085 122 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w
2086 122 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_t_processed
2087 122 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2088 : TYPE(qs_scf_env_type), POINTER :: scf_env
2089 :
2090 122 : CALL timeset(routineN, handle)
2091 :
2092 : ! store matrices to speed up the next scf run
2093 122 : CALL almo_scf_store_extrapolation_data(almo_scf_env)
2094 :
2095 : ! orthogonalize orbitals before returning them to QS
2096 494 : ALLOCATE (matrix_t_processed(almo_scf_env%nspins))
2097 : !ALLOCATE (matrix_v_processed(almo_scf_env%nspins))
2098 :
2099 250 : DO ispin = 1, almo_scf_env%nspins
2100 :
2101 : CALL dbcsr_create(matrix_t_processed(ispin), &
2102 : template=almo_scf_env%matrix_t(ispin), &
2103 128 : matrix_type=dbcsr_type_no_symmetry)
2104 :
2105 : CALL dbcsr_copy(matrix_t_processed(ispin), &
2106 128 : almo_scf_env%matrix_t(ispin))
2107 :
2108 250 : IF (almo_scf_env%return_orthogonalized_mos) THEN
2109 :
2110 : CALL orthogonalize_mos(ket=matrix_t_processed(ispin), &
2111 : overlap=almo_scf_env%matrix_sigma(ispin), &
2112 : metric=almo_scf_env%matrix_s(1), &
2113 : retain_locality=.FALSE., &
2114 : only_normalize=.FALSE., &
2115 : nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
2116 : eps_filter=almo_scf_env%eps_filter, &
2117 : order_lanczos=almo_scf_env%order_lanczos, &
2118 : eps_lanczos=almo_scf_env%eps_lanczos, &
2119 : max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
2120 106 : smear=almo_scf_env%smear)
2121 : END IF
2122 :
2123 : END DO
2124 :
2125 : !! RS-WARNING: If smearing ALMO is requested, rescaled fully-occupied orbitals are returned to QS
2126 : !! RS-WARNING: Beware that QS will not be informed about electronic entropy.
2127 : !! If you want a quick and dirty transfer to QS energy, uncomment the following hack:
2128 : !! IF (almo_scf_env%smear) THEN
2129 : !! qs_env%energy%kTS = 0.0_dp
2130 : !! DO ispin = 1, almo_scf_env%nspins
2131 : !! qs_env%energy%kTS = qs_env%energy%kTS + almo_scf_env%kTS(ispin)
2132 : !! END DO
2133 : !! END IF
2134 :
2135 : ! return orbitals to QS
2136 122 : NULLIFY (mos, mo_coeff, scf_env)
2137 :
2138 122 : CALL get_qs_env(qs_env, mos=mos, scf_env=scf_env)
2139 :
2140 250 : DO ispin = 1, almo_scf_env%nspins
2141 :
2142 : ! Currently only fm version of mo_set is usable.
2143 : ! First transform the matrix_t to fm version
2144 128 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
2145 128 : CALL copy_dbcsr_to_fm(matrix_t_processed(ispin), mo_coeff)
2146 250 : CALL dbcsr_release(matrix_t_processed(ispin))
2147 : END DO
2148 250 : DO ispin = 1, almo_scf_env%nspins
2149 250 : CALL dbcsr_release(matrix_t_processed(ispin))
2150 : END DO
2151 122 : DEALLOCATE (matrix_t_processed)
2152 :
2153 : ! calculate post scf properties
2154 :
2155 122 : CALL almo_post_scf_compute_properties(qs_env)
2156 :
2157 : ! compute the W matrix if associated
2158 122 : IF (almo_scf_env%calc_forces) THEN
2159 66 : CALL get_qs_env(qs_env, matrix_w=matrix_w)
2160 66 : IF (ASSOCIATED(matrix_w)) THEN
2161 66 : CALL calculate_w_matrix_almo(matrix_w, almo_scf_env)
2162 : ELSE
2163 0 : CPABORT("Matrix W is needed but not associated")
2164 : END IF
2165 : END IF
2166 :
2167 122 : CALL timestop(handle)
2168 :
2169 122 : END SUBROUTINE almo_scf_post
2170 :
2171 : ! **************************************************************************************************
2172 : !> \brief create various matrices
2173 : !> \param almo_scf_env ...
2174 : !> \param matrix_s0 ...
2175 : !> \par History
2176 : !> 2011.07 created [Rustam Z Khaliullin]
2177 : !> \author Rustam Z Khaliullin
2178 : ! **************************************************************************************************
2179 122 : SUBROUTINE almo_scf_env_create_matrices(almo_scf_env, matrix_s0)
2180 :
2181 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
2182 : TYPE(dbcsr_type), INTENT(IN) :: matrix_s0
2183 :
2184 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_env_create_matrices'
2185 :
2186 : INTEGER :: handle, ispin, nspins
2187 :
2188 122 : CALL timeset(routineN, handle)
2189 :
2190 122 : nspins = almo_scf_env%nspins
2191 :
2192 : ! AO overlap matrix and its various functions
2193 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s(1), &
2194 : matrix_qs=matrix_s0, &
2195 : almo_scf_env=almo_scf_env, &
2196 : name_new="S", &
2197 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_aobasis], &
2198 : symmetry_new=dbcsr_type_symmetric, &
2199 : spin_key=0, &
2200 122 : init_domains=.FALSE.)
2201 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk(1), &
2202 : matrix_qs=matrix_s0, &
2203 : almo_scf_env=almo_scf_env, &
2204 : name_new="S_BLK", &
2205 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_aobasis], &
2206 : symmetry_new=dbcsr_type_symmetric, &
2207 : spin_key=0, &
2208 122 : init_domains=.TRUE.)
2209 122 : IF (almo_scf_env%almo_update_algorithm == almo_scf_diag) THEN
2210 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk_sqrt_inv(1), &
2211 : matrix_qs=matrix_s0, &
2212 : almo_scf_env=almo_scf_env, &
2213 : name_new="S_BLK_SQRT_INV", &
2214 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_aobasis], &
2215 : symmetry_new=dbcsr_type_symmetric, &
2216 : spin_key=0, &
2217 76 : init_domains=.TRUE.)
2218 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk_sqrt(1), &
2219 : matrix_qs=matrix_s0, &
2220 : almo_scf_env=almo_scf_env, &
2221 : name_new="S_BLK_SQRT", &
2222 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_aobasis], &
2223 : symmetry_new=dbcsr_type_symmetric, &
2224 : spin_key=0, &
2225 76 : init_domains=.TRUE.)
2226 46 : ELSE IF (almo_scf_env%almo_update_algorithm == almo_scf_dm_sign) THEN
2227 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_s_blk_inv(1), &
2228 : matrix_qs=matrix_s0, &
2229 : almo_scf_env=almo_scf_env, &
2230 : name_new="S_BLK_INV", &
2231 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_aobasis], &
2232 : symmetry_new=dbcsr_type_symmetric, &
2233 : spin_key=0, &
2234 0 : init_domains=.TRUE.)
2235 : END IF
2236 :
2237 : ! MO coeff matrices and their derivatives
2238 494 : ALLOCATE (almo_scf_env%matrix_t_blk(nspins))
2239 372 : ALLOCATE (almo_scf_env%quench_t_blk(nspins))
2240 372 : ALLOCATE (almo_scf_env%matrix_err_blk(nspins))
2241 372 : ALLOCATE (almo_scf_env%matrix_err_xx(nspins))
2242 372 : ALLOCATE (almo_scf_env%matrix_sigma(nspins))
2243 372 : ALLOCATE (almo_scf_env%matrix_sigma_inv(nspins))
2244 372 : ALLOCATE (almo_scf_env%matrix_sigma_sqrt(nspins))
2245 372 : ALLOCATE (almo_scf_env%matrix_sigma_sqrt_inv(nspins))
2246 372 : ALLOCATE (almo_scf_env%matrix_sigma_blk(nspins))
2247 372 : ALLOCATE (almo_scf_env%matrix_sigma_inv_0deloc(nspins))
2248 372 : ALLOCATE (almo_scf_env%matrix_t(nspins))
2249 372 : ALLOCATE (almo_scf_env%matrix_t_tr(nspins))
2250 250 : DO ispin = 1, nspins
2251 : ! create the blocked quencher
2252 : CALL matrix_almo_create(matrix_new=almo_scf_env%quench_t_blk(ispin), &
2253 : matrix_qs=matrix_s0, &
2254 : almo_scf_env=almo_scf_env, &
2255 : name_new="Q_BLK", &
2256 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_occ], &
2257 : symmetry_new=dbcsr_type_no_symmetry, &
2258 : spin_key=ispin, &
2259 128 : init_domains=.TRUE.)
2260 : ! create ALMO coefficient matrix
2261 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_t_blk(ispin), &
2262 : matrix_qs=matrix_s0, &
2263 : almo_scf_env=almo_scf_env, &
2264 : name_new="T_BLK", &
2265 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_occ], &
2266 : symmetry_new=dbcsr_type_no_symmetry, &
2267 : spin_key=ispin, &
2268 128 : init_domains=.TRUE.)
2269 : ! create the error matrix
2270 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_err_blk(ispin), &
2271 : matrix_qs=matrix_s0, &
2272 : almo_scf_env=almo_scf_env, &
2273 : name_new="ERR_BLK", &
2274 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_aobasis], &
2275 : symmetry_new=dbcsr_type_no_symmetry, &
2276 : spin_key=ispin, &
2277 128 : init_domains=.TRUE.)
2278 : ! create the error matrix for the quenched ALMOs
2279 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_err_xx(ispin), &
2280 : matrix_qs=matrix_s0, &
2281 : almo_scf_env=almo_scf_env, &
2282 : name_new="ERR_XX", &
2283 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_occ], &
2284 : symmetry_new=dbcsr_type_no_symmetry, &
2285 : spin_key=ispin, &
2286 128 : init_domains=.FALSE.)
2287 : ! create a matrix with dimensions of a transposed mo coefficient matrix
2288 : ! it might be necessary to perform the correction step using cayley
2289 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_t_tr(ispin), &
2290 : matrix_qs=matrix_s0, &
2291 : almo_scf_env=almo_scf_env, &
2292 : name_new="T_TR", &
2293 : size_keys=[almo_mat_dim_occ, almo_mat_dim_aobasis], &
2294 : symmetry_new=dbcsr_type_no_symmetry, &
2295 : spin_key=ispin, &
2296 128 : init_domains=.FALSE.)
2297 : ! create mo overlap matrix
2298 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma(ispin), &
2299 : matrix_qs=matrix_s0, &
2300 : almo_scf_env=almo_scf_env, &
2301 : name_new="SIG", &
2302 : size_keys=[almo_mat_dim_occ, almo_mat_dim_occ], &
2303 : symmetry_new=dbcsr_type_symmetric, &
2304 : spin_key=ispin, &
2305 128 : init_domains=.FALSE.)
2306 : ! create blocked mo overlap matrix
2307 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_blk(ispin), &
2308 : matrix_qs=matrix_s0, &
2309 : almo_scf_env=almo_scf_env, &
2310 : name_new="SIG_BLK", &
2311 : size_keys=[almo_mat_dim_occ, almo_mat_dim_occ], &
2312 : symmetry_new=dbcsr_type_symmetric, &
2313 : spin_key=ispin, &
2314 128 : init_domains=.TRUE.)
2315 : ! create blocked inverse mo overlap matrix
2316 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
2317 : matrix_qs=matrix_s0, &
2318 : almo_scf_env=almo_scf_env, &
2319 : name_new="SIGINV_BLK", &
2320 : size_keys=[almo_mat_dim_occ, almo_mat_dim_occ], &
2321 : symmetry_new=dbcsr_type_symmetric, &
2322 : spin_key=ispin, &
2323 128 : init_domains=.TRUE.)
2324 : ! create inverse mo overlap matrix
2325 : CALL matrix_almo_create( &
2326 : matrix_new=almo_scf_env%matrix_sigma_inv(ispin), &
2327 : matrix_qs=matrix_s0, &
2328 : almo_scf_env=almo_scf_env, &
2329 : name_new="SIGINV", &
2330 : size_keys=[almo_mat_dim_occ, almo_mat_dim_occ], &
2331 : symmetry_new=dbcsr_type_symmetric, &
2332 : spin_key=ispin, &
2333 128 : init_domains=.FALSE.)
2334 : ! create various templates that will be necessary later
2335 : CALL matrix_almo_create( &
2336 : matrix_new=almo_scf_env%matrix_t(ispin), &
2337 : matrix_qs=matrix_s0, &
2338 : almo_scf_env=almo_scf_env, &
2339 : name_new="T", &
2340 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_occ], &
2341 : symmetry_new=dbcsr_type_no_symmetry, &
2342 : spin_key=ispin, &
2343 128 : init_domains=.FALSE.)
2344 : CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt(ispin), &
2345 : template=almo_scf_env%matrix_sigma(ispin), &
2346 128 : matrix_type=dbcsr_type_no_symmetry)
2347 : CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
2348 : template=almo_scf_env%matrix_sigma(ispin), &
2349 250 : matrix_type=dbcsr_type_no_symmetry)
2350 : END DO
2351 :
2352 : ! create virtual orbitals if necessary
2353 122 : IF (almo_scf_env%need_virtuals) THEN
2354 372 : ALLOCATE (almo_scf_env%matrix_v_blk(nspins))
2355 372 : ALLOCATE (almo_scf_env%matrix_v_full_blk(nspins))
2356 372 : ALLOCATE (almo_scf_env%matrix_v(nspins))
2357 372 : ALLOCATE (almo_scf_env%matrix_vo(nspins))
2358 372 : ALLOCATE (almo_scf_env%matrix_x(nspins))
2359 372 : ALLOCATE (almo_scf_env%matrix_ov(nspins))
2360 372 : ALLOCATE (almo_scf_env%matrix_ov_full(nspins))
2361 372 : ALLOCATE (almo_scf_env%matrix_sigma_vv(nspins))
2362 372 : ALLOCATE (almo_scf_env%matrix_sigma_vv_blk(nspins))
2363 372 : ALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt(nspins))
2364 372 : ALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt_inv(nspins))
2365 372 : ALLOCATE (almo_scf_env%matrix_vv_full_blk(nspins))
2366 :
2367 122 : IF (almo_scf_env%deloc_truncate_virt /= virt_full) THEN
2368 0 : ALLOCATE (almo_scf_env%matrix_k_blk(nspins))
2369 0 : ALLOCATE (almo_scf_env%matrix_k_blk_ones(nspins))
2370 0 : ALLOCATE (almo_scf_env%matrix_k_tr(nspins))
2371 0 : ALLOCATE (almo_scf_env%matrix_v_disc(nspins))
2372 0 : ALLOCATE (almo_scf_env%matrix_v_disc_blk(nspins))
2373 0 : ALLOCATE (almo_scf_env%matrix_ov_disc(nspins))
2374 0 : ALLOCATE (almo_scf_env%matrix_vv_disc_blk(nspins))
2375 0 : ALLOCATE (almo_scf_env%matrix_vv_disc(nspins))
2376 0 : ALLOCATE (almo_scf_env%opt_k_t_dd(nspins))
2377 0 : ALLOCATE (almo_scf_env%opt_k_t_rr(nspins))
2378 0 : ALLOCATE (almo_scf_env%opt_k_denom(nspins))
2379 : END IF
2380 :
2381 250 : DO ispin = 1, nspins
2382 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_full_blk(ispin), &
2383 : matrix_qs=matrix_s0, &
2384 : almo_scf_env=almo_scf_env, &
2385 : name_new="V_FULL_BLK", &
2386 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_virt_full], &
2387 : symmetry_new=dbcsr_type_no_symmetry, &
2388 : spin_key=ispin, &
2389 128 : init_domains=.FALSE.)
2390 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_blk(ispin), &
2391 : matrix_qs=matrix_s0, &
2392 : almo_scf_env=almo_scf_env, &
2393 : name_new="V_BLK", &
2394 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_virt], &
2395 : symmetry_new=dbcsr_type_no_symmetry, &
2396 : spin_key=ispin, &
2397 128 : init_domains=.FALSE.)
2398 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v(ispin), &
2399 : matrix_qs=matrix_s0, &
2400 : almo_scf_env=almo_scf_env, &
2401 : name_new="V", &
2402 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_virt], &
2403 : symmetry_new=dbcsr_type_no_symmetry, &
2404 : spin_key=ispin, &
2405 128 : init_domains=.FALSE.)
2406 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ov_full(ispin), &
2407 : matrix_qs=matrix_s0, &
2408 : almo_scf_env=almo_scf_env, &
2409 : name_new="OV_FULL", &
2410 : size_keys=[almo_mat_dim_occ, almo_mat_dim_virt_full], &
2411 : symmetry_new=dbcsr_type_no_symmetry, &
2412 : spin_key=ispin, &
2413 128 : init_domains=.FALSE.)
2414 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ov(ispin), &
2415 : matrix_qs=matrix_s0, &
2416 : almo_scf_env=almo_scf_env, &
2417 : name_new="OV", &
2418 : size_keys=[almo_mat_dim_occ, almo_mat_dim_virt], &
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_vo(ispin), &
2423 : matrix_qs=matrix_s0, &
2424 : almo_scf_env=almo_scf_env, &
2425 : name_new="VO", &
2426 : size_keys=[almo_mat_dim_virt, almo_mat_dim_occ], &
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_x(ispin), &
2431 : matrix_qs=matrix_s0, &
2432 : almo_scf_env=almo_scf_env, &
2433 : name_new="VO", &
2434 : size_keys=[almo_mat_dim_virt, almo_mat_dim_occ], &
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_sigma_vv(ispin), &
2439 : matrix_qs=matrix_s0, &
2440 : almo_scf_env=almo_scf_env, &
2441 : name_new="SIG_VV", &
2442 : size_keys=[almo_mat_dim_virt, almo_mat_dim_virt], &
2443 : symmetry_new=dbcsr_type_symmetric, &
2444 : spin_key=ispin, &
2445 128 : init_domains=.FALSE.)
2446 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vv_full_blk(ispin), &
2447 : matrix_qs=matrix_s0, &
2448 : almo_scf_env=almo_scf_env, &
2449 : name_new="VV_FULL_BLK", &
2450 : size_keys=[almo_mat_dim_virt_full, almo_mat_dim_virt_full], &
2451 : symmetry_new=dbcsr_type_no_symmetry, &
2452 : spin_key=ispin, &
2453 128 : init_domains=.TRUE.)
2454 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_sigma_vv_blk(ispin), &
2455 : matrix_qs=matrix_s0, &
2456 : almo_scf_env=almo_scf_env, &
2457 : name_new="SIG_VV_BLK", &
2458 : size_keys=[almo_mat_dim_virt, almo_mat_dim_virt], &
2459 : symmetry_new=dbcsr_type_symmetric, &
2460 : spin_key=ispin, &
2461 128 : init_domains=.TRUE.)
2462 : CALL dbcsr_create(almo_scf_env%matrix_sigma_vv_sqrt(ispin), &
2463 : template=almo_scf_env%matrix_sigma_vv(ispin), &
2464 128 : matrix_type=dbcsr_type_no_symmetry)
2465 : CALL dbcsr_create(almo_scf_env%matrix_sigma_vv_sqrt_inv(ispin), &
2466 : template=almo_scf_env%matrix_sigma_vv(ispin), &
2467 128 : matrix_type=dbcsr_type_no_symmetry)
2468 :
2469 250 : IF (almo_scf_env%deloc_truncate_virt /= virt_full) THEN
2470 : CALL matrix_almo_create(matrix_new=almo_scf_env%opt_k_t_rr(ispin), &
2471 : matrix_qs=matrix_s0, &
2472 : almo_scf_env=almo_scf_env, &
2473 : name_new="OPT_K_U_RR", &
2474 : size_keys=[almo_mat_dim_virt, almo_mat_dim_virt], &
2475 : symmetry_new=dbcsr_type_no_symmetry, &
2476 : spin_key=ispin, &
2477 0 : init_domains=.FALSE.)
2478 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vv_disc(ispin), &
2479 : matrix_qs=matrix_s0, &
2480 : almo_scf_env=almo_scf_env, &
2481 : name_new="VV_DISC", &
2482 : size_keys=[almo_mat_dim_virt_disc, almo_mat_dim_virt_disc], &
2483 : symmetry_new=dbcsr_type_symmetric, &
2484 : spin_key=ispin, &
2485 0 : init_domains=.FALSE.)
2486 : CALL matrix_almo_create(matrix_new=almo_scf_env%opt_k_t_dd(ispin), &
2487 : matrix_qs=matrix_s0, &
2488 : almo_scf_env=almo_scf_env, &
2489 : name_new="OPT_K_U_DD", &
2490 : size_keys=[almo_mat_dim_virt_disc, almo_mat_dim_virt_disc], &
2491 : symmetry_new=dbcsr_type_no_symmetry, &
2492 : spin_key=ispin, &
2493 0 : init_domains=.FALSE.)
2494 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_vv_disc_blk(ispin), &
2495 : matrix_qs=matrix_s0, &
2496 : almo_scf_env=almo_scf_env, &
2497 : name_new="VV_DISC_BLK", &
2498 : size_keys=[almo_mat_dim_virt_disc, almo_mat_dim_virt_disc], &
2499 : symmetry_new=dbcsr_type_symmetric, &
2500 : spin_key=ispin, &
2501 0 : init_domains=.TRUE.)
2502 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_k_blk(ispin), &
2503 : matrix_qs=matrix_s0, &
2504 : almo_scf_env=almo_scf_env, &
2505 : name_new="K_BLK", &
2506 : size_keys=[almo_mat_dim_virt_disc, almo_mat_dim_virt], &
2507 : symmetry_new=dbcsr_type_no_symmetry, &
2508 : spin_key=ispin, &
2509 0 : init_domains=.TRUE.)
2510 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_k_blk_ones(ispin), &
2511 : matrix_qs=matrix_s0, &
2512 : almo_scf_env=almo_scf_env, &
2513 : name_new="K_BLK_1", &
2514 : size_keys=[almo_mat_dim_virt_disc, almo_mat_dim_virt], &
2515 : symmetry_new=dbcsr_type_no_symmetry, &
2516 : spin_key=ispin, &
2517 0 : init_domains=.TRUE.)
2518 : CALL matrix_almo_create(matrix_new=almo_scf_env%opt_k_denom(ispin), &
2519 : matrix_qs=matrix_s0, &
2520 : almo_scf_env=almo_scf_env, &
2521 : name_new="OPT_K_DENOM", &
2522 : size_keys=[almo_mat_dim_virt_disc, almo_mat_dim_virt], &
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_k_tr(ispin), &
2527 : matrix_qs=matrix_s0, &
2528 : almo_scf_env=almo_scf_env, &
2529 : name_new="K_TR", &
2530 : size_keys=[almo_mat_dim_virt, almo_mat_dim_virt_disc], &
2531 : symmetry_new=dbcsr_type_no_symmetry, &
2532 : spin_key=ispin, &
2533 0 : init_domains=.FALSE.)
2534 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_disc_blk(ispin), &
2535 : matrix_qs=matrix_s0, &
2536 : almo_scf_env=almo_scf_env, &
2537 : name_new="V_DISC_BLK", &
2538 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_virt_disc], &
2539 : symmetry_new=dbcsr_type_no_symmetry, &
2540 : spin_key=ispin, &
2541 0 : init_domains=.FALSE.)
2542 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_v_disc(ispin), &
2543 : matrix_qs=matrix_s0, &
2544 : almo_scf_env=almo_scf_env, &
2545 : name_new="V_DISC", &
2546 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_virt_disc], &
2547 : symmetry_new=dbcsr_type_no_symmetry, &
2548 : spin_key=ispin, &
2549 0 : init_domains=.FALSE.)
2550 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ov_disc(ispin), &
2551 : matrix_qs=matrix_s0, &
2552 : almo_scf_env=almo_scf_env, &
2553 : name_new="OV_DISC", &
2554 : size_keys=[almo_mat_dim_occ, almo_mat_dim_virt_disc], &
2555 : symmetry_new=dbcsr_type_no_symmetry, &
2556 : spin_key=ispin, &
2557 0 : init_domains=.FALSE.)
2558 :
2559 : END IF ! end need_discarded_virtuals
2560 :
2561 : END DO ! spin
2562 : END IF
2563 :
2564 : ! create matrices of orbital energies if necessary
2565 122 : IF (almo_scf_env%need_orbital_energies) THEN
2566 372 : ALLOCATE (almo_scf_env%matrix_eoo(nspins))
2567 372 : ALLOCATE (almo_scf_env%matrix_evv_full(nspins))
2568 250 : DO ispin = 1, nspins
2569 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_eoo(ispin), &
2570 : matrix_qs=matrix_s0, &
2571 : almo_scf_env=almo_scf_env, &
2572 : name_new="E_OCC", &
2573 : size_keys=[almo_mat_dim_occ, almo_mat_dim_occ], &
2574 : symmetry_new=dbcsr_type_no_symmetry, &
2575 : spin_key=ispin, &
2576 128 : init_domains=.FALSE.)
2577 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_evv_full(ispin), &
2578 : matrix_qs=matrix_s0, &
2579 : almo_scf_env=almo_scf_env, &
2580 : name_new="E_VIRT", &
2581 : size_keys=[almo_mat_dim_virt_full, almo_mat_dim_virt_full], &
2582 : symmetry_new=dbcsr_type_no_symmetry, &
2583 : spin_key=ispin, &
2584 250 : init_domains=.FALSE.)
2585 : END DO
2586 : END IF
2587 :
2588 : ! Density and KS matrices
2589 372 : ALLOCATE (almo_scf_env%matrix_p(nspins))
2590 372 : ALLOCATE (almo_scf_env%matrix_p_blk(nspins))
2591 372 : ALLOCATE (almo_scf_env%matrix_ks(nspins))
2592 372 : ALLOCATE (almo_scf_env%matrix_ks_blk(nspins))
2593 122 : IF (almo_scf_env%need_previous_ks) THEN
2594 372 : ALLOCATE (almo_scf_env%matrix_ks_0deloc(nspins))
2595 : END IF
2596 250 : DO ispin = 1, nspins
2597 : ! RZK-warning copy with symmery but remember that this might cause problems
2598 : CALL dbcsr_create(almo_scf_env%matrix_p(ispin), &
2599 : template=almo_scf_env%matrix_s(1), &
2600 128 : matrix_type=dbcsr_type_symmetric)
2601 : CALL dbcsr_create(almo_scf_env%matrix_ks(ispin), &
2602 : template=almo_scf_env%matrix_s(1), &
2603 128 : matrix_type=dbcsr_type_symmetric)
2604 128 : IF (almo_scf_env%need_previous_ks) THEN
2605 : CALL dbcsr_create(almo_scf_env%matrix_ks_0deloc(ispin), &
2606 : template=almo_scf_env%matrix_s(1), &
2607 128 : matrix_type=dbcsr_type_symmetric)
2608 : END IF
2609 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_p_blk(ispin), &
2610 : matrix_qs=matrix_s0, &
2611 : almo_scf_env=almo_scf_env, &
2612 : name_new="P_BLK", &
2613 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_aobasis], &
2614 : symmetry_new=dbcsr_type_symmetric, &
2615 : spin_key=ispin, &
2616 128 : init_domains=.TRUE.)
2617 : CALL matrix_almo_create(matrix_new=almo_scf_env%matrix_ks_blk(ispin), &
2618 : matrix_qs=matrix_s0, &
2619 : almo_scf_env=almo_scf_env, &
2620 : name_new="KS_BLK", &
2621 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_aobasis], &
2622 : symmetry_new=dbcsr_type_symmetric, &
2623 : spin_key=ispin, &
2624 250 : init_domains=.TRUE.)
2625 : END DO
2626 :
2627 122 : CALL timestop(handle)
2628 :
2629 122 : END SUBROUTINE almo_scf_env_create_matrices
2630 :
2631 : ! **************************************************************************************************
2632 : !> \brief clean up procedures for almo scf
2633 : !> \param almo_scf_env ...
2634 : !> \par History
2635 : !> 2011.06 created [Rustam Z Khaliullin]
2636 : !> 2018.09 smearing support [Ruben Staub]
2637 : !> \author Rustam Z Khaliullin
2638 : ! **************************************************************************************************
2639 122 : SUBROUTINE almo_scf_clean_up(almo_scf_env)
2640 :
2641 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
2642 :
2643 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_clean_up'
2644 :
2645 : INTEGER :: handle, ispin, unit_nr
2646 : TYPE(cp_logger_type), POINTER :: logger
2647 :
2648 122 : CALL timeset(routineN, handle)
2649 :
2650 : ! get a useful output_unit
2651 122 : logger => cp_get_default_logger()
2652 122 : IF (logger%para_env%is_source()) THEN
2653 61 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
2654 : ELSE
2655 : unit_nr = -1
2656 : END IF
2657 :
2658 : ! release matrices
2659 122 : CALL dbcsr_release(almo_scf_env%matrix_s(1))
2660 122 : CALL dbcsr_release(almo_scf_env%matrix_s_blk(1))
2661 122 : IF (almo_scf_env%almo_update_algorithm == almo_scf_diag) THEN
2662 76 : CALL dbcsr_release(almo_scf_env%matrix_s_blk_sqrt_inv(1))
2663 76 : CALL dbcsr_release(almo_scf_env%matrix_s_blk_sqrt(1))
2664 46 : ELSE IF (almo_scf_env%almo_update_algorithm == almo_scf_dm_sign) THEN
2665 0 : CALL dbcsr_release(almo_scf_env%matrix_s_blk_inv(1))
2666 : END IF
2667 250 : DO ispin = 1, almo_scf_env%nspins
2668 128 : CALL dbcsr_release(almo_scf_env%quench_t(ispin))
2669 128 : CALL dbcsr_release(almo_scf_env%quench_t_blk(ispin))
2670 128 : CALL dbcsr_release(almo_scf_env%matrix_t_blk(ispin))
2671 128 : CALL dbcsr_release(almo_scf_env%matrix_err_blk(ispin))
2672 128 : CALL dbcsr_release(almo_scf_env%matrix_err_xx(ispin))
2673 128 : CALL dbcsr_release(almo_scf_env%matrix_t_tr(ispin))
2674 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma(ispin))
2675 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma_blk(ispin))
2676 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma_inv_0deloc(ispin))
2677 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma_inv(ispin))
2678 128 : CALL dbcsr_release(almo_scf_env%matrix_t(ispin))
2679 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma_sqrt(ispin))
2680 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma_sqrt_inv(ispin))
2681 128 : CALL dbcsr_release(almo_scf_env%matrix_p(ispin))
2682 128 : CALL dbcsr_release(almo_scf_env%matrix_ks(ispin))
2683 128 : CALL dbcsr_release(almo_scf_env%matrix_p_blk(ispin))
2684 128 : CALL dbcsr_release(almo_scf_env%matrix_ks_blk(ispin))
2685 128 : IF (almo_scf_env%need_previous_ks) THEN
2686 128 : CALL dbcsr_release(almo_scf_env%matrix_ks_0deloc(ispin))
2687 : END IF
2688 128 : IF (almo_scf_env%need_virtuals) THEN
2689 128 : CALL dbcsr_release(almo_scf_env%matrix_v_blk(ispin))
2690 128 : CALL dbcsr_release(almo_scf_env%matrix_v_full_blk(ispin))
2691 128 : CALL dbcsr_release(almo_scf_env%matrix_v(ispin))
2692 128 : CALL dbcsr_release(almo_scf_env%matrix_vo(ispin))
2693 128 : CALL dbcsr_release(almo_scf_env%matrix_x(ispin))
2694 128 : CALL dbcsr_release(almo_scf_env%matrix_ov(ispin))
2695 128 : CALL dbcsr_release(almo_scf_env%matrix_ov_full(ispin))
2696 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma_vv(ispin))
2697 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma_vv_blk(ispin))
2698 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma_vv_sqrt(ispin))
2699 128 : CALL dbcsr_release(almo_scf_env%matrix_sigma_vv_sqrt_inv(ispin))
2700 128 : CALL dbcsr_release(almo_scf_env%matrix_vv_full_blk(ispin))
2701 128 : IF (almo_scf_env%deloc_truncate_virt /= virt_full) THEN
2702 0 : CALL dbcsr_release(almo_scf_env%matrix_k_tr(ispin))
2703 0 : CALL dbcsr_release(almo_scf_env%matrix_k_blk(ispin))
2704 0 : CALL dbcsr_release(almo_scf_env%matrix_k_blk_ones(ispin))
2705 0 : CALL dbcsr_release(almo_scf_env%matrix_v_disc(ispin))
2706 0 : CALL dbcsr_release(almo_scf_env%matrix_v_disc_blk(ispin))
2707 0 : CALL dbcsr_release(almo_scf_env%matrix_ov_disc(ispin))
2708 0 : CALL dbcsr_release(almo_scf_env%matrix_vv_disc_blk(ispin))
2709 0 : CALL dbcsr_release(almo_scf_env%matrix_vv_disc(ispin))
2710 0 : CALL dbcsr_release(almo_scf_env%opt_k_t_dd(ispin))
2711 0 : CALL dbcsr_release(almo_scf_env%opt_k_t_rr(ispin))
2712 0 : CALL dbcsr_release(almo_scf_env%opt_k_denom(ispin))
2713 : END IF
2714 : END IF
2715 250 : IF (almo_scf_env%need_orbital_energies) THEN
2716 128 : CALL dbcsr_release(almo_scf_env%matrix_eoo(ispin))
2717 128 : CALL dbcsr_release(almo_scf_env%matrix_evv_full(ispin))
2718 : END IF
2719 : END DO
2720 :
2721 : ! deallocate matrices
2722 122 : DEALLOCATE (almo_scf_env%matrix_p)
2723 122 : DEALLOCATE (almo_scf_env%matrix_p_blk)
2724 122 : DEALLOCATE (almo_scf_env%matrix_ks)
2725 122 : DEALLOCATE (almo_scf_env%matrix_ks_blk)
2726 122 : DEALLOCATE (almo_scf_env%matrix_t_blk)
2727 122 : DEALLOCATE (almo_scf_env%matrix_err_blk)
2728 122 : DEALLOCATE (almo_scf_env%matrix_err_xx)
2729 122 : DEALLOCATE (almo_scf_env%matrix_t)
2730 122 : DEALLOCATE (almo_scf_env%matrix_t_tr)
2731 122 : DEALLOCATE (almo_scf_env%matrix_sigma)
2732 122 : DEALLOCATE (almo_scf_env%matrix_sigma_blk)
2733 122 : DEALLOCATE (almo_scf_env%matrix_sigma_inv_0deloc)
2734 122 : DEALLOCATE (almo_scf_env%matrix_sigma_sqrt)
2735 122 : DEALLOCATE (almo_scf_env%matrix_sigma_sqrt_inv)
2736 122 : DEALLOCATE (almo_scf_env%matrix_sigma_inv)
2737 122 : DEALLOCATE (almo_scf_env%quench_t)
2738 122 : DEALLOCATE (almo_scf_env%quench_t_blk)
2739 122 : IF (almo_scf_env%need_virtuals) THEN
2740 122 : DEALLOCATE (almo_scf_env%matrix_v_blk)
2741 122 : DEALLOCATE (almo_scf_env%matrix_v_full_blk)
2742 122 : DEALLOCATE (almo_scf_env%matrix_v)
2743 122 : DEALLOCATE (almo_scf_env%matrix_vo)
2744 122 : DEALLOCATE (almo_scf_env%matrix_x)
2745 122 : DEALLOCATE (almo_scf_env%matrix_ov)
2746 122 : DEALLOCATE (almo_scf_env%matrix_ov_full)
2747 122 : DEALLOCATE (almo_scf_env%matrix_sigma_vv)
2748 122 : DEALLOCATE (almo_scf_env%matrix_sigma_vv_blk)
2749 122 : DEALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt)
2750 122 : DEALLOCATE (almo_scf_env%matrix_sigma_vv_sqrt_inv)
2751 122 : DEALLOCATE (almo_scf_env%matrix_vv_full_blk)
2752 122 : IF (almo_scf_env%deloc_truncate_virt /= virt_full) THEN
2753 0 : DEALLOCATE (almo_scf_env%matrix_k_tr)
2754 0 : DEALLOCATE (almo_scf_env%matrix_k_blk)
2755 0 : DEALLOCATE (almo_scf_env%matrix_v_disc)
2756 0 : DEALLOCATE (almo_scf_env%matrix_v_disc_blk)
2757 0 : DEALLOCATE (almo_scf_env%matrix_ov_disc)
2758 0 : DEALLOCATE (almo_scf_env%matrix_vv_disc_blk)
2759 0 : DEALLOCATE (almo_scf_env%matrix_vv_disc)
2760 0 : DEALLOCATE (almo_scf_env%matrix_k_blk_ones)
2761 0 : DEALLOCATE (almo_scf_env%opt_k_t_dd)
2762 0 : DEALLOCATE (almo_scf_env%opt_k_t_rr)
2763 0 : DEALLOCATE (almo_scf_env%opt_k_denom)
2764 : END IF
2765 : END IF
2766 122 : IF (almo_scf_env%need_previous_ks) THEN
2767 122 : DEALLOCATE (almo_scf_env%matrix_ks_0deloc)
2768 : END IF
2769 122 : IF (almo_scf_env%need_orbital_energies) THEN
2770 122 : DEALLOCATE (almo_scf_env%matrix_eoo)
2771 122 : DEALLOCATE (almo_scf_env%matrix_evv_full)
2772 : END IF
2773 :
2774 : ! clean up other variables
2775 250 : DO ispin = 1, almo_scf_env%nspins
2776 : CALL release_submatrices( &
2777 128 : almo_scf_env%domain_preconditioner(:, ispin))
2778 128 : CALL release_submatrices(almo_scf_env%domain_s_inv(:, ispin))
2779 128 : CALL release_submatrices(almo_scf_env%domain_s_sqrt_inv(:, ispin))
2780 128 : CALL release_submatrices(almo_scf_env%domain_s_sqrt(:, ispin))
2781 128 : CALL release_submatrices(almo_scf_env%domain_ks_xx(:, ispin))
2782 128 : CALL release_submatrices(almo_scf_env%domain_t(:, ispin))
2783 128 : CALL release_submatrices(almo_scf_env%domain_err(:, ispin))
2784 250 : CALL release_submatrices(almo_scf_env%domain_r_down_up(:, ispin))
2785 : END DO
2786 956 : DEALLOCATE (almo_scf_env%domain_preconditioner)
2787 956 : DEALLOCATE (almo_scf_env%domain_s_inv)
2788 956 : DEALLOCATE (almo_scf_env%domain_s_sqrt_inv)
2789 956 : DEALLOCATE (almo_scf_env%domain_s_sqrt)
2790 956 : DEALLOCATE (almo_scf_env%domain_ks_xx)
2791 956 : DEALLOCATE (almo_scf_env%domain_t)
2792 956 : DEALLOCATE (almo_scf_env%domain_err)
2793 956 : DEALLOCATE (almo_scf_env%domain_r_down_up)
2794 250 : DO ispin = 1, almo_scf_env%nspins
2795 128 : DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs)
2796 250 : DEALLOCATE (almo_scf_env%domain_map(ispin)%index1)
2797 : END DO
2798 250 : DEALLOCATE (almo_scf_env%domain_map)
2799 122 : DEALLOCATE (almo_scf_env%domain_index_of_ao)
2800 122 : DEALLOCATE (almo_scf_env%domain_index_of_atom)
2801 122 : DEALLOCATE (almo_scf_env%first_atom_of_domain)
2802 122 : DEALLOCATE (almo_scf_env%last_atom_of_domain)
2803 122 : DEALLOCATE (almo_scf_env%nbasis_of_domain)
2804 122 : IF (ALLOCATED(almo_scf_env%nocc_of_domain)) THEN
2805 122 : DEALLOCATE (almo_scf_env%nocc_of_domain)
2806 : END IF
2807 122 : DEALLOCATE (almo_scf_env%real_ne_of_domain)
2808 122 : DEALLOCATE (almo_scf_env%nvirt_full_of_domain)
2809 122 : DEALLOCATE (almo_scf_env%nvirt_of_domain)
2810 122 : DEALLOCATE (almo_scf_env%nvirt_disc_of_domain)
2811 122 : DEALLOCATE (almo_scf_env%mu_of_domain)
2812 122 : DEALLOCATE (almo_scf_env%cpu_of_domain)
2813 122 : DEALLOCATE (almo_scf_env%charge_of_domain)
2814 122 : DEALLOCATE (almo_scf_env%multiplicity_of_domain)
2815 122 : DEALLOCATE (almo_scf_env%activate)
2816 122 : IF (almo_scf_env%smear) THEN
2817 4 : DEALLOCATE (almo_scf_env%mo_energies)
2818 4 : DEALLOCATE (almo_scf_env%kTS)
2819 : END IF
2820 :
2821 122 : DEALLOCATE (almo_scf_env%domain_index_of_ao_block)
2822 122 : DEALLOCATE (almo_scf_env%domain_index_of_mo_block)
2823 :
2824 122 : CALL mp_para_env_release(almo_scf_env%para_env)
2825 122 : CALL cp_blacs_env_release(almo_scf_env%blacs_env)
2826 :
2827 122 : CALL timestop(handle)
2828 :
2829 122 : END SUBROUTINE almo_scf_clean_up
2830 :
2831 : ! **************************************************************************************************
2832 : !> \brief Do post scf calculations with ALMO
2833 : !> WARNING: ALMO post scf calculation may not work for certain quantities,
2834 : !> like forces, since ALMO wave function is only 'partially' optimized
2835 : !> \param qs_env ...
2836 : !> \par History
2837 : !> 2016.12 created [Yifei Shi]
2838 : !> \author Yifei Shi
2839 : ! **************************************************************************************************
2840 122 : SUBROUTINE almo_post_scf_compute_properties(qs_env)
2841 : TYPE(qs_environment_type), POINTER :: qs_env
2842 :
2843 122 : CALL qs_scf_compute_properties(qs_env)
2844 :
2845 122 : END SUBROUTINE almo_post_scf_compute_properties
2846 :
2847 : END MODULE almo_scf
2848 :
|
