Skip to content

Commit

Permalink
Add keyword and refpoint for non local commutator
Browse files Browse the repository at this point in the history
  • Loading branch information
@mattiatj @dev-zero
mattiatj authored and dev-zero committed Feb 24, 2020
1 parent 1202fe0 commit 5a454a7
Show file tree
Hide file tree
Showing 4 changed files with 85 additions and 26 deletions.
82 changes: 62 additions & 20 deletions src/commutator_rpnl.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@ MODULE commutator_rpnl
USE basis_set_types, ONLY: gto_basis_set_p_type,&
gto_basis_set_type
USE block_p_types, ONLY: block_p_type
USE cell_types, ONLY: cell_type,&
pbc
USE dbcsr_api, ONLY: dbcsr_get_block_p,&
dbcsr_p_type
USE external_potential_types, ONLY: gth_potential_p_type,&
Expand All @@ -25,6 +27,7 @@ MODULE commutator_rpnl
USE orbital_pointers, ONLY: init_orbital_pointers,&
nco,&
ncoset
USE particle_types, ONLY: particle_type
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
Expand Down Expand Up @@ -409,8 +412,11 @@ END SUBROUTINE build_com_rpnl
!> \param matrix_rxrv ...
!> \param matrix_rrv ...
!> \param ref_point ...
!> \param particle_set ...
!> \param cell ...
! **************************************************************************************************
SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl, matrix_rxrv, matrix_rrv, ref_point)
SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl, matrix_rxrv, matrix_rrv, ref_point, &
particle_set, cell)

TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rv
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
Expand All @@ -420,6 +426,9 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: matrix_rxrv, matrix_rrv
REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: ref_point
TYPE(particle_type), DIMENSION(:), OPTIONAL, &
POINTER :: particle_set
TYPE(cell_type), OPTIONAL, POINTER :: cell

CHARACTER(LEN=*), PARAMETER :: routineN = 'build_com_mom_nl', &
routineP = moduleN//':'//routineN
Expand All @@ -433,11 +442,12 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
nsgf_seta
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
LOGICAL :: asso_rrv, asso_rv, asso_rxrv, found, go, &
gpot, ppnl_present, spot
gpot, my_ref, my_rrv, my_rxrv, &
ppnl_present, spot
REAL(KIND=dp) :: dac, ppnl_radius
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work, sab, work
REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc
REAL(KIND=dp), DIMENSION(3) :: rab, rac, rcc
REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, raf, rc, rcf, rf
REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_ppnl, set_radius_a
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, rpgfa, sphi_a, vprj_ppnl, zeta
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
Expand Down Expand Up @@ -505,7 +515,12 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
END IF
END DO

IF (PRESENT(matrix_rxrv) .OR. PRESENT(matrix_rrv)) THEN
my_rxrv = .FALSE.
my_rrv = .FALSE.
IF (PRESENT(matrix_rxrv)) my_rxrv = .TRUE.
IF (PRESENT(matrix_rrv)) my_rrv = .TRUE.

IF (my_rxrv .OR. my_rrv) THEN
maxder = 10
order = 2
CPASSERT(PRESENT(ref_point))
Expand All @@ -514,8 +529,13 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
order = 1
ENDIF

my_ref = .FALSE.
IF (PRESENT(ref_point)) THEN
rcc = ref_point
CPASSERT((PRESENT(cell) .AND. PRESENT(particle_set)))
rf = ref_point
my_ref = .TRUE.
ELSE
rf(:) = 0._dp
ENDIF

nthread = 1
Expand All @@ -526,12 +546,12 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP SHARED (nl_iterator, basis_set, spotential, gpotential, maxder, ncoset, &
!$OMP sap_int, nkind, ldsab, ldai, nco, order) &
!$OMP sap_int, nkind, ldsab, ldai, nco, cell, particle_set, order, rf, my_ref) &
!$OMP PRIVATE (mepos, ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
!$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, &
!$OMP sphi_a, zeta, cprj, lppnl, nppnl, nprj_ppnl, &
!$OMP clist, iset, ncoa, sgfa, prjc, work, sab, ai_work, nprjc, ppnl_radius, &
!$OMP ncoc, rpgfa, vprj_ppnl, i, l, gpot, spot, &
!$OMP ncoc, rpgfa, vprj_ppnl, i, l, gpot, spot, raf, rcf, ra, rc, &
!$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl)
mepos = 0
!$ mepos = omp_get_thread_num()
Expand Down Expand Up @@ -607,6 +627,14 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
clist%achint = 0._dp
clist%nsgf_cnt = 0
NULLIFY (clist%sgf_list)

