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 localize wavefunctions
10 : !> linear response scf
11 : !> \par History
12 : !> created 07-2005 [MI]
13 : !> \author MI
14 : ! **************************************************************************************************
15 : MODULE qs_linres_methods
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
18 : dbcsr_p_type,&
19 : dbcsr_set,&
20 : dbcsr_type
21 : USE cp_dbcsr_contrib, ONLY: dbcsr_checksum
22 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
23 : USE cp_external_control, ONLY: external_control
24 : USE cp_files, ONLY: close_file,&
25 : open_file
26 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
27 : cp_fm_trace
28 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
29 : cp_fm_struct_release,&
30 : cp_fm_struct_type
31 : USE cp_fm_types, ONLY: cp_fm_create,&
32 : cp_fm_get_info,&
33 : cp_fm_get_submatrix,&
34 : cp_fm_release,&
35 : cp_fm_set_submatrix,&
36 : cp_fm_to_fm,&
37 : cp_fm_type
38 : USE cp_log_handling, ONLY: cp_get_default_logger,&
39 : cp_logger_type,&
40 : cp_to_string
41 : USE cp_output_handling, ONLY: cp_p_file,&
42 : cp_print_key_finished_output,&
43 : cp_print_key_generate_filename,&
44 : cp_print_key_should_output,&
45 : cp_print_key_unit_nr
46 : USE input_constants, ONLY: do_loc_none,&
47 : op_loc_berry,&
48 : ot_precond_none,&
49 : ot_precond_solver_default,&
50 : state_loc_all
51 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
52 : section_vals_type,&
53 : section_vals_val_get
54 : USE kinds, ONLY: default_path_length,&
55 : default_string_length,&
56 : dp
57 : USE machine, ONLY: m_flush,&
58 : m_walltime
59 : USE message_passing, ONLY: mp_para_env_type
60 : USE parallel_gemm_api, ONLY: parallel_gemm
61 : USE preconditioner, ONLY: apply_preconditioner,&
62 : make_preconditioner
63 : USE qs_2nd_kernel_ao, ONLY: build_dm_response
64 : USE qs_environment_types, ONLY: get_qs_env,&
65 : qs_environment_type
66 : USE qs_gapw_densities, ONLY: prepare_gapw_den
67 : USE qs_linres_kernel, ONLY: apply_op_2
68 : USE qs_linres_types, ONLY: linres_control_type
69 : USE qs_loc_main, ONLY: qs_loc_driver
70 : USE qs_loc_types, ONLY: get_qs_loc_env,&
71 : localized_wfn_control_type,&
72 : qs_loc_env_create,&
73 : qs_loc_env_type
74 : USE qs_loc_utils, ONLY: loc_write_restart,&
75 : qs_loc_control_init,&
76 : qs_loc_init
77 : USE qs_mo_types, ONLY: get_mo_set,&
78 : mo_set_type
79 : USE qs_p_env_methods, ONLY: p_env_check_i_alloc,&
80 : p_env_update_rho
81 : USE qs_p_env_types, ONLY: qs_p_env_type
82 : USE qs_rho_methods, ONLY: qs_rho_update_rho
83 : USE qs_rho_types, ONLY: qs_rho_type
84 : USE string_utilities, ONLY: xstring
85 : #include "./base/base_uses.f90"
86 :
87 : IMPLICIT NONE
88 :
89 : PRIVATE
90 :
91 : ! *** Public subroutines ***
92 : PUBLIC :: linres_localize, linres_solver
93 : PUBLIC :: linres_write_restart, linres_read_restart
94 : PUBLIC :: build_dm_response
95 :
96 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_methods'
97 :
98 : ! **************************************************************************************************
99 :
100 : CONTAINS
101 :
102 : ! **************************************************************************************************
103 : !> \brief Find the centers and spreads of the wfn,
104 : !> if required apply a localization algorithm
105 : !> \param qs_env ...
106 : !> \param linres_control ...
107 : !> \param nspins ...
108 : !> \param centers_only ...
109 : !> \par History
110 : !> 07.2005 created [MI]
111 : !> \author MI
112 : ! **************************************************************************************************
113 190 : SUBROUTINE linres_localize(qs_env, linres_control, nspins, centers_only)
114 :
115 : TYPE(qs_environment_type), POINTER :: qs_env
116 : TYPE(linres_control_type), POINTER :: linres_control
117 : INTEGER, INTENT(IN) :: nspins
118 : LOGICAL, INTENT(IN), OPTIONAL :: centers_only
119 :
120 : INTEGER :: iounit, ispin, istate, nmoloc(2)
121 : LOGICAL :: my_centers_only
122 190 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mos_localized
123 : TYPE(cp_fm_type), POINTER :: mo_coeff
124 : TYPE(cp_logger_type), POINTER :: logger
125 : TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
126 190 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
127 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
128 : TYPE(section_vals_type), POINTER :: loc_print_section, loc_section, &
129 : lr_section
130 :
131 190 : NULLIFY (logger, lr_section, loc_section, loc_print_section, localized_wfn_control)
132 380 : logger => cp_get_default_logger()
133 190 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
134 190 : loc_section => section_vals_get_subs_vals(lr_section, "LOCALIZE")
135 190 : loc_print_section => section_vals_get_subs_vals(lr_section, "LOCALIZE%PRINT")
136 : iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
137 190 : extension=".linresLog")
138 190 : my_centers_only = .FALSE.
139 190 : IF (PRESENT(centers_only)) my_centers_only = centers_only
140 :
141 190 : NULLIFY (mos, mo_coeff, qs_loc_env)
142 190 : CALL get_qs_env(qs_env=qs_env, mos=mos)
143 1520 : ALLOCATE (qs_loc_env)
144 190 : CALL qs_loc_env_create(qs_loc_env)
145 190 : linres_control%qs_loc_env => qs_loc_env
146 190 : CALL qs_loc_control_init(qs_loc_env, loc_section, do_homo=.TRUE.)
147 190 : CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
148 :
149 838 : ALLOCATE (mos_localized(nspins))
150 458 : DO ispin = 1, nspins
151 268 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
152 268 : CALL cp_fm_create(mos_localized(ispin), mo_coeff%matrix_struct)
153 458 : CALL cp_fm_to_fm(mo_coeff, mos_localized(ispin))
154 : END DO
155 :
156 : nmoloc(1:2) = 0
157 190 : IF (my_centers_only) THEN
158 0 : localized_wfn_control%set_of_states = state_loc_all
159 0 : localized_wfn_control%localization_method = do_loc_none
160 0 : localized_wfn_control%operator_type = op_loc_berry
161 : END IF
162 :
163 : CALL qs_loc_init(qs_env, qs_loc_env, loc_section, mos_localized=mos_localized, &
164 190 : do_homo=.TRUE.)
165 :
166 : ! The orbital centers are stored in linres_control%localized_wfn_control
167 458 : DO ispin = 1, nspins
168 : CALL qs_loc_driver(qs_env, qs_loc_env, loc_print_section, myspin=ispin, &
169 268 : ext_mo_coeff=mos_localized(ispin))
170 268 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
171 458 : CALL cp_fm_to_fm(mos_localized(ispin), mo_coeff)
172 : END DO
173 :
174 : CALL loc_write_restart(qs_loc_env, loc_print_section, mos, &
175 190 : mos_localized, do_homo=.TRUE.)
176 190 : CALL cp_fm_release(mos_localized)
177 :
178 : ! Write Centers and Spreads on std out
179 190 : IF (iounit > 0) THEN
180 229 : DO ispin = 1, nspins
181 : WRITE (iounit, "(/,T2,A,I2)") &
182 134 : "WANNIER CENTERS for spin ", ispin
183 : WRITE (iounit, "(/,T18,A,3X,A)") &
184 134 : "--------------- Centers --------------- ", &
185 268 : "--- Spreads ---"
186 1102 : DO istate = 1, SIZE(localized_wfn_control%centers_set(ispin)%array, 2)
187 : WRITE (iounit, "(T5,A6,I6,2X,3f12.6,5X,f12.6)") &
188 3492 : 'state ', istate, localized_wfn_control%centers_set(ispin)%array(1:3, istate), &
189 1880 : localized_wfn_control%centers_set(ispin)%array(4, istate)
190 : END DO
191 : END DO
192 95 : CALL m_flush(iounit)
193 : END IF
194 :
195 570 : END SUBROUTINE linres_localize
196 :
197 : ! **************************************************************************************************
198 : !> \brief scf loop to optimize the first order wavefunctions (psi1)
199 : !> given a perturbation as an operator applied to the ground
200 : !> state orbitals (h1_psi0)
201 : !> psi1 is defined wrt psi0_order (can be a subset of the occupied space)
202 : !> The reference ground state is defined through qs_env (density and ground state MOs)
203 : !> psi1 is orthogonal to the occupied orbitals in the ground state MO set (qs_env%mos)
204 : !> \param p_env ...
205 : !> \param qs_env ...
206 : !> \param psi1 ...
207 : !> \param h1_psi0 ...
208 : !> \param psi0_order ...
209 : !> \param iounit ...
210 : !> \param should_stop ...
211 : !> \param silent ...
212 : !> \par History
213 : !> 07.2005 created [MI]
214 : !> \author MI
215 : ! **************************************************************************************************
216 5100 : SUBROUTINE linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, &
217 : should_stop, silent)
218 : TYPE(qs_p_env_type) :: p_env
219 : TYPE(qs_environment_type), POINTER :: qs_env
220 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: psi1, h1_psi0, psi0_order
221 : INTEGER, INTENT(IN) :: iounit
222 : LOGICAL, INTENT(OUT) :: should_stop
223 : LOGICAL, INTENT(IN), OPTIONAL :: silent
224 :
225 : CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_solver'
226 :
227 : INTEGER :: handle, ispin, iter, maxnmo, nao, ncol, &
228 : nmo, nocc, nspins
229 : LOGICAL :: my_silent, restart
230 : REAL(dp) :: alpha, beta, norm_res, t1, t2
231 5100 : REAL(dp), DIMENSION(:), POINTER :: tr_pAp, tr_rz0, tr_rz00, tr_rz1
232 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
233 : TYPE(cp_fm_type) :: buf
234 5100 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: Ap, chc, mo_coeff_array, p, r, z
235 5100 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: Sc
236 : TYPE(cp_fm_type), POINTER :: mo_coeff
237 5100 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_t
238 : TYPE(dbcsr_type), POINTER :: matrix_x
239 : TYPE(dft_control_type), POINTER :: dft_control
240 : TYPE(linres_control_type), POINTER :: linres_control
241 5100 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
242 : TYPE(mp_para_env_type), POINTER :: para_env
243 :
244 5100 : CALL timeset(routineN, handle)
245 :
246 5100 : my_silent = .FALSE.
247 5100 : IF (PRESENT(silent)) my_silent = silent
248 :
249 5100 : NULLIFY (dft_control, linres_control, matrix_s, matrix_t, matrix_ks, para_env)
250 5100 : NULLIFY (mos, tmp_fm_struct, mo_coeff)
251 :
252 5100 : t1 = m_walltime()
253 :
254 : CALL get_qs_env(qs_env=qs_env, &
255 : matrix_ks=matrix_ks, &
256 : matrix_s=matrix_s, &
257 : kinetic=matrix_t, &
258 : dft_control=dft_control, &
259 : linres_control=linres_control, &
260 : para_env=para_env, &
261 5100 : mos=mos)
262 :
263 5100 : nspins = dft_control%nspins
264 5100 : CALL cp_fm_get_info(psi1(1), nrow_global=nao)
265 5100 : maxnmo = 0
266 12264 : DO ispin = 1, nspins
267 7164 : CALL cp_fm_get_info(psi1(ispin), ncol_global=ncol)
268 12264 : maxnmo = MAX(maxnmo, ncol)
269 : END DO
270 : !
271 5100 : CALL check_p_env_init(p_env, linres_control, nspins)
272 : !
273 : ! allocate the vectors
274 : ALLOCATE (tr_pAp(nspins), tr_rz0(nspins), tr_rz00(nspins), tr_rz1(nspins), &
275 84756 : r(nspins), p(nspins), z(nspins), Ap(nspins))
276 : !
277 12264 : DO ispin = 1, nspins
278 7164 : CALL cp_fm_create(r(ispin), psi1(ispin)%matrix_struct)
279 7164 : CALL cp_fm_create(p(ispin), psi1(ispin)%matrix_struct)
280 7164 : CALL cp_fm_create(z(ispin), psi1(ispin)%matrix_struct)
281 12264 : CALL cp_fm_create(Ap(ispin), psi1(ispin)%matrix_struct)
282 : END DO
283 : !
284 : ! build C0 occupied vectors and S*C0 matrix
285 34728 : ALLOCATE (Sc(nspins), mo_coeff_array(nspins))
286 12264 : DO ispin = 1, nspins
287 7164 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, homo=nocc)
288 7164 : NULLIFY (tmp_fm_struct)
289 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
290 : ncol_global=nocc, para_env=para_env, &
291 7164 : context=mo_coeff%matrix_struct%context)
292 7164 : CALL cp_fm_create(mo_coeff_array(ispin), tmp_fm_struct)
293 7164 : CALL cp_fm_to_fm(mo_coeff, mo_coeff_array(ispin), nocc)
294 7164 : CALL cp_fm_create(Sc(ispin), tmp_fm_struct)
295 19428 : CALL cp_fm_struct_release(tmp_fm_struct)
296 : END DO
297 : !
298 : ! Allocate C0_order'*H*C0_order
299 22464 : ALLOCATE (chc(nspins))
300 12264 : DO ispin = 1, nspins
301 7164 : CALL cp_fm_get_info(psi1(ispin), ncol_global=nmo)
302 7164 : NULLIFY (tmp_fm_struct)
303 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
304 : ncol_global=nmo, para_env=para_env, &
305 7164 : context=mo_coeff%matrix_struct%context)
306 7164 : CALL cp_fm_create(chc(ispin), tmp_fm_struct, set_zero=.TRUE.)
307 19428 : CALL cp_fm_struct_release(tmp_fm_struct)
308 : END DO
309 : !
310 12264 : DO ispin = 1, nspins
311 : !
312 : ! C0_order' * H * C0_order
313 : ASSOCIATE (mo_coeff => psi0_order(ispin))
314 7164 : CALL cp_fm_create(buf, mo_coeff%matrix_struct)
315 7164 : CALL cp_fm_get_info(mo_coeff, ncol_global=ncol)
316 7164 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, mo_coeff, buf, ncol)
317 7164 : CALL parallel_gemm('T', 'N', ncol, ncol, nao, -1.0_dp, mo_coeff, buf, 0.0_dp, chc(ispin))
318 7164 : CALL cp_fm_release(buf)
319 : END ASSOCIATE
320 : !
321 : ! S * C0
322 7164 : CALL cp_fm_get_info(mo_coeff_array(ispin), ncol_global=ncol)
323 19428 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff_array(ispin), Sc(ispin), ncol)
324 : END DO
325 : !
326 : ! header
327 5100 : IF (iounit > 0 .AND. .NOT. my_silent) THEN
328 : WRITE (iounit, "(/,T3,A,T16,A,T25,A,T38,A,T52,A,T72,A,/,T3,A)") &
329 2537 : "Iteration", "Method", "Restart", "Stepsize", "Convergence", "Time", &
330 5074 : REPEAT("-", 78)
331 : END IF
332 : !
333 : ! orthogonalize x with respect to the psi0
334 5100 : CALL preortho(psi1, mo_coeff_array, Sc)
335 : !
336 : ! build the preconditioner
337 5100 : IF (linres_control%preconditioner_type /= ot_precond_none) THEN
338 5100 : IF (p_env%new_preconditioner) THEN
339 2796 : DO ispin = 1, nspins
340 2796 : IF (ASSOCIATED(matrix_t)) THEN
341 : CALL make_preconditioner(p_env%preconditioner(ispin), &
342 : linres_control%preconditioner_type, ot_precond_solver_default, &
343 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, matrix_t(1)%matrix, &
344 1478 : mos(ispin), linres_control%energy_gap)
345 : ELSE
346 24 : NULLIFY (matrix_x)
347 : CALL make_preconditioner(p_env%preconditioner(ispin), &
348 : linres_control%preconditioner_type, ot_precond_solver_default, &
349 : matrix_ks(ispin)%matrix, matrix_s(1)%matrix, matrix_x, &
350 24 : mos(ispin), linres_control%energy_gap)
351 : END IF
352 : END DO
353 1294 : p_env%new_preconditioner = .FALSE.
354 : END IF
355 : END IF
356 : !
357 : ! initialization of the linear solver
358 : !
359 : ! A * x0
360 5100 : CALL apply_op(qs_env, p_env, psi0_order, psi1, Ap, chc)
361 : !
362 : !
363 : ! r_0 = b - Ax0
364 12264 : DO ispin = 1, nspins
365 7164 : CALL cp_fm_to_fm(h1_psi0(ispin), r(ispin))
366 12264 : CALL cp_fm_scale_and_add(-1.0_dp, r(ispin), -1.0_dp, Ap(ispin))
367 : END DO
368 : !
369 : ! proj r
370 5100 : CALL postortho(r, mo_coeff_array, Sc)
371 : !
372 : ! preconditioner
373 5100 : linres_control%flag = ""
374 5100 : IF (linres_control%preconditioner_type == ot_precond_none) THEN
375 : !
376 : ! z_0 = r_0
377 0 : DO ispin = 1, nspins
378 0 : CALL cp_fm_to_fm(r(ispin), z(ispin))
379 : END DO
380 0 : linres_control%flag = "CG"
381 : ELSE
382 : !
383 : ! z_0 = M * r_0
384 12264 : DO ispin = 1, nspins
385 : CALL apply_preconditioner(p_env%preconditioner(ispin), r(ispin), &
386 12264 : z(ispin))
387 : END DO
388 5100 : linres_control%flag = "PCG"
389 : END IF
390 : !
391 12264 : DO ispin = 1, nspins
392 : !
393 : ! p_0 = z_0
394 7164 : CALL cp_fm_to_fm(z(ispin), p(ispin))
395 : !
396 : ! trace(r_0 * z_0)
397 12264 : CALL cp_fm_trace(r(ispin), z(ispin), tr_rz0(ispin))
398 : END DO
399 12264 : IF (SUM(tr_rz0) < 0.0_dp) CPABORT("tr(r_j*z_j) < 0")
400 12264 : norm_res = ABS(SUM(tr_rz0))/SQRT(REAL(nspins*nao*maxnmo, dp))
401 : !
402 5100 : alpha = 0.0_dp
403 5100 : restart = .FALSE.
404 5100 : should_stop = .FALSE.
405 33746 : iteration: DO iter = 1, linres_control%max_iter
406 : !
407 : ! check convergence
408 32426 : linres_control%converged = .FALSE.
409 32426 : IF (norm_res < linres_control%eps) THEN
410 3780 : linres_control%converged = .TRUE.
411 : END IF
412 : !
413 32426 : t2 = m_walltime()
414 : IF (iter == 1 .OR. MOD(iter, 1) == 0 .OR. linres_control%converged &
415 : .OR. restart .OR. should_stop) THEN
416 32426 : IF (iounit > 0 .AND. .NOT. my_silent) THEN
417 : WRITE (iounit, "(T5,I5,T18,A3,T28,L1,T38,1E8.2,T48,F16.10,T68,F8.2)") &
418 16111 : iter, linres_control%flag, restart, alpha, norm_res, t2 - t1
419 16111 : CALL m_flush(iounit)
420 : END IF
421 : END IF
422 : !
423 32426 : IF (linres_control%converged) THEN
424 3780 : IF (iounit > 0) THEN
425 1882 : WRITE (iounit, "(T3,A,I4,A,T73,F8.2)") "The linear solver converged in ", &
426 3764 : iter, " iterations.", t2 - t1
427 1882 : CALL m_flush(iounit)
428 : END IF
429 : EXIT iteration
430 28646 : ELSE IF (should_stop) THEN
431 0 : IF (iounit > 0) THEN
432 0 : WRITE (iounit, "(T3,A,I4,A)") "The linear solver did NOT converge! External stop"
433 0 : CALL m_flush(iounit)
434 : END IF
435 : EXIT iteration
436 : END IF
437 : !
438 : ! Max number of iteration reached
439 28646 : IF (iter == linres_control%max_iter) THEN
440 1320 : IF (iounit > 0) THEN
441 : CALL cp_warn(__LOCATION__, &
442 660 : "The linear solver didn't converge! Maximum number of iterations reached.")
443 660 : CALL m_flush(iounit)
444 : END IF
445 1320 : linres_control%converged = .FALSE.
446 : END IF
447 : !
448 : ! Apply the operators that do not depend on the perturbation
449 28646 : CALL apply_op(qs_env, p_env, psi0_order, p, Ap, chc)
450 : !
451 : ! proj Ap onto the virtual subspace
452 28646 : CALL postortho(Ap, mo_coeff_array, Sc)
453 : !
454 70224 : DO ispin = 1, nspins
455 : !
456 : ! tr(Ap_j*p_j)
457 70224 : CALL cp_fm_trace(Ap(ispin), p(ispin), tr_pAp(ispin))
458 : END DO
459 : !
460 : ! alpha = tr(r_j*z_j) / tr(Ap_j*p_j)
461 70224 : IF (SUM(tr_pAp) < 1.0e-10_dp) THEN
462 176 : alpha = 1.0_dp
463 : ELSE
464 111274 : alpha = SUM(tr_rz0)/SUM(tr_pAp)
465 : END IF
466 70224 : DO ispin = 1, nspins
467 : !
468 : ! x_j+1 = x_j + alpha * p_j
469 70224 : CALL cp_fm_scale_and_add(1.0_dp, psi1(ispin), alpha, p(ispin))
470 : END DO
471 : !
472 : ! need to recompute the residue
473 28646 : restart = .FALSE.
474 28646 : IF (MOD(iter, linres_control%restart_every) == 0) THEN
475 : !
476 : ! r_j+1 = b - A * x_j+1
477 206 : CALL apply_op(qs_env, p_env, psi0_order, psi1, Ap, chc)
478 : !
479 446 : DO ispin = 1, nspins
480 240 : CALL cp_fm_to_fm(h1_psi0(ispin), r(ispin))
481 446 : CALL cp_fm_scale_and_add(-1.0_dp, r(ispin), -1.0_dp, Ap(ispin))
482 : END DO
483 206 : CALL postortho(r, mo_coeff_array, Sc)
484 : !
485 206 : restart = .TRUE.
486 : ELSE
487 : ! proj Ap onto the virtual subspace
488 28440 : CALL postortho(Ap, mo_coeff_array, Sc)
489 : !
490 : ! r_j+1 = r_j - alpha * Ap_j
491 69778 : DO ispin = 1, nspins
492 69778 : CALL cp_fm_scale_and_add(1.0_dp, r(ispin), -alpha, Ap(ispin))
493 : END DO
494 28440 : restart = .FALSE.
495 : END IF
496 : !
497 : ! preconditioner
498 28646 : linres_control%flag = ""
499 28646 : IF (linres_control%preconditioner_type == ot_precond_none) THEN
500 : !
501 : ! z_j+1 = r_j+1
502 0 : DO ispin = 1, nspins
503 0 : CALL cp_fm_to_fm(r(ispin), z(ispin))
504 : END DO
505 0 : linres_control%flag = "CG"
506 : ELSE
507 : !
508 : ! z_j+1 = M * r_j+1
509 70224 : DO ispin = 1, nspins
510 : CALL apply_preconditioner(p_env%preconditioner(ispin), r(ispin), &
511 70224 : z(ispin))
512 : END DO
513 28646 : linres_control%flag = "PCG"
514 : END IF
515 : !
516 70224 : DO ispin = 1, nspins
517 : !
518 : ! tr(r_j+1*z_j+1)
519 70224 : CALL cp_fm_trace(r(ispin), z(ispin), tr_rz1(ispin))
520 : END DO
521 70224 : IF (SUM(tr_rz1) < 0.0_dp) CPABORT("tr(r_j+1*z_j+1) < 0")
522 70224 : norm_res = SUM(tr_rz1)/SQRT(REAL(nspins*nao*maxnmo, dp))
523 : !
524 : ! beta = tr(r_j+1*z_j+1) / tr(r_j*z_j)
525 70224 : IF (SUM(tr_rz0) < 1.0e-10_dp) THEN
526 218 : beta = 0.0_dp
527 : ELSE
528 111148 : beta = SUM(tr_rz1)/SUM(tr_rz0)
529 : END IF
530 70224 : DO ispin = 1, nspins
531 : !
532 : ! p_j+1 = z_j+1 + beta * p_j
533 41578 : CALL cp_fm_scale_and_add(beta, p(ispin), 1.0_dp, z(ispin))
534 41578 : tr_rz00(ispin) = tr_rz0(ispin)
535 70224 : tr_rz0(ispin) = tr_rz1(ispin)
536 : END DO
537 : !
538 : ! Can we exit the SCF loop?
539 : CALL external_control(should_stop, "LINRES", target_time=qs_env%target_time, &
540 29966 : start_time=qs_env%start_time)
541 :
542 : END DO iteration
543 : !
544 : ! proj psi1
545 5100 : CALL preortho(psi1, mo_coeff_array, Sc)
546 : !
547 : !
548 : ! clean up
549 5100 : CALL cp_fm_release(r)
550 5100 : CALL cp_fm_release(p)
551 5100 : CALL cp_fm_release(z)
552 5100 : CALL cp_fm_release(Ap)
553 : !
554 5100 : CALL cp_fm_release(mo_coeff_array)
555 5100 : CALL cp_fm_release(Sc)
556 5100 : CALL cp_fm_release(chc)
557 : !
558 5100 : DEALLOCATE (tr_pAp, tr_rz0, tr_rz00, tr_rz1)
559 : !
560 5100 : CALL timestop(handle)
561 : !
562 15300 : END SUBROUTINE linres_solver
563 :
564 : ! **************************************************************************************************
565 : !> \brief ...
566 : !> \param qs_env ...
567 : !> \param p_env ...
568 : !> \param c0 ...
569 : !> \param v ...
570 : !> \param Av ...
571 : !> \param chc ...
572 : ! **************************************************************************************************
573 33952 : SUBROUTINE apply_op(qs_env, p_env, c0, v, Av, chc)
574 : !
575 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
576 : TYPE(qs_p_env_type) :: p_env
577 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: c0, v
578 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: Av
579 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: chc
580 :
581 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op'
582 :
583 : INTEGER :: handle, ispin, nc1, nc2, nc3, nc4, nr1, &
584 : nr2, nr3, nr4, nspins
585 : REAL(dp) :: chksum
586 33952 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
587 : TYPE(dft_control_type), POINTER :: dft_control
588 : TYPE(linres_control_type), POINTER :: linres_control
589 : TYPE(qs_rho_type), POINTER :: rho
590 :
591 33952 : CALL timeset(routineN, handle)
592 :
593 33952 : NULLIFY (dft_control, matrix_ks, matrix_s, linres_control)
594 : CALL get_qs_env(qs_env=qs_env, &
595 : matrix_ks=matrix_ks, &
596 : matrix_s=matrix_s, &
597 : dft_control=dft_control, &
598 33952 : linres_control=linres_control)
599 :
600 33952 : nspins = dft_control%nspins
601 :
602 82934 : DO ispin = 1, nspins
603 : !c0, v, Av, chc
604 48982 : CALL cp_fm_get_info(c0(ispin), ncol_global=nc1, nrow_global=nr1)
605 48982 : CALL cp_fm_get_info(v(ispin), ncol_global=nc2, nrow_global=nr2)
606 48982 : CALL cp_fm_get_info(Av(ispin), ncol_global=nc3, nrow_global=nr3)
607 48982 : CALL cp_fm_get_info(chc(ispin), ncol_global=nc4, nrow_global=nr4)
608 48982 : IF (.NOT. (nc1 == nc2 .AND. nr1 == nr2 .AND. nc1 == nc3 .AND. nr1 == nr3 &
609 82934 : .AND. nc4 == nr4 .AND. nc1 <= nc4)) THEN
610 : CALL cp_abort(__LOCATION__, &
611 0 : "Number of vectors inconsistent or CHC matrix too small")
612 : END IF
613 : END DO
614 :
615 : ! apply the uncoupled operator
616 82934 : DO ispin = 1, nspins
617 : CALL apply_op_1(v(ispin), Av(ispin), matrix_ks(ispin)%matrix, &
618 82934 : matrix_s(1)%matrix, chc(ispin))
619 : END DO
620 :
621 33952 : IF (linres_control%do_kernel) THEN
622 :
623 : ! build DM, refill p1, build_dm_response keeps sparse structure
624 20382 : DO ispin = 1, nspins
625 20382 : CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
626 : END DO
627 9622 : CALL build_dm_response(c0, v, p_env%p1)
628 :
629 9622 : chksum = 0.0_dp
630 20382 : DO ispin = 1, nspins
631 20382 : chksum = chksum + dbcsr_checksum(p_env%p1(ispin)%matrix)
632 : END DO
633 :
634 : ! skip the kernel if the DM is very small
635 9622 : IF (chksum > 1.0E-14_dp) THEN
636 :
637 8134 : CALL p_env_check_i_alloc(p_env, qs_env)
638 :
639 8134 : CALL p_env_update_rho(p_env, qs_env)
640 :
641 8134 : CALL get_qs_env(qs_env, rho=rho) ! that could be called before
642 8134 : CALL qs_rho_update_rho(rho, qs_env=qs_env) ! that could be called before
643 8134 : IF (dft_control%qs_control%gapw) THEN
644 1014 : CALL prepare_gapw_den(qs_env)
645 7120 : ELSEIF (dft_control%qs_control%gapw_xc) THEN
646 414 : CALL prepare_gapw_den(qs_env, do_rho0=.FALSE.)
647 : END IF
648 :
649 17244 : DO ispin = 1, nspins
650 9110 : CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
651 17244 : IF (ASSOCIATED(p_env%kpp1_admm)) CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
652 : END DO
653 :
654 8134 : CALL apply_op_2(qs_env, p_env, c0, Av)
655 :
656 : END IF
657 :
658 : END IF
659 :
660 33952 : CALL timestop(handle)
661 :
662 33952 : END SUBROUTINE apply_op
663 :
664 : ! **************************************************************************************************
665 : !> \brief ...
666 : !> \param v ...
667 : !> \param Av ...
668 : !> \param matrix_ks ...
669 : !> \param matrix_s ...
670 : !> \param chc ...
671 : ! **************************************************************************************************
672 146946 : SUBROUTINE apply_op_1(v, Av, matrix_ks, matrix_s, chc)
673 : !
674 : TYPE(cp_fm_type), INTENT(IN) :: v
675 : TYPE(cp_fm_type), INTENT(INOUT) :: Av
676 : TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s
677 : TYPE(cp_fm_type), INTENT(IN) :: chc
678 :
679 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op_1'
680 :
681 : INTEGER :: handle, ncol, nrow
682 : TYPE(cp_fm_type) :: buf
683 :
684 48982 : CALL timeset(routineN, handle)
685 : !
686 48982 : CALL cp_fm_create(buf, v%matrix_struct)
687 : !
688 48982 : CALL cp_fm_get_info(v, ncol_global=ncol, nrow_global=nrow)
689 : ! H * v
690 48982 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks, v, Av, ncol)
691 : ! v * e (chc already multiplied by -1)
692 48982 : CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, v, chc, 0.0_dp, buf)
693 : ! S * ve
694 48982 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, buf, Av, ncol, alpha=1.0_dp, beta=1.0_dp)
695 : !Results is H*C1 - S*<iHj>*C1
696 : !
697 48982 : CALL cp_fm_release(buf)
698 : !
699 48982 : CALL timestop(handle)
700 : !
701 48982 : END SUBROUTINE apply_op_1
702 :
703 : !MERGE
704 : ! **************************************************************************************************
705 : !> \brief ...
706 : !> \param v ...
707 : !> \param psi0 ...
708 : !> \param S_psi0 ...
709 : ! **************************************************************************************************
710 10200 : SUBROUTINE preortho(v, psi0, S_psi0)
711 : !v = (I-PS)v
712 : !
713 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: v, psi0, S_psi0
714 :
715 : CHARACTER(LEN=*), PARAMETER :: routineN = 'preortho'
716 :
717 : INTEGER :: handle, ispin, mp, mt, mv, np, nspins, &
718 : nt, nv
719 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
720 : TYPE(cp_fm_type) :: buf
721 :
722 10200 : CALL timeset(routineN, handle)
723 : !
724 10200 : NULLIFY (tmp_fm_struct)
725 : !
726 10200 : nspins = SIZE(v, 1)
727 : !
728 24528 : DO ispin = 1, nspins
729 14328 : CALL cp_fm_get_info(v(ispin), ncol_global=mv, nrow_global=nv)
730 14328 : CALL cp_fm_get_info(psi0(ispin), ncol_global=mp, nrow_global=np)
731 : !
732 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nv, ncol_global=mp, &
733 : para_env=v(ispin)%matrix_struct%para_env, &
734 14328 : context=v(ispin)%matrix_struct%context)
735 14328 : CALL cp_fm_create(buf, tmp_fm_struct)
736 14328 : CALL cp_fm_struct_release(tmp_fm_struct)
737 : !
738 14328 : CALL cp_fm_get_info(buf, ncol_global=mt, nrow_global=nt)
739 14328 : CPASSERT(nv == np)
740 14328 : CPASSERT(mt >= mv)
741 14328 : CPASSERT(mt >= mp)
742 14328 : CPASSERT(nt == nv)
743 : !
744 : ! buf = v' * S_psi0
745 14328 : CALL parallel_gemm('T', 'N', mv, mp, nv, 1.0_dp, v(ispin), S_psi0(ispin), 0.0_dp, buf)
746 : ! v = v - psi0 * buf'
747 14328 : CALL parallel_gemm('N', 'T', nv, mv, mp, -1.0_dp, psi0(ispin), buf, 1.0_dp, v(ispin))
748 : !
749 67512 : CALL cp_fm_release(buf)
750 : END DO
751 : !
752 10200 : CALL timestop(handle)
753 : !
754 10200 : END SUBROUTINE preortho
755 :
756 : ! **************************************************************************************************
757 : !> \brief projects first index of v onto the virtual subspace
758 : !> \param v matrix to be projected
759 : !> \param psi0 matrix with occupied orbitals
760 : !> \param S_psi0 matrix containing product of metric and occupied orbitals
761 : ! **************************************************************************************************
762 62392 : SUBROUTINE postortho(v, psi0, S_psi0)
763 : !v = (I-SP)v
764 : !
765 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: v, psi0, S_psi0
766 :
767 : CHARACTER(LEN=*), PARAMETER :: routineN = 'postortho'
768 :
769 : INTEGER :: handle, ispin, mp, mt, mv, np, nspins, &
770 : nt, nv
771 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
772 : TYPE(cp_fm_type) :: buf
773 :
774 62392 : CALL timeset(routineN, handle)
775 : !
776 62392 : NULLIFY (tmp_fm_struct)
777 : !
778 62392 : nspins = SIZE(v, 1)
779 : !
780 152712 : DO ispin = 1, nspins
781 90320 : CALL cp_fm_get_info(v(ispin), ncol_global=mv, nrow_global=nv)
782 90320 : CALL cp_fm_get_info(psi0(ispin), ncol_global=mp, nrow_global=np)
783 : !
784 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nv, ncol_global=mp, &
785 : para_env=v(ispin)%matrix_struct%para_env, &
786 90320 : context=v(ispin)%matrix_struct%context)
787 90320 : CALL cp_fm_create(buf, tmp_fm_struct)
788 90320 : CALL cp_fm_struct_release(tmp_fm_struct)
789 : !
790 90320 : CALL cp_fm_get_info(buf, ncol_global=mt, nrow_global=nt)
791 90320 : CPASSERT(nv == np)
792 90320 : CPASSERT(mt >= mv)
793 90320 : CPASSERT(mt >= mp)
794 90320 : CPASSERT(nt == nv)
795 : !
796 : ! buf = v' * psi0
797 90320 : CALL parallel_gemm('T', 'N', mv, mp, nv, 1.0_dp, v(ispin), psi0(ispin), 0.0_dp, buf)
798 : ! v = v - S_psi0 * buf'
799 90320 : CALL parallel_gemm('N', 'T', nv, mv, mp, -1.0_dp, S_psi0(ispin), buf, 1.0_dp, v(ispin))
800 : !
801 423672 : CALL cp_fm_release(buf)
802 : END DO
803 : !
804 62392 : CALL timestop(handle)
805 : !
806 62392 : END SUBROUTINE postortho
807 :
808 : ! **************************************************************************************************
809 : !> \brief ...
810 : !> \param qs_env ...
811 : !> \param linres_section ...
812 : !> \param vec ...
813 : !> \param ivec ...
814 : !> \param tag ...
815 : !> \param ind ...
816 : ! **************************************************************************************************
817 3522 : SUBROUTINE linres_write_restart(qs_env, linres_section, vec, ivec, tag, ind)
818 : TYPE(qs_environment_type), POINTER :: qs_env
819 : TYPE(section_vals_type), POINTER :: linres_section
820 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
821 : INTEGER, INTENT(IN) :: ivec
822 : CHARACTER(LEN=*) :: tag
823 : INTEGER, INTENT(IN), OPTIONAL :: ind
824 :
825 : CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_write_restart'
826 :
827 : CHARACTER(LEN=default_path_length) :: filename
828 : CHARACTER(LEN=default_string_length) :: my_middle, my_pos, my_status
829 : INTEGER :: handle, i, i_block, ia, ie, iounit, &
830 : ispin, j, max_block, nao, nmo, nspins, &
831 : rst_unit
832 3522 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
833 : TYPE(cp_fm_type), POINTER :: mo_coeff
834 : TYPE(cp_logger_type), POINTER :: logger
835 3522 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
836 : TYPE(mp_para_env_type), POINTER :: para_env
837 : TYPE(section_vals_type), POINTER :: print_key
838 :
839 3522 : NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
840 :
841 3522 : CALL timeset(routineN, handle)
842 :
843 3522 : logger => cp_get_default_logger()
844 :
845 3522 : IF (BTEST(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", &
846 : used_print_key=print_key), &
847 : cp_p_file)) THEN
848 :
849 : iounit = cp_print_key_unit_nr(logger, linres_section, &
850 3522 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
851 :
852 : CALL get_qs_env(qs_env=qs_env, &
853 : mos=mos, &
854 3522 : para_env=para_env)
855 :
856 3522 : nspins = SIZE(mos)
857 :
858 3522 : my_status = "REPLACE"
859 3522 : my_pos = "REWIND"
860 3522 : CALL XSTRING(tag, ia, ie)
861 3522 : IF (PRESENT(ind)) THEN
862 2460 : my_middle = "RESTART-"//tag(ia:ie)//TRIM(ADJUSTL(cp_to_string(ivec)))
863 : ELSE
864 1062 : my_middle = "RESTART-"//tag(ia:ie)
865 1062 : IF (ivec > 1) THEN
866 708 : my_status = "OLD"
867 708 : my_pos = "APPEND"
868 : END IF
869 : END IF
870 : rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", &
871 : extension=".lr", middle_name=TRIM(my_middle), file_status=TRIM(my_status), &
872 3522 : file_position=TRIM(my_pos), file_action="WRITE", file_form="UNFORMATTED")
873 :
874 : filename = cp_print_key_generate_filename(logger, print_key, &
875 3522 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
876 :
877 3522 : IF (iounit > 0) THEN
878 : WRITE (UNIT=iounit, FMT="(T2,A)") &
879 1761 : "LINRES| Writing response functions to the restart file <"//TRIM(ADJUSTL(filename))//">"
880 : END IF
881 :
882 : !
883 : ! write data to file
884 : ! use the scalapack block size as a default for buffering columns
885 3522 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
886 3522 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
887 14088 : ALLOCATE (vecbuffer(nao, max_block))
888 :
889 3522 : IF (PRESENT(ind)) THEN
890 2460 : IF (rst_unit > 0) WRITE (rst_unit) ind, ivec, nspins, nao
891 : ELSE
892 1062 : IF (rst_unit > 0) WRITE (rst_unit) ivec, nspins, nao
893 : END IF
894 :
895 8934 : DO ispin = 1, nspins
896 5412 : CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
897 :
898 5412 : IF (rst_unit > 0) WRITE (rst_unit) nmo
899 :
900 22494 : DO i = 1, nmo, MAX(max_block, 1)
901 8148 : i_block = MIN(max_block, nmo - i + 1)
902 8148 : CALL cp_fm_get_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
903 : ! doing this in one write would increase efficiency, but breaks RESTART compatibility.
904 : ! to old ones, and in cases where max_block is different between runs, as might happen during
905 : ! restarts with a different number of CPUs
906 50352 : DO j = 1, i_block
907 44940 : IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
908 : END DO
909 : END DO
910 : END DO
911 :
912 3522 : DEALLOCATE (vecbuffer)
913 :
914 : CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
915 7044 : "PRINT%RESTART")
916 : END IF
917 :
918 3522 : CALL timestop(handle)
919 :
920 3522 : END SUBROUTINE linres_write_restart
921 :
922 : ! **************************************************************************************************
923 : !> \brief ...
924 : !> \param qs_env ...
925 : !> \param linres_section ...
926 : !> \param vec ...
927 : !> \param ivec ...
928 : !> \param tag ...
929 : !> \param ind ...
930 : ! **************************************************************************************************
931 72 : SUBROUTINE linres_read_restart(qs_env, linres_section, vec, ivec, tag, ind)
932 : TYPE(qs_environment_type), POINTER :: qs_env
933 : TYPE(section_vals_type), POINTER :: linres_section
934 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
935 : INTEGER, INTENT(IN) :: ivec
936 : CHARACTER(LEN=*) :: tag
937 : INTEGER, INTENT(INOUT), OPTIONAL :: ind
938 :
939 : CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_read_restart'
940 :
941 : CHARACTER(LEN=default_path_length) :: filename
942 : CHARACTER(LEN=default_string_length) :: my_middle
943 : INTEGER :: handle, i, i_block, ia, ie, iostat, iounit, ispin, iv, iv1, ivec_tmp, j, &
944 : max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, nspins, nspins_tmp, rst_unit
945 : LOGICAL :: file_exists
946 72 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
947 : TYPE(cp_fm_type), POINTER :: mo_coeff
948 : TYPE(cp_logger_type), POINTER :: logger
949 72 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
950 : TYPE(mp_para_env_type), POINTER :: para_env
951 : TYPE(section_vals_type), POINTER :: print_key
952 :
953 72 : file_exists = .FALSE.
954 :
955 72 : CALL timeset(routineN, handle)
956 :
957 72 : NULLIFY (mos, para_env, logger, print_key, vecbuffer)
958 72 : logger => cp_get_default_logger()
959 :
960 : iounit = cp_print_key_unit_nr(logger, linres_section, &
961 72 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
962 :
963 : CALL get_qs_env(qs_env=qs_env, &
964 : para_env=para_env, &
965 72 : mos=mos)
966 :
967 72 : nspins = SIZE(mos)
968 :
969 72 : rst_unit = -1
970 72 : IF (para_env%is_source()) THEN
971 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
972 36 : n_rep_val=n_rep_val)
973 :
974 36 : CALL XSTRING(tag, ia, ie)
975 36 : IF (PRESENT(ind)) THEN
976 9 : my_middle = "RESTART-"//tag(ia:ie)//TRIM(ADJUSTL(cp_to_string(ivec)))
977 : ELSE
978 27 : my_middle = "RESTART-"//tag(ia:ie)
979 : END IF
980 :
981 36 : IF (n_rep_val > 0) THEN
982 18 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
983 18 : CALL xstring(filename, ia, ie)
984 18 : filename = filename(ia:ie)//TRIM(my_middle)//".lr"
985 : ELSE
986 : ! try to read from the filename that is generated automatically from the printkey
987 18 : print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
988 : filename = cp_print_key_generate_filename(logger, print_key, &
989 18 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
990 : END IF
991 36 : INQUIRE (FILE=filename, exist=file_exists)
992 : !
993 : ! open file
994 36 : IF (file_exists) THEN
995 : CALL open_file(file_name=TRIM(filename), &
996 : file_action="READ", &
997 : file_form="UNFORMATTED", &
998 : file_position="REWIND", &
999 : file_status="OLD", &
1000 15 : unit_number=rst_unit)
1001 :
1002 15 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
1003 15 : "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">"
1004 : ELSE
1005 21 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
1006 21 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
1007 : END IF
1008 : END IF
1009 :
1010 72 : CALL para_env%bcast(file_exists)
1011 :
1012 72 : IF (file_exists) THEN
1013 :
1014 30 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
1015 30 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
1016 :
1017 120 : ALLOCATE (vecbuffer(nao, max_block))
1018 : !
1019 : ! read headers
1020 30 : IF (PRESENT(ind)) THEN
1021 6 : iv1 = ivec
1022 : ELSE
1023 : iv1 = 1
1024 : END IF
1025 54 : DO iv = iv1, ivec
1026 :
1027 54 : IF (PRESENT(ind)) THEN
1028 6 : IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) ind, ivec_tmp, nspins_tmp, nao_tmp
1029 6 : CALL para_env%bcast(iostat)
1030 6 : CALL para_env%bcast(ind)
1031 : ELSE
1032 48 : IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) ivec_tmp, nspins_tmp, nao_tmp
1033 48 : CALL para_env%bcast(iostat)
1034 : END IF
1035 :
1036 54 : IF (iostat /= 0) EXIT
1037 54 : CALL para_env%bcast(ivec_tmp)
1038 54 : CALL para_env%bcast(nspins_tmp)
1039 54 : CALL para_env%bcast(nao_tmp)
1040 :
1041 : ! check that the number nao, nmo and nspins are
1042 : ! the same as in the current mos
1043 54 : IF (nspins_tmp /= nspins) CPABORT("nspins not consistent")
1044 54 : IF (nao_tmp /= nao) CPABORT("nao not consistent")
1045 : !
1046 108 : DO ispin = 1, nspins
1047 54 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
1048 54 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
1049 : !
1050 54 : IF (rst_unit > 0) READ (rst_unit) nmo_tmp
1051 54 : CALL para_env%bcast(nmo_tmp)
1052 54 : IF (nmo_tmp /= nmo) CPABORT("nmo not consistent")
1053 : !
1054 : ! read the response
1055 216 : DO i = 1, nmo, MAX(max_block, 1)
1056 54 : i_block = MIN(max_block, nmo - i + 1)
1057 270 : DO j = 1, i_block
1058 270 : IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
1059 : END DO
1060 108 : IF (iv == ivec_tmp) THEN
1061 10422 : CALL para_env%bcast(vecbuffer)
1062 54 : CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
1063 : END IF
1064 : END DO
1065 : END DO
1066 54 : IF (ivec == ivec_tmp) EXIT
1067 : END DO
1068 :
1069 30 : IF (iostat /= 0) THEN
1070 0 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
1071 0 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
1072 : END IF
1073 :
1074 60 : DEALLOCATE (vecbuffer)
1075 :
1076 : END IF
1077 :
1078 72 : IF (para_env%is_source()) THEN
1079 36 : IF (file_exists) CALL close_file(unit_number=rst_unit)
1080 : END IF
1081 :
1082 72 : CALL timestop(handle)
1083 :
1084 72 : END SUBROUTINE linres_read_restart
1085 :
1086 : ! **************************************************************************************************
1087 : !> \brief ...
1088 : !> \param p_env ...
1089 : !> \param linres_control ...
1090 : !> \param nspins ...
1091 : ! **************************************************************************************************
1092 5100 : SUBROUTINE check_p_env_init(p_env, linres_control, nspins)
1093 : !
1094 : TYPE(qs_p_env_type) :: p_env
1095 : TYPE(linres_control_type), INTENT(IN) :: linres_control
1096 : INTEGER, INTENT(IN) :: nspins
1097 :
1098 : INTEGER :: ispin, ncol, nrow
1099 :
1100 5100 : IF (linres_control%preconditioner_type /= ot_precond_none) THEN
1101 5100 : CPASSERT(ASSOCIATED(p_env%preconditioner))
1102 12264 : DO ispin = 1, nspins
1103 7164 : CALL cp_fm_get_info(p_env%PS_psi0(ispin), nrow_global=nrow, ncol_global=ncol)
1104 7164 : CPASSERT(nrow == p_env%n_ao(ispin))
1105 19428 : CPASSERT(ncol == p_env%n_mo(ispin))
1106 : END DO
1107 : END IF
1108 :
1109 5100 : END SUBROUTINE check_p_env_init
1110 :
1111 : END MODULE qs_linres_methods
|
