Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief orbital transformations
10 : !> \par History
11 : !> Added Taylor expansion based computation of the matrix functions (01.2004)
12 : !> added additional rotation variables for non-equivalent occupied orbs (08.2004)
13 : !> \author Joost VandeVondele (06.2002)
14 : ! **************************************************************************************************
15 : MODULE qs_ot_types
16 : USE bibliography, ONLY: VandeVondele2003,&
17 : Weber2008,&
18 : cite_reference
19 : USE cp_blacs_env, ONLY: cp_blacs_env_release,&
20 : cp_blacs_env_type
21 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
22 : dbcsr_get_info,&
23 : dbcsr_init_p,&
24 : dbcsr_p_type,&
25 : dbcsr_release_p,&
26 : dbcsr_set,&
27 : dbcsr_type,&
28 : dbcsr_type_no_symmetry
29 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag
30 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_m_by_n_from_row_template,&
31 : cp_dbcsr_m_by_n_from_template,&
32 : dbcsr_allocate_matrix_set,&
33 : dbcsr_deallocate_matrix_set
34 : USE cp_fm_struct, ONLY: cp_fm_struct_get,&
35 : cp_fm_struct_type
36 : USE cp_fm_types, ONLY: cp_fm_release,&
37 : cp_fm_type
38 : USE input_constants, ONLY: &
39 : cholesky_reduce, ls_2pnt, ls_3pnt, ls_adapt, ls_gold, ls_none, ot_algo_irac, &
40 : ot_algo_taylor_or_diag, ot_chol_irac, ot_lwdn_irac, ot_mini_broyden, ot_mini_cg, &
41 : ot_mini_diis, ot_mini_lbfgs, ot_mini_sd, ot_poly_irac, ot_precond_full_all, &
42 : ot_precond_full_kinetic, ot_precond_full_single, ot_precond_full_single_inverse, &
43 : ot_precond_none, ot_precond_s_inverse, ot_precond_solver_chebyshev, &
44 : ot_precond_solver_default, ot_precond_solver_direct, ot_precond_solver_inv_chol, &
45 : ot_precond_solver_update
46 : USE input_section_types, ONLY: section_vals_type,&
47 : section_vals_val_get
48 : USE kinds, ONLY: dp
49 : USE message_passing, ONLY: mp_para_env_release,&
50 : mp_para_env_type
51 : USE preconditioner_types, ONLY: preconditioner_type
52 : #include "./base/base_uses.f90"
53 :
54 : IMPLICIT NONE
55 :
56 : PRIVATE
57 :
58 : PUBLIC :: qs_ot_type
59 : PUBLIC :: qs_ot_settings_type
60 : PUBLIC :: qs_ot_destroy
61 : PUBLIC :: qs_ot_allocate
62 : PUBLIC :: qs_ot_allocate_complex_state
63 : PUBLIC :: qs_ot_channel_index
64 : PUBLIC :: qs_ot_check_channel_context
65 : PUBLIC :: qs_ot_init
66 : PUBLIC :: qs_ot_kpoint_preconditioner_solver_supported
67 : PUBLIC :: qs_ot_kpoint_preconditioner_scale
68 : PUBLIC :: qs_ot_kpoint_preconditioner_supported
69 : PUBLIC :: qs_ot_number_of_channels
70 : PUBLIC :: qs_ot_settings_init
71 : PUBLIC :: qs_ot_set_context
72 : PUBLIC :: ot_readwrite_input
73 :
74 : ! **************************************************************************************************
75 : !> \brief notice, this variable needs to be copyable, needed for spins as e.g. in qs_ot_scf
76 : ! **************************************************************************************************
77 : TYPE qs_ot_settings_type
78 : LOGICAL :: do_rotation = .FALSE., do_ener = .FALSE.
79 : LOGICAL :: ks = .FALSE.
80 : CHARACTER(LEN=4) :: ot_method = ""
81 : CHARACTER(LEN=3) :: ot_algorithm = ""
82 : CHARACTER(LEN=4) :: line_search_method = ""
83 : CHARACTER(LEN=20) :: preconditioner_name = ""
84 : INTEGER :: preconditioner_type = -1
85 : INTEGER :: cholesky_type = -1
86 : INTEGER :: ot_state = 0
87 : CHARACTER(LEN=20) :: precond_solver_name = ""
88 : INTEGER :: precond_solver_type = -1
89 : INTEGER :: chebyshev_degree = 8
90 : LOGICAL :: safer_diis = .FALSE.
91 : REAL(KIND=dp) :: ds_min = -1.0_dp
92 : REAL(KIND=dp) :: energy_gap = -1.0_dp
93 : INTEGER :: diis_m = -1
94 : REAL(KIND=dp) :: lbfgs_curvature_tol = 1.0E-4_dp
95 : LOGICAL :: lbfgs_damping = .TRUE.
96 : INTEGER :: max_scf_diis = 0
97 : REAL(KIND=dp) :: gold_target = -1.0_dp
98 : REAL(KIND=dp) :: eps_taylor = -1.0_dp ! minimum accuracy of Taylor expansion
99 : INTEGER :: max_taylor = -1 ! maximum order of Taylor expansion before switching to diagonalization
100 : INTEGER :: irac_degree = -1 ! used to control the refinement polynomial degree
101 : INTEGER :: max_irac = -1 ! maximum number of iteration for refinement
102 : REAL(KIND=dp) :: eps_irac = -1.0_dp ! target accuracy for refinement
103 : REAL(KIND=dp) :: eps_irac_quick_exit = -1.0_dp
104 : REAL(KIND=dp) :: eps_irac_filter_matrix = -1.0_dp
105 : REAL(KIND=dp) :: eps_irac_switch = -1.0_dp
106 : LOGICAL :: on_the_fly_loc = .FALSE.
107 : CHARACTER(LEN=4) :: ortho_irac = ""
108 : LOGICAL :: occupation_preconditioner = .FALSE., add_nondiag_energy = .FALSE.
109 : REAL(KIND=dp) :: nondiag_energy_strength = -1.0_dp
110 : REAL(KIND=dp) :: broyden_beta = -1.0_dp, broyden_gamma = -1.0_dp, broyden_sigma = -1.0_dp
111 : REAL(KIND=dp) :: broyden_eta = -1.0_dp, broyden_omega = -1.0_dp, broyden_sigma_decrease = -1.0_dp
112 : REAL(KIND=dp) :: broyden_sigma_min = -1.0_dp
113 : LOGICAL :: broyden_forget_history = .FALSE., broyden_adaptive_sigma = .FALSE.
114 : LOGICAL :: broyden_enable_flip = .FALSE.
115 : END TYPE qs_ot_settings_type
116 :
117 : ! **************************************************************************************************
118 : TYPE qs_ot_type
119 : ! this sets the method to be used
120 : TYPE(qs_ot_settings_type) :: settings = qs_ot_settings_type()
121 : LOGICAL :: restricted = .FALSE.
122 : INTEGER :: spin_index = 0, kpoint_index = 0, local_kpoint_index = 0
123 : REAL(KIND=dp) :: kpoint_weight = 1.0_dp
124 : LOGICAL :: has_kpoint_context = .FALSE.
125 : LOGICAL :: state_allocated = .FALSE.
126 : LOGICAL :: has_complex_kpoint_state = .FALSE.
127 :
128 : ! first part of the variables, for occupied subspace invariant optimisation
129 :
130 : ! add a preconditioner matrix. should be symmetric and positive definite
131 : ! the type of this matrix might change in the future
132 : TYPE(preconditioner_type), POINTER :: preconditioner => NULL()
133 :
134 : ! these will/might change during iterations
135 :
136 : ! OT / TOD
137 : TYPE(dbcsr_type), POINTER :: matrix_p => NULL(), matrix_p_im => NULL()
138 : TYPE(dbcsr_type), POINTER :: matrix_r => NULL(), matrix_r_im => NULL()
139 : TYPE(dbcsr_type), POINTER :: matrix_sinp => NULL(), matrix_sinp_im => NULL()
140 : TYPE(dbcsr_type), POINTER :: matrix_cosp => NULL(), matrix_cosp_im => NULL()
141 : TYPE(dbcsr_type), POINTER :: matrix_sinp_b => NULL()
142 : TYPE(dbcsr_type), POINTER :: matrix_cosp_b => NULL()
143 : TYPE(dbcsr_type), POINTER :: matrix_buf1 => NULL(), matrix_buf1_im => NULL()
144 : TYPE(dbcsr_type), POINTER :: matrix_buf2 => NULL(), matrix_buf2_im => NULL()
145 : TYPE(dbcsr_type), POINTER :: matrix_buf3 => NULL(), matrix_buf3_im => NULL()
146 : TYPE(dbcsr_type), POINTER :: matrix_buf4 => NULL(), matrix_buf4_im => NULL()
147 : TYPE(dbcsr_type), POINTER :: matrix_os => NULL(), matrix_os_im => NULL()
148 : TYPE(dbcsr_type), POINTER :: matrix_buf1_ortho => NULL(), matrix_buf1_ortho_im => NULL()
149 : TYPE(dbcsr_type), POINTER :: matrix_buf2_ortho => NULL(), matrix_buf2_ortho_im => NULL()
150 : TYPE(dbcsr_type), POINTER :: matrix_tmp_ortho => NULL()
151 : TYPE(dbcsr_type), POINTER :: matrix_buf_nk => NULL(), matrix_buf_nk_im => NULL(), &
152 : matrix_tmp_nk => NULL()
153 :
154 : REAL(KIND=dp), DIMENSION(:), POINTER :: evals => NULL()
155 : REAL(KIND=dp), DIMENSION(:), POINTER :: dum => NULL()
156 :
157 : ! matrix os valid
158 : LOGICAL :: os_valid = .FALSE.
159 :
160 : ! for efficient/parallel writing to the blacs_matrix
161 : TYPE(mp_para_env_type), POINTER :: para_env => NULL()
162 : TYPE(cp_blacs_env_type), POINTER :: blacs_env => NULL()
163 :
164 : ! mo-like vectors
165 : TYPE(dbcsr_type), POINTER :: matrix_c0 => NULL(), matrix_sc0 => NULL(), matrix_psc0 => NULL()
166 : TYPE(dbcsr_type), POINTER :: matrix_c0_im => NULL(), matrix_sc0_im => NULL(), matrix_psc0_im => NULL()
167 :
168 : ! OT / IR
169 : TYPE(dbcsr_type), POINTER :: buf1_k_k_nosym => NULL(), buf2_k_k_nosym => NULL(), &
170 : buf3_k_k_nosym => NULL(), buf4_k_k_nosym => NULL(), &
171 : buf1_k_k_sym => NULL(), buf2_k_k_sym => NULL(), &
172 : buf3_k_k_sym => NULL(), buf4_k_k_sym => NULL(), &
173 : p_k_k_sym => NULL(), buf1_n_k => NULL(), buf1_n_k_dp => NULL()
174 :
175 : ! only here for the ease of programming. These will have to be supplied
176 : ! explicitly at all times
177 : TYPE(dbcsr_type), POINTER :: matrix_x => NULL(), matrix_sx => NULL(), matrix_gx => NULL()
178 : TYPE(dbcsr_type), POINTER :: matrix_x_im => NULL(), matrix_sx_im => NULL(), matrix_gx_im => NULL()
179 : TYPE(dbcsr_type), POINTER :: matrix_preconditioned_gx => NULL(), &
180 : matrix_preconditioned_gx_im => NULL()
181 : TYPE(dbcsr_type), POINTER :: matrix_response_gx => NULL(), matrix_response_gx_im => NULL()
182 : TYPE(dbcsr_type), POINTER :: matrix_mermin_g0 => NULL(), matrix_mermin_g0_im => NULL()
183 : TYPE(dbcsr_type), POINTER :: matrix_ref_inv_sqrt => NULL(), matrix_ref_inv_sqrt_im => NULL()
184 : ! Owned physical endpoints for the bounded, gauge-invariant Hxc response secant.
185 : TYPE(cp_fm_type), POINTER :: fm_mermin_c0 => NULL(), fm_mermin_c0_im => NULL()
186 : TYPE(cp_fm_type), POINTER :: fm_mermin_h0 => NULL(), fm_mermin_h0_im => NULL()
187 : TYPE(cp_fm_type), POINTER :: fm_mermin_c_previous => NULL(), &
188 : fm_mermin_c_previous_im => NULL()
189 : TYPE(cp_fm_type), POINTER :: fm_mermin_y_previous => NULL(), &
190 : fm_mermin_y_previous_im => NULL()
191 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mermin_occupation0, &
192 : mermin_occupation_previous, &
193 : ener_mermin_g0
194 : LOGICAL :: mermin_physical_ref_valid = .FALSE., mermin_physical_secant_valid = .FALSE.
195 : LOGICAL :: mermin_gradient_ref_valid = .FALSE.
196 : TYPE(dbcsr_type), POINTER :: matrix_dx => NULL(), matrix_gx_old => NULL()
197 : TYPE(dbcsr_type), POINTER :: matrix_dx_im => NULL(), matrix_gx_old_im => NULL()
198 :
199 : LOGICAL :: use_gx_old = .FALSE., use_dx = .FALSE.
200 :
201 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h_e => NULL(), matrix_h_x => NULL()
202 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h_e_im => NULL(), matrix_h_x_im => NULL()
203 : REAL(KIND=dp), DIMENSION(:), POINTER :: lbfgs_rho => NULL(), lbfgs_yy => NULL()
204 : REAL(KIND=dp), DIMENSION(:), POINTER :: lbfgs_sy_rotation => NULL(), &
205 : lbfgs_yy_rotation => NULL()
206 :
207 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ls_diis => NULL()
208 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: lss_diis => NULL()
209 : REAL(KIND=dp), DIMENSION(:), POINTER :: c_diis => NULL()
210 : REAL(KIND=dp), DIMENSION(:), POINTER :: c_broy => NULL()
211 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy_h => NULL()
212 : INTEGER, DIMENSION(:), POINTER :: ipivot => NULL()
213 :
214 : REAL(KIND=dp) :: ot_pos(53) = -1.0_dp, ot_energy(53) = -1.0_dp, ot_grad(53) = -1.0_dp ! HARD LIMIT FOR THE LS
215 : INTEGER :: line_search_left = -1, line_search_right = -1, line_search_mid = -1
216 : INTEGER :: line_search_count = -1
217 : LOGICAL :: line_search_might_be_done = .FALSE.
218 : REAL(KIND=dp) :: delta = -1.0_dp, gnorm = -1.0_dp, gnorm_old = -1.0_dp, etotal = -1.0_dp, gradient = -1.0_dp
219 : LOGICAL :: energy_only = .FALSE.
220 : INTEGER :: diis_iter = -1
221 : CHARACTER(LEN=8) :: OT_METHOD_FULL = ""
222 : INTEGER :: OT_count = -1
223 : REAL(KIND=dp) :: ds_min = -1.0_dp
224 : REAL(KIND=dp) :: broyden_adaptive_sigma = -1.0_dp
225 :
226 : LOGICAL :: do_taylor = .FALSE.
227 : INTEGER :: taylor_order = -1
228 : REAL(KIND=dp) :: largest_eval_upper_bound = -1.0_dp
229 :
230 : ! second part of the variables, if an explicit rotation is required as well
231 : TYPE(dbcsr_type), POINTER :: rot_mat_u => NULL() ! rotation matrix
232 : TYPE(dbcsr_type), POINTER :: rot_mat_u_im => NULL() ! imaginary part at complex k points
233 : TYPE(dbcsr_type), POINTER :: rot_mat_x => NULL() ! antisymmetric matrix that parametrises rot_matrix_u
234 : TYPE(dbcsr_type), POINTER :: rot_mat_x_im => NULL() ! symmetric imaginary generator component
235 : TYPE(dbcsr_type), POINTER :: rot_mat_dedu => NULL() ! derivative of the total energy wrt to u
236 : TYPE(dbcsr_type), POINTER :: rot_mat_dedu_im => NULL()
237 : TYPE(dbcsr_type), POINTER :: rot_mat_chc => NULL() ! for convencience, the matrix c^T H c
238 : TYPE(dbcsr_type), POINTER :: rot_mat_chc_im => NULL()
239 :
240 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rot_mat_h_e => NULL()
241 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rot_mat_h_x => NULL()
242 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rot_mat_h_e_im => NULL()
243 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rot_mat_h_x_im => NULL()
244 : TYPE(dbcsr_type), POINTER :: rot_mat_gx => NULL()
245 : TYPE(dbcsr_type), POINTER :: rot_mat_gx_im => NULL()
246 : TYPE(dbcsr_type), POINTER :: rot_mat_response_gx => NULL()
247 : TYPE(dbcsr_type), POINTER :: rot_mat_response_gx_im => NULL()
248 : TYPE(dbcsr_type), POINTER :: rot_mat_mermin_g0 => NULL()
249 : TYPE(dbcsr_type), POINTER :: rot_mat_mermin_g0_im => NULL()
250 : TYPE(dbcsr_type), POINTER :: rot_mat_gx_old => NULL()
251 : TYPE(dbcsr_type), POINTER :: rot_mat_gx_old_im => NULL()
252 : TYPE(dbcsr_type), POINTER :: rot_mat_dx => NULL()
253 : TYPE(dbcsr_type), POINTER :: rot_mat_dx_im => NULL()
254 :
255 : REAL(KIND=dp), DIMENSION(:), POINTER :: rot_mat_evals => NULL()
256 : TYPE(dbcsr_type), POINTER :: rot_mat_evec_re => NULL()
257 : TYPE(dbcsr_type), POINTER :: rot_mat_evec_im => NULL()
258 : LOGICAL :: rotation_response_valid = .FALSE.
259 : LOGICAL :: response_candidate_pending = .FALSE.
260 : LOGICAL :: response_shadow_pending = .FALSE.
261 : INTEGER :: response_candidate_directions = 0
262 : INTEGER :: response_candidate_good_samples = 0
263 : INTEGER :: response_candidate_cooldown = 0
264 : INTEGER :: response_shadow_good_samples = 0
265 : REAL(KIND=dp) :: response_model_curvature = 0.0_dp
266 : REAL(KIND=dp) :: response_shadow_curvature = 0.0_dp
267 : LOGICAL :: response_hxc_direction_valid = .FALSE.
268 : REAL(KIND=dp) :: response_reference_energy = 0.0_dp
269 : REAL(KIND=dp) :: response_reference_residual = 0.0_dp
270 : REAL(KIND=dp) :: response_predicted_slope = 0.0_dp
271 : REAL(KIND=dp) :: response_predicted_curvature = 0.0_dp
272 :
273 : ! third part of the variables, if we need to optimize orbital energies
274 : REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_x => NULL()
275 : REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_rayleigh => NULL()
276 : REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_dx => NULL()
277 : REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_gx => NULL()
278 : REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_preconditioned_gx => NULL()
279 : REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_response_gx => NULL()
280 : REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_gx_old => NULL()
281 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: ener_h_e => NULL()
282 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: ener_h_x => NULL()
283 : END TYPE qs_ot_type
284 :
285 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_types'
286 :
287 : CONTAINS
288 :
289 : ! **************************************************************************************************
290 : !> \brief sets default values for the settings type
291 : !> \param settings ...
292 : !> \par History
293 : !> 10.2004 created [Joost VandeVondele]
294 : ! **************************************************************************************************
295 17155 : SUBROUTINE qs_ot_settings_init(settings)
296 : TYPE(qs_ot_settings_type) :: settings
297 :
298 17155 : settings%ot_method = "CG"
299 17155 : settings%ot_algorithm = "TOD"
300 17155 : settings%diis_m = 7
301 17155 : settings%lbfgs_curvature_tol = 1.0E-4_dp
302 17155 : settings%lbfgs_damping = .TRUE.
303 17155 : settings%preconditioner_name = "FULL_KINETIC"
304 17155 : settings%preconditioner_type = ot_precond_full_kinetic
305 17155 : settings%cholesky_type = cholesky_reduce
306 17155 : settings%precond_solver_name = "CHOLESKY_INVERSE"
307 17155 : settings%precond_solver_type = ot_precond_solver_inv_chol
308 17155 : settings%chebyshev_degree = 8
309 17155 : settings%line_search_method = "2PNT"
310 17155 : settings%ds_min = 0.15_dp
311 17155 : settings%safer_diis = .TRUE.
312 17155 : settings%energy_gap = 0.2_dp
313 17155 : settings%eps_taylor = 1.0E-16_dp
314 17155 : settings%max_taylor = 4
315 17155 : settings%gold_target = 0.01_dp
316 17155 : settings%do_rotation = .FALSE.
317 17155 : settings%do_ener = .FALSE.
318 17155 : settings%irac_degree = 4
319 17155 : settings%max_irac = 50
320 17155 : settings%max_scf_diis = 0
321 17155 : settings%eps_irac = 1.0E-10_dp
322 17155 : settings%eps_irac_quick_exit = 1.0E-5_dp
323 17155 : settings%eps_irac_switch = 1.0E-2
324 17155 : settings%eps_irac_filter_matrix = 0.0_dp
325 17155 : settings%on_the_fly_loc = .FALSE.
326 17155 : settings%ortho_irac = "CHOL"
327 17155 : settings%ks = .TRUE.
328 17155 : settings%occupation_preconditioner = .FALSE.
329 17155 : settings%add_nondiag_energy = .FALSE.
330 17155 : settings%nondiag_energy_strength = 0.0_dp
331 :
332 17155 : END SUBROUTINE qs_ot_settings_init
333 :
334 : ! **************************************************************************************************
335 : !> \brief label an OT environment by spin and optional irreducible k-point context
336 : !> \param qs_ot_env ...
337 : !> \param spin_index ...
338 : !> \param kpoint_index ...
339 : !> \param local_kpoint_index ...
340 : !> \param kpoint_weight ...
341 : ! **************************************************************************************************
342 17912 : SUBROUTINE qs_ot_set_context(qs_ot_env, spin_index, kpoint_index, local_kpoint_index, kpoint_weight)
343 : TYPE(qs_ot_type) :: qs_ot_env
344 : INTEGER, INTENT(IN), OPTIONAL :: spin_index, kpoint_index, &
345 : local_kpoint_index
346 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kpoint_weight
347 :
348 17912 : IF (PRESENT(spin_index)) qs_ot_env%spin_index = spin_index
349 17912 : IF (PRESENT(kpoint_index)) qs_ot_env%kpoint_index = kpoint_index
350 17912 : IF (PRESENT(local_kpoint_index)) qs_ot_env%local_kpoint_index = local_kpoint_index
351 17912 : IF (PRESENT(kpoint_weight)) qs_ot_env%kpoint_weight = kpoint_weight
352 17912 : qs_ot_env%has_kpoint_context = PRESENT(kpoint_index) .OR. PRESENT(local_kpoint_index)
353 :
354 17912 : END SUBROUTINE qs_ot_set_context
355 :
356 : ! **************************************************************************************************
357 : !> \brief number of OT optimization channels for spin and k-point resolved state
358 : !> \param nspin ...
359 : !> \param nkpoint ...
360 : !> \param restricted ...
361 : !> \return ...
362 : ! **************************************************************************************************
363 17832 : INTEGER FUNCTION qs_ot_number_of_channels(nspin, nkpoint, restricted)
364 : INTEGER, INTENT(IN) :: nspin
365 : INTEGER, INTENT(IN), OPTIONAL :: nkpoint
366 : LOGICAL, INTENT(IN), OPTIONAL :: restricted
367 :
368 : INTEGER :: nkpoint_eff, nspin_eff
369 :
370 17832 : nspin_eff = nspin
371 17832 : IF (PRESENT(restricted)) THEN
372 17832 : IF (restricted) nspin_eff = 1
373 : END IF
374 17832 : nkpoint_eff = 1
375 17832 : IF (PRESENT(nkpoint)) nkpoint_eff = nkpoint
376 :
377 17832 : CPASSERT(nspin_eff >= 1)
378 17832 : CPASSERT(nkpoint_eff >= 1)
379 17832 : qs_ot_number_of_channels = nspin_eff*nkpoint_eff
380 :
381 17832 : END FUNCTION qs_ot_number_of_channels
382 :
383 : ! **************************************************************************************************
384 : !> \brief Return whether a preconditioner is implemented for complex k-point OT.
385 : !> \param preconditioner_type selected OT preconditioner
386 : !> \param use_real_wfn whether the k-point channel stores a real wavefunction
387 : !> \return ...
388 : ! **************************************************************************************************
389 322 : PURE ELEMENTAL FUNCTION qs_ot_kpoint_preconditioner_supported(preconditioner_type, use_real_wfn) &
390 : RESULT(supported)
391 : INTEGER, INTENT(IN) :: preconditioner_type
392 : LOGICAL, INTENT(IN) :: use_real_wfn
393 : LOGICAL :: supported
394 :
395 322 : supported = .FALSE.
396 322 : IF (.NOT. use_real_wfn) THEN
397 : supported = preconditioner_type == ot_precond_none .OR. &
398 : preconditioner_type == ot_precond_full_all .OR. &
399 : preconditioner_type == ot_precond_full_single .OR. &
400 : preconditioner_type == ot_precond_full_single_inverse .OR. &
401 : preconditioner_type == ot_precond_full_kinetic .OR. &
402 320 : preconditioner_type == ot_precond_s_inverse
403 : END IF
404 :
405 322 : END FUNCTION qs_ot_kpoint_preconditioner_supported
406 :
407 : ! **************************************************************************************************
408 : !> \brief Return whether a solver is implemented for a complex k-point OT preconditioner.
409 : !> \param preconditioner_type selected OT preconditioner
410 : !> \param solver_type selected preconditioner solver
411 : !> \param use_real_wfn whether the k-point channel stores a real wavefunction
412 : !> \return ...
413 : ! **************************************************************************************************
414 168 : PURE ELEMENTAL FUNCTION qs_ot_kpoint_preconditioner_solver_supported( &
415 : preconditioner_type, solver_type, use_real_wfn) RESULT(supported)
416 : INTEGER, INTENT(IN) :: preconditioner_type, solver_type
417 : LOGICAL, INTENT(IN) :: use_real_wfn
418 : LOGICAL :: supported
419 :
420 168 : supported = qs_ot_kpoint_preconditioner_supported(preconditioner_type, use_real_wfn)
421 168 : IF (.NOT. supported) RETURN
422 168 : IF (preconditioner_type == ot_precond_full_all .OR. &
423 : preconditioner_type == ot_precond_full_single) THEN
424 74 : supported = solver_type == ot_precond_solver_default
425 : ELSE IF (preconditioner_type == ot_precond_full_single_inverse .OR. &
426 94 : preconditioner_type == ot_precond_full_kinetic .OR. &
427 : preconditioner_type == ot_precond_s_inverse) THEN
428 : supported = solver_type == ot_precond_solver_default .OR. &
429 86 : solver_type == ot_precond_solver_inv_chol
430 : END IF
431 :
432 : END FUNCTION qs_ot_kpoint_preconditioner_solver_supported
433 :
434 : ! **************************************************************************************************
435 : !> \brief Scale an inverse k-point Hessian block consistently with its irreducible weight.
436 : !> \param kpoint_weight positive irreducible k-point weight
437 : !> \return ...
438 : ! **************************************************************************************************
439 4467 : PURE ELEMENTAL REAL(KIND=dp) FUNCTION qs_ot_kpoint_preconditioner_scale(kpoint_weight)
440 : REAL(KIND=dp), INTENT(IN) :: kpoint_weight
441 :
442 4467 : qs_ot_kpoint_preconditioner_scale = 1.0_dp/kpoint_weight
443 :
444 4467 : END FUNCTION qs_ot_kpoint_preconditioner_scale
445 :
446 : ! **************************************************************************************************
447 : !> \brief flat OT channel index for a spin/k-point pair
448 : !> \param ispin ...
449 : !> \param ikpoint ...
450 : !> \param nspin ...
451 : !> \return ...
452 : ! **************************************************************************************************
453 51986 : INTEGER FUNCTION qs_ot_channel_index(ispin, ikpoint, nspin)
454 : INTEGER, INTENT(IN) :: ispin, ikpoint, nspin
455 :
456 51986 : CPASSERT(ispin >= 1)
457 51986 : CPASSERT(ikpoint >= 1)
458 51986 : CPASSERT(nspin >= 1)
459 51986 : CPASSERT(ispin <= nspin)
460 51986 : qs_ot_channel_index = (ikpoint - 1)*nspin + ispin
461 :
462 51986 : END FUNCTION qs_ot_channel_index
463 :
464 : ! **************************************************************************************************
465 : !> \brief validate the flat OT channel identity for spin and optional irreducible k-points
466 : !> \param qs_ot_env ...
467 : !> \param nspin ...
468 : !> \param nkpoint ...
469 : !> \param restricted ...
470 : !> \param require_kpoint ...
471 : !> \param kp_range ...
472 : !> \param wkp ...
473 : !> \param require_local_state ...
474 : !> \param require_complex_state ...
475 : ! **************************************************************************************************
476 10281 : SUBROUTINE qs_ot_check_channel_context(qs_ot_env, nspin, nkpoint, restricted, require_kpoint, &
477 : kp_range, wkp, require_local_state, require_complex_state)
478 : TYPE(qs_ot_type), DIMENSION(:), INTENT(IN) :: qs_ot_env
479 : INTEGER, INTENT(IN) :: nspin
480 : INTEGER, INTENT(IN), OPTIONAL :: nkpoint
481 : LOGICAL, INTENT(IN), OPTIONAL :: restricted, require_kpoint
482 : INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: kp_range
483 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: wkp
484 : LOGICAL, INTENT(IN), OPTIONAL :: require_local_state, &
485 : require_complex_state
486 :
487 : INTEGER :: expected_local_kpoint, ikpoint, ispin, &
488 : nkpoint_eff, nspin_eff, ot_channel
489 : LOGICAL :: check_complex_state, &
490 : check_kpoint_context, &
491 : check_local_state, check_weights, &
492 : context_ok, restricted_eff
493 : REAL(KIND=dp) :: weight_tol
494 :
495 10281 : nkpoint_eff = 1
496 10281 : IF (PRESENT(nkpoint)) nkpoint_eff = nkpoint
497 10281 : nspin_eff = nspin
498 10281 : restricted_eff = .FALSE.
499 10281 : IF (PRESENT(restricted)) THEN
500 10281 : restricted_eff = restricted
501 10281 : IF (restricted_eff) nspin_eff = 1
502 : END IF
503 :
504 10281 : check_kpoint_context = (nkpoint_eff > 1)
505 10281 : IF (PRESENT(require_kpoint)) check_kpoint_context = require_kpoint
506 10281 : check_weights = .FALSE.
507 10281 : IF (PRESENT(wkp)) check_weights = ASSOCIATED(wkp)
508 10281 : check_local_state = .FALSE.
509 10281 : IF (PRESENT(require_local_state)) check_local_state = require_local_state
510 10281 : check_complex_state = .FALSE.
511 10281 : IF (PRESENT(require_complex_state)) check_complex_state = require_complex_state
512 :
513 : context_ok = (SIZE(qs_ot_env) == qs_ot_number_of_channels(nspin, &
514 : nkpoint=nkpoint_eff, &
515 10281 : restricted=restricted_eff))
516 10281 : CPASSERT(context_ok)
517 23556 : DO ikpoint = 1, nkpoint_eff
518 13275 : expected_local_kpoint = 0
519 13275 : IF (PRESENT(kp_range)) THEN
520 13275 : IF (ikpoint >= kp_range(1) .AND. ikpoint <= kp_range(2)) THEN
521 4464 : expected_local_kpoint = ikpoint - kp_range(1) + 1
522 : END IF
523 : END IF
524 38578 : DO ispin = 1, nspin_eff
525 15022 : ot_channel = qs_ot_channel_index(ispin, ikpoint, nspin_eff)
526 15022 : CPASSERT(qs_ot_env(ot_channel)%spin_index == ispin)
527 15022 : IF (check_kpoint_context) THEN
528 6216 : CPASSERT(qs_ot_env(ot_channel)%has_kpoint_context)
529 6216 : CPASSERT(qs_ot_env(ot_channel)%kpoint_index == ikpoint)
530 6216 : IF (PRESENT(kp_range)) THEN
531 6216 : context_ok = (qs_ot_env(ot_channel)%local_kpoint_index == expected_local_kpoint)
532 6216 : CPASSERT(context_ok)
533 : END IF
534 6216 : IF (check_weights) THEN
535 5916 : weight_tol = 1000.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, ABS(wkp(ikpoint)))
536 5916 : context_ok = (ABS(qs_ot_env(ot_channel)%kpoint_weight - wkp(ikpoint)) <= weight_tol)
537 5916 : CPASSERT(context_ok)
538 : END IF
539 : END IF
540 28297 : IF (check_local_state) THEN
541 5916 : IF (expected_local_kpoint > 0 .OR. .NOT. PRESENT(kp_range)) THEN
542 4556 : CPASSERT(qs_ot_env(ot_channel)%state_allocated)
543 4556 : IF (check_complex_state) THEN
544 4556 : CPASSERT(qs_ot_env(ot_channel)%has_complex_kpoint_state)
545 : END IF
546 : ELSE
547 1360 : CPASSERT(.NOT. qs_ot_env(ot_channel)%state_allocated)
548 : END IF
549 : END IF
550 : END DO
551 : END DO
552 :
553 10281 : END SUBROUTINE qs_ot_check_channel_context
554 :
555 : ! **************************************************************************************************
556 : !> \brief init matrices, needs c0 and sc0 so that c0*sc0=1
557 : !> \param qs_ot_env ...
558 : ! **************************************************************************************************
559 9827 : SUBROUTINE qs_ot_init(qs_ot_env)
560 : TYPE(qs_ot_type) :: qs_ot_env
561 :
562 530658 : qs_ot_env%OT_energy(:) = 0.0_dp
563 530658 : qs_ot_env%OT_pos(:) = 0.0_dp
564 530658 : qs_ot_env%OT_grad(:) = 0.0_dp
565 9827 : qs_ot_env%line_search_count = 0
566 :
567 9827 : qs_ot_env%energy_only = .FALSE.
568 9827 : qs_ot_env%gnorm_old = 1.0_dp
569 9827 : qs_ot_env%diis_iter = 0
570 9827 : qs_ot_env%ds_min = qs_ot_env%settings%ds_min
571 9827 : qs_ot_env%os_valid = .FALSE.
572 :
573 9827 : CALL dbcsr_set(qs_ot_env%matrix_gx, 0.0_dp)
574 9827 : IF (ASSOCIATED(qs_ot_env%matrix_gx_im)) THEN
575 347 : CALL dbcsr_set(qs_ot_env%matrix_gx_im, 0.0_dp)
576 : END IF
577 9827 : IF (ASSOCIATED(qs_ot_env%matrix_preconditioned_gx)) THEN
578 112 : CALL dbcsr_set(qs_ot_env%matrix_preconditioned_gx, 0.0_dp)
579 : END IF
580 9827 : IF (ASSOCIATED(qs_ot_env%matrix_preconditioned_gx_im)) THEN
581 112 : CALL dbcsr_set(qs_ot_env%matrix_preconditioned_gx_im, 0.0_dp)
582 : END IF
583 9827 : IF (ASSOCIATED(qs_ot_env%matrix_response_gx)) CALL dbcsr_set(qs_ot_env%matrix_response_gx, 0.0_dp)
584 9827 : IF (ASSOCIATED(qs_ot_env%matrix_response_gx_im)) THEN
585 64 : CALL dbcsr_set(qs_ot_env%matrix_response_gx_im, 0.0_dp)
586 : END IF
587 9827 : IF (ASSOCIATED(qs_ot_env%matrix_mermin_g0)) CALL dbcsr_set(qs_ot_env%matrix_mermin_g0, 0.0_dp)
588 9827 : IF (ASSOCIATED(qs_ot_env%matrix_mermin_g0_im)) THEN
589 64 : CALL dbcsr_set(qs_ot_env%matrix_mermin_g0_im, 0.0_dp)
590 : END IF
591 9827 : qs_ot_env%mermin_gradient_ref_valid = .FALSE.
592 :
593 9827 : IF (qs_ot_env%use_dx) THEN
594 4213 : CALL dbcsr_set(qs_ot_env%matrix_dx, 0.0_dp)
595 : END IF
596 9827 : IF (qs_ot_env%use_dx .AND. ASSOCIATED(qs_ot_env%matrix_dx_im)) THEN
597 212 : CALL dbcsr_set(qs_ot_env%matrix_dx_im, 0.0_dp)
598 : END IF
599 :
600 9827 : IF (qs_ot_env%use_gx_old) THEN
601 4213 : CALL dbcsr_set(qs_ot_env%matrix_gx_old, 0.0_dp)
602 : END IF
603 9827 : IF (qs_ot_env%use_gx_old .AND. ASSOCIATED(qs_ot_env%matrix_gx_old_im)) THEN
604 212 : CALL dbcsr_set(qs_ot_env%matrix_gx_old_im, 0.0_dp)
605 : END IF
606 :
607 9827 : IF (qs_ot_env%settings%ot_method == "LBFG") THEN
608 772 : qs_ot_env%lbfgs_rho = 0.0_dp
609 772 : qs_ot_env%lbfgs_yy = 0.0_dp
610 772 : qs_ot_env%lbfgs_sy_rotation = 0.0_dp
611 772 : qs_ot_env%lbfgs_yy_rotation = 0.0_dp
612 : END IF
613 :
614 9827 : IF (qs_ot_env%settings%do_rotation) THEN
615 386 : CALL dbcsr_set(qs_ot_env%rot_mat_u, 0.0_dp)
616 386 : CALL dbcsr_add_on_diag(qs_ot_env%rot_mat_u, 1.0_dp)
617 386 : CALL dbcsr_set(qs_ot_env%rot_mat_x, 0.0_dp)
618 386 : CALL dbcsr_set(qs_ot_env%rot_mat_dedu, 0.0_dp)
619 386 : CALL dbcsr_set(qs_ot_env%rot_mat_chc, 0.0_dp)
620 386 : CALL dbcsr_set(qs_ot_env%rot_mat_gx, 0.0_dp)
621 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_response_gx)) THEN
622 112 : CALL dbcsr_set(qs_ot_env%rot_mat_response_gx, 0.0_dp)
623 : END IF
624 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_u_im)) CALL dbcsr_set(qs_ot_env%rot_mat_u_im, 0.0_dp)
625 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_x_im)) CALL dbcsr_set(qs_ot_env%rot_mat_x_im, 0.0_dp)
626 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_dedu_im)) CALL dbcsr_set(qs_ot_env%rot_mat_dedu_im, 0.0_dp)
627 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_chc_im)) CALL dbcsr_set(qs_ot_env%rot_mat_chc_im, 0.0_dp)
628 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_gx_im)) CALL dbcsr_set(qs_ot_env%rot_mat_gx_im, 0.0_dp)
629 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_response_gx_im)) THEN
630 112 : CALL dbcsr_set(qs_ot_env%rot_mat_response_gx_im, 0.0_dp)
631 : END IF
632 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_mermin_g0)) THEN
633 64 : CALL dbcsr_set(qs_ot_env%rot_mat_mermin_g0, 0.0_dp)
634 : END IF
635 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_mermin_g0_im)) THEN
636 64 : CALL dbcsr_set(qs_ot_env%rot_mat_mermin_g0_im, 0.0_dp)
637 : END IF
638 386 : qs_ot_env%rotation_response_valid = .FALSE.
639 386 : qs_ot_env%response_candidate_pending = .FALSE.
640 386 : qs_ot_env%response_shadow_pending = .FALSE.
641 386 : qs_ot_env%response_candidate_directions = 0
642 386 : qs_ot_env%response_candidate_good_samples = 0
643 386 : qs_ot_env%response_candidate_cooldown = 0
644 386 : qs_ot_env%response_shadow_good_samples = 0
645 386 : qs_ot_env%response_model_curvature = 0.0_dp
646 386 : qs_ot_env%response_shadow_curvature = 0.0_dp
647 386 : qs_ot_env%response_hxc_direction_valid = .FALSE.
648 386 : qs_ot_env%response_reference_energy = 0.0_dp
649 386 : qs_ot_env%response_reference_residual = 0.0_dp
650 386 : qs_ot_env%response_predicted_slope = 0.0_dp
651 386 : qs_ot_env%response_predicted_curvature = 0.0_dp
652 386 : IF (qs_ot_env%use_dx) THEN
653 300 : CALL dbcsr_set(qs_ot_env%rot_mat_dx, 0.0_dp)
654 300 : IF (ASSOCIATED(qs_ot_env%rot_mat_dx_im)) CALL dbcsr_set(qs_ot_env%rot_mat_dx_im, 0.0_dp)
655 : END IF
656 386 : IF (qs_ot_env%use_gx_old) THEN
657 300 : CALL dbcsr_set(qs_ot_env%rot_mat_gx_old, 0.0_dp)
658 300 : IF (ASSOCIATED(qs_ot_env%rot_mat_gx_old_im)) CALL dbcsr_set(qs_ot_env%rot_mat_gx_old_im, 0.0_dp)
659 : END IF
660 : END IF
661 9827 : IF (qs_ot_env%settings%do_ener) THEN
662 886 : qs_ot_env%ener_rayleigh(:) = 0.0_dp
663 886 : qs_ot_env%ener_gx(:) = 0.0_dp
664 114 : IF (ASSOCIATED(qs_ot_env%ener_preconditioned_gx)) THEN
665 840 : qs_ot_env%ener_preconditioned_gx(:) = 0.0_dp
666 : END IF
667 842 : IF (ASSOCIATED(qs_ot_env%ener_response_gx)) qs_ot_env%ener_response_gx(:) = 0.0_dp
668 506 : IF (ALLOCATED(qs_ot_env%ener_mermin_g0)) qs_ot_env%ener_mermin_g0(:) = 0.0_dp
669 114 : IF (qs_ot_env%use_dx) THEN
670 502 : qs_ot_env%ener_dx(:) = 0.0_dp
671 : END IF
672 114 : IF (qs_ot_env%use_gx_old) THEN
673 502 : qs_ot_env%ener_gx_old(:) = 0.0_dp
674 : END IF
675 : END IF
676 :
677 9827 : END SUBROUTINE qs_ot_init
678 :
679 : ! **************************************************************************************************
680 : !> \brief allocates the data in qs_ot_env, for a calculation with fm_struct_ref
681 : !> ortho_k allows for specifying an additional orthogonal subspace (i.e. c will
682 : !> be kept orthogonal provided c0 was, used in qs_ot_eigensolver)
683 : !> \param qs_ot_env ...
684 : !> \param matrix_s ...
685 : !> \param fm_struct_ref ...
686 : !> \param ortho_k ...
687 : !> \param energy_dimension number of auxiliary energy variables, defaults to the orbital count
688 : ! **************************************************************************************************
689 9827 : SUBROUTINE qs_ot_allocate(qs_ot_env, matrix_s, fm_struct_ref, ortho_k, energy_dimension)
690 : TYPE(qs_ot_type) :: qs_ot_env
691 : TYPE(dbcsr_type), POINTER :: matrix_s
692 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_ref
693 : INTEGER, OPTIONAL :: ortho_k, energy_dimension
694 :
695 : INTEGER :: i, k, m_diis, my_energy_dimension, &
696 : my_ortho_k, n, ncoef, nhistory
697 : TYPE(cp_blacs_env_type), POINTER :: context
698 : TYPE(mp_para_env_type), POINTER :: para_env
699 :
700 9827 : CALL cite_reference(VandeVondele2003)
701 :
702 9827 : CPASSERT(.NOT. qs_ot_env%state_allocated)
703 9827 : qs_ot_env%has_complex_kpoint_state = .FALSE.
704 9827 : NULLIFY (qs_ot_env%preconditioner)
705 9827 : NULLIFY (qs_ot_env%matrix_psc0)
706 9827 : NULLIFY (qs_ot_env%matrix_psc0_im)
707 9827 : NULLIFY (qs_ot_env%para_env)
708 9827 : NULLIFY (qs_ot_env%blacs_env)
709 :
710 : CALL cp_fm_struct_get(fm_struct_ref, nrow_global=n, ncol_global=k, &
711 9827 : para_env=para_env, context=context)
712 :
713 9827 : qs_ot_env%para_env => para_env
714 9827 : qs_ot_env%blacs_env => context
715 9827 : CALL para_env%retain()
716 9827 : CALL context%retain()
717 :
718 9827 : IF (PRESENT(ortho_k)) THEN
719 674 : my_ortho_k = ortho_k
720 : ELSE
721 9153 : my_ortho_k = k
722 : END IF
723 9827 : my_energy_dimension = k
724 9827 : IF (PRESENT(energy_dimension)) my_energy_dimension = energy_dimension
725 9827 : IF (qs_ot_env%settings%do_ener) THEN
726 114 : CPASSERT(my_energy_dimension > 0)
727 : END IF
728 :
729 9827 : m_diis = qs_ot_env%settings%diis_m
730 :
731 9827 : qs_ot_env%use_gx_old = .FALSE.
732 9827 : qs_ot_env%use_dx = .FALSE.
733 :
734 4213 : SELECT CASE (qs_ot_env%settings%ot_method)
735 : CASE ("SD")
736 : ! nothing
737 : CASE ("CG", "LBFG")
738 4213 : qs_ot_env%use_gx_old = .TRUE.
739 4213 : qs_ot_env%use_dx = .TRUE.
740 4213 : IF (qs_ot_env%settings%ot_method == "LBFG" .AND. m_diis < 1) CPABORT("m_diis less than one")
741 : CASE ("DIIS", "BROY")
742 5596 : IF (m_diis < 1) CPABORT("m_diis less than one")
743 : CASE DEFAULT
744 9827 : CPABORT("Unknown option")
745 : END SELECT
746 :
747 9827 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
748 : qs_ot_env%settings%ot_method == "BROY") THEN
749 22384 : ALLOCATE (qs_ot_env%ls_diis(m_diis + 1, m_diis + 1))
750 410582 : qs_ot_env%ls_diis = 0.0_dp
751 16788 : ALLOCATE (qs_ot_env%lss_diis(m_diis + 1, m_diis + 1))
752 16788 : ALLOCATE (qs_ot_env%c_diis(m_diis + 1))
753 16788 : ALLOCATE (qs_ot_env%c_broy(m_diis))
754 11192 : ALLOCATE (qs_ot_env%energy_h(m_diis))
755 16788 : ALLOCATE (qs_ot_env%ipivot(m_diis + 1))
756 : END IF
757 9827 : IF (qs_ot_env%settings%ot_method == "LBFG") THEN
758 : ALLOCATE (qs_ot_env%lbfgs_rho(m_diis), qs_ot_env%lbfgs_yy(m_diis), &
759 672 : qs_ot_env%lbfgs_sy_rotation(m_diis), qs_ot_env%lbfgs_yy_rotation(m_diis))
760 772 : qs_ot_env%lbfgs_rho = 0.0_dp
761 772 : qs_ot_env%lbfgs_yy = 0.0_dp
762 772 : qs_ot_env%lbfgs_sy_rotation = 0.0_dp
763 772 : qs_ot_env%lbfgs_yy_rotation = 0.0_dp
764 : END IF
765 :
766 29175 : ALLOCATE (qs_ot_env%evals(k))
767 19348 : ALLOCATE (qs_ot_env%dum(k))
768 :
769 9827 : NULLIFY (qs_ot_env%matrix_os)
770 9827 : NULLIFY (qs_ot_env%matrix_os_im)
771 9827 : NULLIFY (qs_ot_env%matrix_buf1_ortho)
772 9827 : NULLIFY (qs_ot_env%matrix_buf1_ortho_im)
773 9827 : NULLIFY (qs_ot_env%matrix_buf2_ortho)
774 9827 : NULLIFY (qs_ot_env%matrix_buf2_ortho_im)
775 9827 : NULLIFY (qs_ot_env%matrix_tmp_ortho)
776 9827 : NULLIFY (qs_ot_env%matrix_buf_nk)
777 9827 : NULLIFY (qs_ot_env%matrix_buf_nk_im)
778 9827 : NULLIFY (qs_ot_env%matrix_tmp_nk)
779 9827 : NULLIFY (qs_ot_env%matrix_p)
780 9827 : NULLIFY (qs_ot_env%matrix_p_im)
781 9827 : NULLIFY (qs_ot_env%matrix_r)
782 9827 : NULLIFY (qs_ot_env%matrix_r_im)
783 9827 : NULLIFY (qs_ot_env%matrix_sinp)
784 9827 : NULLIFY (qs_ot_env%matrix_sinp_im)
785 9827 : NULLIFY (qs_ot_env%matrix_cosp)
786 9827 : NULLIFY (qs_ot_env%matrix_cosp_im)
787 9827 : NULLIFY (qs_ot_env%matrix_sinp_b)
788 9827 : NULLIFY (qs_ot_env%matrix_cosp_b)
789 9827 : NULLIFY (qs_ot_env%matrix_buf1)
790 9827 : NULLIFY (qs_ot_env%matrix_buf1_im)
791 9827 : NULLIFY (qs_ot_env%matrix_buf2)
792 9827 : NULLIFY (qs_ot_env%matrix_buf2_im)
793 9827 : NULLIFY (qs_ot_env%matrix_buf3)
794 9827 : NULLIFY (qs_ot_env%matrix_buf3_im)
795 9827 : NULLIFY (qs_ot_env%matrix_buf4)
796 9827 : NULLIFY (qs_ot_env%matrix_buf4_im)
797 9827 : NULLIFY (qs_ot_env%matrix_c0)
798 9827 : NULLIFY (qs_ot_env%matrix_sc0)
799 9827 : NULLIFY (qs_ot_env%matrix_c0_im)
800 9827 : NULLIFY (qs_ot_env%matrix_sc0_im)
801 9827 : NULLIFY (qs_ot_env%matrix_x)
802 9827 : NULLIFY (qs_ot_env%matrix_sx)
803 9827 : NULLIFY (qs_ot_env%matrix_x_im)
804 9827 : NULLIFY (qs_ot_env%matrix_sx_im)
805 9827 : NULLIFY (qs_ot_env%matrix_gx)
806 9827 : NULLIFY (qs_ot_env%matrix_gx_im)
807 9827 : NULLIFY (qs_ot_env%matrix_preconditioned_gx)
808 9827 : NULLIFY (qs_ot_env%matrix_preconditioned_gx_im)
809 9827 : NULLIFY (qs_ot_env%matrix_response_gx)
810 9827 : NULLIFY (qs_ot_env%matrix_response_gx_im)
811 9827 : NULLIFY (qs_ot_env%matrix_mermin_g0)
812 9827 : NULLIFY (qs_ot_env%matrix_mermin_g0_im)
813 9827 : NULLIFY (qs_ot_env%matrix_ref_inv_sqrt)
814 9827 : NULLIFY (qs_ot_env%matrix_ref_inv_sqrt_im)
815 9827 : NULLIFY (qs_ot_env%fm_mermin_c0)
816 9827 : NULLIFY (qs_ot_env%fm_mermin_c0_im)
817 9827 : NULLIFY (qs_ot_env%fm_mermin_h0)
818 9827 : NULLIFY (qs_ot_env%fm_mermin_h0_im)
819 9827 : NULLIFY (qs_ot_env%fm_mermin_c_previous)
820 9827 : NULLIFY (qs_ot_env%fm_mermin_c_previous_im)
821 9827 : NULLIFY (qs_ot_env%fm_mermin_y_previous)
822 9827 : NULLIFY (qs_ot_env%fm_mermin_y_previous_im)
823 9827 : qs_ot_env%mermin_physical_ref_valid = .FALSE.
824 9827 : qs_ot_env%mermin_physical_secant_valid = .FALSE.
825 9827 : qs_ot_env%mermin_gradient_ref_valid = .FALSE.
826 9827 : NULLIFY (qs_ot_env%matrix_gx_old)
827 9827 : NULLIFY (qs_ot_env%matrix_gx_old_im)
828 9827 : NULLIFY (qs_ot_env%matrix_dx)
829 9827 : NULLIFY (qs_ot_env%matrix_dx_im)
830 9827 : NULLIFY (qs_ot_env%buf1_k_k_nosym)
831 9827 : NULLIFY (qs_ot_env%buf2_k_k_nosym)
832 9827 : NULLIFY (qs_ot_env%buf3_k_k_nosym)
833 9827 : NULLIFY (qs_ot_env%buf4_k_k_nosym)
834 9827 : NULLIFY (qs_ot_env%buf1_k_k_sym)
835 9827 : NULLIFY (qs_ot_env%buf2_k_k_sym)
836 9827 : NULLIFY (qs_ot_env%buf3_k_k_sym)
837 9827 : NULLIFY (qs_ot_env%buf4_k_k_sym)
838 9827 : NULLIFY (qs_ot_env%buf1_n_k)
839 9827 : NULLIFY (qs_ot_env%buf1_n_k_dp)
840 9827 : NULLIFY (qs_ot_env%p_k_k_sym)
841 :
842 : ! COMMON MATRICES
843 9827 : CALL dbcsr_init_p(qs_ot_env%matrix_c0)
844 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_c0, template=matrix_s, n=k, &
845 9827 : sym=dbcsr_type_no_symmetry)
846 :
847 9827 : CALL dbcsr_init_p(qs_ot_env%matrix_sc0)
848 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_sc0, template=matrix_s, n=my_ortho_k, &
849 9827 : sym=dbcsr_type_no_symmetry)
850 :
851 9827 : CALL dbcsr_init_p(qs_ot_env%matrix_x)
852 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_x, template=matrix_s, n=k, &
853 9827 : sym=dbcsr_type_no_symmetry)
854 :
855 9827 : CALL dbcsr_init_p(qs_ot_env%matrix_sx)
856 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_sx, template=matrix_s, n=k, &
857 9827 : sym=dbcsr_type_no_symmetry)
858 :
859 9827 : CALL dbcsr_init_p(qs_ot_env%matrix_gx)
860 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_gx, template=matrix_s, n=k, &
861 9827 : sym=dbcsr_type_no_symmetry)
862 :
863 9827 : IF (qs_ot_env%settings%occupation_preconditioner) THEN
864 112 : CALL dbcsr_init_p(qs_ot_env%matrix_preconditioned_gx)
865 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_preconditioned_gx, &
866 : template=matrix_s, n=k, &
867 112 : sym=dbcsr_type_no_symmetry)
868 : END IF
869 :
870 9827 : IF (qs_ot_env%use_dx) THEN
871 4213 : CALL dbcsr_init_p(qs_ot_env%matrix_dx)
872 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_dx, template=matrix_s, n=k, &
873 4213 : sym=dbcsr_type_no_symmetry)
874 : END IF
875 :
876 9827 : IF (qs_ot_env%use_gx_old) THEN
877 4213 : CALL dbcsr_init_p(qs_ot_env%matrix_gx_old)
878 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_gx_old, template=matrix_s, n=k, &
879 4213 : sym=dbcsr_type_no_symmetry)
880 : END IF
881 :
882 8524 : SELECT CASE (qs_ot_env%settings%ot_algorithm)
883 : CASE ("TOD")
884 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_p)
885 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_p, template=matrix_s, m=k, n=k, &
886 8524 : sym=dbcsr_type_no_symmetry)
887 :
888 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_r)
889 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_r, template=matrix_s, m=k, n=k, &
890 8524 : sym=dbcsr_type_no_symmetry)
891 :
892 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_sinp)
893 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_sinp, template=matrix_s, m=k, n=k, &
894 8524 : sym=dbcsr_type_no_symmetry)
895 :
896 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_cosp)
897 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_cosp, template=matrix_s, m=k, n=k, &
898 8524 : sym=dbcsr_type_no_symmetry)
899 :
900 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_sinp_b)
901 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_sinp_b, template=matrix_s, m=k, n=k, &
902 8524 : sym=dbcsr_type_no_symmetry)
903 :
904 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_cosp_b)
905 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_cosp_b, template=matrix_s, m=k, n=k, &
906 8524 : sym=dbcsr_type_no_symmetry)
907 :
908 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_buf1)
909 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf1, template=matrix_s, m=k, n=k, &
910 8524 : sym=dbcsr_type_no_symmetry)
911 :
912 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_buf2)
913 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf2, template=matrix_s, m=k, n=k, &
914 8524 : sym=dbcsr_type_no_symmetry)
915 :
916 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_buf3)
917 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf3, template=matrix_s, m=k, n=k, &
918 8524 : sym=dbcsr_type_no_symmetry)
919 :
920 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_buf4)
921 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf4, template=matrix_s, m=k, n=k, &
922 8524 : sym=dbcsr_type_no_symmetry)
923 :
924 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_os)
925 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_os, template=matrix_s, m=my_ortho_k, n=my_ortho_k, &
926 8524 : sym=dbcsr_type_no_symmetry)
927 :
928 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_buf1_ortho)
929 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf1_ortho, template=matrix_s, m=my_ortho_k, n=k, &
930 8524 : sym=dbcsr_type_no_symmetry)
931 :
932 8524 : CALL dbcsr_init_p(qs_ot_env%matrix_buf2_ortho)
933 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf2_ortho, template=matrix_s, m=my_ortho_k, n=k, &
934 8524 : sym=dbcsr_type_no_symmetry)
935 :
936 : CASE ("REF")
937 1303 : CALL dbcsr_init_p(qs_ot_env%buf1_k_k_nosym)
938 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf1_k_k_nosym, template=matrix_s, m=k, n=k, &
939 1303 : sym=dbcsr_type_no_symmetry)
940 :
941 1303 : CALL dbcsr_init_p(qs_ot_env%buf2_k_k_nosym)
942 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf2_k_k_nosym, template=matrix_s, m=k, n=k, &
943 1303 : sym=dbcsr_type_no_symmetry)
944 :
945 1303 : CALL dbcsr_init_p(qs_ot_env%buf3_k_k_nosym)
946 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf3_k_k_nosym, template=matrix_s, m=k, n=k, &
947 1303 : sym=dbcsr_type_no_symmetry)
948 :
949 1303 : CALL dbcsr_init_p(qs_ot_env%buf4_k_k_nosym)
950 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf4_k_k_nosym, template=matrix_s, m=k, n=k, &
951 1303 : sym=dbcsr_type_no_symmetry)
952 :
953 : ! It claims to be symmetric but to avoid dbcsr confusion nonsymmetric is kept
954 1303 : CALL dbcsr_init_p(qs_ot_env%buf1_k_k_sym)
955 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf1_k_k_sym, template=matrix_s, m=k, n=k, &
956 1303 : sym=dbcsr_type_no_symmetry)
957 :
958 1303 : CALL dbcsr_init_p(qs_ot_env%buf2_k_k_sym)
959 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf2_k_k_sym, template=matrix_s, m=k, n=k, &
960 1303 : sym=dbcsr_type_no_symmetry)
961 :
962 1303 : CALL dbcsr_init_p(qs_ot_env%buf3_k_k_sym)
963 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf3_k_k_sym, template=matrix_s, m=k, n=k, &
964 1303 : sym=dbcsr_type_no_symmetry)
965 : !
966 1303 : CALL dbcsr_init_p(qs_ot_env%buf4_k_k_sym)
967 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf4_k_k_sym, template=matrix_s, m=k, n=k, &
968 1303 : sym=dbcsr_type_no_symmetry)
969 : !
970 1303 : CALL dbcsr_init_p(qs_ot_env%p_k_k_sym)
971 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%p_k_k_sym, template=matrix_s, m=k, n=k, &
972 1303 : sym=dbcsr_type_no_symmetry)
973 : !
974 1303 : CALL dbcsr_init_p(qs_ot_env%buf1_n_k)
975 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%buf1_n_k, template=matrix_s, n=k, &
976 1303 : sym=dbcsr_type_no_symmetry)
977 : !
978 1303 : CALL dbcsr_init_p(qs_ot_env%matrix_buf1)
979 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf1, template=matrix_s, m=k, n=k, &
980 11130 : sym=dbcsr_type_no_symmetry)
981 :
982 : END SELECT
983 :
984 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
985 9827 : qs_ot_env%settings%ot_method == "BROY" .OR. &
986 : qs_ot_env%settings%ot_method == "LBFG") THEN
987 5708 : NULLIFY (qs_ot_env%matrix_h_e)
988 5708 : NULLIFY (qs_ot_env%matrix_h_x)
989 5708 : NULLIFY (qs_ot_env%matrix_h_e_im)
990 5708 : NULLIFY (qs_ot_env%matrix_h_x_im)
991 5708 : ncoef = m_diis
992 5708 : IF (qs_ot_env%settings%ot_method == "LBFG") ncoef = m_diis + 1
993 5708 : CALL dbcsr_allocate_matrix_set(qs_ot_env%matrix_h_e, ncoef)
994 5708 : CALL dbcsr_allocate_matrix_set(qs_ot_env%matrix_h_x, ncoef)
995 45490 : DO i = 1, ncoef
996 39782 : CALL dbcsr_init_p(qs_ot_env%matrix_h_x(i)%matrix)
997 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_h_x(i)%matrix, template=matrix_s, n=k, &
998 39782 : sym=dbcsr_type_no_symmetry)
999 :
1000 39782 : CALL dbcsr_init_p(qs_ot_env%matrix_h_e(i)%matrix)
1001 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_h_e(i)%matrix, template=matrix_s, n=k, &
1002 49609 : sym=dbcsr_type_no_symmetry)
1003 : END DO
1004 : END IF
1005 :
1006 9827 : NULLIFY (qs_ot_env%rot_mat_u, qs_ot_env%rot_mat_u_im, &
1007 9827 : qs_ot_env%rot_mat_x, qs_ot_env%rot_mat_x_im, &
1008 9827 : qs_ot_env%rot_mat_h_e, qs_ot_env%rot_mat_h_x, &
1009 9827 : qs_ot_env%rot_mat_h_e_im, qs_ot_env%rot_mat_h_x_im, qs_ot_env%rot_mat_gx, &
1010 9827 : qs_ot_env%rot_mat_gx_im, qs_ot_env%rot_mat_response_gx, &
1011 9827 : qs_ot_env%rot_mat_response_gx_im, qs_ot_env%rot_mat_mermin_g0, &
1012 9827 : qs_ot_env%rot_mat_mermin_g0_im, qs_ot_env%rot_mat_gx_old, &
1013 9827 : qs_ot_env%rot_mat_gx_old_im, qs_ot_env%rot_mat_dx, qs_ot_env%rot_mat_dx_im, &
1014 9827 : qs_ot_env%rot_mat_evals, qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_dedu_im, &
1015 9827 : qs_ot_env%rot_mat_chc, qs_ot_env%rot_mat_chc_im, &
1016 9827 : qs_ot_env%rot_mat_evec_re, qs_ot_env%rot_mat_evec_im)
1017 :
1018 9827 : IF (qs_ot_env%settings%do_rotation) THEN
1019 386 : CALL dbcsr_init_p(qs_ot_env%rot_mat_u)
1020 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_u, template=matrix_s, m=k, n=k, &
1021 386 : sym=dbcsr_type_no_symmetry)
1022 :
1023 386 : CALL dbcsr_init_p(qs_ot_env%rot_mat_x)
1024 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_x, template=matrix_s, m=k, n=k, &
1025 386 : sym=dbcsr_type_no_symmetry)
1026 :
1027 386 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dedu)
1028 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_dedu, template=matrix_s, m=k, n=k, &
1029 386 : sym=dbcsr_type_no_symmetry)
1030 :
1031 386 : CALL dbcsr_init_p(qs_ot_env%rot_mat_chc)
1032 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_chc, template=matrix_s, m=k, n=k, &
1033 386 : sym=dbcsr_type_no_symmetry)
1034 :
1035 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
1036 386 : qs_ot_env%settings%ot_method == "BROY" .OR. &
1037 : qs_ot_env%settings%ot_method == "LBFG") THEN
1038 134 : ncoef = m_diis
1039 134 : IF (qs_ot_env%settings%ot_method == "LBFG") ncoef = m_diis + 1
1040 134 : CALL dbcsr_allocate_matrix_set(qs_ot_env%rot_mat_h_e, ncoef)
1041 134 : CALL dbcsr_allocate_matrix_set(qs_ot_env%rot_mat_h_x, ncoef)
1042 1016 : DO i = 1, ncoef
1043 882 : CALL dbcsr_init_p(qs_ot_env%rot_mat_h_e(i)%matrix)
1044 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_h_e(i)%matrix, template=matrix_s, m=k, n=k, &
1045 882 : sym=dbcsr_type_no_symmetry)
1046 :
1047 882 : CALL dbcsr_init_p(qs_ot_env%rot_mat_h_x(i)%matrix)
1048 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_h_x(i)%matrix, template=matrix_s, m=k, n=k, &
1049 1268 : sym=dbcsr_type_no_symmetry)
1050 : END DO
1051 : END IF
1052 :
1053 1154 : ALLOCATE (qs_ot_env%rot_mat_evals(k))
1054 386 : CALL dbcsr_init_p(qs_ot_env%rot_mat_evec_re)
1055 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_evec_re, template=matrix_s, m=k, n=k, &
1056 386 : sym=dbcsr_type_no_symmetry)
1057 386 : CALL dbcsr_init_p(qs_ot_env%rot_mat_evec_im)
1058 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_evec_im, template=matrix_s, m=k, n=k, &
1059 386 : sym=dbcsr_type_no_symmetry)
1060 :
1061 386 : CALL dbcsr_init_p(qs_ot_env%rot_mat_gx)
1062 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_gx, template=matrix_s, m=k, n=k, &
1063 386 : sym=dbcsr_type_no_symmetry)
1064 :
1065 386 : IF (qs_ot_env%settings%occupation_preconditioner) THEN
1066 112 : CALL dbcsr_init_p(qs_ot_env%rot_mat_response_gx)
1067 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_response_gx, &
1068 : template=matrix_s, m=k, n=k, &
1069 112 : sym=dbcsr_type_no_symmetry)
1070 : END IF
1071 :
1072 386 : IF (qs_ot_env%use_gx_old) THEN
1073 300 : CALL dbcsr_init_p(qs_ot_env%rot_mat_gx_old)
1074 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_gx_old, template=matrix_s, m=k, n=k, &
1075 300 : sym=dbcsr_type_no_symmetry)
1076 : END IF
1077 :
1078 386 : IF (qs_ot_env%use_dx) THEN
1079 300 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dx)
1080 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_dx, template=matrix_s, m=k, n=k, &
1081 300 : sym=dbcsr_type_no_symmetry)
1082 : END IF
1083 :
1084 : END IF
1085 :
1086 9827 : IF (qs_ot_env%settings%do_ener) THEN
1087 : ncoef = my_energy_dimension
1088 342 : ALLOCATE (qs_ot_env%ener_x(ncoef))
1089 228 : ALLOCATE (qs_ot_env%ener_rayleigh(ncoef))
1090 :
1091 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
1092 114 : qs_ot_env%settings%ot_method == "BROY" .OR. &
1093 : qs_ot_env%settings%ot_method == "LBFG") THEN
1094 54 : nhistory = m_diis
1095 54 : IF (qs_ot_env%settings%ot_method == "LBFG") nhistory = m_diis + 1
1096 216 : ALLOCATE (qs_ot_env%ener_h_e(nhistory, ncoef))
1097 162 : ALLOCATE (qs_ot_env%ener_h_x(nhistory, ncoef))
1098 3478 : qs_ot_env%ener_h_e = 0.0_dp
1099 3478 : qs_ot_env%ener_h_x = 0.0_dp
1100 : END IF
1101 :
1102 228 : ALLOCATE (qs_ot_env%ener_gx(ncoef))
1103 114 : IF (qs_ot_env%settings%occupation_preconditioner) THEN
1104 224 : ALLOCATE (qs_ot_env%ener_preconditioned_gx(ncoef))
1105 224 : ALLOCATE (qs_ot_env%ener_response_gx(ncoef))
1106 : END IF
1107 :
1108 114 : IF (qs_ot_env%use_gx_old) THEN
1109 132 : ALLOCATE (qs_ot_env%ener_gx_old(ncoef))
1110 : END IF
1111 :
1112 114 : IF (qs_ot_env%use_dx) THEN
1113 132 : ALLOCATE (qs_ot_env%ener_dx(ncoef))
1114 502 : qs_ot_env%ener_dx = 0.0_dp
1115 : END IF
1116 : END IF
1117 :
1118 9827 : qs_ot_env%state_allocated = .TRUE.
1119 :
1120 9827 : END SUBROUTINE qs_ot_allocate
1121 :
1122 : ! **************************************************************************************************
1123 : !> \brief ...
1124 : !> \param qs_ot_env ...
1125 : !> \param matrix_s ...
1126 : ! **************************************************************************************************
1127 1041 : SUBROUTINE qs_ot_allocate_complex_state(qs_ot_env, matrix_s)
1128 : TYPE(qs_ot_type) :: qs_ot_env
1129 : TYPE(dbcsr_type), POINTER :: matrix_s
1130 :
1131 : INTEGER :: i, ncoef, nmo, nmo_ortho
1132 :
1133 347 : CPASSERT(qs_ot_env%state_allocated)
1134 347 : CPASSERT(.NOT. qs_ot_env%has_complex_kpoint_state)
1135 347 : CPASSERT(ASSOCIATED(matrix_s))
1136 :
1137 347 : CALL dbcsr_get_info(qs_ot_env%matrix_c0, nfullcols_total=nmo)
1138 347 : CALL dbcsr_get_info(qs_ot_env%matrix_sc0, nfullcols_total=nmo_ortho)
1139 :
1140 347 : CALL dbcsr_init_p(qs_ot_env%matrix_c0_im)
1141 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_c0_im, template=matrix_s, n=nmo, &
1142 347 : sym=dbcsr_type_no_symmetry)
1143 347 : CALL dbcsr_init_p(qs_ot_env%matrix_sc0_im)
1144 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_sc0_im, template=matrix_s, n=nmo_ortho, &
1145 347 : sym=dbcsr_type_no_symmetry)
1146 347 : CALL dbcsr_init_p(qs_ot_env%matrix_x_im)
1147 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_x_im, template=matrix_s, n=nmo, &
1148 347 : sym=dbcsr_type_no_symmetry)
1149 347 : CALL dbcsr_init_p(qs_ot_env%matrix_sx_im)
1150 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_sx_im, template=matrix_s, n=nmo, &
1151 347 : sym=dbcsr_type_no_symmetry)
1152 347 : CALL dbcsr_init_p(qs_ot_env%matrix_gx_im)
1153 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_gx_im, template=matrix_s, n=nmo, &
1154 347 : sym=dbcsr_type_no_symmetry)
1155 347 : IF (qs_ot_env%settings%occupation_preconditioner) THEN
1156 112 : CALL dbcsr_init_p(qs_ot_env%matrix_preconditioned_gx_im)
1157 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_preconditioned_gx_im, &
1158 : template=matrix_s, n=nmo, &
1159 112 : sym=dbcsr_type_no_symmetry)
1160 112 : IF (qs_ot_env%use_dx) THEN
1161 64 : CALL dbcsr_init_p(qs_ot_env%matrix_response_gx)
1162 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_response_gx, &
1163 : template=matrix_s, n=nmo, &
1164 64 : sym=dbcsr_type_no_symmetry)
1165 64 : CALL dbcsr_init_p(qs_ot_env%matrix_response_gx_im)
1166 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_response_gx_im, &
1167 : template=matrix_s, n=nmo, &
1168 64 : sym=dbcsr_type_no_symmetry)
1169 64 : CALL dbcsr_init_p(qs_ot_env%matrix_mermin_g0)
1170 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_mermin_g0, &
1171 : template=matrix_s, n=nmo, &
1172 64 : sym=dbcsr_type_no_symmetry)
1173 64 : CALL dbcsr_init_p(qs_ot_env%matrix_mermin_g0_im)
1174 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_mermin_g0_im, &
1175 : template=matrix_s, n=nmo, &
1176 64 : sym=dbcsr_type_no_symmetry)
1177 64 : IF (qs_ot_env%settings%do_ener) THEN
1178 64 : CPASSERT(ASSOCIATED(qs_ot_env%ener_x))
1179 192 : ALLOCATE (qs_ot_env%ener_mermin_g0(SIZE(qs_ot_env%ener_x)))
1180 : END IF
1181 : END IF
1182 : END IF
1183 152 : SELECT CASE (qs_ot_env%settings%ot_algorithm)
1184 : CASE ("TOD")
1185 152 : CPASSERT(nmo_ortho == nmo)
1186 152 : CALL allocate_complex_copy(qs_ot_env%matrix_p_im, qs_ot_env%matrix_p, "matrix_p_im")
1187 152 : CALL allocate_complex_copy(qs_ot_env%matrix_r_im, qs_ot_env%matrix_r, "matrix_r_im")
1188 152 : CALL allocate_complex_copy(qs_ot_env%matrix_sinp_im, qs_ot_env%matrix_sinp, "matrix_sinp_im")
1189 152 : CALL allocate_complex_copy(qs_ot_env%matrix_cosp_im, qs_ot_env%matrix_cosp, "matrix_cosp_im")
1190 152 : CALL allocate_complex_copy(qs_ot_env%matrix_buf1_im, qs_ot_env%matrix_buf1, "matrix_buf1_im")
1191 152 : CALL allocate_complex_copy(qs_ot_env%matrix_buf2_im, qs_ot_env%matrix_buf2, "matrix_buf2_im")
1192 152 : CALL allocate_complex_copy(qs_ot_env%matrix_buf3_im, qs_ot_env%matrix_buf3, "matrix_buf3_im")
1193 152 : CALL allocate_complex_copy(qs_ot_env%matrix_buf4_im, qs_ot_env%matrix_buf4, "matrix_buf4_im")
1194 152 : CALL allocate_complex_copy(qs_ot_env%matrix_os_im, qs_ot_env%matrix_os, "matrix_os_im")
1195 : CALL allocate_complex_copy(qs_ot_env%matrix_buf1_ortho_im, qs_ot_env%matrix_buf1_ortho, &
1196 152 : "matrix_buf1_ortho_im")
1197 : CALL allocate_complex_copy(qs_ot_env%matrix_buf2_ortho_im, qs_ot_env%matrix_buf2_ortho, &
1198 152 : "matrix_buf2_ortho_im")
1199 : CALL allocate_complex_copy(qs_ot_env%matrix_tmp_ortho, qs_ot_env%matrix_buf1_ortho, &
1200 152 : "matrix_tmp_ortho")
1201 152 : CALL allocate_complex_copy(qs_ot_env%matrix_buf_nk, qs_ot_env%matrix_x, "matrix_buf_nk")
1202 152 : CALL allocate_complex_copy(qs_ot_env%matrix_buf_nk_im, qs_ot_env%matrix_x, "matrix_buf_nk_im")
1203 152 : CALL allocate_complex_copy(qs_ot_env%matrix_tmp_nk, qs_ot_env%matrix_x, "matrix_tmp_nk")
1204 : CASE ("REF")
1205 195 : CALL dbcsr_init_p(qs_ot_env%matrix_ref_inv_sqrt)
1206 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_ref_inv_sqrt, template=matrix_s, m=nmo, n=nmo, &
1207 195 : sym=dbcsr_type_no_symmetry)
1208 195 : CALL dbcsr_init_p(qs_ot_env%matrix_ref_inv_sqrt_im)
1209 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_ref_inv_sqrt_im, template=matrix_s, m=nmo, n=nmo, &
1210 195 : sym=dbcsr_type_no_symmetry)
1211 195 : CALL dbcsr_init_p(qs_ot_env%buf1_n_k_dp)
1212 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%buf1_n_k_dp, template=matrix_s, n=nmo, &
1213 195 : sym=dbcsr_type_no_symmetry)
1214 : CASE DEFAULT
1215 347 : CPABORT("Complex K-point OT state requires ALGORITHM STRICT or IRAC")
1216 : END SELECT
1217 347 : IF (qs_ot_env%use_dx) THEN
1218 212 : CALL dbcsr_init_p(qs_ot_env%matrix_dx_im)
1219 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_dx_im, template=matrix_s, n=nmo, &
1220 212 : sym=dbcsr_type_no_symmetry)
1221 : END IF
1222 347 : IF (qs_ot_env%use_gx_old) THEN
1223 212 : CALL dbcsr_init_p(qs_ot_env%matrix_gx_old_im)
1224 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_gx_old_im, template=matrix_s, n=nmo, &
1225 212 : sym=dbcsr_type_no_symmetry)
1226 : END IF
1227 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
1228 347 : qs_ot_env%settings%ot_method == "BROY" .OR. &
1229 : qs_ot_env%settings%ot_method == "LBFG") THEN
1230 227 : ncoef = qs_ot_env%settings%diis_m
1231 227 : IF (qs_ot_env%settings%ot_method == "LBFG") ncoef = ncoef + 1
1232 227 : CALL dbcsr_allocate_matrix_set(qs_ot_env%matrix_h_e_im, ncoef)
1233 227 : CALL dbcsr_allocate_matrix_set(qs_ot_env%matrix_h_x_im, ncoef)
1234 1692 : DO i = 1, ncoef
1235 1465 : CALL dbcsr_init_p(qs_ot_env%matrix_h_e_im(i)%matrix)
1236 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_h_e_im(i)%matrix, &
1237 : template=matrix_s, n=nmo, &
1238 1465 : sym=dbcsr_type_no_symmetry)
1239 1465 : CALL dbcsr_init_p(qs_ot_env%matrix_h_x_im(i)%matrix)
1240 : CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_h_x_im(i)%matrix, &
1241 : template=matrix_s, n=nmo, &
1242 1812 : sym=dbcsr_type_no_symmetry)
1243 : END DO
1244 : END IF
1245 :
1246 347 : IF (qs_ot_env%settings%do_rotation) THEN
1247 166 : CALL dbcsr_init_p(qs_ot_env%rot_mat_u_im)
1248 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_u_im, template=matrix_s, m=nmo, n=nmo, &
1249 166 : sym=dbcsr_type_no_symmetry)
1250 166 : CALL dbcsr_init_p(qs_ot_env%rot_mat_x_im)
1251 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_x_im, template=matrix_s, m=nmo, n=nmo, &
1252 166 : sym=dbcsr_type_no_symmetry)
1253 166 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dedu_im)
1254 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_dedu_im, template=matrix_s, m=nmo, n=nmo, &
1255 166 : sym=dbcsr_type_no_symmetry)
1256 166 : CALL dbcsr_init_p(qs_ot_env%rot_mat_chc_im)
1257 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_chc_im, template=matrix_s, m=nmo, n=nmo, &
1258 166 : sym=dbcsr_type_no_symmetry)
1259 166 : CALL dbcsr_init_p(qs_ot_env%rot_mat_gx_im)
1260 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_gx_im, template=matrix_s, m=nmo, n=nmo, &
1261 166 : sym=dbcsr_type_no_symmetry)
1262 166 : IF (qs_ot_env%settings%occupation_preconditioner) THEN
1263 112 : CALL dbcsr_init_p(qs_ot_env%rot_mat_response_gx_im)
1264 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_response_gx_im, &
1265 : template=matrix_s, m=nmo, n=nmo, &
1266 112 : sym=dbcsr_type_no_symmetry)
1267 112 : IF (qs_ot_env%use_dx) THEN
1268 64 : CALL dbcsr_init_p(qs_ot_env%rot_mat_mermin_g0)
1269 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_mermin_g0, &
1270 : template=matrix_s, m=nmo, n=nmo, &
1271 64 : sym=dbcsr_type_no_symmetry)
1272 64 : CALL dbcsr_init_p(qs_ot_env%rot_mat_mermin_g0_im)
1273 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_mermin_g0_im, &
1274 : template=matrix_s, m=nmo, n=nmo, &
1275 64 : sym=dbcsr_type_no_symmetry)
1276 : END IF
1277 : END IF
1278 166 : IF (qs_ot_env%use_dx) THEN
1279 110 : CALL dbcsr_init_p(qs_ot_env%rot_mat_dx_im)
1280 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_dx_im, template=matrix_s, m=nmo, n=nmo, &
1281 110 : sym=dbcsr_type_no_symmetry)
1282 : END IF
1283 166 : IF (qs_ot_env%use_gx_old) THEN
1284 110 : CALL dbcsr_init_p(qs_ot_env%rot_mat_gx_old_im)
1285 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_gx_old_im, template=matrix_s, m=nmo, n=nmo, &
1286 110 : sym=dbcsr_type_no_symmetry)
1287 : END IF
1288 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
1289 166 : qs_ot_env%settings%ot_method == "BROY" .OR. &
1290 : qs_ot_env%settings%ot_method == "LBFG") THEN
1291 100 : ncoef = qs_ot_env%settings%diis_m
1292 100 : IF (qs_ot_env%settings%ot_method == "LBFG") ncoef = ncoef + 1
1293 100 : CALL dbcsr_allocate_matrix_set(qs_ot_env%rot_mat_h_e_im, ncoef)
1294 100 : CALL dbcsr_allocate_matrix_set(qs_ot_env%rot_mat_h_x_im, ncoef)
1295 740 : DO i = 1, ncoef
1296 640 : CALL dbcsr_init_p(qs_ot_env%rot_mat_h_e_im(i)%matrix)
1297 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_h_e_im(i)%matrix, &
1298 : template=matrix_s, m=nmo, n=nmo, &
1299 640 : sym=dbcsr_type_no_symmetry)
1300 640 : CALL dbcsr_init_p(qs_ot_env%rot_mat_h_x_im(i)%matrix)
1301 : CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_h_x_im(i)%matrix, &
1302 : template=matrix_s, m=nmo, n=nmo, &
1303 806 : sym=dbcsr_type_no_symmetry)
1304 : END DO
1305 : END IF
1306 : END IF
1307 :
1308 347 : qs_ot_env%has_complex_kpoint_state = .TRUE.
1309 :
1310 : CONTAINS
1311 :
1312 : ! **************************************************************************************************
1313 : !> \brief allocate a zeroed matrix with the distribution and sparsity of a real companion
1314 : !> \param matrix output matrix
1315 : !> \param template real companion matrix
1316 : !> \param name DBCSR matrix name
1317 : ! **************************************************************************************************
1318 2280 : SUBROUTINE allocate_complex_copy(matrix, template, name)
1319 : TYPE(dbcsr_type), POINTER :: matrix, template
1320 : CHARACTER(LEN=*), INTENT(IN) :: name
1321 :
1322 2280 : CALL dbcsr_init_p(matrix)
1323 2280 : CALL dbcsr_copy(matrix, template, name=name)
1324 2280 : CALL dbcsr_set(matrix, 0.0_dp)
1325 2280 : END SUBROUTINE allocate_complex_copy
1326 :
1327 : END SUBROUTINE qs_ot_allocate_complex_state
1328 :
1329 : ! **************************************************************************************************
1330 : !> \brief deallocates data
1331 : !> \param qs_ot_env ...
1332 : ! **************************************************************************************************
1333 9867 : SUBROUTINE qs_ot_destroy(qs_ot_env)
1334 : TYPE(qs_ot_type) :: qs_ot_env
1335 :
1336 9867 : IF (.NOT. qs_ot_env%state_allocated) RETURN
1337 :
1338 9827 : CALL mp_para_env_release(qs_ot_env%para_env)
1339 9827 : CALL cp_blacs_env_release(qs_ot_env%blacs_env)
1340 :
1341 9827 : IF (ASSOCIATED(qs_ot_env%evals)) DEALLOCATE (qs_ot_env%evals)
1342 9827 : IF (ASSOCIATED(qs_ot_env%dum)) DEALLOCATE (qs_ot_env%dum)
1343 :
1344 9827 : IF (ASSOCIATED(qs_ot_env%matrix_os)) CALL dbcsr_release_p(qs_ot_env%matrix_os)
1345 9827 : IF (ASSOCIATED(qs_ot_env%matrix_os_im)) CALL dbcsr_release_p(qs_ot_env%matrix_os_im)
1346 9827 : IF (ASSOCIATED(qs_ot_env%matrix_p)) CALL dbcsr_release_p(qs_ot_env%matrix_p)
1347 9827 : IF (ASSOCIATED(qs_ot_env%matrix_p_im)) CALL dbcsr_release_p(qs_ot_env%matrix_p_im)
1348 9827 : IF (ASSOCIATED(qs_ot_env%matrix_cosp)) CALL dbcsr_release_p(qs_ot_env%matrix_cosp)
1349 9827 : IF (ASSOCIATED(qs_ot_env%matrix_cosp_im)) CALL dbcsr_release_p(qs_ot_env%matrix_cosp_im)
1350 9827 : IF (ASSOCIATED(qs_ot_env%matrix_sinp)) CALL dbcsr_release_p(qs_ot_env%matrix_sinp)
1351 9827 : IF (ASSOCIATED(qs_ot_env%matrix_sinp_im)) CALL dbcsr_release_p(qs_ot_env%matrix_sinp_im)
1352 9827 : IF (ASSOCIATED(qs_ot_env%matrix_r)) CALL dbcsr_release_p(qs_ot_env%matrix_r)
1353 9827 : IF (ASSOCIATED(qs_ot_env%matrix_r_im)) CALL dbcsr_release_p(qs_ot_env%matrix_r_im)
1354 9827 : IF (ASSOCIATED(qs_ot_env%matrix_cosp_b)) CALL dbcsr_release_p(qs_ot_env%matrix_cosp_b)
1355 9827 : IF (ASSOCIATED(qs_ot_env%matrix_sinp_b)) CALL dbcsr_release_p(qs_ot_env%matrix_sinp_b)
1356 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf1)) CALL dbcsr_release_p(qs_ot_env%matrix_buf1)
1357 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf1_im)) CALL dbcsr_release_p(qs_ot_env%matrix_buf1_im)
1358 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf2)) CALL dbcsr_release_p(qs_ot_env%matrix_buf2)
1359 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf2_im)) CALL dbcsr_release_p(qs_ot_env%matrix_buf2_im)
1360 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf3)) CALL dbcsr_release_p(qs_ot_env%matrix_buf3)
1361 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf3_im)) CALL dbcsr_release_p(qs_ot_env%matrix_buf3_im)
1362 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf4)) CALL dbcsr_release_p(qs_ot_env%matrix_buf4)
1363 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf4_im)) CALL dbcsr_release_p(qs_ot_env%matrix_buf4_im)
1364 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf1_ortho)) CALL dbcsr_release_p(qs_ot_env%matrix_buf1_ortho)
1365 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf1_ortho_im)) CALL dbcsr_release_p(qs_ot_env%matrix_buf1_ortho_im)
1366 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf2_ortho)) CALL dbcsr_release_p(qs_ot_env%matrix_buf2_ortho)
1367 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf2_ortho_im)) CALL dbcsr_release_p(qs_ot_env%matrix_buf2_ortho_im)
1368 9827 : IF (ASSOCIATED(qs_ot_env%matrix_tmp_ortho)) CALL dbcsr_release_p(qs_ot_env%matrix_tmp_ortho)
1369 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf_nk)) CALL dbcsr_release_p(qs_ot_env%matrix_buf_nk)
1370 9827 : IF (ASSOCIATED(qs_ot_env%matrix_buf_nk_im)) CALL dbcsr_release_p(qs_ot_env%matrix_buf_nk_im)
1371 9827 : IF (ASSOCIATED(qs_ot_env%matrix_tmp_nk)) CALL dbcsr_release_p(qs_ot_env%matrix_tmp_nk)
1372 9827 : IF (ASSOCIATED(qs_ot_env%matrix_c0)) CALL dbcsr_release_p(qs_ot_env%matrix_c0)
1373 9827 : IF (ASSOCIATED(qs_ot_env%matrix_sc0)) CALL dbcsr_release_p(qs_ot_env%matrix_sc0)
1374 9827 : IF (ASSOCIATED(qs_ot_env%matrix_psc0)) CALL dbcsr_release_p(qs_ot_env%matrix_psc0)
1375 9827 : IF (ASSOCIATED(qs_ot_env%matrix_psc0_im)) CALL dbcsr_release_p(qs_ot_env%matrix_psc0_im)
1376 9827 : IF (ASSOCIATED(qs_ot_env%matrix_c0_im)) CALL dbcsr_release_p(qs_ot_env%matrix_c0_im)
1377 9827 : IF (ASSOCIATED(qs_ot_env%matrix_sc0_im)) CALL dbcsr_release_p(qs_ot_env%matrix_sc0_im)
1378 9827 : IF (ASSOCIATED(qs_ot_env%matrix_x)) CALL dbcsr_release_p(qs_ot_env%matrix_x)
1379 9827 : IF (ASSOCIATED(qs_ot_env%matrix_sx)) CALL dbcsr_release_p(qs_ot_env%matrix_sx)
1380 9827 : IF (ASSOCIATED(qs_ot_env%matrix_x_im)) CALL dbcsr_release_p(qs_ot_env%matrix_x_im)
1381 9827 : IF (ASSOCIATED(qs_ot_env%matrix_sx_im)) CALL dbcsr_release_p(qs_ot_env%matrix_sx_im)
1382 9827 : IF (ASSOCIATED(qs_ot_env%matrix_gx)) CALL dbcsr_release_p(qs_ot_env%matrix_gx)
1383 9827 : IF (ASSOCIATED(qs_ot_env%matrix_gx_im)) CALL dbcsr_release_p(qs_ot_env%matrix_gx_im)
1384 9827 : IF (ASSOCIATED(qs_ot_env%matrix_preconditioned_gx)) THEN
1385 112 : CALL dbcsr_release_p(qs_ot_env%matrix_preconditioned_gx)
1386 : END IF
1387 9827 : IF (ASSOCIATED(qs_ot_env%matrix_preconditioned_gx_im)) THEN
1388 112 : CALL dbcsr_release_p(qs_ot_env%matrix_preconditioned_gx_im)
1389 : END IF
1390 9827 : IF (ASSOCIATED(qs_ot_env%matrix_response_gx)) CALL dbcsr_release_p(qs_ot_env%matrix_response_gx)
1391 9827 : IF (ASSOCIATED(qs_ot_env%matrix_response_gx_im)) THEN
1392 64 : CALL dbcsr_release_p(qs_ot_env%matrix_response_gx_im)
1393 : END IF
1394 9827 : IF (ASSOCIATED(qs_ot_env%matrix_mermin_g0)) CALL dbcsr_release_p(qs_ot_env%matrix_mermin_g0)
1395 9827 : IF (ASSOCIATED(qs_ot_env%matrix_mermin_g0_im)) THEN
1396 64 : CALL dbcsr_release_p(qs_ot_env%matrix_mermin_g0_im)
1397 : END IF
1398 9827 : IF (ASSOCIATED(qs_ot_env%matrix_ref_inv_sqrt)) CALL dbcsr_release_p(qs_ot_env%matrix_ref_inv_sqrt)
1399 9827 : IF (ASSOCIATED(qs_ot_env%matrix_ref_inv_sqrt_im)) CALL dbcsr_release_p(qs_ot_env%matrix_ref_inv_sqrt_im)
1400 9827 : IF (ASSOCIATED(qs_ot_env%fm_mermin_c0)) THEN
1401 104 : CALL cp_fm_release(qs_ot_env%fm_mermin_c0)
1402 104 : DEALLOCATE (qs_ot_env%fm_mermin_c0)
1403 : END IF
1404 9827 : IF (ASSOCIATED(qs_ot_env%fm_mermin_c0_im)) THEN
1405 104 : CALL cp_fm_release(qs_ot_env%fm_mermin_c0_im)
1406 104 : DEALLOCATE (qs_ot_env%fm_mermin_c0_im)
1407 : END IF
1408 9827 : IF (ASSOCIATED(qs_ot_env%fm_mermin_h0)) THEN
1409 104 : CALL cp_fm_release(qs_ot_env%fm_mermin_h0)
1410 104 : DEALLOCATE (qs_ot_env%fm_mermin_h0)
1411 : END IF
1412 9827 : IF (ASSOCIATED(qs_ot_env%fm_mermin_h0_im)) THEN
1413 104 : CALL cp_fm_release(qs_ot_env%fm_mermin_h0_im)
1414 104 : DEALLOCATE (qs_ot_env%fm_mermin_h0_im)
1415 : END IF
1416 9827 : IF (ASSOCIATED(qs_ot_env%fm_mermin_c_previous)) THEN
1417 96 : CALL cp_fm_release(qs_ot_env%fm_mermin_c_previous)
1418 96 : DEALLOCATE (qs_ot_env%fm_mermin_c_previous)
1419 : END IF
1420 9827 : IF (ASSOCIATED(qs_ot_env%fm_mermin_c_previous_im)) THEN
1421 96 : CALL cp_fm_release(qs_ot_env%fm_mermin_c_previous_im)
1422 96 : DEALLOCATE (qs_ot_env%fm_mermin_c_previous_im)
1423 : END IF
1424 9827 : IF (ASSOCIATED(qs_ot_env%fm_mermin_y_previous)) THEN
1425 96 : CALL cp_fm_release(qs_ot_env%fm_mermin_y_previous)
1426 96 : DEALLOCATE (qs_ot_env%fm_mermin_y_previous)
1427 : END IF
1428 9827 : IF (ASSOCIATED(qs_ot_env%fm_mermin_y_previous_im)) THEN
1429 96 : CALL cp_fm_release(qs_ot_env%fm_mermin_y_previous_im)
1430 96 : DEALLOCATE (qs_ot_env%fm_mermin_y_previous_im)
1431 : END IF
1432 9827 : IF (ALLOCATED(qs_ot_env%mermin_occupation0)) DEALLOCATE (qs_ot_env%mermin_occupation0)
1433 9827 : IF (ALLOCATED(qs_ot_env%mermin_occupation_previous)) THEN
1434 96 : DEALLOCATE (qs_ot_env%mermin_occupation_previous)
1435 : END IF
1436 9827 : IF (ALLOCATED(qs_ot_env%ener_mermin_g0)) DEALLOCATE (qs_ot_env%ener_mermin_g0)
1437 9827 : qs_ot_env%mermin_physical_ref_valid = .FALSE.
1438 9827 : qs_ot_env%mermin_physical_secant_valid = .FALSE.
1439 9827 : qs_ot_env%mermin_gradient_ref_valid = .FALSE.
1440 9827 : IF (ASSOCIATED(qs_ot_env%matrix_dx)) CALL dbcsr_release_p(qs_ot_env%matrix_dx)
1441 9827 : IF (ASSOCIATED(qs_ot_env%matrix_dx_im)) CALL dbcsr_release_p(qs_ot_env%matrix_dx_im)
1442 9827 : IF (ASSOCIATED(qs_ot_env%matrix_gx_old)) CALL dbcsr_release_p(qs_ot_env%matrix_gx_old)
1443 9827 : IF (ASSOCIATED(qs_ot_env%matrix_gx_old_im)) CALL dbcsr_release_p(qs_ot_env%matrix_gx_old_im)
1444 9827 : IF (ASSOCIATED(qs_ot_env%buf1_k_k_nosym)) CALL dbcsr_release_p(qs_ot_env%buf1_k_k_nosym)
1445 9827 : IF (ASSOCIATED(qs_ot_env%buf2_k_k_nosym)) CALL dbcsr_release_p(qs_ot_env%buf2_k_k_nosym)
1446 9827 : IF (ASSOCIATED(qs_ot_env%buf3_k_k_nosym)) CALL dbcsr_release_p(qs_ot_env%buf3_k_k_nosym)
1447 9827 : IF (ASSOCIATED(qs_ot_env%buf4_k_k_nosym)) CALL dbcsr_release_p(qs_ot_env%buf4_k_k_nosym)
1448 9827 : IF (ASSOCIATED(qs_ot_env%p_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%p_k_k_sym)
1449 9827 : IF (ASSOCIATED(qs_ot_env%buf1_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%buf1_k_k_sym)
1450 9827 : IF (ASSOCIATED(qs_ot_env%buf2_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%buf2_k_k_sym)
1451 9827 : IF (ASSOCIATED(qs_ot_env%buf3_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%buf3_k_k_sym)
1452 9827 : IF (ASSOCIATED(qs_ot_env%buf4_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%buf4_k_k_sym)
1453 9827 : IF (ASSOCIATED(qs_ot_env%buf1_n_k)) CALL dbcsr_release_p(qs_ot_env%buf1_n_k)
1454 9827 : IF (ASSOCIATED(qs_ot_env%buf1_n_k_dp)) CALL dbcsr_release_p(qs_ot_env%buf1_n_k_dp)
1455 :
1456 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
1457 9827 : qs_ot_env%settings%ot_method == "BROY" .OR. &
1458 : qs_ot_env%settings%ot_method == "LBFG") THEN
1459 5708 : CALL dbcsr_deallocate_matrix_set(qs_ot_env%matrix_h_x)
1460 5708 : CALL dbcsr_deallocate_matrix_set(qs_ot_env%matrix_h_e)
1461 5708 : IF (ASSOCIATED(qs_ot_env%matrix_h_x_im)) THEN
1462 227 : CALL dbcsr_deallocate_matrix_set(qs_ot_env%matrix_h_x_im)
1463 : END IF
1464 5708 : IF (ASSOCIATED(qs_ot_env%matrix_h_e_im)) THEN
1465 227 : CALL dbcsr_deallocate_matrix_set(qs_ot_env%matrix_h_e_im)
1466 : END IF
1467 : END IF
1468 9827 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
1469 : qs_ot_env%settings%ot_method == "BROY") THEN
1470 5596 : DEALLOCATE (qs_ot_env%ls_diis)
1471 5596 : DEALLOCATE (qs_ot_env%lss_diis)
1472 5596 : DEALLOCATE (qs_ot_env%c_diis)
1473 5596 : DEALLOCATE (qs_ot_env%c_broy)
1474 5596 : DEALLOCATE (qs_ot_env%energy_h)
1475 5596 : DEALLOCATE (qs_ot_env%ipivot)
1476 : END IF
1477 9827 : IF (qs_ot_env%settings%ot_method == "LBFG") THEN
1478 0 : DEALLOCATE (qs_ot_env%lbfgs_rho, qs_ot_env%lbfgs_yy, &
1479 112 : qs_ot_env%lbfgs_sy_rotation, qs_ot_env%lbfgs_yy_rotation)
1480 : END IF
1481 :
1482 9827 : IF (qs_ot_env%settings%do_rotation) THEN
1483 :
1484 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_u)) CALL dbcsr_release_p(qs_ot_env%rot_mat_u)
1485 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_u_im)) CALL dbcsr_release_p(qs_ot_env%rot_mat_u_im)
1486 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_x)) CALL dbcsr_release_p(qs_ot_env%rot_mat_x)
1487 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_x_im)) CALL dbcsr_release_p(qs_ot_env%rot_mat_x_im)
1488 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_dedu)) CALL dbcsr_release_p(qs_ot_env%rot_mat_dedu)
1489 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_dedu_im)) CALL dbcsr_release_p(qs_ot_env%rot_mat_dedu_im)
1490 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_chc)) CALL dbcsr_release_p(qs_ot_env%rot_mat_chc)
1491 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_chc_im)) CALL dbcsr_release_p(qs_ot_env%rot_mat_chc_im)
1492 :
1493 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
1494 386 : qs_ot_env%settings%ot_method == "BROY" .OR. &
1495 : qs_ot_env%settings%ot_method == "LBFG") THEN
1496 134 : CALL dbcsr_deallocate_matrix_set(qs_ot_env%rot_mat_h_x)
1497 134 : CALL dbcsr_deallocate_matrix_set(qs_ot_env%rot_mat_h_e)
1498 134 : IF (ASSOCIATED(qs_ot_env%rot_mat_h_x_im)) THEN
1499 100 : CALL dbcsr_deallocate_matrix_set(qs_ot_env%rot_mat_h_x_im)
1500 : END IF
1501 134 : IF (ASSOCIATED(qs_ot_env%rot_mat_h_e_im)) THEN
1502 100 : CALL dbcsr_deallocate_matrix_set(qs_ot_env%rot_mat_h_e_im)
1503 : END IF
1504 : END IF
1505 :
1506 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_evals)) DEALLOCATE (qs_ot_env%rot_mat_evals)
1507 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_evec_re)) CALL dbcsr_release_p(qs_ot_env%rot_mat_evec_re)
1508 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_evec_im)) CALL dbcsr_release_p(qs_ot_env%rot_mat_evec_im)
1509 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_gx)) CALL dbcsr_release_p(qs_ot_env%rot_mat_gx)
1510 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_gx_im)) CALL dbcsr_release_p(qs_ot_env%rot_mat_gx_im)
1511 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_response_gx)) THEN
1512 112 : CALL dbcsr_release_p(qs_ot_env%rot_mat_response_gx)
1513 : END IF
1514 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_response_gx_im)) THEN
1515 112 : CALL dbcsr_release_p(qs_ot_env%rot_mat_response_gx_im)
1516 : END IF
1517 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_mermin_g0)) THEN
1518 64 : CALL dbcsr_release_p(qs_ot_env%rot_mat_mermin_g0)
1519 : END IF
1520 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_mermin_g0_im)) THEN
1521 64 : CALL dbcsr_release_p(qs_ot_env%rot_mat_mermin_g0_im)
1522 : END IF
1523 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_gx_old)) CALL dbcsr_release_p(qs_ot_env%rot_mat_gx_old)
1524 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_gx_old_im)) CALL dbcsr_release_p(qs_ot_env%rot_mat_gx_old_im)
1525 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_dx)) CALL dbcsr_release_p(qs_ot_env%rot_mat_dx)
1526 386 : IF (ASSOCIATED(qs_ot_env%rot_mat_dx_im)) CALL dbcsr_release_p(qs_ot_env%rot_mat_dx_im)
1527 : END IF
1528 :
1529 9827 : IF (qs_ot_env%settings%do_ener) THEN
1530 114 : IF (ASSOCIATED(qs_ot_env%ener_x)) DEALLOCATE (qs_ot_env%ener_x)
1531 114 : IF (ASSOCIATED(qs_ot_env%ener_rayleigh)) DEALLOCATE (qs_ot_env%ener_rayleigh)
1532 114 : IF (ASSOCIATED(qs_ot_env%ener_gx)) DEALLOCATE (qs_ot_env%ener_gx)
1533 114 : IF (ASSOCIATED(qs_ot_env%ener_preconditioned_gx)) THEN
1534 112 : DEALLOCATE (qs_ot_env%ener_preconditioned_gx)
1535 : END IF
1536 114 : IF (ASSOCIATED(qs_ot_env%ener_response_gx)) DEALLOCATE (qs_ot_env%ener_response_gx)
1537 : IF (qs_ot_env%settings%ot_method == "DIIS" .OR. &
1538 114 : qs_ot_env%settings%ot_method == "BROY" .OR. &
1539 : qs_ot_env%settings%ot_method == "LBFG") THEN
1540 54 : IF (ASSOCIATED(qs_ot_env%ener_h_x)) DEALLOCATE (qs_ot_env%ener_h_x)
1541 54 : IF (ASSOCIATED(qs_ot_env%ener_h_e)) DEALLOCATE (qs_ot_env%ener_h_e)
1542 : END IF
1543 114 : IF (qs_ot_env%use_dx) THEN
1544 114 : IF (ASSOCIATED(qs_ot_env%ener_dx)) DEALLOCATE (qs_ot_env%ener_dx)
1545 : END IF
1546 114 : IF (qs_ot_env%use_gx_old) THEN
1547 66 : IF (ASSOCIATED(qs_ot_env%ener_gx_old)) DEALLOCATE (qs_ot_env%ener_gx_old)
1548 : END IF
1549 : END IF
1550 :
1551 9827 : qs_ot_env%state_allocated = .FALSE.
1552 9827 : qs_ot_env%has_complex_kpoint_state = .FALSE.
1553 :
1554 : END SUBROUTINE qs_ot_destroy
1555 :
1556 : ! **************************************************************************************************
1557 : !> \brief ...
1558 : !> \param settings ...
1559 : !> \param ot_section ...
1560 : !> \param output_unit ...
1561 : !> \param complex_kpoints whether defaults are for complex k-point OT
1562 : !> \param eigensolver whether OT minimizes a fixed-H orbital subspace for DIAGONALIZATION
1563 : ! **************************************************************************************************
1564 45426 : SUBROUTINE ot_readwrite_input(settings, ot_section, output_unit, complex_kpoints, eigensolver)
1565 : TYPE(qs_ot_settings_type) :: settings
1566 : TYPE(section_vals_type), POINTER :: ot_section
1567 : INTEGER, INTENT(IN) :: output_unit
1568 : LOGICAL, INTENT(IN), OPTIONAL :: complex_kpoints, eigensolver
1569 :
1570 : CHARACTER(len=*), PARAMETER :: routineN = 'ot_readwrite_input'
1571 :
1572 : INTEGER :: handle, ls_method, ot_algorithm, &
1573 : ot_method, ot_ortho_irac
1574 : LOGICAL :: energy_gap_explicit, &
1575 : use_complex_kpoints, use_eigensolver
1576 :
1577 7571 : CALL timeset(routineN, handle)
1578 :
1579 7571 : use_complex_kpoints = .FALSE.
1580 7571 : IF (PRESENT(complex_kpoints)) use_complex_kpoints = complex_kpoints
1581 7571 : use_eigensolver = .FALSE.
1582 7571 : IF (PRESENT(eigensolver)) use_eigensolver = eigensolver
1583 :
1584 : ! choose algorithm
1585 7571 : CALL section_vals_val_get(ot_section, "ALGORITHM", i_val=ot_algorithm)
1586 6911 : SELECT CASE (ot_algorithm)
1587 : CASE (ot_algo_taylor_or_diag)
1588 6911 : settings%ot_algorithm = "TOD"
1589 : CASE (ot_algo_irac)
1590 660 : CALL cite_reference(Weber2008)
1591 660 : settings%ot_algorithm = "REF"
1592 : CASE DEFAULT
1593 7571 : CPABORT("Value unknown")
1594 : END SELECT
1595 :
1596 : ! irac input
1597 7571 : CALL section_vals_val_get(ot_section, "IRAC_DEGREE", i_val=settings%irac_degree)
1598 7571 : IF (settings%irac_degree < 2 .OR. settings%irac_degree > 4) THEN
1599 0 : CPABORT("READ OT IRAC_DEGREE: Value unknown")
1600 : END IF
1601 7571 : CALL section_vals_val_get(ot_section, "MAX_IRAC", i_val=settings%max_irac)
1602 7571 : IF (settings%max_irac < 1) THEN
1603 0 : CPABORT("READ OT MAX_IRAC: VALUE MUST BE GREATER THAN ZERO")
1604 : END IF
1605 7571 : CALL section_vals_val_get(ot_section, "EPS_IRAC_FILTER_MATRIX", r_val=settings%eps_irac_filter_matrix)
1606 7571 : CALL section_vals_val_get(ot_section, "EPS_IRAC", r_val=settings%eps_irac)
1607 7571 : IF (settings%eps_irac < 0.0_dp) THEN
1608 0 : CPABORT("READ OT EPS_IRAC: VALUE MUST BE GREATER THAN ZERO")
1609 : END IF
1610 7571 : CALL section_vals_val_get(ot_section, "EPS_IRAC_QUICK_EXIT", r_val=settings%eps_irac_quick_exit)
1611 7571 : IF (settings%eps_irac_quick_exit < 0.0_dp) THEN
1612 0 : CPABORT("READ OT EPS_IRAC_QUICK_EXIT: VALUE MUST BE GREATER THAN ZERO")
1613 : END IF
1614 :
1615 7571 : CALL section_vals_val_get(ot_section, "EPS_IRAC_SWITCH", r_val=settings%eps_irac_switch)
1616 7571 : IF (settings%eps_irac_switch < 0.0_dp) THEN
1617 0 : CPABORT("READ OT EPS_IRAC_SWITCH: VALUE MUST BE GREATER THAN ZERO")
1618 : END IF
1619 :
1620 7571 : CALL section_vals_val_get(ot_section, "ORTHO_IRAC", i_val=ot_ortho_irac)
1621 7487 : SELECT CASE (ot_ortho_irac)
1622 : CASE (ot_chol_irac)
1623 7487 : settings%ortho_irac = "CHOL"
1624 : CASE (ot_poly_irac)
1625 34 : settings%ortho_irac = "POLY"
1626 : CASE (ot_lwdn_irac)
1627 50 : settings%ortho_irac = "LWDN"
1628 : CASE DEFAULT
1629 7571 : CPABORT("READ OT ORTHO_IRAC: Value unknown")
1630 : END SELECT
1631 :
1632 7571 : CALL section_vals_val_get(ot_section, "ON_THE_FLY_LOC", l_val=settings%on_the_fly_loc)
1633 :
1634 7571 : CALL section_vals_val_get(ot_section, "MINIMIZER", i_val=ot_method)
1635 : ! overwrite input if ot_state is set
1636 7571 : IF (settings%ot_state == 1) THEN
1637 6 : ot_method = ot_mini_diis
1638 : END IF
1639 : ! compatibility
1640 14 : SELECT CASE (ot_method)
1641 : CASE (ot_mini_sd)
1642 14 : settings%ot_method = "SD"
1643 : CASE (ot_mini_cg)
1644 2492 : settings%ot_method = "CG"
1645 : CASE (ot_mini_diis)
1646 5007 : settings%ot_method = "DIIS"
1647 5007 : CALL section_vals_val_get(ot_section, "N_HISTORY_VEC", i_val=settings%diis_m)
1648 : CASE (ot_mini_broyden)
1649 16 : CALL section_vals_val_get(ot_section, "N_HISTORY_VEC", i_val=settings%diis_m)
1650 16 : CALL section_vals_val_get(ot_section, "BROYDEN_BETA", r_val=settings%broyden_beta)
1651 16 : CALL section_vals_val_get(ot_section, "BROYDEN_GAMMA", r_val=settings%broyden_gamma)
1652 16 : CALL section_vals_val_get(ot_section, "BROYDEN_SIGMA", r_val=settings%broyden_sigma)
1653 16 : CALL section_vals_val_get(ot_section, "BROYDEN_ETA", r_val=settings%broyden_eta)
1654 16 : CALL section_vals_val_get(ot_section, "BROYDEN_OMEGA", r_val=settings%broyden_omega)
1655 16 : CALL section_vals_val_get(ot_section, "BROYDEN_SIGMA_DECREASE", r_val=settings%broyden_sigma_decrease)
1656 16 : CALL section_vals_val_get(ot_section, "BROYDEN_SIGMA_MIN", r_val=settings%broyden_sigma_min)
1657 16 : CALL section_vals_val_get(ot_section, "BROYDEN_FORGET_HISTORY", l_val=settings%broyden_forget_history)
1658 16 : CALL section_vals_val_get(ot_section, "BROYDEN_ADAPTIVE_SIGMA", l_val=settings%broyden_adaptive_sigma)
1659 16 : CALL section_vals_val_get(ot_section, "BROYDEN_ENABLE_FLIP", l_val=settings%broyden_enable_flip)
1660 16 : settings%ot_method = "BROY"
1661 : CASE (ot_mini_lbfgs)
1662 42 : CALL section_vals_val_get(ot_section, "N_HISTORY_VEC", i_val=settings%diis_m)
1663 : CALL section_vals_val_get(ot_section, "LBFGS_CURVATURE_TOL", &
1664 42 : r_val=settings%lbfgs_curvature_tol)
1665 : CALL section_vals_val_get(ot_section, "LBFGS_DAMPING", &
1666 42 : l_val=settings%lbfgs_damping)
1667 42 : IF (settings%lbfgs_curvature_tol < 0.0_dp .OR. &
1668 : settings%lbfgs_curvature_tol >= 1.0_dp) THEN
1669 0 : CPABORT("READ OT LBFGS_CURVATURE_TOL: Value must be in [0, 1)")
1670 : END IF
1671 42 : settings%ot_method = "LBFG"
1672 : CASE DEFAULT
1673 7571 : CPABORT("READ OTSCF MINIMIZER: Value unknown")
1674 : END SELECT
1675 7571 : CALL section_vals_val_get(ot_section, "SAFER_DIIS", l_val=settings%safer_diis)
1676 7571 : CALL section_vals_val_get(ot_section, "LINESEARCH", i_val=ls_method)
1677 2 : SELECT CASE (ls_method)
1678 : CASE (ls_none)
1679 2 : settings%line_search_method = "NONE"
1680 : CASE (ls_2pnt)
1681 7429 : settings%line_search_method = "2PNT"
1682 : CASE (ls_3pnt)
1683 42 : settings%line_search_method = "3PNT"
1684 : CASE (ls_adapt)
1685 86 : settings%line_search_method = "ADPT"
1686 : CASE (ls_gold)
1687 12 : settings%line_search_method = "GOLD"
1688 12 : CALL section_vals_val_get(ot_section, "GOLD_TARGET", r_val=settings%gold_target)
1689 : CASE DEFAULT
1690 7571 : CPABORT("READ OTSCF LS: Value unknown")
1691 : END SELECT
1692 :
1693 7571 : CALL section_vals_val_get(ot_section, "PRECOND_SOLVER", i_val=settings%precond_solver_type)
1694 15118 : SELECT CASE (settings%precond_solver_type)
1695 : CASE (ot_precond_solver_default)
1696 7547 : settings%precond_solver_name = "DEFAULT"
1697 : CASE (ot_precond_solver_inv_chol)
1698 16 : settings%precond_solver_name = "INVERSE_CHOLESKY"
1699 : CASE (ot_precond_solver_direct)
1700 0 : settings%precond_solver_name = "DIRECT"
1701 : CASE (ot_precond_solver_update)
1702 6 : settings%precond_solver_name = "INVERSE_UPDATE"
1703 : CASE (ot_precond_solver_chebyshev)
1704 2 : settings%precond_solver_name = "CHEBYSHEV"
1705 : CASE DEFAULT
1706 7571 : CPABORT("READ OTSCF SOLVER: Value unknown")
1707 : END SELECT
1708 7571 : CALL section_vals_val_get(ot_section, "CHEBYSHEV_DEGREE", i_val=settings%chebyshev_degree)
1709 7571 : IF (settings%chebyshev_degree < 1 .OR. settings%chebyshev_degree > 100) THEN
1710 0 : CPABORT("READ OT CHEBYSHEV_DEGREE: Value must be between 1 and 100")
1711 : END IF
1712 :
1713 7571 : CALL section_vals_val_get(ot_section, "MAX_SCF_DIIS", i_val=settings%max_scf_diis)
1714 :
1715 : !If these values are negative we will set them "optimal" for a given precondtioner below
1716 7571 : CALL section_vals_val_get(ot_section, "STEPSIZE", r_val=settings%ds_min)
1717 : CALL section_vals_val_get(ot_section, "ENERGY_GAP", r_val=settings%energy_gap, &
1718 7571 : explicit=energy_gap_explicit)
1719 :
1720 7571 : CALL section_vals_val_get(ot_section, "PRECONDITIONER", i_val=settings%preconditioner_type)
1721 8303 : SELECT CASE (settings%preconditioner_type)
1722 : CASE (ot_precond_none)
1723 732 : settings%preconditioner_name = "NONE"
1724 732 : IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
1725 732 : IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.2_dp
1726 : CASE (ot_precond_full_single)
1727 34 : settings%preconditioner_name = "FULL_SINGLE"
1728 34 : IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
1729 34 : IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.2_dp
1730 : CASE (ot_precond_full_single_inverse)
1731 2890 : settings%preconditioner_name = "FULL_SINGLE_INVERSE"
1732 2890 : IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.08_dp
1733 2890 : IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.08_dp
1734 : CASE (ot_precond_full_all)
1735 2572 : settings%preconditioner_name = "FULL_ALL"
1736 2572 : IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
1737 2572 : IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.08_dp
1738 : CASE (ot_precond_full_kinetic)
1739 1283 : settings%preconditioner_name = "FULL_KINETIC"
1740 1283 : IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
1741 1283 : IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.2_dp
1742 : CASE (ot_precond_s_inverse)
1743 60 : settings%preconditioner_name = "FULL_S_INVERSE"
1744 60 : IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
1745 60 : IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.2_dp
1746 : CASE DEFAULT
1747 7571 : CPABORT("READ OTSCF PRECONDITIONER: Value unknown")
1748 : END SELECT
1749 : IF (use_complex_kpoints .AND. &
1750 7571 : settings%preconditioner_type == ot_precond_full_single_inverse .AND. &
1751 2 : .NOT. energy_gap_explicit) settings%energy_gap = 0.20_dp
1752 7571 : CALL section_vals_val_get(ot_section, "CHOLESKY", i_val=settings%cholesky_type)
1753 7571 : CALL section_vals_val_get(ot_section, "EPS_TAYLOR", r_val=settings%eps_taylor)
1754 7571 : CALL section_vals_val_get(ot_section, "MAX_TAYLOR", i_val=settings%max_taylor)
1755 7571 : CALL section_vals_val_get(ot_section, "ROTATION", l_val=settings%do_rotation)
1756 7571 : CALL section_vals_val_get(ot_section, "ENERGIES", l_val=settings%do_ener)
1757 : CALL section_vals_val_get(ot_section, "OCCUPATION_PRECONDITIONER", &
1758 7571 : l_val=settings%occupation_preconditioner)
1759 7571 : CALL section_vals_val_get(ot_section, "NONDIAG_ENERGY", l_val=settings%add_nondiag_energy)
1760 : CALL section_vals_val_get(ot_section, "NONDIAG_ENERGY_STRENGTH", &
1761 7571 : r_val=settings%nondiag_energy_strength)
1762 7571 : IF (settings%ot_method == "LBFG") THEN
1763 42 : IF (settings%ot_algorithm /= "TOD" .AND. settings%ot_algorithm /= "REF") THEN
1764 0 : CPABORT("MINIMIZER LBFGS requires ALGORITHM STRICT or IRAC")
1765 : END IF
1766 : END IF
1767 7571 : IF (settings%do_ener .AND. .NOT. use_complex_kpoints .AND. .NOT. use_eigensolver) THEN
1768 0 : CPABORT("OT%ENERGIES is currently supported only by the complex K-point Mermin path.")
1769 : END IF
1770 7571 : IF (use_eigensolver .AND. &
1771 : (settings%do_rotation .OR. settings%do_ener .OR. &
1772 : settings%occupation_preconditioner .OR. settings%add_nondiag_energy)) THEN
1773 : CALL cp_warn(__LOCATION__, &
1774 : "DIAGONALIZATION%OT ignores ROTATION, ENERGIES, "// &
1775 : "OCCUPATION_PRECONDITIONER, and NONDIAG_ENERGY because occupations "// &
1776 2 : "are assigned after eigenspace minimization.")
1777 2 : settings%do_rotation = .FALSE.
1778 2 : settings%do_ener = .FALSE.
1779 2 : settings%occupation_preconditioner = .FALSE.
1780 2 : settings%add_nondiag_energy = .FALSE.
1781 : END IF
1782 :
1783 : ! write OT output
1784 :
1785 7571 : IF (output_unit > 0) THEN
1786 3902 : WRITE (output_unit, '(/,A)') " ----------------------------------- OT ---------------------------------------"
1787 3902 : IF (settings%do_rotation) THEN
1788 154 : WRITE (output_unit, '(A)') " Allowing for rotations "
1789 : END IF
1790 3902 : IF (settings%do_ener) THEN
1791 41 : WRITE (output_unit, '(A,L2)') " Optimizing orbital energies "
1792 : END IF
1793 7 : SELECT CASE (settings%OT_METHOD)
1794 : CASE ("SD")
1795 7 : WRITE (output_unit, '(A)') " Minimizer : SD : steepest descent"
1796 : CASE ("CG")
1797 1275 : WRITE (output_unit, '(A)') " Minimizer : CG : conjugate gradient"
1798 : CASE ("DIIS")
1799 2591 : WRITE (output_unit, '(A)') " Minimizer : DIIS : direct inversion"
1800 2591 : WRITE (output_unit, '(A)') " in the iterative subspace"
1801 2591 : WRITE (output_unit, '(A,I3,A)') " using ", settings%diis_m, " DIIS vectors"
1802 2591 : IF (settings%safer_diis) THEN
1803 2591 : WRITE (output_unit, '(A,I3,A)') " safer DIIS on"
1804 : ELSE
1805 0 : WRITE (output_unit, '(A,I3,A)') " safer DIIS off"
1806 : END IF
1807 : CASE ("BROY")
1808 8 : WRITE (output_unit, '(A)') " Minimizer : BROYDEN : Broyden "
1809 8 : WRITE (output_unit, '(A,F16.8)') " BETA : ", settings%broyden_beta
1810 8 : WRITE (output_unit, '(A,F16.8)') " GAMMA : ", settings%broyden_gamma
1811 8 : WRITE (output_unit, '(A,F16.8)') " SIGMA : ", settings%broyden_sigma
1812 8 : WRITE (output_unit, '(A,I3,A)') " using : - ", &
1813 16 : settings%diis_m, " BROYDEN vectors"
1814 : CASE ("LBFG")
1815 21 : WRITE (output_unit, '(A)') " Minimizer : LBFGS : limited-memory BFGS"
1816 21 : WRITE (output_unit, '(A,I3,A)') " using : ", &
1817 42 : settings%diis_m, " secant pairs"
1818 21 : WRITE (output_unit, '(A,ES12.4)') " curvature tolerance : ", &
1819 42 : settings%lbfgs_curvature_tol
1820 21 : WRITE (output_unit, '(A,L7)') " damp weak curvature : ", &
1821 42 : settings%lbfgs_damping
1822 : CASE DEFAULT
1823 3902 : WRITE (output_unit, '(3A)') " Minimizer : ", settings%OT_METHOD, " : UNKNOWN"
1824 : END SELECT
1825 17 : SELECT CASE (settings%preconditioner_name)
1826 : CASE ("FULL_SINGLE")
1827 17 : WRITE (output_unit, '(A)') " Preconditioner : FULL_SINGLE : diagonalization based"
1828 : CASE ("FULL_SINGLE_INVERSE")
1829 1512 : WRITE (output_unit, '(A,/,A)') " Preconditioner : FULL_SINGLE_INVERSE : inversion of ", &
1830 3024 : " H + eS - 2*(Sc)(c^T*H*c+const)(Sc)^T"
1831 : CASE ("FULL_ALL")
1832 1316 : WRITE (output_unit, '(A)') " Preconditioner : FULL_ALL : diagonalization, state selective"
1833 : CASE ("FULL_KINETIC")
1834 661 : WRITE (output_unit, '(A)') " Preconditioner : FULL_KINETIC : inversion of T + eS"
1835 : CASE ("FULL_S_INVERSE")
1836 30 : WRITE (output_unit, '(A)') " Preconditioner : FULL_S_INVERSE : cholesky inversion of S"
1837 : CASE ("SPARSE_DIAG")
1838 : WRITE (output_unit, '(A)') &
1839 0 : " Preconditioner : SPARSE_DIAG : diagonal atomic block diagonalization"
1840 : CASE ("SPARSE_KINETIC")
1841 0 : WRITE (output_unit, '(A)') " Preconditioner : SPARSE_KINETIC : sparse linear solver for T + eS"
1842 : CASE ("NONE")
1843 366 : WRITE (output_unit, '(A)') " Preconditioner : NONE"
1844 : CASE DEFAULT
1845 3902 : WRITE (output_unit, '(3A)') " Preconditioner : ", settings%preconditioner_name, " : UNKNOWN"
1846 : END SELECT
1847 :
1848 3902 : WRITE (output_unit, '(A)') " Precond_solver : "//TRIM(settings%precond_solver_name)
1849 3902 : IF (settings%precond_solver_type == ot_precond_solver_chebyshev) THEN
1850 1 : WRITE (output_unit, '(A,I14)') " Chebyshev degree:", settings%chebyshev_degree
1851 : END IF
1852 :
1853 3902 : IF (settings%OT_METHOD == "SD" .OR. settings%OT_METHOD == "CG" .OR. &
1854 : settings%OT_METHOD == "LBFG") THEN
1855 1271 : SELECT CASE (settings%line_search_method)
1856 : CASE ("2PNT")
1857 1271 : WRITE (output_unit, '(A)') " Line search : 2PNT : 2 energies, one gradient"
1858 : CASE ("3PNT")
1859 19 : WRITE (output_unit, '(A)') " Line search : 3PNT : 3 energies"
1860 : CASE ("GOLD")
1861 6 : WRITE (output_unit, '(A)') " Line search : GOLD : bracketing and golden section search"
1862 6 : WRITE (output_unit, '(A,F14.8)') " target rel accuracy : ", settings%gold_target
1863 : CASE ("NONE")
1864 1 : WRITE (output_unit, '(A)') " Line search : NONE"
1865 : CASE DEFAULT
1866 3902 : WRITE (output_unit, '(3A)') " Line search : ", settings%line_search_method, " : UNKNOWN"
1867 : END SELECT
1868 : END IF
1869 3902 : WRITE (output_unit, '(A,F14.8,T49,A,F14.8)') " stepsize :", settings%ds_min, &
1870 7804 : " energy_gap :", settings%energy_gap
1871 3902 : IF (settings%ot_algorithm == 'TOD') THEN
1872 3554 : WRITE (output_unit, '(A,E14.5,T49,A,I14)') " eps_taylor :", settings%eps_taylor, &
1873 7108 : " max_taylor :", settings%max_taylor
1874 : END IF
1875 3902 : IF (settings%ot_algorithm == 'REF') THEN
1876 348 : WRITE (output_unit, '(A,1X,A,T49,A,I14)') " ortho_irac :", settings%ortho_irac, &
1877 696 : " irac_degree :", settings%irac_degree
1878 348 : WRITE (output_unit, '(A,I14,T49,A,E14.5)') " max_irac :", settings%max_irac, &
1879 696 : " eps_irac :", settings%eps_irac
1880 348 : WRITE (output_unit, '(A,E14.5,T49,A,E10.3)') " eps_irac_switch:", settings%eps_irac_switch, &
1881 696 : " eps_irac_quick_exit:", settings%eps_irac_quick_exit
1882 348 : WRITE (output_unit, '(A,L2)') " on_the_fly_loc :", settings%on_the_fly_loc
1883 : END IF
1884 3902 : WRITE (output_unit, '(A)') " ----------------------------------- OT ---------------------------------------"
1885 : WRITE (UNIT=output_unit, &
1886 : FMT="(/,T3,A,T12,A,T31,A,T39,A,T59,A,T75,A,/,T3,A)") &
1887 3902 : "Step", "Update method", "Time", "Convergence", "Total energy", "Change", &
1888 7804 : REPEAT("-", 78)
1889 : END IF
1890 :
1891 7571 : CALL timestop(handle)
1892 :
1893 7571 : END SUBROUTINE ot_readwrite_input
1894 :
1895 0 : END MODULE qs_ot_types
|
