From 78ff48db3aa7bc62878951469cf43e962e79d2f7 Mon Sep 17 00:00:00 2001 From: glb96 <118453102+glb96@users.noreply.github.com> Date: Fri, 24 Feb 2023 13:16:34 +0100 Subject: [PATCH] RTP: Add velocity gauge --- src/commutator_rpnl.F | 181 ++++++++- src/cp_control_types.F | 2 + src/cp_control_utils.F | 8 +- src/efield_utils.F | 101 +++-- src/emd/rt_delta_pulse.F | 8 + src/emd/rt_propagation_methods.F | 26 +- src/emd/rt_propagation_output.F | 66 ++- src/input_cp2k_dft.F | 10 + src/motion/rt_propagation.F | 1 - src/qs_ks_methods.F | 34 +- src/qs_ks_utils.F | 16 +- src/qs_mo_io.F | 1 + src/qs_moments.F | 382 +++++++++++------- src/rt_propagation_velocity_gauge.F | 37 +- .../QS/regtest-rtp-2/H2-emd-efield-custom.inp | 2 +- tests/QS/regtest-rtp-2/H2-emd-efield-ramp.inp | 2 +- tests/QS/regtest-rtp-2/H2-emd-efield.inp | 2 +- tests/QS/regtest-rtp-2/H2-rtp-efield.inp | 2 +- tests/QS/regtest-rtp-2/TEST_FILES | 10 +- tests/QS/regtest-rtp-4/H2O-vg.inp | 14 +- tests/QS/regtest-rtp-4/TEST_FILES | 4 +- 21 files changed, 616 insertions(+), 293 deletions(-) diff --git a/src/commutator_rpnl.F b/src/commutator_rpnl.F index adf5df37e7..4d209c2cd2 100644 --- a/src/commutator_rpnl.F +++ b/src/commutator_rpnl.F @@ -413,8 +413,12 @@ CONTAINS END SUBROUTINE build_com_rpnl ! ************************************************************************************************** -!> \brief Calculate [r,Vnl], r x [r,Vnl] or [rr,Vnl] in AO basis -!> reference point is required for the two latter options +!> \brief Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv) +!> or [rr,Vnl] (matrix_rrv) in AO basis. +!> Reference point is required for the two latter options +!> Update: Calculate rxVnlxr (matrix_rvr) and rxrxVnl + Vnlxrxr (matrix_rrv_vrr) +!> in AO basis. Added in the first place for current correction in +!> the VG formalism (first order wrt vector potential). !> \param qs_kind_set ... !> \param sab_all ... !> \param sap_ppnl ... @@ -424,6 +428,8 @@ CONTAINS !> \param matrix_rv ... !> \param matrix_rxrv ... !> \param matrix_rrv ... +!> \param matrix_rvr ... +!> \param matrix_rrv_vrr ... !> \param matrix_r_rxvr ... !> \param matrix_rxvr_r ... !> \param matrix_r_doublecom ... @@ -431,7 +437,7 @@ CONTAINS !> \param ref_point ... ! ************************************************************************************************** SUBROUTINE build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, & - matrix_rrv, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point) + matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point) TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), & POINTER :: qs_kind_set @@ -442,7 +448,8 @@ CONTAINS POINTER :: particle_set TYPE(cell_type), INTENT(IN), POINTER :: cell TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), & - OPTIONAL :: matrix_rv, matrix_rxrv, matrix_rrv + OPTIONAL :: matrix_rv, matrix_rxrv, matrix_rrv, & + matrix_rvr, matrix_rrv_vrr TYPE(dbcsr_p_type), DIMENSION(:, :), & INTENT(INOUT), OPTIONAL :: matrix_r_rxvr, matrix_rxvr_r, & matrix_r_doublecom @@ -458,12 +465,13 @@ CONTAINS kac, kbc, kkind, na, natom, nb, nkind, & np, order, slot INTEGER, DIMENSION(3) :: cell_b - LOGICAL :: asso_r_doublecom, asso_r_rxvr, asso_rrv, asso_rv, asso_rxrv, asso_rxvr_r, & - do_symmetric, found, go, my_r_doublecom, my_r_rxvr, my_ref, my_rrv, my_rv, my_rxrv, & - my_rxvr_r, periodic, ppnl_present + LOGICAL :: asso_r_doublecom, asso_r_rxvr, asso_rrv, asso_rrv_vrr, asso_rv, asso_rvr, & + asso_rxrv, asso_rxvr_r, do_symmetric, found, go, my_r_doublecom, my_r_rxvr, my_ref, & + my_rrv, my_rrv_vrr, my_rv, my_rvr, my_rxrv, my_rxvr_r, periodic, ppnl_present REAL(KIND=dp), DIMENSION(3) :: rab, rf TYPE(alist_type), POINTER :: alist_ac, alist_bc - TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_rrv, blocks_rv, blocks_rxrv + TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_rrv, blocks_rrv_vrr, blocks_rv, & + blocks_rvr, blocks_rxrv TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_r_doublecom, blocks_r_rxvr, & blocks_rxvr_r TYPE(gto_basis_set_p_type), ALLOCATABLE, & @@ -487,13 +495,17 @@ CONTAINS my_rxrv = .FALSE. my_rrv = .FALSE. my_rv = .FALSE. + my_rvr = .FALSE. + my_rrv_vrr = .FALSE. IF (PRESENT(matrix_r_doublecom)) my_r_doublecom = .TRUE. IF (PRESENT(matrix_r_rxvr)) my_r_rxvr = .TRUE. IF (PRESENT(matrix_rxvr_r)) my_rxvr_r = .TRUE. IF (PRESENT(matrix_rxrv)) my_rxrv = .TRUE. IF (PRESENT(matrix_rrv)) my_rrv = .TRUE. IF (PRESENT(matrix_rv)) my_rv = .TRUE. - IF (.NOT. (my_rv .OR. my_rxrv .OR. my_rrv .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom)) THEN + IF (PRESENT(matrix_rvr)) my_rvr = .TRUE. + IF (PRESENT(matrix_rrv_vrr)) my_rrv_vrr = .TRUE. + IF (.NOT. (my_rv .OR. my_rxrv .OR. my_rrv .OR. my_rvr .OR. my_rrv_vrr .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom)) THEN CPABORT('No dbcsr matrix provided for commutator calculation!') END IF @@ -502,6 +514,8 @@ CONTAINS IF (my_rxrv .OR. my_rrv .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom) THEN order = 2 CPASSERT(PRESENT(ref_point)) ! need reference point for r x [r,Vnl] and [rr,Vnl] + ELSE IF (my_rvr .OR. my_rrv_vrr) THEN + order = 2 ELSE order = 1 END IF @@ -561,15 +575,21 @@ CONTAINS !$OMP PARALLEL & !$OMP DEFAULT (NONE) & -!$OMP SHARED (basis_set, matrix_rv, matrix_rxrv, matrix_rrv, matrix_r_doublecom, & -!$OMP sap_int, natom, nkind, eps_ppnl, my_rv, my_rxrv, my_rrv, my_r_doublecom, locks, sab_all, & +!$OMP SHARED (basis_set, matrix_rv, matrix_rxrv, matrix_rrv, & +!$OMP matrix_rvr, matrix_rrv_vrr, matrix_r_doublecom, & +!$OMP sap_int, natom, nkind, eps_ppnl, locks, sab_all, & +!$OMP my_rv, my_rxrv, my_rrv, my_rvr, my_rrv_vrr, & +!$OMP my_r_doublecom, & !$OMP matrix_r_rxvr, matrix_rxvr_r, my_r_rxvr, my_rxvr_r, & !$OMP pseudoatom, do_symmetric) & !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, & -!$OMP iab, irow, icol, blocks_rv, blocks_rxrv, blocks_rrv, blocks_r_rxvr, blocks_rxvr_r, blocks_r_doublecom, & +!$OMP iab, irow, icol, & +!$OMP blocks_rv, blocks_rxrv, blocks_rrv, blocks_rvr, blocks_rrv_vrr, & +!$OMP blocks_r_rxvr, blocks_rxvr_r, blocks_r_doublecom, & !$OMP found, iac, ibc, alist_ac, alist_bc, & !$OMP na, np, nb, kkind, kac, kbc, i, & -!$OMP go, asso_rv, asso_rxrv, asso_rrv, asso_r_rxvr, asso_rxvr_r, asso_r_doublecom, hash) +!$OMP go, asso_rv, asso_rxrv, asso_rrv, asso_rvr, asso_rrv_vrr, & +!$OMP asso_r_rxvr, asso_rxvr_r, asso_r_doublecom, hash) !$OMP SINGLE !$ ALLOCATE (locks(nlock)) @@ -619,7 +639,12 @@ CONTAINS IF (my_rrv) THEN ALLOCATE (blocks_rrv(6)) END IF - + IF (my_rvr) THEN + ALLOCATE (blocks_rvr(6)) + END IF + IF (my_rrv_vrr) THEN + ALLOCATE (blocks_rrv_vrr(6)) + END IF IF (my_r_rxvr) THEN ALLOCATE (blocks_r_rxvr(3, 3)) END IF @@ -652,6 +677,18 @@ CONTAINS END DO END IF + IF (my_rvr) THEN + DO ind = 1, 6 + CALL dbcsr_get_block_p(matrix_rvr(ind)%matrix, irow, icol, blocks_rvr(ind)%block, found) + END DO + END IF + + IF (my_rrv_vrr) THEN + DO ind = 1, 6 + CALL dbcsr_get_block_p(matrix_rrv_vrr(ind)%matrix, irow, icol, blocks_rrv_vrr(ind)%block, found) + END DO + END IF + IF (my_r_rxvr) THEN DO ind = 1, 3 DO ind2 = 1, 3 @@ -703,6 +740,20 @@ CONTAINS go = go .AND. asso_rrv END IF + IF (my_rvr) THEN + asso_rvr = (ASSOCIATED(blocks_rvr(1)%block) .AND. ASSOCIATED(blocks_rvr(2)%block) .AND. & + ASSOCIATED(blocks_rvr(3)%block) .AND. ASSOCIATED(blocks_rvr(4)%block) .AND. & + ASSOCIATED(blocks_rvr(5)%block) .AND. ASSOCIATED(blocks_rvr(6)%block)) + go = go .AND. asso_rvr + END IF + + IF (my_rrv_vrr) THEN + asso_rrv_vrr = (ASSOCIATED(blocks_rrv_vrr(1)%block) .AND. ASSOCIATED(blocks_rrv_vrr(2)%block) .AND. & + ASSOCIATED(blocks_rrv_vrr(3)%block) .AND. ASSOCIATED(blocks_rrv_vrr(4)%block) .AND. & + ASSOCIATED(blocks_rrv_vrr(5)%block) .AND. ASSOCIATED(blocks_rrv_vrr(6)%block)) + go = go .AND. asso_rrv_vrr + END IF + IF (my_r_rxvr) THEN asso_r_rxvr = .TRUE. DO ind = 1, 3 @@ -998,6 +1049,96 @@ CONTAINS END IF END IF + IF (my_rvr) THEN + ! r_alpha * Vnl * r_beta + IF (iatom <= jatom) THEN + blocks_rvr(1)%block(1:na, 1:nb) = blocks_rvr(1)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 2))) ! xVx + blocks_rvr(2)%block(1:na, 1:nb) = blocks_rvr(2)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! xVy + blocks_rvr(3)%block(1:na, 1:nb) = blocks_rvr(3)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! xVz + blocks_rvr(4)%block(1:na, 1:nb) = blocks_rvr(4)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! yVy + blocks_rvr(5)%block(1:na, 1:nb) = blocks_rvr(5)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! yVz + blocks_rvr(6)%block(1:na, 1:nb) = blocks_rvr(6)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! zVz + ELSE + blocks_rvr(1)%block(1:nb, 1:na) = blocks_rvr(1)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 2))) ! xVx + blocks_rvr(2)%block(1:nb, 1:na) = blocks_rvr(2)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 3))) ! xVy + blocks_rvr(3)%block(1:nb, 1:na) = blocks_rvr(3)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 4))) ! xVz + blocks_rvr(4)%block(1:nb, 1:na) = blocks_rvr(4)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 3))) ! yVy + blocks_rvr(5)%block(1:nb, 1:na) = blocks_rvr(5)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 4))) ! yVz + blocks_rvr(6)%block(1:nb, 1:na) = blocks_rvr(6)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 4))) ! zVz + END IF + END IF + + IF (my_rrv_vrr) THEN + ! r_alpha * r_beta * Vnl + IF (iatom <= jatom) THEN + blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 5), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xxV + blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xyV + blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xzV + blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 8), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yyV + blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yzV + blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 10), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! zzV + ELSE + blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 5), TRANSPOSE(acint(1:na, 1:np, 1))) ! xxV + blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1))) ! xyV + blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1))) ! xzV + blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 8), TRANSPOSE(acint(1:na, 1:np, 1))) ! yyV + blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1))) ! yzV + blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 10), TRANSPOSE(acint(1:na, 1:np, 1))) ! zzV + END IF + ! + Vnl * r_alpha * r_beta + IF (iatom <= jatom) THEN + blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 5))) ! +Vxx + blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 6))) ! +Vxy + blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 7))) ! +Vxz + blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 8))) ! +Vyy + blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 9))) ! +Vyz + blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + & + MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 10))) ! +Vzz + ELSE + blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 5))) ! +Vxx + blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 6))) ! +Vxy + blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 7))) ! +Vxz + blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 8))) ! +Vyy + blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 9))) ! +Vyz + blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + & + MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 10))) ! +Vzz + END IF + END IF + ! The indices are stored in i_1, i_x, ..., i_zzz ! matrix_r_rxvr(alpha, beta) ! = sum_(gamma delta) epsilon_(alpha gamma delta) @@ -1325,6 +1466,18 @@ CONTAINS END DO DEALLOCATE (blocks_rrv) END IF + IF (my_rvr) THEN + DO ind = 1, 6 + NULLIFY (blocks_rvr(ind)%block) + END DO + DEALLOCATE (blocks_rvr) + END IF + IF (my_rrv_vrr) THEN + DO ind = 1, 6 + NULLIFY (blocks_rrv_vrr(ind)%block) + END DO + DEALLOCATE (blocks_rrv_vrr) + END IF IF (my_r_rxvr) THEN DO ind = 1, 3 DO ind2 = 1, 3 diff --git a/src/cp_control_types.F b/src/cp_control_types.F index b62ba184c3..4d4b10c05e 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -66,6 +66,7 @@ MODULE cp_control_types INTEGER, DIMENSION(3) :: delta_pulse_direction REAL(KIND=dp) :: delta_pulse_scale LOGICAL :: velocity_gauge + REAL(KIND=dp), DIMENSION(3) :: vec_pot LOGICAL :: nl_gauge_transform END TYPE rtp_control_type ! ************************************************************************************************** @@ -565,6 +566,7 @@ MODULE cp_control_types dft_plus_u, & apply_efield, & apply_efield_field, & + apply_vector_potential, & apply_period_efield, & apply_external_potential, & eval_external_potential, & diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index bf37d64084..1cdeb9b818 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -399,13 +399,19 @@ CONTAINS ! Read the finite field input section dft_control%apply_efield = .FALSE. dft_control%apply_efield_field = .FALSE. !this is for RTP + dft_control%apply_vector_potential = .FALSE. !this is for RTP tmp_section => section_vals_get_subs_vals(dft_section, "EFIELD") CALL section_vals_get(tmp_section, n_repetition=nrep, explicit=is_present) IF (is_present) THEN ALLOCATE (dft_control%efield_fields(nrep)) CALL read_efield_sections(dft_control, tmp_section) IF (do_rtp) THEN - dft_control%apply_efield_field = .TRUE. + IF (.NOT. dft_control%rtp_control%velocity_gauge) THEN + dft_control%apply_efield_field = .TRUE. + ELSE + dft_control%apply_vector_potential = .TRUE. + dft_control%rtp_control%vec_pot = 0.0_dp + END IF ELSE dft_control%apply_efield = .TRUE. CPASSERT(nrep == 1) diff --git a/src/efield_utils.F b/src/efield_utils.F index f659233d45..4166740b69 100644 --- a/src/efield_utils.F +++ b/src/efield_utils.F @@ -16,24 +16,26 @@ MODULE efield_utils pbc USE cp_control_types, ONLY: dft_control_type,& efield_type - USE cp_para_types, ONLY: cp_para_env_type + USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& + dbcsr_deallocate_matrix_set + USE dbcsr_api, ONLY: dbcsr_add,& + dbcsr_copy,& + dbcsr_p_type,& + dbcsr_set USE input_constants, ONLY: constant_env,& custom_env,& gaussian_env,& ramp_env - USE kahan_sum, ONLY: accurate_dot_product USE kinds, ONLY: dp USE mathconstants, ONLY: pi USE particle_types, ONLY: particle_type - USE pw_types, ONLY: pw_type USE qs_energy_types, ONLY: qs_energy_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,& qs_kind_type - USE qs_rho_types, ONLY: qs_rho_get,& - qs_rho_type + USE qs_moments, ONLY: build_local_moment_matrix #include "./base/base_uses.f90" IMPLICIT NONE @@ -44,75 +46,64 @@ MODULE efield_utils ! *** Public subroutines *** - PUBLIC :: efield_potential, & - calculate_ecore_efield + PUBLIC :: efield_potential_lengh_gauge, & + calculate_ecore_efield, & + make_field CONTAINS ! ************************************************************************************************** -!> \brief computes the time dependend potential on the grid +!> \brief Replace the original implementation of the electric-electronic +!> interaction in the lenght gauge. This calculation is no longer done in +!> the grid but using matrices to match the velocity gauge implementation. +!> Note: The energy is store in energy%core and computed later on. !> \param qs_env ... -!> \param v_efield_rspace ... -!> \author Florian Schiffmann (02.09) +!> \author Guillaume Le Breton (02.23) ! ************************************************************************************************** - SUBROUTINE efield_potential(qs_env, v_efield_rspace) + SUBROUTINE efield_potential_lengh_gauge(qs_env) TYPE(qs_environment_type), POINTER :: qs_env - TYPE(pw_type) :: v_efield_rspace - CHARACTER(len=*), PARAMETER :: routineN = 'efield_potential' + CHARACTER(len=*), PARAMETER :: routineN = 'efield_potential_lengh_gauge' - INTEGER :: handle, i, j, k - INTEGER, DIMENSION(2, 3) :: bo_global, bo_local - REAL(kind=dp) :: dvol, efield_ener, field(3) - REAL(kind=dp), DIMENSION(3) :: dr, grid_p - TYPE(cp_para_env_type), POINTER :: para_env + INTEGER :: handle, i, image + REAL(kind=dp) :: field(3) + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h TYPE(dft_control_type), POINTER :: dft_control - TYPE(pw_type), DIMENSION(:), POINTER :: rho_r - TYPE(qs_energy_type), POINTER :: energy - TYPE(qs_rho_type), POINTER :: rho - NULLIFY (dft_control, para_env, rho_r) + NULLIFY (dft_control) CALL timeset(routineN, handle) CALL get_qs_env(qs_env, & - energy=energy, & - rho=rho, & dft_control=dft_control, & - para_env=para_env) + matrix_h_kp=matrix_h, & + matrix_s=matrix_s) - CALL qs_rho_get(rho, rho_r=rho_r) + NULLIFY (moments) + CALL dbcsr_allocate_matrix_set(moments, 3) + DO i = 1, 3 + ALLOCATE (moments(i)%matrix) + CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments") + CALL dbcsr_set(moments(i)%matrix, 0.0_dp) + END DO - v_efield_rspace%cr3d = 0.0_dp - - bo_local = v_efield_rspace%pw_grid%bounds_local - bo_global = v_efield_rspace%pw_grid%bounds - - dvol = v_efield_rspace%pw_grid%dvol - dr = v_efield_rspace%pw_grid%dr + CALL build_local_moment_matrix(qs_env, moments, 1) CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time) - DO k = bo_local(1, 3), bo_local(2, 3) - DO j = bo_local(1, 2), bo_local(2, 2) - DO i = bo_local(1, 1), bo_local(2, 1) - grid_p(1) = (i - bo_global(1, 1))*dr(1) - grid_p(2) = (j - bo_global(1, 2))*dr(2) - grid_p(3) = (k - bo_global(1, 3))*dr(3) - v_efield_rspace%cr3d(i, j, k) = v_efield_rspace%cr3d(i, j, k) + DOT_PRODUCT(field(:), grid_p(:)) - END DO + DO i = 1, 3 + DO image = 1, dft_control%nimages + CALL dbcsr_add(matrix_h(1, image)%matrix, moments(i)%matrix, 1.0_dp, field(i)) END DO END DO - efield_ener = 0.0_dp - DO i = 1, dft_control%nspins - efield_ener = efield_ener + accurate_dot_product(v_efield_rspace%cr3d, rho_r(i)%cr3d)*dvol - END DO - CALL para_env%group%sum(efield_ener) - energy%efield = efield_ener + + CALL dbcsr_deallocate_matrix_set(moments) + CALL timestop(handle) - END SUBROUTINE efield_potential + END SUBROUTINE efield_potential_lengh_gauge ! ************************************************************************************************** !> \brief computes the amplitude of the efield within a given envelop @@ -124,10 +115,10 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE make_field(dft_control, field, sim_step, sim_time) - TYPE(dft_control_type) :: dft_control - REAL(dp) :: field(3) - INTEGER :: sim_step - REAL(KIND=dp) :: sim_time + TYPE(dft_control_type), INTENT(IN) :: dft_control + REAL(dp), INTENT(OUT) :: field(3) + INTEGER, INTENT(IN) :: sim_step + REAL(KIND=dp), INTENT(IN) :: sim_time INTEGER :: i, lower, nfield, upper REAL(dp) :: c, env, nu, pol(3), strength @@ -182,7 +173,10 @@ CONTAINS END SUBROUTINE make_field ! ************************************************************************************************** -!> \brief computes the force and the energy due to a efield on the cores +!> \brief Computes the force and the energy due to a efield on the cores +!> Note: In the velocity gauge, the energy term is not added because +!> it would lead to an unbalanced energy (center of negative charge not +!> computed in velocity gauge). !> \param qs_env ... !> \param calculate_forces ... !> \author Florian Schiffmann (02.09) @@ -212,7 +206,6 @@ CONTAINS CALL timeset(routineN, handle) CALL get_qs_env(qs_env, dft_control=dft_control) IF (dft_control%apply_efield_field) THEN - my_force = .FALSE. IF (PRESENT(calculate_forces)) my_force = calculate_forces diff --git a/src/emd/rt_delta_pulse.F b/src/emd/rt_delta_pulse.F index 17804ef2c1..7174b01537 100644 --- a/src/emd/rt_delta_pulse.F +++ b/src/emd/rt_delta_pulse.F @@ -273,6 +273,10 @@ CONTAINS cell%h_inv(3, :)*rtp_control%delta_pulse_direction(3) kvec = -kvec*twopi*rtp_control%delta_pulse_scale + IF (rtp_control%velocity_gauge) THEN + rtp_control%vec_pot = rtp_control%vec_pot + kvec + END IF + ! struct for fm matrices fm_struct => rtp%ao_ao_fmstruct @@ -417,6 +421,10 @@ CONTAINS ! 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 + IF (rtp_control%velocity_gauge) THEN + rtp_control%vec_pot = rtp_control%vec_pot + kvec + END IF + ! calculate exponentials (= Berry moments) NULLIFY (cosmat, sinmat) ALLOCATE (cosmat, sinmat) diff --git a/src/emd/rt_propagation_methods.F b/src/emd/rt_propagation_methods.F index 204ec2d2e3..a2e109e308 100644 --- a/src/emd/rt_propagation_methods.F +++ b/src/emd/rt_propagation_methods.F @@ -18,7 +18,8 @@ MODULE rt_propagation_methods USE cp_cfm_types, ONLY: cp_cfm_create,& cp_cfm_release,& cp_cfm_type - USE cp_control_types, ONLY: rtp_control_type + USE cp_control_types, ONLY: dft_control_type,& + rtp_control_type USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,& cp_dbcsr_cholesky_invert USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,& @@ -45,6 +46,7 @@ MODULE rt_propagation_methods 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_type_antisymmetric + USE efield_utils, ONLY: efield_potential_lengh_gauge USE input_constants, ONLY: do_arnoldi,& do_bch,& do_em,& @@ -73,7 +75,8 @@ 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 + USE rt_propagation_velocity_gauge, ONLY: update_vector_potential,& + velocity_gauge_ks_matrix #include "../base/base_uses.f90" IMPLICIT NONE @@ -112,6 +115,7 @@ CONTAINS TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: delta_P, H_last_iter, ks_mix, ks_mix_im, & matrix_ks, matrix_ks_im, matrix_s, & rho_new + TYPE(dft_control_type), POINTER :: dft_control CALL timeset(routineN, handle) @@ -125,15 +129,21 @@ CONTAINS NULLIFY (delta_P, rho_new, delta_mos, mos_new) NULLIFY (ks_mix, ks_mix_im) ! get everything needed and set some values - CALL get_qs_env(qs_env, matrix_s=matrix_s) + CALL get_qs_env(qs_env, matrix_s=matrix_s, dft_control=dft_control) + IF (rtp%iter == 1) THEN CALL qs_energies_init(qs_env, .FALSE.) !the above recalculates matrix_s, but matrix not changed if ions are fixed IF (rtp_control%fixed_ions) CALL set_ks_env(qs_env%ks_env, s_mstruct_changed=.FALSE.) - ! add additional terms for the velocity gauge to matrix_h and matrix_h_im + ! add additional terms to matrix_h and matrix_h_im in the case of applied electric field, + ! either in the lengh or velocity gauge. ! should be called imediately after qs_energies_init and before qs_ks_update_qs_env + IF (dft_control%apply_efield_field) & + CALL efield_potential_lengh_gauge(qs_env) IF (rtp_control%velocity_gauge) THEN + IF (dft_control%apply_vector_potential) & + CALL update_vector_potential(qs_env, dft_control) CALL velocity_gauge_ks_matrix(qs_env, subtract_nl_term=.FALSE.) END IF @@ -288,7 +298,8 @@ CONTAINS CALL timeset(routineN, handle) - CALL get_qs_env(qs_env=qs_env, rtp=rtp, matrix_s=s_mat, matrix_ks=matrix_ks, matrix_ks_im=matrix_ks_im, energy=energy) + CALL get_qs_env(qs_env=qs_env, rtp=rtp, matrix_s=s_mat, & + matrix_ks=matrix_ks, matrix_ks_im=matrix_ks_im, energy=energy) CALL get_rtp(rtp=rtp, exp_H_old=exp_H_old, exp_H_new=exp_H_new) IF (rtp_control%sc_check_start .LT. rtp%iter) THEN @@ -425,8 +436,9 @@ CONTAINS CALL dbcsr_desymmetrize(matrix_ks_im(ispin)%matrix, matrix_ks_nosym) ! take care of the EMD case and add the velocity scaled S_derivative - IF (.NOT. rtp_control%fixed_ions) & + IF (.NOT. rtp_control%fixed_ions) THEN CALL dbcsr_add(matrix_ks_nosym, B_mat, 1.0_dp, -1.0_dp) + END IF CALL dbcsr_multiply("N", "N", -one, S_inv, matrix_ks_nosym, zero, exp_H(re)%matrix, & filter_eps=rtp%filter_eps) @@ -437,7 +449,7 @@ CONTAINS ELSE ! in case of pure EMD its only needed once as B is the same for both spins CALL dbcsr_multiply("N", "N", one, S_inv, B_mat, zero, exp_H(1)%matrix, filter_eps=rtp%filter_eps) - + CALL get_rtp(rtp=rtp, B_mat=B_mat, SinvB=SinvB) CALL dbcsr_copy(SinvB(1)%matrix, exp_H(1)%matrix) IF (SIZE(matrix_ks) == 2) CALL dbcsr_copy(exp_H(3)%matrix, exp_H(1)%matrix) diff --git a/src/emd/rt_propagation_output.F b/src/emd/rt_propagation_output.F index 7f1b7c67d6..2e5bbe7c8d 100644 --- a/src/emd/rt_propagation_output.F +++ b/src/emd/rt_propagation_output.F @@ -42,6 +42,7 @@ MODULE rt_propagation_output dbcsr_get_occupation, dbcsr_init_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, & dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_scale, & dbcsr_set, dbcsr_type + USE efield_utils, ONLY: make_field USE input_constants, ONLY: ehrenfest,& real_time_propagation USE input_section_types, ONLY: section_get_ivals,& @@ -198,11 +199,14 @@ CONTAINS "PRINT%PROGRAM_RUN_INFO") IF (rtp%converged) THEN + dft_section => section_vals_get_subs_vals(input, "DFT") + IF (BTEST(cp_print_key_should_output(logger%iter_info, & + dft_section, "REAL_TIME_PROPAGATION%PRINT%FIELD"), cp_p_file)) & + CALL print_field_applied(qs_env, dft_section) CALL make_moment(qs_env) IF (.NOT. dft_control%qs_control%dftb) THEN CALL write_available_results(qs_env=qs_env, rtp=rtp) END IF - dft_section => section_vals_get_subs_vals(input, "DFT") IF (rtp%linear_scaling) THEN CALL get_rtp(rtp=rtp, rho_new=rho_new) IF (BTEST(cp_print_key_should_output(logger%iter_info, & @@ -798,4 +802,64 @@ CONTAINS END SUBROUTINE write_available_results +! ************************************************************************************************** +!> \brief Print the field applied to the system. Either the electric +!> field or the vector potential depending on the gauge used +!> \param qs_env ... +!> \param dft_section ... +!> \par History +!> 2023-01 Created [Guillaume Le Breton] +! ************************************************************************************************** + SUBROUTINE print_field_applied(qs_env, dft_section) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(section_vals_type), POINTER :: dft_section + + CHARACTER(LEN=3), DIMENSION(3) :: rlab + CHARACTER(LEN=default_path_length) :: filename + INTEGER :: i, output_unit, unit_nr + REAL(kind=dp) :: field(3), to_write(3) + TYPE(cp_logger_type), POINTER :: logger + TYPE(dft_control_type), POINTER :: dft_control + + NULLIFY (dft_control) + + logger => cp_get_default_logger() + output_unit = cp_logger_get_default_io_unit(logger) + + CALL get_qs_env(qs_env, dft_control=dft_control) + + unit_nr = cp_print_key_unit_nr(logger, dft_section, & + "REAL_TIME_PROPAGATION%PRINT%FIELD", extension=".dat") + + IF (output_unit > 0) THEN + IF (unit_nr /= output_unit) THEN + INQUIRE (UNIT=unit_nr, NAME=filename) + WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") & + "FIELD", "The field applied is written to the file:", & + TRIM(filename) + WRITE (UNIT=unit_nr, FMT="(/,(T2,A,T40,I6))") & + "Real time propagation step:", qs_env%sim_step + ELSE + WRITE (UNIT=output_unit, FMT="(/,T2,A)") "FIELD APPLIED" + END IF + END IF + + rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"] + + IF (dft_control%apply_efield_field) THEN + WRITE (unit_nr, "(T3,A)") "Electric Field in atomic units:" + CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time) + to_write = field + ELSE IF (dft_control%apply_vector_potential) THEN + WRITE (unit_nr, "(T3,A)") "Vector potential in atomic units:" + to_write = dft_control%rtp_control%vec_pot + ELSE + to_write = 0._dp + END IF + + WRITE (unit_nr, "(T5,3(A,A,E16.8,1X))") & + (TRIM(rlab(i)), "=", to_write(i), i=1, 3) + + END SUBROUTINE print_field_applied + END MODULE rt_propagation_output diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index e5255aeac9..65dc7e79d3 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -8482,6 +8482,16 @@ CONTAINS CALL section_add_subsection(print_section, print_key) CALL section_release(print_key) + CALL cp_print_key_section_create(print_key, __LOCATION__, "FIELD", & + description="Print the time-dependent field applied during an EMD simulation in "// & + "atomic unit.", & + print_level=high_print_level, common_iter_levels=-1, & + each_iter_names=s2a("MD"), & + each_iter_values=(/1/), & + filename="applied_field") + CALL section_add_subsection(print_section, print_key) + CALL section_release(print_key) + CALL cp_print_key_section_create(print_key, __LOCATION__, "CURRENT", & description="Print the current during an EMD simulation to cube files.", & print_level=high_print_level, common_iter_levels=0, & diff --git a/src/motion/rt_propagation.F b/src/motion/rt_propagation.F index c93f0c7adb..0d0b205889 100644 --- a/src/motion/rt_propagation.F +++ b/src/motion/rt_propagation.F @@ -349,7 +349,6 @@ CONTAINS qs_env%sim_step = i_step rtp%istep = i_step - rtp%i_start CALL calculate_ecore_efield(qs_env, .FALSE.) - ! IF (dft_control%apply_external_potential) THEN IF (.NOT. dft_control%expot_control%static) THEN dft_control%eval_external_potential = .TRUE. diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F index 43b2a191ad..765894f467 100644 --- a/src/qs_ks_methods.F +++ b/src/qs_ks_methods.F @@ -47,7 +47,6 @@ MODULE qs_ks_methods dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, & dbcsr_type_symmetric USE dft_plus_u, ONLY: plus_u - USE efield_utils, ONLY: efield_potential USE hartree_local_methods, ONLY: Vh_1c_gg_integrals USE hartree_local_types, ONLY: ecoul_1center_type USE hfx_admm_utils, ONLY: hfx_admm_init,& @@ -209,8 +208,8 @@ CONTAINS TYPE(pw_type), DIMENSION(:), POINTER :: rho_r, v_rspace_embed, v_rspace_new, & v_rspace_new_aux_fit, v_tau_rspace, & v_tau_rspace_aux_fit - TYPE(pw_type), POINTER :: rho0_s_rs, rho_core, rho_nlcc, v_efield_rspace, v_hartree_rspace, & - v_sccs_rspace, v_sic_rspace, v_spin_ddapc_rest_r, vee, vppl_rspace + TYPE(pw_type), POINTER :: rho0_s_rs, rho_core, rho_nlcc, v_hartree_rspace, v_sccs_rspace, & + v_sic_rspace, v_spin_ddapc_rest_r, vee, vppl_rspace TYPE(qs_energy_type), POINTER :: energy TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(qs_rho_type), POINTER :: rho, rho_struct, rho_xc @@ -424,17 +423,6 @@ CONTAINS END IF CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace) - IF (dft_control%apply_efield_field) THEN - NULLIFY (v_efield_rspace) - ALLOCATE (v_efield_rspace) - CALL pw_pool_create_pw(auxbas_pw_pool, & - v_efield_rspace, & - use_data=REALDATA3D, & - in_space=REALSPACE) - CALL efield_potential(qs_env, v_efield_rspace) - CALL pw_scale(v_efield_rspace, v_efield_rspace%pw_grid%dvol) - END IF - IF (dft_control%correct_surf_dip) THEN IF (dft_control%surf_dip_correct_switch) THEN CALL calc_dipsurf_potential(qs_env, energy) @@ -691,7 +679,7 @@ CONTAINS CALL sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, & v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, & - v_efield_rspace, v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, & + v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, & cdft_control, calculate_forces) IF (use_virial .AND. calculate_forces) THEN @@ -736,11 +724,6 @@ CONTAINS DEALLOCATE (v_spin_ddapc_rest_r) END IF - IF (dft_control%apply_efield_field) THEN - CALL pw_pool_give_back_pw(auxbas_pw_pool, v_efield_rspace) - DEALLOCATE (v_efield_rspace) - END IF - IF (calculate_forces .AND. dft_control%qs_control%cdft) THEN IF (.NOT. cdft_control%transfer_pot) THEN DO iatom = 1, SIZE(cdft_control%group) @@ -1136,6 +1119,7 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'evaluate_core_matrix_traces' INTEGER :: handle + REAL(KIND=dp) :: energy_core_im TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_t, rho_ao_kp TYPE(dft_control_type), POINTER :: dft_control TYPE(qs_energy_type), POINTER :: energy @@ -1155,6 +1139,16 @@ CONTAINS CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy%core, dft_control%nspins) + ! Add the imaginary part in the RTP case + IF (qs_env%run_rtp) THEN + IF (dft_control%rtp_control%velocity_gauge) THEN + CALL get_qs_env(qs_env, matrix_h_im_kp=matrixkp_h) + CALL qs_rho_get(rho, rho_ao_im_kp=rho_ao_kp) + CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy_core_im, dft_control%nspins) + energy%core = energy%core - energy_core_im + END IF + END IF + ! kinetic energy IF (ASSOCIATED(matrixkp_t)) & CALL calculate_ptrace(matrixkp_t, rho_ao_kp, energy%kinetic, dft_control%nspins) diff --git a/src/qs_ks_utils.F b/src/qs_ks_utils.F index 0495f4ef28..1cf3bbcedd 100644 --- a/src/qs_ks_utils.F +++ b/src/qs_ks_utils.F @@ -1254,7 +1254,6 @@ CONTAINS !> \param v_rspace_new_aux_fit ... !> \param v_tau_rspace ... !> \param v_tau_rspace_aux_fit ... -!> \param v_efield_rspace ... !> \param v_sic_rspace ... !> \param v_spin_ddapc_rest_r ... !> \param v_sccs_rspace ... @@ -1269,7 +1268,7 @@ CONTAINS SUBROUTINE sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, & vppl_rspace, v_rspace_new, & v_rspace_new_aux_fit, v_tau_rspace, & - v_tau_rspace_aux_fit, v_efield_rspace, & + v_tau_rspace_aux_fit, & v_sic_rspace, v_spin_ddapc_rest_r, & v_sccs_rspace, v_rspace_embed, cdft_control, & calculate_forces) @@ -1281,8 +1280,8 @@ CONTAINS TYPE(pw_type), POINTER :: vppl_rspace TYPE(pw_type), DIMENSION(:), POINTER :: v_rspace_new, v_rspace_new_aux_fit, & v_tau_rspace, v_tau_rspace_aux_fit - TYPE(pw_type), POINTER :: v_efield_rspace, v_sic_rspace, & - v_spin_ddapc_rest_r, v_sccs_rspace + TYPE(pw_type), POINTER :: v_sic_rspace, v_spin_ddapc_rest_r, & + v_sccs_rspace TYPE(pw_type), DIMENSION(:), POINTER :: v_rspace_embed TYPE(cdft_control_type), POINTER :: cdft_control LOGICAL, INTENT(in) :: calculate_forces @@ -1413,11 +1412,6 @@ CONTAINS cdft_control%strength(igroup) END DO END IF - ! The efield contribution - IF (dft_control%apply_efield_field) THEN - v_rspace_new(ispin)%cr3d = v_rspace_new(ispin)%cr3d + & - v_efield_rspace%cr3d - END IF ! functional derivative of the Hartree energy wrt the density in the presence of dielectric ! (vhartree + v_eps); v_eps is nonzero only if the dielectric constant is defind as a function ! of the charge density @@ -1522,10 +1516,6 @@ CONTAINS CPASSERT(dft_control%sic_method_id == sic_none) CPASSERT(.NOT. dft_control%qs_control%ddapc_restraint_is_spin) DO ispin = 1, nspins - ! the efield contribution - IF (dft_control%apply_efield_field) THEN - v_rspace%cr3d = v_rspace%cr3d + v_efield_rspace%cr3d - END IF ! extra contribution attributed to the dependency of the dielectric constant to the charge density IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol diff --git a/src/qs_mo_io.F b/src/qs_mo_io.F index 73c27f9a8c..80544afabd 100644 --- a/src/qs_mo_io.F +++ b/src/qs_mo_io.F @@ -820,6 +820,7 @@ CONTAINS CALL group%bcast(mos(ispin)%eigenvalues, source) CALL group%bcast(mos(ispin)%uniform_occupation, source) CALL group%bcast(mos(ispin)%occupation_numbers, source) + IF (PRESENT(rt_mos)) THEN DO imat = 2*ispin - 1, 2*ispin DO i = 1, nmo diff --git a/src/qs_moments.F b/src/qs_moments.F index 2453d58eb4..b47d282b7e 100644 --- a/src/qs_moments.F +++ b/src/qs_moments.F @@ -2249,7 +2249,8 @@ CONTAINS order LOGICAL :: my_com_nl, my_velreprs REAL(dp) :: charge, dd, strace, trace - REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rv, nlcom_rxrv, & + REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, & + nlcom_rv, nlcom_rvr, nlcom_rxrv, & qupole_der, rmom_vel REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom REAL(dp), DIMENSION(3) :: rcc, ria @@ -2305,13 +2306,17 @@ CONTAINS IF ((nmom >= 2) .AND. my_velreprs) THEN ALLOCATE (nlcom_rrv(6)) nlcom_rrv(:) = 0._dp + ALLOCATE (nlcom_rvr(6)) + nlcom_rvr(:) = 0._dp + ALLOCATE (nlcom_rrv_vrr(6)) + nlcom_rrv_vrr(:) = 0._dp END IF IF (magnetic) THEN ALLOCATE (nlcom_rxrv(3)) nlcom_rxrv = 0._dp END IF ! Calculate non local correction terms - CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, rcc) + CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, rcc) END IF NULLIFY (moments) @@ -2578,10 +2583,14 @@ CONTAINS mmom(:) = nlcom_rxrv(:) END IF IF (my_velreprs .AND. (nmom >= 1)) THEN + DEALLOCATE (rmom_vel) + ALLOCATE (rmom_vel(21)) rmom_vel(1:3) = nlcom_rv END IF IF (my_velreprs .AND. (nmom >= 2)) THEN rmom_vel(4:9) = nlcom_rrv + rmom_vel(10:15) = nlcom_rvr + rmom_vel(16:21) = nlcom_rrv_vrr END IF IF (magnetic .AND. .NOT. my_velreprs) THEN CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom) @@ -2595,7 +2604,11 @@ CONTAINS IF (my_com_nl) THEN IF (nmom >= 1 .AND. my_velreprs) DEALLOCATE (nlcom_rv) - IF (nmom >= 2 .AND. my_velreprs) DEALLOCATE (nlcom_rrv) + IF (nmom >= 2 .AND. my_velreprs) THEN + DEALLOCATE (nlcom_rrv) + DEALLOCATE (nlcom_rvr) + DEALLOCATE (nlcom_rrv_vrr) + END IF IF (magnetic) DEALLOCATE (nlcom_rxrv) END IF @@ -2789,6 +2802,16 @@ CONTAINS (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=4, 6) WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") & (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=7, 9) + WRITE (unit_number, "(T3,A)") "Expectation value of r x V_nl x r [a. u.]" + WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") & + (TRIM(rlab(i + 1 - 6)), "=", rmom_vel(i), i=10, 12) + WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") & + (TRIM(rlab(i + 1 - 6)), "=", rmom_vel(i), i=13, 15) + WRITE (unit_number, "(T3,A)") "Expectation value of r x r x V_nl + V_nl x r x r [a. u.]" + WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") & + (TRIM(rlab(i + 1 - 12)), "=", rmom_vel(i), i=16, 18) + WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") & + (TRIM(rlab(i + 1 - 12)), "=", rmom_vel(i), i=19, 21) CASE DEFAULT END SELECT END DO @@ -2798,27 +2821,41 @@ CONTAINS END SUBROUTINE print_moments_nl ! ************************************************************************************************** -!> \brief ... +!> \brief Calculate the expectation value of operators related to non-local potential: +!> [r, Vnl], noted rv +!> r x [r,Vnl], noted rxrv +!> [rr,Vnl], noted rrv +!> r x Vnl x r, noted rvr +!> r x r x Vnl + Vnl x r x r, noted rrv_vrr +!> Note that the 3 first operator are commutator while the 2 last +!> are not. For reading clarity the same notation is used for all 5 +!> operators. !> \param qs_env ... !> \param nlcom_rv ... !> \param nlcom_rxrv ... !> \param nlcom_rrv ... +!> \param nlcom_rvr ... +!> \param nlcom_rrv_vrr ... !> \param ref_point ... ! ************************************************************************************************** - SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, ref_point) + SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, & + nlcom_rrv_vrr, ref_point) TYPE(qs_environment_type), POINTER :: qs_env - REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: nlcom_rv, nlcom_rxrv, nlcom_rrv + REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: nlcom_rv, nlcom_rxrv, nlcom_rrv, & + nlcom_rvr, nlcom_rrv_vrr REAL(dp), DIMENSION(3) :: ref_point CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_commutator_nl_terms' INTEGER :: handle, ind, ispin - LOGICAL :: calc_rrv, calc_rv, calc_rxrv + LOGICAL :: calc_rrv, calc_rrv_vrr, calc_rv, & + calc_rvr, calc_rxrv REAL(dp) :: eps_ppnl, strace, trace TYPE(cell_type), POINTER :: cell - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rrv, matrix_rv, matrix_rxrv, & - matrix_s, rho_ao + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rrv, matrix_rrv_vrr, matrix_rv, & + matrix_rvr, matrix_rxrv, matrix_s, & + rho_ao TYPE(dbcsr_type), POINTER :: tmp_ao TYPE(dft_control_type), POINTER :: dft_control TYPE(neighbor_list_set_p_type), DIMENSION(:), & @@ -2832,162 +2869,211 @@ CONTAINS calc_rv = .FALSE. calc_rxrv = .FALSE. calc_rrv = .FALSE. + calc_rvr = .FALSE. + calc_rrv_vrr = .FALSE. - IF (ALLOCATED(nlcom_rv)) calc_rv = .TRUE. - IF (ALLOCATED(nlcom_rrv)) calc_rrv = .TRUE. - IF (ALLOCATED(nlcom_rxrv)) calc_rxrv = .TRUE. + ! rv, rxrv and rrv are commutator matrices: anti-symmetric. + ! The real part of the density matrix rho_ao is symmetric so that + ! the expectation value of real density matrix is zero. Hence, if + ! the density matrix is real, no need to compute these quantities. + ! This is not the case for rvr and rrv_vrr which are symmetric. - IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv)) RETURN - - ! calculate expectation values - ! real part - ! This evaluates to zero, because all commutator matrices are anti-symmetric while - ! the real part of the density matrix rho_ao is symmetric. Tr[A*S] = 0 - IF (calc_rv) THEN + IF (ALLOCATED(nlcom_rv)) THEN nlcom_rv(:) = 0._dp + IF (qs_env%run_rtp) calc_rv = .TRUE. + END IF + IF (ALLOCATED(nlcom_rxrv)) THEN + nlcom_rxrv(:) = 0._dp + IF (qs_env%run_rtp) calc_rxrv = .TRUE. + END IF + IF (ALLOCATED(nlcom_rrv)) THEN + nlcom_rrv(:) = 0._dp + IF (qs_env%run_rtp) calc_rrv = .TRUE. + END IF + IF (ALLOCATED(nlcom_rvr)) THEN + nlcom_rvr(:) = 0._dp + calc_rvr = .TRUE. + END IF + IF (ALLOCATED(nlcom_rrv_vrr)) THEN + nlcom_rrv_vrr(:) = 0._dp + calc_rrv_vrr = .TRUE. END IF - IF (calc_rrv) THEN - nlcom_rrv(:) = 0._dp + IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv & + .OR. calc_rvr .OR. calc_rrv_vrr)) RETURN + + 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, & + matrix_s=matrix_s, & + particle_set=particle_set, & + 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 + + ! Allocate storage + NULLIFY (matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr) + IF (calc_rv) THEN + CALL dbcsr_allocate_matrix_set(matrix_rv, 3) + DO ind = 1, 3 + CALL dbcsr_init_p(matrix_rv(ind)%matrix) + CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, & + matrix_type=dbcsr_type_antisymmetric) + CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb) + CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp) + END DO END IF IF (calc_rxrv) THEN - nlcom_rxrv(:) = 0._dp + CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3) + DO ind = 1, 3 + CALL dbcsr_init_p(matrix_rxrv(ind)%matrix) + CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, & + matrix_type=dbcsr_type_antisymmetric) + CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb) + CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp) + END DO + END IF + + IF (calc_rrv) THEN + CALL dbcsr_allocate_matrix_set(matrix_rrv, 6) + DO ind = 1, 6 + CALL dbcsr_init_p(matrix_rrv(ind)%matrix) + CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, & + matrix_type=dbcsr_type_antisymmetric) + CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb) + CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp) + END DO + END IF + + IF (calc_rvr) THEN + CALL dbcsr_allocate_matrix_set(matrix_rvr, 6) + DO ind = 1, 6 + CALL dbcsr_init_p(matrix_rvr(ind)%matrix) + CALL dbcsr_create(matrix_rvr(ind)%matrix, template=matrix_s(1)%matrix, & + matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(matrix_rvr(ind)%matrix, sab_orb) + CALL dbcsr_set(matrix_rvr(ind)%matrix, 0._dp) + END DO + END IF + IF (calc_rrv_vrr) THEN + CALL dbcsr_allocate_matrix_set(matrix_rrv_vrr, 6) + DO ind = 1, 6 + CALL dbcsr_init_p(matrix_rrv_vrr(ind)%matrix) + CALL dbcsr_create(matrix_rrv_vrr(ind)%matrix, template=matrix_s(1)%matrix, & + matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv_vrr(ind)%matrix, sab_orb) + CALL dbcsr_set(matrix_rrv_vrr(ind)%matrix, 0._dp) + END DO + END IF + + ! calculate evaluation of operators in AO basis set + 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, matrix_rvr=matrix_rvr, & + matrix_rrv_vrr=matrix_rrv_vrr, ref_point=ref_point) + + ! Calculate expectation values + ! Real part + NULLIFY (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_rvr .OR. calc_rrv_vrr) THEN + NULLIFY (rho_ao) + CALL qs_rho_get(rho, rho_ao=rho_ao) + + IF (calc_rvr) THEN + trace = 0._dp + DO ind = 1, SIZE(matrix_rvr) + strace = 0._dp + DO ispin = 1, dft_control%nspins + CALL dbcsr_set(tmp_ao, 0.0_dp) + CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rvr(ind)%matrix, & + 0.0_dp, tmp_ao) + CALL dbcsr_trace(tmp_ao, trace) + strace = strace + trace + END DO + nlcom_rvr(ind) = nlcom_rvr(ind) + strace + END DO + END IF + + IF (calc_rrv_vrr) THEN + trace = 0._dp + DO ind = 1, SIZE(matrix_rrv_vrr) + strace = 0._dp + DO ispin = 1, dft_control%nspins + CALL dbcsr_set(tmp_ao, 0.0_dp) + CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv_vrr(ind)%matrix, & + 0.0_dp, tmp_ao) + CALL dbcsr_trace(tmp_ao, trace) + strace = strace + trace + END DO + nlcom_rrv_vrr(ind) = nlcom_rrv_vrr(ind) + strace + END DO + END IF END IF ! imagninary part of the density matrix - IF (qs_env%run_rtp) THEN - 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, & - matrix_s=matrix_s, & - particle_set=particle_set, & - qs_kind_set=qs_kind_set, & - rho=rho, & - sab_orb=sab_orb, & - sab_all=sab_all, & - sap_ppnl=sap_ppnl) + NULLIFY (rho_ao) + CALL qs_rho_get(rho, rho_ao_im=rho_ao) - eps_ppnl = dft_control%qs_control%eps_ppnl - ! Calculate commutators, only needed of there is an imaginary part of the density matrix - ! Allocate storage - NULLIFY (matrix_rv, matrix_rrv, matrix_rxrv) - IF (calc_rv) THEN - CALL dbcsr_allocate_matrix_set(matrix_rv, 3) - DO ind = 1, 3 - CALL dbcsr_init_p(matrix_rv(ind)%matrix) - CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, & - matrix_type=dbcsr_type_antisymmetric) - CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb) - CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp) + IF (calc_rv) THEN + trace = 0._dp + DO ind = 1, SIZE(matrix_rv) + strace = 0._dp + DO ispin = 1, dft_control%nspins + CALL dbcsr_set(tmp_ao, 0.0_dp) + CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, & + 0.0_dp, tmp_ao) + CALL dbcsr_trace(tmp_ao, trace) + strace = strace + trace END DO - END IF + nlcom_rv(ind) = nlcom_rv(ind) + strace + END DO + END IF - IF (calc_rrv) THEN - CALL dbcsr_allocate_matrix_set(matrix_rrv, 6) - DO ind = 1, 6 - CALL dbcsr_init_p(matrix_rrv(ind)%matrix) - CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, & - matrix_type=dbcsr_type_antisymmetric) - CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb) - CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp) + IF (calc_rrv) THEN + trace = 0._dp + DO ind = 1, SIZE(matrix_rrv) + strace = 0._dp + DO ispin = 1, dft_control%nspins + CALL dbcsr_set(tmp_ao, 0.0_dp) + CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, & + 0.0_dp, tmp_ao) + CALL dbcsr_trace(tmp_ao, trace) + strace = strace + trace END DO - END IF + nlcom_rrv(ind) = nlcom_rrv(ind) + strace + END DO + END IF - IF (calc_rxrv) THEN - CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3) - DO ind = 1, 3 - CALL dbcsr_init_p(matrix_rxrv(ind)%matrix) - CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, & - matrix_type=dbcsr_type_antisymmetric) - CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb) - CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp) + IF (calc_rxrv) THEN + trace = 0._dp + DO ind = 1, SIZE(matrix_rxrv) + strace = 0._dp + DO ispin = 1, dft_control%nspins + CALL dbcsr_set(tmp_ao, 0.0_dp) + CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, & + 0.0_dp, tmp_ao) + CALL dbcsr_trace(tmp_ao, trace) + strace = strace + trace END DO - END IF - - ! 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, 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, 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, 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, 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, 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, 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, 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) - 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 - DO ind = 1, SIZE(matrix_rv) - strace = 0._dp - DO ispin = 1, dft_control%nspins - CALL dbcsr_set(tmp_ao, 0.0_dp) - CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, & - 0.0_dp, tmp_ao) - CALL dbcsr_trace(tmp_ao, trace) - strace = strace + trace - END DO - nlcom_rv(ind) = nlcom_rv(ind) + strace - END DO - END IF - - IF (calc_rrv) THEN - trace = 0._dp - DO ind = 1, SIZE(matrix_rrv) - strace = 0._dp - DO ispin = 1, dft_control%nspins - CALL dbcsr_set(tmp_ao, 0.0_dp) - CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, & - 0.0_dp, tmp_ao) - CALL dbcsr_trace(tmp_ao, trace) - strace = strace + trace - END DO - nlcom_rrv(ind) = nlcom_rrv(ind) + strace - END DO - END IF - - IF (calc_rxrv) THEN - trace = 0._dp - DO ind = 1, SIZE(matrix_rxrv) - strace = 0._dp - DO ispin = 1, dft_control%nspins - CALL dbcsr_set(tmp_ao, 0.0_dp) - CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, & - 0.0_dp, tmp_ao) - CALL dbcsr_trace(tmp_ao, trace) - strace = strace + trace - END DO - nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace - END DO - END IF - CALL dbcsr_deallocate_matrix(tmp_ao) - IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv) - IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv) - IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv) - END IF ! qs_env%run_rtp + nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace + END DO + END IF + CALL dbcsr_deallocate_matrix(tmp_ao) + IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv) + IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv) + IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv) + IF (calc_rvr) CALL dbcsr_deallocate_matrix_set(matrix_rvr) + IF (calc_rrv_vrr) CALL dbcsr_deallocate_matrix_set(matrix_rrv_vrr) CALL timestop(handle) END SUBROUTINE calculate_commutator_nl_terms diff --git a/src/rt_propagation_velocity_gauge.F b/src/rt_propagation_velocity_gauge.F index a49515aef7..4891fe95f5 100644 --- a/src/rt_propagation_velocity_gauge.F +++ b/src/rt_propagation_velocity_gauge.F @@ -31,6 +31,7 @@ MODULE rt_propagation_velocity_gauge dbcsr_set,& dbcsr_type_antisymmetric,& dbcsr_type_symmetric + USE efield_utils, ONLY: make_field USE external_potential_types, ONLY: gth_potential_p_type,& gth_potential_type,& sgp_potential_p_type,& @@ -41,7 +42,6 @@ MODULE rt_propagation_velocity_gauge USE kpoint_types, ONLY: get_kpoint_info,& kpoint_type USE mathconstants, ONLY: one,& - twopi,& zero USE orbital_pointers, ONLY: init_orbital_pointers,& nco,& @@ -79,7 +79,7 @@ MODULE rt_propagation_velocity_gauge CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_velocity_gauge' - PUBLIC :: velocity_gauge_ks_matrix + PUBLIC :: velocity_gauge_ks_matrix, update_vector_potential CONTAINS @@ -128,6 +128,7 @@ CONTAINS sap_ppnl, particle_set, qs_kind_set, atomic_kind_set, virial, force, matrix_p, rho, matrix_nl) CALL get_qs_env(qs_env, & + rho=rho, & dft_control=dft_control, & sab_orb=sab_orb, & sap_ppnl=sap_ppnl, & @@ -182,8 +183,7 @@ CONTAINS END IF !get vector potential - vec_pot(:) = 0._dp - CALL get_vector_potential(dft_control, cell, vec_pot) + vec_pot = dft_control%rtp_control%vec_pot ! allocate and build matrices for linear momentum term NULLIFY (momentum) @@ -216,10 +216,12 @@ CONTAINS 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 + ! add Non local term IF (ppnl_present) THEN IF (dft_control%rtp_control%nl_gauge_transform) THEN NULLIFY (nl_term) @@ -252,26 +254,21 @@ CONTAINS END SUBROUTINE velocity_gauge_ks_matrix ! ************************************************************************************************** -!> \brief ... +!> \brief Update the vector potential in the case where a time-dependant +!> electric field is apply. +!> \param qs_env ... !> \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 + SUBROUTINE update_vector_potential(qs_env, dft_control) + TYPE(qs_environment_type), INTENT(INOUT), POINTER :: qs_env + TYPE(dft_control_type), INTENT(INOUT), POINTER :: dft_control - ! 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 + REAL(kind=dp) :: field(3) - END SUBROUTINE get_vector_potential + CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time) + dft_control%rtp_control%vec_pot = dft_control%rtp_control%vec_pot - field*qs_env%rtp%dt + + END SUBROUTINE update_vector_potential ! ************************************************************************************************** !> \brief ... diff --git a/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp b/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp index d4110d0c42..a7a72bafad 100644 --- a/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp +++ b/tests/QS/regtest-rtp-2/H2-emd-efield-custom.inp @@ -2,7 +2,7 @@ METHOD QUICKSTEP &DFT &REAL_TIME_PROPAGATION - MAX_ITER 5 + MAX_ITER 25 MAT_EXP TAYLOR EPS_ITER 1.0E-9 INITIAL_WFN SCF_WFN diff --git a/tests/QS/regtest-rtp-2/H2-emd-efield-ramp.inp b/tests/QS/regtest-rtp-2/H2-emd-efield-ramp.inp index 6a36b4387b..9e6d87a451 100644 --- a/tests/QS/regtest-rtp-2/H2-emd-efield-ramp.inp +++ b/tests/QS/regtest-rtp-2/H2-emd-efield-ramp.inp @@ -2,7 +2,7 @@ METHOD QUICKSTEP &DFT &REAL_TIME_PROPAGATION - MAX_ITER 5 + MAX_ITER 25 MAT_EXP TAYLOR EPS_ITER 1.0E-9 INITIAL_WFN SCF_WFN diff --git a/tests/QS/regtest-rtp-2/H2-emd-efield.inp b/tests/QS/regtest-rtp-2/H2-emd-efield.inp index 86367dc022..5afeb128b7 100644 --- a/tests/QS/regtest-rtp-2/H2-emd-efield.inp +++ b/tests/QS/regtest-rtp-2/H2-emd-efield.inp @@ -2,7 +2,7 @@ METHOD QUICKSTEP &DFT &REAL_TIME_PROPAGATION - MAX_ITER 5 + MAX_ITER 25 MAT_EXP TAYLOR EPS_ITER 1.0E-9 INITIAL_WFN SCF_WFN diff --git a/tests/QS/regtest-rtp-2/H2-rtp-efield.inp b/tests/QS/regtest-rtp-2/H2-rtp-efield.inp index 47af471b92..d4f6475cb0 100644 --- a/tests/QS/regtest-rtp-2/H2-rtp-efield.inp +++ b/tests/QS/regtest-rtp-2/H2-rtp-efield.inp @@ -2,7 +2,7 @@ METHOD QUICKSTEP &DFT &REAL_TIME_PROPAGATION - MAX_ITER 5 + MAX_ITER 25 MAT_EXP TAYLOR EPS_ITER 1.0E-9 INITIAL_WFN SCF_WFN diff --git a/tests/QS/regtest-rtp-2/TEST_FILES b/tests/QS/regtest-rtp-2/TEST_FILES index c48eaf8ca4..330f9c6f7a 100644 --- a/tests/QS/regtest-rtp-2/TEST_FILES +++ b/tests/QS/regtest-rtp-2/TEST_FILES @@ -7,10 +7,10 @@ H2-emd_restart.inp 2 1.0E-14 - H2-emd_restart-1.inp 2 1.0E-14 -0.902241027707E+00 H2-rtp_restart.inp 1 3e-13 -0.90223968349591 H2-rtp_restart-1.inp 1 3e-13 -0.90223968361437 -H2-rtp-efield.inp 1 2e-10 -0.8296879654 -H2-emd-efield.inp 2 4e-11 -0.894866854822E+00 -H2-emd-efield-ramp.inp 2 5e-12 -0.886674608727E+00 -H2-emd-efield-custom.inp 2 1e-14 -0.902238072012E+00 +H2-rtp-efield.inp 1 2e-10 -0.82613324527108 +H2-emd-efield.inp 2 4e-11 -0.902242711172 +H2-emd-efield-ramp.inp 2 5e-12 -0.902242710949 +H2-emd-efield-custom.inp 2 1e-14 -0.902242709393 H2-rtp_ETRS_ARNOLDI.inp 1 3e-13 -0.90223968349550 -H2-emd_ETRS_ARNOLDI.inp 1 1e-10 -17.085427015 +H2-emd_ETRS_ARNOLDI.inp 1 1e-10 -17.08557183219799 #EOF diff --git a/tests/QS/regtest-rtp-4/H2O-vg.inp b/tests/QS/regtest-rtp-4/H2O-vg.inp index c56f6539ef..107f892cef 100644 --- a/tests/QS/regtest-rtp-4/H2O-vg.inp +++ b/tests/QS/regtest-rtp-4/H2O-vg.inp @@ -27,11 +27,8 @@ &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 @@ -39,6 +36,17 @@ PROPAGATOR ETRS INITIAL_WFN SCF_WFN ! RT_RESTART &END + &EFIELD + #Equivalent to a pulse scale of 0.001 with this cell size (7 7 7) + INTENSITY 7.917741E+11 + POLARISATION 1 0 0 + ENVELOP CONSTANT + WAVELENGTH [nm] 1000 + &CONSTANT_ENV + START_STEP 0 + END_STEP 1 + &END CONSTANT_ENV + &END EFIELD &PRINT &MOMENTS PERIODIC .FALSE. diff --git a/tests/QS/regtest-rtp-4/TEST_FILES b/tests/QS/regtest-rtp-4/TEST_FILES index 762a477a68..5a4c724f78 100644 --- a/tests/QS/regtest-rtp-4/TEST_FILES +++ b/tests/QS/regtest-rtp-4/TEST_FILES @@ -7,5 +7,5 @@ # 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 +H2O-vg.inp 1 1e-13 -17.17816522919518 +#EOF