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

RTP gauges (#2343)

Browse source 
* Add bare routines for velocity gauge rtp

* RTP-VG for all-electron (no pseudo potentials yet)

* add routine for nl term

* allocate, deallocate exp. integrals

* add .fortls to .gitignore

* calculate cos and sine sap integrals for exponential

* OMP parallelization of integrals

* assemble gauge term, something still wrong

* fixed indexing for integrals

* add nl gauge term to ks matrix

* omp parallelization

* only 1 dbcsr_addition for quadratic term

* only 1 dbcsr_addition for quadratic term

* fix some inconsistency

* fix out of bounds error

* debug options for the individual terms

* change reference point for nl term to pseudized atom

* adapt some factors

* new banner, correct sign, do not use non-gauge corrected term by default

* add nl gauge term to matrix_h_kp instead of matrix_ks

* anti-symmetric initialization of imaginary part of ks matrix also during mixing

* make matrix_ks_im a kpoint_transitional_type

* make rho_ao_im a kpoint_transitional_type and adapt allocations for rho_ao_im and matrix_ks_im

* add matrix_h_im as kpoint_transitional_type and add to matrix_ks_im if necessary; subtract non-gauge transformed term if necessary

* initiate eps_ppnl

* simplify doing the gauge transform of nl term or not

* change sign of gauge transformation

* add delta pulse with DFPT in length representation

* print expectation values of nl commutators separately

* change format specifiers for multipole moments to D16.8

* make pretty

* rebase, make pretty, calculate mag mom as anti-symmetric

* make pretty

* release matrix_ks_im as k-point transitional type

* re-calculate and subtract core_ppnl only if necessary, make gauge transformation of non-local part default

* refactor delta-pulse code and add warnings for PBC

* remove debug input parameters, adapt .gitignore

* fix deallocation

* add linear magnetic delta pulse

* correct sign of magnetic dipole moment contributino

* use distributed refpoint for [r,V_nl] in periodic systems

* add regetests for lenrep/magnetic delta pulse and velocity gauge

* add counter for rtp step in output and keep maxsteps consistent after restarts

* fixes after rebase

* add COM_NL .FALSE. to delta-pulse regtests for consistency with previous result

* add references

* change format specifier to exponential notation

* fix allocation for non-local moment integrals

* make pretty
This commit is contained in:
mattiatj 2022-10-11 16:05:52 +02:00 • committed by GitHub
parent 5277e49418
commit e64daf41ce
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
25 changed files with 2269 additions and 384 deletions
Download patch file Download diff file Expand all files Collapse all files

34
src/common/bibliography.F
View file

@ -89,7 +89,7 @@ MODULE bibliography
Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, Lass2018, cp2kqs2020, &
Behler2007, Behler2011, Schran2020a, Schran2020b, &
Rycroft2009, Thomas2015, Brehm2018, Brehm2020, Shigeta2001, Heinecke2016, &
Brehm2021, Bussy2021a, Bussy2021b, Ditler2021
Brehm2021, Bussy2021a, Bussy2021b, Ditler2021, Mattiat2019, Mattiat2022
CONTAINS
@ -4719,6 +4719,38 @@ CONTAINS
"ER"), &
DOI="10.1063/5.0041056")
CALL add_reference(key=Mattiat2019, ISI_record=s2a( &
"TY JOUR", &
"PT J", &
"AU Mattiat, Johann", &
" Luber, Sandra", &
"TI Vibrational (resonance) Raman optical activity with real time time dependent", &
" density functional theory", &
"SO The Journal of Chemical Physics", &
"PY 2019", &
"VL 151", &
"AR 234110", &
"DI 10.1063/1.5132294", &
"ER"), &
DOI="10.1063/1.5132294")
CALL add_reference(key=Mattiat2022, ISI_record=s2a( &
"TY JOUR", &
"PT J", &
"AU Mattiat, Johann", &
" Luber, Sandra", &
"TI Comparison of Length, Velocity, and Symmetric Gauges for the Calculation of", &
" Absorption and Electric Circular Dichroism Spectra with Real-Time Time-Dependent", &
" Density Functional Theory", &
"SO The Journal of Chemical Theory and Computation", &
"PY 2022", &
"VL 18", &
"BP 5513", &
"EP 5526", &
"DI 10.1021/acs.jctc.2c00644", &
"ER"), &
DOI="10.1021/acs.jctc.2c00644")
END SUBROUTINE add_all_references
END MODULE bibliography

261
src/commutator_rpnl.F
View file

@ -59,7 +59,7 @@ MODULE commutator_rpnl
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'commutator_rpnl'
PUBLIC :: build_com_rpnl, build_com_mom_nl
PUBLIC :: build_com_rpnl, build_com_mom_nl, build_com_nl_mag
CONTAINS
@ -415,14 +415,14 @@ CONTAINS
!> \param sap_ppnl ...
!> \param eps_ppnl ...
!> \param particle_set ...
!> \param cell ...
!> \param matrix_rv ...
!> \param matrix_rxrv ...
!> \param matrix_rrv ...
!> \param ref_point ...
!> \param cell ...
! **************************************************************************************************
SUBROUTINE build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_rv, matrix_rxrv, &
matrix_rrv, ref_point, cell)
SUBROUTINE build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, &
matrix_rrv, ref_point)
TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
POINTER :: qs_kind_set
@ -431,10 +431,10 @@ CONTAINS
REAL(KIND=dp), INTENT(IN) :: eps_ppnl
TYPE(particle_type), DIMENSION(:), INTENT(IN), &
POINTER :: particle_set
TYPE(cell_type), INTENT(IN), POINTER :: cell
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
OPTIONAL, POINTER :: matrix_rv, matrix_rxrv, matrix_rrv
REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: ref_point
TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_com_mom_nl'
@ -445,7 +445,7 @@ CONTAINS
INTEGER, DIMENSION(3) :: cell_b
LOGICAL :: asso_rrv, asso_rv, asso_rxrv, found, go, &
my_ref, my_rrv, my_rv, my_rxrv, &
ppnl_present
periodic, ppnl_present
REAL(KIND=dp), DIMENSION(3) :: rab, rf
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
TYPE(alist_type), POINTER :: alist_ac, alist_bc
@ -484,11 +484,17 @@ CONTAINS
order = 1
END IF
periodic = ANY(cell%perd > 0)
my_ref = .FALSE.
IF (PRESENT(ref_point)) THEN
CPASSERT(PRESENT(cell)) ! need cell as well if refpoint is provided
rf = ref_point
my_ref = .TRUE.
IF (.NOT. periodic) THEN
rf = ref_point
my_ref = .TRUE.
ELSE ! use my_ref = False in periodic case, corresponds to distributed ref point
IF (order .GT. 1) THEN
CPWARN("Not clear how to define reference point for order > 1 in periodic cells.")
END IF
END IF
END IF
nkind = SIZE(qs_kind_set)
@ -924,4 +930,241 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE build_com_mom_nl
! **************************************************************************************************
!> \brief calculate \sum_R_ps (R_ps - R_nu) x [V_nl, r] summing over all pseudized atoms R
!> \param qs_kind_set ...
!> \param sab_all ...
!> \param sap_ppnl ...
!> \param eps_ppnl ...
!> \param particle_set ...
!> \param matrix_mag_nl ...
!> \param refpoint ...
!> \param cell ...
! **************************************************************************************************
SUBROUTINE build_com_nl_mag(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_mag_nl, refpoint, cell)
TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
POINTER :: qs_kind_set
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
INTENT(IN), POINTER :: sab_all, sap_ppnl
REAL(KIND=dp), INTENT(IN) :: eps_ppnl
TYPE(particle_type), DIMENSION(:), INTENT(IN), &
POINTER :: particle_set
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
POINTER :: matrix_mag_nl
REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: refpoint
TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_com_nl_mag'
INTEGER :: handle, iab, iac, iatom, ibc, icol, &
ikind, ind, irow, jatom, jkind, kac, &
kbc, kkind, na, natom, nb, nkind, np, &
order, slot
INTEGER, DIMENSION(3) :: cell_b
LOGICAL :: found, go, my_ref, ppnl_present
REAL(KIND=dp), DIMENSION(3) :: r_b, r_ps, rab
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
TYPE(alist_type), POINTER :: alist_ac, alist_bc
TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_mag
TYPE(gto_basis_set_p_type), ALLOCATABLE, &
DIMENSION(:) :: basis_set
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
!$ INTEGER(kind=omp_lock_kind), &
!$ ALLOCATABLE, DIMENSION(:) :: locks
!$ INTEGER :: lock_num, hash
!$ INTEGER, PARAMETER :: nlock = 501
ppnl_present = ASSOCIATED(sap_ppnl)
IF (.NOT. ppnl_present) RETURN
CALL timeset(routineN, handle)
my_ref = .FALSE.
IF (PRESENT(refpoint)) THEN
my_ref = .TRUE.
CPASSERT(PRESENT(cell))
END IF
natom = SIZE(particle_set)
nkind = SIZE(qs_kind_set)
! allocate integral storage
NULLIFY (sap_int)
ALLOCATE (sap_int(nkind*nkind))
DO ind = 1, nkind*nkind
NULLIFY (sap_int(ind)%alist, sap_int(ind)%asort, sap_int(ind)%aindex)
sap_int(ind)%nalist = 0
END DO
! build integrals over GTO + projector functions, refpoint actually
order = 1 ! only need first moments (x, y, z)
! refpoint actually does not matter in this case, i. e. (order = 1 .and. commutator)
IF (my_ref) THEN
CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE., refpoint=refpoint, &
particle_set=particle_set, cell=cell)
ELSE
CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE.)
END IF
CALL sap_sort(sap_int)
! get access to basis sets
ALLOCATE (basis_set(nkind))
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
IF (ASSOCIATED(orb_basis_set)) THEN
basis_set(ikind)%gto_basis_set => orb_basis_set
ELSE
NULLIFY (basis_set(ikind)%gto_basis_set)
END IF
END DO
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP SHARED (basis_set, matrix_mag_nl, sap_int, natom, nkind, eps_ppnl, locks, sab_all, &
!$OMP particle_set, my_ref, refpoint) &
!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
!$OMP iab, irow, icol, blocks_mag, r_ps, r_b, go, hash, &
!$OMP found, iac, ibc, alist_ac, alist_bc, acint, bcint, &
!$OMP achint, bchint, na, np, nb, kkind, kac, kbc)
!$OMP SINGLE
!$ ALLOCATE (locks(nlock))
!$OMP END SINGLE
!$OMP DO
!$ DO lock_num = 1, nlock
!$ call omp_init_lock(locks(lock_num))
!$ END DO
!$OMP END DO
!$OMP DO SCHEDULE(GUIDED)
DO slot = 1, sab_all(1)%nl_size
! get indices
ikind = sab_all(1)%nlist_task(slot)%ikind
jkind = sab_all(1)%nlist_task(slot)%jkind
iatom = sab_all(1)%nlist_task(slot)%iatom
jatom = sab_all(1)%nlist_task(slot)%jatom
cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
iab = ikind + nkind*(jkind - 1)
IF (iatom <= jatom) THEN
irow = iatom
icol = jatom
ELSE
irow = jatom
icol = iatom
END IF
! get blocks
ALLOCATE (blocks_mag(3))
DO ind = 1, 3
CALL dbcsr_get_block_p(matrix_mag_nl(ind)%matrix, irow, icol, blocks_mag(ind)%block, found)
END DO
go = (ASSOCIATED(blocks_mag(1)%block) .AND. ASSOCIATED(blocks_mag(2)%block) .AND. ASSOCIATED(blocks_mag(3)%block))
IF (go) THEN
DO kkind = 1, nkind
iac = ikind + nkind*(kkind - 1)
ibc = jkind + nkind*(kkind - 1)
IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
CALL get_alist(sap_int(iac), alist_ac, iatom)
CALL get_alist(sap_int(ibc), alist_bc, jatom)
IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
DO kac = 1, alist_ac%nclist
DO kbc = 1, alist_bc%nclist
IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
acint => alist_ac%clist(kac)%acint
bcint => alist_bc%clist(kbc)%acint
achint => alist_ac%clist(kac)%achint
bchint => alist_bc%clist(kbc)%achint
na = SIZE(acint, 1)
np = SIZE(acint, 2)
nb = SIZE(bcint, 1)
! Position of the pseudized atom
r_ps = particle_set(alist_ac%clist(kac)%catom)%r
r_b = refpoint
!$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
!$ CALL omp_set_lock(locks(hash))
! assemble integrals
IF (iatom <= jatom) THEN
blocks_mag(1)%block(1:na, 1:nb) = blocks_mag(1)%block(1:na, 1:nb) + &
(r_ps(2) - r_b(2))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 4))) - &
MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 1)))) & ! R_y [V_nl, z]
- (r_ps(3) - r_b(3))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 3))) - &
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 1)))) ! - R_z [V_nl, y]
blocks_mag(2)%block(1:na, 1:nb) = blocks_mag(2)%block(1:na, 1:nb) + &
(r_ps(3) - r_b(3))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 2))) - &
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1)))) & ! R_z [V_nl, x]
- (r_ps(1) - r_b(1))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 4))) - &
MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 1)))) ! - R_x [V_nl, z]
blocks_mag(3)%block(1:na, 1:nb) = blocks_mag(3)%block(1:na, 1:nb) + &
(r_ps(1) - r_b(1))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 3))) - &
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 1)))) & ! R_x [V_nl, y]
- (r_ps(2) - r_b(2))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 2))) - &
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1)))) ! - R_y [V_nl, x]
ELSE
blocks_mag(1)%block(1:nb, 1:na) = blocks_mag(1)%block(1:nb, 1:na) + &
(r_ps(2) - r_b(2))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 4))) - &
MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 1)))) & ! R_y [V_nl, z]
- (r_ps(3) - r_b(3))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 3))) - &
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 1)))) ! - R_z [V_nl, y]
blocks_mag(2)%block(1:nb, 1:na) = blocks_mag(2)%block(1:nb, 1:na) + &
(r_ps(3) - r_b(3))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 2))) - &
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 1)))) & ! R_z [V_nl, x]
- (r_ps(1) - r_b(1))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 4))) - &
MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 1)))) ! - R_x [V_nl, z]
blocks_mag(3)%block(1:nb, 1:na) = blocks_mag(3)%block(1:nb, 1:na) + &
(r_ps(1) - r_b(1))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 3))) - &
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 1)))) & ! R_x [V_nl, y]
- (r_ps(2) - r_b(2))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 2))) - &
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 1)))) ! - R_y [V_nl, x]
END IF
!$ CALL omp_unset_lock(locks(hash))
EXIT ! We have found a match and there can be only one single match
END IF
END DO
END DO
END DO
END IF
DO ind = 1, 3
NULLIFY (blocks_mag(ind)%block)
END DO
DEALLOCATE (blocks_mag)
END DO
!$OMP DO
!$ DO lock_num = 1, nlock
!$ call omp_destroy_lock(locks(lock_num))
!$ END DO
!$OMP END DO
!$OMP SINGLE
!$ DEALLOCATE (locks)
!$OMP END SINGLE NOWAIT
!$OMP END PARALLEL
DEALLOCATE (basis_set)
CALL release_sap_int(sap_int)
CALL timestop(handle)
END SUBROUTINE build_com_nl_mag
END MODULE commutator_rpnl

3
src/cp_control_types.F
View file

@ -56,6 +56,7 @@ MODULE cp_control_types
INTEGER :: aspc_order
INTEGER :: sc_check_start
LOGICAL :: apply_delta_pulse
LOGICAL :: apply_delta_pulse_mag
LOGICAL :: periodic
LOGICAL :: linear_scaling
LOGICAL :: write_restart
@ -64,6 +65,8 @@ MODULE cp_control_types
REAL(dp) :: mcweeny_eps
INTEGER, DIMENSION(3) :: delta_pulse_direction
REAL(KIND=dp) :: delta_pulse_scale
LOGICAL :: velocity_gauge
LOGICAL :: nl_gauge_transform
END TYPE rtp_control_type
! **************************************************************************************************
! \brief Control parameters for DFTB calculations

6
src/cp_control_utils.F
View file

