subroutine wann_main_gamma(kmesh_info, wann_control, omega, print_output, wannier_data, &
m_matrix, u_matrix, real_lattice, num_kpts, num_wann, stdout, timer, &
error, comm)
!================================================!
!
! Calculate the Unitary Rotations to give
! Maximally Localised Wannier Functions
! Gamma version
!================================================
use w90_constants, only: dp, cmplx_1, cmplx_0
use w90_io, only: io_time, io_stopwatch_start, io_stopwatch_stop
use w90_wannier90_types, only: wann_control_type, wann_omega_type
use w90_types, only: kmesh_info_type, print_output_type, &
wannier_data_type, timer_list_type
use w90_wannier90_readwrite, only: w90_wannier90_readwrite_write_chkpt
use w90_utility, only: utility_frac_to_cart, utility_zgemm
use w90_comms, only: w90_comm_type, mpirank
implicit none
! JJ this function is entirely serial
! note that the scaling may be inferior to the non-gamma case
! so this is not necessarily the quicker branch
! JJ surely numkpts==1 identically???
! arguments
type(wannier_data_type), intent(inout) :: wannier_data
type(w90_comm_type), intent(in) :: comm
type(wann_control_type), intent(inout) :: wann_control
type(wann_omega_type), intent(inout) :: omega
type(print_output_type), intent(in) :: print_output
type(kmesh_info_type), intent(in) :: kmesh_info
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), allocatable, intent(out) :: error
integer, intent(in) :: stdout
integer, intent(in) :: num_wann
integer, intent(in) :: num_kpts
real(kind=dp), intent(in) :: real_lattice(3, 3)
complex(kind=dp), intent(inout) :: u_matrix(:, :, :)
complex(kind=dp), intent(inout) :: m_matrix(:, :, :, :)
! local variables
type(localisation_vars_type) :: old_spread
type(localisation_vars_type) :: wann_spread
integer :: counts(0:1)
integer :: global_k(1)
real(kind=dp), allocatable :: rnkb(:, :, :)
real(kind=dp), allocatable :: ln_tmp(:, :, :)
logical :: first_pass
! guiding centres
real(kind=dp), allocatable :: rguide(:, :)
integer :: irguide
! local arrays used and passed in subroutines
real(kind=dp), allocatable :: m_w(:, :, :)
complex(kind=dp), allocatable :: csheet(:, :, :)
real(kind=dp), allocatable :: sheet(:, :, :)
real(kind=dp), allocatable :: rave(:, :), r2ave(:), rave2(:)
!local arrays not passed into subroutines
complex(kind=dp), allocatable :: u0(:, :, :)
complex(kind=dp), allocatable :: uc_rot(:, :)
real(kind=dp), allocatable :: ur_rot(:, :)
complex(kind=dp), allocatable :: cz(:, :)
real(kind=dp) :: sqwb
integer :: i, n, nn, iter, ind, ierr, iw
integer :: tnntot
logical :: lprint, ldump
real(kind=dp), allocatable :: history(:)
logical :: lconverged, lrandom, lfirst
real(kind=dp) :: save_spread
integer :: conv_count, noise_count
if (mpirank(comm) > 0) then
! this cannot happen under ordinary circumstances
call set_error_alloc(error, &
'wann_main_gamma called by non-root rank (but the algorithm is serial)', comm)
return
end if
if (print_output%timing_level > 0) call io_stopwatch_start('wann: main_gamma', timer)
first_pass = .true.
allocate (history(wann_control%conv_window), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error allocating history in wann_main_gamma', comm)
return
end if
allocate (rnkb(num_wann, kmesh_info%nntot, num_kpts), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating rnkb in wann_main_gamma', comm)
return
end if
allocate (ln_tmp(num_wann, kmesh_info%nntot, num_kpts), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating ln_tmp in wann_main_gamma', comm)
return
end if
rnkb = 0.0_dp
tnntot = 2*kmesh_info%nntot
allocate (m_w(num_wann, num_wann, tnntot), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating m_w in wann_main_gamma', comm)
return
end if
allocate (csheet(num_wann, kmesh_info%nntot, num_kpts), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating csheet in wann_main_gamma', comm)
return
end if
allocate (sheet(num_wann, kmesh_info%nntot, num_kpts), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sheet in wann_main_gamma', comm)
return
end if
allocate (rave(3, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating rave in wann_main_gamma', comm)
return
end if
allocate (r2ave(num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating r2ave in wann_main_gamma', comm)
return
end if
allocate (rave2(num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating rave2 in wann_main_gamma', comm)
return
end if
allocate (rguide(3, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating rguide in wann_main_gamma', comm)
return
end if
csheet = cmplx_1
sheet = 0.0_dp; rave = 0.0_dp; r2ave = 0.0_dp; rave2 = 0.0_dp; rguide = 0.0_dp
allocate (u0(num_wann, num_wann, num_kpts), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating u0 in wann_main_gamma', comm)
return
end if
allocate (uc_rot(num_wann, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating uc_rot in wann_main_gamma', comm)
return
end if
allocate (ur_rot(num_wann, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating ur_rot in wann_main_gamma', comm)
return
end if
allocate (cz(num_wann, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating cz in wann_main_gamma', comm)
return
end if
cz = cmplx_0
! Set up the MPI arrays for a serial run.
counts(0) = 1; global_k(1) = 1
! store original U before rotating
!~ ! phase factor ph_g is applied to u_matrix
!~ ! NB: ph_g is applied to u_matrix_opt if (have_disentangled)
!~ if (have_disentangled) then
!~ u0=u_matrix
!~ else
!~ do iw=1,num_wann
!~ u0(iw,:,:)= conjg(ph_g(iw))*u_matrix(iw,:,:)
!~ end do
!~ endif
u0 = u_matrix
! guiding centres are not neede for orthorhombic systems
if (kmesh_info%nntot .eq. 3) wann_control%guiding_centres%enable = .false.
if (wann_control%guiding_centres%enable) then
do n = 1, num_wann
call utility_frac_to_cart(wann_control%guiding_centres%centres(:, n), rguide(:, n), &
real_lattice)
end do
end if
write (stdout, *)
write (stdout, '(1x,a)') '*------------------------------- WANNIERISE ---------------------------------*'
write (stdout, '(1x,a)') '+--------------------------------------------------------------------+<-- CONV'
if (trim(print_output%length_unit) == 'Ang') then
write (stdout, '(1x,a)') '| Iter Delta Spread RMS Gradient Spread (Ang^2) Time |<-- CONV'
else
write (stdout, '(1x,a)') '| Iter Delta Spread RMS Gradient Spread (Bohr^2) Time |<-- CONV'
end if
write (stdout, '(1x,a)') '+--------------------------------------------------------------------+<-- CONV'
write (stdout, *)
irguide = 0
if (wann_control%guiding_centres%enable .and. (wann_control%guiding_centres%num_no_guide_iter .le. 0)) then
call wann_phases(csheet, sheet, rguide, irguide, num_wann, kmesh_info, num_kpts, &
wann_control%use_ss_functional, m_matrix, print_output%timing_level, &
print_output%iprint, timer, num_kpts, global_k, error, comm) ! no plellisation so num_kpts_local = num_kpts
if (allocated(error)) return
irguide = 1
end if
! weight m_matrix first to reduce number of operations
! m_w : weighted real matrix
do nn = 1, kmesh_info%nntot
sqwb = sqrt(kmesh_info%wb(nn))
m_w(:, :, 2*nn - 1) = sqwb*real(m_matrix(1:num_wann, 1:num_wann, nn, 1), dp)
m_w(:, :, 2*nn) = sqwb*aimag(m_matrix(1:num_wann, 1:num_wann, nn, 1))
end do
! calculate initial centers and spread
call wann_omega_gamma(m_w, csheet, sheet, rave, r2ave, rave2, wann_spread, num_wann, &
kmesh_info%nntot, kmesh_info%wbtot, kmesh_info%wb, kmesh_info%bk, &
omega%invariant, ln_tmp, first_pass, &
print_output%timing_level, timer, error, comm)
if (allocated(error)) return
! public variables
omega%total = wann_spread%om_tot
omega%invariant = wann_spread%om_i
omega%tilde = wann_spread%om_d + wann_spread%om_od
! Public array of Wannier centres and spreads
wannier_data%centres = rave
wannier_data%spreads = r2ave - rave2
iter = 0
old_spread%om_tot = 0.0_dp
! print initial state
write (stdout, '(1x,a78)') repeat('-', 78)
write (stdout, '(1x,a)') 'Initial State'
do iw = 1, num_wann
write (stdout, 1000) iw, (rave(ind, iw)*print_output%lenconfac, ind=1, 3), &
(r2ave(iw) - rave2(iw))*print_output%lenconfac**2
end do
write (stdout, 1001) (sum(rave(ind, :))*print_output%lenconfac, ind=1, 3), (sum(r2ave) - sum(rave2))*print_output%lenconfac**2
write (stdout, *)
write (stdout, '(1x,i6,2x,E12.3,19x,F18.10,3x,F8.2,2x,a)') &
iter, (wann_spread%om_tot - old_spread%om_tot)*print_output%lenconfac**2, &
wann_spread%om_tot*print_output%lenconfac**2, io_time(), '<-- CONV'
write (stdout, '(8x,a,F15.7,a,F15.7,a,F15.7,a)') &
'O_D=', wann_spread%om_d*print_output%lenconfac**2, ' O_OD=', wann_spread%om_od*print_output%lenconfac**2, &
' O_TOT=', wann_spread%om_tot*print_output%lenconfac**2, ' <-- SPRD'
write (stdout, '(1x,a78)') repeat('-', 78)
lconverged = .false.
! initialize ur_rot
ur_rot = 0.0_dp
do i = 1, num_wann
ur_rot(i, i) = 1.0_dp
end do
! main iteration loop
do iter = 1, wann_control%num_iter
lprint = .false.
if ((mod(iter, wann_control%num_print_cycles) .eq. 0) .or. (iter .eq. 1) &
.or. (iter .eq. wann_control%num_iter)) lprint = .true.
ldump = .false.
if ((wann_control%num_dump_cycles .gt. 0) .and. &
(mod(iter, wann_control%num_dump_cycles) .eq. 0)) ldump = .true.
if (lprint .and. print_output%iprint > 0) write (stdout, '(1x,a,i6)') 'Cycle: ', iter
if (wann_control%guiding_centres%enable .and. &
(iter .gt. wann_control%guiding_centres%num_no_guide_iter) &
.and. (mod(iter, wann_control%guiding_centres%num_guide_cycles) .eq. 0)) then
call wann_phases(csheet, sheet, rguide, irguide, num_wann, kmesh_info, num_kpts, &
wann_control%use_ss_functional, m_matrix, print_output%timing_level, &
print_output%iprint, timer, num_kpts, global_k, error, comm, m_w) ! num_kpts_loc == num_kpts here
if (allocated(error)) return
irguide = 1
end if
call internal_new_u_and_m_gamma(m_w, ur_rot, tnntot, num_wann, print_output%timing_level, &
timer)
call wann_spread_copy(wann_spread, old_spread)
! calculate the new centers and spread
call wann_omega_gamma(m_w, csheet, sheet, rave, r2ave, rave2, wann_spread, num_wann, &
kmesh_info%nntot, kmesh_info%wbtot, kmesh_info%wb, kmesh_info%bk, &
omega%invariant, ln_tmp, first_pass, print_output%timing_level, &
timer, error, comm)
if (allocated(error)) return
! print the new centers and spreads
if (lprint) then
do iw = 1, num_wann
write (stdout, 1000) iw, (rave(ind, iw)*print_output%lenconfac, ind=1, 3), &
(r2ave(iw) - rave2(iw))*print_output%lenconfac**2
end do
write (stdout, 1001) (sum(rave(ind, :))*print_output%lenconfac, ind=1, 3), &
(sum(r2ave) - sum(rave2))*print_output%lenconfac**2
write (stdout, *)
write (stdout, '(1x,i6,2x,E12.3,19x,F18.10,3x,F8.2,2x,a)') &
iter, (wann_spread%om_tot - old_spread%om_tot)*print_output%lenconfac**2, &
wann_spread%om_tot*print_output%lenconfac**2, io_time(), '<-- CONV'
write (stdout, '(8x,a,F15.7,a,F15.7,a,F15.7,a)') &
'O_D=', wann_spread%om_d*print_output%lenconfac**2, &
' O_OD=', wann_spread%om_od*print_output%lenconfac**2, &
' O_TOT=', wann_spread%om_tot*print_output%lenconfac**2, ' <-- SPRD'
write (stdout, '(1x,a,E15.7,a,E15.7,a,E15.7,a)') &
'Delta: O_D=', (wann_spread%om_d - old_spread%om_d)*print_output%lenconfac**2, &
' O_OD=', (wann_spread%om_od - old_spread%om_od)*print_output%lenconfac**2, &
' O_TOT=', (wann_spread%om_tot - old_spread%om_tot)*print_output%lenconfac**2, ' <-- DLTA'
write (stdout, '(1x,a78)') repeat('-', 78)
end if
! Public array of Wannier centres and spreads
wannier_data%centres = rave
wannier_data%spreads = r2ave - rave2
! Public variables
omega%total = wann_spread%om_tot
omega%tilde = wann_spread%om_d + wann_spread%om_od
! (Jerome Jackson) Removing checkpoint from WF optimisation loop because benefit is limited
! if (ldump) then
! uc_rot(:, :) = cmplx(ur_rot(:, :), 0.0_dp, dp)
! call utility_zgemm(u_matrix, u0, 'N', uc_rot, 'N', num_wann)
! call w90_wannier90_readwrite_write_chkpt('postdis', exclude_bands, wannier_data, &
! kmesh_info, kpt_latt, num_kpts, dis_manifold, &
! num_bands, num_wann, u_matrix, u_matrix_opt, &
! m_matrix, mp_grid, real_lattice, omega%invariant, &
! have_disentangled, stdout, seedname)
! endif
if (wann_control%conv_window .gt. 1) then
call internal_test_convergence(old_spread, wann_spread, history, save_spread, iter, &
conv_count, noise_count, lconverged, lrandom, lfirst, &
wann_control, error, comm)
if (allocated(error)) return
end if
if (lconverged) then
write (stdout, '(/13x,a,es10.3,a,i2,a)') &
'<<< Delta <', wann_control%conv_tol, &
' over ', wann_control%conv_window, ' iterations >>>'
write (stdout, '(13x,a/)') '<<< Wannierisation convergence criteria satisfied >>>'
exit
end if
end do
! end of the minimization loop
! update M
do nn = 1, kmesh_info%nntot
sqwb = 1.0_dp/sqrt(kmesh_info%wb(nn))
m_matrix(1:num_wann, 1:num_wann, nn, 1) = sqwb*cmplx(m_w(:, :, 2*nn - 1), m_w(:, :, 2*nn), dp)
end do
! update U
uc_rot(:, :) = cmplx(ur_rot(:, :), 0.0_dp, dp)
call utility_zgemm(u_matrix, u0, 'N', uc_rot, 'N', num_wann)
write (stdout, '(1x,a)') 'Final State'
do iw = 1, num_wann
write (stdout, 1000) iw, (rave(ind, iw)*print_output%lenconfac, ind=1, 3), &
(r2ave(iw) - rave2(iw))*print_output%lenconfac**2
end do
write (stdout, 1001) (sum(rave(ind, :))*print_output%lenconfac, ind=1, 3), &
(sum(r2ave) - sum(rave2))*print_output%lenconfac**2
write (stdout, *)
write (stdout, '(3x,a21,a,f15.9)') ' Spreads ('//trim(print_output%length_unit)//'^2)', &
' Omega I = ', wann_spread%om_i*print_output%lenconfac**2
write (stdout, '(3x,a,f15.9)') ' ================ Omega D = ', &
wann_spread%om_d*print_output%lenconfac**2
write (stdout, '(3x,a,f15.9)') ' Omega OD = ', &
wann_spread%om_od*print_output%lenconfac**2
write (stdout, '(3x,a21,a,f15.9)') 'Final Spread ('//trim(print_output%length_unit)//'^2)', &
' Omega Total = ', wann_spread%om_tot*print_output%lenconfac**2
write (stdout, '(1x,a78)') repeat('-', 78)
if (wann_control%guiding_centres%enable) then
call wann_phases(csheet, sheet, rguide, irguide, num_wann, kmesh_info, num_kpts, &
wann_control%use_ss_functional, m_matrix, print_output%timing_level, &
print_output%iprint, timer, num_kpts, global_k, error, comm) ! num_kpts_loc == num_kpts here
if (allocated(error)) return
end if
! unitarity is checked
call wann_check_unitarity(num_kpts, num_wann, u_matrix, print_output%timing_level, &
print_output%iprint, stdout, timer, error, comm)
if (allocated(error)) return
! deallocate sub vars not passed into other subs
deallocate (cz, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating cz in wann_main_gamma', comm)
return
end if
deallocate (ur_rot, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating ur_rot in wann_main_gamma', comm)
return
end if
deallocate (uc_rot, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating uc_rot in wann_main_gamma', comm)
return
end if
deallocate (u0, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating u0 in wann_main_gamma', comm)
return
end if
! deallocate sub vars passed into other subs
deallocate (rguide, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating rguide in wann_main_gamma', comm)
return
end if
deallocate (rave2, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating rave2 in wann_main_gamma', comm)
return
end if
deallocate (rave, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating rave in wann_main_gamma', comm)
return
end if
deallocate (sheet, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating sheet in wann_main_gamma', comm)
return
end if
deallocate (csheet, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating csheet in wann_main_gamma', comm)
return
end if
deallocate (m_w, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating m_w in wann_main_gamma', comm)
return
end if
! deallocate module data
deallocate (ln_tmp, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating ln_tmp in wann_main_gamma', comm)
return
end if
deallocate (history, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error deallocating history in wann_main_gamma', comm)
return
end if
if (print_output%timing_level > 0) call io_stopwatch_stop('wann: main_gamma', timer)
return
1000 format(2x, 'WF centre and spread', i5, 2x, '(', f10.6, ',', f10.6, ',', f10.6, ' )', f15.8)
1001 format(2x, 'Sum of centres and spreads', 1x, '(', f10.6, ',', f10.6, ',', f10.6, ' )', f15.8)
contains
!================================================!
subroutine internal_new_u_and_m_gamma(m_w, ur_rot, tnntot, num_wann, timing_level, timer)
!================================================!
use w90_constants, only: pi, eps10
use w90_io, only: io_stopwatch_start, io_stopwatch_stop
use w90_types, only: timer_list_type
implicit none
! arguments
real(kind=dp), intent(inout) :: m_w(:, :, :)
real(kind=dp), intent(inout) :: ur_rot(:, :)
integer, intent(in) :: tnntot
integer, intent(in) :: num_wann
integer, intent(in) :: timing_level
type(timer_list_type), intent(inout) :: timer
! local variables
real(kind=dp) :: theta, twotheta
real(kind=dp) :: a11, a12, a21, a22
real(kind=dp) :: cc, ss, rtmp1, rtmp2
real(kind=dp), parameter :: pifour = 0.25_dp*pi
integer :: nn, nw1, nw2, nw3
if (timing_level > 1) call io_stopwatch_start('wann: main_gamma: new_u_and_m_gamma', timer)
loop_nw1: do nw1 = 1, num_wann
loop_nw2: do nw2 = nw1 + 1, num_wann
a11 = 0.0_dp; a12 = 0.0_dp; a22 = 0.0_dp
do nn = 1, tnntot
a11 = a11 + (m_w(nw1, nw1, nn) - m_w(nw2, nw2, nn))**2
a12 = a12 + m_w(nw1, nw2, nn)*(m_w(nw1, nw1, nn) - m_w(nw2, nw2, nn))
a22 = a22 + m_w(nw1, nw2, nn)**2
end do
a12 = 2.0_dp*a12
a22 = 4.0_dp*a22
a21 = a22 - a11
if (abs(a12) .gt. eps10) then
twotheta = 0.5_dp*(a21 + sqrt(a21**2 + 4.0_dp*a12**2))/a12
theta = 0.5_dp*atan(twotheta)
elseif (a21 .lt. eps10) then
theta = 0.0_dp
else
theta = pifour
end if
cc = cos(theta)
ss = sin(theta)
! update M
do nn = 1, tnntot
! MR
do nw3 = 1, num_wann
rtmp1 = m_w(nw3, nw1, nn)*cc + m_w(nw3, nw2, nn)*ss
rtmp2 = -m_w(nw3, nw1, nn)*ss + m_w(nw3, nw2, nn)*cc
m_w(nw3, nw1, nn) = rtmp1
m_w(nw3, nw2, nn) = rtmp2
end do
! R^+ M R
do nw3 = 1, num_wann
rtmp1 = cc*m_w(nw1, nw3, nn) + ss*m_w(nw2, nw3, nn)
rtmp2 = -ss*m_w(nw1, nw3, nn) + cc*m_w(nw2, nw3, nn)
m_w(nw1, nw3, nn) = rtmp1
m_w(nw2, nw3, nn) = rtmp2
end do
end do
! update U : U=UR
do nw3 = 1, num_wann
rtmp1 = ur_rot(nw3, nw1)*cc + ur_rot(nw3, nw2)*ss
rtmp2 = -ur_rot(nw3, nw1)*ss + ur_rot(nw3, nw2)*cc
ur_rot(nw3, nw1) = rtmp1
ur_rot(nw3, nw2) = rtmp2
end do
end do loop_nw2
end do loop_nw1
if (timing_level > 1) call io_stopwatch_stop('wann: main_gamma: new_u_and_m_gamma', timer)
return
end subroutine internal_new_u_and_m_gamma
end subroutine wann_main_gamma
