From e64daf41ce02fb1dd7f72831acc068f4c092e916 Mon Sep 17 00:00:00 2001 From: mattiatj <44839849+mattiatj@users.noreply.github.com> Date: Tue, 11 Oct 2022 16:05:52 +0200 Subject: [PATCH] RTP gauges (#2343) * 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 --- src/common/bibliography.F | 34 +- src/commutator_rpnl.F | 261 +++++- src/cp_control_types.F | 3 + src/cp_control_utils.F | 6 + src/emd/rt_delta_pulse.F | 610 +++++++++----- src/emd/rt_hfx_utils.F | 141 +++- src/emd/rt_propagation_methods.F | 11 +- src/hfx_admm_utils.F | 15 +- src/input_cp2k_dft.F | 61 +- src/motion/rt_propagation.F | 48 +- src/qs_core_hamiltonian.F | 22 +- src/qs_environment_types.F | 24 +- src/qs_ks_methods.F | 23 +- src/qs_ks_types.F | 71 +- src/qs_moments.F | 182 +++-- src/qs_rho_types.F | 25 +- src/rt_propagation_velocity_gauge.F | 832 ++++++++++++++++++++ src/sap_kind_types.F | 23 +- tests/QS/regtest-rtp-3/H2O-delta-01.inp | 1 + tests/QS/regtest-rtp-3/H2O-delta-02.inp | 1 + tests/QS/regtest-rtp-4/H2O-delta-lenrep.inp | 82 ++ tests/QS/regtest-rtp-4/H2O-delta-mag.inp | 82 ++ tests/QS/regtest-rtp-4/H2O-vg.inp | 82 ++ tests/QS/regtest-rtp-4/TEST_FILES | 11 + tests/TEST_DIRS | 2 +- 25 files changed, 2269 insertions(+), 384 deletions(-) create mode 100644 src/rt_propagation_velocity_gauge.F create mode 100644 tests/QS/regtest-rtp-4/H2O-delta-lenrep.inp create mode 100644 tests/QS/regtest-rtp-4/H2O-delta-mag.inp create mode 100644 tests/QS/regtest-rtp-4/H2O-vg.inp create mode 100644 tests/QS/regtest-rtp-4/TEST_FILES diff --git a/src/common/bibliography.F b/src/common/bibliography.F index 23b859fcfd..4692840584 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -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 diff --git a/src/commutator_rpnl.F b/src/commutator_rpnl.F index 356be95a97..ef12a1973f 100644 --- a/src/commutator_rpnl.F +++ b/src/commutator_rpnl.F @@ -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 diff --git a/src/cp_control_types.F b/src/cp_control_types.F index b9a802906b..23260d5bc3 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -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 diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index 42f5bce2a7..006c549c43 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -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", & diff --git a/src/emd/rt_delta_pulse.F b/src/emd/rt_delta_pulse.F index 6fd3c0fbab..6a58e51e65 100644 --- a/src/emd/rt_delta_pulse.F +++ b/src/emd/rt_delta_pulse.F @@ -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 diff --git a/src/emd/rt_hfx_utils.F b/src/emd/rt_hfx_utils.F index dc0860dfb4..2d35a3fa8f 100644 --- a/src/emd/rt_hfx_utils.F +++ b/src/emd/rt_hfx_utils.F @@ -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 diff --git a/src/emd/rt_propagation_methods.F b/src/emd/rt_propagation_methods.F index b26810aa29..76b8e153ef 100644 --- a/src/emd/rt_propagation_methods.F +++ b/src/emd/rt_propagation_methods.F @@ -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 diff --git a/src/hfx_admm_utils.F b/src/hfx_admm_utils.F index 37772d3a00..84c6428386 100644 --- a/src/hfx_admm_utils.F +++ b/src/hfx_admm_utils.F @@ -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 diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 39508e0af7..989c5b80f7 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -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) diff --git a/src/motion/rt_propagation.F b/src/motion/rt_propagation.F index 1ba0b07f68..126078f3e7 100644 --- a/src/motion/rt_propagation.F +++ b/src/motion/rt_propagation.F @@ -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) diff --git a/src/qs_core_hamiltonian.F b/src/qs_core_hamiltonian.F index 4d0036670f..ac18e709b7 100644 --- a/src/qs_core_hamiltonian.F +++ b/src/qs_core_hamiltonian.F @@ -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 diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F index c5702cd1d8..bacf627ee4 100644 --- a/src/qs_environment_types.F +++ b/src/qs_environment_types.F @@ -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, & diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F index 0c618bd606..e6dcda79d3 100644 --- a/src/qs_ks_methods.F +++ b/src/qs_ks_methods.F @@ -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 diff --git a/src/qs_ks_types.F b/src/qs_ks_types.F index 48dc9550e9..ba8e777cf1 100644 --- a/src/qs_ks_types.F +++ b/src/qs_ks_types.F @@ -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) diff --git a/src/qs_moments.F b/src/qs_moments.F index 5b30ab5b90..d42a1755ad 100644 --- a/src/qs_moments.F +++ b/src/qs_moments.F @@ -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 diff --git a/src/qs_rho_types.F b/src/qs_rho_types.F index 5d42d5d7f1..73bf5e4ef2 100644 --- a/src/qs_rho_types.F +++ b/src/qs_rho_types.F @@ -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 diff --git a/src/rt_propagation_velocity_gauge.F b/src/rt_propagation_velocity_gauge.F new file mode 100644 index 0000000000..efa42397d7 --- /dev/null +++ b/src/rt_propagation_velocity_gauge.F @@ -0,0 +1,832 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2022 CP2K developers group ! +! ! +! 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 h_ij 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 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 and + ! 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 diff --git a/src/sap_kind_types.F b/src/sap_kind_types.F index 7ebf27d9bf..b4c40287d5 100644 --- a/src/sap_kind_types.F +++ b/src/sap_kind_types.F @@ -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 !$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) diff --git a/tests/QS/regtest-rtp-3/H2O-delta-01.inp b/tests/QS/regtest-rtp-3/H2O-delta-01.inp index 63b4625fdf..e8acd48e77 100644 --- a/tests/QS/regtest-rtp-3/H2O-delta-01.inp +++ b/tests/QS/regtest-rtp-3/H2O-delta-01.inp @@ -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 diff --git a/tests/QS/regtest-rtp-3/H2O-delta-02.inp b/tests/QS/regtest-rtp-3/H2O-delta-02.inp index 63aa5b51b2..88cc44774f 100644 --- a/tests/QS/regtest-rtp-3/H2O-delta-02.inp +++ b/tests/QS/regtest-rtp-3/H2O-delta-02.inp @@ -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 diff --git a/tests/QS/regtest-rtp-4/H2O-delta-lenrep.inp b/tests/QS/regtest-rtp-4/H2O-delta-lenrep.inp new file mode 100644 index 0000000000..5a5d80ec43 --- /dev/null +++ b/tests/QS/regtest-rtp-4/H2O-delta-lenrep.inp @@ -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 + diff --git a/tests/QS/regtest-rtp-4/H2O-delta-mag.inp b/tests/QS/regtest-rtp-4/H2O-delta-mag.inp new file mode 100644 index 0000000000..120e139ecc --- /dev/null +++ b/tests/QS/regtest-rtp-4/H2O-delta-mag.inp @@ -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 + diff --git a/tests/QS/regtest-rtp-4/H2O-vg.inp b/tests/QS/regtest-rtp-4/H2O-vg.inp new file mode 100644 index 0000000000..c56f6539ef --- /dev/null +++ b/tests/QS/regtest-rtp-4/H2O-vg.inp @@ -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 + diff --git a/tests/QS/regtest-rtp-4/TEST_FILES b/tests/QS/regtest-rtp-4/TEST_FILES new file mode 100644 index 0000000000..762a477a68 --- /dev/null +++ b/tests/QS/regtest-rtp-4/TEST_FILES @@ -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 \ No newline at end of file diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 0bc5d52a2d..c2cf011ecf 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -327,4 +327,4 @@ Fist/regtest-spgr spglib Fist/regtest-plumed2 plumed2 Fist/regtest-quip quip Fist/regtest-16 - +QS/regtest-rtp-4