mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-29 06:35:28 -04:00
RTP: Add velocity gauge
This commit is contained in:
glb96
GitHub
parent
90e3ec868d
commit
78ff48db3a
21 changed files with 616 additions and 293 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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, &
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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, &
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
382
src/qs_moments.F
382
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
|
||||
|
|
|
|||
|
|
@ -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 ...
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
|
|
@ -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
|
||||
H2O-vg.inp 1 1e-13 -17.17816522919518
|
||||
#EOF
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue