mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-27 13:45:19 -04:00
Infrastructure for Kpoint symmetries (#3482)
This commit is contained in:
Juerg Hutter
GitHub
parent
a391c5a87d
commit
17f893d634
10 changed files with 3860 additions and 120 deletions
|
|
@ -304,6 +304,7 @@ list(
|
|||
kpoint_methods.F
|
||||
kpoint_transitional.F
|
||||
kpoint_types.F
|
||||
kpsym.F
|
||||
libgrpp_integrals.F
|
||||
libint_2c_3c.F
|
||||
libint_wrapper.F
|
||||
|
|
@ -492,6 +493,7 @@ list(
|
|||
qs_atomic_block.F
|
||||
qs_band_structure.F
|
||||
qs_basis_gradient.F
|
||||
qs_basis_rotation_methods.F
|
||||
qs_block_davidson_types.F
|
||||
qs_cdft_methods.F
|
||||
qs_cdft_opt_types.F
|
||||
|
|
|
|||
|
|
@ -46,16 +46,15 @@ MODULE orbital_transformation_matrices
|
|||
|
||||
TYPE(orbtramat_type), DIMENSION(:), POINTER :: orbtramat => NULL()
|
||||
|
||||
TYPE orbrotmat_type
|
||||
REAL(KIND=dp), DIMENSION(:, :), POINTER :: mat => NULL()
|
||||
END TYPE orbrotmat_type
|
||||
|
||||
INTEGER, SAVE :: current_maxl = -1
|
||||
|
||||
! Public variables
|
||||
|
||||
PUBLIC :: orbtramat
|
||||
|
||||
! Public subroutines
|
||||
|
||||
PUBLIC :: deallocate_spherical_harmonics, &
|
||||
init_spherical_harmonics
|
||||
PUBLIC :: orbrotmat_type, calculate_rotmat, release_rotmat
|
||||
PUBLIC :: deallocate_spherical_harmonics, init_spherical_harmonics
|
||||
|
||||
CONTAINS
|
||||
|
||||
|
|
@ -307,6 +306,274 @@ CONTAINS
|
|||
|
||||
END SUBROUTINE init_spherical_harmonics
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Calculate rotation matrices for spherical harmonics up to value l
|
||||
!> Joseph Ivanic and Klaus Ruedenberg
|
||||
!> Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion
|
||||
!> J. Phys. Chem. 1996, 100, 6342-6347
|
||||
!> \param orbrotmat ...
|
||||
!> \param rotmat ...
|
||||
!> \param lval ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE calculate_rotmat(orbrotmat, rotmat, lval)
|
||||
TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrotmat
|
||||
REAL(KIND=dp), DIMENSION(3, 3) :: rotmat
|
||||
INTEGER, INTENT(IN) :: lval
|
||||
|
||||
INTEGER :: l, m1, m2, ns
|
||||
REAL(KIND=dp) :: s3, u, uf, v, vf, w, wf
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r
|
||||
REAL(KIND=dp), DIMENSION(-1:1, -1:1) :: r1
|
||||
REAL(KIND=dp), DIMENSION(-2:2, -2:2) :: r2
|
||||
REAL(KIND=dp), DIMENSION(3, 3) :: t
|
||||
|
||||
MARK_USED(rotmat)
|
||||
|
||||
CALL release_rotmat(orbrotmat)
|
||||
|
||||
ALLOCATE (orbrotmat(0:lval))
|
||||
DO l = 0, lval
|
||||
ns = nso(l)
|
||||
ALLOCATE (orbrotmat(l)%mat(ns, ns))
|
||||
END DO
|
||||
|
||||
IF (lval >= 0) THEN
|
||||
orbrotmat(0)%mat = 1.0_dp
|
||||
END IF
|
||||
IF (lval >= 1) THEN
|
||||
t(1, 1:3) = rotmat(2, 1:3)
|
||||
t(2, 1:3) = rotmat(3, 1:3)
|
||||
t(3, 1:3) = rotmat(1, 1:3)
|
||||
r1(-1:1, -1) = t(1:3, 2)
|
||||
r1(-1:1, 0) = t(1:3, 3)
|
||||
r1(-1:1, 1) = t(1:3, 1)
|
||||
orbrotmat(1)%mat(1:3, 1:3) = r1(-1:1, -1:1)
|
||||
END IF
|
||||
IF (lval >= 2) THEN
|
||||
s3 = SQRT(3.0_dp)
|
||||
! Table 4
|
||||
r2(0, 0) = (3.0_dp*r1(0, 0)**2 - 1.0_dp)*0.5_dp
|
||||
r2(0, 1) = s3*r1(0, 1)*r1(0, 0)
|
||||
r2(0, -1) = s3*r1(0, -1)*r1(0, 0)
|
||||
r2(0, -2) = s3*(r1(0, 1)**2 - r1(0, -1)**2)*0.5_dp
|
||||
r2(0, -2) = s3*r1(0, 1)*r1(0, -1)
|
||||
!
|
||||
r2(1, 0) = s3*r1(1, 0)*r1(0, 0)
|
||||
r2(1, 1) = r1(1, 1)*r1(0, 0) + r1(1, 0)*r1(0, 1)
|
||||
r2(1, -1) = r1(1, -1)*r1(0, 0) + r1(1, 0)*r1(0, -1)
|
||||
r2(1, 2) = r1(1, 1)*r1(0, 1) - r1(1, -1)*r1(0, -1)
|
||||
r2(1, -2) = r1(1, 1)*r1(0, -1) + r1(1, -1)*r1(0, 1)
|
||||
!
|
||||
r2(-1, 0) = s3*r1(-1, 0)*r1(0, 0)
|
||||
r2(-1, 1) = r1(-1, 1)*r1(0, 0) + r1(-1, 0)*r1(0, 1)
|
||||
r2(-1, -1) = r1(-1, -1)*r1(0, 0) + r1(-1, 0)*r1(0, -1)
|
||||
r2(-1, 2) = r1(-1, 1)*r1(0, 1) - r1(-1, -1)*r1(0, -1)
|
||||
r2(-1, -2) = r1(-1, 1)*r1(0, -1) + r1(-1, -1)*r1(0, 1)
|
||||
!
|
||||
r2(2, 0) = s3*(r1(1, 0)**2 - r1(-1, 0)**2)*0.5_dp
|
||||
r2(2, 1) = r1(1, 1)*r1(1, 0) - r1(-1, 1)*r1(-1, 0)
|
||||
r2(2, -1) = r1(1, -1)*r1(1, 0) - r1(-1, -1)*r1(-1, 0)
|
||||
r2(2, 2) = (r1(1, 1)**2 - r1(1, -1)**2 - r1(-1, 1)**2 + r1(-1, -1)**2)*0.5_dp
|
||||
r2(2, -2) = r1(1, 1)*r1(1, -1) - r1(-1, 1)*r1(-1, -1)
|
||||
!
|
||||
r2(-2, 0) = s3*r1(1, 0)*r1(-1, 0)
|
||||
r2(2, 1) = r1(1, 1)*r1(-1, 0) + r1(1, 0)*r1(-1, 1)
|
||||
r2(2, -1) = r1(1, -1)*r1(-1, 0) + r1(1, 0)*r1(-1, -1)
|
||||
r2(2, 2) = r1(1, 1)*r1(-1, 1) - r1(1, -1)*r1(-1, -1)
|
||||
r2(2, -2) = r1(1, 1)*r1(-1, -1) + r1(1, -1)*r1(-1, 1)
|
||||
!
|
||||
orbrotmat(2)%mat(1:5, 1:5) = r2(-2:2, -2:2)
|
||||
END IF
|
||||
IF (lval >= 3) THEN
|
||||
! General recursion
|
||||
ALLOCATE (r(0:lval, -lval:lval, -lval:lval))
|
||||
r = 0.0_dp
|
||||
r(0, 0, 0) = 1.0_dp
|
||||
r(1, -1:1, -1:1) = r1(-1:1, -1:1)
|
||||
r(2, -2:2, -2:2) = r2(-2:2, -2:2)
|
||||
DO l = 3, lval
|
||||
DO m1 = -l, l
|
||||
DO m2 = -l, l
|
||||
u = u_func(l, m1, m2)
|
||||
v = v_func(l, m1, m2)
|
||||
w = w_func(l, m1, m2)
|
||||
CALL r_recursion(l, m1, m2, r1, r, lval, uf, vf, wf)
|
||||
r(l, m1, m2) = u*uf + v*vf + w*wf
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
DO l = 3, lval
|
||||
ns = nso(l)
|
||||
orbrotmat(l)%mat(1:ns, 1:ns) = r(l, -l:l, -l:l)
|
||||
END DO
|
||||
DEALLOCATE (r)
|
||||
END IF
|
||||
|
||||
END SUBROUTINE calculate_rotmat
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param l ...
|
||||
!> \param ma ...
|
||||
!> \param mb ...
|
||||
!> \return ...
|
||||
! **************************************************************************************************
|
||||
FUNCTION u_func(l, ma, mb) RESULT(u)
|
||||
INTEGER :: l, ma, mb
|
||||
REAL(KIND=dp) :: u
|
||||
|
||||
IF (ABS(mb) == l) THEN
|
||||
u = REAL((l + ma)*(l - ma), KIND=dp)/REAL(2*l*(2*l - 1), KIND=dp)
|
||||
u = SQRT(u)
|
||||
ELSE IF (ABS(mb) < l) THEN
|
||||
u = REAL((l + ma)*(l - ma), KIND=dp)/REAL((l + mb)*(l - mb), KIND=dp)
|
||||
u = SQRT(u)
|
||||
ELSE
|
||||
CPABORT("Illegal m value")
|
||||
END IF
|
||||
END FUNCTION u_func
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param l ...
|
||||
!> \param ma ...
|
||||
!> \param mb ...
|
||||
!> \return ...
|
||||
! **************************************************************************************************
|
||||
FUNCTION v_func(l, ma, mb) RESULT(v)
|
||||
INTEGER :: l, ma, mb
|
||||
REAL(KIND=dp) :: v
|
||||
|
||||
INTEGER :: a1, a2, dm0
|
||||
|
||||
dm0 = 0
|
||||
IF (ma == 0) dm0 = 1
|
||||
IF (ABS(mb) == l) THEN
|
||||
a1 = (1 + dm0)*(l + ABS(ma) - 1)*(l + ABS(ma))
|
||||
a2 = 2*l*(2*l - 1)
|
||||
ELSE IF (ABS(mb) < l) THEN
|
||||
a1 = (1 + dm0)*(l + ABS(ma) - 1)*(l + ABS(ma))
|
||||
a2 = (l + mb)*(l - mb)
|
||||
ELSE
|
||||
CPABORT("Illegal m value")
|
||||
END IF
|
||||
v = 0.5_dp*SQRT(REAL(a1, KIND=dp)/REAL(a2, KIND=dp))*REAL(1 - 2*dm0, KIND=dp)
|
||||
END FUNCTION v_func
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param l ...
|
||||
!> \param ma ...
|
||||
!> \param mb ...
|
||||
!> \return ...
|
||||
! **************************************************************************************************
|
||||
FUNCTION w_func(l, ma, mb) RESULT(w)
|
||||
INTEGER :: l, ma, mb
|
||||
REAL(KIND=dp) :: w
|
||||
|
||||
INTEGER :: a1, a2, dm0
|
||||
|
||||
dm0 = 0
|
||||
IF (ma == 0) dm0 = 1
|
||||
IF (ABS(mb) == l) THEN
|
||||
a1 = (l - ABS(ma) - 1)*(l - ABS(ma))
|
||||
a2 = 2*l*(2*l - 1)
|
||||
ELSE IF (ABS(mb) < l) THEN
|
||||
a1 = (l - ABS(ma) - 1)*(l - ABS(ma))
|
||||
a2 = (l + mb)*(l - mb)
|
||||
ELSE
|
||||
CPABORT("Illegal m value")
|
||||
END IF
|
||||
w = -0.5_dp*SQRT(REAL(a1, KIND=dp)/REAL(a2, KIND=dp))*REAL(1 - dm0, KIND=dp)
|
||||
END FUNCTION w_func
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param i ...
|
||||
!> \param l ...
|
||||
!> \param mu ...
|
||||
!> \param m ...
|
||||
!> \param r1 ...
|
||||
!> \param r ...
|
||||
!> \param lr ...
|
||||
!> \return ...
|
||||
! **************************************************************************************************
|
||||
FUNCTION p_func(i, l, mu, m, r1, r, lr) RESULT(p)
|
||||
INTEGER :: i, l, mu, m
|
||||
REAL(KIND=dp), DIMENSION(-1:1, -1:1) :: r1
|
||||
INTEGER :: lr
|
||||
REAL(KIND=dp), DIMENSION(0:lr, -lr:lr, -lr:lr) :: r
|
||||
REAL(KIND=dp) :: p
|
||||
|
||||
IF (m == l) THEN
|
||||
p = r1(i, 1)*r(l - 1, mu, m) - r1(i, -1)*r(l - 1, mu, -m)
|
||||
ELSE IF (m == -l) THEN
|
||||
p = r1(i, 1)*r(l - 1, mu, m) + r1(i, -1)*r(l - 1, mu, -m)
|
||||
ELSE IF (ABS(m) < l) THEN
|
||||
p = r1(i, 0)*r(l - 1, mu, m)
|
||||
ELSE
|
||||
CPABORT("Illegal m value")
|
||||
END IF
|
||||
END FUNCTION p_func
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param l ...
|
||||
!> \param ma ...
|
||||
!> \param mb ...
|
||||
!> \param r1 ...
|
||||
!> \param r ...
|
||||
!> \param lr ...
|
||||
!> \param u ...
|
||||
!> \param v ...
|
||||
!> \param w ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE r_recursion(l, ma, mb, r1, r, lr, u, v, w)
|
||||
INTEGER :: l, ma, mb
|
||||
REAL(KIND=dp), DIMENSION(-1:1, -1:1) :: r1
|
||||
INTEGER :: lr
|
||||
REAL(KIND=dp), DIMENSION(0:lr, -lr:lr, -lr:lr) :: r
|
||||
REAL(KIND=dp) :: u, v, w
|
||||
|
||||
REAL(KIND=dp) :: dd
|
||||
|
||||
IF (ma == 0) THEN
|
||||
u = p_func(0, l, 0, mb, r1, r, lr)
|
||||
v = p_func(1, l, 1, mb, r1, r, lr) + p_func(-1, l, -1, mb, r1, r, lr)
|
||||
w = 0.0_dp
|
||||
ELSE IF (ma > 0) THEN
|
||||
dd = 1.0_dp
|
||||
IF (ma == 1) dd = 1.0_dp
|
||||
u = p_func(0, l, ma, mb, r1, r, lr)
|
||||
v = p_func(1, l, ma - 1, mb, r1, r, lr)*SQRT(1.0_dp + dd) - p_func(-1, l, -ma + 1, mb, r1, r, lr)*(1.0_dp - dd)
|
||||
w = p_func(1, l, ma + 1, mb, r1, r, lr) + p_func(-1, l, -1 - ma, mb, r1, r, lr)
|
||||
ELSE IF (ma < 0) THEN
|
||||
dd = 1.0_dp
|
||||
IF (ma == -1) dd = 1.0_dp
|
||||
u = p_func(0, l, ma, mb, r1, r, lr)
|
||||
v = p_func(1, l, ma + 1, mb, r1, r, lr)*(1.0_dp - dd) + p_func(-1, l, -ma - 1, mb, r1, r, lr)*SQRT(1.0_dp + dd)
|
||||
w = p_func(1, l, ma - 1, mb, r1, r, lr) - p_func(-1, l, -ma + 1, mb, r1, r, lr)
|
||||
END IF
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Release rotation matrices
|
||||
!> \param orbrotmat ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE release_rotmat(orbrotmat)
|
||||
TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrotmat
|
||||
|
||||
INTEGER :: i
|
||||
|
||||
IF (ASSOCIATED(orbrotmat)) THEN
|
||||
DO i = LBOUND(orbrotmat, 1), UBOUND(orbrotmat, 1)
|
||||
IF (ASSOCIATED(orbrotmat(i)%mat)) DEALLOCATE (orbrotmat(i)%mat)
|
||||
END DO
|
||||
DEALLOCATE (orbrotmat)
|
||||
END IF
|
||||
|
||||
END SUBROUTINE release_rotmat
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Print a matrix to the logical unit "lunit".
|
||||
!> \param matrix ...
|
||||
|
|
|
|||
252
src/cryssym.F
252
src/cryssym.F
|
|
@ -14,17 +14,23 @@ MODULE cryssym
|
|||
USE bibliography, ONLY: Togo2018,&
|
||||
cite_reference
|
||||
USE kinds, ONLY: dp
|
||||
USE spglib_f08, ONLY: &
|
||||
spg_get_international, spg_get_ir_reciprocal_mesh, spg_get_major_version, &
|
||||
spg_get_micro_version, spg_get_minor_version, spg_get_multiplicity, spg_get_pointgroup, &
|
||||
spg_get_schoenflies, spg_get_symmetry
|
||||
USE kpsym, ONLY: group1s,&
|
||||
k290s
|
||||
USE spglib_f08, ONLY: spg_get_international,&
|
||||
spg_get_major_version,&
|
||||
spg_get_micro_version,&
|
||||
spg_get_minor_version,&
|
||||
spg_get_multiplicity,&
|
||||
spg_get_pointgroup,&
|
||||
spg_get_schoenflies,&
|
||||
spg_get_symmetry
|
||||
USE string_utilities, ONLY: strip_control_codes
|
||||
#include "./base/base_uses.f90"
|
||||
|
||||
IMPLICIT NONE
|
||||
PRIVATE
|
||||
PUBLIC :: csym_type, release_csym_type, print_crys_symmetry, print_kp_symmetry
|
||||
PUBLIC :: crys_sym_gen, kpoint_gen, apply_rotation_coord
|
||||
PUBLIC :: crys_sym_gen, kpoint_gen
|
||||
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cryssym'
|
||||
|
||||
|
|
@ -60,6 +66,12 @@ MODULE cryssym
|
|||
INTEGER :: n_operations = 0
|
||||
INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: rotations
|
||||
REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: translations
|
||||
!K290
|
||||
REAL(KIND=dp), DIMENSION(3, 3, 48) :: rt = 0.0_dp
|
||||
REAL(KIND=dp), DIMENSION(3, 48) :: vt = 0.0_dp
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0
|
||||
INTEGER :: nrtot = 0
|
||||
INTEGER, DIMENSION(48) :: ibrot = 1
|
||||
END TYPE csym_type
|
||||
|
||||
CONTAINS
|
||||
|
|
@ -98,6 +110,9 @@ CONTAINS
|
|||
IF (ALLOCATED(csym%kpop)) THEN
|
||||
DEALLOCATE (csym%kpop)
|
||||
END IF
|
||||
IF (ALLOCATED(csym%f0)) THEN
|
||||
DEALLOCATE (csym%f0)
|
||||
END IF
|
||||
|
||||
END SUBROUTINE release_csym_type
|
||||
|
||||
|
|
@ -214,14 +229,11 @@ CONTAINS
|
|||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_gen'
|
||||
|
||||
INTEGER :: handle, i, ik, is_shift(3), &
|
||||
is_time_reversal, j, nkp, nkpts
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:) :: iwkp, xptr
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ixkp
|
||||
INTEGER :: handle, i, ik, j, nkp, nkpts
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:) :: kpop, xptr
|
||||
LOGICAL :: fullmesh
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wkp
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp
|
||||
REAL(KIND=dp), DIMENSION(3) :: rxkp
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
|
|
@ -247,7 +259,7 @@ CONTAINS
|
|||
csym%mesh(1:3) = nk(1:3)
|
||||
|
||||
nkpts = nk(1)*nk(2)*nk(3)
|
||||
ALLOCATE (xkp(3, nkpts), wkp(nkpts))
|
||||
ALLOCATE (xkp(3, nkpts), wkp(nkpts), kpop(nkpts))
|
||||
! kp: link
|
||||
ALLOCATE (csym%kplink(2, nkpts))
|
||||
csym%kplink = 0
|
||||
|
|
@ -274,20 +286,10 @@ CONTAINS
|
|||
CALL cp_abort(__LOCATION__, "MacDonald shifted k-point meshes are only "// &
|
||||
"possible without symmetrization.")
|
||||
END IF
|
||||
ALLOCATE (ixkp(3, nkpts), iwkp(nkpts))
|
||||
is_shift(1:3) = MOD(nk(1:3) + 1, 2)
|
||||
is_time_reversal = 1
|
||||
nkp = spg_get_ir_reciprocal_mesh(ixkp, iwkp, nk, is_shift, is_time_reversal, &
|
||||
TRANSPOSE(csym%hmat), csym%scoord, csym%atype, &
|
||||
csym%nat, csym%delta)
|
||||
wkp = 0.0_dp
|
||||
DO i = 1, nkpts
|
||||
xkp(:, i) = REAL(is_shift + 2*ixkp(:, i), KIND=dp)/REAL(2*nk(:), KIND=dp)
|
||||
j = iwkp(i) + 1
|
||||
wkp(j) = wkp(j) + 1.0_dp
|
||||
END DO
|
||||
csym%kplink(1, 1:nkpts) = iwkp(1:nkpts) + 1
|
||||
DEALLOCATE (ixkp, iwkp)
|
||||
|
||||
CALL full_grid_gen(nk, xkp, wkp, shift)
|
||||
CALL kp_symmetry(csym, xkp, wkp, kpop)
|
||||
|
||||
END IF
|
||||
ELSE
|
||||
! no symmetry library is available
|
||||
|
|
@ -343,24 +345,9 @@ CONTAINS
|
|||
ALLOCATE (csym%kpop(nkpts))
|
||||
IF (csym%symlib .AND. csym%istriz == 1 .AND. .NOT. fullmesh) THEN
|
||||
! atomic symmetry operations possible
|
||||
csym%kpop(1:nkpts) = kpop(1:nkpts)
|
||||
DO ik = 1, nkpts
|
||||
IF (wkp(ik) > 0.0_dp) THEN
|
||||
csym%kpop(ik) = 1
|
||||
ELSE
|
||||
csym%kpop(ik) = 0
|
||||
j = csym%kplink(2, ik)
|
||||
DO i = 1, csym%n_operations
|
||||
IF (SUM(ABS(csym%translations(:, i))) > 1.e-10_dp) CYCLE
|
||||
rxkp(1:3) = kp_apply_operation(csym%xkpoint(1:3, j), csym%rotations(:, :, i))
|
||||
rxkp(1:3) = ABS(xkp(1:3, ik) - rxkp(1:3))
|
||||
rxkp(1:3) = rxkp(1:3) - REAL(NINT(rxkp(1:3)), KIND=dp)
|
||||
IF (ALL((rxkp(1:3)) < 1.e-12_dp)) THEN
|
||||
csym%kpop(ik) = i
|
||||
EXIT
|
||||
END IF
|
||||
END DO
|
||||
CPASSERT(csym%kpop(ik) /= 0)
|
||||
END IF
|
||||
CPASSERT(csym%kpop(ik) /= 0)
|
||||
END DO
|
||||
ELSE
|
||||
! only time reversal symmetry
|
||||
|
|
@ -373,7 +360,7 @@ CONTAINS
|
|||
END DO
|
||||
END IF
|
||||
|
||||
DEALLOCATE (xkp, wkp)
|
||||
DEALLOCATE (xkp, wkp, kpop)
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
|
|
@ -381,6 +368,128 @@ CONTAINS
|
|||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param csym ...
|
||||
!> \param xkp ...
|
||||
!> \param wkp ...
|
||||
!> \param kpop ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE kp_symmetry(csym, xkp, wkp, kpop)
|
||||
TYPE(csym_type) :: csym
|
||||
REAL(KIND=dp), DIMENSION(:, :) :: xkp
|
||||
REAL(KIND=dp), DIMENSION(:) :: wkp
|
||||
INTEGER, DIMENSION(:) :: kpop
|
||||
|
||||
INTEGER :: i, ihc, ihg, ik, indpg, iou, iq1, iq2, &
|
||||
iq3, istriz, isy, itoj, j, kr, li, lr, &
|
||||
nat, nc, nhash, nkm, nkp, nkpoint, &
|
||||
nsp, ntvec
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:) :: includ, isc, list, lwght, ty
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0, lrot
|
||||
INTEGER, DIMENSION(48) :: ib
|
||||
REAL(KIND=dp) :: alat
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rlist, rx, tvec, wvkl, xkapa
|
||||
REAL(KIND=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3, origin, rr, wvk0
|
||||
REAL(KIND=dp), DIMENSION(3, 3) :: hmat, strain
|
||||
REAL(KIND=dp), DIMENSION(3, 3, 48) :: r
|
||||
REAL(KIND=dp), DIMENSION(3, 48) :: vt
|
||||
|
||||
iou = csym%punit
|
||||
hmat = csym%hmat
|
||||
nat = csym%nat
|
||||
iq1 = csym%mesh(1)
|
||||
iq2 = csym%mesh(2)
|
||||
iq3 = csym%mesh(3)
|
||||
nkpoint = 10*iq1*iq2*iq3
|
||||
nkpoint = 2*MAX(iq1, iq2, iq3)**3
|
||||
wvk0 = csym%wvk0
|
||||
istriz = csym%istriz
|
||||
a1(1:3) = hmat(1:3, 1)
|
||||
a2(1:3) = hmat(1:3, 2)
|
||||
a3(1:3) = hmat(1:3, 3)
|
||||
alat = hmat(1, 1)
|
||||
strain = 0.0_dp
|
||||
ALLOCATE (xkapa(3, nat), rx(3, nat), tvec(3, 200), ty(nat), isc(nat), f0(49, nat))
|
||||
ty(1:nat) = csym%atype(1:nat)
|
||||
nsp = MAXVAL(ty)
|
||||
DO i = 1, nat
|
||||
xkapa(1:3, i) = MATMUL(hmat, csym%scoord(1:3, i))
|
||||
END DO
|
||||
nhash = 1000
|
||||
ALLOCATE (wvkl(3, nkpoint), rlist(3, nkpoint), includ(nkpoint), list(nhash + nkpoint))
|
||||
ALLOCATE (lrot(48, nkpoint), lwght(nkpoint))
|
||||
|
||||
IF (iou > 0) THEN
|
||||
WRITE (iou, '(/,(T2,A79))') &
|
||||
"*******************************************************************************", &
|
||||
"** Special K-Point Generation by K290 **", &
|
||||
"*******************************************************************************"
|
||||
END IF
|
||||
|
||||
CALL K290s(iou, nat, nkpoint, nsp, iq1, iq2, iq3, istriz, &
|
||||
a1, a2, a3, alat, strain, xkapa, rx, tvec, &
|
||||
ty, isc, f0, ntvec, wvk0, wvkl, lwght, lrot, &
|
||||
nhash, includ, list, rlist, csym%delta)
|
||||
|
||||
CALL GROUP1s(0, a1, a2, a3, nat, ty, xkapa, b1, b2, b3, &
|
||||
ihg, ihc, isy, li, nc, indpg, ib, ntvec, &
|
||||
vt, f0, r, tvec, origin, rx, isc, csym%delta)
|
||||
|
||||
IF (iou > 0) THEN
|
||||
WRITE (iou, '((T2,A79))') &
|
||||
"*******************************************************************************", &
|
||||
"** Finished K290 **", &
|
||||
"*******************************************************************************"
|
||||
END IF
|
||||
|
||||
csym%rt = r
|
||||
csym%vt = vt
|
||||
csym%nrtot = nc
|
||||
ALLOCATE (csym%f0(nat, nc))
|
||||
DO i = 1, nc
|
||||
csym%f0(1:nat, i) = f0(i, 1:nat)
|
||||
END DO
|
||||
csym%ibrot = 0
|
||||
csym%ibrot(1:nc) = ib(1:nc)
|
||||
|
||||
kpop = 0
|
||||
nkm = iq1*iq2*iq3
|
||||
nkp = 0
|
||||
DO i = 1, nkm
|
||||
IF (lwght(i) == 0) EXIT
|
||||
nkp = nkp + 1
|
||||
END DO
|
||||
wkp = 0
|
||||
ik = 0
|
||||
DO i = 1, nkp
|
||||
DO j = 1, nkm
|
||||
wvk0(1:3) = xkp(1:3, j) - wvkl(1:3, i)
|
||||
IF (ALL(ABS(wvk0(1:3)) < 1.e-12_dp)) THEN
|
||||
wkp(j) = lwght(i)
|
||||
itoj = j
|
||||
EXIT
|
||||
END IF
|
||||
END DO
|
||||
DO lr = 1, lwght(i)
|
||||
kr = lrot(lr, i)
|
||||
rr(1:3) = kp_apply_operation(wvkl(1:3, i), r(1:3, 1:3, ABS(kr)))
|
||||
IF (kr < 0) rr(1:3) = -rr(1:3)
|
||||
DO j = 1, nkm
|
||||
wvk0(1:3) = xkp(1:3, j) - rr(1:3)
|
||||
IF (ALL(ABS(wvk0(1:3)) < 1.e-12_dp)) THEN
|
||||
csym%kplink(1, j) = itoj
|
||||
kpop(j) = kr
|
||||
EXIT
|
||||
END IF
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
DEALLOCATE (xkapa, rx, tvec, ty, isc, f0)
|
||||
DEALLOCATE (wvkl, rlist, includ, list)
|
||||
DEALLOCATE (lrot, lwght)
|
||||
|
||||
END SUBROUTINE kp_symmetry
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param nk ...
|
||||
!> \param xkp ...
|
||||
!> \param wkp ...
|
||||
|
|
@ -453,55 +562,15 @@ CONTAINS
|
|||
!> \return ...
|
||||
! **************************************************************************************************
|
||||
FUNCTION kp_apply_operation(x, r) RESULT(y)
|
||||
REAL(KIND=dp), INTENT(IN) :: x(3)
|
||||
INTEGER, INTENT(IN) :: r(3, 3)
|
||||
REAL(KIND=dp), INTENT(IN) :: x(3), r(3, 3)
|
||||
REAL(KIND=dp) :: y(3)
|
||||
|
||||
y(1) = REAL(r(1, 1), dp)*x(1) + REAL(r(1, 2), dp)*x(2) + REAL(r(1, 3), dp)*x(3)
|
||||
y(2) = REAL(r(2, 1), dp)*x(1) + REAL(r(2, 2), dp)*x(2) + REAL(r(2, 3), dp)*x(3)
|
||||
y(3) = REAL(r(3, 1), dp)*x(1) + REAL(r(3, 2), dp)*x(2) + REAL(r(3, 3), dp)*x(3)
|
||||
y(1) = r(1, 1)*x(1) + r(1, 2)*x(2) + r(1, 3)*x(3)
|
||||
y(2) = r(2, 1)*x(1) + r(2, 2)*x(2) + r(2, 3)*x(3)
|
||||
y(3) = r(3, 1)*x(1) + r(3, 2)*x(2) + r(3, 3)*x(3)
|
||||
|
||||
END FUNCTION kp_apply_operation
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param f0 ...
|
||||
!> \param csym ...
|
||||
!> \param ir ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE apply_rotation_coord(f0, csym, ir)
|
||||
INTEGER, DIMENSION(:), INTENT(INOUT) :: f0
|
||||
TYPE(csym_type) :: csym
|
||||
INTEGER, INTENT(IN) :: ir
|
||||
|
||||
INTEGER :: ia, ib, natom
|
||||
REAL(KIND=dp) :: diff
|
||||
REAL(KIND=dp), DIMENSION(3) :: rb, ri, ro, tr
|
||||
REAL(KIND=dp), DIMENSION(3, 3) :: rot
|
||||
|
||||
natom = csym%nat
|
||||
rot(1:3, 1:3) = csym%rotations(1:3, 1:3, ir)
|
||||
tr(1:3) = csym%translations(1:3, ir)
|
||||
|
||||
f0 = 0
|
||||
DO ia = 1, natom
|
||||
ri(1:3) = csym%scoord(1:3, ia)
|
||||
ro(1) = REAL(rot(1, 1), dp)*ri(1) + REAL(rot(2, 1), dp)*ri(2) + REAL(rot(3, 1), dp)*ri(3) + tr(1)
|
||||
ro(2) = REAL(rot(1, 2), dp)*ri(1) + REAL(rot(2, 2), dp)*ri(2) + REAL(rot(3, 2), dp)*ri(3) + tr(2)
|
||||
ro(3) = REAL(rot(1, 3), dp)*ri(1) + REAL(rot(2, 3), dp)*ri(2) + REAL(rot(3, 3), dp)*ri(3) + tr(3)
|
||||
DO ib = 1, natom
|
||||
rb(1:3) = csym%scoord(1:3, ib)
|
||||
diff = SQRT(SUM((ri(:) - rb(:))**2))
|
||||
IF (diff < csym%delta) THEN
|
||||
f0(ia) = ib
|
||||
EXIT
|
||||
END IF
|
||||
END DO
|
||||
CPASSERT(f0(ia) /= 0)
|
||||
END DO
|
||||
|
||||
END SUBROUTINE apply_rotation_coord
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param csym ...
|
||||
|
|
@ -542,7 +611,7 @@ CONTAINS
|
|||
SUBROUTINE print_kp_symmetry(csym)
|
||||
TYPE(csym_type), INTENT(IN) :: csym
|
||||
|
||||
INTEGER :: i, iunit, plevel
|
||||
INTEGER :: i, iunit, nat, plevel
|
||||
|
||||
iunit = csym%punit
|
||||
IF (iunit >= 0) THEN
|
||||
|
|
@ -559,6 +628,13 @@ CONTAINS
|
|||
WRITE (iunit, '(T2,i10,3F10.5,T45,3I12)') i, csym%kpmesh(1:3, i), &
|
||||
csym%kplink(1:2, i), csym%kpop(i)
|
||||
END DO
|
||||
IF (csym%nrtot > 0) THEN
|
||||
WRITE (iunit, '(/,A)') " Atom Transformation Table"
|
||||
nat = SIZE(csym%f0, 1)
|
||||
DO i = 1, csym%nrtot
|
||||
WRITE (iunit, '(T10,A,I3,(T21,12I5))') " Rot=", csym%ibrot(i), csym%f0(1:nat, i)
|
||||
END DO
|
||||
END IF
|
||||
END IF
|
||||
|
||||
END SUBROUTINE print_kp_symmetry
|
||||
|
|
|
|||
|
|
@ -36,8 +36,7 @@ MODULE kpoint_methods
|
|||
copy_info_type, cp_fm_cleanup_copy_general, cp_fm_create, cp_fm_finish_copy_general, &
|
||||
cp_fm_get_info, cp_fm_release, cp_fm_start_copy_general, cp_fm_to_fm, cp_fm_type
|
||||
USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
|
||||
USE cryssym, ONLY: apply_rotation_coord,&
|
||||
crys_sym_gen,&
|
||||
USE cryssym, ONLY: crys_sym_gen,&
|
||||
csym_type,&
|
||||
kpoint_gen,&
|
||||
print_crys_symmetry,&
|
||||
|
|
@ -48,6 +47,7 @@ MODULE kpoint_methods
|
|||
USE input_constants, ONLY: smear_fermi_dirac
|
||||
USE kinds, ONLY: dp
|
||||
USE kpoint_types, ONLY: get_kpoint_info,&
|
||||
kind_rotmat_type,&
|
||||
kpoint_env_create,&
|
||||
kpoint_env_p_type,&
|
||||
kpoint_env_type,&
|
||||
|
|
@ -109,8 +109,8 @@ CONTAINS
|
|||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize'
|
||||
|
||||
INTEGER :: handle, i, ik, iounit, ir, is, natom, &
|
||||
nr, ns
|
||||
INTEGER :: handle, i, ic, ik, iounit, ir, ira, is, &
|
||||
j, natom, nkind, nr, ns
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:) :: atype
|
||||
LOGICAL :: spez
|
||||
REAL(KIND=dp) :: wsum
|
||||
|
|
@ -159,6 +159,9 @@ CONTAINS
|
|||
ELSE
|
||||
iounit = -1
|
||||
END IF
|
||||
! kind type list
|
||||
ALLOCATE (kpoint%atype(natom))
|
||||
kpoint%atype = atype
|
||||
|
||||
CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=kpoint%eps_geo, iounit=iounit)
|
||||
CALL kpoint_gen(crys_sym, kpoint%nkp_grid, symm=kpoint%symmetry, shift=kpoint%kp_shift, &
|
||||
|
|
@ -185,26 +188,50 @@ CONTAINS
|
|||
! set up the symmetrization information
|
||||
kpsym%nwght = NINT(crys_sym%wkpoint(ik))
|
||||
ns = kpsym%nwght
|
||||
! to be done correctly
|
||||
!
|
||||
IF (ns > 1) THEN
|
||||
kpsym%apply_symmetry = .TRUE.
|
||||
natom = SIZE(particle_set)
|
||||
ALLOCATE (kpsym%rot(3, 3, ns))
|
||||
ALLOCATE (kpsym%xkp(3, ns))
|
||||
ALLOCATE (kpsym%rotp(ns))
|
||||
ALLOCATE (kpsym%f0(natom, ns))
|
||||
nr = 0
|
||||
DO is = 1, SIZE(crys_sym%kplink, 2)
|
||||
IF (crys_sym%kplink(2, is) == ik) THEN
|
||||
nr = nr + 1
|
||||
ir = crys_sym%kpop(is)
|
||||
kpsym%rot(1:3, 1:3, nr) = crys_sym%rotations(1:3, 1:3, ir)
|
||||
ira = ABS(ir)
|
||||
kpsym%rotp(nr) = ir
|
||||
IF (ir > 0) THEN
|
||||
kpsym%rot(1:3, 1:3, nr) = crys_sym%rt(1:3, 1:3, ira)
|
||||
ELSE
|
||||
kpsym%rot(1:3, 1:3, nr) = -crys_sym%rt(1:3, 1:3, ira)
|
||||
END IF
|
||||
kpsym%xkp(1:3, nr) = crys_sym%kpmesh(1:3, is)
|
||||
CALL apply_rotation_coord(kpsym%f0(1:natom, nr), crys_sym, ir)
|
||||
DO ic = 1, crys_sym%nrtot
|
||||
IF (crys_sym%ibrot(ic) == ira) THEN
|
||||
kpsym%f0(1:natom, nr) = crys_sym%f0(1:natom, ic)
|
||||
EXIT
|
||||
END IF
|
||||
END DO
|
||||
END IF
|
||||
END DO
|
||||
! Reduce inversion?
|
||||
kpsym%nwred = nr
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
IF (kpoint%symmetry) THEN
|
||||
nkind = MAXVAL(atype)
|
||||
ns = crys_sym%nrtot
|
||||
ALLOCATE (kpoint%kind_rotmat(ns, nkind))
|
||||
DO i = 1, ns
|
||||
DO j = 1, nkind
|
||||
NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
|
||||
END DO
|
||||
END DO
|
||||
END IF
|
||||
|
||||
CALL release_csym_type(crys_sym)
|
||||
DEALLOCATE (scoord, atype)
|
||||
|
|
@ -1125,8 +1152,8 @@ CONTAINS
|
|||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_transform'
|
||||
|
||||
INTEGER :: handle, ic, ik, ikk, indx, is, ispin, &
|
||||
nc, nimg, nkp, nspin
|
||||
INTEGER :: handle, ic, ik, ikk, indx, ir, is, &
|
||||
ispin, nc, nimg, nkp, nspin
|
||||
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
|
||||
LOGICAL :: aux_fit, do_ext, do_symmetric, my_kpgrp, &
|
||||
real_only
|
||||
|
|
@ -1136,6 +1163,7 @@ CONTAINS
|
|||
TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:) :: info
|
||||
TYPE(cp_fm_type) :: fmdummy
|
||||
TYPE(dbcsr_type), POINTER :: cpmat, rpmat, scpmat, srpmat
|
||||
TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
|
||||
TYPE(kpoint_env_type), POINTER :: kp
|
||||
TYPE(kpoint_sym_type), POINTER :: kpsym
|
||||
TYPE(mp_para_env_type), POINTER :: para_env
|
||||
|
|
@ -1162,17 +1190,21 @@ CONTAINS
|
|||
CALL dbcsr_create(rpmat, template=tempmat, &
|
||||
matrix_type=MERGE(dbcsr_type_symmetric, dbcsr_type_no_symmetry, do_symmetric))
|
||||
CALL cp_dbcsr_alloc_block_from_nbl(rpmat, sab_nl)
|
||||
CALL dbcsr_set(rpmat, 0.0_dp)
|
||||
ALLOCATE (cpmat)
|
||||
CALL dbcsr_create(cpmat, template=tempmat, &
|
||||
matrix_type=MERGE(dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, do_symmetric))
|
||||
CALL cp_dbcsr_alloc_block_from_nbl(cpmat, sab_nl)
|
||||
CALL dbcsr_set(cpmat, 0.0_dp)
|
||||
IF (.NOT. kpoint%full_grid) THEN
|
||||
ALLOCATE (srpmat)
|
||||
CALL dbcsr_create(srpmat, template=rpmat)
|
||||
CALL cp_dbcsr_alloc_block_from_nbl(srpmat, sab_nl)
|
||||
CALL dbcsr_set(srpmat, 0.0_dp)
|
||||
ALLOCATE (scpmat)
|
||||
CALL dbcsr_create(scpmat, template=cpmat)
|
||||
CALL cp_dbcsr_alloc_block_from_nbl(scpmat, sab_nl)
|
||||
CALL dbcsr_set(scpmat, 0.0_dp)
|
||||
END IF
|
||||
|
||||
CALL get_kpoint_info(kpoint, nkp=nkp, xkp=xkp, wkp=wkp, &
|
||||
|
|
@ -1259,11 +1291,16 @@ CONTAINS
|
|||
IF (kpsym%apply_symmetry) THEN
|
||||
wkpx = wkp(ik)/REAL(kpsym%nwght, KIND=dp)
|
||||
DO is = 1, kpsym%nwght
|
||||
ir = ABS(kpsym%rotp(is))
|
||||
kind_rot => kpoint%kind_rotmat(ir, :)
|
||||
IF (real_only) THEN
|
||||
CALL symtrans(srpmat, rpmat, kpsym%rot(1:3, 1:3, is), kpsym%f0(:, is), symmetric=.TRUE.)
|
||||
CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
|
||||
kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
|
||||
ELSE
|
||||
CALL symtrans(srpmat, rpmat, kpsym%rot(1:3, 1:3, is), kpsym%f0(:, is), symmetric=.TRUE.)
|
||||
CALL symtrans(scpmat, cpmat, kpsym%rot(1:3, 1:3, is), kpsym%f0(:, is), antisymmetric=.TRUE.)
|
||||
CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
|
||||
kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
|
||||
CALL symtrans(scpmat, cpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
|
||||
kpsym%f0(:, is), kpoint%atype, antisymmetric=.TRUE.)
|
||||
END IF
|
||||
CALL transform_dmat(denmat, srpmat, scpmat, ispin, real_only, sab_nl, &
|
||||
cell_to_index, kpsym%xkp(1:3, is), wkpx)
|
||||
|
|
@ -1363,7 +1400,7 @@ CONTAINS
|
|||
!We have a FT from KP to real-space: S(R) = sum_k S(k)*exp(-i*k*R), with S(k) a complex number
|
||||
!Therefore, we have: S(R) = sum_k Re(S(k))*cos(k*R) -i^2*Im(S(k))*sin(k*R)
|
||||
! = sum_k Re(S(k))*cos(k*R) + Im(S(k))*sin(k*R)
|
||||
!fc = +- 1 is due to the usual non-symmetric real-sapce matrices stored as symmetric ones
|
||||
!fc = +- 1 is due to the usual non-symmetric real-space matrices stored as symmetric ones
|
||||
|
||||
irow = iatom
|
||||
icol = jatom
|
||||
|
|
@ -1408,24 +1445,27 @@ CONTAINS
|
|||
!> \brief Symmetrization of density matrix - transform to new k-point
|
||||
!> \param smat density matrix at new kpoint
|
||||
!> \param pmat reference density matrix
|
||||
!> \param kmat Kind type rotation matrix
|
||||
!> \param rot Rotation matrix
|
||||
!> \param f0 Permutation of atoms under transformation
|
||||
!> \param atype Atom to kind pointer
|
||||
!> \param symmetric Symmetric matrix
|
||||
!> \param antisymmetric Anti-Symmetric matrix
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE symtrans(smat, pmat, rot, f0, symmetric, antisymmetric)
|
||||
SUBROUTINE symtrans(smat, pmat, kmat, rot, f0, atype, symmetric, antisymmetric)
|
||||
TYPE(dbcsr_type), POINTER :: smat, pmat
|
||||
TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kmat
|
||||
REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: rot
|
||||
INTEGER, DIMENSION(:), INTENT(IN) :: f0
|
||||
INTEGER, DIMENSION(:), INTENT(IN) :: f0, atype
|
||||
LOGICAL, INTENT(IN), OPTIONAL :: symmetric, antisymmetric
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'symtrans'
|
||||
|
||||
INTEGER :: handle, iatom, icol, ip, irow, jcol, jp, &
|
||||
jrow, natom, numnodes
|
||||
INTEGER :: handle, iatom, icol, ikind, ip, irow, &
|
||||
jcol, jkind, jp, jrow, natom, numnodes
|
||||
LOGICAL :: asym, dorot, found, perm, sym, trans
|
||||
REAL(KIND=dp) :: dr, fsign
|
||||
REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, sblock
|
||||
REAL(KIND=dp), DIMENSION(:, :), POINTER :: kroti, krotj, pblock, sblock
|
||||
TYPE(dbcsr_distribution_type) :: dist
|
||||
TYPE(dbcsr_iterator_type) :: iter
|
||||
|
||||
|
|
@ -1459,6 +1499,8 @@ CONTAINS
|
|||
IF (asym) fsign = -1.0_dp
|
||||
|
||||
IF (dorot .OR. perm) THEN
|
||||
CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
|
||||
"Reduced grids not yet working correctly")
|
||||
CALL dbcsr_set(smat, 0.0_dp)
|
||||
IF (perm) THEN
|
||||
CALL dbcsr_get_info(pmat, distribution=dist)
|
||||
|
|
@ -1480,11 +1522,16 @@ CONTAINS
|
|||
trans = .TRUE.
|
||||
END IF
|
||||
CALL dbcsr_get_block_p(matrix=smat, row=jrow, col=jcol, BLOCK=sblock, found=found)
|
||||
IF (.NOT. found) CYCLE
|
||||
CPASSERT(found)
|
||||
ikind = atype(irow)
|
||||
jkind = atype(icol)
|
||||
kroti => kmat(ikind)%rmat
|
||||
krotj => kmat(jkind)%rmat
|
||||
! rotation
|
||||
IF (trans) THEN
|
||||
sblock = fsign*TRANSPOSE(pblock)
|
||||
sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(pblock), kroti))
|
||||
ELSE
|
||||
sblock = pblock
|
||||
sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(pblock, krotj))
|
||||
END IF
|
||||
END DO
|
||||
CALL dbcsr_iterator_stop(iter)
|
||||
|
|
@ -1496,10 +1543,34 @@ CONTAINS
|
|||
END IF
|
||||
ELSE
|
||||
! no atom permutations, this is always local
|
||||
! ignore rotations for now
|
||||
CALL dbcsr_copy(smat, pmat)
|
||||
CALL cp_abort(__LOCATION__, "k-points need FULL_GRID currently. "// &
|
||||
"Reduced grids not yet working correctly")
|
||||
CALL dbcsr_iterator_start(iter, smat)
|
||||
DO WHILE (dbcsr_iterator_blocks_left(iter))
|
||||
CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
|
||||
ip = f0(irow)
|
||||
jp = f0(icol)
|
||||
IF (ip <= jp) THEN
|
||||
jrow = ip
|
||||
jcol = jp
|
||||
trans = .FALSE.
|
||||
ELSE
|
||||
jrow = jp
|
||||
jcol = ip
|
||||
trans = .TRUE.
|
||||
END IF
|
||||
ikind = atype(irow)
|
||||
jkind = atype(icol)
|
||||
kroti => kmat(ikind)%rmat
|
||||
krotj => kmat(jkind)%rmat
|
||||
! rotation
|
||||
IF (trans) THEN
|
||||
sblock = fsign*MATMUL(TRANSPOSE(krotj), MATMUL(TRANSPOSE(sblock), kroti))
|
||||
ELSE
|
||||
sblock = fsign*MATMUL(TRANSPOSE(kroti), MATMUL(sblock, krotj))
|
||||
END IF
|
||||
END DO
|
||||
CALL dbcsr_iterator_stop(iter)
|
||||
!
|
||||
END IF
|
||||
ELSE
|
||||
! this is the identity operation, just copy the matrix
|
||||
|
|
|
|||
|
|
@ -53,6 +53,7 @@ MODULE kpoint_types
|
|||
PUBLIC :: read_kpoint_section, write_kpoint_info
|
||||
PUBLIC :: kpoint_env_type, kpoint_env_p_type
|
||||
PUBLIC :: kpoint_env_create, get_kpoint_env
|
||||
PUBLIC :: kind_rotmat_type
|
||||
PUBLIC :: kpoint_sym_type
|
||||
PUBLIC :: kpoint_sym_create
|
||||
|
||||
|
|
@ -85,6 +86,15 @@ MODULE kpoint_types
|
|||
TYPE(kpoint_env_type), POINTER :: kpoint_env => NULL()
|
||||
END TYPE kpoint_env_p_type
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Rotation matrices for basis sets
|
||||
!> \param rmat atom basis function rotation matrix
|
||||
!> \author JGH
|
||||
! **************************************************************************************************
|
||||
TYPE kind_rotmat_type
|
||||
REAL(KIND=dp), DIMENSION(:, :), POINTER :: rmat
|
||||
END TYPE kind_rotmat_type
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Keeps symmetry information about a specific k-point
|
||||
!> \param apply_symmetry ...
|
||||
|
|
@ -97,8 +107,10 @@ MODULE kpoint_types
|
|||
TYPE kpoint_sym_type
|
||||
LOGICAL :: apply_symmetry = .FALSE.
|
||||
INTEGER :: nwght = -1
|
||||
INTEGER :: nwred = -1
|
||||
REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp => NULL()
|
||||
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rot => NULL()
|
||||
INTEGER, DIMENSION(:), POINTER :: rotp => NULL()
|
||||
INTEGER, DIMENSION(:, :), POINTER :: f0 => NULL()
|
||||
END TYPE kpoint_sym_type
|
||||
|
||||
|
|
@ -168,8 +180,13 @@ MODULE kpoint_types
|
|||
POINTER :: kp_env => Null()
|
||||
TYPE(kpoint_env_p_type), DIMENSION(:), &
|
||||
POINTER :: kp_aux_env => Null()
|
||||
! symmetry
|
||||
TYPE(kpoint_sym_p_type), DIMENSION(:), &
|
||||
POINTER :: kp_sym => Null()
|
||||
INTEGER, DIMENSION(:), POINTER :: atype => Null()
|
||||
TYPE(kind_rotmat_type), DIMENSION(:, :), &
|
||||
POINTER :: kind_rotmat => Null()
|
||||
! pools
|
||||
TYPE(qs_matrix_pools_type), POINTER :: mpools => Null()
|
||||
TYPE(qs_diis_buffer_type_kp), POINTER :: scf_diis_buffer => Null()
|
||||
TYPE(qs_matrix_pools_type), POINTER :: mpools_aux_fit => Null()
|
||||
|
|
@ -233,7 +250,7 @@ CONTAINS
|
|||
SUBROUTINE kpoint_release(kpoint)
|
||||
TYPE(kpoint_type), POINTER :: kpoint
|
||||
|
||||
INTEGER :: ik
|
||||
INTEGER :: i, ik, j
|
||||
|
||||
IF (ASSOCIATED(kpoint)) THEN
|
||||
|
||||
|
|
@ -281,6 +298,19 @@ CONTAINS
|
|||
DEALLOCATE (kpoint%kp_sym)
|
||||
END IF
|
||||
|
||||
IF (ASSOCIATED(kpoint%atype)) DEALLOCATE (kpoint%atype)
|
||||
|
||||
IF (ASSOCIATED(kpoint%kind_rotmat)) THEN
|
||||
DO i = 1, SIZE(kpoint%kind_rotmat, 1)
|
||||
DO j = 1, SIZE(kpoint%kind_rotmat, 2)
|
||||
IF (ASSOCIATED(kpoint%kind_rotmat(i, j)%rmat)) THEN
|
||||
DEALLOCATE (kpoint%kind_rotmat(i, j)%rmat)
|
||||
END IF
|
||||
END DO
|
||||
END DO
|
||||
DEALLOCATE (kpoint%kind_rotmat)
|
||||
END IF
|
||||
|
||||
IF (ASSOCIATED(kpoint%scf_diis_buffer)) THEN
|
||||
CALL qs_diis_b_release_kp(kpoint%scf_diis_buffer)
|
||||
DEALLOCATE (kpoint%scf_diis_buffer)
|
||||
|
|
@ -773,10 +803,12 @@ CONTAINS
|
|||
ALLOCATE (kp_sym)
|
||||
|
||||
kp_sym%nwght = 0
|
||||
kp_sym%nwred = 0
|
||||
kp_sym%apply_symmetry = .FALSE.
|
||||
|
||||
NULLIFY (kp_sym%rot)
|
||||
NULLIFY (kp_sym%xkp)
|
||||
NULLIFY (kp_sym%rotp)
|
||||
NULLIFY (kp_sym%f0)
|
||||
|
||||
END SUBROUTINE kpoint_sym_create
|
||||
|
|
@ -800,6 +832,9 @@ CONTAINS
|
|||
IF (ASSOCIATED(kp_sym%f0)) THEN
|
||||
DEALLOCATE (kp_sym%f0)
|
||||
END IF
|
||||
IF (ASSOCIATED(kp_sym%rotp)) THEN
|
||||
DEALLOCATE (kp_sym%rotp)
|
||||
END IF
|
||||
|
||||
DEALLOCATE (kp_sym)
|
||||
|
||||
|
|
|
|||
3070
src/kpsym.F
Normal file
3070
src/kpsym.F
Normal file
File diff suppressed because it is too large
Load diff
145
src/qs_basis_rotation_methods.F
Normal file
145
src/qs_basis_rotation_methods.F
Normal file
|
|
@ -0,0 +1,145 @@
|
|||
!--------------------------------------------------------------------------------------------------!
|
||||
! CP2K: A general program to perform molecular dynamics simulations !
|
||||
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
|
||||
! !
|
||||
! SPDX-License-Identifier: GPL-2.0-or-later !
|
||||
!--------------------------------------------------------------------------------------------------!
|
||||
|
||||
! **************************************************************************************************
|
||||
MODULE qs_basis_rotation_methods
|
||||
USE basis_set_types, ONLY: get_gto_basis_set,&
|
||||
gto_basis_set_type
|
||||
USE cell_types, ONLY: cell_type
|
||||
USE cp_control_types, ONLY: dft_control_type
|
||||
USE input_constants, ONLY: do_method_dftb
|
||||
USE kinds, ONLY: dp
|
||||
USE kpoint_types, ONLY: kpoint_sym_type,&
|
||||
kpoint_type
|
||||
USE orbital_pointers, ONLY: nso
|
||||
USE orbital_transformation_matrices, ONLY: calculate_rotmat,&
|
||||
orbrotmat_type,&
|
||||
release_rotmat
|
||||
USE qs_environment_types, ONLY: get_qs_env,&
|
||||
qs_environment_type
|
||||
USE qs_kind_types, ONLY: get_qs_kind,&
|
||||
get_qs_kind_set,&
|
||||
qs_kind_type
|
||||
#include "./base/base_uses.f90"
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
PRIVATE
|
||||
|
||||
! Global parameters (only in this module)
|
||||
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_basis_rotation_methods'
|
||||
|
||||
! Public subroutines
|
||||
|
||||
PUBLIC :: qs_basis_rotation
|
||||
|
||||
CONTAINS
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Construct basis set rotation matrices
|
||||
!> \param qs_env ...
|
||||
!> \param kpoints ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE qs_basis_rotation(qs_env, kpoints)
|
||||
|
||||
TYPE(qs_environment_type), POINTER :: qs_env
|
||||
TYPE(kpoint_type), POINTER :: kpoints
|
||||
|
||||
INTEGER :: ik, ikind, ir, ira, irot, lval, nkind, &
|
||||
nrot
|
||||
REAL(KIND=dp), DIMENSION(3, 3) :: rotmat
|
||||
TYPE(cell_type), POINTER :: cell
|
||||
TYPE(dft_control_type), POINTER :: dft_control
|
||||
TYPE(gto_basis_set_type), POINTER :: orb_basis
|
||||
TYPE(kpoint_sym_type), POINTER :: kpsym
|
||||
TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrot
|
||||
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
|
||||
|
||||
CPASSERT(ASSOCIATED(qs_env))
|
||||
CPASSERT(ASSOCIATED(kpoints))
|
||||
IF (ASSOCIATED(kpoints%kind_rotmat)) THEN
|
||||
CALL get_qs_env(qs_env, cell=cell)
|
||||
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
|
||||
CALL get_qs_kind_set(qs_kind_set, maxlgto=lval)
|
||||
nrot = SIZE(kpoints%kind_rotmat, 1)
|
||||
nkind = SIZE(kpoints%kind_rotmat, 2)
|
||||
! remove possible old rotation matrices
|
||||
DO irot = 1, nrot
|
||||
DO ikind = 1, nkind
|
||||
IF (ASSOCIATED(kpoints%kind_rotmat(irot, ikind)%rmat)) THEN
|
||||
DEALLOCATE (kpoints%kind_rotmat(irot, ikind)%rmat)
|
||||
END IF
|
||||
END DO
|
||||
END DO
|
||||
! check all rotations needed
|
||||
NULLIFY (orbrot)
|
||||
CALL get_qs_env(qs_env, dft_control=dft_control)
|
||||
DO ik = 1, kpoints%nkp
|
||||
kpsym => kpoints%kp_sym(ik)%kpoint_sym
|
||||
IF (kpsym%apply_symmetry) THEN
|
||||
DO irot = 1, SIZE(kpsym%rotp)
|
||||
ir = kpsym%rotp(irot)
|
||||
ira = ABS(ir)
|
||||
IF (ira > 0 .AND. ira <= nrot) THEN
|
||||
IF (.NOT. ASSOCIATED(kpoints%kind_rotmat(ira, 1)%rmat)) THEN
|
||||
rotmat(1:3, 1:3) = MATMUL(cell%h_inv, &
|
||||
MATMUL(kpsym%rot(:, :, irot), cell%hmat))
|
||||
NULLIFY (orbrot)
|
||||
CALL calculate_rotmat(orbrot, rotmat, lval)
|
||||
IF (dft_control%qs_control%method_id == do_method_dftb) THEN
|
||||
CPABORT("ROTMAT")
|
||||
ELSE
|
||||
DO ikind = 1, nkind
|
||||
CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis)
|
||||
NULLIFY (kpoints%kind_rotmat(ira, ikind)%rmat)
|
||||
CALL set_rotmat_basis(kpoints%kind_rotmat(ira, ikind)%rmat, orbrot, orb_basis)
|
||||
END DO
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
END IF
|
||||
END DO
|
||||
CALL release_rotmat(orbrot)
|
||||
END IF
|
||||
|
||||
END SUBROUTINE qs_basis_rotation
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param rmat ...
|
||||
!> \param orbrot ...
|
||||
!> \param basis ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE set_rotmat_basis(rmat, orbrot, basis)
|
||||
REAL(KIND=dp), DIMENSION(:, :), POINTER :: rmat
|
||||
TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrot
|
||||
TYPE(gto_basis_set_type), POINTER :: basis
|
||||
|
||||
INTEGER :: fs1, fs2, iset, ishell, l, nset, nsgf
|
||||
INTEGER, DIMENSION(:), POINTER :: nshell
|
||||
INTEGER, DIMENSION(:, :), POINTER :: first_sgf, lshell
|
||||
|
||||
CALL get_gto_basis_set(gto_basis_set=basis, nsgf=nsgf)
|
||||
ALLOCATE (rmat(nsgf, nsgf))
|
||||
rmat = 0.0_dp
|
||||
|
||||
CALL get_gto_basis_set(gto_basis_set=basis, nset=nset, nshell=nshell, l=lshell, &
|
||||
first_sgf=first_sgf)
|
||||
DO iset = 1, nset
|
||||
DO ishell = 1, nshell(iset)
|
||||
l = lshell(ishell, iset)
|
||||
fs1 = first_sgf(ishell, iset)
|
||||
fs2 = fs1 + nso(l) - 1
|
||||
rmat(fs1:fs2, fs1:fs2) = orbrot(l)%mat(1:nso(l), 1:nso(l))
|
||||
END DO
|
||||
END DO
|
||||
|
||||
END SUBROUTINE set_rotmat_basis
|
||||
|
||||
END MODULE qs_basis_rotation_methods
|
||||
|
|
@ -132,6 +132,7 @@ MODULE qs_environment
|
|||
USE particle_types, ONLY: particle_type
|
||||
USE pw_env_types, ONLY: pw_env_type
|
||||
USE qmmm_types_low, ONLY: qmmm_env_qm_type
|
||||
USE qs_basis_rotation_methods, ONLY: qs_basis_rotation
|
||||
USE qs_dftb_parameters, ONLY: qs_dftb_param_init
|
||||
USE qs_dftb_types, ONLY: qs_dftb_pairpot_type
|
||||
USE qs_dispersion_nonloc, ONLY: qs_dispersion_nonloc_init
|
||||
|
|
@ -391,6 +392,10 @@ CONTAINS
|
|||
CALL get_qs_env(qs_env=qs_env, wf_history=wf_history)
|
||||
CALL wfi_create_for_kp(wf_history)
|
||||
END IF
|
||||
! basis set symmetry rotations
|
||||
IF (do_kpoints) THEN
|
||||
CALL qs_basis_rotation(qs_env, kpoints)
|
||||
END IF
|
||||
|
||||
do_hfx = .FALSE.
|
||||
hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
|
||||
|
|
|
|||
|
|
@ -8,4 +8,5 @@ bn1.inp 94 1.0E-14
|
|||
#cc2.inp 1 1.0E-12 -45.38989062944852
|
||||
#cc3.inp 1 1.0E-12 -45.38990176407365
|
||||
#cc4.inp 1 1.0E-12 -45.38975726647867
|
||||
#cc5.inp 1 1.0E-12 -45.38975726647867
|
||||
#EOF
|
||||
|
|
|
|||
68
tests/QS/regtest-kp-3/cc5.inp
Normal file
68
tests/QS/regtest-kp-3/cc5.inp
Normal file
|
|
@ -0,0 +1,68 @@
|
|||
# Reference energy -45.2700189877
|
||||
&GLOBAL
|
||||
PRINT_LEVEL LOW
|
||||
PROJECT C
|
||||
RUN_TYPE ENERGY
|
||||
&END GLOBAL
|
||||
|
||||
&FORCE_EVAL
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME GTH_BASIS_SETS
|
||||
POTENTIAL_FILE_NAME POTENTIAL
|
||||
&KPOINTS
|
||||
EPS_GEO 1.e-8
|
||||
FULL_GRID OFF
|
||||
PARALLEL_GROUP_SIZE 0
|
||||
SCHEME MONKHORST-PACK 2 2 2
|
||||
SYMMETRY ON
|
||||
VERBOSE T
|
||||
&END KPOINTS
|
||||
&MGRID
|
||||
CUTOFF 280
|
||||
REL_CUTOFF 50
|
||||
&END MGRID
|
||||
&QS
|
||||
EPS_DEFAULT 1.0E-14
|
||||
METHOD GPW
|
||||
&END QS
|
||||
&SCF
|
||||
CHOLESKY OFF
|
||||
EPS_EIGVAL 1.e-8
|
||||
EPS_SCF 1.0E-8
|
||||
MAX_SCF 20
|
||||
SCF_GUESS ATOMIC
|
||||
&MIXING
|
||||
ALPHA 0.50
|
||||
METHOD BROYDEN_MIXING
|
||||
&END MIXING
|
||||
&PRINT
|
||||
&RESTART off
|
||||
&END RESTART
|
||||
&END PRINT
|
||||
&END SCF
|
||||
&XC
|
||||
&XC_FUNCTIONAL PADE
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 3.56683 3.56683 3.56683
|
||||
&END CELL
|
||||
&COORD
|
||||
SCALED
|
||||
C 0.100000 0.000000 0.000000
|
||||
C 0.500000 0.500000 0.000000
|
||||
C 0.500000 0.000000 0.500000
|
||||
C 0.000000 0.500000 0.500000
|
||||
C 0.250000 0.250000 0.250000
|
||||
C 0.250000 0.750000 0.750000
|
||||
C 0.750000 0.250000 0.750000
|
||||
C 0.750000 0.750000 0.250000
|
||||
&END COORD
|
||||
&KIND C
|
||||
BASIS_SET SZV-GTH
|
||||
POTENTIAL GTH-PADE-q4
|
||||
&END KIND
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
Loading…
Add table
Add a link
Reference in a new issue