@ -2189,6 +2189,12 @@ CONTAINS
l_val=dft_control%rtp_control%hfx_redistribute)
CALL section_vals_val_get(rtp_section, "APPLY_DELTA_PULSE", &
l_val=dft_control%rtp_control%apply_delta_pulse)
CALL section_vals_val_get(rtp_section, "APPLY_DELTA_PULSE_MAG", &
l_val=dft_control%rtp_control%apply_delta_pulse_mag)
CALL section_vals_val_get(rtp_section, "VELOCITY_GAUGE", &
l_val=dft_control%rtp_control%velocity_gauge)
CALL section_vals_val_get(rtp_section, "VG_COM_NL", &
l_val=dft_control%rtp_control%nl_gauge_transform)
CALL section_vals_val_get(rtp_section, "PERIODIC", &
l_val=dft_control%rtp_control%periodic)
CALL section_vals_val_get(rtp_section, "DENSITY_PROPAGATION", &

610
src/emd/rt_delta_pulse.F
View file

@ -10,22 +10,28 @@
! **************************************************************************************************
MODULE rt_delta_pulse
USE bibliography, ONLY: Mattiat2019,&
Mattiat2022,&
cite_reference
USE cell_types, ONLY: cell_type
USE commutator_rpnl, ONLY: build_com_mom_nl
USE commutator_rpnl, ONLY: build_com_mom_nl,&
build_com_nl_mag
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale
USE cp_cfm_diag, ONLY: cp_cfm_heevd
USE cp_cfm_types, ONLY: cp_cfm_create,&
cp_cfm_release,&
cp_cfm_to_cfm,&
cp_cfm_type
USE cp_control_types, ONLY: dft_control_type
USE cp_control_types, ONLY: dft_control_type,&
rtp_control_type
USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
copy_fm_to_dbcsr,&
cp_dbcsr_sm_fm_multiply,&
dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
cp_fm_triangular_multiply,&
cp_fm_upper_to_full
USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
cp_fm_cholesky_invert,&
@ -42,13 +48,21 @@ MODULE rt_delta_pulse
cp_fm_set_all,&
cp_fm_to_fm,&
cp_fm_type
USE cp_fm_vect, ONLY: cp_fm_vect_dealloc
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: &
dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_filter, dbcsr_init_p, &
dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric
USE input_section_types, ONLY: section_get_lval,&
section_vals_type
dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_info, &
dbcsr_init_p, dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
dbcsr_type_symmetric
USE input_section_types, ONLY: section_get_ival,&
section_get_lval,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE mathconstants, ONLY: twopi
USE mathconstants, ONLY: one,&
twopi,&
zero
USE moments_utils, ONLY: get_reference_point
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_types, ONLY: particle_type
USE qs_environment_types, ONLY: get_qs_env,&
@ -56,10 +70,11 @@ MODULE rt_delta_pulse
USE qs_kind_types, ONLY: qs_kind_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_p_type
USE qs_moments, ONLY: build_berry_moment_matrix
USE qs_moments, ONLY: build_berry_moment_matrix,&
build_local_magmom_matrix,&
build_local_moment_matrix
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE rt_propagation_types, ONLY: get_rtp,&
rt_prop_type
USE rt_propagation_types, ONLY: rt_prop_type
#include "../base/base_uses.f90"
IMPLICIT NONE
@ -69,12 +84,14 @@ MODULE rt_delta_pulse
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_delta_pulse'
PUBLIC :: apply_delta_pulse_periodic, &
apply_delta_pulse
apply_delta_pulse, &
apply_delta_pulse_mag
CONTAINS
! **************************************************************************************************
!> \brief uses perturbation theory to get the proper initial conditions
!> The len_rep option is NOT compatible with periodic boundary conditions!
!> \param qs_env ...
!> \param mos_old ...
!> \param mos_new ...
@ -87,25 +104,21 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_periodic'
COMPLEX(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_sqrt
INTEGER :: handle, icol, idir, irow, ispin, nao, &
ncol_local, nmo, nrow_global, &
nrow_local, nvirt
ncol_local, nmo, nrow_local, nvirt, &
reference
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
LOGICAL :: com_nl
LOGICAL :: com_nl, len_rep, periodic_cell
REAL(KIND=dp) :: eps_ppnl, factor
REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
POINTER :: local_data
REAL(KIND=dp), DIMENSION(3) :: kvec
REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
REAL(KIND=dp), DIMENSION(3) :: kvec, rcc
REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, ref_point
TYPE(cell_type), POINTER :: cell
TYPE(cp_cfm_type) :: oo_c, oo_v, oo_vt
TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_tmp
TYPE(cp_fm_type) :: eigenvectors, mat_ks, mat_tmp, momentum, &
oo_1, oo_2, S_chol, S_inv_fm, tmpS, &
virtuals
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_rv, matrix_s
TYPE(dbcsr_type), POINTER :: S_inv
S_chol, virtuals
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_r, matrix_rv, matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
@ -113,11 +126,12 @@ CONTAINS
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(rt_prop_type), POINTER :: rtp
TYPE(rtp_control_type), POINTER :: rtp_control
TYPE(section_vals_type), POINTER :: input
CALL timeset(routineN, handle)
NULLIFY (cell, mos, rtp, matrix_s, matrix_ks, input, dft_control, particle_set)
NULLIFY (cell, mos, rtp, matrix_s, matrix_ks, input, dft_control, particle_set, fm_struct)
! we need the overlap and ks matrix for a full diagionalization
CALL get_qs_env(qs_env, &
cell=cell, &
@ -128,21 +142,17 @@ CONTAINS
dft_control=dft_control, &
input=input, &
particle_set=particle_set)
com_nl = section_get_lval(section_vals=input, keyword_name="DFT%REAL_TIME_PROPAGATION%COM_NL")
CALL get_rtp(rtp=rtp, S_inv=S_inv)
CALL cp_fm_create(S_chol, matrix_struct=rtp%ao_ao_fmstruct, name="S_chol")
CALL cp_fm_create(S_inv_fm, matrix_struct=rtp%ao_ao_fmstruct, name="S_inv_fm")
CALL cp_fm_create(tmpS, matrix_struct=rtp%ao_ao_fmstruct)
CALL copy_dbcsr_to_fm(S_inv, S_inv_fm)
CALL cp_fm_upper_to_full(S_inv_fm, tmpS)
CALL cp_fm_get_info(S_inv_fm, nrow_global=nrow_global)
CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_chol)
CALL cp_fm_cholesky_decompose(S_chol)
CALL cp_fm_create(mat_ks, matrix_struct=S_inv_fm%matrix_struct, name="mat_ks")
CALL cp_fm_create(eigenvectors, matrix_struct=S_inv_fm%matrix_struct, name="eigenvectors")
! calculate non-local commutator
rtp_control => dft_control%rtp_control
periodic_cell = ANY(cell%perd > 0)
! relevant input parameters
com_nl = section_get_lval(section_vals=input, keyword_name="DFT%REAL_TIME_PROPAGATION%COM_NL")
len_rep = section_get_lval(section_vals=input, keyword_name="DFT%REAL_TIME_PROPAGATION%LEN_REP")
! calculate non-local commutator if necessary
IF (com_nl) THEN
CALL cite_reference(Mattiat2019)
NULLIFY (qs_kind_set, sab_orb, sap_ppnl)
CALL get_qs_env(qs_env, &
sap_ppnl=sap_ppnl, &
@ -159,23 +169,69 @@ CONTAINS
CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(idir)%matrix, sab_orb)
CALL dbcsr_set(matrix_rv(idir)%matrix, 0._dp)
END DO
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_rv=matrix_rv)
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv)
END IF
DO ispin = 1, SIZE(matrix_ks)
ALLOCATE (eigenvalues(nrow_global))
CALL cp_fm_create(mat_tmp, matrix_struct=S_inv_fm%matrix_struct, name="mat_tmp")
! calculate dipole moment matrix if required, NOT for periodic boundary conditions!
IF (len_rep) THEN
CALL cite_reference(Mattiat2022)
IF (periodic_cell) THEN
CPWARN("This application of the delta pulse is not compatible with PBC!")
END IF
! get reference point
reference = section_get_ival(section_vals=input, &
keyword_name="DFT%PRINT%MOMENTS%REFERENCE")
NULLIFY (ref_point)
CALL section_vals_val_get(input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
NULLIFY (sab_orb)
CALL get_qs_env(qs_env, sab_orb=sab_orb)
! calculate dipole moment operator
NULLIFY (matrix_r)
CALL dbcsr_allocate_matrix_set(matrix_r, 3)
DO idir = 1, 3
CALL dbcsr_init_p(matrix_r(idir)%matrix)
CALL dbcsr_create(matrix_r(idir)%matrix, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_r(idir)%matrix, sab_orb)
CALL dbcsr_set(matrix_r(idir)%matrix, 0._dp)
END DO
CALL build_local_moment_matrix(qs_env, matrix_r, 1, rcc)
END IF
! the prefactor (strength of the electric field)
kvec(:) = cell%h_inv(1, :)*rtp_control%delta_pulse_direction(1) + &
cell%h_inv(2, :)*rtp_control%delta_pulse_direction(2) + &
cell%h_inv(3, :)*rtp_control%delta_pulse_direction(3)
kvec = -kvec*twopi*rtp_control%delta_pulse_scale
! struct for fm matrices
fm_struct => rtp%ao_ao_fmstruct
! create matrices and get Cholesky decomposition of S
CALL cp_fm_create(mat_ks, matrix_struct=fm_struct, name="mat_ks")
CALL cp_fm_create(eigenvectors, matrix_struct=fm_struct, name="eigenvectors")
CALL cp_fm_create(S_chol, matrix_struct=fm_struct, name="S_chol")
CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_chol)
CALL cp_fm_cholesky_decompose(S_chol)
! get number of atomic orbitals
CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
DO ispin = 1, SIZE(matrix_ks)
! diagonalize KS matrix to get occ and virt mos
ALLOCATE (eigenvalues(nao))
CALL cp_fm_create(mat_tmp, matrix_struct=fm_struct, name="mat_tmp")
CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, mat_ks)
CALL cp_fm_cholesky_reduce(mat_ks, S_chol)
CALL cp_fm_syevd(mat_ks, mat_tmp, eigenvalues)
CALL cp_fm_cholesky_restore(mat_tmp, nrow_global, S_chol, eigenvectors, "SOLVE")
CALL cp_fm_cholesky_restore(mat_tmp, nao, S_chol, eigenvectors, "SOLVE")
! virtuals
CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo)
CALL get_mo_set(mo_set=mos(ispin)%mo_set, nmo=nmo)
nvirt = nao - nmo
CALL cp_fm_struct_create(fm_struct_tmp, para_env=S_inv_fm%matrix_struct%para_env, context=S_inv_fm%matrix_struct%context, &
nrow_global=nrow_global, ncol_global=nvirt)
CALL cp_fm_struct_create(fm_struct_tmp, para_env=fm_struct%para_env, context=fm_struct%context, &
nrow_global=nao, ncol_global=nvirt)
CALL cp_fm_create(virtuals, matrix_struct=fm_struct_tmp, name="virtuals")
CALL cp_fm_struct_release(fm_struct_tmp)
CALL cp_fm_to_fm(eigenvectors, virtuals, nvirt, nmo + 1, 1)
@ -183,7 +239,7 @@ CONTAINS
! occupied
CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1)%matrix, nmo, 1, 1)
CALL cp_fm_struct_create(fm_struct_tmp, para_env=S_inv_fm%matrix_struct%para_env, context=S_inv_fm%matrix_struct%context, &
CALL cp_fm_struct_create(fm_struct_tmp, para_env=fm_struct%para_env, context=fm_struct%context, &
nrow_global=nvirt, ncol_global=nmo)
CALL cp_fm_create(momentum, matrix_struct=fm_struct_tmp, name="momentum")
CALL cp_fm_struct_release(fm_struct_tmp)
@ -191,17 +247,17 @@ CONTAINS
! the momentum operator (in a given direction)
CALL cp_fm_set_all(mos_new(2*ispin - 1)%matrix, 0.0_dp)
! the prefactor (strength of the electric field)
kvec(:) = cell%h_inv(1, :)*dft_control%rtp_control%delta_pulse_direction(1) + &
cell%h_inv(2, :)*dft_control%rtp_control%delta_pulse_direction(2) + &
cell%h_inv(3, :)*dft_control%rtp_control%delta_pulse_direction(3)
kvec = -kvec*twopi*dft_control%rtp_control%delta_pulse_scale
DO idir = 1, 3
factor = kvec(idir)
IF (factor .NE. 0.0_dp) THEN
CALL cp_dbcsr_sm_fm_multiply(matrix_s(idir + 1)%matrix, mos_old(2*ispin - 1)%matrix, &
mos_old(2*ispin)%matrix, ncol=nmo)
IF (.NOT. len_rep) THEN
CALL cp_dbcsr_sm_fm_multiply(matrix_s(idir + 1)%matrix, mos_old(2*ispin - 1)%matrix, &
mos_old(2*ispin)%matrix, ncol=nmo)
ELSE
CALL cp_dbcsr_sm_fm_multiply(matrix_r(idir)%matrix, mos_old(2*ispin - 1)%matrix, &
mos_old(2*ispin)%matrix, ncol=nmo)
END IF
CALL cp_fm_scale_and_add(1.0_dp, mos_new(2*ispin - 1)%matrix, factor, mos_old(2*ispin)%matrix)
IF (com_nl) THEN
CALL cp_fm_set_all(mos_old(2*ispin)%matrix, 0.0_dp)
@ -215,14 +271,16 @@ CONTAINS
CALL parallel_gemm('T', 'N', nvirt, nmo, nao, 1.0_dp, virtuals, mos_new(2*ispin - 1)%matrix, 0.0_dp, momentum)
! the tricky bit ... rescale by the eigenvalue difference
CALL cp_fm_get_info(momentum, nrow_local=nrow_local, ncol_local=ncol_local, &
row_indices=row_indices, col_indices=col_indices, local_data=local_data)
DO icol = 1, ncol_local
DO irow = 1, nrow_local
factor = 1/(eigenvalues(col_indices(icol)) - eigenvalues(nmo + row_indices(irow)))
local_data(irow, icol) = factor*local_data(irow, icol)
IF (.NOT. len_rep) THEN
CALL cp_fm_get_info(momentum, nrow_local=nrow_local, ncol_local=ncol_local, &
row_indices=row_indices, col_indices=col_indices, local_data=local_data)
DO icol = 1, ncol_local
DO irow = 1, nrow_local
factor = 1/(eigenvalues(col_indices(icol)) - eigenvalues(nmo + row_indices(irow)))
local_data(irow, icol) = factor*local_data(irow, icol)
END DO
END DO
END DO
END IF
CALL cp_fm_release(mat_tmp)
DEALLOCATE (eigenvalues)
@ -232,76 +290,18 @@ CONTAINS
CALL cp_fm_release(virtuals)
CALL cp_fm_release(momentum)
! orthonormalize afterwards
CALL cp_fm_struct_create(fm_struct_tmp, para_env=S_inv_fm%matrix_struct%para_env, context=S_inv_fm%matrix_struct%context, &
nrow_global=nmo, ncol_global=nmo)
CALL cp_fm_create(oo_1, matrix_struct=fm_struct_tmp, name="oo_1")
CALL cp_fm_create(oo_2, matrix_struct=fm_struct_tmp, name="oo_2")
CALL cp_fm_struct_release(fm_struct_tmp)
CALL cp_fm_create(mat_tmp, matrix_struct=mos_old(2*ispin - 1)%matrix%matrix_struct, name="tmp_mat")
! get the complex overlap matrix
! x^T S x + y^T S y + i (-y^TS x+x^T S y)
CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mos_old(2*ispin - 1)%matrix, &
mat_tmp, ncol=nmo)
CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*ispin - 1)%matrix, mat_tmp, 0.0_dp, oo_1)
CALL parallel_gemm("T", "N", nmo, nmo, nao, -1.0_dp, mos_old(2*ispin)%matrix, mat_tmp, 0.0_dp, oo_2)
CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mos_old(2*ispin)%matrix, &
mat_tmp, ncol=nmo)
CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*ispin)%matrix, mat_tmp, 1.0_dp, oo_1)
CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*ispin - 1)%matrix, mat_tmp, 1.0_dp, oo_2)
CALL cp_fm_release(mat_tmp)
CALL cp_cfm_create(oo_c, oo_1%matrix_struct)
CALL cp_cfm_create(oo_v, oo_1%matrix_struct)
CALL cp_cfm_create(oo_vt, oo_1%matrix_struct)
oo_c%local_data = CMPLX(oo_1%local_data, oo_2%local_data, KIND=dp)
! compute inv(sqrt(overlap))
ALLOCATE (eigenvalues(nmo))
ALLOCATE (eigenvalues_sqrt(nmo))
CALL cp_cfm_heevd(oo_c, oo_v, eigenvalues)
eigenvalues_sqrt = CMPLX(1.0_dp/SQRT(eigenvalues), 0.0_dp, dp)
CALL cp_cfm_to_cfm(oo_v, oo_vt)
CALL cp_cfm_column_scale(oo_v, eigenvalues_sqrt)
DEALLOCATE (eigenvalues)
DEALLOCATE (eigenvalues_sqrt)
CALL parallel_gemm('N', 'C', nmo, nmo, nmo, (1.0_dp, 0.0_dp), &
oo_v, oo_vt, (0.0_dp, 0.0_dp), oo_c)
oo_1%local_data = REAL(oo_c%local_data, KIND=dp)
oo_2%local_data = AIMAG(oo_c%local_data)
CALL cp_cfm_release(oo_c)
CALL cp_cfm_release(oo_v)
CALL cp_cfm_release(oo_vt)
! use this to compute the orthonormal vectors
CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*ispin - 1)%matrix, oo_1, 0.0_dp, mos_new(2*ispin - 1)%matrix)
CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*ispin - 1)%matrix, oo_2, 0.0_dp, mos_new(2*ispin)%matrix)
CALL parallel_gemm("N", "N", nao, nmo, nmo, -1.0_dp, mos_old(2*ispin)%matrix, oo_2, 0.0_dp, mos_old(2*ispin - 1)%matrix)
CALL cp_fm_scale_and_add(1.0_dp, mos_old(2*ispin - 1)%matrix, 1.0_dp, mos_new(2*ispin - 1)%matrix)
CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*ispin)%matrix, oo_1, 1.0_dp, mos_new(2*ispin)%matrix)
CALL cp_fm_to_fm(mos_new(2*ispin)%matrix, mos_old(2*ispin)%matrix)
CALL cp_fm_release(oo_1)
CALL cp_fm_release(oo_2)
END DO
! release matrices
CALL cp_fm_release(S_chol)
CALL cp_fm_release(mat_ks)
CALL cp_fm_release(eigenvectors)
IF (com_nl) CALL dbcsr_deallocate_matrix_set(matrix_rv)
IF (len_rep) CALL dbcsr_deallocate_matrix_set(matrix_r)
! orthonormalize afterwards
CALL orthonormalize_complex_mos(qs_env, mos_old)
!***************************************************************
!remove later
CALL cp_fm_release(S_inv_fm)
CALL cp_fm_release(tmpS)
!**************************************************************
CALL timestop(handle)
END SUBROUTINE apply_delta_pulse_periodic
@ -320,102 +320,321 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse'
COMPLEX(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_sqrt
INTEGER :: handle, i, nao, nmo
REAL(KIND=dp), DIMENSION(3) :: kvec
REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
TYPE(cell_type), POINTER :: cell
TYPE(cp_cfm_type) :: oo_c, oo_v, oo_vt
TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
TYPE(cp_fm_type) :: mat_S, oo_1, oo_2, S_inv_fm, tmp
TYPE(cp_fm_type) :: S_inv_fm, tmp
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
TYPE(dbcsr_type), POINTER :: cosmat, S_inv, sinmat
TYPE(dbcsr_type), POINTER :: cosmat, sinmat
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos
TYPE(rt_prop_type), POINTER :: rtp
NULLIFY (dft_control)
TYPE(rtp_control_type), POINTER :: rtp_control
CALL timeset(routineN, handle)
! we need the inverse overlap
NULLIFY (cell, dft_control, matrix_s, mos, rtp, rtp_control)
CALL get_qs_env(qs_env, &
mos=mos, &
rtp=rtp, &
cell=cell, &
dft_control=dft_control, &
matrix_s=matrix_s, &
dft_control=dft_control)
CALL get_rtp(rtp=rtp, S_inv=S_inv)
CALL cp_fm_create(S_inv_fm, matrix_struct=rtp%ao_ao_fmstruct, name="tmp_mat")
CALL cp_fm_create(tmp, matrix_struct=rtp%ao_ao_fmstruct, name="tmp_mat")
CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_inv_fm)
CALL cp_fm_cholesky_decompose(S_inv_fm)
CALL cp_fm_cholesky_invert(S_inv_fm)
CALL cp_fm_upper_to_full(S_inv_fm, tmp)
CALL cp_fm_create(mat_S, matrix_struct=S_inv_fm%matrix_struct, name="mat_S")
CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, mat_S)
CALL cp_fm_upper_to_full(mat_S, tmp)
CALL cp_fm_release(tmp)
! we need the berry matrix
CALL get_qs_env(qs_env, cell=cell)
mos=mos, &
rtp=rtp)
rtp_control => dft_control%rtp_control
! direction ... unscaled, this will yield a exp(ikr) that is periodic with the cell
kvec(:) = cell%h_inv(1, :)*dft_control%rtp_control%delta_pulse_direction(1) + &
cell%h_inv(2, :)*dft_control%rtp_control%delta_pulse_direction(2) + &
cell%h_inv(3, :)*dft_control%rtp_control%delta_pulse_direction(3)
kvec(:) = cell%h_inv(1, :)*rtp_control%delta_pulse_direction(1) + &
cell%h_inv(2, :)*rtp_control%delta_pulse_direction(2) + &
cell%h_inv(3, :)*rtp_control%delta_pulse_direction(3)
kvec = -kvec*twopi
! scaling will make the things not periodic with the cell, which would only be good for gas phase systems ?
kvec(:) = dft_control%rtp_control%delta_pulse_scale*kvec
! calculate exponentials (= Berry moments)
NULLIFY (cosmat, sinmat)
ALLOCATE (cosmat, sinmat)
CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
! apply inv(S)*operator to C
! need inverse of overlap matrix
CALL cp_fm_create(S_inv_fm, matrix_struct=rtp%ao_ao_fmstruct, name="S_inv_fm")
CALL cp_fm_create(tmp, matrix_struct=rtp%ao_ao_fmstruct, name="tmp_mat")
CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_inv_fm)
CALL cp_fm_cholesky_decompose(S_inv_fm)
CALL cp_fm_cholesky_invert(S_inv_fm)
CALL cp_fm_upper_to_full(S_inv_fm, tmp)
CALL cp_fm_release(tmp)
DO i = 1, SIZE(mos)
! apply exponentials to mo coefficients
CALL get_mo_set(mos(i)%mo_set, nao=nao, nmo=nmo)
CALL cp_dbcsr_sm_fm_multiply(cosmat, mos(i)%mo_set%mo_coeff, mos_new(2*i - 1)%matrix, ncol=nmo)
CALL cp_dbcsr_sm_fm_multiply(sinmat, mos(i)%mo_set%mo_coeff, mos_new(2*i)%matrix, ncol=nmo)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_inv_fm, mos_new(2*i - 1)%matrix, 0.0_dp, mos_old(2*i - 1)%matrix)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_inv_fm, mos_new(2*i)%matrix, 0.0_dp, mos_old(2*i)%matrix)
END DO
! in a finite basis, unfortunately, inv(S)*operator is not unitary, so orthonormalize afterwards
CALL cp_fm_struct_create(fm_struct_tmp, para_env=S_inv_fm%matrix_struct%para_env, context=S_inv_fm%matrix_struct%context, &
CALL cp_fm_release(S_inv_fm)
CALL dbcsr_deallocate_matrix(cosmat)
CALL dbcsr_deallocate_matrix(sinmat)
! orthonormalize afterwards
CALL orthonormalize_complex_mos(qs_env, mos_old)
CALL timestop(handle)
END SUBROUTINE apply_delta_pulse
! **************************************************************************************************
!> \brief apply magnetic delta pulse to linear order
!> \param qs_env ...
!> \param mos_old ...
!> \param mos_new ...
! **************************************************************************************************
SUBROUTINE apply_delta_pulse_mag(qs_env, mos_old, mos_new)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_old, mos_new
CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_mag'
INTEGER :: gauge_orig, handle, idir, ispin, nao, &
nmo, nrow_global, nvirt
REAL(KIND=dp) :: eps_ppnl, factor
REAL(KIND=dp), DIMENSION(3) :: kvec, rcc
REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, ref_point
TYPE(cell_type), POINTER :: cell
TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
TYPE(cp_fm_type) :: eigenvectors, mat_ks, perturbation, &
S_chol, virtuals
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_mag, matrix_nl, &
matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_all, sab_orb, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(rt_prop_type), POINTER :: rtp
TYPE(section_vals_type), POINTER :: input
CALL timeset(routineN, handle)
CALL cite_reference(Mattiat2022)
NULLIFY (rtp, dft_control, matrix_ks, matrix_s, input, mos, cell, sab_orb, sab_all, sap_ppnl, &
qs_kind_set, particle_set)
CALL get_qs_env(qs_env, &
rtp=rtp, &
dft_control=dft_control, &
mos=mos, &
matrix_ks=matrix_ks, &
matrix_s=matrix_s, &
input=input, &
cell=cell, &
sab_orb=sab_orb, &
sab_all=sab_all, &
sap_ppnl=sap_ppnl)
gauge_orig = section_get_ival(section_vals=input, &
keyword_name="DFT%REAL_TIME_PROPAGATION%GAUGE_ORIG")
NULLIFY (ref_point)
CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%GAUGE_ORIG_MANUAL", r_vals=ref_point)
CALL get_reference_point(rcc, qs_env=qs_env, reference=gauge_orig, ref_point=ref_point)
! Create fm matrices
CALL cp_fm_create(S_chol, matrix_struct=rtp%ao_ao_fmstruct, name='Cholesky S')
CALL cp_fm_create(eigenvectors, matrix_struct=rtp%ao_ao_fmstruct, name="gs evecs fm")
CALL cp_fm_create(mat_ks, matrix_struct=rtp%ao_ao_fmstruct, name='KS matrix')
! get nrows_global
CALL cp_fm_get_info(mat_ks, nrow_global=nrow_global)
! cholesky decomposition of overlap matrix
CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_chol)
CALL cp_fm_cholesky_decompose(S_chol)
! initiate perturbation matrix
NULLIFY (matrix_mag)
CALL dbcsr_allocate_matrix_set(matrix_mag, 3)
DO idir = 1, 3
CALL dbcsr_init_p(matrix_mag(idir)%matrix)
CALL dbcsr_create(matrix_mag(idir)%matrix, template=matrix_s(1)%matrix, &
matrix_type=dbcsr_type_antisymmetric)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_mag(idir)%matrix, sab_orb)
CALL dbcsr_set(matrix_mag(idir)%matrix, 0._dp)
END DO
! construct magnetic dipole moment matrix
CALL build_local_magmom_matrix(qs_env, matrix_mag, 1, ref_point=rcc)
! work matrix for non-local potential part if necessary
NULLIFY (matrix_nl)
IF (ASSOCIATED(sap_ppnl)) THEN
CALL dbcsr_allocate_matrix_set(matrix_nl, 3)
DO idir = 1, 3
CALL dbcsr_init_p(matrix_nl(idir)%matrix)
CALL dbcsr_create(matrix_nl(idir)%matrix, template=matrix_s(1)%matrix, &
matrix_type=dbcsr_type_antisymmetric)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_nl(idir)%matrix, sab_orb)
CALL dbcsr_set(matrix_nl(idir)%matrix, 0._dp)
END DO
! construct non-local contribution
CALL get_qs_env(qs_env, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set)
eps_ppnl = dft_control%qs_control%eps_ppnl
CALL build_com_nl_mag(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_nl, rcc, cell)
DO idir = 1, 3
CALL dbcsr_add(matrix_mag(idir)%matrix, matrix_nl(idir)%matrix, -one, one)
END DO
CALL dbcsr_deallocate_matrix_set(matrix_nl)
END IF
DO ispin = 1, dft_control%nspins
! allocate eigenvalues
NULLIFY (eigenvalues)
ALLOCATE (eigenvalues(nrow_global))
! diagonalize KS matrix in AO basis using Cholesky decomp. of S
CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, mat_ks)
CALL cp_fm_cholesky_reduce(mat_ks, S_chol)
CALL cp_fm_syevd(mat_ks, eigenvectors, eigenvalues)
CALL cp_fm_triangular_multiply(S_chol, eigenvectors, invert_tr=.TRUE.)
! virtuals
CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo)
nvirt = nao - nmo
CALL cp_fm_struct_create(fm_struct_tmp, para_env=rtp%ao_ao_fmstruct%para_env, context=rtp%ao_ao_fmstruct%context, &
nrow_global=nrow_global, ncol_global=nvirt)
CALL cp_fm_create(virtuals, matrix_struct=fm_struct_tmp, name="virtuals")
CALL cp_fm_struct_release(fm_struct_tmp)
CALL cp_fm_to_fm(eigenvectors, virtuals, nvirt, nmo + 1, 1)
! occupied
CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1)%matrix, nmo, 1, 1)
CALL cp_fm_struct_create(fm_struct_tmp, para_env=rtp%ao_ao_fmstruct%para_env, context=rtp%ao_ao_fmstruct%context, &
nrow_global=nvirt, ncol_global=nmo)
CALL cp_fm_create(perturbation, matrix_struct=fm_struct_tmp, name="perturbation")
CALL cp_fm_struct_release(fm_struct_tmp)
! apply perturbation
CALL cp_fm_set_all(mos_new(2*ispin - 1)%matrix, 0.0_dp)
! the prefactor (strength of the magnetic field, should be divided by 2c)
kvec(:) = cell%h_inv(1, :)*dft_control%rtp_control%delta_pulse_direction(1) + &
cell%h_inv(2, :)*dft_control%rtp_control%delta_pulse_direction(2) + &
cell%h_inv(3, :)*dft_control%rtp_control%delta_pulse_direction(3)
kvec = -kvec*twopi*dft_control%rtp_control%delta_pulse_scale
DO idir = 1, 3
factor = kvec(idir)/2 ! divide by two b.c. of pre-factor for magnetic term
IF (factor .NE. 0.0_dp) THEN
CALL cp_dbcsr_sm_fm_multiply(matrix_mag(idir)%matrix, mos_old(2*ispin - 1)%matrix, &
mos_old(2*ispin)%matrix, ncol=nmo)
CALL cp_fm_scale_and_add(1.0_dp, mos_new(2*ispin - 1)%matrix, factor, mos_old(2*ispin)%matrix)
END IF
END DO
CALL parallel_gemm('T', 'N', nvirt, nmo, nao, 1.0_dp, virtuals, mos_new(2*ispin - 1)%matrix, 0.0_dp, perturbation)
DEALLOCATE (eigenvalues)
! now obtain the initial condition in mos_old
CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1)%matrix, nmo, 1, 1)
CALL parallel_gemm("N", "N", nao, nmo, nvirt, 1.0_dp, virtuals, perturbation, 0.0_dp, mos_old(2*ispin)%matrix)
CALL cp_fm_release(virtuals)
CALL cp_fm_release(perturbation)
END DO
! deallocations
CALL cp_fm_release(S_chol)
CALL cp_fm_release(mat_ks)
CALL cp_fm_release(eigenvectors)
CALL dbcsr_deallocate_matrix_set(matrix_mag)
! orthonormalize afterwards
CALL orthonormalize_complex_mos(qs_env, mos_old)
CALL timestop(handle)
END SUBROUTINE apply_delta_pulse_mag
! **************************************************************************************************
!> \brief orthonormalize complex mos, e. g. after non-unitary transformations using Löwdin's algorithm
!> \param qs_env ...
!> \param coeffs ...
! **************************************************************************************************
SUBROUTINE orthonormalize_complex_mos(qs_env, coeffs)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cp_fm_p_type), DIMENSION(:), INTENT(INOUT), &
POINTER :: coeffs
COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_sqrt
INTEGER :: im, ispin, j, nao, nmo, nspins, re
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_cfm_type) :: oo_c, oo_v, oo_vt
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: coeffs_tmp
TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
TYPE(cp_fm_type) :: oo_1, oo_2, S_fm, tmp
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos
NULLIFY (para_env, blacs_env, dft_control, matrix_s, mos)
CALL get_qs_env(qs_env, &
blacs_env=blacs_env, &
dft_control=dft_control, &
matrix_s=matrix_s, &
mos=mos, &
para_env=para_env)
nspins = dft_control%nspins
CALL cp_fm_get_info(coeffs(1)%matrix, nrow_global=nao)
! get overlap matrix
CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nao, &
context=blacs_env, para_env=para_env)
CALL cp_fm_create(S_fm, matrix_struct=fm_struct_tmp, name="overlap fm")
CALL cp_fm_struct_release(fm_struct_tmp)
! copy overlap matrix
CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_fm)
DO ispin = 1, nspins
CALL get_mo_set(mos(ispin)%mo_set, nmo=nmo)
CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
nrow_global=nmo, ncol_global=nmo)
CALL cp_fm_create(oo_1, matrix_struct=fm_struct_tmp, name="oo_1")
CALL cp_fm_create(oo_2, matrix_struct=fm_struct_tmp, name="oo_2")
CALL cp_fm_struct_release(fm_struct_tmp)
CALL cp_fm_create(tmp, matrix_struct=mos_old(2*i - 1)%matrix%matrix_struct, name="tmp_mat")
! get the complex overlap matrix
CALL cp_fm_create(tmp, matrix_struct=coeffs(2*ispin - 1)%matrix%matrix_struct, name="tmp_mat")
! get the complex overlap matrix in MO basis
! x^T S x + y^T S y + i (-y^TS x+x^T S y)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, mat_S, mos_old(2*i - 1)%matrix, 0.0_dp, tmp)
CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*i - 1)%matrix, tmp, 0.0_dp, oo_1)
CALL parallel_gemm("T", "N", nmo, nmo, nao, -1.0_dp, mos_old(2*i)%matrix, tmp, 0.0_dp, oo_2)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_fm, coeffs(2*ispin - 1)%matrix, 0.0_dp, tmp)
CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, coeffs(2*ispin - 1)%matrix, tmp, 0.0_dp, oo_1)
CALL parallel_gemm("T", "N", nmo, nmo, nao, -1.0_dp, coeffs(2*ispin)%matrix, tmp, 0.0_dp, oo_2)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, mat_S, mos_old(2*i)%matrix, 0.0_dp, tmp)
CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*i)%matrix, tmp, 1.0_dp, oo_1)
CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mos_old(2*i - 1)%matrix, tmp, 1.0_dp, oo_2)
CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_fm, coeffs(2*ispin)%matrix, 0.0_dp, tmp)
CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, coeffs(2*ispin)%matrix, tmp, 1.0_dp, oo_1)
CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, coeffs(2*ispin - 1)%matrix, tmp, 1.0_dp, oo_2)
CALL cp_fm_release(tmp)
! complex Löwdin
CALL cp_cfm_create(oo_c, oo_1%matrix_struct)
CALL cp_cfm_create(oo_v, oo_1%matrix_struct)
CALL cp_cfm_create(oo_vt, oo_1%matrix_struct)
oo_c%local_data = CMPLX(oo_1%local_data, oo_2%local_data, KIND=dp)
! compute inv(sqrt(overlap))
ALLOCATE (eigenvalues(nmo))
ALLOCATE (eigenvalues_sqrt(nmo))
CALL cp_cfm_heevd(oo_c, oo_v, eigenvalues)
eigenvalues_sqrt = CMPLX(1.0_dp/SQRT(eigenvalues), 0.0_dp, dp)
eigenvalues_sqrt(:) = CMPLX(one/SQRT(eigenvalues(:)), zero, dp)
CALL cp_cfm_to_cfm(oo_v, oo_vt)
CALL cp_cfm_column_scale(oo_v, eigenvalues_sqrt)
DEALLOCATE (eigenvalues)
@ -428,34 +647,33 @@ CONTAINS
CALL cp_cfm_release(oo_v)
CALL cp_cfm_release(oo_vt)
! use this to compute the orthonormal vectors
CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*i - 1)%matrix, oo_1, 0.0_dp, mos_new(2*i - 1)%matrix)
CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*i - 1)%matrix, oo_2, 0.0_dp, mos_new(2*i)%matrix)
! transform coefficients accordingly
NULLIFY (coeffs_tmp)
ALLOCATE (coeffs_tmp(2))
DO j = 1, 2
NULLIFY (coeffs_tmp(j)%matrix)
ALLOCATE (coeffs_tmp(j)%matrix)
CALL cp_fm_create(coeffs_tmp(j)%matrix, matrix_struct=coeffs(2*(ispin - 1) + j)%matrix%matrix_struct)
END DO
CALL parallel_gemm("N", "N", nao, nmo, nmo, -1.0_dp, mos_old(2*i)%matrix, oo_2, 0.0_dp, mos_old(2*i - 1)%matrix)
CALL cp_fm_scale_and_add(1.0_dp, mos_old(2*i - 1)%matrix, 1.0_dp, mos_new(2*i - 1)%matrix)
! indices for coeffs_tmp
re = 1
im = 2
CALL parallel_gemm("N", "N", nao, nmo, nmo, one, coeffs(2*ispin - 1)%matrix, oo_1, zero, coeffs_tmp(re)%matrix)
CALL parallel_gemm("N", "N", nao, nmo, nmo, one, coeffs(2*ispin - 1)%matrix, oo_2, zero, coeffs_tmp(im)%matrix)
CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mos_old(2*i)%matrix, oo_1, 1.0_dp, mos_new(2*i)%matrix)
CALL cp_fm_to_fm(mos_new(2*i)%matrix, mos_old(2*i)%matrix)
CALL parallel_gemm("N", "N", nao, nmo, nmo, -one, coeffs(2*ispin)%matrix, oo_2, zero, coeffs(2*ispin - 1)%matrix)
CALL cp_fm_scale_and_add(one, coeffs(2*ispin - 1)%matrix, one, coeffs_tmp(re)%matrix)
CALL parallel_gemm("N", "N", nao, nmo, nmo, one, coeffs(2*ispin)%matrix, oo_1, one, coeffs_tmp(im)%matrix)
CALL cp_fm_to_fm(coeffs_tmp(im)%matrix, coeffs(2*ispin)%matrix)
CALL cp_fm_vect_dealloc(coeffs_tmp)
CALL cp_fm_release(oo_1)
CALL cp_fm_release(oo_2)
END DO
CALL cp_fm_release(S_fm)
CALL cp_fm_release(mat_S)
CALL dbcsr_deallocate_matrix(cosmat)
CALL dbcsr_deallocate_matrix(sinmat)
!***************************************************************
!remove later
CALL copy_fm_to_dbcsr(S_inv_fm, S_inv)
CALL dbcsr_filter(S_inv, rtp%filter_eps)
CALL cp_fm_release(S_inv_fm)
!**************************************************************
CALL timestop(handle)
END SUBROUTINE apply_delta_pulse
END SUBROUTINE orthonormalize_complex_mos
END MODULE rt_delta_pulse

