subroutine wann_main_gamma
!==================================================================!
! !
! 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: stdout, io_error, io_time, io_stopwatch
use w90_parameters, only: num_wann, num_iter, wb, &
nntot, u_matrix, m_matrix, num_kpts, iprint, &
num_print_cycles, num_dump_cycles, omega_invariant, &
param_write_chkpt, length_unit, lenconfac, &
proj_site, real_lattice, write_r2mn, guiding_centres, &
num_guide_cycles, num_no_guide_iter, timing_level, &
write_proj, have_disentangled, conv_tol, conv_window, &
wannier_centres, write_xyz, wannier_spreads, omega_total, &
omega_tilde, write_vdw_data
use w90_utility, only: utility_frac_to_cart, utility_zgemm
implicit none
type(localisation_vars) :: old_spread
type(localisation_vars) :: wann_spread
! 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
if (timing_level > 0) call io_stopwatch('wann: main_gamma', 1)
first_pass = .true.
! Allocate stuff
allocate (history(conv_window), stat=ierr)
if (ierr /= 0) call io_error('Error allocating history in wann_main_gamma')
!~ if (.not.allocated(ph_g)) then
!~ allocate( ph_g(num_wann),stat=ierr )
!~ if (ierr/=0) call io_error('Error in allocating ph_g in wann_main_gamma')
!~ ph_g = cmplx_1
!~ endif
! module data
allocate (rnkb(num_wann, nntot, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating rnkb in wann_main_gamma')
allocate (ln_tmp(num_wann, nntot, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ln_tmp in wann_main_gamma')
rnkb = 0.0_dp
tnntot = 2*nntot
! sub vars passed into other subs
allocate (m_w(num_wann, num_wann, tnntot), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating m_w in wann_main_gamma')
allocate (csheet(num_wann, nntot, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating csheet in wann_main_gamma')
allocate (sheet(num_wann, nntot, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating sheet in wann_main_gamma')
allocate (rave(3, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating rave in wann_main_gamma')
allocate (r2ave(num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating r2ave in wann_main_gamma')
allocate (rave2(num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating rave2 in wann_main_gamma')
allocate (rguide(3, num_wann))
if (ierr /= 0) call io_error('Error in allocating rguide in wann_main_gamma')
csheet = cmplx_1
sheet = 0.0_dp; rave = 0.0_dp; r2ave = 0.0_dp; rave2 = 0.0_dp; rguide = 0.0_dp
! sub vars not passed into other subs
allocate (u0(num_wann, num_wann, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating u0 in wann_main_gamma')
allocate (uc_rot(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating uc_rot in wann_main_gamma')
allocate (ur_rot(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating ur_rot in wann_main_gamma')
allocate (cz(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating cz in wann_main_gamma')
cz = cmplx_0
! Set up the MPI arrays for a serial run.
allocate (counts(0:0), displs(0:0), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating counts and displs in wann_main_gamma')
counts(0) = 1; displs(0) = 0
! 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
!~ lguide = .false.
! guiding centres are not neede for orthorhombic systems
if (nntot .eq. 3) guiding_centres = .false.
if (guiding_centres) then
! initialise rguide to projection centres (Cartesians in units of Ang)
!~ if ( use_bloch_phases) then
!~ lguide = .true.
!~ else
do n = 1, num_wann
call utility_frac_to_cart(proj_site(:, n), rguide(:, n), real_lattice)
enddo
!~ endif
endif
write (stdout, *)
write (stdout, '(1x,a)') '*------------------------------- WANNIERISE ---------------------------------*'
write (stdout, '(1x,a)') '+--------------------------------------------------------------------+<-- CONV'
if (lenconfac .eq. 1.0_dp) 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'
endif
write (stdout, '(1x,a)') '+--------------------------------------------------------------------+<-- CONV'
write (stdout, *)
irguide = 0
!~ if (guiding_centres.and.(num_no_guide_iter.le.0)) then
!~ if (nntot.gt.3) call wann_phases(csheet,sheet,rguide,irguide)
!~ irguide=1
!~ endif
if (guiding_centres .and. (num_no_guide_iter .le. 0)) then
call wann_phases(csheet, sheet, rguide, irguide)
irguide = 1
endif
! weight m_matrix first to reduce number of operations
! m_w : weighted real matrix
do nn = 1, nntot
sqwb = sqrt(wb(nn))
m_w(:, :, 2*nn - 1) = sqwb*real(m_matrix(:, :, nn, 1), dp)
m_w(:, :, 2*nn) = sqwb*aimag(m_matrix(:, :, nn, 1))
end do
! calculate initial centers and spread
call wann_omega_gamma(m_w, csheet, sheet, rave, r2ave, rave2, wann_spread)
! 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_centres = rave
wannier_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)*lenconfac, ind=1, 3), &
(r2ave(iw) - rave2(iw))*lenconfac**2
end do
write (stdout, 1001) (sum(rave(ind, :))*lenconfac, ind=1, 3), (sum(r2ave) - sum(rave2))*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)*lenconfac**2, &
wann_spread%om_tot*lenconfac**2, io_time(), '<-- CONV'
write (stdout, '(8x,a,F15.7,a,F15.7,a,F15.7,a)') &
'O_D=', wann_spread%om_d*lenconfac**2, ' O_OD=', wann_spread%om_od*lenconfac**2, &
' O_TOT=', wann_spread%om_tot*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, num_iter
lprint = .false.
if ((mod(iter, num_print_cycles) .eq. 0) .or. (iter .eq. 1) &
.or. (iter .eq. num_iter)) lprint = .true.
ldump = .false.
if ((num_dump_cycles .gt. 0) .and. (mod(iter, num_dump_cycles) .eq. 0)) ldump = .true.
if (lprint .and. on_root) write (stdout, '(1x,a,i6)') 'Cycle: ', iter
!~ ! initialize rguide as rave for use_bloch_phases
!~ if ( (iter.gt.num_no_guide_iter) .and. lguide ) then
!~ rguide(:,:) = rave(:,:)
!~ lguide = .false.
!~ endif
!~ if ( guiding_centres.and.(iter.gt.num_no_guide_iter) &
!~ .and.(mod(iter,num_guide_cycles).eq.0) ) then
!~ if(nntot.gt.3) call wann_phases(csheet,sheet,rguide,irguide)
!~ irguide=1
!~ endif
if (guiding_centres .and. (iter .gt. num_no_guide_iter) &
.and. (mod(iter, num_guide_cycles) .eq. 0)) then
call wann_phases(csheet, sheet, rguide, irguide, m_w)
irguide = 1
endif
call internal_new_u_and_m_gamma()
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)
! print the new centers and spreads
if (lprint) then
do iw = 1, num_wann
write (stdout, 1000) iw, (rave(ind, iw)*lenconfac, ind=1, 3), &
(r2ave(iw) - rave2(iw))*lenconfac**2
end do
write (stdout, 1001) (sum(rave(ind, :))*lenconfac, ind=1, 3), &
(sum(r2ave) - sum(rave2))*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)*lenconfac**2, &
wann_spread%om_tot*lenconfac**2, io_time(), '<-- CONV'
write (stdout, '(8x,a,F15.7,a,F15.7,a,F15.7,a)') &
'O_D=', wann_spread%om_d*lenconfac**2, &
' O_OD=', wann_spread%om_od*lenconfac**2, &
' O_TOT=', wann_spread%om_tot*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)*lenconfac**2, &
' O_OD=', (wann_spread%om_od - old_spread%om_od)*lenconfac**2, &
' O_TOT=', (wann_spread%om_tot - old_spread%om_tot)*lenconfac**2, ' <-- DLTA'
write (stdout, '(1x,a78)') repeat('-', 78)
end if
! Public array of Wannier centres and spreads
wannier_centres = rave
wannier_spreads = r2ave - rave2
! Public variables
omega_total = wann_spread%om_tot
omega_tilde = wann_spread%om_d + wann_spread%om_od
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 param_write_chkpt('postdis')
endif
if (conv_window .gt. 1) call internal_test_convergence_gamma()
if (lconverged) then
write (stdout, '(/13x,a,es10.3,a,i2,a)') &
'<<< Delta <', conv_tol, &
' over ', conv_window, ' iterations >>>'
write (stdout, '(13x,a/)') '<<< Wannierisation convergence criteria satisfied >>>'
exit
endif
enddo
! end of the minimization loop
! update M
do nn = 1, nntot
sqwb = 1.0_dp/sqrt(wb(nn))
m_matrix(:, :, 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)*lenconfac, ind=1, 3), &
(r2ave(iw) - rave2(iw))*lenconfac**2
end do
write (stdout, 1001) (sum(rave(ind, :))*lenconfac, ind=1, 3), &
(sum(r2ave) - sum(rave2))*lenconfac**2
write (stdout, *)
write (stdout, '(3x,a21,a,f15.9)') ' Spreads ('//trim(length_unit)//'^2)', &
' Omega I = ', wann_spread%om_i*lenconfac**2
write (stdout, '(3x,a,f15.9)') ' ================ Omega D = ', &
wann_spread%om_d*lenconfac**2
write (stdout, '(3x,a,f15.9)') ' Omega OD = ', &
wann_spread%om_od*lenconfac**2
write (stdout, '(3x,a21,a,f15.9)') 'Final Spread ('//trim(length_unit)//'^2)', &
' Omega Total = ', wann_spread%om_tot*lenconfac**2
write (stdout, '(1x,a78)') repeat('-', 78)
if (write_xyz .and. on_root) call wann_write_xyz()
if (guiding_centres) call wann_phases(csheet, sheet, rguide, irguide)
! unitarity is checked
!~ call internal_check_unitarity()
call wann_check_unitarity()
! write extra info regarding omega_invariant
!~ if (iprint>2) call internal_svd_omega_i()
if (iprint > 2) call wann_svd_omega_i()
! write matrix elements <m|r^2|n> to file
!~ if (write_r2mn) call internal_write_r2mn()
if (write_r2mn) call wann_write_r2mn()
! calculate and write projection of WFs on original bands in outer window
if (have_disentangled .and. write_proj) call wann_calc_projection()
! aam: write data required for vdW utility
if (write_vdw_data) call wann_write_vdw_data()
! deallocate sub vars not passed into other subs
deallocate (cz, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating cz in wann_main_gamma')
deallocate (ur_rot, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating ur_rot in wann_main_gamma')
deallocate (uc_rot, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating uc_rot in wann_main_gamma')
deallocate (u0, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating u0 in wann_main_gamma')
! deallocate sub vars passed into other subs
deallocate (rguide, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating rguide in wann_main_gamma')
deallocate (rave2, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating rave2 in wann_main_gamma')
deallocate (rave, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating rave in wann_main_gamma')
deallocate (sheet, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating sheet in wann_main_gamma')
deallocate (csheet, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating csheet in wann_main_gamma')
deallocate (m_w, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating m_w in wann_main_gamma')
! deallocate module data
deallocate (ln_tmp, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating ln_tmp in wann_main_gamma')
deallocate (rnkb, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating rnkb in wann_main_gamma')
deallocate (history, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating history in wann_main_gamma')
if (timing_level > 0) call io_stopwatch('wann: main_gamma', 2)
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()
!===============================================!
use w90_constants, only: pi, eps10
implicit none
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('wann: main_gamma: new_u_and_m_gamma', 1)
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
endif
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('wann: main_gamma: new_u_and_m_gamma', 2)
return
end subroutine internal_new_u_and_m_gamma
!===============================================!
subroutine internal_test_convergence_gamma()
!===============================================!
! !
! Determine whether minimisation of non-gauge- !
! invariant spread is converged !
! !
!===============================================!
implicit none
real(kind=dp) :: delta_omega
integer :: j, ierr
real(kind=dp), allocatable :: temp_hist(:)
allocate (temp_hist(conv_window), stat=ierr)
if (ierr /= 0) call io_error('Error allocating temp_hist in wann_main')
delta_omega = wann_spread%om_tot - old_spread%om_tot
if (iter .le. conv_window) then
history(iter) = delta_omega
else
temp_hist = eoshift(history, 1, delta_omega)
history = temp_hist
endif
lconverged = .false.
if (iter .ge. conv_window) then
!~ write(stdout,*) (history(j),j=1,conv_window)
do j = 1, conv_window
if (abs(history(j)) .gt. conv_tol) exit
lconverged = .true.
enddo
endif
deallocate (temp_hist, stat=ierr)
if (ierr /= 0) call io_error('Error deallocating temp_hist in wann_main_gamma')
return
end subroutine internal_test_convergence_gamma
!~ !========================================!
!~ subroutine internal_check_unitarity()
!~ !========================================!
!~
!~ implicit none
!~
!~ integer :: nkp,i,j,m
!~ complex(kind=dp) :: ctmp1,ctmp2
!~
!~ if (timing_level>1) call io_stopwatch('wann: main: check_unitarity',1)
!~
!~ do nkp = 1, num_kpts
!~ do i = 1, num_wann
!~ do j = 1, num_wann
!~ ctmp1 = cmplx_0
!~ ctmp2 = cmplx_0
!~ do m = 1, num_wann
!~ ctmp1 = ctmp1 + u_matrix (i, m, nkp) * conjg (u_matrix (j, m, nkp) )
!~ ctmp2 = ctmp2 + u_matrix (m, j, nkp) * conjg (u_matrix (m, i, nkp) )
!~ enddo
!~ if ( (i.eq.j) .and. (abs (ctmp1 - cmplx_1 ) .gt. eps5) ) &
!~ then
!~ write ( stdout , * ) ' ERROR: unitariety of final U', nkp, i, j, &
!~ ctmp1
!~ call io_error('wann_main: unitariety error 1')
!~ endif
!~ if ( (i.eq.j) .and. (abs (ctmp2 - cmplx_1 ) .gt. eps5) ) &
!~ then
!~ write ( stdout , * ) ' ERROR: unitariety of final U', nkp, i, j, &
!~ ctmp2
!~ call io_error('wann_main: unitariety error 2')
!~ endif
!~ if ( (i.ne.j) .and. (abs (ctmp1) .gt. eps5) ) then
!~ write ( stdout , * ) ' ERROR: unitariety of final U', nkp, i, j, &
!~ ctmp1
!~ call io_error('wann_main: unitariety error 3')
!~ endif
!~ if ( (i.ne.j) .and. (abs (ctmp2) .gt. eps5) ) then
!~ write ( stdout , * ) ' ERROR: unitariety of final U', nkp, i, j, &
!~ ctmp2
!~ call io_error('wann_main: unitariety error 4')
!~ endif
!~ enddo
!~ enddo
!~ enddo
!~
!~ if (timing_level>1) call io_stopwatch('wann: main: check_unitarity',2)
!~
!~ return
!~
!~ end subroutine internal_check_unitarity
!~ !========================================!
!~ subroutine internal_write_r2mn()
!~ !========================================!
!~ ! !
!~ ! Write seedname.r2mn file !
!~ ! !
!~ !========================================!
!~ use w90_io, only: seedname,io_file_unit,io_error
!~
!~ implicit none
!~
!~ integer :: r2mnunit,nw1,nw2,nkp,nn
!~ real(kind=dp) :: r2ave_mn,delta
!~
!~ ! note that here I use formulas analogue to Eq. 23, and not to the
!~ ! shift-invariant Eq. 32 .
!~ r2mnunit=io_file_unit()
!~ open(r2mnunit,file=trim(seedname)//'.r2mn',form='formatted',err=158)
!~ do nw1 = 1, num_wann
!~ do nw2 = 1, num_wann
!~ r2ave_mn = 0.0_dp
!~ delta = 0.0_dp
!~ if (nw1.eq.nw2) delta = 1.0_dp
!~ do nkp = 1, num_kpts
!~ do nn = 1, nntot
!~ r2ave_mn = r2ave_mn + wb(nn) * &
!~ ! [GP-begin, Apr13, 2012: corrected sign inside "real"]
!~ ( 2.0_dp * delta - real(m_matrix(nw1,nw2,nn,nkp) + &
!~ conjg(m_matrix(nw2,nw1,nn,nkp)),kind=dp) )
!~ enddo
!~ enddo
!~ r2ave_mn = r2ave_mn / real(num_kpts,dp)
!~ write (r2mnunit, '(2i6,f20.12)') nw1, nw2, r2ave_mn
!~ enddo
!~ enddo
!~ close(r2mnunit)
!~
!~ return
!~
!~158 call io_error('Error opening file '//trim(seedname)//'.r2mn in wann_main')
!~
!~ end subroutine internal_write_r2mn
!~ !========================================!
!~ subroutine internal_svd_omega_i()
!~ !========================================!
!~
!~ implicit none
!~
!~ complex(kind=dp), allocatable :: cv1(:,:),cv2(:,:)
!~ complex(kind=dp), allocatable :: cw1(:),cw2(:)
!~ complex(kind=dp), allocatable :: cpad1 (:)
!~ real(kind=dp), allocatable :: singvd (:)
!~
!~ integer :: nkp,nn,nb,na,ind
!~ real(kind=dp) :: omt1,omt2,omt3
!~
!~ if (timing_level>1) call io_stopwatch('wann: main: svd_omega_i',1)
!~
!~ allocate( cw1 (10 * num_wann),stat=ierr )
!~ if (ierr/=0) call io_error('Error in allocating cw1 in wann_main')
!~ allocate( cw2 (10 * num_wann),stat=ierr )
!~ if (ierr/=0) call io_error('Error in allocating cw2 in wann_main')
!~ allocate( cv1 (num_wann, num_wann),stat=ierr )
!~ if (ierr/=0) call io_error('Error in allocating cv1 in wann_main')
!~ allocate( cv2 (num_wann, num_wann),stat=ierr )
!~ if (ierr/=0) call io_error('Error in allocating cv2 in wann_main')
!~ allocate( singvd (num_wann),stat=ierr )
!~ if (ierr/=0) call io_error('Error in allocating singvd in wann_main')
!~ allocate( cpad1 (num_wann * num_wann),stat=ierr )
!~ if (ierr/=0) call io_error('Error in allocating cpad1 in wann_main')
!~
!~ cw1=cmplx_0; cw2=cmplx_0; cv1=cmplx_0; cv2=cmplx_0; cpad1=cmplx_0
!~ singvd=0.0_dp
!~
!~ ! singular value decomposition
!~ omt1 = 0.0_dp ; omt2 = 0.0_dp ; omt3 = 0.0_dp
!~ do nkp = 1, num_kpts
!~ do nn = 1, nntot
!~ ind = 1
!~ do nb = 1, num_wann
!~ do na = 1, num_wann
!~ cpad1 (ind) = m_matrix (na, nb, nn, nkp)
!~ ind = ind+1
!~ enddo
!~ enddo
!~ call zgesvd ('A', 'A', num_wann, num_wann, cpad1, num_wann, singvd, cv1, &
!~ num_wann, cv2, num_wann, cw1, 10 * num_wann, cw2, info)
!~ if (info.ne.0) then
!~ call io_error('ERROR: Singular value decomp. zgesvd failed')
!~ endif
!~
!~ do nb = 1, num_wann
!~ omt1 = omt1 + wb(nn) * (1.0_dp - singvd (nb) **2)
!~ omt2 = omt2 - wb(nn) * (2.0_dp * log (singvd (nb) ) )
!~ omt3 = omt3 + wb(nn) * (acos (singvd (nb) ) **2)
!~ enddo
!~ enddo
!~ enddo
!~ omt1 = omt1 / real(num_kpts,dp)
!~ omt2 = omt2 / real(num_kpts,dp)
!~ omt3 = omt3 / real(num_kpts,dp)
!~ write ( stdout , * ) ' '
!~ write(stdout,'(2x,a,f15.9,1x,a)') 'Omega Invariant: 1-s^2 = ',&
!~ omt1*lenconfac**2,'('//trim(length_unit)//'^2)'
!~ write(stdout,'(2x,a,f15.9,1x,a)') ' -2log s = ',&
!~ omt2*lenconfac**2,'('//trim(length_unit)//'^2)'
!~ write(stdout,'(2x,a,f15.9,1x,a)') ' acos^2 = ',&
!~ omt3*lenconfac**2,'('//trim(length_unit)//'^2)'
!~
!~ deallocate(cpad1,stat=ierr)
!~ if (ierr/=0) call io_error('Error in deallocating cpad1 in wann_main')
!~ deallocate(singvd,stat=ierr)
!~ if (ierr/=0) call io_error('Error in deallocating singvd in wann_main')
!~ deallocate(cv2,stat=ierr)
!~ if (ierr/=0) call io_error('Error in deallocating cv2 in wann_main')
!~ deallocate(cv1,stat=ierr)
!~ if (ierr/=0) call io_error('Error in deallocating cv1 in wann_main')
!~ deallocate(cw2,stat=ierr)
!~ if (ierr/=0) call io_error('Error in deallocating cw2 in wann_main')
!~ deallocate(cw1,stat=ierr)
!~ if (ierr/=0) call io_error('Error in deallocating cw1 in wann_main')
!~
!~ if (timing_level>1) call io_stopwatch('wann: main: svd_omega_i',2)
!~
!~ return
!~
!~ end subroutine internal_svd_omega_i
end subroutine wann_main_gamma