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