diff --git a/src/emd/rt_make_propagators.F b/src/emd/rt_make_propagators.F index 992d5e8fe4..e174ac65f3 100644 --- a/src/emd/rt_make_propagators.F +++ b/src/emd/rt_make_propagators.F @@ -150,7 +150,7 @@ CONTAINS CALL get_rtp(rtp=rtp, propagator_matrix=propagator_matrix, exp_H_new=exp_H_new, & exp_H_old=exp_H_old, rho_old=rho_old, rho_new=rho_new, rho_next=rho_next) - CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp) + CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp_control, rtp) !I could store these matrices in the type NULLIFY (tmp_re) @@ -235,7 +235,7 @@ CONTAINS DO ispin = 1, nspin re = ispin*2 - 1 im = ispin*2 - IF (.NOT. rtp%track_imag_density) THEN + IF (rtp_control%fixed_ions .AND. .NOT. rtp%propagate_complex_ks) THEN CALL arnoldi(mos_old(re:im), mos_new(re:im), & eps_arnoldi, Him=propagator_matrix_fm(im), & mos_next=mos_next(re:im), narn_old=rtp%narn_old) @@ -296,7 +296,7 @@ CONTAINS CALL get_rtp(rtp=rtp, exp_H_old=exp_H_old, dt=dt) prefac = -0.5_dp*dt CALL dbcsr_scale(exp_H_old(im)%matrix, prefac) - IF (.NOT. rtp%track_imag_density) THEN + IF (rtp_control%fixed_ions .AND. .NOT. rtp%propagate_complex_ks) THEN CALL bch_expansion_imaginary_propagator( & exp_H_old(im)%matrix, rho_next(re)%matrix, rho_next(im)%matrix, & rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp) @@ -310,7 +310,7 @@ CONTAINS END IF CALL dbcsr_copy(rho_new(re)%matrix, rho_next(re)%matrix) CALL dbcsr_copy(rho_new(im)%matrix, rho_next(im)%matrix) - IF (.NOT. rtp%track_imag_density) THEN + IF (rtp_control%fixed_ions .AND. .NOT. rtp%propagate_complex_ks) THEN CALL bch_expansion_imaginary_propagator( & propagator_matrix(im)%matrix, rho_new(re)%matrix, rho_new(im)%matrix, & rtp%filter_eps, rtp%filter_eps_small, rtp_control%eps_exp) @@ -363,7 +363,7 @@ CONTAINS SELECT CASE (rtp_control%mat_exp) CASE (do_taylor) - IF (.NOT. rtp%track_imag_density) THEN + IF (rtp_control%fixed_ions .AND. .NOT. rtp%propagate_complex_ks) THEN CALL taylor_only_imaginary(propagator_fm(re:im), propagator_matrix_fm(im), & rtp%orders(1, ispin), rtp%orders(2, ispin)) ELSE @@ -371,7 +371,7 @@ CONTAINS rtp%orders(1, ispin), rtp%orders(2, ispin)) END IF CASE (do_pade) - IF (.NOT. rtp%track_imag_density) THEN + IF (rtp_control%fixed_ions .AND. .NOT. rtp%propagate_complex_ks) THEN CALL exp_pade_only_imaginary(propagator_fm(re:im), propagator_matrix_fm(im), & rtp%orders(1, ispin), rtp%orders(2, ispin)) ELSE @@ -394,12 +394,14 @@ CONTAINS !> \brief Sparse versions of the matrix exponentials !> \param propagator ... !> \param propagator_matrix ... +!> \param rtp_control ... !> \param rtp ... !> \author Samuel Andermatt (02.14) ! ************************************************************************************************** - SUBROUTINE compute_exponential_sparse(propagator, propagator_matrix, rtp) + SUBROUTINE compute_exponential_sparse(propagator, propagator_matrix, rtp_control, rtp) TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: propagator, propagator_matrix + TYPE(rtp_control_type), POINTER :: rtp_control TYPE(rt_prop_type), POINTER :: rtp CHARACTER(len=*), PARAMETER :: routineN = 'compute_exponential_sparse' @@ -411,7 +413,7 @@ CONTAINS DO ispin = 1, SIZE(propagator)/2 re = 2*ispin - 1 im = 2*ispin - IF (.NOT. rtp%track_imag_density) THEN + IF (rtp_control%fixed_ions .AND. .NOT. rtp%propagate_complex_ks) THEN CALL taylor_only_imaginary_dbcsr(propagator(re:im), propagator_matrix(im)%matrix, & rtp%orders(1, ispin), rtp%orders(2, ispin), rtp%filter_eps) ELSE diff --git a/src/emd/rt_propagation_methods.F b/src/emd/rt_propagation_methods.F index a2e109e308..0c61e8678c 100644 --- a/src/emd/rt_propagation_methods.F +++ b/src/emd/rt_propagation_methods.F @@ -212,7 +212,7 @@ CONTAINS im = 2*i CALL dbcsr_add(ks_mix(i)%matrix, matrix_ks(i)%matrix, 0.0_dp, rtp%mixing_factor) CALL dbcsr_add(ks_mix(i)%matrix, H_last_iter(re)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor) - IF (rtp%do_hfx .OR. rtp_control%velocity_gauge) THEN + IF (rtp%propagate_complex_ks) THEN CALL dbcsr_add(ks_mix_im(i)%matrix, matrix_ks_im(i)%matrix, 0.0_dp, rtp%mixing_factor) CALL dbcsr_add(ks_mix_im(i)%matrix, H_last_iter(im)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor) END IF @@ -222,7 +222,7 @@ CONTAINS re = 2*i - 1 im = 2*i CALL dbcsr_copy(H_last_iter(re)%matrix, ks_mix(i)%matrix) - IF (rtp%do_hfx .OR. rtp_control%velocity_gauge) THEN + IF (rtp%propagate_complex_ks) THEN CALL dbcsr_copy(H_last_iter(im)%matrix, ks_mix_im(i)%matrix) END IF END DO @@ -234,7 +234,7 @@ CONTAINS re = 2*i - 1 im = 2*i CALL dbcsr_copy(H_last_iter(re)%matrix, matrix_ks(i)%matrix) - IF (rtp%do_hfx .OR. rtp_control%velocity_gauge) THEN + IF (rtp%propagate_complex_ks) THEN CALL dbcsr_copy(H_last_iter(im)%matrix, matrix_ks_im(i)%matrix) END IF END DO @@ -426,9 +426,11 @@ CONTAINS CALL dbcsr_copy(SinvH(ispin)%matrix, exp_H(im)%matrix) END IF END DO - IF (rtp%track_imag_density) THEN + + ! in these two cases the argument in the matrix exponentials has a real part + IF (.NOT. rtp_control%fixed_ions .OR. rtp%propagate_complex_ks) THEN CALL get_rtp(rtp=rtp, B_mat=B_mat, SinvB=SinvB) - IF (rtp%do_hfx .OR. rtp_control%velocity_gauge) THEN + IF (rtp%propagate_complex_ks) THEN DO ispin = 1, SIZE(matrix_ks) re = ispin*2 - 1 im = ispin*2 diff --git a/src/emd/rt_propagator_init.F b/src/emd/rt_propagator_init.F index 543ad6749a..912b985ad0 100644 --- a/src/emd/rt_propagator_init.F +++ b/src/emd/rt_propagator_init.F @@ -152,7 +152,7 @@ CONTAINS ELSEIF (rtp_control%mat_exp == do_bch) THEN ELSE IF (rtp%linear_scaling) THEN - CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp) + CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp_control, rtp) ELSE CALL compute_exponential(exp_H_new, propagator_matrix, rtp_control, rtp) END IF diff --git a/src/motion/rt_propagation.F b/src/motion/rt_propagation.F index 0d0b205889..d0a08d6d42 100644 --- a/src/motion/rt_propagation.F +++ b/src/motion/rt_propagation.F @@ -110,7 +110,7 @@ CONTAINS TYPE(force_env_type), POINTER :: force_env INTEGER :: aspc_order - LOGICAL :: ks_is_complex, magnetic, vel_reprs + LOGICAL :: magnetic, vel_reprs TYPE(dft_control_type), POINTER :: dft_control TYPE(global_environment_type), POINTER :: globenv TYPE(qs_energy_type), POINTER :: energy @@ -160,16 +160,16 @@ CONTAINS CALL section_vals_get(hfx_sections, explicit=rtp%do_hfx) CALL section_vals_val_get(print_moments_section, "MAGNETIC", l_val=magnetic) CALL section_vals_val_get(print_moments_section, "VEL_REPRS", l_val=vel_reprs) - rtp%track_imag_density = magnetic .OR. vel_reprs .OR. rtp_control%velocity_gauge & - .OR. rtp%do_hfx .OR. (.NOT. rtp_control%fixed_ions) - ks_is_complex = rtp%do_hfx .OR. rtp_control%velocity_gauge + + rtp%track_imag_density = magnetic .OR. vel_reprs .OR. rtp_control%velocity_gauge .OR. rtp%do_hfx + rtp%propagate_complex_ks = rtp%do_hfx .OR. rtp_control%velocity_gauge CALL rt_init_complex_quantities(qs_env, imag_p=rtp%track_imag_density, & - imag_ks=ks_is_complex, imag_h=rtp_control%velocity_gauge) + imag_ks=rtp%propagate_complex_ks, imag_h=rtp_control%velocity_gauge) ! Hmm, not really like to initialize with the structure of S but I reckon it is ! done everywhere like this - IF (rtp%track_imag_density) CALL rtp_hfx_rebuild(qs_env) + IF (rtp%do_hfx) CALL rtp_hfx_rebuild(qs_env) CALL init_propagation_run(qs_env) IF (.NOT. rtp_control%fixed_ions) THEN diff --git a/src/rt_propagation_types.F b/src/rt_propagation_types.F index 8199cd4e86..350656e8de 100644 --- a/src/rt_propagation_types.F +++ b/src/rt_propagation_types.F @@ -99,6 +99,7 @@ MODULE rt_propagation_types REAL(KIND=dp) :: mixing_factor = 0.0_dp LOGICAL :: mixing = .FALSE. LOGICAL :: do_hfx = .FALSE. + LOGICAL :: propagate_complex_ks = .FALSE. LOGICAL :: track_imag_density = .FALSE. INTEGER, DIMENSION(:, :), ALLOCATABLE :: orders INTEGER :: nsteps = -1 diff --git a/tests/QS/regtest-rtp-4/H2O-delta-velrepr.inp b/tests/QS/regtest-rtp-4/H2O-delta-velrepr.inp new file mode 100644 index 0000000000..057deacf66 --- /dev/null +++ b/tests/QS/regtest-rtp-4/H2O-delta-velrepr.inp @@ -0,0 +1,85 @@ +&GLOBAL + PROJECT H2O_rtp + RUN_TYPE RT_PROPAGATION + PRINT_LEVEL LOW +&END GLOBAL +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + POTENTIAL_FILE_NAME POTENTIAL + !RESTART_FILE_NAME H2O_rtp-RESTART.rtpwfn + &MGRID + CUTOFF 200 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + &END QS + &SCF + MAX_SCF 20 + SCF_GUESS ATOMIC + &END SCF + &XC + &XC_FUNCTIONAL PADE + &END XC_FUNCTIONAL + &END XC + &REAL_TIME_PROPAGATION + APPLY_DELTA_PULSE T + DELTA_PULSE_DIRECTION 1 0 0 + DELTA_PULSE_SCALE 0.001 + PERIODIC .TRUE. + MAX_ITER 50 + MAT_EXP PADE + EXP_ACCURACY 1.0E-10 + EPS_ITER 1.0E-9 + PROPAGATOR ETRS + INITIAL_WFN SCF_WFN ! RT_RESTART + &END + &PRINT + &MOMENTS low + PERIODIC F + MAGNETIC T + MAX_MOMENT 2 + VEL_REPRS T + COM_NL T + REFERENCE COM + &END MOMENTS + &END + &END DFT + &SUBSYS + &CELL + ABC 7 7 7 + PERIODIC NONE + &END CELL + &TOPOLOGY + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &COORD + O 0.000000 0.000000 -0.065587 H2O + H 0.000000 -0.757136 0.520545 H2O + H 0.000000 0.757136 0.520545 H2O + &END COORD + &KIND H + BASIS_SET DZVP-MOLOPT-GTH + POTENTIAL GTH-PADE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-MOLOPT-GTH + POTENTIAL GTH-PADE-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&MOTION + &MD + STEPS 2 + TIMESTEP [au_t] 0.1 + &END +&END MOTION +!&EXT_RESTART +! RESTART_FILE_NAME H2O_rtp-1.restart +!&END + diff --git a/tests/QS/regtest-rtp-4/TEST_FILES b/tests/QS/regtest-rtp-4/TEST_FILES index 5a4c724f78..eec53a6525 100644 --- a/tests/QS/regtest-rtp-4/TEST_FILES +++ b/tests/QS/regtest-rtp-4/TEST_FILES @@ -8,4 +8,5 @@ H2O-delta-mag.inp 1 1e-13 -17.17816619814588 H2O-delta-lenrep.inp 1 1e-13 -17.17816511743962 H2O-vg.inp 1 1e-13 -17.17816522919518 +H2O-delta-velrepr.inp 1 1e-13 -17.17816551662581 #EOF