subroutine ModifyMomentum(FDet)
integer :: i, j ! Loop variables
integer, intent(inout) :: FDet(NEl)
integer :: k_total(3) ! Stores the total momentum of FDet
integer :: delta_k(3) ! Stores the difference between the current FDet and the FDet we're aiming for
integer, allocatable :: kPointToBasisFn(:, :, :, :) ! Look up table for kPoints to basis functions
integer :: kmaxX, kminX, kmaxY, kminY, kmaxZ, kminZ, iSpinIndex ! Stores the limits of kPointToBasisFn
integer :: det_sorted(NEl), e_store ! Storage for the sorting routine
logical :: sorted ! As above
integer :: wrapped_index
integer :: k_new
IF (.not. tUEG) call stop_all("ModifyMomentum", "Only works for UEG")
! Finds current momentum, and finds the difference between this and the target momentum
! most commonly will be zero to start with
k_total(1) = 0
k_total(2) = 0
k_total(3) = 0
do j = 1, NEl
k_total(1) = k_total(1) + G1(FDet(j))%k(1)
k_total(2) = k_total(2) + G1(FDet(j))%k(2)
k_total(3) = k_total(3) + G1(FDet(j))%k(3)
end do
delta_k = k_momentum - k_total
if (delta_k(1) == 0 .and. delta_k(2) == 0 .and. delta_k(3) == 0) write (stdout, *) "WARNING: specified momentum is ground state"
! Creates a look-up table for k-points (this was taken from symrandexcit2.F90)
kmaxX = 0
kminX = 0
kmaxY = 0
kminY = 0
kminZ = 0
kmaxZ = 0
do i = 1, nBasis
IF (G1(i)%k(1) > kmaxX) kmaxX = G1(i)%k(1)
IF (G1(i)%k(1) < kminX) kminX = G1(i)%k(1)
IF (G1(i)%k(2) > kmaxY) kmaxY = G1(i)%k(2)
IF (G1(i)%k(2) < kminY) kminY = G1(i)%k(2)
IF (G1(i)%k(3) > kmaxZ) kmaxZ = G1(i)%k(3)
IF (G1(i)%k(3) < kminZ) kminZ = G1(i)%k(3)
end do
allocate (kPointToBasisFn(kminX:kmaxX, kminY:kmaxY, kminZ:kmaxZ, 2))
do i = 1, nBasis
iSpinIndex = (G1(i)%Ms + 1) / 2 + 1 ! iSpinIndex equals 1 for a beta spin (ms=-1), and 2 for an alpha spin (ms=1)
kPointToBasisFn(G1(i)%k(1), G1(i)%k(2), G1(i)%k(3), iSpinIndex) = i
end do
! For each of the three dimensions, nudge electrons one at a time by one momentum unit until delta_k is reached
det_sorted = FDet
! Bubble sort to order det_sorted in order of kx of the corresponding electron
do
sorted = .true.
do i = 1, NEl - 1
j = i + 1
if (G1(det_sorted(j))%k(1) > G1(det_sorted(i))%k(1)) then
sorted = .false.
e_store = det_sorted(i)
det_sorted(i) = det_sorted(j)
det_sorted(j) = e_store
end if
end do
if (sorted) exit
end do
! Nudge momenta one at a time
if (delta_k(1) > 0) then
do i = 1, delta_k(1)
wrapped_index = mod(i, NEl) ! Take the modulus to know which electron to nudge
if (wrapped_index == 0) then ! Deal with the i=NEl case
wrapped_index = NEl
end if
j = wrapped_index ! For convenience asign this to j
k_new = G1(det_sorted(j))%k(1) + 1 ! Find the new momentum of this electron
if (k_new > kmaxX) then ! Check that this momentum isn't outside the cell
call stop_all("ModifyMomentum", "Electron moved outside of the cell limits")
end if
iSpinIndex = (G1(j)%Ms + 1) / 2 + 1 ! Spin of the new orbital is the same as the old
! Finds basis number for the new momentum
det_sorted(j) = kPointToBasisFn(k_new, G1(det_sorted(j))%k(2), G1(det_sorted(j))%k(3), iSpinIndex)
end do
else if (delta_k(1) < 0) then ! For the negative case, i must run through negative numbers
do i = -1, delta_k(1), -1
wrapped_index = mod(i, NEl)
if (wrapped_index == 0) then
wrapped_index = -NEl
end if
j = NEl + wrapped_index + 1 ! Now this goes through the list backward (wrapped_index is negative)
k_new = G1(det_sorted(j))%k(1) - 1 ! Find the new momentum of this electron, this time in the opposite direction
if (k_new < kminX) then ! Check the limits of the cell again
call stop_all("ModifyMomentum", "Electron moved outside of the cell limits")
end if
iSpinIndex = (G1(j)%Ms + 1) / 2 + 1 ! Spin of the new orbital is the same as the old
! Finds basis number for the new momentum
det_sorted(j) = kPointToBasisFn(k_new, G1(det_sorted(j))%k(2), G1(det_sorted(j))%k(3), iSpinIndex)
end do
end if
FDet = det_sorted
!====ky treated as kx above
do
sorted = .true.
do i = 1, NEl - 1
j = i + 1
if (G1(det_sorted(j))%k(2) > G1(det_sorted(i))%k(2)) then
sorted = .false.
e_store = det_sorted(i)
det_sorted(i) = det_sorted(j)
det_sorted(j) = e_store
end if
end do
if (sorted) exit
end do
! Nudge momenta one at a time
if (delta_k(2) > 0) then
do i = 1, delta_k(2)
wrapped_index = mod(i, NEl) ! Take the modulus to know which electron to nudge
if (wrapped_index == 0) then ! Deal with the i=NEl case
wrapped_index = NEl
end if
j = wrapped_index ! For convenience asign this to j
k_new = G1(det_sorted(j))%k(2) + 1 ! Find the new momentum of this electron
if (k_new > kmaxY) then ! Check that this momentum isn't outside the cell
call stop_all("ModifyMomentum", "Electron moved outside of the cell limits")
end if
iSpinIndex = (G1(j)%Ms + 1) / 2 + 1 ! Spin of the new orbital is the same as the old
! Finds basis number for the new momentum
det_sorted(j) = kPointToBasisFn(G1(det_sorted(j))%k(1), k_new, G1(det_sorted(j))%k(3), iSpinIndex)
end do
else if (delta_k(2) < 0) then ! For the negative case, i must run through negative numbers
do i = -1, delta_k(2), -1
wrapped_index = mod(i, NEl)
if (wrapped_index == 0) then
wrapped_index = -NEl
end if
j = NEl + wrapped_index + 1 ! Now this goes through the list backward (wrapped_index is negative)
k_new = G1(det_sorted(j))%k(2) - 1 ! Find the new momentum of this electron, this time in the opposite direction
if (k_new < kminY) then ! Check the limits of the cell again
call stop_all("ModifyMomentum", "Electron moved outside of the cell limits")
end if
iSpinIndex = (G1(j)%Ms + 1) / 2 + 1 ! Spin of the new orbital is the same as the old
! Finds basis number for the new momentum
det_sorted(j) = kPointToBasisFn(G1(det_sorted(j))%k(1), k_new, G1(det_sorted(j))%k(3), iSpinIndex)
end do
end if
FDet = det_sorted
!====kz treated as kx and ky above
do
sorted = .true.
do i = 1, NEl - 1
j = i + 1
if (G1(det_sorted(j))%k(3) > G1(det_sorted(i))%k(3)) then
sorted = .false.
e_store = det_sorted(i)
det_sorted(i) = det_sorted(j)
det_sorted(j) = e_store
end if
end do
if (sorted) exit
end do
! Nudge momenta one at a time
if (delta_k(3) > 0) then
do i = 1, delta_k(3)
wrapped_index = mod(i, NEl) ! Take the modulus to know which electron to nudge
if (wrapped_index == 0) then ! Deal with the i=NEl case
wrapped_index = NEl
end if
j = wrapped_index ! For convenience asign this to j
k_new = G1(det_sorted(j))%k(3) + 1 ! Find the new momentum of this electron
if (k_new > kmaxZ) then ! Check that this momentum isn't outside the cell
call stop_all("ModifyMomentum", "Electron moved outside of the cell limits")
end if
iSpinIndex = (G1(j)%Ms + 1) / 2 + 1 ! Spin of the new orbital is the same as the old
! Finds basis number for the new momentum
det_sorted(j) = kPointToBasisFn(G1(det_sorted(j))%k(1), G1(det_sorted(j))%k(2), k_new, iSpinIndex)
end do
else if (delta_k(3) < 0) then ! For the negative case, i must run through negative numbers
do i = -1, delta_k(3), -1
wrapped_index = mod(i, NEl)
if (wrapped_index == 0) then
wrapped_index = -NEl
end if
j = NEl + wrapped_index + 1 ! Now this goes through the list backward (wrapped_index is negative)
k_new = G1(det_sorted(j))%k(3) - 1 ! Find the new momentum of this electron, this time in the opposite direction
if (k_new < kminZ) then ! Check the limits of the cell again
call stop_all("ModifyMomentum", "Electron moved outside of the cell limits")
end if
iSpinIndex = (G1(j)%Ms + 1) / 2 + 1 ! Spin of the new orbital is the same as the old
! Finds basis number for the new momentum
det_sorted(j) = kPointToBasisFn(G1(det_sorted(j))%k(1), G1(det_sorted(j))%k(2), k_new, iSpinIndex)
end do
end if
FDet = det_sorted
! Bubble sort to order FDet back into increasing order by number
do
sorted = .true.
do i = 1, NEl - 1
j = i + 1
if (FDet(j) < FDet(i)) then
sorted = .false.
e_store = FDet(i)
FDet(i) = FDet(j)
FDet(j) = e_store
end if
end do
if (sorted) exit
end do
write (stdout, *) "Total momentum set to", k_momentum
end subroutine ModifyMomentum