141
src/emd/rt_hfx_utils.F
View file

@ -22,8 +22,10 @@ MODULE rt_hfx_utils
USE dbcsr_api, ONLY: dbcsr_create,&
dbcsr_p_type,&
dbcsr_set,&
dbcsr_type,&
dbcsr_type_antisymmetric
USE kinds, ONLY: dp
USE kinds, ONLY: default_string_length,&
dp
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_ks_types, ONLY: qs_ks_env_type,&
@ -32,7 +34,6 @@ MODULE rt_hfx_utils
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_set,&
qs_rho_type
USE rt_propagation_types, ONLY: rt_prop_type
#include "../base/base_uses.f90"
IMPLICIT NONE
@ -52,33 +53,27 @@ CONTAINS
SUBROUTINE rtp_hfx_rebuild(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_im, matrix_ks_im, &
matrix_s, matrix_s_aux, rho_ao_im, &
LOGICAL :: need_h_im = .FALSE.
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_im, matrix_s_aux, &
rho_aux_ao_im
TYPE(dft_control_type), POINTER :: dft_control
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_aux, sab_orb
POINTER :: sab_aux
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho, rho_aux
TYPE(rt_prop_type), POINTER :: rtp
NULLIFY (ks_env, matrix_s, matrix_ks_im, dft_control, sab_orb, rtp, rho, rho_ao_im)
NULLIFY (ks_env, dft_control, rho)
NULLIFY (sab_aux, rho_aux, rho_aux_ao_im, matrix_ks_aux_im, matrix_s_aux)
CALL get_qs_env(qs_env, &
ks_env=ks_env, &
dft_control=dft_control, &
rtp=rtp, &
matrix_s=matrix_s, &
rho=rho, &
matrix_ks_im=matrix_ks_im, &
sab_orb=sab_orb)
rho=rho)
CALL qs_rho_get(rho, rho_ao_im=rho_ao_im)
CALL rebuild_matrices(rho_ao_im, matrix_ks_im, sab_orb, matrix_s, &
dft_control%nspins)
CALL set_ks_env(ks_env, matrix_ks_im=matrix_ks_im)
CALL qs_rho_set(rho, rho_ao_im=rho_ao_im)
need_h_im = dft_control%rtp_control%velocity_gauge
CALL rebuild_matrices_im_kp(qs_env, rho, ks_env, dft_control%nspins, dft_control%nimages, &
rebuild_matrix_p=.TRUE., rebuild_matrix_ks=.TRUE., rebuild_matrix_h=need_h_im)
IF (dft_control%do_admm) THEN
CALL get_admm_env(qs_env%admm_env, &
@ -142,4 +137,116 @@ CONTAINS
END SUBROUTINE rebuild_matrices
! **************************************************************************************************
!> \brief (Re)-allocate imaginary parts of density matrix and KS-matrix as kpoint_transitional_types
!> \param qs_env ...
!> \param rho ...
!> \param ks_env ...
!> \param nspins ...
!> \param nimages ...
!> \param rebuild_matrix_p ...
!> \param rebuild_matrix_ks ...
!> \param rebuild_matrix_h ...
! **************************************************************************************************
SUBROUTINE rebuild_matrices_im_kp(qs_env, rho, ks_env, nspins, nimages, rebuild_matrix_p, rebuild_matrix_ks, rebuild_matrix_h)
TYPE(qs_environment_type), INTENT(inout), POINTER :: qs_env
TYPE(qs_rho_type), INTENT(INOUT), POINTER :: rho
TYPE(qs_ks_env_type), INTENT(INOUT), POINTER :: ks_env
INTEGER, INTENT(in), OPTIONAL :: nspins, nimages
LOGICAL, INTENT(in), OPTIONAL :: rebuild_matrix_p, rebuild_matrix_ks, &
rebuild_matrix_h
CHARACTER(LEN=default_string_length) :: headline
INTEGER :: image, ispin
LOGICAL :: my_rebuild_h = .FALSE., &
my_rebuild_ks = .FALSE., &
my_rebuild_p = .FALSE.
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h_im, matrix_h_kp, matrix_ks_im, &
matrix_ks_kp, matrix_p_im, rho_ao_kp
TYPE(dbcsr_type), POINTER :: refmatrix
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb
IF (PRESENT(rebuild_matrix_p)) my_rebuild_p = rebuild_matrix_p
IF (PRESENT(rebuild_matrix_ks)) my_rebuild_ks = rebuild_matrix_ks
IF (PRESENT(rebuild_matrix_h)) my_rebuild_h = rebuild_matrix_h
NULLIFY (matrix_ks_kp, rho_ao_kp, sab_orb, matrix_h_kp)
CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks_kp, matrix_h_kp=matrix_h_kp, sab_orb=sab_orb)
CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
IF (my_rebuild_p) THEN
NULLIFY (matrix_p_im)
CALL qs_rho_get(rho, rho_ao_im_kp=matrix_p_im)
CALL dbcsr_allocate_matrix_set(matrix_p_im, nspins, nimages)
refmatrix => rho_ao_kp(1, 1)%matrix
DO ispin = 1, nspins
DO image = 1, nimages
IF (nspins > 1) THEN
IF (ispin == 1) THEN
headline = "IMAGINARY DENSITY MATRIX FOR ALPHA SPIN"
ELSE
headline = "IMAGINARY DENSITY MATRIX FOR BETA SPIN"
END IF
ELSE
headline = "IMAGINARY DENSITY MATRIX"
END IF
ALLOCATE (matrix_p_im(ispin, image)%matrix)
CALL dbcsr_create(matrix=matrix_p_im(ispin, image)%matrix, template=refmatrix, &
name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric, nze=0)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_p_im(ispin, image)%matrix, sab_orb)
CALL dbcsr_set(matrix_p_im(ispin, image)%matrix, 0.0_dp)
END DO
END DO
CALL qs_rho_set(rho, rho_ao_im_kp=matrix_p_im)
END IF
IF (my_rebuild_ks) THEN
NULLIFY (matrix_ks_im)
CALL get_qs_env(qs_env, matrix_ks_im_kp=matrix_ks_im)
CALL dbcsr_allocate_matrix_set(matrix_ks_im, nspins, nimages)
refmatrix => matrix_ks_kp(1, 1)%matrix
DO ispin = 1, nspins
DO image = 1, nimages
IF (nspins > 1) THEN
IF (ispin == 1) THEN
headline = "IMAGINARY KOHN-SHAM MATRIX FOR ALPHA SPIN"
ELSE
headline = "IMAGINARY KOHN-SHAM MATRIX FOR BETA SPIN"
END IF
ELSE
headline = "IMAGINARY KOHN-SHAM MATRIX"
END IF
ALLOCATE (matrix_ks_im(ispin, image)%matrix)
CALL dbcsr_create(matrix=matrix_ks_im(ispin, image)%matrix, template=refmatrix, &
name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric, nze=0)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_im(ispin, image)%matrix, sab_orb)
CALL dbcsr_set(matrix_ks_im(ispin, image)%matrix, 0.0_dp)
END DO
END DO
CALL set_ks_env(ks_env, matrix_ks_im_kp=matrix_ks_im)
END IF
IF (my_rebuild_h) THEN
NULLIFY (matrix_h_im)
CALL get_qs_env(qs_env, matrix_h_im_kp=matrix_h_im)
CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
refmatrix => matrix_h_kp(1, 1)%matrix
DO image = 1, nimages
headline = "IMAGINARY CORE HAMILTONIAN MATRIX"
ALLOCATE (matrix_h_im(1, image)%matrix)
CALL dbcsr_create(matrix=matrix_h_im(1, image)%matrix, template=refmatrix, &
name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric, nze=0)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, image)%matrix, sab_orb)
CALL dbcsr_set(matrix_h_im(1, image)%matrix, 0.0_dp)
END DO
CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
END IF
END SUBROUTINE
END MODULE rt_hfx_utils

