comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Fix logicals for RTP/Ehrenfest

Browse source
This commit is contained in:
mattiatj 2023-02-27 13:58:25 +01:00 • committed by GitHub
parent 26d600e12b
commit e1cb72ec33
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
7 changed files with 111 additions and 20 deletions
Download patch file Download diff file Expand all files Collapse all files

18
src/emd/rt_make_propagators.F
View file

@ -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

12
src/emd/rt_propagation_methods.F
View file

@ -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

2
src/emd/rt_propagator_init.F
View file

@ -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

12
src/motion/rt_propagation.F
View file

@ -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

1
src/rt_propagation_types.F
View file

@ -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

85
tests/QS/regtest-rtp-4/H2O-delta-velrepr.inp Normal file
View file

@ -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

1
tests/QS/regtest-rtp-4/TEST_FILES
View file

@ -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