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