11
src/emd/rt_propagation_methods.F
View file

@ -45,7 +45,7 @@ MODULE rt_propagation_methods
dbcsr_filter, dbcsr_frobenius_norm, dbcsr_get_block_p, dbcsr_init_p, &
dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type
dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_antisymmetric
USE input_constants, ONLY: do_arnoldi,&
do_bch,&
do_em,&
@ -73,6 +73,7 @@ MODULE rt_propagation_methods
USE rt_propagation_utils, ONLY: calc_S_derivs,&
calc_update_rho,&
calc_update_rho_sparse
USE rt_propagation_velocity_gauge, ONLY: velocity_gauge_ks_matrix
#include "../base/base_uses.f90"
IMPLICIT NONE
@ -127,6 +128,12 @@ CONTAINS
CALL get_qs_env(qs_env, matrix_s=matrix_s)
IF (rtp%iter == 1) THEN
CALL qs_energies_init(qs_env, .FALSE.)
! add additional terms for the velocity gauge to matrix_h and matrix_h_im
! should be called imediately after qs_energies_init and before qs_ks_update_qs_env
IF (rtp_control%velocity_gauge) THEN
CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.FALSE.)
END IF
CALL get_qs_env(qs_env, matrix_s=matrix_s)
IF (.NOT. rtp_control%fixed_ions) THEN
CALL s_matrices_create(matrix_s, rtp)
@ -186,7 +193,7 @@ CONTAINS
CALL dbcsr_init_p(ks_mix(i)%matrix)
CALL dbcsr_create(ks_mix(i)%matrix, template=matrix_ks(1)%matrix)
CALL dbcsr_init_p(ks_mix_im(i)%matrix)
CALL dbcsr_create(ks_mix_im(i)%matrix, template=matrix_ks(1)%matrix)
CALL dbcsr_create(ks_mix_im(i)%matrix, template=matrix_ks(1)%matrix, matrix_type=dbcsr_type_antisymmetric)
END DO
DO i = 1, SIZE(matrix_ks)
re = 2*i - 1

15
src/hfx_admm_utils.F
View file

@ -738,7 +738,8 @@ CONTAINS
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_1d, matrix_ks_aux_fit, &
matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_im, matrix_ks_im, rho_ao_1d, rho_ao_resp
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks_orb, rho_ao_orb
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im, matrix_ks_orb, &
rho_ao_orb
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array
@ -755,13 +756,14 @@ CONTAINS
NULLIFY (auxbas_pw_pool, dft_control, hfx_sections, input, &
para_env, poisson_env, pw_env, virial, matrix_ks_im, &
matrix_ks_orb, rho_ao_orb, matrix_h, matrix_ks_aux_fit, &
matrix_ks_orb, rho_ao_orb, matrix_h, matrix_h_im, matrix_ks_aux_fit, &
matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx)
CALL get_qs_env(qs_env=qs_env, &
dft_control=dft_control, &
input=input, &
matrix_h_kp=matrix_h, &
matrix_h_im_kp=matrix_h_im, &
para_env=para_env, &
pw_env=pw_env, &
virial=virial, &
@ -953,6 +955,15 @@ CONTAINS
!! If required, the calculation of the forces will be done later with adiabatic rescaling
IF (do_adiabatic_rescaling) hf_energy(irep) = ehfx + ehfxrt
IF (dft_control%rtp_control%velocity_gauge) THEN
CPASSERT(ASSOCIATED(matrix_h_im))
DO ispin = 1, nspins
CALL dbcsr_add(matrix_ks_im(ispin)%matrix, matrix_h_im(1, 1)%matrix, &
1.0_dp, 1.0_dp)
END DO
END IF
END IF
IF (.NOT. qs_env%run_rtp) THEN

61
src/input_cp2k_dft.F
View file

@ -8277,9 +8277,68 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="APPLY_DELTA_PULSE_MAG", &
description="Applies a magnetic delta kick to the initial wfn (only RTP for now - the EMD "// &
" case is not yet implemented).", &
usage="APPLY_DELTA_PULSE_MAG", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="VELOCITY_GAUGE", &
description="Perform propagation in the velocity gauge; only impulse at the moment (delta-pulse) "// &
"only works for all-electron up to now."// &
"uses DELTA_PULSE_SCALE and DELTA_PULSE_DIRECTION", &
usage="VELOCITY_GAUGE T", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="GAUGE_ORIG", &
description="Define gauge origin for magnetic perturbation", &
usage="GAUGE_ORIG COM", &
enum_c_vals=s2a("COM", "COAC", "USER_DEFINED", "ZERO"), &
enum_desc=s2a("Use Center of Mass", &
"Use Center of Atomic Charges", &
"Use User Defined Point (Keyword:REF_POINT)", &
"Use Origin of Coordinate System"), &
enum_i_vals=(/use_mom_ref_com, &
use_mom_ref_coac, &
use_mom_ref_user, &
use_mom_ref_zero/), &
default_i_val=use_mom_ref_com)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="GAUGE_ORIG_MANUAL", &
description="Manually defined gauge origin for magnetic perturbation [in Bohr!]", &
usage="SPATIAL_REF_POINT x y z", &
repeats=.FALSE., &
n_var=3, default_r_vals=(/0._dp, 0._dp, 0._dp/), &
type_of_var=real_t, &
unit_str='bohr')
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="VG_COM_NL", &
description="apply gauge transformed non-local potential term", &
usage="VG_COM_NL T", &
default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="COM_NL", &
description="Include non local commutator for periodic delta pulse.", &
description="Include non-local commutator for periodic delta pulse.", &
usage="COM_NL", &
default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="LEN_REP", &
description="Use length representation delta pulse (in conjunction with PERIODIC T). "// &
"Note that this is NOT compatible with PBC!!!"// &
" Uses the reference point defined in DFT%PRINT%MOMENTS ", &
usage="LEN_REP T", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

48
src/motion/rt_propagation.F
View file

@ -13,6 +13,7 @@
MODULE rt_propagation
USE bibliography, ONLY: Andermatt2016,&
cite_reference
USE cell_types, ONLY: cell_type
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type,&
rtp_control_type
@ -29,6 +30,7 @@ MODULE rt_propagation
cp_to_string
USE cp_output_handling, ONLY: cp_add_iter_level,&
cp_iterate,&
cp_print_key_unit_nr,&
cp_rm_iter_level
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: dbcsr_copy,&
@ -67,6 +69,7 @@ MODULE rt_propagation
init_mo_set,&
mo_set_p_type
USE rt_delta_pulse, ONLY: apply_delta_pulse,&
apply_delta_pulse_mag,&
apply_delta_pulse_periodic
USE rt_hfx_utils, ONLY: rtp_hfx_rebuild
USE rt_propagation_methods, ONLY: propagation_step
@ -79,6 +82,7 @@ MODULE rt_propagation
calc_update_rho,&
calc_update_rho_sparse,&
get_restart_wfn
USE rt_propagation_velocity_gauge, ONLY: velocity_gauge_ks_matrix
USE rt_propagator_init, ONLY: init_propagators,&
rt_initialize_rho_from_mos
#include "../base/base_uses.f90"
@ -153,7 +157,10 @@ 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
rtp%track_imag_density = magnetic .OR. vel_reprs .OR. rtp_control%velocity_gauge
IF (rtp_control%velocity_gauge) THEN
rtp%do_hfx = .TRUE.
END IF
! Hmm, not really like to initialize with the structure of S but I reckon it is
! done everywhere like this
IF (rtp%do_hfx .OR. rtp%track_imag_density) &
@ -191,6 +198,8 @@ CONTAINS
INTEGER :: i, ispin, re
INTEGER, DIMENSION(2) :: nelectron_spin
LOGICAL :: periodic_cell
TYPE(cell_type), POINTER :: cell
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mos_new, mos_old
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_new, rho_old
TYPE(dft_control_type), POINTER :: dft_control
@ -198,16 +207,19 @@ CONTAINS
TYPE(rt_prop_type), POINTER :: rtp
TYPE(rtp_control_type), POINTER :: rtp_control
NULLIFY (matrix_s, dft_control)
NULLIFY (matrix_s, dft_control, cell, rtp, rtp_control)
CALL cite_reference(Andermatt2016)
CALL get_qs_env(qs_env, &
cell=cell, &
rtp=rtp, &
matrix_s=matrix_s, &
dft_control=dft_control)
rtp_control => dft_control%rtp_control
periodic_cell = ANY(cell%perd > 0)
IF (.NOT. rtp%linear_scaling) THEN
CALL get_qs_env(qs_env, mos=mos)
CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
@ -218,8 +230,16 @@ CONTAINS
IF (rtp_control%periodic) THEN
CALL apply_delta_pulse_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(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
DO i = 1, SIZE(mos)
CALL cp_fm_to_fm(mos(i)%mo_set%mo_coeff, mos_old(2*i - 1)%matrix)
@ -254,6 +274,12 @@ CONTAINS
CALL calc_update_rho_sparse(qs_env)
END IF
! Modify KS matrix to include the additional terms in the velocity gauge
IF (rtp_control%velocity_gauge) THEN
! As matrix_h and matrix_h_im are not updated by qs_ks_update_qs_env()
! the non-gauge transformed non-local part has to be subtracted here
CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.TRUE.)
END IF
CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE.)
END SUBROUTINE init_propagation_run
@ -275,7 +301,7 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'run_propagation'
INTEGER :: aspc_order, handle, i_iter, i_step, &
max_iter, max_steps
max_iter, max_steps, output_unit
LOGICAL :: should_stop
REAL(Kind=dp) :: eps_ener, time_iter_start, &
time_iter_stop, used_time
@ -285,19 +311,20 @@ CONTAINS
TYPE(qs_energy_type), POINTER :: energy
TYPE(rt_prop_type), POINTER :: rtp
TYPE(rtp_control_type), POINTER :: rtp_control
TYPE(section_vals_type), POINTER :: input, rtp_section
should_stop = .FALSE.
CALL timeset(routineN, handle)
CALL cite_reference(Andermatt2016)
NULLIFY (logger, dft_control, energy, rtp, rtp_control)
NULLIFY (logger, dft_control, energy, rtp, rtp_control, input, rtp_section)
logger => cp_get_default_logger()
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, rtp=rtp, energy=energy)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, rtp=rtp, energy=energy, input=input)
rtp_control => dft_control%rtp_control
max_steps = MIN(rtp%i_start + rtp%nsteps, rtp%max_steps)
max_steps = MIN(rtp%nsteps, rtp%max_steps)
max_iter = rtp_control%max_iter
eps_ener = rtp_control%eps_ener
@ -309,7 +336,16 @@ CONTAINS
CALL cp_iterate(logger%iter_info, iter_nr=0)
IF (rtp%i_start >= max_steps) CALL cp_abort(__LOCATION__, &
"maximum step number smaller than initial step value")
rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
output_unit = cp_print_key_unit_nr(logger, rtp_section, "PRINT%PROGRAM_RUN_INFO", &
extension=".scfLog")
DO i_step = rtp%i_start + 1, max_steps
IF (output_unit > 0) THEN
WRITE (output_unit, FMT="(/,(T2,A,T40,I6))") &
"Real time propagation step:", i_step
END IF
energy%efield_core = 0.0_dp
qs_env%sim_time = REAL(i_step, dp)*rtp%dt
CALL get_qs_env(qs_env, pw_env=pw_env)

