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