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