22
src/qs_core_hamiltonian.F
View file

@ -422,8 +422,8 @@ CONTAINS
INTEGER :: nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: all_present, ppl_present, ppnl_present, &
use_virial
LOGICAL :: all_present, my_gt_nl, ppl_present, &
ppnl_present, use_virial
REAL(KIND=dp) :: eps_ppnl
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dft_control_type), POINTER :: dft_control
@ -440,6 +440,16 @@ CONTAINS
NULLIFY (dft_control)
CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control)
nimages = dft_control%nimages
! check whether a gauge transformed version of the non-local potential part has to be used
my_gt_nl = .FALSE.
IF (qs_env%run_rtp) THEN
CPASSERT(ASSOCIATED(dft_control%rtp_control))
IF (dft_control%rtp_control%velocity_gauge) THEN
my_gt_nl = dft_control%rtp_control%nl_gauge_transform
END IF
END IF
! prepare for k-points
NULLIFY (cell_to_index)
IF (nimages > 1) THEN
@ -504,9 +514,11 @@ CONTAINS
eps_ppnl = dft_control%qs_control%eps_ppnl
ppnl_present = ASSOCIATED(sap_ppnl)
IF (ppnl_present) THEN
CALL build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
nimages, cell_to_index, "ORB")
IF (.NOT. my_gt_nl) THEN
CALL build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
nimages, cell_to_index, "ORB")
END IF
END IF
END SUBROUTINE core_matrices

24
src/qs_environment_types.F
View file

@ -371,13 +371,16 @@ CONTAINS
!> \param energy ...
!> \param force ...
!> \param matrix_h ...
!> \param matrix_h_im ...
!> \param matrix_ks ...
!> \param matrix_ks_im ...
!> \param matrix_vxc ...
!> \param run_rtp ...
!> \param rtp ...
!> \param matrix_h_kp ...
!> \param matrix_h_im_kp ...
!> \param matrix_ks_kp ...
!> \param matrix_ks_im_kp ...
!> \param matrix_vxc_kp ...
!> \param kinetic_kp ...
!> \param matrix_s_kp ...
@ -497,9 +500,9 @@ CONTAINS
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, particle_set, energy, force, &
matrix_h, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, &
matrix_h_kp, matrix_ks_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_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_buffer, rho_xc, pw_env, ewald_env, ewald_pw, active_space, &
mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, &
@ -542,15 +545,13 @@ CONTAINS
TYPE(qs_force_type), DIMENSION(:), OPTIONAL, &
POINTER :: force
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: matrix_h, matrix_ks, matrix_ks_im, &
matrix_vxc
POINTER :: matrix_h, matrix_h_im, matrix_ks, &
matrix_ks_im, matrix_vxc
LOGICAL, OPTIONAL :: run_rtp
TYPE(rt_prop_type), OPTIONAL, POINTER :: rtp
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
POINTER :: matrix_h_kp, matrix_ks_kp, &
matrix_vxc_kp, kinetic_kp, &
matrix_s_kp, matrix_w_kp, &
matrix_s_RI_aux_kp
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, POINTER :: 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
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: matrix_s, matrix_s_RI_aux, matrix_w, &
matrix_p_mp2, matrix_p_mp2_admm
@ -763,16 +764,19 @@ CONTAINS
potential_changed=potential_changed, &
forces_up_to_date=forces_up_to_date, &
matrix_h=matrix_h, &
matrix_h_im=matrix_h_im, &
matrix_ks=matrix_ks, &
matrix_ks_im=matrix_ks_im, &
matrix_vxc=matrix_vxc, &
kinetic=kinetic, &
matrix_s=matrix_s, &
matrix_s_RI_aux=matrix_s_RI_aux, &
matrix_ks_im_kp=matrix_ks_im_kp, &
matrix_w=matrix_w, &
matrix_p_mp2=matrix_p_mp2, &
matrix_p_mp2_admm=matrix_p_mp2_admm, &
matrix_h_kp=matrix_h_kp, &
matrix_h_im_kp=matrix_h_im_kp, &
matrix_ks_kp=matrix_ks_kp, &
matrix_vxc_kp=matrix_vxc_kp, &
kinetic_kp=kinetic_kp, &

23
src/qs_ks_methods.F
View file

@ -195,8 +195,8 @@ CONTAINS
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_vxc, mo_derivs
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h, matrix_s, my_rho, &
rho_ao
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, ks_matrix_im, matrix_h, &
matrix_h_im, matrix_s, my_rho, rho_ao
TYPE(dft_control_type), POINTER :: dft_control
TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
TYPE(local_rho_type), POINTER :: local_rho_set
@ -224,10 +224,10 @@ CONTAINS
NULLIFY (admm_env, cell, dft_control, logger, mo_derivs, my_rho, &
rho_struct, para_env, pw_env, virial, vppl_rspace, &
adiabatic_rescaling_section, hfx_sections, &
input, scf_section, xc_section, matrix_h, matrix_s, &
input, scf_section, xc_section, matrix_h, matrix_h_im, matrix_s, &
auxbas_pw_pool, poisson_env, v_rspace_new, v_rspace_new_aux_fit, &
v_tau_rspace, v_tau_rspace_aux_fit, matrix_vxc, vee, rho_nlcc, &
ks_env, ks_matrix, rho, energy, rho_xc, rho_r, rho_ao, rho_core)
ks_env, ks_matrix, ks_matrix_im, rho, energy, rho_xc, rho_r, rho_ao, rho_core)
CPASSERT(ASSOCIATED(qs_env))
@ -239,8 +239,10 @@ CONTAINS
ks_env=ks_env, &
dft_control=dft_control, &
matrix_h_kp=matrix_h, &
matrix_h_im_kp=matrix_h_im, &
matrix_s_kp=matrix_s, &
matrix_ks_kp=ks_matrix, &
matrix_ks_im_kp=ks_matrix_im, &
matrix_vxc=matrix_vxc, &
pw_env=pw_env, &
cell=cell, &
@ -666,6 +668,19 @@ CONTAINS
CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix, name=name)
END DO
END DO
! imaginary part if required
IF (qs_env%run_rtp) THEN
IF (dft_control%rtp_control%velocity_gauge) THEN
CPASSERT(ASSOCIATED(matrix_h_im))
CPASSERT(ASSOCIATED(ks_matrix_im))
DO ispin = 1, nspins
DO img = 1, nimages
CALL dbcsr_get_info(ks_matrix_im(ispin, img)%matrix, name=name) ! keep the name
CALL dbcsr_copy(ks_matrix_im(ispin, img)%matrix, matrix_h_im(1, img)%matrix, name=name)
END DO
END DO
END IF
END IF
END IF
IF (use_virial .AND. calculate_forces) THEN

71
src/qs_ks_types.F
View file

