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