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