@ -155,10 +155,11 @@ MODULE qs_ks_types
matrix_h, &
matrix_w, &
matrix_vxc, &
matrix_s_RI_aux
matrix_s_RI_aux, &
matrix_ks_im, &
matrix_h_im
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_im => Null(), &
matrix_p_mp2 => Null(), &
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p_mp2 => Null(), &
matrix_p_mp2_admm => Null()
TYPE(qs_rho_type), POINTER :: rho => Null(), &
@ -242,6 +243,7 @@ CONTAINS
!> \param potential_changed ...
!> \param forces_up_to_date ...
!> \param matrix_h ...
!> \param matrix_h_im ...
!> \param matrix_ks ...
!> \param matrix_ks_im ...
!> \param matrix_vxc ...
@ -252,12 +254,14 @@ CONTAINS
!> \param matrix_p_mp2 ...
!> \param matrix_p_mp2_admm ...
!> \param matrix_h_kp ...
!> \param matrix_h_im_kp ...
!> \param matrix_ks_kp ...
!> \param matrix_vxc_kp ...
!> \param kinetic_kp ...
!> \param matrix_s_kp ...
!> \param matrix_w_kp ...
!> \param matrix_s_RI_aux_kp ...
!> \param matrix_ks_im_kp ...
!> \param rho ...
!> \param rho_buffer ...
!> \param rho_xc ...
@ -320,11 +324,11 @@ CONTAINS
SUBROUTINE get_ks_env(ks_env, v_hartree_rspace, &
s_mstruct_changed, rho_changed, &
potential_changed, forces_up_to_date, &
matrix_h, matrix_ks, matrix_ks_im, matrix_vxc, &
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_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, &
matrix_s_RI_aux_kp, &
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, &
rho, &
rho_buffer, rho_xc, &
vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, &
@ -343,13 +347,12 @@ CONTAINS
TYPE(pw_type), OPTIONAL, POINTER :: v_hartree_rspace
LOGICAL, OPTIONAL :: s_mstruct_changed, rho_changed, &
potential_changed, forces_up_to_date
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, POINTER :: matrix_h, matrix_ks, matrix_ks_im, &
matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
POINTER :: matrix_h_kp, matrix_ks_kp, &
matrix_vxc_kp, kinetic_kp, &
matrix_s_kp, matrix_w_kp, &
matrix_s_RI_aux_kp
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, POINTER :: 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
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, POINTER :: 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
TYPE(qs_rho_type), OPTIONAL, POINTER :: rho, rho_buffer, rho_xc
TYPE(pw_type), OPTIONAL, POINTER :: vppl, rho_core, rho_nlcc, rho_nlcc_g, vee
INTEGER, OPTIONAL :: neighbor_list_id
@ -404,18 +407,21 @@ CONTAINS
IF (PRESENT(matrix_s_RI_aux)) matrix_s_RI_aux => get_1d_pointer(ks_env%matrix_s_RI_aux)
IF (PRESENT(kinetic)) kinetic => get_1d_pointer(ks_env%kinetic)
IF (PRESENT(matrix_h)) matrix_h => get_1d_pointer(ks_env%matrix_h)
IF (PRESENT(matrix_h_im)) matrix_h_im => get_1d_pointer(ks_env%matrix_h_im)
IF (PRESENT(matrix_w)) matrix_w => get_1d_pointer(ks_env%matrix_w)
IF (PRESENT(matrix_vxc)) matrix_vxc => get_1d_pointer(ks_env%matrix_vxc)
IF (PRESENT(matrix_ks_kp)) matrix_ks_kp => get_2d_pointer(ks_env%matrix_ks)
IF (PRESENT(matrix_ks_im_kp)) matrix_ks_im_kp => get_2d_pointer(ks_env%matrix_ks_im)
IF (PRESENT(matrix_s_kp)) matrix_s_kp => get_2d_pointer(ks_env%matrix_s)
IF (PRESENT(matrix_s_RI_aux_kp)) matrix_s_RI_aux_kp => get_2d_pointer(ks_env%matrix_s_RI_aux)
IF (PRESENT(matrix_w_kp)) matrix_w_kp => get_2d_pointer(ks_env%matrix_w)
IF (PRESENT(kinetic_kp)) kinetic_kp => get_2d_pointer(ks_env%kinetic)
IF (PRESENT(matrix_h_kp)) matrix_h_kp => get_2d_pointer(ks_env%matrix_h)
IF (PRESENT(matrix_h_im_kp)) matrix_h_im_kp => get_2d_pointer(ks_env%matrix_h_im)
IF (PRESENT(matrix_vxc_kp)) matrix_vxc_kp => get_2d_pointer(ks_env%matrix_vxc)
IF (PRESENT(matrix_ks_im)) matrix_ks_im => get_1d_pointer(ks_env%matrix_ks_im)
IF (PRESENT(matrix_ks_im)) matrix_ks_im => ks_env%matrix_ks_im
IF (PRESENT(matrix_p_mp2)) matrix_p_mp2 => ks_env%matrix_p_mp2
IF (PRESENT(matrix_p_mp2_admm)) matrix_p_mp2_admm => ks_env%matrix_p_mp2_admm
IF (PRESENT(rho)) rho => ks_env%rho
@ -495,6 +501,7 @@ CONTAINS
!> \param potential_changed ...
!> \param forces_up_to_date ...
!> \param matrix_h ...
!> \param matrix_h_im ...
!> \param matrix_ks ...
!> \param matrix_ks_im ...
!> \param matrix_vxc ...
@ -505,12 +512,14 @@ CONTAINS
!> \param matrix_p_mp2 ...
!> \param matrix_p_mp2_admm ...
!> \param matrix_h_kp ...
!> \param matrix_h_im_kp ...
!> \param matrix_ks_kp ...
!> \param matrix_vxc_kp ...
!> \param kinetic_kp ...
!> \param matrix_s_kp ...
!> \param matrix_w_kp ...
!> \param matrix_s_RI_aux_kp ...
!> \param matrix_ks_im_kp ...
!> \param vppl ...
!> \param rho_core ...
!> \param rho_nlcc ...
@ -549,11 +558,11 @@ CONTAINS
SUBROUTINE set_ks_env(ks_env, v_hartree_rspace, &
s_mstruct_changed, rho_changed, &
potential_changed, forces_up_to_date, &
matrix_h, matrix_ks, matrix_ks_im, matrix_vxc, &
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_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, &
matrix_s_RI_aux_kp, &
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, &
@ -567,13 +576,12 @@ CONTAINS
TYPE(pw_type), OPTIONAL, POINTER :: v_hartree_rspace
LOGICAL, OPTIONAL :: s_mstruct_changed, rho_changed, &
potential_changed, forces_up_to_date
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, POINTER :: matrix_h, matrix_ks, matrix_ks_im, &
matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
POINTER :: matrix_h_kp, matrix_ks_kp, &
matrix_vxc_kp, kinetic_kp, &
matrix_s_kp, matrix_w_kp, &
matrix_s_RI_aux_kp
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, POINTER :: 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
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, POINTER :: 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
TYPE(pw_type), OPTIONAL, POINTER :: vppl, rho_core, rho_nlcc, rho_nlcc_g, vee
INTEGER, OPTIONAL :: neighbor_list_id
TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
@ -602,6 +610,9 @@ CONTAINS
IF (PRESENT(matrix_h)) CALL set_1d_pointer(ks_env%matrix_h, matrix_h)
IF (PRESENT(matrix_h_kp)) CALL set_2d_pointer(ks_env%matrix_h, matrix_h_kp)
IF (PRESENT(matrix_h_im)) CALL set_1d_pointer(ks_env%matrix_h_im, matrix_h_im)
IF (PRESENT(matrix_h_im_kp)) CALL set_2d_pointer(ks_env%matrix_h_im, matrix_h_im_kp)
IF (PRESENT(matrix_ks)) CALL set_1d_pointer(ks_env%matrix_ks, matrix_ks)
IF (PRESENT(matrix_ks_kp)) CALL set_2d_pointer(ks_env%matrix_ks, matrix_ks_kp)
@ -620,7 +631,9 @@ CONTAINS
IF (PRESENT(matrix_s_RI_aux)) CALL set_1d_pointer(ks_env%matrix_s_RI_aux, matrix_s_RI_aux)
IF (PRESENT(matrix_s_RI_aux_kp)) CALL set_2d_pointer(ks_env%matrix_s_RI_aux, matrix_s_RI_aux_kp)
IF (PRESENT(matrix_ks_im)) ks_env%matrix_ks_im => matrix_ks_im
IF (PRESENT(matrix_ks_im)) CALL set_1d_pointer(ks_env%matrix_ks_im, matrix_ks_im)
IF (PRESENT(matrix_ks_im_kp)) CALL set_2d_pointer(ks_env%matrix_ks_im, matrix_ks_im_kp)
IF (PRESENT(matrix_p_mp2)) ks_env%matrix_p_mp2 => matrix_p_mp2
IF (PRESENT(matrix_p_mp2_admm)) ks_env%matrix_p_mp2_admm => matrix_p_mp2_admm
IF (PRESENT(rho_core)) ks_env%rho_core => rho_core
@ -711,11 +724,11 @@ CONTAINS
CALL pw_release(ks_env%v_hartree_rspace)
DEALLOCATE (ks_env%v_hartree_rspace)
END IF
IF (ASSOCIATED(ks_env%matrix_ks_im)) &
CALL dbcsr_deallocate_matrix_set(ks_env%matrix_ks_im)
CALL kpoint_transitional_release(ks_env%matrix_ks)
CALL kpoint_transitional_release(ks_env%matrix_ks_im)
CALL kpoint_transitional_release(ks_env%matrix_h)
CALL kpoint_transitional_release(ks_env%matrix_h_im)
CALL kpoint_transitional_release(ks_env%matrix_vxc)
CALL kpoint_transitional_release(ks_env%matrix_s)
CALL kpoint_transitional_release(ks_env%matrix_w)
@ -812,10 +825,10 @@ CONTAINS
CALL pw_release(ks_env%v_hartree_rspace)
DEALLOCATE (ks_env%v_hartree_rspace)
END IF
IF (ASSOCIATED(ks_env%matrix_ks_im)) &
CALL dbcsr_deallocate_matrix_set(ks_env%matrix_ks_im)
CALL kpoint_transitional_release(ks_env%matrix_h)
CALL kpoint_transitional_release(ks_env%matrix_h_im)
CALL kpoint_transitional_release(ks_env%matrix_ks_im)
CALL kpoint_transitional_release(ks_env%matrix_vxc)
CALL kpoint_transitional_release(ks_env%matrix_w)
CALL kpoint_transitional_release(ks_env%kinetic)

182
src/qs_moments.F
View file

@ -22,6 +22,8 @@ MODULE qs_moments
get_atomic_kind
USE basis_set_types, ONLY: gto_basis_set_p_type,&
gto_basis_set_type
USE bibliography, ONLY: Mattiat2019,&
cite_reference
USE block_p_types, ONLY: block_p_type
USE cell_types, ONLY: cell_type,&
pbc
@ -51,10 +53,9 @@ MODULE qs_moments
put_results
USE cp_result_types, ONLY: cp_result_type
USE dbcsr_api, ONLY: &
dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
dbcsr_distribution_type, dbcsr_dot, dbcsr_get_block_p, dbcsr_init_p, dbcsr_multiply, &
dbcsr_p_type, dbcsr_set, dbcsr_trace, dbcsr_type, dbcsr_type_antisymmetric, &
dbcsr_type_no_symmetry, dbcsr_type_symmetric
dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_dot, &
dbcsr_get_block_p, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_set, dbcsr_trace, &
dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
USE distribution_1d_types, ONLY: distribution_1d_type
USE kinds, ONLY: default_string_length,&
dp
@ -102,7 +103,7 @@ MODULE qs_moments
! Public subroutines
PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix
PUBLIC :: build_berry_kpoint_matrix
PUBLIC :: build_berry_kpoint_matrix, build_local_magmom_matrix
PUBLIC :: qs_moment_berry_phase, qs_moment_locop
PUBLIC :: dipole_deriv_ao
@ -610,6 +611,8 @@ CONTAINS
DEALLOCATE (work)
DO i = 1, nm
NULLIFY (mom_block(i)%block)
END DO
DO i = 1, M_dim
DO idir = 1, 3
NULLIFY (mom_block_der(i, idir)%block)
END DO
@ -642,8 +645,8 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_magmom_matrix'
INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, maxco, &
maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, npgfa, npgfb, &
nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
@ -662,7 +665,7 @@ CONTAINS
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_all
POINTER :: sab_orb
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
@ -671,26 +674,19 @@ CONTAINS
CALL timeset(routineN, handle)
NULLIFY (qs_kind_set, cell, particle_set, sab_all)
NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
NULLIFY (matrix_s)
CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
! magnetic dipoles/angular moments only
nm = 3
CALL dbcsr_allocate_matrix_set(magmom, nm)
DO i = 1, nm
ALLOCATE (magmom(i)%matrix)
CALL dbcsr_desymmetrize(matrix_s(1)%matrix, magmom(i)%matrix)
CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
END DO
NULLIFY (qs_kind_set, particle_set, sab_all, cell)
NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
CALL get_qs_env(qs_env=qs_env, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, cell=cell, &
sab_all=sab_all)
sab_orb=sab_orb)
nkind = SIZE(qs_kind_set)
natom = SIZE(particle_set)
@ -720,7 +716,7 @@ CONTAINS
NULLIFY (basis_set_list(ikind)%gto_basis_set)
END IF
END DO
CALL neighbor_list_iterator_create(nl_iterator, sab_all)
CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
iatom=iatom, jatom=jatom, r=rab)
@ -750,19 +746,14 @@ CONTAINS
sphi_b => basis_set_b%sphi
zetb => basis_set_b%zet
IF (inode == 1) last_jatom = 0
! this guarantees minimum image convention
! anything else would not make sense
IF (jatom == last_jatom) THEN
CYCLE
IF (iatom <= jatom) THEN
irow = iatom
icol = jatom
ELSE
irow = jatom
icol = iatom
END IF
last_jatom = jatom
irow = iatom
icol = jatom
DO i = 1, nm
NULLIFY (mint(i)%block)
CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
@ -809,17 +800,24 @@ CONTAINS
! Contraction step
DO i = 1, nm
CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
sphi_b(1, sgfb), SIZE(sphi_b, 1), &
0.0_dp, work(1, 1), SIZE(work, 1))
CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
work(1, 1), SIZE(work, 1), &
1.0_dp, mint(i)%block(sgfa, sgfb), &
SIZE(mint(i)%block, 1))
IF (iatom <= jatom) THEN
CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
work(1, 1), SIZE(work, 1), &
1.0_dp, mint(i)%block(sgfa, sgfb), &
SIZE(mint(i)%block, 1))
ELSE
CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
-1.0_dp, work(1, 1), SIZE(work, 1), &
sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1.0_dp, mint(i)%block(sgfb, sgfa), &
SIZE(mint(i)%block, 1))
END IF
END DO
@ -1773,8 +1771,8 @@ CONTAINS
CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
CHARACTER(LEN=default_string_length) :: description
INTEGER :: akind, handle, i, ia, iatom, idir, &
ikind, ispin, ix, iy, iz, l, nm, nmm, &
nmom, order
ikind, ispin, ix, iy, iz, l, nm, nmom, &
order
LOGICAL :: my_com_nl, my_velreprs
REAL(dp) :: charge, dd, strace, trace
REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rv, nlcom_rxrv, &
@ -1792,7 +1790,7 @@ CONTAINS
TYPE(dft_control_type), POINTER :: dft_control
TYPE(distribution_1d_type), POINTER :: local_particles
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb
POINTER :: sab_all, sab_orb
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_rho_type), POINTER :: rho
@ -1804,6 +1802,7 @@ CONTAINS
my_velreprs = .FALSE.
IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
IF (PRESENT(com_nl)) my_com_nl = com_nl
IF (my_velreprs) CALL cite_reference(Mattiat2019)
! reference point
CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
@ -1811,7 +1810,7 @@ CONTAINS
! only allow for moments up to maxl set by basis
nmom = MIN(nmoments, current_maxl)
! electronic contribution
NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_orb)
NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_all, sab_orb)
CALL get_qs_env(qs_env, &
dft_control=dft_control, &
rho=rho, &
@ -1821,7 +1820,8 @@ CONTAINS
qs_kind_set=qs_kind_set, &
para_env=para_env, &
matrix_s=matrix_s, &
sab_all=sab_orb)
sab_all=sab_all, &
sab_orb=sab_orb)
IF (my_com_nl) THEN
IF ((nmom >= 1) .AND. my_velreprs) THEN
@ -1848,7 +1848,7 @@ CONTAINS
IF (my_velreprs .AND. (nmom >= 2)) THEN
CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
matrix_type=dbcsr_type_no_symmetry)
CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_all)
ELSE
CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
END IF
@ -1864,7 +1864,7 @@ CONTAINS
CALL dbcsr_init_p(moments_der(i, idir)%matrix)
CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
matrix_type=dbcsr_type_no_symmetry)
CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_all)
CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
END DO
END DO
@ -1884,8 +1884,10 @@ CONTAINS
! Allocate matrix to store the matrix product to be traced (dbcsr_dot only works for products of
! symmetric matrices)
NULLIFY (tmp_ao)
ALLOCATE (tmp_ao)
CALL dbcsr_desymmetrize(matrix_s(1)%matrix, tmp_ao)
CALL dbcsr_init_p(tmp_ao)
CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
CALL dbcsr_set(tmp_ao, 0.0_dp)
END IF
trace = 0.0_dp
@ -1944,10 +1946,19 @@ CONTAINS
! magnetic moments
IF (magnetic) THEN
NULLIFY (magmom)
CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
nmm = SIZE(magmom)
ALLOCATE (mmom(nmm))
CALL dbcsr_allocate_matrix_set(magmom, 3)
DO i = 1, SIZE(magmom)
CALL dbcsr_init_p(magmom(i)%matrix)
CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
matrix_type=dbcsr_type_antisymmetric)
CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
END DO
CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
ALLOCATE (mmom(SIZE(magmom)))
mmom(:) = 0.0_dp
IF (qs_env%run_rtp) THEN
! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
@ -1980,8 +1991,6 @@ CONTAINS
SELECT CASE (order)
CASE (1) ! expectation value of momentum
NULLIFY (sab_orb)
CALL get_qs_env(qs_env, sab_orb=sab_orb)
NULLIFY (momentum)
CALL dbcsr_allocate_matrix_set(momentum, 3)
DO i = 1, 3
@ -2092,13 +2101,13 @@ CONTAINS
IF (my_com_nl) THEN
IF (magnetic) THEN
mmom(:) = mmom(:) + nlcom_rxrv(:)
mmom(:) = nlcom_rxrv(:)
END IF
IF (my_velreprs .AND. (nmom >= 1)) THEN
rmom_vel(1:3) = rmom_vel(1:3) + nlcom_rv
rmom_vel(1:3) = nlcom_rv
END IF
IF (my_velreprs .AND. (nmom >= 2)) THEN
rmom_vel(4:9) = rmom_vel(4:9) + nlcom_rrv
rmom_vel(4:9) = nlcom_rrv
END IF
IF (magnetic .AND. .NOT. my_velreprs) THEN
CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom)
@ -2199,13 +2208,13 @@ CONTAINS
END IF
dd = SQRT(SUM(rmom(2:4, 3)**2))*debye
WRITE (unit_number, "(T3,A)") "Dipole moment [Debye]"
WRITE (unit_number, "(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
(TRIM(rlab(i)), "=", rmom(i, 3)*debye, i=2, 4), "Total=", dd
CASE (2)
WRITE (unit_number, "(T3,A)") "Quadrupole moment [Debye*Angstrom]"
WRITE (unit_number, "(T17,3(A,A,F14.8,9X))") &
WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
(TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=5, 7)
WRITE (unit_number, "(T17,3(A,A,F14.8,9X))") &
WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
(TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=8, 10)
CASE (3)
WRITE (unit_number, "(T3,A)") "Octapole moment [Debye*Angstrom**2]"
@ -2241,7 +2250,7 @@ CONTAINS
IF (nmom >= 1) THEN
dd = SQRT(SUM(mmom(1:3)**2))
WRITE (unit_number, "(T3,A)") "Orbital angular momentum [a. u.]"
WRITE (unit_number, "(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
(TRIM(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
END IF
END IF
@ -2251,13 +2260,13 @@ CONTAINS
CASE (1)
dd = SQRT(SUM(rmom_vel(1:3)**2))
WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator [a. u.]"
WRITE (unit_number, "(T5,3(A,A,D14.8,1X),T60,A,T67,F14.8)") &
WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
(TRIM(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
CASE (2)
WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. [a. u.]"
WRITE (unit_number, "(T17,3(A,A,D14.8,9X))") &
WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
(TRIM(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
WRITE (unit_number, "(T17,3(A,A,D14.8,9X))") &
WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
(TRIM(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
CASE DEFAULT
END SELECT
@ -2287,8 +2296,8 @@ CONTAINS
IF (PRESENT(mmom)) THEN
IF (nmom >= 1) THEN
dd = SQRT(SUM(mmom(1:3)**2))
WRITE (unit_number, "(T3,A)") "Orbital angular momentum incl. com_nl [a. u.]"
WRITE (unit_number, "(T5,3(A,A,D14.8,1X),T60,A,T67,F14.8)") &
WRITE (unit_number, "(T3,A)") "Expectation value of rx[r,V_nl] [a. u.]"
WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
(TRIM(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
END IF
END IF
@ -2297,14 +2306,14 @@ CONTAINS
SELECT CASE (l)
CASE (1)
dd = SQRT(SUM(rmom_vel(1:3)**2))
WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator incl. com_nl [a. u.]"
WRITE (unit_number, "(T5,3(A,A,D14.8,1X),T60,A,T67,F14.8)") &
WRITE (unit_number, "(T3,A)") "Expectation value of [r,V_nl] [a. u.]"
WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
(TRIM(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
CASE (2)
WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. incl. com_nl [a. u.]"
WRITE (unit_number, "(T17,3(A,A,D14.8,9X))") &
WRITE (unit_number, "(T3,A)") "Expectation value of [rr,V_nl] [a. u.]"
WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
(TRIM(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
WRITE (unit_number, "(T17,3(A,A,D14.8,9X))") &
WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
(TRIM(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
CASE DEFAULT
END SELECT
@ -2339,7 +2348,7 @@ CONTAINS
TYPE(dbcsr_type), POINTER :: tmp_ao
TYPE(dft_control_type), POINTER :: dft_control
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sap_ppnl
POINTER :: sab_all, sab_orb, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_rho_type), POINTER :: rho
@ -2374,7 +2383,7 @@ CONTAINS
! imagninary part of the density matrix
IF (qs_env%run_rtp) THEN
NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_orb, sap_ppnl)
NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_all, sab_orb, sap_ppnl)
CALL get_qs_env(qs_env, &
cell=cell, &
dft_control=dft_control, &
@ -2383,6 +2392,7 @@ CONTAINS
qs_kind_set=qs_kind_set, &
rho=rho, &
sab_orb=sab_orb, &
sab_all=sab_all, &
sap_ppnl=sap_ppnl)
eps_ppnl = dft_control%qs_control%eps_ppnl
@ -2424,34 +2434,36 @@ CONTAINS
! calculate commutators in AO basis
IF (calc_rv .AND. calc_rrv .AND. .NOT. calc_rxrv) THEN
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_rv=matrix_rv, &
matrix_rrv=matrix_rrv, ref_point=ref_point, cell=cell)
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
matrix_rrv=matrix_rrv, ref_point=ref_point)
ELSE IF (calc_rv .AND. calc_rxrv .AND. .NOT. calc_rrv) THEN
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_rv=matrix_rv, &
matrix_rxrv=matrix_rxrv, ref_point=ref_point, cell=cell)
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
matrix_rxrv=matrix_rxrv, ref_point=ref_point)
ELSE IF (calc_rrv .AND. calc_rxrv .AND. .NOT. calc_rv) THEN
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_rxrv=matrix_rxrv, &
matrix_rrv=matrix_rrv, ref_point=ref_point, cell=cell)
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rxrv=matrix_rxrv, &
matrix_rrv=matrix_rrv, ref_point=ref_point)
ELSE IF (calc_rv .AND. calc_rxrv .AND. calc_rrv) THEN
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_rv=matrix_rv, &
matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, ref_point=ref_point, cell=cell)
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, ref_point=ref_point)
ELSE IF (calc_rv .AND. .NOT. calc_rrv .AND. .NOT. calc_rxrv) THEN
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_rv=matrix_rv, &
ref_point=ref_point, cell=cell)
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
ref_point=ref_point)
ELSE IF (calc_rrv .AND. .NOT. calc_rv .AND. .NOT. calc_rxrv) THEN
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_rrv=matrix_rrv, &
ref_point=ref_point, cell=cell)
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rrv=matrix_rrv, &
ref_point=ref_point)
ELSE IF (calc_rxrv .AND. .NOT. calc_rv .AND. .NOT. calc_rrv) THEN
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_rxrv=matrix_rxrv, &
ref_point=ref_point, cell=cell)
CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rxrv=matrix_rxrv, &
ref_point=ref_point)
END IF
NULLIFY (rho_ao)
CALL qs_rho_get(rho, rho_ao_im=rho_ao)
NULLIFY (tmp_ao)
ALLOCATE (tmp_ao)
CALL dbcsr_desymmetrize(matrix_s(1)%matrix, tmp_ao)
CALL dbcsr_init_p(tmp_ao)
CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
CALL dbcsr_set(tmp_ao, 0.0_dp)
IF (calc_rv) THEN
trace = 0._dp

25
src/qs_rho_types.F
View file

@ -16,7 +16,6 @@
!> \author Fawzi Mohamed
! **************************************************************************************************
MODULE qs_rho_types
USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
USE dbcsr_api, ONLY: dbcsr_p_type
USE kinds, ONLY: dp
USE kpoint_transitional, ONLY: get_1d_pointer,&
@ -69,7 +68,7 @@ MODULE qs_rho_types
TYPE qs_rho_type
PRIVATE
TYPE(kpoint_transitional_type) :: rho_ao
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_im => Null()
TYPE(kpoint_transitional_type) :: rho_ao_im
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g => Null(), &
rho_r => Null(), &
tau_g => Null(), &
@ -219,8 +218,8 @@ CONTAINS
CALL kpoint_transitional_release(rho_struct%rho_ao)
IF (ASSOCIATED(rho_struct%rho_ao_im)) &
CALL dbcsr_deallocate_matrix_set(rho_struct%rho_ao_im)
CALL kpoint_transitional_release(rho_struct%rho_ao_im)
IF (ASSOCIATED(rho_struct%tot_rho_r)) &
DEALLOCATE (rho_struct%tot_rho_r)
IF (ASSOCIATED(rho_struct%tot_rho_g)) &
@ -246,6 +245,7 @@ CONTAINS
!> \param rho_ao ...
!> \param rho_ao_im ...
!> \param rho_ao_kp ...
!> \param rho_ao_im_kp ...
!> \param rho_r ...
!> \param drho_r ...
!> \param rho_g ...
@ -267,7 +267,7 @@ CONTAINS
!> 08.2002 created [fawzi]
!> \author Fawzi Mohamed
! **************************************************************************************************
SUBROUTINE qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_r, drho_r, &
SUBROUTINE qs_rho_get(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, rebuild_each, tot_rho_r, tot_rho_g, &
rho_r_sccs, soft_valid)
@ -275,7 +275,7 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: rho_ao, rho_ao_im
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
POINTER :: rho_ao_kp
POINTER :: rho_ao_kp, rho_ao_im_kp
TYPE(pw_p_type), DIMENSION(:), OPTIONAL, POINTER :: rho_r
TYPE(pw_type), DIMENSION(:, :), OPTIONAL, POINTER :: drho_r
TYPE(pw_p_type), DIMENSION(:), OPTIONAL, POINTER :: rho_g
@ -294,7 +294,9 @@ CONTAINS
IF (PRESENT(rho_ao)) rho_ao => get_1d_pointer(rho_struct%rho_ao)
IF (PRESENT(rho_ao_kp)) rho_ao_kp => get_2d_pointer(rho_struct%rho_ao)
IF (PRESENT(rho_ao_im)) rho_ao_im => rho_struct%rho_ao_im
IF (PRESENT(rho_ao_im)) rho_ao_im => get_1d_pointer(rho_struct%rho_ao_im)
IF (PRESENT(rho_ao_im_kp)) rho_ao_im_kp => get_2d_pointer(rho_struct%rho_ao_im)
IF (PRESENT(rho_r)) rho_r => rho_struct%rho_r
IF (PRESENT(drho_r)) drho_r => rho_struct%drho_r
IF (PRESENT(rho_g)) rho_g => rho_struct%rho_g
@ -321,6 +323,7 @@ CONTAINS
!> \param rho_ao ...
!> \param rho_ao_im ...
!> \param rho_ao_kp ...
!> \param rho_ao_im_kp ...
!> \param rho_r ...
!> \param drho_r ...
!> \param rho_g ...
@ -340,7 +343,7 @@ CONTAINS
!> \param soft_valid ...
!> \author Ole Schuett
! **************************************************************************************************
SUBROUTINE qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_r, drho_r, &
SUBROUTINE 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, rebuild_each, tot_rho_r, tot_rho_g, &
rho_r_sccs, soft_valid)
@ -348,7 +351,7 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: rho_ao, rho_ao_im
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
POINTER :: rho_ao_kp
POINTER :: rho_ao_kp, rho_ao_im_kp
TYPE(pw_p_type), DIMENSION(:), OPTIONAL, POINTER :: rho_r
TYPE(pw_type), DIMENSION(:, :), OPTIONAL, POINTER :: drho_r
TYPE(pw_p_type), DIMENSION(:), OPTIONAL, POINTER :: rho_g
@ -367,7 +370,9 @@ CONTAINS
IF (PRESENT(rho_ao)) CALL set_1d_pointer(rho_struct%rho_ao, rho_ao)
IF (PRESENT(rho_ao_kp)) CALL set_2d_pointer(rho_struct%rho_ao, rho_ao_kp)
IF (PRESENT(rho_ao_im)) rho_struct%rho_ao_im => rho_ao_im
IF (PRESENT(rho_ao_im)) CALL set_1d_pointer(rho_struct%rho_ao_im, rho_ao_im)
IF (PRESENT(rho_ao_im_kp)) CALL set_2d_pointer(rho_struct%rho_ao_im, rho_ao_im_kp)
IF (PRESENT(rho_r)) rho_struct%rho_r => rho_r
IF (PRESENT(rho_g)) rho_struct%rho_g => rho_g
IF (PRESENT(drho_r)) rho_struct%drho_r => drho_r