IF (my_ref) THEN
ra(:) = pbc(particle_set(iatom)%r(:) - rf, cell) + rf
rc(:) = pbc(particle_set(katom)%r(:) - rf, cell) + rf
raf(:) = ra(:) - rf(:)
rcf(:) = rc(:) - rf(:)
ENDIF

DO iset = 1, nseta
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
Expand All @@ -624,8 +652,13 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
! Calculate the primitive overlap and dipole moment integrals
CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab(:, :, 1), 0, .FALSE., ai_work, ldai)
CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
lc_max, 1, zetc, rprjc, order, rac, (/0._dp, 0._dp, 0._dp/), sab(:, :, 2:maxder))
IF (.NOT. my_ref) THEN
CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
lc_max, 1, zetc, rprjc, order, -rac, (/0._dp, 0._dp, 0._dp/), sab(:, :, 2:maxder))
ELSE
CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
lc_max, 1, zetc, rprjc, order, raf, rcf, sab(:, :, 2:maxder))
ENDIF
! *** Transformation step projector functions (cartesian->spherical) ***
DO i = 1, maxder
CALL dgemm("N", "N", ncoa, nprjc, ncoc, 1.0_dp, sab(1, 1, i), ldsab, &
Expand Down Expand Up @@ -659,7 +692,7 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP SHARED (nl_iterator, basis_set, matrix_rv, matrix_rxrv, matrix_rrv, &
!$OMP sap_int, nkind, eps_ppnl ) &
!$OMP sap_int, nkind, eps_ppnl, my_rxrv, my_rrv ) &
!$OMP PRIVATE (mepos, ikind, jkind, iatom, jatom, nlist, ilist, nnode, inode, cell_b, rab, &
!$OMP iab, irow, icol, blocks_rv, blocks_rxrv, blocks_rrv, &
!$OMP found, iac, ibc, alist_ac, alist_bc, acint, bcint, &
Expand All @@ -672,10 +705,10 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
NULLIFY (blocks_rv, blocks_rxrv, blocks_rrv)
ALLOCATE (blocks_rv(3))

IF (PRESENT(matrix_rxrv)) THEN
IF (my_rxrv) THEN
ALLOCATE (blocks_rxrv(3))
ENDIF
IF (PRESENT(matrix_rxrv)) THEN
IF (my_rxrv) THEN
ALLOCATE (blocks_rrv(6))
ENDIF

Expand All @@ -694,13 +727,13 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
CALL dbcsr_get_block_p(matrix_rv(ind)%matrix, irow, icol, blocks_rv(ind)%block, found)
ENDDO

IF (PRESENT(matrix_rxrv)) THEN
IF (my_rxrv) THEN
DO ind = 1, 3
CALL dbcsr_get_block_p(matrix_rxrv(ind)%matrix, irow, icol, blocks_rxrv(ind)%block, found)
ENDDO
ENDIF

IF (PRESENT(matrix_rrv)) THEN
IF (my_rrv) THEN
DO ind = 1, 3
CALL dbcsr_get_block_p(matrix_rrv(ind)%matrix, irow, icol, blocks_rrv(ind)%block, found)
ENDDO
Expand All @@ -711,13 +744,13 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
ASSOCIATED(blocks_rv(3)%block))
go = asso_rv

IF (PRESENT(matrix_rxrv)) THEN
IF (my_rxrv) THEN
asso_rxrv = (ASSOCIATED(blocks_rxrv(1)%block) .AND. ASSOCIATED(blocks_rxrv(2)%block) .AND. &
ASSOCIATED(blocks_rxrv(1)%block))
go = go .AND. asso_rxrv
ENDIF

IF (PRESENT(matrix_rxrv)) THEN
IF (my_rxrv) THEN
asso_rrv = (ASSOCIATED(blocks_rrv(1)%block) .AND. ASSOCIATED(blocks_rrv(2)%block) .AND. &
ASSOCIATED(blocks_rrv(1)%block) .AND. ASSOCIATED(blocks_rrv(2)%block) .AND. &
ASSOCIATED(blocks_rrv(1)%block) .AND. ASSOCIATED(blocks_rrv(2)%block))
Expand Down Expand Up @@ -763,7 +796,7 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
bcint(1, 1, 4), nb, 1.0_dp, blocks_rv(3)%block, SIZE(blocks_rv(3)%block, 1)) ! -Vz

IF (PRESENT(matrix_rxrv)) THEN
IF (my_rxrv) THEN
! x-component (y [z,Vnl] - z [y, Vnl])
CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 9), na, &
bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! yzV
Expand Down Expand Up @@ -795,7 +828,7 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
bcint(1, 1, 2), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! yVx
ENDIF

