Skip to content

Commit

Permalink
Add linear density delta kick
Browse files Browse the repository at this point in the history
  • Loading branch information
@glb96
glb96 committed Feb 3, 2023
1 parent 3114f69 commit 427e198
Show file tree
Hide file tree
Showing 36 changed files with 2,345 additions and 147 deletions.
90 changes: 80 additions & 10 deletions src/emd/rt_delta_pulse.F
View file
Expand Up @@ -63,6 +63,7 @@ MODULE rt_delta_pulse
USE moments_utils, ONLY: get_reference_point
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_types, ONLY: particle_type
USE qs_dftb_matrices, ONLY: build_dftb_overlap
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_kind_types, ONLY: qs_kind_type
Expand All @@ -72,7 +73,9 @@ MODULE rt_delta_pulse
build_local_magmom_matrix,&
build_local_moment_matrix
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE rt_propagation_types, ONLY: rt_prop_type
USE rt_propagation_types, ONLY: get_rtp,&
rt_prop_create_mos,&
rt_prop_type
#include "../base/base_uses.f90"

IMPLICIT NONE
Expand All @@ -81,12 +84,79 @@ MODULE rt_delta_pulse

CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_delta_pulse'

PUBLIC :: apply_delta_pulse_periodic, &
apply_delta_pulse, &
apply_delta_pulse_mag
PUBLIC :: apply_delta_pulse

CONTAINS

! **************************************************************************************************
!> \brief Interface to call the delta pulse depending on the type of calculation.
!> \param qs_env ...
!> \param rtp ...
!> \param rtp_control ...
!> \author Update: Guillaume Le Breton (2023.01)
! **************************************************************************************************

SUBROUTINE apply_delta_pulse(qs_env, rtp, rtp_control)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(rt_prop_type), POINTER :: rtp
TYPE(rtp_control_type), POINTER :: rtp_control

LOGICAL :: my_apply_pulse, periodic_cell
TYPE(cell_type), POINTER :: cell
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_old
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos

CALL get_qs_env(qs_env, &
cell=cell, &
dft_control=dft_control, &
matrix_s=matrix_s)
periodic_cell = ANY(cell%perd > 0)
my_apply_pulse = .TRUE.
CALL get_qs_env(qs_env, mos=mos)
IF (rtp%linear_scaling) THEN
IF (.NOT. ASSOCIATED(mos)) THEN
CALL cp_warn(__LOCATION__, "Delta Pulse not implemented for Linear-Scaling based ground "// &
"state calculation. If you want to perform a Linear-Scaling RTP from a "// &
"Linear-Scaling GS calculation you can do the following: (i) LSCF froms "// &
"scratch, (ii) MO-based SCF (for 1 SCF loop for instance) with the LSCF "// &
"result as a restart and (iii) linear scaling RTP + delta kick (for 1 "// &
"SCF loop for instance).")
my_apply_pulse = .FALSE.
ELSE
! create temporary mos_old and mos_new to use delta kick routine designed for MOs-based RTP
CALL rt_prop_create_mos(rtp, mos, qs_env%mpools, dft_control, &
init_mos_old=.TRUE., init_mos_new=.TRUE., &
init_mos_next=.FALSE., init_mos_admn=.FALSE.)
END IF
END IF

IF (my_apply_pulse) THEN
CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
IF (rtp_control%apply_delta_pulse) THEN
IF (dft_control%qs_control%dftb) &
CALL build_dftb_overlap(qs_env, 1, matrix_s)
IF (rtp_control%periodic) THEN
CALL apply_delta_pulse_electric_periodic(qs_env, mos_old, mos_new)
ELSE
IF (periodic_cell) THEN
CPWARN("This application of the delta pulse is not compatible with PBC!")
END IF
CALL apply_delta_pulse_electric(qs_env, mos_old, mos_new)
END IF
ELSE IF (rtp_control%apply_delta_pulse_mag) THEN
IF (periodic_cell) THEN
CPWARN("This application of the delta pulse is not compatible with PBC!")
END IF
CALL apply_delta_pulse_mag(qs_env, mos_old, mos_new)
ELSE
CPABORT("Code error: this case should not happen!")
END IF
END IF

