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