Commit 676c6978 authored by Daniel Wortmann's avatar Daniel Wortmann

More fixes

Introduced module types_rcmat with a general type for real/complex matrices
Introduced general IO for these matrices
parent 09b6c67a
...@@ -59,9 +59,11 @@ CONTAINS ...@@ -59,9 +59,11 @@ CONTAINS
!---> in a spin-spiral calculation different basis sets are used for !---> in a spin-spiral calculation different basis sets are used for
!---> the two spin directions, because the cutoff radius is defined !---> the two spin directions, because the cutoff radius is defined
!---> by |G + k +/- qss/2| < rkmax. !---> by |G + k +/- qss/2| < rkmax.
IF (nk>kpts%nkpt) THEN
bkpt(:) = kpts%bk(:,nk) bkpt(:)=kpts%bkf(:,nk)
ELSE
bkpt(:) = kpts%bk(:,nk)
ENDIF
!---> Determine rkmax box of size mk1, mk2, mk3, !---> Determine rkmax box of size mk1, mk2, mk3,
! for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax ! for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax
! arltv(i) length of reciprical lattice vector along direction (i) ! arltv(i) length of reciprical lattice vector along direction (i)
......
...@@ -40,6 +40,7 @@ CONTAINS ...@@ -40,6 +40,7 @@ CONTAINS
USE m_icorrkeys USE m_icorrkeys
USE m_eig66_io, ONLY : open_eig, write_eig, close_eig,read_eig USE m_eig66_io, ONLY : open_eig, write_eig, close_eig,read_eig
USE m_xmlOutput USE m_xmlOutput
USE m_io_hybrid
#ifdef CPP_MPI #ifdef CPP_MPI
USE m_mpi_bc_pot USE m_mpi_bc_pot
#endif #endif
...@@ -98,6 +99,7 @@ CONTAINS ...@@ -98,6 +99,7 @@ CONTAINS
TYPE(t_enpara) :: enpara TYPE(t_enpara) :: enpara
TYPE(t_zMat) :: zMat TYPE(t_zMat) :: zMat
TYPE(t_hamOvlp) :: hamOvlp TYPE(t_hamOvlp) :: hamOvlp
TYPE(T_mat) :: olap
! !
INTEGER fh,nn,n INTEGER fh,nn,n
INTEGER ierr(3) INTEGER ierr(3)
...@@ -179,7 +181,7 @@ CONTAINS ...@@ -179,7 +181,7 @@ CONTAINS
input, v%mt,v%vacz, enpara_in, enpara) input, v%mt,v%vacz, enpara_in, enpara)
! !
CALL eigen_hf_init(hybrid,kpts,atoms,input,dimension,hybdat,irank2,isize2) CALL eigen_hf_init(hybrid,kpts,atoms,input,dimension,hybdat,irank2,isize2,l_real)
!---> set up and solve the eigenvalue problem !---> set up and solve the eigenvalue problem
!---> loop over energy windows !---> loop over energy windows
...@@ -267,7 +269,7 @@ CONTAINS ...@@ -267,7 +269,7 @@ CONTAINS
DO jsp = 1,nspins DO jsp = 1,nspins
CALL eigen_HF_setup(hybrid,input,sym,kpts,dimension,atoms,mpi,noco,cell,oneD,results,jsp,eig_id,& CALL eigen_HF_setup(hybrid,input,sym,kpts,dimension,atoms,mpi,noco,cell,oneD,results,jsp,eig_id,&
hybdat,irank2,it,l_real,vr0) hybdat,irank2,it,l_real,vr0,eig_irr)
! !
!---> set up k-point independent t(l'm',lm) matrices !---> set up k-point independent t(l'm',lm) matrices
...@@ -353,6 +355,11 @@ CONTAINS ...@@ -353,6 +355,11 @@ CONTAINS
ENDIF ENDIF
! !
IF( hybrid%l_hybrid ) THEN IF( hybrid%l_hybrid ) THEN
!write overlap matrix b to direct access file olap
print *,"Wrong overlap matrix used, fix this later"
call olap%from_packed(dimension%nbasfcn,l_real,hamovlp%b_r,hamovlp%b_c)
call write_olap(olap,nrec)
CALL hsfock(nk,atoms,hybrid,lapw,DIMENSION,kpts,jsp,input,hybdat,eig_irr,& CALL hsfock(nk,atoms,hybrid,lapw,DIMENSION,kpts,jsp,input,hybdat,eig_irr,&
sym,cell,noco,results,it,maxval(hybdat%nobd),xcpot,& sym,cell,noco,results,it,maxval(hybdat%nobd),xcpot,&
...@@ -383,21 +390,9 @@ CONTAINS ...@@ -383,21 +390,9 @@ CONTAINS
IF (noco%l_noco) CLOSE (25) IF (noco%l_noco) CLOSE (25)
!write overlap matrix b to direct access file olap
INQUIRE(file='olap',exist=l_file)
IF (l_file) THEN
IF (l_real) THEN
OPEN(88,file='olap',form='unformatted',access='direct', recl=matsize*8)
WRITE(88,rec=nrec) hamOvlp%b_r
CLOSE(88)
ELSE
OPEN(88,file='olap',form='unformatted',access='direct', recl=matsize*16)
WRITE(88,rec=nrec) hamOvlp%b_c
CLOSE(88)
ENDIF
ENDIF
CALL eigen_diag(jsp,eig_id,it,atoms,DIMENSION,matsize,mpi, mpi%n_rank,mpi%n_size,ne,nk,lapw,input,& CALL eigen_diag(jsp,eig_id,it,atoms,DIMENSION,matsize,mpi, mpi%n_rank,mpi%n_size,ne,nk,lapw,input,&
nred,mpi%sub_comm, sym,l_zref,matind,kveclo, noco,cell,bkpt,enpara%el0,jij,l_wu,& nred,mpi%sub_comm, sym,l_zref,matind,kveclo, noco,cell,bkpt,enpara%el0,jij,l_wu,&
oneD,td,ud, eig,ne_found,hamOvlp,zMat) oneD,td,ud, eig,ne_found,hamOvlp,zMat)
......
...@@ -27,5 +27,6 @@ global/phasy1.f90 ...@@ -27,5 +27,6 @@ global/phasy1.f90
global/qfix.f90 global/qfix.f90
global/radflo.F90 global/radflo.F90
global/types.F90 global/types.F90
global/types_rcmat.F90
global/utility.F90 global/utility.F90
) )
...@@ -5,6 +5,7 @@ ...@@ -5,6 +5,7 @@
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
MODULE m_types MODULE m_types
USE m_types_rcmat
!************************************************************* !*************************************************************
! This module contains definitions for all kind of types ! This module contains definitions for all kind of types
!************************************************************* !*************************************************************
...@@ -83,7 +84,7 @@ MODULE m_types ...@@ -83,7 +84,7 @@ MODULE m_types
REAL,ALLOCATABLE,DIMENSION(:,:,:) :: dulon !(nlod,ntype,jspd) REAL,ALLOCATABLE,DIMENSION(:,:,:) :: dulon !(nlod,ntype,jspd)
REAL,ALLOCATABLE,DIMENSION(:,:,:,:) :: uloulopn! (nlod,nlod,ntypd,jspd) REAL,ALLOCATABLE,DIMENSION(:,:,:,:) :: uloulopn! (nlod,nlod,ntypd,jspd)
CONTAINS CONTAINS
PROCEDURE init => usdus_init PROCEDURE :: init => usdus_init
END TYPE t_usdus END TYPE t_usdus
...@@ -765,14 +766,7 @@ MODULE m_types ...@@ -765,14 +766,7 @@ MODULE m_types
COMPLEX, ALLOCATABLE :: a_c(:), b_c(:) COMPLEX, ALLOCATABLE :: a_c(:), b_c(:)
END TYPE t_hamOvlp END TYPE t_hamOvlp
TYPE t_lapwmat
LOGICAL :: l_real
INTEGER :: matsize
REAL, ALLOCATABLE :: mat_r(:)
COMPLEX, ALLOCATABLE :: mat_c(:)
CONTAINS
PROCEDURE allocate_space =>t_lapwmat_allocate
END TYPE t_lapwmat
! !
! type for wannier-functions ! type for wannier-functions
! !
...@@ -947,19 +941,5 @@ CONTAINS ...@@ -947,19 +941,5 @@ CONTAINS
ENDIF ENDIF
END SUBROUTINE init_potden_simple END SUBROUTINE init_potden_simple
SUBROUTINE t_lapwmat_allocate(mat)
USE m_judft
IMPLICIT NONE
CLASS(t_lapwmat),INTENT(INOUT) :: mat
INTEGER:: err
IF (mat%l_real) THEN
if (allocated(mat%mat_r)) CALL juDFT_error("Matrix already allocated",hint="this is a bug in the code,please report")
allocate(mat%mat_r(mat%matsize),stat=err)
ELSE
if (allocated(mat%mat_c)) CALL juDFT_error("Matrix already allocated",hint="this is a bug in the code,please report")
allocate(mat%mat_c(mat%matsize),stat=err)
endif
IF (err>0) CALL judft_error("Not enough memory allocating a lapw-matrix")
end SUBROUTINE t_lapwmat_allocate
END MODULE m_types END MODULE m_types
module m_types_rcmat
use m_judft
IMPLICIT NONE
TYPE t_mat
LOGICAL :: l_real
INTEGER :: matsize1=-1
INTEGER :: matsize2=-1
REAL,ALLOCATABLE :: data_r(:,:)
COMPLEX,ALLOCATABLE :: data_c(:,:)
CONTAINS
PROCEDURE :: alloc => t_mat_alloc
PROCEDURE :: multiply=>t_mat_multiply
PROCEDURE :: transpose=>t_mat_transpose
PROCEDURE :: from_packed=>t_mat_from_packed
END type t_mat
CONTAINS
SUBROUTINE t_mat_alloc(mat,l_real,matsize1,matsize2)
CLASS(t_mat) :: mat
LOGICAL,INTENT(IN),OPTIONAL:: l_real
INTEGER,INTENT(IN),OPTIONAL:: matsize1,matsize2
INTEGER:: err
IF (present(l_real)) mat%l_real=l_real
IF (present(matsize1)) mat%matsize1=matsize1
IF (present(matsize2)) mat%matsize2=matsize2
IF (mat%matsize1<0) CALL judft_error("Cannot allocate memory for mat datatype that is not initialized",hint="This is a BUG in FLEUR, please report")
IF (mat%matsize2<0) mat%matsize2=mat%matsize1 !Square by default
IF (allocated(mat%data_r)) deallocate(mat%data_r)
IF (allocated(mat%data_c)) deallocate(mat%data_c)
IF (mat%l_real) THEN
ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2),STAT=err)
ELSE
ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2),STAT=err)
ENDIF
IF (err>0) CALL judft_error("Allocation of memmory failed for mat datatype",hint="You probably run out of memory")
END SUBROUTINE t_mat_alloc
SUBROUTINE t_mat_multiply(mat1,mat2,res)
CLASS(t_mat),INTENT(INOUT) ::mat1
TYPE(t_mat),INTENT(IN) ::mat2
TYPE(t_mat),INTENT(OUT),OPTIONAL ::res
if (mat1%matsize2.ne.mat2%matsize1) CALL judft_error("Cannot multiply matrices because of non-matching dimensions",hint="This is a BUG in FLEUR, please report")
IF (present(res)) THEN
call res%alloc(mat1%l_real,mat1%matsize1,mat2%matsize2)
IF (mat1%l_real) THEN
res%data_r=matmul(mat1%data_r(:mat1%matsize1,:mat1%matsize2),mat2%data_r(:mat2%matsize1,:mat2%matsize2))
ELSE
res%data_c=matmul(mat1%data_c(:mat1%matsize1,:mat1%matsize2),mat2%data_c(:mat2%matsize1,:mat2%matsize2))
ENDIF
else
if (mat1%matsize1.ne.mat1%matsize2) CALL judft_error("Cannot multiply matrices inplace because of non-matching dimensions",hint="This is a BUG in FLEUR, please report")
if (mat1%l_real) THEN
mat1%data_r(:mat1%matsize1,:mat1%matsize2)=matmul(mat1%data_r(:mat1%matsize1,:mat1%matsize2),mat2%data_r(:mat2%matsize1,:mat2%matsize2))
ELSE
mat1%data_c(:mat1%matsize1,:mat1%matsize2)=matmul(mat1%data_c(:mat1%matsize1,:mat1%matsize2),mat2%data_c(:mat2%matsize1,:mat2%matsize2))
ENDIF
end IF
end SUBROUTINE t_mat_multiply
SUBROUTINE t_mat_transpose(mat1,res)
CLASS(t_mat),INTENT(INOUT) ::mat1
TYPE(t_mat),INTENT(OUT),OPTIONAL ::res
IF (present(res)) THEN
call res%alloc(mat1%l_real,mat1%matsize2,mat1%matsize1)
IF (mat1%l_real) THEN
res%data_r=transpose(mat1%data_r(:mat1%matsize1,:mat1%matsize2))
ELSE
res%data_c=transpose(mat1%data_c(:mat1%matsize1,:mat1%matsize2))
ENDIF
else
if (mat1%matsize1.ne.mat1%matsize2) CALL judft_error("Cannot transpose matrices inplace because of non-matching dimensions",hint="This is a BUG in FLEUR, please report")
IF (mat1%l_real) THEN
mat1%data_r(:mat1%matsize1,:mat1%matsize2)=transpose(mat1%data_r(:mat1%matsize1,:mat1%matsize2))
ELSE
mat1%data_c(:mat1%matsize1,:mat1%matsize2)=transpose(mat1%data_c(:mat1%matsize1,:mat1%matsize2))
ENDIF
end IF
end SUBROUTINE t_mat_transpose
SUBROUTINE t_mat_from_packed(mat1,matsize,l_real,packed_r,packed_c)
CLASS(t_mat),INTENT(INOUT) :: mat1
INTEGER,INTENT(IN) :: matsize
LOGICAL,INTENT(IN) :: l_real
REAL,INTENT(IN) :: packed_r(:)
COMPLEX,INTENT(IN) :: packed_c(:)
INTEGER:: n,nn,i
call mat1%alloc(l_real,matsize)
i=1
DO n=1,matsize
DO nn=1,n
if (l_real) THEN
mat1%data_r(n,nn)=packed_r(i)
else
mat1%data_c(n,nn)=packed_c(i)
end if
i=i+1
end DO
end DO
end SUBROUTINE t_mat_from_packed
end module m_types_rcmat
MODULE m_abcrot MODULE m_abcrot
CONTAINS CONTAINS
SUBROUTINE abcrot(atoms,neig,sym,cell,oneD,& SUBROUTINE abcrot(hybrid,atoms,neig,sym,cell,oneD,&
& acof,bcof,ccof) & acof,bcof,ccof)
! *************************************************************** ! ***************************************************************
! * This routine transforms a/b/cof which are given wrt rotated * ! * This routine transforms a/b/cof which are given wrt rotated *
...@@ -13,6 +13,7 @@ MODULE m_abcrot ...@@ -13,6 +13,7 @@ MODULE m_abcrot
USE m_dwigner USE m_dwigner
USE m_types USE m_types
IMPLICIT NONE IMPLICIT NONE
TYPE(t_hybrid),INTENT(IN) :: hybrid
TYPE(t_oneD),INTENT(IN) :: oneD TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_sym),INTENT(IN) :: sym TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell TYPE(t_cell),INTENT(IN) :: cell
...@@ -34,7 +35,7 @@ MODULE m_abcrot ...@@ -34,7 +35,7 @@ MODULE m_abcrot
!***** COMPLEX, ALLOCATABLE :: d_wgn(:,:,:,:) !put into module m_savewigner !***** COMPLEX, ALLOCATABLE :: d_wgn(:,:,:,:) !put into module m_savewigner
! !
IF ( .NOT.ALLOCATED(sym%d_wgn) ) THEN !calculate sym%d_wgn only once IF ( .NOT.ALLOCATED(hybrid%d_wgn2) ) THEN !calculate sym%d_wgn only once
#ifndef CPP_MPI #ifndef CPP_MPI
PRINT*,"calculate wigner-matrix" PRINT*,"calculate wigner-matrix"
#endif #endif
...@@ -72,15 +73,15 @@ MODULE m_abcrot ...@@ -72,15 +73,15 @@ MODULE m_abcrot
DO l=1,atoms%lmax(itype) DO l=1,atoms%lmax(itype)
! replaced d_wgn by conjg(d_wgn),FF October 2006 ! replaced d_wgn by conjg(d_wgn),FF October 2006
DO i=1,neig DO i=1,neig
acof(i,l**2:l*(l+2),iatom) = ifac**l * matmul(conjg(sym%d_wgn(-l:l,-l:l,l,iop)), acof(i,l**2:l*(l+2),iatom)) acof(i,l**2:l*(l+2),iatom) = ifac**l * matmul(conjg(hybrid%d_wgn2(-l:l,-l:l,l,iop)), acof(i,l**2:l*(l+2),iatom))
bcof(i,l**2:l*(l+2),iatom) = ifac**l * matmul(conjg(sym%d_wgn(-l:l,-l:l,l,iop)), bcof(i,l**2:l*(l+2),iatom)) bcof(i,l**2:l*(l+2),iatom) = ifac**l * matmul(conjg(hybrid%d_wgn2(-l:l,-l:l,l,iop)), bcof(i,l**2:l*(l+2),iatom))
ENDDO ENDDO
ENDDO ENDDO
DO ilo=1,atoms%nlo(itype) DO ilo=1,atoms%nlo(itype)
l=atoms%llo(ilo,itype) l=atoms%llo(ilo,itype)
IF(l.gt.0) THEN IF(l.gt.0) THEN
DO i=1,neig DO i=1,neig
ccof(-l:l,i,ilo,iatom) = ifac**l * matmul(conjg(sym%d_wgn(-l:l,-l:l,l,iop)), ccof(-l:l,i,ilo,iatom)) ccof(-l:l,i,ilo,iatom) = ifac**l * matmul(conjg(hybrid%d_wgn2(-l:l,-l:l,l,iop)), ccof(-l:l,i,ilo,iatom))
ENDDO ENDDO
ENDIF ENDIF
ENDDO ENDDO
......
...@@ -12,6 +12,7 @@ ...@@ -12,6 +12,7 @@
USE m_constants USE m_constants
USE m_apws USE m_apws
USE m_types USE m_types
USE m_io_hybrid
IMPLICIT NONE IMPLICIT NONE
TYPE(t_hybdat),INTENT(IN) :: hybdat TYPE(t_hybdat),INTENT(IN) :: hybdat
...@@ -89,21 +90,17 @@ ...@@ -89,21 +90,17 @@
cmt = 0 cmt = 0
! initialize gridf ! initialize gridf -> was done in eigen_HF_init
CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf) !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
! read in cmt ! read in cmt
irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16
OPEN(unit=777,file='cmt',form='unformatted',access='direct',&
& recl=irecl_cmt)
DO ikpt = 1,nkpti DO ikpt = 1,nkpti
#ifdef CPP_MPI #ifdef CPP_MPI
IF ( skip_kpt(ikpt) ) CYCLE IF ( skip_kpt(ikpt) ) CYCLE
#endif #endif
READ(777,rec=ikpt) cmt(:,:,:,ikpt) call read_cmt(cmt(:,:,:,ikpt),ikpt)
END DO END DO
CLOSE(777)
IF ( mpi%irank == 0 ) WRITE(6,'(/A)') ' Overlap <core|core>' IF ( mpi%irank == 0 ) WRITE(6,'(/A)') ' Overlap <core|core>'
DO itype=1,atoms%ntype DO itype=1,atoms%ntype
IF(atoms%ntype.gt.1 .AND. mpi%irank==0) & IF(atoms%ntype.gt.1 .AND. mpi%irank==0) &
......
...@@ -2064,12 +2064,12 @@ ...@@ -2064,12 +2064,12 @@
# else # else
!write coulomb matrix to direct access file coulomb !write coulomb matrix to direct access file coulomb
OPEN(unit=777,file=fname,form='unformatted',access='direct',& OPEN(unit=677,file=fname,form='unformatted',access='direct',&
& recl=irecl_coulomb) & recl=irecl_coulomb)
DO i=1,kpts%nkpt DO i=1,kpts%nkpt
WRITE(777,rec=i) coulomb(:,i) WRITE(677,rec=i) coulomb(:,i)
END DO END DO
CLOSE(777) CLOSE(677)
# endif # endif
......
...@@ -4,11 +4,11 @@ module m_eigen_hf_init ...@@ -4,11 +4,11 @@ module m_eigen_hf_init
! preparations for HF and hybrid functional calculation ! preparations for HF and hybrid functional calculation
! !
contains contains
subroutine eigen_hf_init(hybrid,kpts,atoms,input,dimension,hybdat,irank2,isize2) subroutine eigen_hf_init(hybrid,kpts,atoms,input,dimension,hybdat,irank2,isize2,l_real)
USE m_types USE m_types
USE m_read_core USE m_read_core
USE m_util USE m_util
USE m_io_hybrid
implicit none implicit none
TYPE(t_hybrid),INTENT(INOUT) :: hybrid TYPE(t_hybrid),INTENT(INOUT) :: hybrid
TYPE(t_kpts),INTENT(IN) :: kpts TYPE(t_kpts),INTENT(IN) :: kpts
...@@ -17,6 +17,7 @@ contains ...@@ -17,6 +17,7 @@ contains
TYPE(t_dimension),INTENT(IN) :: dimension TYPE(t_dimension),INTENT(IN) :: dimension
INTEGER,INTENT(OUT) :: irank2(:),isize2(:) INTEGER,INTENT(OUT) :: irank2(:),isize2(:)
TYPE(t_hybdat),INTENT(OUT) :: hybdat TYPE(t_hybdat),INTENT(OUT) :: hybdat
LOGICAL,INTENT(IN) :: l_real
...@@ -42,7 +43,9 @@ contains ...@@ -42,7 +43,9 @@ contains
END IF END IF
IF( hybrid%l_calhf ) THEN IF( hybrid%l_calhf ) THEN
hybrid%maxlmindx = MAXVAL((/ ( SUM( (/ (hybrid%nindx(l,itype)*(2*l+1), l=0,atoms%lmax(itype)) /) ),itype=1,atoms%ntype) /) )
call open_hybrid_io(hybrid,dimension,atoms,l_real)
!initialize hybdat%gridf for radial integration !initialize hybdat%gridf for radial integration
CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf) CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
......
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
MODULE m_eigen_hf_setup MODULE m_eigen_hf_setup
CONTAINS CONTAINS
SUBROUTINE eigen_hf_setup(hybrid,input,sym,kpts,DIMENSION,atoms,mpi,noco,cell,oneD,results,jsp,eig_id_hf,& SUBROUTINE eigen_hf_setup(hybrid,input,sym,kpts,DIMENSION,atoms,mpi,noco,cell,oneD,results,jsp,eig_id_hf,&
hybdat,irank2,it,l_real,vr0) hybdat,irank2,it,l_real,vr0,eig_irr)
USE m_types USE m_types
USE m_eig66_io USE m_eig66_io
USE m_util USE m_util
...@@ -28,13 +28,13 @@ CONTAINS ...@@ -28,13 +28,13 @@ CONTAINS
REAL, INTENT(IN) :: vr0(:,:,:) REAL, INTENT(IN) :: vr0(:,:,:)
LOGICAL,INTENT(IN) :: l_real LOGICAL,INTENT(IN) :: l_real
TYPE(t_hybdat),INTENT(INOUT):: hybdat TYPE(t_hybdat),INTENT(INOUT):: hybdat
REAL,ALLOCATABLE,INTENT(OUT):: eig_irr(:,:)
INTEGER:: ok,nk,nrec1,i,j,ll,l1,l2,ng,itype,n,l,n1,n2,nn INTEGER:: ok,nk,nrec1,i,j,ll,l1,l2,ng,itype,n,l,n1,n2,nn
TYPE(t_zmat),ALLOCATABLE :: zmat(:) TYPE(t_zmat),ALLOCATABLE :: zmat(:)
REAL, ALLOCATABLE :: eig_irr(:,:)
REAL, ALLOCATABLE :: basprod(:) REAL, ALLOCATABLE :: basprod(:)
REAL :: el_eig(0:atoms%lmaxd,atoms%ntype), ello_eig(atoms%nlod,atoms%ntype),bk(3) REAL :: el_eig(0:atoms%lmaxd,atoms%ntype), ello_eig(atoms%nlod,atoms%ntype),bk(3)
INTEGER :: degenerat(DIMENSION%neigd2+1,kpts%nkpt) INTEGER :: degenerat(DIMENSION%neigd2+1,kpts%nkpt)
...@@ -54,19 +54,20 @@ CONTAINS ...@@ -54,19 +54,20 @@ CONTAINS
ALLOCATE(zmat(kpts%nkptf),stat=ok) ALLOCATE(zmat(kpts%nkptf),stat=ok)
IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation z_c' IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation z_c'
ALLOCATE ( eig_irr(DIMENSION%neigd2,kpts%nkpt) ,stat=ok ) ALLOCATE ( eig_irr(DIMENSION%neigd2,kpts%nkpt) ,stat=ok )
IF( ok .NE. 0 ) STOP'eigen_hf: failure allocation eig_irr' IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation eig_irr'
ALLOCATE ( hybdat%kveclo_eig(atoms%nlotot,kpts%nkpt) ,stat=ok ) ALLOCATE ( hybdat%kveclo_eig(atoms%nlotot,kpts%nkpt) ,stat=ok )
IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation hybdat%kveclo_eig' IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation hybdat%kveclo_eig'
eig_irr = 0 ; hybdat%kveclo_eig = 0 eig_irr = 0 ; hybdat%kveclo_eig = 0
! Reading the eig file
zmat(:)%l_real=l_real
! Reading the eig file
DO nk = 1,kpts%nkpt DO nk = 1,kpts%nkpt
# ifdef CPP_MPI # ifdef CPP_MPI
! jump to next k-point if this process is not present in communicator ! jump to next k-point if this process is not present in communicator
IF ( skip_kpt(nk) ) CYCLE IF ( skip_kpt(nk) ) CYCLE
# endif # endif
nrec1 = kpts%nkpt*(jsp-1) + nk nrec1 = kpts%nkpt*(jsp-1) + nk
zmat(nk)%l_real=l_real
zmat(nk)%nbasfcn=dimension%nbasfcn zmat(nk)%nbasfcn=dimension%nbasfcn
zmat(nk)%nbands=dimension%neigd2 zmat(nk)%nbands=dimension%neigd2
if (l_real) THEN if (l_real) THEN
...@@ -80,7 +81,16 @@ CONTAINS ...@@ -80,7 +81,16 @@ CONTAINS
print *,"Done" print *,"Done"
END DO END DO
!Allocate further space
DO nk=kpts%nkpt+1,kpts%nkptf
zmat(nk)%nbasfcn=dimension%nbasfcn
zmat(nk)%nbands=dimension%neigd2
if (l_real) THEN
ALLOCATE(zmat(nk)%z_r(dimension%nbasfcn,dimension%neigd2))
else
ALLOCATE(zmat(nk)%z_c(dimension%nbasfcn,dimension%neigd2))
endif
Enddo
! !
!determine degenerate states at each k-point !determine degenerate states at each k-point
! !
...@@ -209,6 +219,7 @@ CONTAINS ...@@ -209,6 +219,7 @@ CONTAINS
! !
! setup dimension of pntgpt ! setup dimension of pntgpt
ALLOCATE(hybdat%pntgptd(3))
hybdat%pntgptd = 0 hybdat%pntgptd = 0
DO nk = 1,kpts%nkptf DO nk = 1,kpts%nkptf
CALL apws(DIMENSION,input,noco, kpts,nk,cell,sym%zrfs,& CALL apws(DIMENSION,input,noco, kpts,nk,cell,sym%zrfs,&
...@@ -239,7 +250,7 @@ CONTAINS ...@@ -239,7 +250,7 @@ CONTAINS
IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation hybdat%prod' IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation hybdat%prod'
basprod = 0 ; hybdat%prodm = 0 ; hybdat%prod%l1 = 0 ; hybdat%prod%l2 = 0 basprod = 0 ; hybdat%prodm = 0 ; hybdat%prod%l1 = 0 ; hybdat%prod%l2 = 0
hybdat%prod%n1 = 0 ; hybdat%prod%n2 = 0 hybdat%prod%n1 = 0 ; hybdat%prod%n2 = 0
ALLOCATE(hybrid%nindxp1(0:hybrid%maxlcutm1,atoms%ntype))
hybrid%nindxp1 = 0 hybrid%nindxp1 = 0
DO itype = 1,atoms%ntype DO itype = 1,atoms%ntype
ng = atoms%jri(itype) ng = atoms%jri(itype)
...@@ -281,16 +292,12 @@ CONTAINS ...@@ -281,16 +292,12 @@ CONTAINS
ELSE IF (hybrid%l_hybrid ) THEN ! hybrid%l_calhf is false ELSE IF (hybrid%l_hybrid ) THEN ! hybrid%l_calhf is false
! Reading the eig file ! Reading the eig file
ALLOCATE ( eig_irr(DIMENSION%neigd2,1) )
ALLOCATE ( hybdat%kveclo_eig(atoms%nlotot,1) )
!DO nk = n_start,kpts%nkpt,n_stride !DO nk = n_start,kpts%nkpt,n_stride
DO nk = 1,kpts%nkpt,1 DO nk = 1,kpts%nkpt,1
nrec1 = kpts%nkpt*(jsp-1) + nk
CALL read_eig(eig_id_hf,nk,jsp,el=el_eig, ello=ello_eig,neig=hybdat%ne_eig(nk)) CALL read_eig(eig_id_hf,nk,jsp,el=el_eig, ello=ello_eig,neig=hybdat%ne_eig(nk))
hybdat%nobd(nk) = COUNT( results%w_iks(:hybdat%ne_eig(nk),nk,jsp) > 0.0 ) hybdat%nobd(nk) = COUNT( results%w_iks(:hybdat%ne_eig(nk),nk,jsp) > 0.0 )
END DO END DO
DEALLOCATE ( eig_irr , hybdat%kveclo_eig )
hybrid%maxlmindx = MAXVAL((/ ( SUM( (/ (hybrid%nindx(l,itype)*(2*l+1), l=0,atoms%lmax(itype)) /) ),itype=1,atoms%ntype) /) ) hybrid%maxlmindx = MAXVAL((/ ( SUM( (/ (hybrid%nindx(l,itype)*(2*l+1), l=0,atoms%lmax(itype)) /) ),itype=1,atoms%ntype) /) )
hybdat%nbands = MIN