END SUBROUTINE apply_delta_pulse

! **************************************************************************************************
!> \brief uses perturbation theory to get the proper initial conditions
!> The len_rep option is NOT compatible with periodic boundary conditions!
Expand All @@ -96,11 +166,11 @@ MODULE rt_delta_pulse
!> \author Joost & Martin (2011)
! **************************************************************************************************

SUBROUTINE apply_delta_pulse_periodic(qs_env, mos_old, mos_new)
SUBROUTINE apply_delta_pulse_electric_periodic(qs_env, mos_old, mos_new)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_old, mos_new

CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_periodic'
CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_electric_periodic'

INTEGER :: handle, icol, idir, irow, ispin, nao, &
ncol_local, nmo, nrow_local, nvirt, &
Expand Down Expand Up @@ -302,7 +372,7 @@ SUBROUTINE apply_delta_pulse_periodic(qs_env, mos_old, mos_new)

CALL timestop(handle)

END SUBROUTINE apply_delta_pulse_periodic
END SUBROUTINE apply_delta_pulse_electric_periodic

! **************************************************************************************************
!> \brief applies exp(ikr) to the wavefunction.... stored in mos_old...
Expand All @@ -312,11 +382,11 @@ END SUBROUTINE apply_delta_pulse_periodic
!> \author Joost & Martin (2011)
! **************************************************************************************************

SUBROUTINE apply_delta_pulse(qs_env, mos_old, mos_new)
SUBROUTINE apply_delta_pulse_electric(qs_env, mos_old, mos_new)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_old, mos_new

CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse'
CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_electric'

INTEGER :: handle, i, nao, nmo
REAL(KIND=dp), DIMENSION(3) :: kvec
Expand Down Expand Up @@ -382,7 +452,7 @@ SUBROUTINE apply_delta_pulse(qs_env, mos_old, mos_new)

CALL timestop(handle)

END SUBROUTINE apply_delta_pulse
END SUBROUTINE apply_delta_pulse_electric

