pure SUBROUTINE MOMPBCSYM(K1, NBASISMAX)
! NB the third column of NBASISMAX tells us whether it is tilted
INTEGER, intent(inout) :: K1(3)
INTEGER, intent(in) :: nBasisMax(5, *)
INTEGER J, LDIM, AX, AY, LENX, LENY, KK2, T1, T2
real(dp) R1, R2, NORM
! [W.D]
! in case of the new k-space hubbard implementation use the
! built-in k-vec mapping
if (t_k_space_hubbard) then
! this should actually never be called anymore!
k1 = lat%map_k_vec(k1)
return
end if
! (AX,AY) is the tilt of the lattice, and corresponds to the lattice vector of the cell. The other lattice vector is (-AY,AX).
!These are expressed in terms of the primitive Hubbard lattice vectors
AX = NBASISMAX(1, 4)
AY = NBASISMAX(2, 4)
! LENX is the length (i.e. number of lattice vectors) in direction (AX,AY). LENY is in the other lattice vector direction
LENX = NBASISMAX(1, 5)
LENY = NBASISMAX(2, 5)
IF (NBASISMAX(1, 3) == 0 .OR. NBASISMAX(1, 3) == 0) THEN
! A non-tilted lattice with PBC
DO J = 1, 3
! non-tilted
KK2 = K1(J)
LDIM = NBASISMAX(J, 2) - NBASISMAX(J, 1) + 1
KK2 = MOD(KK2, LDIM)
IF (KK2 < NBASISMAX(J, 1)) KK2 = KK2 + LDIM
IF (KK2 > NBASISMAX(J, 2)) KK2 = KK2 - LDIM
K1(J) = KK2
end do
else if (NBASISMAX(1, 3) == 1) THEN
! we have a tilted lattice with PBC
! we want the a1,a2 components of k
NORM = AX * AX + AY * AY
R1 = (AX * K1(1) + AY * K1(2)) / NORM
R2 = (AX * K1(2) - AY * K1(1)) / NORM
R1 = R1 / LENX + 0.5_dp
R2 = R2 / LENY + 0.5_dp
! R1 is now in terms of the shifted extended unit cell. We want 0<R1<=1
! R2 is now in terms of the shifted extended unit cell. We want 0<=R2<1
! T1 and T2 will be the vectors to translate the point K1 back into our cell
T1 = INT(ABS(R1))
IF (R1 < 0.0_dp) THEN
T1 = -T1
IF (abs(t1 - r1) > 1.0e-12_dp) T1 = T1 - 1
end if
! Now T1= highest integer less than R1 (i.e. FLOOR)
! We want to include R1=1, so we have one fewer translations in that case.
IF (abs(r1 - t1) < 1.0e-12_dp) T1 = T1 - 1
T2 = INT(ABS(R2))
IF (R2 < 0.0_dp) THEN
T2 = -T2
IF (abs(t2 - r2) > 1.0e-12_dp) T2 = T2 - 1
end if
! The following line conserves the top-right rule for edge effects
IF (abs(r1 - t1) < 1.0e-12_dp) T1 = T1 - 1
! Now T2= highest integer less than R2 (i.e. FLOOR)
! Do the translation
IF (R1 > 1.0_dp .OR. R1 <= 0.0_dp) R1 = R1 - T1
IF (R2 >= 1.0_dp .OR. R2 < 0.0_dp) R2 = R2 - T2
! Now convert back to the the Coords in the extended brillouin zone (ish)
R1 = (R1 - 0.5_dp) * LENX
R2 = (R2 - 0.5_dp) * LENY
! Now K1 is defined as the re-scaled (R1,R2) vector
K1(1) = NINT(R1 * AX - R2 * AY)
K1(2) = NINT(R1 * AY + R2 * AX)
else if (NBASISMAX(1, 3) == -1) THEN
if ((abs(k1(1)) > NBASISMAX(1, 2)) .and. tperiodicinmom) then
if (k1(1) > 0) then
k1(1) = k1(1) - (2 * NBASISMAX(1, 2) + 1)
else
k1(1) = k1(1) + (2 * NBASISMAX(1, 2) + 1)
end if
end if
end if
END SUBROUTINE MOMPBCSYM