832
src/rt_propagation_velocity_gauge.F Normal file
View file

@ -0,0 +1,832 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2022 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Routines to perform the RTP in the velocity gauge
! **************************************************************************************************
MODULE rt_propagation_velocity_gauge
USE ai_moments, ONLY: cossin
USE atomic_kind_types, ONLY: atomic_kind_type
USE basis_set_types, ONLY: gto_basis_set_p_type,&
gto_basis_set_type
USE bibliography, ONLY: Mattiat2022,&
cite_reference
USE cell_types, ONLY: cell_type,&
pbc
USE core_ppnl, ONLY: build_core_ppnl
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE dbcsr_api, ONLY: dbcsr_add,&
dbcsr_create,&
dbcsr_get_block_p,&
dbcsr_init_p,&
dbcsr_p_type,&
dbcsr_set,&
dbcsr_type_antisymmetric,&
dbcsr_type_symmetric
USE external_potential_types, ONLY: gth_potential_p_type,&
gth_potential_type,&
sgp_potential_p_type,&
sgp_potential_type
USE input_section_types, ONLY: section_vals_type
USE kinds, ONLY: dp,&
int_8
USE kpoint_types, ONLY: get_kpoint_info,&
kpoint_type
USE mathconstants, ONLY: one,&
twopi,&
zero
USE orbital_pointers, ONLY: init_orbital_pointers,&
nco,&
ncoset
USE particle_types, ONLY: particle_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
USE qs_ks_types, ONLY: get_ks_env,&
qs_ks_env_type
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE qs_operators_ao, ONLY: build_lin_mom_matrix
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE sap_kind_types, ONLY: alist_type,&
clist_type,&
get_alist,&
release_sap_int,&
sap_int_type,&
sap_sort
USE virial_types, ONLY: virial_type
!$ USE OMP_LIB, ONLY: omp_lock_kind, &
!$ omp_init_lock, omp_set_lock, &
!$ omp_unset_lock, omp_destroy_lock
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_velocity_gauge'
PUBLIC :: velocity_gauge_ks_matrix
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param subtract_nl_term ...
! **************************************************************************************************
SUBROUTINE velocity_gauge_ks_matrix(qs_env, subtract_nl_term)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN), OPTIONAL :: subtract_nl_term
CHARACTER(len=*), PARAMETER :: routineN = 'velocity_gauge_ks_matrix'
INTEGER :: handle, idir, image, nder, nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: calculate_forces, my_subtract_nl_term, &
ppnl_present, use_virial
REAL(KIND=dp) :: eps_ppnl, factor
REAL(KIND=dp), DIMENSION(3) :: vec_pot
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: momentum, nl_term
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im, matrix_nl, &
matrix_p, matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kpoint_type), POINTER :: kpoints
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(section_vals_type), POINTER :: input
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
CALL cite_reference(Mattiat2022)
my_subtract_nl_term = .FALSE.
IF (PRESENT(subtract_nl_term)) my_subtract_nl_term = subtract_nl_term
NULLIFY (dft_control, matrix_s, sab_orb, matrix_h, cell, input, matrix_h_im, kpoints, cell_to_index, &
sap_ppnl, particle_set, qs_kind_set, atomic_kind_set, virial, force, matrix_p, rho, matrix_nl)
CALL get_qs_env(qs_env, &
dft_control=dft_control, &
sab_orb=sab_orb, &
sap_ppnl=sap_ppnl, &
matrix_s_kp=matrix_s, &
matrix_h_kp=matrix_h, &
cell=cell, &
input=input, &
matrix_h_im_kp=matrix_h_im)
nimages = dft_control%nimages
ppnl_present = ASSOCIATED(sap_ppnl)
IF (nimages > 1) THEN
CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
END IF
IF (my_subtract_nl_term) THEN
IF (ppnl_present) THEN
CALL get_qs_env(qs_env, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, &
atomic_kind_set=atomic_kind_set, &
virial=virial, &
rho=rho, &
force=force)
CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
calculate_forces = .FALSE.
use_virial = .FALSE.
nder = 1
eps_ppnl = dft_control%qs_control%eps_ppnl
CALL dbcsr_allocate_matrix_set(matrix_nl, 1, nimages)
DO image = 1, nimages
ALLOCATE (matrix_nl(1, image)%matrix)
CALL dbcsr_create(matrix_nl(1, image)%matrix, template=matrix_s(1, 1)%matrix)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_nl(1, image)%matrix, sab_orb)
CALL dbcsr_set(matrix_nl(1, image)%matrix, zero)
END DO
CALL build_core_ppnl(matrix_nl, matrix_p, force, virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
nimages, cell_to_index, "ORB")
DO image = 1, nimages
CALL dbcsr_add(matrix_h(1, image)%matrix, matrix_nl(1, image)%matrix, one, -one)
END DO
CALL dbcsr_deallocate_matrix_set(matrix_nl)
END IF
END IF
!get vector potential
vec_pot(:) = 0._dp
CALL get_vector_potential(dft_control, cell, vec_pot)
! allocate and build matrices for linear momentum term
NULLIFY (momentum)
CALL dbcsr_allocate_matrix_set(momentum, 3)
DO idir = 1, 3
CALL dbcsr_init_p(momentum(idir)%matrix)
CALL dbcsr_create(momentum(idir)%matrix, template=matrix_s(1, 1)%matrix, &
matrix_type=dbcsr_type_antisymmetric)
CALL cp_dbcsr_alloc_block_from_nbl(momentum(idir)%matrix, sab_orb)
CALL dbcsr_set(momentum(idir)%matrix, zero)
END DO
CALL build_lin_mom_matrix(qs_env, momentum)
! set imaginary part of KS matrix to zero
DO image = 1, nimages
CALL dbcsr_set(matrix_h_im(1, image)%matrix, zero)
END DO
! add linear term in vector potential to imaginary part of KS-matrix
DO image = 1, nimages
DO idir = 1, 3
CALL dbcsr_add(matrix_h_im(1, image)%matrix, momentum(idir)%matrix, one, -vec_pot(idir))
END DO
END DO
CALL dbcsr_deallocate_matrix_set(momentum)
! add quadratic term to real part of KS matrix
factor = 0._dp
DO idir = 1, 3
factor = factor + vec_pot(idir)**2
END DO
DO image = 1, nimages
CALL dbcsr_add(matrix_h(1, image)%matrix, matrix_s(1, image)%matrix, one, 0.5*factor)
END DO
IF (ppnl_present) THEN
IF (dft_control%rtp_control%nl_gauge_transform) THEN
NULLIFY (nl_term)
CALL dbcsr_allocate_matrix_set(nl_term, 2)
CALL dbcsr_init_p(nl_term(1)%matrix)
CALL dbcsr_create(nl_term(1)%matrix, template=matrix_s(1, 1)%matrix, &
matrix_type=dbcsr_type_symmetric, name="nl gauge term real part")
CALL cp_dbcsr_alloc_block_from_nbl(nl_term(1)%matrix, sab_orb)
CALL dbcsr_set(nl_term(1)%matrix, zero)
CALL dbcsr_init_p(nl_term(2)%matrix)
CALL dbcsr_create(nl_term(2)%matrix, template=matrix_s(1, 1)%matrix, &
matrix_type=dbcsr_type_antisymmetric, name="nl gauge term imaginary part")
CALL cp_dbcsr_alloc_block_from_nbl(nl_term(2)%matrix, sab_orb)
CALL dbcsr_set(nl_term(2)%matrix, zero)
CALL velocity_gauge_nl_term(qs_env, nl_term, vec_pot)
DO image = 1, nimages
CALL dbcsr_add(matrix_h(1, image)%matrix, nl_term(1)%matrix, one, one)
CALL dbcsr_add(matrix_h_im(1, image)%matrix, nl_term(2)%matrix, one, one)
END DO
CALL dbcsr_deallocate_matrix_set(nl_term)
END IF
END IF
CALL timestop(handle)
END SUBROUTINE velocity_gauge_ks_matrix
! **************************************************************************************************
!> \brief ...
!> \param dft_control ...
!> \param cell ...
!> \param vec_pot ...
! **************************************************************************************************
PURE SUBROUTINE get_vector_potential(dft_control, cell, vec_pot)
TYPE(dft_control_type), POINTER :: dft_control
TYPE(cell_type), POINTER :: cell
REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: vec_pot
CHARACTER(len=*), PARAMETER :: routineN = 'get_vector_potential'
! fieldstrength
! for a delta pulse: vec_pot = - c * kvec
! the speed of light cancels in the terms in the velocity gauge Hamiltonian and is not
! taken into account here
vec_pot(:) = cell%h_inv(1, :)*dft_control%rtp_control%delta_pulse_direction(1) + &
cell%h_inv(2, :)*dft_control%rtp_control%delta_pulse_direction(2) + &
cell%h_inv(3, :)*dft_control%rtp_control%delta_pulse_direction(3)
vec_pot = -vec_pot*twopi*dft_control%rtp_control%delta_pulse_scale
END SUBROUTINE get_vector_potential
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param nl_term ...
!> \param vec_pot ...
! **************************************************************************************************
SUBROUTINE velocity_gauge_nl_term(qs_env, nl_term, vec_pot)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
POINTER :: nl_term
REAL(KIND=dp), DIMENSION(3), INTENT(in) :: vec_pot
CHARACTER(len=*), PARAMETER :: routiuneN = "velocity_gauge_nl_term"
INTEGER :: handle, i, iac, iatom, ibc, icol, ikind, &
irow, jatom, jkind, kac, kbc, kkind, &
maxl, maxlgto, maxlppnl, na, natom, &
nb, nkind, np, slot
INTEGER, DIMENSION(3) :: cell_b
LOGICAL :: found
REAL(dp) :: eps_ppnl
REAL(KIND=dp), DIMENSION(3) :: rab
REAL(KIND=dp), DIMENSION(:, :), POINTER :: imag_block, real_block
REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint_cos, achint_sin, acint_cos, &
acint_sin, bchint_cos, bchint_sin, &
bcint_cos, bcint_sin
TYPE(alist_type), POINTER :: alist_cos_ac, alist_cos_bc, &
alist_sin_ac, alist_sin_bc
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_p_type), ALLOCATABLE, &
DIMENSION(:) :: basis_set
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int_cos, sap_int_sin
!$ INTEGER(kind=omp_lock_kind), &
!$ ALLOCATABLE, DIMENSION(:) :: locks
!$ INTEGER(KIND=int_8) :: iatom8
!$ INTEGER :: lock_num, hash
!$ INTEGER, PARAMETER :: nlock = 501
MARK_USED(int_8)
CALL timeset(routiuneN, handle)
NULLIFY (sap_ppnl, sab_orb)
CALL get_qs_env(qs_env, &
sap_ppnl=sap_ppnl, &
sab_orb=sab_orb)
IF (ASSOCIATED(sap_ppnl)) THEN
NULLIFY (qs_kind_set, particle_set, cell, dft_control)
CALL get_qs_env(qs_env, &
dft_control=dft_control, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, &
cell=cell, &
atomic_kind_set=atomic_kind_set)
nkind = SIZE(atomic_kind_set)
natom = SIZE(particle_set)
eps_ppnl = dft_control%qs_control%eps_ppnl
CALL get_qs_kind_set(qs_kind_set, &
maxlgto=maxlgto, &
maxlppnl=maxlppnl)
maxl = MAX(maxlppnl, maxlgto)
CALL init_orbital_pointers(maxl + 1)
! initalize sab_int types to store the integrals
NULLIFY (sap_int_cos, sap_int_sin)
ALLOCATE (sap_int_cos(nkind*nkind), sap_int_sin(nkind*nkind))
DO i = 1, SIZE(sap_int_cos)
NULLIFY (sap_int_cos(i)%alist, sap_int_cos(i)%asort, sap_int_cos(i)%aindex)
sap_int_cos(i)%nalist = 0
NULLIFY (sap_int_sin(i)%alist, sap_int_sin(i)%asort, sap_int_sin(i)%aindex)
sap_int_sin(i)%nalist = 0
END DO
! get basis set
ALLOCATE (basis_set(nkind))
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
IF (ASSOCIATED(orb_basis_set)) THEN
basis_set(ikind)%gto_basis_set => orb_basis_set
ELSE
NULLIFY (basis_set(ikind)%gto_basis_set)
END IF
END DO
! calculate exponential integrals
CALL build_sap_exp_ints(sap_int_cos, sap_int_sin, sap_ppnl, qs_kind_set, particle_set, &
cell, kvec=vec_pot, basis_set=basis_set, nkind=nkind)
CALL sap_sort(sap_int_cos)
CALL sap_sort(sap_int_sin)
! assemble the integrals for the gauge term
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP SHARED (basis_set, nl_term, sab_orb, sap_int_cos, sap_int_sin, eps_ppnl, locks, nkind, natom) &
!$OMP PRIVATE (real_block, imag_block, acint_cos, achint_cos, bcint_cos, bchint_cos, acint_sin,&
!$OMP achint_sin, bcint_sin, bchint_sin, slot, ikind, jkind, iatom, jatom, cell_b, rab, irow, icol,&
!$OMP found, kkind, iac, ibc, alist_cos_ac, alist_cos_bc, alist_sin_ac, alist_sin_bc, kac, kbc, &
!$OMP na, np, nb, iatom8, hash)
!$OMP SINGLE
!$ ALLOCATE (locks(nlock))
!$OMP END SINGLE
!$OMP DO
!$ DO lock_num = 1, nlock
!$ call omp_init_lock(locks(lock_num))
!$ END DO
!$OMP END DO
NULLIFY (real_block, imag_block)
NULLIFY (acint_cos, achint_cos, achint_cos, bchint_cos)
NULLIFY (acint_sin, achint_sin, achint_sin, bchint_sin)
! loop over atom pairs
!$OMP DO SCHEDULE(GUIDED)
DO slot = 1, sab_orb(1)%nl_size
ikind = sab_orb(1)%nlist_task(slot)%ikind
jkind = sab_orb(1)%nlist_task(slot)%jkind
iatom = sab_orb(1)%nlist_task(slot)%iatom
jatom = sab_orb(1)%nlist_task(slot)%jatom
cell_b(:) = sab_orb(1)%nlist_task(slot)%cell
rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
IF (iatom <= jatom) THEN
irow = iatom
icol = jatom
ELSE
irow = jatom
icol = iatom
END IF
CALL dbcsr_get_block_p(nl_term(1)%matrix, irow, icol, real_block, found)
CALL dbcsr_get_block_p(nl_term(2)%matrix, irow, icol, imag_block, found)
IF (ASSOCIATED(real_block) .AND. ASSOCIATED(imag_block)) THEN
! loop over the <gto_a|ppln_c>h_ij<ppnl_c|gto_b> pairs
DO kkind = 1, nkind
iac = ikind + nkind*(kkind - 1)
ibc = jkind + nkind*(kkind - 1)
IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist)) CYCLE
IF (.NOT. ASSOCIATED(sap_int_cos(ibc)%alist)) CYCLE
IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist)) CYCLE
IF (.NOT. ASSOCIATED(sap_int_sin(ibc)%alist)) CYCLE
CALL get_alist(sap_int_cos(iac), alist_cos_ac, iatom)
CALL get_alist(sap_int_cos(ibc), alist_cos_bc, jatom)
CALL get_alist(sap_int_sin(iac), alist_sin_ac, iatom)
CALL get_alist(sap_int_sin(ibc), alist_sin_bc, jatom)
IF (.NOT. ASSOCIATED(alist_cos_ac)) CYCLE
IF (.NOT. ASSOCIATED(alist_cos_bc)) CYCLE
IF (.NOT. ASSOCIATED(alist_sin_ac)) CYCLE
IF (.NOT. ASSOCIATED(alist_sin_bc)) CYCLE
! only use cos for indexing, as cos and sin integrals are constructed by the same routine
! in the same way
DO kac = 1, alist_cos_ac%nclist
DO kbc = 1, alist_cos_bc%nclist
! the next two ifs should be the same for sine integrals
IF (alist_cos_ac%clist(kac)%catom /= alist_cos_bc%clist(kbc)%catom) CYCLE
IF (ALL(cell_b + alist_cos_bc%clist(kbc)%cell - alist_cos_ac%clist(kac)%cell == 0)) THEN
! screening
IF (alist_cos_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
.AND. alist_cos_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl &
.AND. alist_sin_ac%clist(kac)%maxac*alist_cos_bc%clist(kbc)%maxach < eps_ppnl &
.AND. alist_sin_ac%clist(kac)%maxac*alist_sin_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
acint_cos => alist_cos_ac%clist(kac)%acint
bcint_cos => alist_cos_bc%clist(kbc)%acint
achint_cos => alist_cos_ac%clist(kac)%achint
bchint_cos => alist_cos_bc%clist(kbc)%achint
acint_sin => alist_sin_ac%clist(kac)%acint
bcint_sin => alist_sin_bc%clist(kbc)%acint
achint_sin => alist_sin_ac%clist(kac)%achint
bchint_sin => alist_sin_bc%clist(kbc)%achint
na = SIZE(acint_cos, 1)
np = SIZE(acint_cos, 2)
nb = SIZE(bcint_cos, 1)
!$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
!$ hash = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
!$ CALL omp_set_lock(locks(hash))
IF (iatom <= jatom) THEN
! cos*cos + sin*sin
real_block(1:na, 1:nb) = real_block(1:na, 1:nb) + &
MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1))) + &
MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1)))
! sin * cos - cos * sin
imag_block(1:na, 1:nb) = imag_block(1:na, 1:nb) - &
MATMUL(achint_sin(1:na, 1:np, 1), TRANSPOSE(bcint_cos(1:nb, 1:np, 1))) + &
MATMUL(achint_cos(1:na, 1:np, 1), TRANSPOSE(bcint_sin(1:nb, 1:np, 1)))
ELSE
! cos*cos + sin*sin
real_block(1:nb, 1:na) = real_block(1:nb, 1:na) + &
MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1))) + &
MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1)))
! sin * cos - cos * sin
imag_block(1:nb, 1:na) = imag_block(1:nb, 1:na) - &
MATMUL(bchint_sin(1:nb, 1:np, 1), TRANSPOSE(acint_cos(1:na, 1:np, 1))) + &
MATMUL(bchint_cos(1:nb, 1:np, 1), TRANSPOSE(acint_sin(1:na, 1:np, 1)))
END IF
!$ CALL omp_unset_lock(locks(hash))
EXIT
END IF
END DO
END DO
END DO
END IF
END DO
!$OMP DO
!$ DO lock_num = 1, nlock
!$ call omp_destroy_lock(locks(lock_num))
!$ END DO
!$OMP END DO
!$OMP SINGLE
!$ DEALLOCATE (locks)
!$OMP END SINGLE NOWAIT
!$OMP END PARALLEL
CALL release_sap_int(sap_int_cos)
CALL release_sap_int(sap_int_sin)
DEALLOCATE (basis_set)
END IF
CALL timestop(handle)
END SUBROUTINE velocity_gauge_nl_term
! **************************************************************************************************
!> \brief calculate <a|sin/cos|p> integrals and store in sap_int_type
!> adapted from build_sap_ints
!> Do this on each MPI task as the integrals need to be available globally.
!> Might be faster than communicating as the integrals are obtained analytically.
!> \param sap_int_cos ...
!> \param sap_int_sin ...
!> \param sap_ppnl ...
!> \param qs_kind_set ...
!> \param particle_set ...
!> \param cell ...
!> \param kvec ...
!> \param basis_set ...
!> \param nkind ...
! **************************************************************************************************
SUBROUTINE build_sap_exp_ints(sap_int_cos, sap_int_sin, sap_ppnl, qs_kind_set, particle_set, cell, kvec, basis_set, nkind)
TYPE(sap_int_type), DIMENSION(:), INTENT(INOUT), &
POINTER :: sap_int_cos, sap_int_sin
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
INTENT(IN), POINTER :: sap_ppnl
TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
POINTER :: qs_kind_set
TYPE(particle_type), DIMENSION(:), INTENT(IN), &
POINTER :: particle_set
TYPE(cell_type), INTENT(IN), POINTER :: cell
REAL(KIND=dp), DIMENSION(3), INTENT(in) :: kvec
TYPE(gto_basis_set_p_type), DIMENSION(:), &
INTENT(IN) :: basis_set
INTEGER, INTENT(IN) :: nkind
CHARACTER(len=*), PARAMETER :: routiuneN = "build_sap_exp_ints"
INTEGER :: handle, i, iac, iatom, ikind, ilist, iset, jneighbor, katom, kkind, l, lc_max, &
lc_min, ldai, ldints, lppnl, maxco, maxl, maxlgto, maxlppnl, maxppnl, maxsgf, na, nb, &
ncoa, ncoc, nlist, nneighbor, np, nppnl, nprjc, nseta, nsgfa, prjc, sgfa, slot
INTEGER, DIMENSION(3) :: cell_c
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
nsgf_seta
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
LOGICAL :: dogth
REAL(KIND=dp) :: dac, ppnl_radius
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ai_work_cos, ai_work_sin, work_cos, &
work_sin
REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc
REAL(KIND=dp), DIMENSION(3) :: ra, rac, raf, rc, rcf
REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_ppnl, set_radius_a
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, rpgfa, sphi_a, vprj_ppnl, zeta
TYPE(clist_type), POINTER :: clist, clist_sin
TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
TYPE(sgp_potential_type), POINTER :: sgp_potential
CALL timeset(routiuneN, handle)
CALL get_qs_kind_set(qs_kind_set, &
maxco=maxco, &
maxlppnl=maxlppnl, &
maxppnl=maxppnl, &
maxsgf=maxsgf, &
maxlgto=maxlgto)
! maximum dimensions for allocations
maxl = MAX(maxlppnl, maxlgto)
ldints = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
ldai = ncoset(maxl + 1)
!set up direct access to basis and potential
NULLIFY (gpotential, spotential)
ALLOCATE (gpotential(nkind), spotential(nkind))
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
NULLIFY (gpotential(ikind)%gth_potential)
NULLIFY (spotential(ikind)%sgp_potential)
IF (ASSOCIATED(gth_potential)) THEN
gpotential(ikind)%gth_potential => gth_potential
ELSE IF (ASSOCIATED(sgp_potential)) THEN
spotential(ikind)%sgp_potential => sgp_potential
END IF
END DO
!allocate sap int
NULLIFY (clist)
DO slot = 1, sap_ppnl(1)%nl_size
ikind = sap_ppnl(1)%nlist_task(slot)%ikind
kkind = sap_ppnl(1)%nlist_task(slot)%jkind
iatom = sap_ppnl(1)%nlist_task(slot)%iatom
katom = sap_ppnl(1)%nlist_task(slot)%jatom
nlist = sap_ppnl(1)%nlist_task(slot)%nlist
ilist = sap_ppnl(1)%nlist_task(slot)%ilist
nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
iac = ikind + nkind*(kkind - 1)
IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
IF (.NOT. ASSOCIATED(gpotential(kkind)%gth_potential) .AND. &
.NOT. ASSOCIATED(spotential(kkind)%sgp_potential)) CYCLE
IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist)) THEN
sap_int_cos(iac)%a_kind = ikind
sap_int_cos(iac)%p_kind = kkind
sap_int_cos(iac)%nalist = nlist
ALLOCATE (sap_int_cos(iac)%alist(nlist))
DO i = 1, nlist
NULLIFY (sap_int_cos(iac)%alist(i)%clist)
sap_int_cos(iac)%alist(i)%aatom = 0
sap_int_cos(iac)%alist(i)%nclist = 0
END DO
END IF
IF (.NOT. ASSOCIATED(sap_int_cos(iac)%alist(ilist)%clist)) THEN
sap_int_cos(iac)%alist(ilist)%aatom = iatom
sap_int_cos(iac)%alist(ilist)%nclist = nneighbor
ALLOCATE (sap_int_cos(iac)%alist(ilist)%clist(nneighbor))
DO i = 1, nneighbor
clist => sap_int_cos(iac)%alist(ilist)%clist(i)
clist%catom = 0
NULLIFY (clist%acint)
NULLIFY (clist%achint)
NULLIFY (clist%sgf_list)
END DO
END IF
IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist)) THEN
sap_int_sin(iac)%a_kind = ikind
sap_int_sin(iac)%p_kind = kkind
sap_int_sin(iac)%nalist = nlist
ALLOCATE (sap_int_sin(iac)%alist(nlist))
DO i = 1, nlist
NULLIFY (sap_int_sin(iac)%alist(i)%clist)
sap_int_sin(iac)%alist(i)%aatom = 0
sap_int_sin(iac)%alist(i)%nclist = 0
END DO
END IF
IF (.NOT. ASSOCIATED(sap_int_sin(iac)%alist(ilist)%clist)) THEN
sap_int_sin(iac)%alist(ilist)%aatom = iatom
sap_int_sin(iac)%alist(ilist)%nclist = nneighbor
ALLOCATE (sap_int_sin(iac)%alist(ilist)%clist(nneighbor))
DO i = 1, nneighbor
clist => sap_int_sin(iac)%alist(ilist)%clist(i)
clist%catom = 0
NULLIFY (clist%acint)
NULLIFY (clist%achint)
NULLIFY (clist%sgf_list)
END DO
END IF
END DO
! actual calculation of the integrals <a|cos|p> and <a|sin|p>
! allocate temporary storage using maximum dimensions
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP SHARED (basis_set, gpotential, ncoset, sap_ppnl, sap_int_cos, sap_int_sin, nkind, &
!$OMP ldints, maxco, nco, cell, particle_set, kvec) &
!$OMP PRIVATE (slot, ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
!$OMP cell_c, rac, dac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta,&
!$OMP rpgfa, set_radius_a, sphi_a, zeta, alpha_ppnl, cprj, lppnl, nppnl, nprj_ppnl,&
!$OMP ppnl_radius, vprj_ppnl, clist, clist_sin, ra, rc, ncoa, sgfa, prjc, work_cos, work_sin,&
!$OMP nprjc, rprjc, lc_max, lc_min, zetc, ncoc, ai_work_sin, ai_work_cos, na, nb, np, dogth, raf, rcf)
ALLOCATE (work_cos(ldints, ldints), work_sin(ldints, ldints))
ALLOCATE (ai_work_cos(maxco, maxco), ai_work_sin(maxco, maxco))
work_cos = 0.0_dp
work_sin = 0.0_dp
ai_work_cos = 0.0_dp
ai_work_sin = 0.0_dp
dogth = .FALSE.
NULLIFY (first_sgfa, la_max, la_min, npgfa, nsgf_seta, rpgfa, set_radius_a, sphi_a, zeta)
NULLIFY (alpha_ppnl, cprj, nprj_ppnl, vprj_ppnl)
NULLIFY (clist, clist_sin)
!$OMP DO SCHEDULE(GUIDED)
DO slot = 1, sap_ppnl(1)%nl_size
ikind = sap_ppnl(1)%nlist_task(slot)%ikind
kkind = sap_ppnl(1)%nlist_task(slot)%jkind
iatom = sap_ppnl(1)%nlist_task(slot)%iatom
katom = sap_ppnl(1)%nlist_task(slot)%jatom
nlist = sap_ppnl(1)%nlist_task(slot)%nlist
ilist = sap_ppnl(1)%nlist_task(slot)%ilist
nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
jneighbor = sap_ppnl(1)%nlist_task(slot)%inode
cell_c(:) = sap_ppnl(1)%nlist_task(slot)%cell(:)
rac(1:3) = sap_ppnl(1)%nlist_task(slot)%r(1:3)
dac = NORM2(rac)
iac = ikind + nkind*(kkind - 1)
IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
! get definition of gto basis set
first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf
la_max => basis_set(ikind)%gto_basis_set%lmax
la_min => basis_set(ikind)%gto_basis_set%lmin
npgfa => basis_set(ikind)%gto_basis_set%npgf
nseta = basis_set(ikind)%gto_basis_set%nset
nsgfa = basis_set(ikind)%gto_basis_set%nsgf
nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set
rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius
set_radius_a => basis_set(ikind)%gto_basis_set%set_radius
sphi_a => basis_set(ikind)%gto_basis_set%sphi
zeta => basis_set(ikind)%gto_basis_set%zet
IF (ASSOCIATED(gpotential(kkind)%gth_potential)) THEN
! GTH potential
dogth = .TRUE.
alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl
cprj => gpotential(kkind)%gth_potential%cprj
lppnl = gpotential(kkind)%gth_potential%lppnl
nppnl = gpotential(kkind)%gth_potential%nppnl
nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
ELSE
CYCLE
END IF
clist => sap_int_cos(iac)%alist(ilist)%clist(jneighbor)
clist_sin => sap_int_sin(iac)%alist(ilist)%clist(jneighbor)
clist%catom = katom
clist%cell = cell_c
clist%rac = rac
clist_sin%catom = katom
clist_sin%cell = cell_c
clist_sin%rac = rac
ALLOCATE (clist%acint(nsgfa, nppnl, 1), clist%achint(nsgfa, nppnl, 1))
clist%acint = 0.0_dp
clist%achint = 0.0_dp
clist%nsgf_cnt = 0
ALLOCATE (clist_sin%acint(nsgfa, nppnl, 1), clist_sin%achint(nsgfa, nppnl, 1))
clist_sin%acint = 0.0_dp
clist_sin%achint = 0.0_dp
clist_sin%nsgf_cnt = 0
! reference point at zero
ra(:) = pbc(particle_set(iatom)%r(:), cell)
rc(:) = ra + rac
! reference point at pseudized atom
raf(:) = ra - rc
rcf(:) = 0._dp
DO iset = 1, nseta
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
IF (dogth) THEN
prjc = 1
work_cos = 0.0_dp
work_sin = 0.0_dp
DO l = 0, lppnl
nprjc = nprj_ppnl(l)*nco(l)
IF (nprjc == 0) CYCLE
rprjc(1) = ppnl_radius
IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE
lc_max = l + 2*(nprj_ppnl(l) - 1)
lc_min = l
zetc(1) = alpha_ppnl(l)
ncoc = ncoset(lc_max)
CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
lc_max, 1, zetc, rprjc, lc_min, raf, rcf, kvec, ai_work_cos, ai_work_sin)
! projector functions: Cartesian -> spherical
na = ncoa
nb = nprjc
np = ncoc
work_cos(1:na, prjc:prjc + nb - 1) = &
MATMUL(ai_work_cos(1:na, 1:np), cprj(1:np, prjc:prjc + nb - 1))
work_sin(1:na, prjc:prjc + nb - 1) = &
MATMUL(ai_work_sin(1:na, 1:np), cprj(1:np, prjc:prjc + nb - 1))
prjc = prjc + nprjc
END DO
! contract gto basis set into acint
na = nsgf_seta(iset)
nb = nppnl
np = ncoa
clist%acint(sgfa:sgfa + na - 1, 1:nb, 1) = &
MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_cos(1:np, 1:nb))
clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1) = &
MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_sin(1:np, 1:nb))
! multiply with interaction matrix h_ij of the nl pp
clist%achint(sgfa:sgfa + na - 1, 1:nb, 1) = &
MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, 1), vprj_ppnl(1:nb, 1:nb))
clist_sin%achint(sgfa:sgfa + na - 1, 1:nb, 1) = &
MATMUL(clist_sin%acint(sgfa:sgfa + na - 1, 1:nb, 1), vprj_ppnl(1:nb, 1:nb))
END IF
END DO
clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
clist_sin%maxac = MAXVAL(ABS(clist_sin%acint(:, :, 1)))
clist_sin%maxach = MAXVAL(ABS(clist_sin%achint(:, :, 1)))
END DO
DEALLOCATE (work_cos, work_sin, ai_work_cos, ai_work_sin)
!$OMP END PARALLEL
DEALLOCATE (gpotential, spotential)
CALL timestop(handle)
END SUBROUTINE build_sap_exp_ints
END MODULE rt_propagation_velocity_gauge