! **************************************************************************************************
!> \brief apply magnetic delta pulse to linear order
Expand Down
33 changes: 24 additions & 9 deletions src/emd/rt_propagation_utils.F
View file
Expand Up @@ -266,6 +266,7 @@ SUBROUTINE get_restart_wfn(qs_env)
unit_nr
REAL(KIND=dp) :: alpha, cs_pos
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_fm_type) :: mos_occ
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_old
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_para_env_type), POINTER :: para_env
Expand Down Expand Up @@ -315,15 +316,29 @@ SUBROUTINE get_restart_wfn(qs_env)
re = 2*ispin - 1
im = 2*ispin
CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, ncol_global=ncol)
alpha = 1.0_dp
IF (SIZE(mo_array) == 1) alpha = 2*alpha
CALL cp_dbcsr_plus_fm_fm_t( &
sparse_matrix=rho_old(re)%matrix, &
matrix_v=mo_array(ispin)%mo_coeff, matrix_g=mo_array(ispin)%mo_coeff, ncol=ncol, &
keep_sparsity=.FALSE., alpha=alpha)
END DO
DO i = 1, nspin
CALL dbcsr_copy(rho_new(i)%matrix, rho_old(i)%matrix)
CALL cp_fm_create(mos_occ, &
matrix_struct=mo_array(ispin)%mo_coeff%matrix_struct, &
name="mos_occ")
CALL cp_fm_to_fm(mo_array(ispin)%mo_coeff, mos_occ)
IF (mo_array(ispin)%uniform_occupation) THEN
alpha = 3.0_dp - REAL(nspin, dp)
CALL cp_fm_column_scale(mos_occ, mo_array(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
ELSE
alpha = 1.0_dp
CALL cp_fm_column_scale(mos_occ, mo_array(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mo_array(ispin)%mo_coeff, &
matrix_g=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
END IF
CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
CALL cp_fm_release(mos_occ)
END DO
CALL calc_update_rho_sparse(qs_env)
ELSE
Expand Down
138 changes: 108 additions & 30 deletions src/emd/rt_propagator_init.F
View file
Expand Up @@ -54,6 +54,7 @@ MODULE rt_propagator_init
put_data_to_history,&
s_matrices_create
USE rt_propagation_types, ONLY: get_rtp,&
rt_prop_release_mos,&
rt_prop_type
#include "../base/base_uses.f90"

Expand Down Expand Up @@ -406,54 +407,131 @@ END SUBROUTINE backtransform_matrix

! **************************************************************************************************
!> \brief Computes the density matrix from the mos
!> Update: Initialized the density matrix from complex mos (for
!> instance after delta kick)
!> \param rtp ...
!> \param mos ...
!> \param mos_old ...
!> \author Samuel Andermatt (08.15)
!> \author Guillaume Le Breton (01.23)
! **************************************************************************************************

SUBROUTINE rt_initialize_rho_from_mos(rtp, mos)
SUBROUTINE rt_initialize_rho_from_mos(rtp, mos, mos_old)

TYPE(rt_prop_type), POINTER :: rtp
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mos_old

INTEGER :: ispin, ncol, re
INTEGER :: im, ispin, ncol, re
REAL(KIND=dp) :: alpha
TYPE(cp_fm_type) :: mos_occ
TYPE(cp_fm_type) :: mos_occ, mos_occ_im
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new, rho_old

CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)

DO ispin = 1, SIZE(mos)
re = 2*ispin - 1
CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)

CALL cp_fm_create(mos_occ, &
matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
name="mos_occ")
CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
IF (mos(ispin)%uniform_occupation) THEN
alpha = 3.0_dp - REAL(SIZE(mos), dp)
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
ELSE
alpha = 1.0_dp
IF (PRESENT(mos_old)) THEN
! Used the mos from delta kick. Initialize both real and im part
DO ispin = 1, SIZE(mos_old)/2
re = 2*ispin - 1; im = 2*ispin
CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
CALL cp_fm_create(mos_occ, &
matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
name="mos_occ")
!Real part of rho
CALL cp_fm_to_fm(mos_old(re), mos_occ)
IF (mos(ispin)%uniform_occupation) THEN
alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
ELSE
alpha = 1.0_dp
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mos_old(re), &
matrix_g=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
END IF

! Add complex part of MOs, i*i=-1
CALL cp_fm_to_fm(mos_old(im), mos_occ)
IF (mos(ispin)%uniform_occupation) THEN
alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
ELSE
alpha = 1.0_dp
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mos_old(im), &
matrix_g=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
END IF
CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)

! Imaginary part of rho
CALL cp_fm_create(mos_occ_im, &
matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
name="mos_occ_im")
alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
CALL cp_fm_to_fm(mos_old(re), mos_occ)
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mos(ispin)%mo_coeff, &
CALL cp_fm_to_fm(mos_old(im), mos_occ_im)
CALL cp_fm_column_scale(mos_occ_im, mos(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(im)%matrix, &
matrix_v=mos_occ_im, &
matrix_g=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
END IF
alpha=2.0_dp*alpha, &
symmetry_mode=-1, keep_sparsity=.FALSE.)

CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
CALL cp_fm_release(mos_occ)

END DO
CALL dbcsr_filter(rho_old(im)%matrix, rtp%filter_eps)
CALL dbcsr_copy(rho_new(im)%matrix, rho_old(im)%matrix)
CALL cp_fm_release(mos_occ_im)
CALL cp_fm_release(mos_occ)
END DO
! Release the mos used to apply the delta kick, no longer required
CALL rt_prop_release_mos(rtp)
ELSE
DO ispin = 1, SIZE(mos)
re = 2*ispin - 1
CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)

CALL cp_fm_create(mos_occ, &
matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
name="mos_occ")
CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
IF (mos(ispin)%uniform_occupation) THEN
alpha = 3.0_dp - REAL(SIZE(mos), dp)
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
ELSE
alpha = 1.0_dp
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
matrix_v=mos(ispin)%mo_coeff, &
matrix_g=mos_occ, &
ncol=ncol, &
alpha=alpha, keep_sparsity=.FALSE.)
END IF
CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
CALL cp_fm_release(mos_occ)
END DO
END IF

END SUBROUTINE rt_initialize_rho_from_mos

Expand Down

0 comments on commit 427e198

Please sign in to comment.