IF (PRESENT(matrix_rrv)) THEN
IF (my_rrv) THEN
! r_alpha * r_beta * Vnl
CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 5), na, &
bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(1)%block, SIZE(blocks_rrv(1)%block, 1)) ! xxV
Expand Down Expand Up @@ -832,11 +865,20 @@ SUBROUTINE build_com_mom_nl(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl,
END DO
ENDIF
END DO
DO ind = 1, 3
NULLIFY (blocks_rv(ind)%block)
ENDDO
DEALLOCATE (blocks_rv)
IF (PRESENT(matrix_rxrv)) THEN
IF (my_rxrv) THEN
DO ind = 1, 3
NULLIFY (blocks_rxrv(ind)%block)
ENDDO
DEALLOCATE (blocks_rxrv)
ENDIF
IF (PRESENT(matrix_rxrv)) THEN
IF (my_rrv) THEN
DO ind = 1, 6
NULLIFY (blocks_rrv(ind)%block)
ENDDO
DEALLOCATE (blocks_rrv)
ENDIF
!$OMP END PARALLEL
Expand Down
13 changes: 12 additions & 1 deletion src/input_cp2k_dft.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1698,14 +1698,25 @@ SUBROUTINE create_print_dft_section(section)
CALL keyword_create(keyword, __LOCATION__, &
name="VEL_REPRS", &
description="Calculate expectation values the el. multipole moments in their velocity"// &
"representation during RTP. Implemented up to el. quadrupole moment", &
"representation during RTP. Implemented up to el. quadrupole moment.", &
usage="VEL_REPS yes", &
repeats=.FALSE., &
n_var=1, &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, &
name="COM_NL", &
description="Include non local commutator for velocity representations."// &
"Necessary for origin independent results.", &
usage="COM_NL yes", &
repeats=.FALSE., &
n_var=1, &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)

Expand Down
8 changes: 5 additions & 3 deletions src/qs_moments.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1755,15 +1755,16 @@ END SUBROUTINE op_orbbas_rtp
!> \param ref_point ...
!> \param unit_number ...
!> \param vel_reprs ...
!> \param com_nl ...
! **************************************************************************************************
SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs)
SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN) :: magnetic
INTEGER, INTENT(IN) :: nmoments, reference
REAL(dp), DIMENSION(:), POINTER :: ref_point
INTEGER, INTENT(IN) :: unit_number
LOGICAL, INTENT(IN), OPTIONAL :: vel_reprs
LOGICAL, INTENT(IN), OPTIONAL :: vel_reprs, com_nl
CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_moment_locop', &
routineP = moduleN//':'//routineN
Expand All @@ -1773,7 +1774,7 @@ SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, uni
INTEGER :: akind, handle, i, ia, iatom, idir, &
ikind, ispin, ix, iy, iz, l, nm, nmm, &
nmom, order
LOGICAL :: my_velreprs
LOGICAL :: my_com_nl, my_velreprs
REAL(dp) :: charge, dd, strace, trace
REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, qupole_der, rmom_vel
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
Expand All @@ -1800,6 +1801,7 @@ SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, uni
my_velreprs = .FALSE.
IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
IF (PRESENT(com_nl)) my_com_nl = com_nl
! reference point
CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
Expand Down
8 changes: 6 additions & 2 deletions src/qs_scf_post_gpw.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1103,7 +1103,8 @@ SUBROUTINE qs_scf_post_moments(input, logger, qs_env, output_unit)

CHARACTER(LEN=default_path_length) :: filename
INTEGER :: handle, maxmom, reference, unit_nr
LOGICAL :: do_kpoints, magnetic, periodic, vel_reprs
LOGICAL :: com_nl, do_kpoints, magnetic, periodic, &
vel_reprs
REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
TYPE(section_vals_type), POINTER :: print_key

Expand All @@ -1124,6 +1125,9 @@ SUBROUTINE qs_scf_post_moments(input, logger, qs_env, output_unit)
keyword_name="DFT%PRINT%MOMENTS%MAGNETIC")
vel_reprs = section_get_lval(section_vals=input, &
keyword_name="DFT%PRINT%MOMENTS%VEL_REPRS")
com_nl = section_get_lval(section_vals=input, &
keyword_name="DFT%PRINT%MOMENTS%COM_NL")

NULLIFY (ref_point)
CALL section_vals_val_get(input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=input, &
Expand All @@ -1149,7 +1153,7 @@ SUBROUTINE qs_scf_post_moments(input, logger, qs_env, output_unit)
IF (periodic) THEN
CALL qs_moment_berry_phase(qs_env, magnetic, maxmom, reference, ref_point, unit_nr)
ELSE
CALL qs_moment_locop(qs_env, magnetic, maxmom, reference, ref_point, unit_nr, vel_reprs)
CALL qs_moment_locop(qs_env, magnetic, maxmom, reference, ref_point, unit_nr, vel_reprs, com_nl)
END IF
END IF

Expand Down

0 comments on commit 5a454a7

Please sign in to comment.