Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for the real time propagation.
10 : !> \author Florian Schiffmann (02.09)
11 : ! **************************************************************************************************
12 :
13 : MODULE rt_propagation
14 : USE bibliography, ONLY: Andermatt2016,&
15 : cite_reference
16 : USE cp_control_types, ONLY: dft_control_type,&
17 : rtp_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
19 : dbcsr_p_type
20 : USE cp_external_control, ONLY: external_control
21 : USE cp_fm_types, ONLY: cp_fm_set_all,&
22 : cp_fm_to_fm,&
23 : cp_fm_type
24 : USE cp_log_handling, ONLY: cp_get_default_logger,&
25 : cp_logger_type,&
26 : cp_to_string
27 : USE cp_output_handling, ONLY: cp_add_iter_level,&
28 : cp_iterate,&
29 : cp_p_file,&
30 : cp_print_key_generate_filename,&
31 : cp_print_key_should_output,&
32 : cp_print_key_unit_nr,&
33 : cp_rm_iter_level
34 : USE efield_utils, ONLY: calculate_ecore_efield
35 : USE force_env_methods, ONLY: force_env_calc_energy_force
36 : USE force_env_types, ONLY: force_env_get,&
37 : force_env_type
38 : USE global_types, ONLY: global_environment_type
39 : USE hfx_admm_utils, ONLY: hfx_admm_init
40 : USE input_constants, ONLY: real_time_propagation,&
41 : use_restart_wfn,&
42 : use_rt_restart,&
43 : use_scf_wfn
44 : USE input_cp2k_restarts, ONLY: write_restart
45 : USE input_section_types, ONLY: section_vals_get,&
46 : section_vals_get_subs_vals,&
47 : section_vals_type,&
48 : section_vals_val_get,&
49 : section_vals_val_set
50 : USE kinds, ONLY: default_path_length,&
51 : dp
52 : USE machine, ONLY: m_walltime
53 : USE md_environment_types, ONLY: md_environment_type
54 : USE pw_env_types, ONLY: pw_env_type
55 : USE qs_core_hamiltonian, ONLY: qs_matrix_h_allocate_imag_from_real
56 : USE qs_energy_init, ONLY: qs_energies_init
57 : USE qs_energy_types, ONLY: qs_energy_type
58 : USE qs_environment_types, ONLY: get_qs_env,&
59 : qs_environment_type
60 : USE qs_external_potential, ONLY: external_c_potential,&
61 : external_e_potential
62 : USE qs_ks_methods, ONLY: qs_ks_allocate_basics,&
63 : qs_ks_update_qs_env
64 : USE qs_ks_types, ONLY: qs_ks_did_change,&
65 : qs_ks_env_type,&
66 : set_ks_env
67 : USE qs_mo_io, ONLY: wfn_restart_file_name
68 : USE qs_mo_types, ONLY: get_mo_set,&
69 : init_mo_set,&
70 : mo_set_type
71 : USE qs_rho_methods, ONLY: allocate_rho_ao_imag_from_real
72 : USE qs_rho_types, ONLY: qs_rho_set,&
73 : qs_rho_type
74 : USE rt_delta_pulse, ONLY: apply_delta_pulse
75 : USE rt_hfx_utils, ONLY: rtp_hfx_rebuild
76 : USE rt_projection_mo_utils, ONLY: init_mo_projection
77 : USE rt_propagation_methods, ONLY: propagation_step
78 : USE rt_propagation_output, ONLY: rt_prop_output
79 : USE rt_propagation_types, ONLY: get_rtp,&
80 : rt_prop_create,&
81 : rt_prop_type,&
82 : rtp_create_SinvH_imag,&
83 : rtp_history_create
84 : USE rt_propagation_utils, ONLY: calc_S_derivs,&
85 : calc_update_rho,&
86 : calc_update_rho_sparse,&
87 : get_restart_wfn,&
88 : warn_section_unused
89 : USE rt_propagation_velocity_gauge, ONLY: velocity_gauge_ks_matrix
90 : USE rt_propagator_init, ONLY: init_propagators,&
91 : rt_initialize_rho_from_mos
92 : #include "../base/base_uses.f90"
93 :
94 : IMPLICIT NONE
95 :
96 : PRIVATE
97 :
98 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation'
99 :
100 : PUBLIC :: rt_prop_setup
101 :
102 : CONTAINS
103 :
104 : ! **************************************************************************************************
105 : !> \brief creates rtp_type, gets the initial state, either by reading MO's
106 : !> from file or calling SCF run
107 : !> \param force_env ...
108 : !> \author Florian Schiffmann (02.09)
109 : ! **************************************************************************************************
110 :
111 594 : SUBROUTINE rt_prop_setup(force_env)
112 : TYPE(force_env_type), POINTER :: force_env
113 :
114 : INTEGER :: aspc_order
115 : LOGICAL :: magnetic, vel_reprs
116 : TYPE(dft_control_type), POINTER :: dft_control
117 : TYPE(global_environment_type), POINTER :: globenv
118 : TYPE(qs_energy_type), POINTER :: energy
119 : TYPE(qs_environment_type), POINTER :: qs_env
120 : TYPE(rt_prop_type), POINTER :: rtp
121 : TYPE(rtp_control_type), POINTER :: rtp_control
122 : TYPE(section_vals_type), POINTER :: hfx_sections, input, ls_scf_section, &
123 : md_section, motion_section, &
124 : print_moments_section, &
125 : rtp_print_section
126 :
127 198 : NULLIFY (qs_env, rtp_control, dft_control)
128 :
129 198 : CALL cite_reference(Andermatt2016)
130 :
131 198 : CALL force_env_get(force_env=force_env, qs_env=qs_env, globenv=globenv)
132 198 : CALL get_qs_env(qs_env, dft_control=dft_control, energy=energy)
133 198 : rtp_control => dft_control%rtp_control
134 :
135 : ! Takes care that an initial wavefunction/density is available
136 : ! Can either be by performing an scf loop or reading a restart
137 198 : CALL rt_initial_guess(qs_env, force_env, rtp_control)
138 :
139 : ! Initializes the extrapolation
140 198 : NULLIFY (rtp)
141 198 : CALL get_qs_env(qs_env=qs_env, rtp=rtp, input=input)
142 198 : aspc_order = rtp_control%aspc_order
143 198 : CALL rtp_history_create(rtp, aspc_order)
144 :
145 : ! Reads the simulation parameters from the input
146 198 : motion_section => section_vals_get_subs_vals(force_env%root_section, "MOTION")
147 198 : md_section => section_vals_get_subs_vals(motion_section, "MD")
148 198 : hfx_sections => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%DFT%XC%HF")
149 198 : print_moments_section => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%DFT%PRINT%MOMENTS")
150 198 : CALL section_vals_val_get(md_section, "TIMESTEP", r_val=qs_env%rtp%dt)
151 198 : CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=qs_env%rtp%i_start)
152 198 : CALL section_vals_val_get(md_section, "STEPS", i_val=rtp%nsteps)
153 198 : CALL section_vals_val_get(md_section, "MAX_STEPS", i_val=rtp%max_steps)
154 :
155 198 : ls_scf_section => section_vals_get_subs_vals(input, "DFT%LS_SCF")
156 198 : CALL section_vals_val_get(ls_scf_section, "EPS_FILTER", r_val=rtp%filter_eps)
157 198 : IF (.NOT. qs_env%rtp%linear_scaling) rtp%filter_eps = 0.0_dp
158 198 : IF (rtp_control%acc_ref < 1) rtp_control%acc_ref = 1
159 198 : rtp%filter_eps_small = rtp%filter_eps/rtp_control%acc_ref
160 198 : CALL section_vals_val_get(ls_scf_section, "EPS_LANCZOS", r_val=rtp%lanzcos_threshold)
161 198 : CALL section_vals_val_get(ls_scf_section, "MAX_ITER_LANCZOS", i_val=rtp%lanzcos_max_iter)
162 198 : CALL section_vals_val_get(ls_scf_section, "SIGN_SQRT_ORDER", i_val=rtp%newton_schulz_order)
163 198 : CALL section_vals_get(hfx_sections, explicit=rtp%do_hfx)
164 198 : CALL section_vals_val_get(print_moments_section, "MAGNETIC", l_val=magnetic)
165 198 : CALL section_vals_val_get(print_moments_section, "VEL_REPRS", l_val=vel_reprs)
166 :
167 : rtp%track_imag_density = (magnetic) .OR. (vel_reprs) .OR. (rtp_control%velocity_gauge) &
168 198 : .OR. (rtp%do_hfx) .OR. (.NOT. rtp_control%fixed_ions)
169 198 : rtp%propagate_complex_ks = rtp%do_hfx .OR. rtp_control%velocity_gauge
170 :
171 : ! Marek : In case some print sections that apply so far only to RTBSE are present,
172 : ! warn the user that the quantities will not be in fact printed out
173 : rtp_print_section => section_vals_get_subs_vals(force_env%root_section, &
174 198 : "FORCE_EVAL%DFT%REAL_TIME_PROPAGATION%PRINT")
175 : CALL warn_section_unused(rtp_print_section, "POLARIZABILITY", &
176 198 : "POLARIZABILITY printing not implemented for non-RTBSE code.")
177 : CALL warn_section_unused(rtp_print_section, "MOMENTS", &
178 198 : "MOMENTS not implemented: Use DFT%PRINT%MOMENTS for moments printing with TDDFT.")
179 : CALL warn_section_unused(rtp_print_section, "MOMENTS_FT", &
180 198 : "MOMENTS_FT printing not implemented for non-RTBSE code.")
181 : CALL warn_section_unused(rtp_print_section, "DENSITY_MATRIX", &
182 198 : "DENSITY_MATRIX printing not implemented for non-RTBSE code.")
183 :
184 : CALL rt_init_complex_quantities(qs_env, imag_p=rtp%track_imag_density, &
185 198 : imag_ks=rtp%propagate_complex_ks, imag_h=rtp_control%velocity_gauge)
186 :
187 : ! Hmm, not really like to initialize with the structure of S but I reckon it is
188 : ! done everywhere like this
189 198 : IF (rtp%do_hfx) CALL rtp_hfx_rebuild(qs_env)
190 :
191 : ! Setup the MO projection environment if required
192 198 : IF (rtp_control%is_proj_mo) CALL init_mo_projection(qs_env, rtp_control)
193 :
194 198 : CALL init_propagation_run(qs_env)
195 198 : IF (.NOT. rtp_control%fixed_ions) THEN
196 : !derivativs of the overlap needed for EMD
197 72 : CALL calc_S_derivs(qs_env)
198 : ! a bit hidden, but computes SinvH and SinvB (calc_SinvH for CN,EM and ARNOLDI)
199 : ! make_etrs_exp in case of ETRS in combination with TAYLOR and PADE
200 : END IF
201 198 : CALL init_propagators(qs_env)
202 198 : IF (rtp_control%fixed_ions) THEN
203 126 : CALL run_propagation(qs_env, force_env, globenv)
204 : ELSE
205 72 : rtp_control%initial_step = .TRUE.
206 72 : CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
207 72 : rtp_control%initial_step = .FALSE.
208 72 : rtp%energy_old = energy%total
209 : END IF
210 :
211 198 : END SUBROUTINE rt_prop_setup
212 :
213 : ! **************************************************************************************************
214 : !> \brief calculates the matrices needed in the first step of EMD/RTP
215 : !> \param qs_env ...
216 : !> \author Florian Schiffmann (02.09)
217 : ! **************************************************************************************************
218 :
219 198 : SUBROUTINE init_propagation_run(qs_env)
220 : TYPE(qs_environment_type), POINTER :: qs_env
221 :
222 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
223 :
224 : INTEGER :: i, ispin, re
225 : INTEGER, DIMENSION(2) :: nelectron_spin
226 198 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_old
227 198 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, rho_new, rho_old
228 : TYPE(dft_control_type), POINTER :: dft_control
229 198 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
230 : TYPE(rt_prop_type), POINTER :: rtp
231 : TYPE(rtp_control_type), POINTER :: rtp_control
232 :
233 198 : NULLIFY (dft_control, rtp, rtp_control)
234 :
235 198 : CALL cite_reference(Andermatt2016)
236 :
237 : CALL get_qs_env(qs_env, &
238 : rtp=rtp, &
239 198 : dft_control=dft_control)
240 198 : rtp_control => dft_control%rtp_control
241 :
242 198 : IF (rtp_control%initial_wfn == use_scf_wfn) THEN
243 162 : IF (rtp_control%apply_delta_pulse .OR. rtp_control%apply_delta_pulse_mag) THEN
244 54 : CALL apply_delta_pulse(qs_env, rtp, rtp_control)
245 : ELSE
246 108 : IF (.NOT. rtp%linear_scaling) THEN
247 74 : CALL get_rtp(rtp=rtp, mos_old=mos_old)
248 74 : CALL get_qs_env(qs_env, mos=mos)
249 162 : DO i = 1, SIZE(mos)
250 88 : CALL cp_fm_to_fm(mos(i)%mo_coeff, mos_old(2*i - 1))
251 162 : CALL cp_fm_set_all(mos_old(2*i), zero, zero)
252 : END DO
253 : END IF
254 : END IF
255 : END IF
256 :
257 198 : IF (.NOT. rtp%linear_scaling) THEN
258 108 : CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
259 384 : DO i = 1, SIZE(mos_old)
260 384 : CALL cp_fm_to_fm(mos_old(i), mos_new(i))
261 : END DO
262 108 : CALL calc_update_rho(qs_env)
263 : ELSE
264 90 : IF (rtp_control%initial_wfn == use_scf_wfn) THEN
265 : CALL get_qs_env(qs_env, &
266 : matrix_ks=matrix_ks, &
267 : mos=mos, &
268 74 : nelectron_spin=nelectron_spin)
269 74 : IF (ASSOCIATED(mos)) THEN
270 : !The wavefunction was minimized by an mo based algorith. P is therefore calculated from the mos
271 64 : IF (ASSOCIATED(rtp%mos)) THEN
272 40 : IF (ASSOCIATED(rtp%mos%old)) THEN
273 : ! Delta kick was applied and the results is in rtp%mos%old
274 40 : CALL rt_initialize_rho_from_mos(rtp, mos, mos_old=rtp%mos%old)
275 : ELSE
276 0 : CALL rt_initialize_rho_from_mos(rtp, mos)
277 : END IF
278 : ELSE
279 24 : CALL rt_initialize_rho_from_mos(rtp, mos)
280 : END IF
281 : ELSE
282 : !The wavefunction was minimized using a linear scaling method. The density matrix is therefore taken from the ls_scf_env.
283 10 : CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
284 24 : DO ispin = 1, SIZE(rho_old)/2
285 14 : re = 2*ispin - 1
286 14 : CALL dbcsr_copy(rho_old(re)%matrix, qs_env%ls_scf_env%matrix_p(ispin))
287 24 : CALL dbcsr_copy(rho_new(re)%matrix, qs_env%ls_scf_env%matrix_p(ispin))
288 : END DO
289 : END IF
290 74 : CALL calc_update_rho_sparse(qs_env)
291 : END IF
292 : END IF
293 : ! Modify KS matrix to include the additional terms in the velocity gauge
294 198 : IF (rtp_control%velocity_gauge) THEN
295 : ! As matrix_h and matrix_h_im are not updated by qs_ks_update_qs_env()
296 : ! the non-gauge transformed non-local part has to be subtracted here
297 8 : CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.TRUE.)
298 : END IF
299 198 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.)
300 :
301 198 : END SUBROUTINE init_propagation_run
302 :
303 : ! **************************************************************************************************
304 : !> \brief performs the real RTP run, gets information from MD section
305 : !> uses MD as iteration level
306 : !> \param qs_env ...
307 : !> \param force_env ...
308 : !> \param globenv ...
309 : !> \author Florian Schiffmann (02.09)
310 : ! **************************************************************************************************
311 :
312 126 : SUBROUTINE run_propagation(qs_env, force_env, globenv)
313 : TYPE(qs_environment_type), POINTER :: qs_env
314 : TYPE(force_env_type), POINTER :: force_env
315 : TYPE(global_environment_type), POINTER :: globenv
316 :
317 : CHARACTER(len=*), PARAMETER :: routineN = 'run_propagation'
318 :
319 : INTEGER :: aspc_order, handle, i_iter, i_step, &
320 : max_iter, max_steps, output_unit
321 : LOGICAL :: should_stop
322 : REAL(Kind=dp) :: eps_ener, time_iter_start, &
323 : time_iter_stop, used_time
324 : TYPE(cp_logger_type), POINTER :: logger
325 : TYPE(dft_control_type), POINTER :: dft_control
326 : TYPE(pw_env_type), POINTER :: pw_env
327 : TYPE(qs_energy_type), POINTER :: energy
328 : TYPE(rt_prop_type), POINTER :: rtp
329 : TYPE(rtp_control_type), POINTER :: rtp_control
330 : TYPE(section_vals_type), POINTER :: input, rtp_section
331 :
332 126 : should_stop = .FALSE.
333 126 : CALL timeset(routineN, handle)
334 :
335 126 : CALL cite_reference(Andermatt2016)
336 :
337 126 : NULLIFY (logger, dft_control, energy, rtp, rtp_control, input, rtp_section)
338 126 : logger => cp_get_default_logger()
339 :
340 126 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, rtp=rtp, energy=energy, input=input)
341 :
342 126 : rtp_control => dft_control%rtp_control
343 126 : max_steps = MIN(rtp%nsteps, rtp%max_steps)
344 126 : max_iter = rtp_control%max_iter
345 126 : eps_ener = rtp_control%eps_ener
346 :
347 126 : aspc_order = rtp_control%aspc_order
348 :
349 126 : rtp%energy_old = energy%total
350 126 : time_iter_start = m_walltime()
351 126 : CALL cp_add_iter_level(logger%iter_info, "MD")
352 126 : CALL cp_iterate(logger%iter_info, iter_nr=0)
353 126 : IF (rtp%i_start >= max_steps) CALL cp_abort(__LOCATION__, &
354 0 : "maximum step number smaller than initial step value")
355 :
356 126 : rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
357 : output_unit = cp_print_key_unit_nr(logger, rtp_section, "PRINT%PROGRAM_RUN_INFO", &
358 126 : extension=".scfLog")
359 :
360 470 : DO i_step = rtp%i_start + 1, max_steps
361 344 : IF (output_unit > 0) THEN
362 : WRITE (output_unit, FMT="(/,(T2,A,T40,I6))") &
363 172 : "Real time propagation step:", i_step
364 : END IF
365 344 : energy%efield_core = 0.0_dp
366 344 : qs_env%sim_time = REAL(i_step, dp)*rtp%dt
367 344 : CALL get_qs_env(qs_env, pw_env=pw_env)
368 344 : pw_env%poisson_env%parameters%dbc_params%time = qs_env%sim_time
369 344 : qs_env%sim_step = i_step
370 344 : rtp%istep = i_step - rtp%i_start
371 344 : CALL calculate_ecore_efield(qs_env, .FALSE.)
372 344 : IF (dft_control%apply_external_potential) THEN
373 0 : IF (.NOT. dft_control%expot_control%static) THEN
374 0 : dft_control%eval_external_potential = .TRUE.
375 : END IF
376 : END IF
377 344 : CALL external_c_potential(qs_env, calculate_forces=.FALSE.)
378 344 : CALL external_e_potential(qs_env)
379 344 : CALL cp_iterate(logger%iter_info, last=(i_step == max_steps), iter_nr=i_step)
380 344 : rtp%converged = .FALSE.
381 1140 : DO i_iter = 1, max_iter
382 1140 : IF (i_step == rtp%i_start + 1 .AND. i_iter == 2 .AND. rtp_control%hfx_redistribute) &
383 0 : CALL qs_ks_did_change(qs_env%ks_env, s_mstruct_changed=.TRUE.)
384 1140 : rtp%iter = i_iter
385 1140 : CALL propagation_step(qs_env, rtp, rtp_control)
386 1140 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.)
387 1140 : rtp%energy_new = energy%total
388 1140 : IF (rtp%converged) EXIT
389 1140 : CALL rt_prop_output(qs_env, real_time_propagation, rtp%delta_iter)
390 : END DO
391 470 : IF (rtp%converged) THEN
392 344 : CALL external_control(should_stop, "MD", globenv=globenv)
393 344 : IF (should_stop) CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=i_step)
394 344 : time_iter_stop = m_walltime()
395 344 : used_time = time_iter_stop - time_iter_start
396 344 : time_iter_start = time_iter_stop
397 344 : CALL rt_prop_output(qs_env, real_time_propagation, delta_iter=rtp%delta_iter, used_time=used_time)
398 344 : CALL rt_write_input_restart(force_env=force_env, qs_env=qs_env)
399 344 : IF (should_stop) EXIT
400 : ELSE
401 : EXIT
402 : END IF
403 : END DO
404 126 : CALL cp_rm_iter_level(logger%iter_info, "MD")
405 :
406 126 : IF (.NOT. rtp%converged) &
407 : CALL cp_abort(__LOCATION__, "propagation did not converge, "// &
408 0 : "either increase MAX_ITER or use a smaller TIMESTEP")
409 :
410 126 : CALL timestop(handle)
411 :
412 126 : END SUBROUTINE run_propagation
413 :
414 : ! **************************************************************************************************
415 : !> \brief overwrites some values in the input file such that the .restart
416 : !> file will contain the appropriate information
417 : !> \param md_env ...
418 : !> \param qs_env ...
419 : !> \param force_env ...
420 : !> \author Florian Schiffmann (02.09)
421 : ! **************************************************************************************************
422 :
423 344 : SUBROUTINE rt_write_input_restart(md_env, qs_env, force_env)
424 : TYPE(md_environment_type), OPTIONAL, POINTER :: md_env
425 : TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
426 : TYPE(force_env_type), POINTER :: force_env
427 :
428 : CHARACTER(len=default_path_length) :: file_name
429 344 : REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_vals
430 : TYPE(cp_logger_type), POINTER :: logger
431 : TYPE(dft_control_type), POINTER :: dft_control
432 : TYPE(section_vals_type), POINTER :: dft_section, efield_section, &
433 : motion_section, print_key, &
434 : root_section, rt_section
435 :
436 344 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
437 344 : root_section => force_env%root_section
438 344 : motion_section => section_vals_get_subs_vals(root_section, "MOTION")
439 344 : dft_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL%DFT")
440 344 : rt_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL%DFT%REAL_TIME_PROPAGATION")
441 :
442 344 : CALL section_vals_val_set(rt_section, "INITIAL_WFN", i_val=use_rt_restart)
443 344 : CALL section_vals_val_set(rt_section, "APPLY_DELTA_PULSE", l_val=.FALSE.)
444 344 : CALL section_vals_val_set(rt_section, "APPLY_DELTA_PULSE_MAG", l_val=.FALSE.)
445 344 : CALL section_vals_val_set(rt_section, "APPLY_WFN_MIX_INIT_RESTART", l_val=.FALSE.)
446 :
447 344 : logger => cp_get_default_logger()
448 :
449 : ! to continue propagating the TD wavefunction we need to read from the new .rtpwfn
450 344 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
451 : rt_section, "PRINT%RESTART"), cp_p_file)) THEN
452 126 : print_key => section_vals_get_subs_vals(rt_section, "PRINT%RESTART")
453 : file_name = cp_print_key_generate_filename(logger, print_key, &
454 126 : extension=".rtpwfn", my_local=.FALSE.)
455 126 : CALL section_vals_val_set(dft_section, "WFN_RESTART_FILE_NAME", c_val=TRIM(file_name))
456 : END IF
457 :
458 : ! coming from RTP
459 344 : IF (.NOT. PRESENT(md_env)) THEN
460 344 : CALL section_vals_val_set(motion_section, "MD%STEP_START_VAL", i_val=force_env%qs_env%sim_step)
461 : END IF
462 :
463 344 : IF (dft_control%apply_vector_potential) THEN
464 22 : efield_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL%DFT%EFIELD")
465 : NULLIFY (tmp_vals)
466 22 : ALLOCATE (tmp_vals(3))
467 88 : tmp_vals = dft_control%efield_fields(1)%efield%vec_pot_initial
468 : CALL section_vals_val_set(efield_section, "VEC_POT_INITIAL", &
469 : r_vals_ptr=tmp_vals, &
470 22 : i_rep_section=1)
471 : END IF
472 :
473 344 : CALL write_restart(md_env=md_env, root_section=root_section)
474 :
475 344 : END SUBROUTINE rt_write_input_restart
476 :
477 : ! **************************************************************************************************
478 : !> \brief Creates the initial electronic states and allocates the necessary
479 : !> matrices
480 : !> \param qs_env ...
481 : !> \param force_env ...
482 : !> \param rtp_control ...
483 : !> \author Florian Schiffmann (02.09)
484 : ! **************************************************************************************************
485 :
486 198 : SUBROUTINE rt_initial_guess(qs_env, force_env, rtp_control)
487 : TYPE(qs_environment_type), POINTER :: qs_env
488 : TYPE(force_env_type), POINTER :: force_env
489 : TYPE(rtp_control_type), POINTER :: rtp_control
490 :
491 : INTEGER :: homo, ispin
492 : LOGICAL :: energy_consistency
493 : TYPE(cp_fm_type), POINTER :: mo_coeff
494 198 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
495 : TYPE(dft_control_type), POINTER :: dft_control
496 :
497 198 : NULLIFY (matrix_s, dft_control)
498 198 : CALL get_qs_env(qs_env, dft_control=dft_control)
499 :
500 360 : SELECT CASE (rtp_control%initial_wfn)
501 : CASE (use_scf_wfn)
502 162 : qs_env%sim_time = 0.0_dp
503 162 : qs_env%sim_step = 0
504 162 : energy_consistency = .TRUE.
505 : !in the linear scaling case we need a correct kohn-sham matrix, which we cannot get with consistent energies
506 162 : IF (rtp_control%linear_scaling) energy_consistency = .FALSE.
507 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., &
508 162 : consistent_energies=energy_consistency)
509 162 : qs_env%run_rtp = .TRUE.
510 162 : ALLOCATE (qs_env%rtp)
511 162 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
512 162 : IF (dft_control%do_admm) THEN
513 6 : CALL hfx_admm_init(qs_env)
514 : CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
515 6 : rtp_control%linear_scaling, qs_env%admm_env%mos_aux_fit)
516 : ELSE
517 : CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
518 156 : rtp_control%linear_scaling)
519 : END IF
520 :
521 : CASE (use_restart_wfn, use_rt_restart)
522 36 : CALL qs_energies_init(qs_env, .FALSE.)
523 36 : IF (.NOT. rtp_control%linear_scaling .OR. rtp_control%initial_wfn == use_restart_wfn) THEN
524 86 : DO ispin = 1, SIZE(qs_env%mos)
525 52 : CALL get_mo_set(qs_env%mos(ispin), mo_coeff=mo_coeff, homo=homo)
526 86 : IF (.NOT. ASSOCIATED(mo_coeff)) THEN
527 : CALL init_mo_set(qs_env%mos(ispin), &
528 : qs_env%mpools%ao_mo_fm_pools(ispin)%pool, &
529 52 : name="qs_env%mo"//TRIM(ADJUSTL(cp_to_string(ispin))))
530 : END IF
531 : END DO
532 34 : IF (dft_control%do_admm) CALL hfx_admm_init(qs_env)
533 : END IF
534 36 : ALLOCATE (qs_env%rtp)
535 36 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
536 : CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
537 36 : rtp_control%linear_scaling, qs_env%admm_env%mos_aux_fit)
538 36 : CALL get_restart_wfn(qs_env)
539 :
540 234 : qs_env%run_rtp = .TRUE.
541 : END SELECT
542 :
543 198 : END SUBROUTINE rt_initial_guess
544 :
545 : ! **************************************************************************************************
546 : !> \brief ...
547 : !> \param qs_env ...
548 : !> \param imag_p ...
549 : !> \param imag_ks ...
550 : !> \param imag_h ...
551 : ! **************************************************************************************************
552 198 : SUBROUTINE rt_init_complex_quantities(qs_env, imag_p, imag_ks, imag_h)
553 : TYPE(qs_environment_type), POINTER :: qs_env
554 : LOGICAL, INTENT(in) :: imag_p, imag_ks, imag_h
555 :
556 : TYPE(dft_control_type), POINTER :: dft_control
557 : TYPE(qs_ks_env_type), POINTER :: ks_env
558 : TYPE(qs_rho_type), POINTER :: rho
559 : TYPE(rt_prop_type), POINTER :: rtp
560 :
561 198 : NULLIFY (ks_env, rho, dft_control)
562 :
563 : CALL get_qs_env(qs_env, &
564 : dft_control=dft_control, &
565 : ks_env=ks_env, &
566 : rho=rho, &
567 198 : rtp=rtp)
568 :
569 : ! rho
570 198 : CALL qs_rho_set(rho, complex_rho_ao=imag_p)
571 198 : IF (imag_p) CALL allocate_rho_ao_imag_from_real(rho, qs_env)
572 :
573 : ! ks
574 198 : CALL set_ks_env(ks_env, complex_ks=imag_ks)
575 198 : IF (imag_ks) THEN
576 40 : CALL qs_ks_allocate_basics(qs_env, is_complex=imag_ks)
577 40 : IF (.NOT. dft_control%rtp_control%fixed_ions) &
578 22 : CALL rtp_create_SinvH_imag(rtp, dft_control%nspins)
579 : END IF
580 :
581 : ! h
582 198 : IF (imag_h) CALL qs_matrix_h_allocate_imag_from_real(qs_env)
583 :
584 198 : END SUBROUTINE rt_init_complex_quantities
585 :
586 : END MODULE rt_propagation
|
