Read the Mmn and Amn from files Note: one needs to call overlap_allocate first!
subroutine overlap_read()
!%%%%%%%%%%%%%%%%%%%%%
!! Read the Mmn and Amn from files
!! Note: one needs to call overlap_allocate first!
use w90_parameters, only: num_bands, num_wann, num_kpts, nntot, nncell, nnlist, &
num_proj, lselproj, proj2wann_map, &
devel_flag, u_matrix, m_matrix, a_matrix, timing_level, &
m_matrix_orig, u_matrix_opt, cp_pp, use_bloch_phases, gamma_only, & ![ysl]
m_matrix_local, m_matrix_orig_local, lhasproj
use w90_io, only: io_file_unit, io_error, seedname, io_stopwatch
use w90_comms, only: my_node_id, num_nodes, &
comms_array_split, comms_scatterv
implicit none
integer :: nkp, nkp2, inn, nn, n, m, i, j
integer :: mmn_in, amn_in, num_mmn, num_amn
integer :: nnl, nnm, nnn, ncount
integer :: nb_tmp, nkp_tmp, nntot_tmp, np_tmp, ierr
real(kind=dp) :: m_real, m_imag, a_real, a_imag, mu_tmp, sigma_tmp
complex(kind=dp), allocatable :: mmn_tmp(:, :)
character(len=50) :: dummy
logical :: nn_found
! Needed to split an array on different nodes
integer, dimension(0:num_nodes - 1) :: counts
integer, dimension(0:num_nodes - 1) :: displs
if (timing_level > 0) call io_stopwatch('overlap: read', 1)
call comms_array_split(num_kpts, counts, displs)
if (disentanglement) then
if (on_root) then
m_matrix_orig = cmplx_0
endif
m_matrix_orig_local = cmplx_0
a_matrix = cmplx_0
u_matrix_opt = cmplx_0
endif
if (on_root) then
! Read M_matrix_orig from file
mmn_in = io_file_unit()
open (unit=mmn_in, file=trim(seedname)//'.mmn', &
form='formatted', status='old', action='read', err=101)
if (on_root) write (stdout, '(/a)', advance='no') ' Reading overlaps from '//trim(seedname)//'.mmn : '
! Read the comment line
read (mmn_in, '(a)', err=103, end=103) dummy
if (on_root) write (stdout, '(a)') trim(dummy)
! Read the number of bands, k-points and nearest neighbours
read (mmn_in, *, err=103, end=103) nb_tmp, nkp_tmp, nntot_tmp
! Checks
if (nb_tmp .ne. num_bands) &
call io_error(trim(seedname)//'.mmn has not the right number of bands')
if (nkp_tmp .ne. num_kpts) &
call io_error(trim(seedname)//'.mmn has not the right number of k-points')
if (nntot_tmp .ne. nntot) &
call io_error(trim(seedname)//'.mmn has not the right number of nearest neighbours')
! Read the overlaps
num_mmn = num_kpts*nntot
allocate (mmn_tmp(num_bands, num_bands), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating mmn_tmp in overlap_read')
do ncount = 1, num_mmn
read (mmn_in, *, err=103, end=103) nkp, nkp2, nnl, nnm, nnn
do n = 1, num_bands
do m = 1, num_bands
read (mmn_in, *, err=103, end=103) m_real, m_imag
mmn_tmp(m, n) = cmplx(m_real, m_imag, kind=dp)
enddo
enddo
nn = 0
nn_found = .false.
do inn = 1, nntot
if ((nkp2 .eq. nnlist(nkp, inn)) .and. &
(nnl .eq. nncell(1, nkp, inn)) .and. &
(nnm .eq. nncell(2, nkp, inn)) .and. &
(nnn .eq. nncell(3, nkp, inn))) then
if (.not. nn_found) then
nn_found = .true.
nn = inn
else
call io_error('Error reading '//trim(seedname)// &
'.mmn. More than one matching nearest neighbour found')
endif
endif
end do
if (nn .eq. 0) then
if (on_root) write (stdout, '(/a,i8,2i5,i4,2x,3i3)') &
' Error reading '//trim(seedname)//'.mmn:', ncount, nkp, nkp2, nn, nnl, nnm, nnn
call io_error('Neighbour not found')
end if
if (disentanglement) then
m_matrix_orig(:, :, nn, nkp) = mmn_tmp(:, :)
else
! disentanglement=.false. means numbands=numwann, so no the dimensions are the same
m_matrix(:, :, nn, nkp) = mmn_tmp(:, :)
end if
end do
deallocate (mmn_tmp, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating mmn_tmp in overlap_read')
close (mmn_in)
endif
if (disentanglement) then
! call comms_bcast(m_matrix_orig(1,1,1,1),num_bands*num_bands*nntot*num_kpts)
call comms_scatterv(m_matrix_orig_local, num_bands*num_bands*nntot*counts(my_node_id), &
m_matrix_orig, num_bands*num_bands*nntot*counts, num_bands*num_bands*nntot*displs)
else
! call comms_bcast(m_matrix(1,1,1,1),num_wann*num_wann*nntot*num_kpts)
call comms_scatterv(m_matrix_local, num_wann*num_wann*nntot*counts(my_node_id), &
m_matrix, num_wann*num_wann*nntot*counts, num_wann*num_wann*nntot*displs)
endif
if (.not. use_bloch_phases) then
if (on_root) then
! Read A_matrix from file wannier.amn
amn_in = io_file_unit()
open (unit=amn_in, file=trim(seedname)//'.amn', form='formatted', status='old', err=102)
if (on_root) write (stdout, '(/a)', advance='no') ' Reading projections from '//trim(seedname)//'.amn : '
! Read the comment line
read (amn_in, '(a)', err=104, end=104) dummy
if (on_root) write (stdout, '(a)') trim(dummy)
! Read the number of bands, k-points and wannier functions
read (amn_in, *, err=104, end=104) nb_tmp, nkp_tmp, np_tmp
! Checks
if (nb_tmp .ne. num_bands) &
call io_error(trim(seedname)//'.amn has not the right number of bands')
if (nkp_tmp .ne. num_kpts) &
call io_error(trim(seedname)//'.amn has not the right number of k-points')
if (np_tmp .ne. num_proj) &
call io_error(trim(seedname)//'.amn has not the right number of projections')
if (num_proj > num_wann .and. .not. lselproj) &
call io_error(trim(seedname)//'.amn has too many projections to be used without selecting a subset')
! Read the projections
num_amn = num_bands*num_proj*num_kpts
if (disentanglement) then
do ncount = 1, num_amn
read (amn_in, *, err=104, end=104) m, n, nkp, a_real, a_imag
if (proj2wann_map(n) < 0) cycle
a_matrix(m, proj2wann_map(n), nkp) = cmplx(a_real, a_imag, kind=dp)
end do
else
do ncount = 1, num_amn
read (amn_in, *, err=104, end=104) m, n, nkp, a_real, a_imag
if (proj2wann_map(n) < 0) cycle
u_matrix(m, proj2wann_map(n), nkp) = cmplx(a_real, a_imag, kind=dp)
end do
end if
close (amn_in)
endif
if (disentanglement) then
call comms_bcast(a_matrix(1, 1, 1), num_bands*num_wann*num_kpts)
else
call comms_bcast(u_matrix(1, 1, 1), num_wann*num_wann*num_kpts)
endif
else
do n = 1, num_kpts
do m = 1, num_wann
u_matrix(m, m, n) = cmplx_1
end do
end do
end if
! If post-processing a Car-Parinello calculation (gamma only)
! then rotate M and A to the basis of Kohn-Sham eigenstates
if (cp_pp) call overlap_rotate()
! Check Mmn(k,b) is symmetric in m and n for gamma_only case
!~ if (gamma_only) call overlap_check_m_symmetry()
! If we don't need to disentangle we can now convert from A to U
! And rotate M accordingly
![ysl-b]
! if(.not.disentanglement .and. (.not.cp_pp) .and. (.not. use_bloch_phases )) &
! call overlap_project
!~ if((.not.cp_pp) .and. (.not. use_bloch_phases )) then
!~ if (.not.disentanglement) then
!~ if ( .not. gamma_only ) then
!~ call overlap_project
!~ else
!~ call overlap_project_gamma()
!~ endif
!~ else
!~ if (gamma_only) call overlap_symmetrize()
!~ endif
!~ endif
!
!~[aam]
if ((.not. disentanglement) .and. (.not. cp_pp) .and. (.not. use_bloch_phases)) then
if (.not. gamma_only) then
call overlap_project()
else
call overlap_project_gamma()
endif
endif
!~[aam]
!~ if( gamma_only .and. use_bloch_phases ) then
!~ write(stdout,'(1x,"+",76("-"),"+")')
!~ write(stdout,'(3x,a)') 'WARNING: gamma_only and use_bloch_phases '
!~ write(stdout,'(3x,a)') ' M must be calculated from *real* Bloch functions'
!~ write(stdout,'(1x,"+",76("-"),"+")')
!~ end if
![ysl-e]
if (timing_level > 0) call io_stopwatch('overlap: read', 2)
return
101 call io_error('Error: Problem opening input file '//trim(seedname)//'.mmn')
102 call io_error('Error: Problem opening input file '//trim(seedname)//'.amn')
103 call io_error('Error: Problem reading input file '//trim(seedname)//'.mmn')
104 call io_error('Error: Problem reading input file '//trim(seedname)//'.amn')
end subroutine overlap_read