29 USE dbcsr_api,
ONLY: dbcsr_copy,&
86 #include "../base/base_uses.f90"
92 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'rt_propagation'
106 TYPE(force_env_type),
POINTER :: force_env
108 INTEGER :: aspc_order
109 LOGICAL :: magnetic, vel_reprs
110 TYPE(dft_control_type),
POINTER :: dft_control
111 TYPE(global_environment_type),
POINTER :: globenv
112 TYPE(qs_energy_type),
POINTER :: energy
113 TYPE(qs_environment_type),
POINTER :: qs_env
114 TYPE(rt_prop_type),
POINTER :: rtp
115 TYPE(rtp_control_type),
POINTER :: rtp_control
116 TYPE(section_vals_type),
POINTER :: hfx_sections, input, ls_scf_section, &
117 md_section, motion_section, &
118 print_moments_section
120 NULLIFY (qs_env, rtp_control, dft_control)
124 CALL force_env_get(force_env=force_env, qs_env=qs_env, globenv=globenv)
125 CALL get_qs_env(qs_env, dft_control=dft_control, energy=energy)
126 rtp_control => dft_control%rtp_control
130 CALL rt_initial_guess(qs_env, force_env, rtp_control)
134 CALL get_qs_env(qs_env=qs_env, rtp=rtp, input=input)
135 aspc_order = rtp_control%aspc_order
150 IF (.NOT. qs_env%rtp%linear_scaling) rtp%filter_eps = 0.0_dp
151 IF (rtp_control%acc_ref < 1) rtp_control%acc_ref = 1
152 rtp%filter_eps_small = rtp%filter_eps/rtp_control%acc_ref
160 rtp%track_imag_density = (magnetic) .OR. (vel_reprs) .OR. (rtp_control%velocity_gauge) &
161 .OR. (rtp%do_hfx) .OR. (.NOT. rtp_control%fixed_ions)
162 rtp%propagate_complex_ks = rtp%do_hfx .OR. rtp_control%velocity_gauge
164 CALL rt_init_complex_quantities(qs_env, imag_p=rtp%track_imag_density, &
165 imag_ks=rtp%propagate_complex_ks, imag_h=rtp_control%velocity_gauge)
174 CALL init_propagation_run(qs_env)
175 IF (.NOT. rtp_control%fixed_ions)
THEN
182 IF (rtp_control%fixed_ions)
THEN
183 CALL run_propagation(qs_env, force_env, globenv)
185 rtp_control%initial_step = .true.
187 rtp_control%initial_step = .false.
188 rtp%energy_old = energy%total
199 SUBROUTINE init_propagation_run(qs_env)
200 TYPE(qs_environment_type),
POINTER :: qs_env
202 REAL(kind=
dp),
PARAMETER :: zero = 0.0_dp
204 INTEGER :: i, ispin, re
205 INTEGER,
DIMENSION(2) :: nelectron_spin
206 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: mos_new, mos_old
207 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, rho_new, rho_old
208 TYPE(dft_control_type),
POINTER :: dft_control
209 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
210 TYPE(rt_prop_type),
POINTER :: rtp
211 TYPE(rtp_control_type),
POINTER :: rtp_control
213 NULLIFY (dft_control, rtp, rtp_control)
219 dft_control=dft_control)
220 rtp_control => dft_control%rtp_control
223 IF (rtp_control%apply_delta_pulse .OR. rtp_control%apply_delta_pulse_mag)
THEN
226 IF (.NOT. rtp%linear_scaling)
THEN
227 CALL get_rtp(rtp=rtp, mos_old=mos_old)
230 CALL cp_fm_to_fm(mos(i)%mo_coeff, mos_old(2*i - 1))
237 IF (.NOT. rtp%linear_scaling)
THEN
238 CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
239 DO i = 1,
SIZE(mos_old)
240 CALL cp_fm_to_fm(mos_old(i), mos_new(i))
246 matrix_ks=matrix_ks, &
248 nelectron_spin=nelectron_spin)
249 IF (
ASSOCIATED(mos))
THEN
251 IF (
ASSOCIATED(rtp%mos))
THEN
252 IF (
ASSOCIATED(rtp%mos%old))
THEN
263 CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
264 DO ispin = 1,
SIZE(rho_old)/2
266 CALL dbcsr_copy(rho_old(re)%matrix, qs_env%ls_scf_env%matrix_p(ispin))
267 CALL dbcsr_copy(rho_new(re)%matrix, qs_env%ls_scf_env%matrix_p(ispin))
274 IF (rtp_control%velocity_gauge)
THEN
281 END SUBROUTINE init_propagation_run
292 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
297 CHARACTER(len=*),
PARAMETER :: routinen =
'run_propagation'
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
312 should_stop = .false.
313 CALL timeset(routinen, handle)
317 NULLIFY (logger, dft_control, energy, rtp, rtp_control, input, rtp_section)
320 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, rtp=rtp, energy=energy, input=input)
322 rtp_control => dft_control%rtp_control
323 max_steps = min(rtp%nsteps, rtp%max_steps)
324 max_iter = rtp_control%max_iter
325 eps_ener = rtp_control%eps_ener
327 aspc_order = rtp_control%aspc_order
329 rtp%energy_old = energy%total
333 IF (rtp%i_start >= max_steps)
CALL cp_abort(__location__, &
334 "maximum step number smaller than initial step value")
340 DO i_step = rtp%i_start + 1, max_steps
341 IF (output_unit > 0)
THEN
342 WRITE (output_unit, fmt=
"(/,(T2,A,T40,I6))") &
343 "Real time propagation step:", i_step
345 energy%efield_core = 0.0_dp
346 qs_env%sim_time = real(i_step,
dp)*rtp%dt
348 pw_env%poisson_env%parameters%dbc_params%time = qs_env%sim_time
349 qs_env%sim_step = i_step
350 rtp%istep = i_step - rtp%i_start
352 IF (dft_control%apply_external_potential)
THEN
353 IF (.NOT. dft_control%expot_control%static)
THEN
354 dft_control%eval_external_potential = .true.
359 CALL cp_iterate(logger%iter_info, last=(i_step == max_steps), iter_nr=i_step)
360 rtp%converged = .false.
361 DO i_iter = 1, max_iter
362 IF (i_step == rtp%i_start + 1 .AND. i_iter == 2 .AND. rtp_control%hfx_redistribute) &
367 rtp%energy_new = energy%total
368 IF (rtp%converged)
EXIT
371 IF (rtp%converged)
THEN
373 IF (should_stop)
CALL cp_iterate(logger%iter_info, last=.true., iter_nr=i_step)
375 used_time = time_iter_stop - time_iter_start
376 time_iter_start = time_iter_stop
378 CALL rt_write_input_restart(force_env=force_env, qs_env=qs_env)
379 IF (should_stop)
EXIT
386 IF (.NOT. rtp%converged) &
387 CALL cp_abort(__location__,
"propagation did not converge, "// &
388 "either increase MAX_ITER or use a smaller TIMESTEP")
390 CALL timestop(handle)
392 END SUBROUTINE run_propagation
403 SUBROUTINE rt_write_input_restart(md_env, qs_env, force_env)
404 TYPE(md_environment_type),
OPTIONAL,
POINTER :: md_env
405 TYPE(qs_environment_type),
OPTIONAL,
POINTER :: qs_env
406 TYPE(force_env_type),
POINTER :: force_env
408 REAL(kind=
dp),
DIMENSION(:),
POINTER :: tmp_vals
409 TYPE(dft_control_type),
POINTER :: dft_control
410 TYPE(section_vals_type),
POINTER :: efield_section, motion_section, &
411 root_section, rt_section
413 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
414 root_section => force_env%root_section
419 IF (.NOT.
PRESENT(md_env))
THEN
423 IF (dft_control%apply_vector_potential)
THEN
426 ALLOCATE (tmp_vals(3))
427 tmp_vals = dft_control%efield_fields(1)%efield%vec_pot_initial
429 r_vals_ptr=tmp_vals, &
435 END SUBROUTINE rt_write_input_restart
446 SUBROUTINE rt_initial_guess(qs_env, force_env, rtp_control)
447 TYPE(qs_environment_type),
POINTER :: qs_env
448 TYPE(force_env_type),
POINTER :: force_env
449 TYPE(rtp_control_type),
POINTER :: rtp_control
451 INTEGER :: homo, ispin
452 LOGICAL :: energy_consistency
453 TYPE(cp_fm_type),
POINTER :: mo_coeff
454 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
455 TYPE(dft_control_type),
POINTER :: dft_control
457 NULLIFY (matrix_s, dft_control)
458 CALL get_qs_env(qs_env, dft_control=dft_control)
460 SELECT CASE (rtp_control%initial_wfn)
462 qs_env%sim_time = 0.0_dp
464 energy_consistency = .true.
466 IF (rtp_control%linear_scaling) energy_consistency = .false.
468 consistent_energies=energy_consistency)
469 qs_env%run_rtp = .true.
470 ALLOCATE (qs_env%rtp)
472 IF (dft_control%do_admm)
THEN
474 CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
475 rtp_control%linear_scaling, qs_env%admm_env%mos_aux_fit)
477 CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
478 rtp_control%linear_scaling)
483 IF (.NOT. rtp_control%linear_scaling .OR. rtp_control%initial_wfn ==
use_restart_wfn)
THEN
484 DO ispin = 1,
SIZE(qs_env%mos)
485 CALL get_mo_set(qs_env%mos(ispin), mo_coeff=mo_coeff, homo=homo)
486 IF (.NOT.
ASSOCIATED(mo_coeff))
THEN
488 qs_env%mpools%ao_mo_fm_pools(ispin)%pool, &
489 name=
"qs_env%mo"//trim(adjustl(cp_to_string(ispin))))
494 ALLOCATE (qs_env%rtp)
496 CALL rt_prop_create(qs_env%rtp, qs_env%mos, qs_env%mpools, dft_control, matrix_s(1)%matrix, &
497 rtp_control%linear_scaling, qs_env%admm_env%mos_aux_fit)
500 qs_env%run_rtp = .true.
503 END SUBROUTINE rt_initial_guess
512 SUBROUTINE rt_init_complex_quantities(qs_env, imag_p, imag_ks, imag_h)
513 TYPE(qs_environment_type),
POINTER :: qs_env
514 LOGICAL,
INTENT(in) :: imag_p, imag_ks, imag_h
516 TYPE(dft_control_type),
POINTER :: dft_control
517 TYPE(qs_ks_env_type),
POINTER :: ks_env
518 TYPE(qs_rho_type),
POINTER :: rho
519 TYPE(rt_prop_type),
POINTER :: rtp
521 NULLIFY (ks_env, rho, dft_control)
524 dft_control=dft_control, &
537 IF (.NOT. dft_control%rtp_control%fixed_ions) &
544 END SUBROUTINE rt_init_complex_quantities
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public andermatt2016
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Routines to handle the external control of CP2K.
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
represent a full matrix distributed on many processors
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
subroutine, public cp_rm_iter_level(iteration_info, level_name, n_rlevel_att)
Removes an iteration level.
subroutine, public cp_add_iter_level(iteration_info, level_name, n_rlevel_new)
Adds an iteration level.
all routins needed for a nonperiodic electric field
subroutine, public calculate_ecore_efield(qs_env, calculate_forces)
Computes the force and the energy due to a efield on the cores Note: In the velocity gauge,...
Interface for the force calculations.
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
Define type storing the global information of a run. Keep the amount of stored data small....
Utilities for hfx and admm methods.
subroutine, public hfx_admm_init(qs_env, calculate_forces)
...
Defines the basic variable types.
integer, parameter, public dp
Machine interface based on Fortran 2003 and POSIX.
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
container for various plainwaves related things
Calculation of the core Hamiltonian integral matrix <a|H|b> over Cartesian Gaussian-type functions.
subroutine, public qs_matrix_h_allocate_imag_from_real(qs_env)
(Re-)allocates matrix_h_im from matrix_h
Utility subroutine for qs energy calculation.
subroutine, public qs_energies_init(qs_env, calc_forces)
Refactoring of qs_energies_scf. Driver routine for the initial setup and calculations for a qs energy...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Routines to handle an external electrostatic field The external field can be generic and is provided ...
subroutine, public external_c_potential(qs_env, calculate_forces)
Computes the force and the energy due to the external potential on the cores.
subroutine, public external_e_potential(qs_env)
Computes the external potential on the grid.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public qs_ks_allocate_basics(qs_env, is_complex)
Allocate ks_matrix if necessary, take current overlap matrix as template.
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition and initialisation of the mo data type.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
methods of the rho structure (defined in qs_rho_types)
subroutine, public allocate_rho_ao_imag_from_real(rho, qs_env)
(Re-)allocates rho_ao_im from real part rho_ao
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
...
Routines to apply a delta pulse for RTP and EMD.
subroutine, public apply_delta_pulse(qs_env, rtp, rtp_control)
Interface to call the delta pulse depending on the type of calculation.
Utility functions that are needed for RTP/EMD in combination with HF or hybrid functionals (needs to ...
subroutine, public rtp_hfx_rebuild(qs_env)
rebuilds the structures of P and KS (imaginary) in case S changed
Function related to MO projection in RTP calculations.
subroutine, public init_mo_projection(qs_env, rtp_control)
Initialize the mo projection objects for time dependent run.
Routines for propagating the orbitals.
subroutine, public propagation_step(qs_env, rtp, rtp_control)
performs a single propagation step a(t+Dt)=U(t+Dt,t)*a(0) and calculates the new exponential
Routine for the real time propagation output.
subroutine, public rt_prop_output(qs_env, run_type, delta_iter, used_time)
...
Types and set_get for real time propagation depending on runtype and diagonalization method different...
subroutine, public get_rtp(rtp, exp_H_old, exp_H_new, H_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, S_inv, S_half, S_minus_half, B_mat, C_mat, propagator_matrix, mixing, mixing_factor, S_der, dt, nsteps, SinvH, SinvH_imag, SinvB, admm_mos)
...
subroutine, public rt_prop_create(rtp, mos, mpools, dft_control, template, linear_scaling, mos_aux)
...
subroutine, public rtp_create_sinvh_imag(rtp, nspins)
Initialize SinvH_imag for rtp.
subroutine, public rtp_history_create(rtp, aspc_order)
...
subroutine, public calc_update_rho_sparse(qs_env)
Copies the density matrix back into the qs_envrhorho_ao.
subroutine, public calc_s_derivs(qs_env)
Calculates dS/dR respectily the velocity weighted derivatves only needed for ehrenfest MD.
subroutine, public calc_update_rho(qs_env)
calculates the density from the complex MOs and passes the density to qs_env.
subroutine, public get_restart_wfn(qs_env)
reads the restart file. At the moment only SCF (means only real)
Routines to perform the RTP in the velocity gauge.
subroutine, public velocity_gauge_ks_matrix(qs_env, subtract_nl_term)
...
Routines for the real time propagation.
subroutine, public rt_prop_setup(force_env)
creates rtp_type, gets the initial state, either by reading MO's from file or calling SCF run
Routines for that prepare rtp and EMD.
subroutine, public init_propagators(qs_env)
prepares the initial matrices for the propagators
subroutine, public rt_initialize_rho_from_mos(rtp, mos, mos_old)
Computes the density matrix from the mos Update: Initialized the density matrix from complex mos (for...