Stengel-Spalding with guiding centres

[JeromeCCP9]

• check wann_phases() implementation correctness for Stengel-Spalding case
• add new benchmark for higher-order-stengel-spalding test case
• remove error exit for this combination

see PR #584 for more detail.

[sjhong6230]

@JeromeCCP9

Dear Jerome,

I looked into the code and found the reason for the failure. There were two mistakes.

The first one is that the Cartesian inner product is missing in the following line. It is my silly mistake :)

wannier90/src/wannierise.F90

Lines 1916 to 1923 in 77b35f7

	if (use_ss_functional) then
	do nn = 1, kmesh_info%nntot
	do loop_wann = 1, num_wann
	do j = 1, 3
	sheet(loop_wann, nn, 1) = kmesh_info%bk(j, nn, 1)*rguide(j, loop_wann)
	enddo
	enddo
	enddo

The second one is trickier. In the following lines in wann_domega subroutine, the sum_mmn matrices are multiplied by csheet, while m_matrix_loc matrices are not. This results in inconsistent phase between m matrices and results in a wrong output.

wannier90/src/wannierise.F90

Lines 2508 to 2533 in 77b35f7

	cdodq_loc = cmplx_0
	do nkp_loc = 1, nkrank
	do n = 1, num_wann
	do m = 1, num_wann
	do nn = 1, kmesh_info%nntot
	nkp = global_k(nkp_loc)
	cnn = kmesh_info%nnord(nn, nkp) ! enforce uniform order of bk vectors
	cnn2 = kmesh_info%nnrev(nn, nkp) ! (b-vector nn2) is opposite to/negative of (b-vector nn)
	cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) + &
	kmesh_info%wb(nn)*m_matrix_loc(m, n, cnn, nkp_loc)* &
	conjg(sum_mnn(n, nn))
	cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - &
	kmesh_info%wb(nn)*conjg(m_matrix_loc(n, m, cnn2, nkp_loc))* &
	conjg(sum_mnn(m, nn))
	cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - &
	kmesh_info%wb(nn)*conjg(m_matrix_loc(n, m, cnn, nkp_loc))* &
	sum_mnn(m, nn)
	cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) + &
	kmesh_info%wb(nn)*m_matrix_loc(m, n, cnn2, nkp_loc)* &
	sum_mnn(n, nn)
	enddo
	enddo
	enddo
	enddo
	cdodq_loc = cdodq_loc/real(num_kpts, dp)

Fixing the line as below gives the same results for KNbO3 with guiding_centres = T and guiding_centres = F.

      cdodq_loc = cmplx_0
      do nkp_loc = 1, nkrank
        do n = 1, num_wann
          do m = 1, num_wann
            do nn = 1, kmesh_info%nntot
              nkp = global_k(nkp_loc)
              cnn = kmesh_info%nnord(nn, nkp) ! enforce uniform order of bk vectors
              cnn2 = kmesh_info%nnrev(nn, nkp) ! (b-vector nn2) is opposite to/negative of (b-vector nn)

              cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) + &
                                         kmesh_info%wb(nn)*m_matrix_loc(m, n, cnn, nkp_loc)* &
                                         conjg(sum_mnn(n, nn)/csheet(n, nn, 1))
              cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - &
                                         kmesh_info%wb(nn)*conjg(m_matrix_loc(n, m, cnn2, nkp_loc))* &
                                         conjg(sum_mnn(m, nn)/csheet(m, nn, 1))
              cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - &
                                         kmesh_info%wb(nn)*conjg(m_matrix_loc(n, m, cnn, nkp_loc))* &
                                         sum_mnn(m, nn)/csheet(m, nn, 1)
              cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) + &
                                         kmesh_info%wb(nn)*m_matrix_loc(m, n, cnn2, nkp_loc)* &
                                         sum_mnn(n, nn)/csheet(n, nn, 1)
            enddo
          enddo
        enddo
      enddo

I also attached the wannierization results.

guiding_false.txt
guiding_true.txt

I checked that all tests are passed, and the code is ready for pull request.

[JeromeCCP9]

Ahh, indeed! This is very nice, thank you also for the explanation. I would not have been able to identify this phase issue without spending a lot of time, so I really appreciate your fixing this.

I implemented the fix (in commit f5640a3) which will soon be merged (I mentioned you in the commit log along with this issue number, I hope that's ok).
J.