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