diff --git a/src/sitesym.F90 b/src/sitesym.F90 index 4b26fd7d..19fa8ea3 100644 --- a/src/sitesym.F90 +++ b/src/sitesym.F90 @@ -560,7 +560,7 @@ subroutine sitesym_dis_extract_symmetry(sitesym, lambda, umat, zmat, ik, n, num_ !================================================! use w90_wannier90_types, only: sitesym_type - use w90_error, only: w90_error_type, set_error_fatal + use w90_error, only: w90_error_type, set_error_fatal, set_error_alloc implicit none @@ -579,15 +579,40 @@ subroutine sitesym_dis_extract_symmetry(sitesym, lambda, umat, zmat, ik, n, num_ complex(kind=dp), intent(inout) :: umat(:, :) !(num_bands, num_wann) ! local variables - complex(kind=dp) :: umatnew(num_bands, num_wann) !jj normally don't we alloc explicitly? - complex(kind=dp) :: ZU(num_bands, num_wann) - complex(kind=dp) :: deltaU(num_bands, num_wann), carr(num_bands) + complex(kind=dp), allocatable :: umatnew(:,:) !(num_bands, num_wann) + complex(kind=dp), allocatable :: ZU(:,:) !(num_bands, num_wann) + complex(kind=dp), allocatable :: deltaU(:,:) !(num_bands, num_wann) + complex(kind=dp), allocatable :: carr(:) !(num_bands) integer :: i, m, INFO, IFAIL(2), IWORK(5*2) complex(kind=dp) :: HP(3), SP(3), V(2, 2), CWORK(2*2) real(kind=dp) :: W(2), RWORK(7*2), sp3 - integer :: iter + integer :: iter, ierr integer, parameter :: niter = 50 + allocate (umatnew(num_bands, num_wann), stat=ierr) + if (ierr /= 0) then + call set_error_alloc(error, 'Error in allocating umatnew in sitesym_dis_extract_symmetry', comm) + return + endif + + allocate (ZU(num_bands, num_wann), stat=ierr) + if (ierr /= 0) then + call set_error_alloc(error, 'Error in allocating ZU in sitesym_dis_extract_symmetry', comm) + return + endif + + allocate (deltaU(num_bands, num_wann), stat=ierr) + if (ierr /= 0) then + call set_error_alloc(error, 'Error in allocating deltaU in sitesym_dis_extract_symmetry', comm) + return + endif + + allocate (carr(num_bands), stat=ierr) + if (ierr /= 0) then + call set_error_alloc(error, 'Error in allocating carr in sitesym_dis_extract_symmetry', comm) + return + endif + do iter = 1, niter ! Z*U call zgemm('N', 'N', n, num_wann, n, cmplx_1, zmat, num_bands, umat, num_bands, cmplx_0, ZU, &