23
src/sap_kind_types.F
View file

@ -282,9 +282,9 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_sap_ints'
INTEGER :: first_col, handle, i, iac, iatom, ikind, ilist, iset, j, jneighbor, katom, kkind, &
l, lc_max, lc_min, ldai, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, maxppnl, &
maxsgf, na, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, nseta, nsgfa, &
prjc, sgfa, slot
l, lc_max, lc_min, ldai, ldai_nl, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, &
maxppnl, maxsgf, na, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, &
nseta, nsgfa, prjc, sgfa, slot
INTEGER, DIMENSION(3) :: cell_c
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
nsgf_seta
@ -339,9 +339,14 @@ CONTAINS
CALL init_orbital_pointers(maxl + nder + 1)
ldsab = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
ldai = ncoset(maxl + nder + 1)
IF (.NOT. my_momentmode) THEN
ldai = ncoset(maxl + nder + 1)
ELSE
ldai = maxco
END IF
!set up direct access to basis and potential
ldai_nl = 0
NULLIFY (gpotential, spotential)
ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind))
DO ikind = 1, nkind
@ -356,8 +361,10 @@ CONTAINS
NULLIFY (spotential(ikind)%sgp_potential)
IF (ASSOCIATED(gth_potential)) THEN
gpotential(ikind)%gth_potential => gth_potential
ldai_nl = MAX(ldai_nl, ncoset(gth_potential%lprj_ppnl_max))
ELSE IF (ASSOCIATED(sgp_potential)) THEN
spotential(ikind)%sgp_potential => sgp_potential
ldai_nl = MAX(ldai_nl, sgp_potential%n_nonlocal*ncoset(sgp_potential%lmax))
END IF
END DO
@ -400,7 +407,7 @@ CONTAINS
!calculate the overlap integrals <a|pp>
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP SHARED (basis_set, gpotential, spotential, maxder, ncoset, my_momentmode, &
!$OMP SHARED (basis_set, gpotential, spotential, maxder, ncoset, my_momentmode, ldai_nl, &
!$OMP sap_ppnl, sap_int, nkind, ldsab, ldai, nder, nco, my_ref, cell, particle_set, rf) &
!$OMP PRIVATE (ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
!$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, &
@ -413,7 +420,11 @@ CONTAINS
! allocate work storage
ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder))
sab = 0.0_dp
ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1)))
IF (.NOT. my_momentmode) THEN
ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1)))
ELSE
ALLOCATE (ai_work(ldai, ldai_nl, ncoset(nder + 1)))
END IF
ai_work = 0.0_dp
!$OMP DO SCHEDULE(GUIDED)

1
tests/QS/regtest-rtp-3/H2O-delta-01.inp
View file

@ -29,6 +29,7 @@
DELTA_PULSE_DIRECTION 1 0 0
DELTA_PULSE_SCALE 0.001
PERIODIC .TRUE.
COM_NL .FALSE.
MAX_ITER 50
MAT_EXP PADE
EXP_ACCURACY 1.0E-10

1
tests/QS/regtest-rtp-3/H2O-delta-02.inp
View file

@ -27,6 +27,7 @@
DELTA_PULSE_DIRECTION 1 0 0
DELTA_PULSE_SCALE 0.001
PERIODIC .TRUE.
COM_NL .FALSE.
MAX_ITER 50
MAT_EXP PADE
EXP_ACCURACY 1.0E-10

82
tests/QS/regtest-rtp-4/H2O-delta-lenrep.inp Normal file
View file

@ -0,0 +1,82 @@
&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
LEN_REP 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
PERIODIC .FALSE.
! FILENAME dipole
&END
&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

82
tests/QS/regtest-rtp-4/H2O-delta-mag.inp Normal file
View file

@ -0,0 +1,82 @@
&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_MAG T
GAUGE_ORIG COM
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
PERIODIC .FALSE.
! FILENAME dipole
&END
&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

82
tests/QS/regtest-rtp-4/H2O-vg.inp Normal file
View file

@ -0,0 +1,82 @@
&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 F
VELOCITY_GAUGE T
VG_COM_NL T
DELTA_PULSE_DIRECTION 1 0 0
DELTA_PULSE_SCALE 0.001
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
PERIODIC .FALSE.
! FILENAME dipole
&END
&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

11
tests/QS/regtest-rtp-4/TEST_FILES Normal file
View file

@ -0,0 +1,11 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
#
# Delta pulse in the length representation, velocity gauge and magnetic delta pulse
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.17816524491995
#EOF

2
tests/TEST_DIRS
View file

@ -327,4 +327,4 @@ Fist/regtest-spgr spglib
Fist/regtest-plumed2 plumed2
Fist/regtest-quip quip
Fist/regtest-16
QS/regtest-rtp-4