From 676c697860f4ecbe2610a0eda1d8e6f0796488ee Mon Sep 17 00:00:00 2001 From: Daniel Wortmann Date: Tue, 11 Jul 2017 10:25:26 +0200 Subject: [PATCH] More fixes Introduced module types_rcmat with a general type for real/complex matrices Introduced general IO for these matrices --- eigen/apws.F90 | 8 +- eigen/eigen.F90 | 29 ++--- global/CMakeLists.txt | 1 + global/types.F90 | 28 +---- global/types_rcmat.F90 | 110 ++++++++++++++++ hybrid/abcrot.F90 | 11 +- hybrid/checkolap.F90 | 15 +-- hybrid/coulombmatrix.F90 | 6 +- hybrid/eigen_HF_init.F90 | 9 +- hybrid/eigen_HF_setup.F90 | 33 +++-- hybrid/exchange_core.F90 | 52 +++----- hybrid/exchange_val_hf.F90 | 19 +-- hybrid/gen_wavf.F90 | 72 ++++------- hybrid/hsfock.F90 | 251 ++++++++----------------------------- hybrid/kp_perturbation.F90 | 13 +- hybrid/symm_hf.F90 | 20 +-- hybrid/symmetrizeh.F90 | 121 +++++++++--------- hybrid/wavefproducts.F90 | 39 +++--- init/gen_bz.F90 | 13 +- init/init_hybrid.f90 | 37 ++++++ init/kpoints.f90 | 61 +++++++++ io/CMakeLists.txt | 2 + io/eig66_da.F90 | 12 +- io/io_hybrid.F90 | 98 +++++++++++++++ io/io_matrix.F90 | 101 +++++++++++++++ math/util.F | 2 +- xc-pot/trafo.F90 | 16 ++- 27 files changed, 687 insertions(+), 492 deletions(-) create mode 100644 global/types_rcmat.F90 create mode 100644 init/init_hybrid.f90 create mode 100644 init/kpoints.f90 create mode 100644 io/io_hybrid.F90 create mode 100644 io/io_matrix.F90 diff --git a/eigen/apws.F90 b/eigen/apws.F90 index a76185c8..60917df1 100644 --- a/eigen/apws.F90 +++ b/eigen/apws.F90 @@ -59,9 +59,11 @@ CONTAINS !---> in a spin-spiral calculation different basis sets are used for !---> the two spin directions, because the cutoff radius is defined !---> by |G + k +/- qss/2| < rkmax. - - bkpt(:) = kpts%bk(:,nk) - + IF (nk>kpts%nkpt) THEN + bkpt(:)=kpts%bkf(:,nk) + ELSE + bkpt(:) = kpts%bk(:,nk) + ENDIF !---> Determine rkmax box of size mk1, mk2, mk3, ! for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax ! arltv(i) length of reciprical lattice vector along direction (i) diff --git a/eigen/eigen.F90 b/eigen/eigen.F90 index 42d5853c..2e503890 100644 --- a/eigen/eigen.F90 +++ b/eigen/eigen.F90 @@ -40,6 +40,7 @@ CONTAINS USE m_icorrkeys USE m_eig66_io, ONLY : open_eig, write_eig, close_eig,read_eig USE m_xmlOutput + USE m_io_hybrid #ifdef CPP_MPI USE m_mpi_bc_pot #endif @@ -98,6 +99,7 @@ CONTAINS TYPE(t_enpara) :: enpara TYPE(t_zMat) :: zMat TYPE(t_hamOvlp) :: hamOvlp + TYPE(T_mat) :: olap ! INTEGER fh,nn,n INTEGER ierr(3) @@ -179,7 +181,7 @@ CONTAINS 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 !---> loop over energy windows @@ -267,7 +269,7 @@ CONTAINS DO jsp = 1,nspins 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 @@ -353,6 +355,11 @@ CONTAINS ENDIF ! 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,& sym,cell,noco,results,it,maxval(hybdat%nobd),xcpot,& @@ -383,21 +390,9 @@ CONTAINS 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,& 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) diff --git a/global/CMakeLists.txt b/global/CMakeLists.txt index 1f2fcea8..a63ebe7a 100644 --- a/global/CMakeLists.txt +++ b/global/CMakeLists.txt @@ -27,5 +27,6 @@ global/phasy1.f90 global/qfix.f90 global/radflo.F90 global/types.F90 +global/types_rcmat.F90 global/utility.F90 ) diff --git a/global/types.F90 b/global/types.F90 index dfd32188..54423f28 100644 --- a/global/types.F90 +++ b/global/types.F90 @@ -5,6 +5,7 @@ !-------------------------------------------------------------------------------- MODULE m_types + USE m_types_rcmat !************************************************************* ! This module contains definitions for all kind of types !************************************************************* @@ -83,7 +84,7 @@ MODULE m_types REAL,ALLOCATABLE,DIMENSION(:,:,:) :: dulon !(nlod,ntype,jspd) REAL,ALLOCATABLE,DIMENSION(:,:,:,:) :: uloulopn! (nlod,nlod,ntypd,jspd) CONTAINS - PROCEDURE init => usdus_init + PROCEDURE :: init => usdus_init END TYPE t_usdus @@ -765,14 +766,7 @@ MODULE m_types COMPLEX, ALLOCATABLE :: a_c(:), b_c(:) 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 ! @@ -947,19 +941,5 @@ CONTAINS ENDIF 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 diff --git a/global/types_rcmat.F90 b/global/types_rcmat.F90 new file mode 100644 index 00000000..18379f99 --- /dev/null +++ b/global/types_rcmat.F90 @@ -0,0 +1,110 @@ +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 diff --git a/hybrid/abcrot.F90 b/hybrid/abcrot.F90 index 4a0cbe89..6272fc58 100644 --- a/hybrid/abcrot.F90 +++ b/hybrid/abcrot.F90 @@ -1,7 +1,7 @@ MODULE m_abcrot CONTAINS - SUBROUTINE abcrot(atoms,neig,sym,cell,oneD,& + SUBROUTINE abcrot(hybrid,atoms,neig,sym,cell,oneD,& & acof,bcof,ccof) ! *************************************************************** ! * This routine transforms a/b/cof which are given wrt rotated * @@ -13,6 +13,7 @@ MODULE m_abcrot USE m_dwigner USE m_types IMPLICIT NONE + TYPE(t_hybrid),INTENT(IN) :: hybrid TYPE(t_oneD),INTENT(IN) :: oneD TYPE(t_sym),INTENT(IN) :: sym TYPE(t_cell),INTENT(IN) :: cell @@ -34,7 +35,7 @@ MODULE m_abcrot !***** 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 PRINT*,"calculate wigner-matrix" #endif @@ -72,15 +73,15 @@ MODULE m_abcrot DO l=1,atoms%lmax(itype) ! replaced d_wgn by conjg(d_wgn),FF October 2006 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)) - 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)) + 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(hybrid%d_wgn2(-l:l,-l:l,l,iop)), bcof(i,l**2:l*(l+2),iatom)) ENDDO ENDDO DO ilo=1,atoms%nlo(itype) l=atoms%llo(ilo,itype) IF(l.gt.0) THEN 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 ENDIF ENDDO diff --git a/hybrid/checkolap.F90 b/hybrid/checkolap.F90 index dd75f15c..148deeb4 100644 --- a/hybrid/checkolap.F90 +++ b/hybrid/checkolap.F90 @@ -12,6 +12,7 @@ USE m_constants USE m_apws USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat @@ -89,21 +90,17 @@ cmt = 0 - ! initialize gridf - CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf) + ! initialize gridf -> was done in eigen_HF_init + !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf) ! 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 #ifdef CPP_MPI IF ( skip_kpt(ikpt) ) CYCLE #endif - READ(777,rec=ikpt) cmt(:,:,:,ikpt) - END DO - CLOSE(777) - + call read_cmt(cmt(:,:,:,ikpt),ikpt) + END DO + IF ( mpi%irank == 0 ) WRITE(6,'(/A)') ' Overlap ' DO itype=1,atoms%ntype IF(atoms%ntype.gt.1 .AND. mpi%irank==0) & diff --git a/hybrid/coulombmatrix.F90 b/hybrid/coulombmatrix.F90 index 32d8faba..b2832932 100644 --- a/hybrid/coulombmatrix.F90 +++ b/hybrid/coulombmatrix.F90 @@ -2064,12 +2064,12 @@ # else !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) DO i=1,kpts%nkpt - WRITE(777,rec=i) coulomb(:,i) + WRITE(677,rec=i) coulomb(:,i) END DO - CLOSE(777) + CLOSE(677) # endif diff --git a/hybrid/eigen_HF_init.F90 b/hybrid/eigen_HF_init.F90 index 1e5502ee..f78b111d 100644 --- a/hybrid/eigen_HF_init.F90 +++ b/hybrid/eigen_HF_init.F90 @@ -4,11 +4,11 @@ module m_eigen_hf_init ! preparations for HF and hybrid functional calculation ! 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_read_core USE m_util - + USE m_io_hybrid implicit none TYPE(t_hybrid),INTENT(INOUT) :: hybrid TYPE(t_kpts),INTENT(IN) :: kpts @@ -17,6 +17,7 @@ contains TYPE(t_dimension),INTENT(IN) :: dimension INTEGER,INTENT(OUT) :: irank2(:),isize2(:) TYPE(t_hybdat),INTENT(OUT) :: hybdat + LOGICAL,INTENT(IN) :: l_real @@ -42,7 +43,9 @@ contains END IF 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 CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf) diff --git a/hybrid/eigen_HF_setup.F90 b/hybrid/eigen_HF_setup.F90 index c89dd665..f718f792 100644 --- a/hybrid/eigen_HF_setup.F90 +++ b/hybrid/eigen_HF_setup.F90 @@ -2,7 +2,7 @@ MODULE m_eigen_hf_setup CONTAINS 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_eig66_io USE m_util @@ -28,13 +28,13 @@ CONTAINS REAL, INTENT(IN) :: vr0(:,:,:) LOGICAL,INTENT(IN) :: l_real 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 TYPE(t_zmat),ALLOCATABLE :: zmat(:) - REAL, ALLOCATABLE :: eig_irr(:,:) REAL, ALLOCATABLE :: basprod(:) REAL :: el_eig(0:atoms%lmaxd,atoms%ntype), ello_eig(atoms%nlod,atoms%ntype),bk(3) INTEGER :: degenerat(DIMENSION%neigd2+1,kpts%nkpt) @@ -54,19 +54,20 @@ CONTAINS ALLOCATE(zmat(kpts%nkptf),stat=ok) IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation z_c' 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 ) IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation hybdat%kveclo_eig' 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 # ifdef CPP_MPI ! jump to next k-point if this process is not present in communicator IF ( skip_kpt(nk) ) CYCLE # endif nrec1 = kpts%nkpt*(jsp-1) + nk - zmat(nk)%l_real=l_real zmat(nk)%nbasfcn=dimension%nbasfcn zmat(nk)%nbands=dimension%neigd2 if (l_real) THEN @@ -80,7 +81,16 @@ CONTAINS print *,"Done" 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 ! @@ -209,6 +219,7 @@ CONTAINS ! ! setup dimension of pntgpt + ALLOCATE(hybdat%pntgptd(3)) hybdat%pntgptd = 0 DO nk = 1,kpts%nkptf CALL apws(DIMENSION,input,noco, kpts,nk,cell,sym%zrfs,& @@ -239,7 +250,7 @@ CONTAINS IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation hybdat%prod' basprod = 0 ; hybdat%prodm = 0 ; hybdat%prod%l1 = 0 ; hybdat%prod%l2 = 0 hybdat%prod%n1 = 0 ; hybdat%prod%n2 = 0 - + ALLOCATE(hybrid%nindxp1(0:hybrid%maxlcutm1,atoms%ntype)) hybrid%nindxp1 = 0 DO itype = 1,atoms%ntype ng = atoms%jri(itype) @@ -281,16 +292,12 @@ CONTAINS ELSE IF (hybrid%l_hybrid ) THEN ! hybrid%l_calhf is false ! 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 = 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)) hybdat%nobd(nk) = COUNT( results%w_iks(:hybdat%ne_eig(nk),nk,jsp) > 0.0 ) 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) /) ) hybdat%nbands = MIN( hybrid%bands1, DIMENSION%neigd ) diff --git a/hybrid/exchange_core.F90 b/hybrid/exchange_core.F90 index cf444de6..b5ca172d 100644 --- a/hybrid/exchange_core.F90 +++ b/hybrid/exchange_core.F90 @@ -33,6 +33,7 @@ USE m_util USE m_wrapper USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat @@ -62,7 +63,7 @@ INTEGER :: iatom,ieq,itype,ic,l,l1,l2,& & ll,lm ,m1,m2,p1,p2,n,n1,n2,i,j INTEGER :: iband1,iband2,ndb1,ndb2,ic1,ic2 - INTEGER :: irecl_cmt,m + INTEGER :: m REAL :: time1,time2 REAL :: rdum @@ -91,11 +92,7 @@ END IF ! read in mt wavefunction coefficients from file cmt - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct', recl=irecl_cmt) - READ(777,rec=nk) cmt(:,:,:) - CLOSE(777) - + call read_cmt(cmt,nk) ALLOCATE ( fprod(atoms%jmtd,5),larr(5),parr(5) ) ! generate ldum(nbands(nk),nbands(nk)), which is true if the corresponding matrix entry is non-zero @@ -261,6 +258,7 @@ USE m_util USE m_wrapper USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat @@ -279,17 +277,12 @@ ! - arays - INTEGER,INTENT(IN) :: nsest(hybdat%nbands(nk)),indx_sest(hybdat%nbands(nk),hybdat%nbands(nk)) - -#ifdef CPP_INVERSION - REAL ,INTENT(INOUT) :: mat_ex(dimension%nbasfcn*(dimension%nbasfcn+1)/2) -#else - COMPLEX ,INTENT(INOUT) :: mat_ex(dimension%nbasfcn*(dimension%nbasfcn+1)/2) -#endif + TYPE(t_mat),INTENT(INOUT):: mat_ex ! - local scalars - INTEGER :: iatom,ieq,itype,ic,l,l1,l2,& & ll,lm ,m1,m2,p1,p2,n,n1,n2,nn2,i,j INTEGER :: iband1,iband2,ndb1,ndb2,ic1,ic2 - INTEGER :: irecl_cmt,m + INTEGER :: m REAL :: time1,time2 REAL :: rdum @@ -314,11 +307,9 @@ ! read in mt wavefunction coefficients from file cmt - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct', recl=irecl_cmt) - READ(777,rec=nk) cmt(:,:,:) - CLOSE(777) - + + CALL read_cmt(cmt,nk) + ALLOCATE ( fprod(atoms%jmtd,5),larr(5),parr(5) ) exchange = 0 @@ -430,13 +421,11 @@ ic = 0 sum_offdia = 0 - DO n1=1,hybdat%nbands(nk) - DO n2=1,n1 - ic = ic + 1 - mat_ex(ic) = mat_ex(ic) + conjg(exchange(n2,n1))/nsymop - END DO - END DO - + if (mat_ex%l_real) THEN + mat_ex%data_r=mat_ex%data_r+exchange/nsymop + else + mat_ex%data_c=mat_ex%data_c+conjg(exchange)/nsymop + end if END SUBROUTINE exchange_vccv1 @@ -451,6 +440,7 @@ USE m_gaunt USE m_trafo USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat TYPE(t_results),INTENT(INOUT) :: results @@ -599,8 +589,8 @@ USE m_wrapper USE m_gaunt USE m_trafo - - USE m_types + USE m_io_hybrid + USE m_types IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat TYPE(t_mpi),INTENT(IN) :: mpi @@ -626,7 +616,6 @@ INTEGER :: lm2,lmp2 INTEGER :: m1,m2 ,mm,m INTEGER :: n1,n2,n,nn - INTEGER :: irecl_cmt REAL :: rdum0,rdum1,rdum2,rdum3,rdum4 COMPLEX :: cdum @@ -670,11 +659,8 @@ END DO ! read in cmt coefficient at k-point nk - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct', recl=irecl_cmt) - READ(777,rec=nk) cmt(:,:,:) - CLOSE(777) - + + call read_cmt(cmt,nk) iatom = 0 DO itype = 1,atoms%ntype DO ieq = 1,atoms%neq(itype) diff --git a/hybrid/exchange_val_hf.F90 b/hybrid/exchange_val_hf.F90 index 22281a08..3f3e4dc4 100644 --- a/hybrid/exchange_val_hf.F90 +++ b/hybrid/exchange_val_hf.F90 @@ -111,11 +111,7 @@ REAL ,INTENT(IN) :: wl_iks(dimension%neigd,kpts%nkptf) REAL ,INTENT(OUT) :: div_vv(hybdat%nbands(nk)) -#ifdef CPP_INVERSION - REAL ,INTENT(OUT) :: mat_ex(dimension%nbasfcn*(dimension%nbasfcn+1)/2) -#else - COMPLEX,INTENT(OUT) :: mat_ex(dimension%nbasfcn*(dimension%nbasfcn+1)/2) -#endif + TYPE(t_mat),INTENT(INOUT):: mat_ex ! - local scalars - INTEGER :: iband,iband1,ibando,ikpt,ikpt0 @@ -217,7 +213,7 @@ vol = cell%omtil svol = sqrt(cell%omtil) - rws = (3*cell%omtil/fpi_const)**(1d0/3) ! Wigner-Seitz radius + rws = (3*cell%omtil/fpi_const)**(1.0/3) ! Wigner-Seitz radius #if ALGORITHM == 1 npot = 3 ! for switching-off function k0 = hybrid%radshmin / 2 ! radius of largest sphere that fits inside the BZ @@ -699,14 +695,9 @@ #endif ! write exch_vv in mat_ex - ic = 0 - DO n1=1,hybdat%nbands(nk) - DO n2=1,n1 - ic = ic + 1 - mat_ex(ic) = mat_ex(ic) + exch_vv(n2,n1) - END DO - END DO - + call mat_ex%alloc(matsize1=hybdat%nbands(nk)) + mat_ex%data_c=exch_vv(n2,n1) + END SUBROUTINE exchange_valence_hf diff --git a/hybrid/gen_wavf.F90 b/hybrid/gen_wavf.F90 index aa32bd27..9cd5fd87 100644 --- a/hybrid/gen_wavf.F90 +++ b/hybrid/gen_wavf.F90 @@ -28,10 +28,10 @@ USE m_abcof USE m_trafo ,ONLY: waveftrafo_genwavf USE m_util ,ONLY: modulo1 - USE m_setabc1locdn USE m_olap USE m_types USE m_abcrot + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(INOUT) :: hybdat @@ -65,16 +65,13 @@ TYPE(t_zmat),INTENT(IN) :: zmat(:) !for all kpoints ! - - local scalars - - - INTEGER :: ilo,idum ,m,maxlmindx - REAL :: rdum,merror - COMPLEX :: cdum,cdum1,cdum2 + INTEGER :: ilo,idum ,m + COMPLEX :: cdum ! local scalars for apws INTEGER :: nred - INTEGER :: ikpt0,ikpt,ikpt1,iband,itype,iop,& - & ispin,ieq,ic,indx,iatom - INTEGER :: i,j,l ,ll,lm,nrkpt,ng,ok - REAL :: time1,time2,time3,time4 + INTEGER :: ikpt0,ikpt,itype,iop,ispin,ieq,indx,iatom + INTEGER :: i,j,l ,ll,lm,ng,ok COMPLEX :: img=(0d0,1d0) ! local scalars for radfun @@ -89,7 +86,6 @@ INTEGER :: rrot(3,3,sym%nsym) INTEGER :: map_lo(atoms%nlod) INTEGER :: iarr(0:atoms%lmaxd,atoms%ntype) - REAL :: rotkpt(3) COMPLEX,ALLOCATABLE :: acof(:,:,:),bcof(:,:,:),ccof(:,:,:,:) REAL,ALLOCATABLE :: zhlp_r(:,:) COMPLEX,ALLOCATABLE :: zhlp_c(:,:) @@ -107,38 +103,21 @@ ! local arrays for radflo REAL :: flo(atoms%jmtd,2,atoms%nlod) - REAL :: uloulopn(atoms%nlod,atoms%nlod,atoms%ntype) REAL :: uuilon(atoms%nlod,atoms%ntype),duilon(atoms%nlod,atoms%ntype) REAL :: ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype) - -! local arrays for setabc1locdn - INTEGER :: nbasf0(atoms%nlod,atoms%nat),nkvec(atoms%nlod,atoms%nat) - INTEGER :: kvec(2*(2*atoms%llod+1) ) - LOGICAL :: enough(atoms%nat) - + ! local arrays for apws INTEGER :: matind(dimension%nbasfcn,2) - INTEGER :: gpthlp1(3,dimension%nvd,dimension%jspd),nvhlp1(dimension%jspd) - INTEGER :: gpthlp2(3,dimension%nvd,dimension%jspd),nvhlp2(dimension%jspd) - INTEGER :: k1hlp(dimension%nvd,dimension%jspd),k2hlp(dimension%nvd,dimension%jspd),& - & k3hlp(dimension%nvd,dimension%jspd),nvhlp(dimension%jspd) - REAL :: bkpt(3) - - REAL,ALLOCATABLE :: olapmt(:,:,:,:) -#if ( defined(CPP_INVERSION) ) - REAL,ALLOCATABLE :: olappw(:,:) -#else - COMPLEX,ALLOCATABLE :: olappw(:,:) -#endif + REAL :: bkpt(3) + INTEGER :: irecl_cmt,irecl_z - INTEGER :: gpt(3,dimension%nvd,dimension%jspd,kpts%nkpt),ngpt(dimension%jspd,kpts%nkpt) - TYPE(t_lapw) :: lapw + TYPE(t_lapw) :: lapw(kpts%nkpt) TYPE(t_usdus):: usdus - CALL CPU_TIME(time1) + !CALL CPU_TIME(time1) call usdus%init(atoms,dimension%jspd) ! setup rotations in reciprocal space @@ -155,7 +134,7 @@ DO ikpt=1,kpts%nkpt CALL apws(dimension,input,noco,& & kpts,ikpt,cell,sym%zrfs,& - & 1,jsp,bkpt,lapw,matind,nred) + & 1,jsp,bkpt,lapw(ikpt),matind,nred) END DO @@ -274,10 +253,7 @@ END IF - maxlmindx = maxval(& - & (/ ( sum( (/ (hybrid%nindx(l,itype)*(2*l+1),l=0,atoms%lmax(itype)) /) ),& - & itype=1,atoms%ntype) /) ) - + ! calculate wavefunction expansion in the the MT region ! (acof,bcof,ccof) and APW-basis coefficients ! (a,b,bascofold_lo) at irred. kpoints @@ -291,9 +267,9 @@ IF( ok .ne. 0 ) STOP 'gen_wavf: failure allocation bcof' ALLOCATE( ccof(-atoms%llod:atoms%llod,dimension%neigd,atoms%nlod,atoms%nat),stat=ok ) IF( ok .ne. 0 ) STOP 'gen_wavf: failure allocation ccof' - ALLOCATE ( cmt(dimension%neigd,maxlmindx,atoms%nat), stat=ok) + ALLOCATE ( cmt(dimension%neigd,hybrid%maxlmindx,atoms%nat), stat=ok) IF( ok .ne. 0 ) STOP 'gen_wavf: Failure allocation cmt' - ALLOCATE ( cmthlp(dimension%neigd,maxlmindx,atoms%nat), stat=ok) + ALLOCATE ( cmthlp(dimension%neigd,hybrid%maxlmindx,atoms%nat), stat=ok) IF( ok .ne. 0) STOP 'gen_wavf: failure allocation cmthlp' if (zmat(1)%l_real) THEN ALLOCATE ( zhlp_r(dimension%nbasfcn,dimension%neigd), stat=ok) @@ -302,9 +278,7 @@ ENDIF IF( ok .ne. 0) STOP 'gen_wavf: failure allocation zhlp' - irecl_cmt = dimension%neigd*maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct',& - & recl=irecl_cmt) + # ifdef CPP_INVERSION irecl_z = dimension%nbasfcn*dimension%neigd*8 @@ -320,9 +294,10 @@ ! abcof calculates the wavefunction coefficients ! stored in acof,bcof,ccof + lapw(ikpt0)%nmat=lapw(ikpt0)%nv(jsp)+atoms%nlotot CALL abcof(& - input,atoms,hybdat%nbands(ikpt0),sym, cell, Kpts%bk(:,ikpt0), lapw, & - ngpt(jsp,ikpt0)+hybdat%nbands(ikpt0),usdus,noco,jsp,hybdat%kveclo_eig(:,ikpt0),& + input,atoms,hybdat%nbands(ikpt0),sym, cell, Kpts%bk(:,ikpt0), lapw(ikpt0), & + hybdat%nbands(ikpt0),usdus,noco,jsp,hybdat%kveclo_eig(:,ikpt0),& oneD,acof(: hybdat%nbands(ikpt0),:,:),bcof(: hybdat%nbands(ikpt0),:,:),ccof(:,: hybdat%nbands(ikpt0),:,:),& zmat(ikpt)) @@ -341,7 +316,7 @@ ! rotate them in the global one CALL abcrot(& - atoms,hybdat%nbands(ikpt0),& + hybrid,atoms,hybdat%nbands(ikpt0),& sym,& cell,oneD,& acof(: hybdat%nbands(ikpt0),:,:),bcof(: hybdat%nbands(ikpt0),:,:),& @@ -407,7 +382,7 @@ ! CALL cpu_time(time3) ! write cmt at irreducible k-points in direct-access file cmt - WRITE(777,rec=ikpt0) cmt + call write_cmt(cmt,ikpt0) ! write z at irreducible k-points in direct-access file z IF (zmat(1)%l_real) THEN @@ -431,11 +406,9 @@ & cmt(:,:,:),zmat(1)%l_real,zmat(ikpt0)%z_r(:,:),zmat(ikpt0)%z_c(:,:),ikpt0,iop,atoms,& & hybrid,kpts,sym,& & jsp,dimension,hybdat%nbands(ikpt0),& - & cell,gpt(:,:ngpt(jsp,ikpt0),jsp,ikpt0),& - & ngpt(:,ikpt0),gpt(:,:ngpt(jsp,ikpt),jsp,ikpt),& - & ngpt(:,ikpt),.true.) + & cell,lapw(ikpt0),lapw(ikpt),.true.) - WRITE(777,rec=ikpt) cmthlp + call write_cmt(cmthlp,ikpt) IF (zmat(1)%l_real) THEN WRITE(778,rec=ikpt) zhlp_r ELSE @@ -449,7 +422,6 @@ DEALLOCATE( cmt,cmthlp) !close file cmt and z - CLOSE(777) CLOSE(778) END SUBROUTINE gen_wavf diff --git a/hybrid/hsfock.F90 b/hybrid/hsfock.F90 index 34573b12..f306f8ed 100644 --- a/hybrid/hsfock.F90 +++ b/hybrid/hsfock.F90 @@ -58,6 +58,7 @@ MODULE m_hsfock USE m_hybridmix USE m_icorrkeys USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat TYPE(t_results),INTENT(INOUT) :: results @@ -85,15 +86,11 @@ MODULE m_hsfock REAL,INTENT(IN) :: eig_irr(dimension%neigd,kpts%nkpt) TYPE(t_hamovlp),INTENT(INOUT)::hamovlp -#if ( !defined(CPP_INVERSION) ) - COMPLEX,ALLOCATABLE :: z(:,:) -#else - REAL,ALLOCATABLE :: z(:,:) -#endif + ! - local scalars - - INTEGER :: i,j,ic,ic1,l,itype + INTEGER :: i,j,ic,ic1,l,itype,n,nn INTEGER :: iband,iband1,iband2 INTEGER :: ikpt,ikpt0 INTEGER :: irec @@ -119,22 +116,7 @@ MODULE m_hsfock REAL :: wl_iks(dimension%neigd,kpts%nkptf) REAL :: div_vv(hybdat%nbands(nk)) - -#ifdef CPP_INVERSION - REAL,ALLOCATABLE :: olap(:,:),olap_p(:) -#else - COMPLEX,ALLOCATABLE :: olap(:,:),olap_p(:) -#endif - -#ifdef CPP_INVERSION - REAL ,ALLOCATABLE :: trafo(:,:),invtrafo(:,:) - REAL ,ALLOCATABLE :: ex(:,:),v(:,:),v_x(:) - REAL ,ALLOCATABLE :: mat_ex(:) -#else - COMPLEX,ALLOCATABLE :: trafo(:,:),invtrafo(:,:) - COMPLEX,ALLOCATABLE :: ex(:,:),v(:,:),v_x(:) - COMPLEX,ALLOCATABLE :: mat_ex(:) -#endif + TYPE(t_mat) :: olap,trafo,invtrafo,ex,tmp,v_x,z COMPLEX :: exch(dimension%neigd,dimension%neigd) COMPLEX,ALLOCATABLE :: carr(:) COMPLEX,ALLOCATABLE :: rep_c(:,:,:,:,:) @@ -149,7 +131,7 @@ MODULE m_hsfock ! ! initialize 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) ! ! initialize weighting factor for HF exchange part @@ -169,49 +151,14 @@ MODULE m_hsfock ! write k1,k2,k3 in gpt DO i=1,lapw%nv(jsp) - gpt(1,i) = lapw%k1(i,jsp) - gpt(2,i) = lapw%k2(i,jsp) - gpt(3,i) = lapw%k3(i,jsp) + gpt(:,i) = (/lapw%k1(i,jsp),lapw%k2(i,jsp),lapw%k3(i,jsp)/) END DO + ! read in lower triangle part of overlap matrix from direct acces file olap - ALLOCATE( olap_p(dimension%nbasfcn*(dimension%nbasfcn+1)/2), stat=ok) - IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation olap_p' - olap_p = 0 - -#ifdef CPP_INVERSION - irecl_olap = dimension%nbasfcn*(dimension%nbasfcn+1)*4 -#else - irecl_olap = dimension%nbasfcn*(dimension%nbasfcn+1)*8 -#endif - irec = kpts%nkpt*(jsp-1) + nk - print *, "Olap read:",irec - OPEN(88,file='olap',form='unformatted',access='direct',& - & recl=irecl_olap) - READ(88,rec=irec) olap_p - CLOSE(88) - - !unpack olap_p into olap - ALLOCATE( olap(dimension%nbasfcn,dimension%nbasfcn), stat=ok) - IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation olap' - olap = 0 - ic = 0 - DO i=1,dimension%nbasfcn - DO j=1,i - ic = ic + 1 - olap(i,j) = olap_p(ic) -#ifdef CPP_INVERSION - olap(j,i) = olap(i,j) -#else - olap(j,i) = conjg(olap_p(ic)) -#endif - END DO - END DO - DEALLOCATE(olap_p) - - ALLOCATE( mat_ex(dimension%nbasfcn*(dimension%nbasfcn+1)),stat=ok ) - IF ( ok .ne. 0) STOP 'mhsfock: failure allocation mat_ex' - mat_ex = 0 + call olap%alloc(hamovlp%l_real,dimension%nbasfcn) + call read_olap(olap, kpts%nkpt*(jsp-1) + nk) + if (.not.olap%l_real) olap%data_c=conjg(olap%data_c) IF( hybrid%l_calhf ) THEN ncstd = sum( (/ ( (hybdat%nindxc(l,itype)*(2*l+1)*atoms%neq(itype),& @@ -258,7 +205,7 @@ MODULE m_hsfock & it,xcpot,& & noco,nsest,indx_sest,& & mpi,irank2,isize2,comm,& - & div_vv,mat_ex) + & div_vv,ex) DEALLOCATE ( rep_c ) CALL timestop("valence exchange calculation") @@ -295,7 +242,7 @@ MODULE m_hsfock & lapw,& & nsymop,nsest,indx_sest,mpi,& & a_ex,results,& - & mat_ex) + & ex) CALL exchange_cccc(& & nk,atoms,hybdat,& & ncstd,& @@ -311,75 +258,28 @@ MODULE m_hsfock IF( dimension%neigd .lt. hybdat%nbands(nk) ) STOP 'mhsfock: neigd < nbands(nk) ; && &trafo from wavefunctions to APW requires at least nbands(nk) ' - ALLOCATE(z(dimension%nbasfcn,dimension%neigd),stat=ok) - IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation z' - z = 0 -#ifdef CPP_INVERSION - irecl_z = dimension%nbasfcn*dimension%neigd*8 -#else - irecl_z = dimension%nbasfcn*dimension%neigd*16 -#endif - - OPEN(unit=778,file='z',form='unformatted',access='direct',& - & recl=irecl_z) - READ(778,rec=nk) z - CLOSE(778) - - !unpack mat_ex - ALLOCATE( ex(hybdat%nbands(nk),hybdat%nbands(nk)), stat = ok ) - IF( ok .ne. 0 ) STOP 'mhsfock: error allocation ex' - ex = 0 - - ic = 0 - DO i = 1,hybdat%nbands(nk) - DO j = 1,i - ic = ic + 1 - ex(j,i) = mat_ex(ic) -#ifdef CPP_INVERSION - ex(i,j) = ex(j,i) -#else - ex(i,j) = conjg(mat_ex(ic)) -#endif - END DO - END DO - + call z%alloc(olap%l_real,dimension%nbasfcn,dimension%neigd) + call read_z(z,nk) !what about spin? + ! calculate trafo ic = lapw%nv(jsp) + atoms%nlotot - - ALLOCATE( trafo(ic,hybdat%nbands(nk)), stat= ok ) - IF( ok .ne. 0 ) STOP 'mhsfock: error allocation trafo' - trafo = matmul(olap(:ic,:ic),z(:ic,:hybdat%nbands(nk))) - - ALLOCATE( invtrafo(hybdat%nbands(nk),ic), stat= ok ) - IF( ok .ne. 0 ) STOP 'mhsfock: error allocation invtrafo' - -#ifdef CPP_INVERSION - invtrafo = transpose(trafo) -#else - invtrafo = conjg(transpose(trafo)) -#endif - - ALLOCATE( v(ic,ic), stat = ok ) - IF( ok .ne. 0 ) STOP 'mhsfock: error allocation v' - - v = matmul(trafo,matmul(ex,invtrafo)) - - DEALLOCATE( mat_ex,ex,trafo,invtrafo,olap ) + z%matsize1=ic + z%matsize2=hybdat%nbands(nk) + olap%matsize1=ic + olap%matsize2=ic + + call olap%multiply(z,trafo) + + call invtrafo%alloc(olap%l_real,hybdat%nbands(nk),ic) + CALL trafo%transpose(invtrafo) + + CALL ex%multiply(invtrafo,tmp) + CALL trafo%multiply(tmp,v_x) + CALL timestop("time for performing T^-1*mat_ex*T^-1*") - !store only lower triangle of v - ALLOCATE( v_x(dimension%nbasfcn*(dimension%nbasfcn+1)/2), stat = ok ) - IF( ok .ne. 0 ) STOP 'mhsfock: error allocation v_x' - v_x = 0 - ic = 0 - DO i=1,lapw%nv(jsp)+atoms%nlotot - DO j = 1,i - ic = ic + 1 - v_x(ic) = v(i,j) - END DO - END DO - + CALL symmetrizeh(atoms,& & kpts%bkf(:,nk),dimension,jsp,lapw,gpt,& & sym,hybdat%kveclo_eig,& @@ -389,84 +289,45 @@ MODULE m_hsfock IF( input%imix .ge. 10 ) THEN -#ifdef CPP_INVERSION - irecl_vx = dimension%nbasfcn*(dimension%nbasfcn+1)*4 -#else - irecl_vx = dimension%nbasfcn*(dimension%nbasfcn+1)*8 -#endif - - irec = kpts%nkpt*(jsp-1) + nk - OPEN(778,file='vex',form='unformatted',access='direct',& - & recl=irecl_vx) -#ifdef CPP_INVERSION - WRITE(778,rec=irec) v_x -#else - WRITE(778,rec=irec) v_x -#endif - CLOSE(778) + CALL write_v_x(v_x,kpts%nkpt*(jsp-1) + nk) END IF ELSE ! not hybrid%l_calhf - - ALLOCATE(z(dimension%nbasfcn,dimension%neigd),stat=ok) - IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation z' - z = 0 -#ifdef CPP_INVERSION - irecl_z = dimension%nbasfcn*dimension%neigd*8 -#else - irecl_z = dimension%nbasfcn*dimension%neigd*16 -#endif - - OPEN(unit=778,file='z',form='unformatted',access='direct',& - & recl=irecl_z) - READ(778,rec=nk) z - CLOSE(778) - - ALLOCATE( v_x(dimension%nbasfcn*(dimension%nbasfcn+1)/2), stat = ok ) - IF( ok .ne. 0 ) STOP 'mhsfock: error allocation v_x' - v_x = 0 -#ifdef CPP_INVERSION - irecl_vx = dimension%nbasfcn*(dimension%nbasfcn+1)*4 -#else - irecl_vx = dimension%nbasfcn*(dimension%nbasfcn+1)*8 -#endif - irec = kpts%nkpt*(jsp-1) + nk - OPEN(778,file='vex',form='unformatted',access='direct',& - & recl=irecl_vx) - READ(778,rec=irec) v_x - CLOSE(778) - - END IF ! hybrid%l_calhf - - ! add non-local x-potential to the hamiltonian a - IF (hamovlp%l_real) THEN - hamovlp%a_r = hamovlp%a_r - a_ex*v_x - ELSE - hamovlp%a_c = hamovlp%a_c - a_ex*v_x - ENDIF + CALL read_v_x(v_x,kpts%nkpt*(jsp-1) + nk) + END IF ! hybrid%l_calhf + ! add non-local x-potential to the hamiltonian a (in packed storage) + ic=1 + DO n=1,v_x%matsize1 + DO nn=1,n + IF (hamovlp%l_real) THEN + hamovlp%a_r(ic) = hamovlp%a_r(ic) - a_ex*v_x%data_r(n,nn) + ELSE + hamovlp%a_c(ic) = hamovlp%a_c(ic) - a_ex*v_x%data_c(n,nn) + ENDIF + ENDDO + END DO ! calculate HF energy IF( hybrid%l_calhf ) THEN - WRITE(6,'(A)') new_line('n')//new_line('n')//' ### '// ' diagonal HF exchange elements (eV) ###' + WRITE(6,'(A)') new_line('n')//new_line('n')//' ### '// ' diagonal HF exchange elements (eV) ###' - WRITE(6,'(A)') new_line('n') // ' k-point '// 'band tail pole input%total(valence+core)' + WRITE(6,'(A)') new_line('n') // ' k-point '// 'band tail pole input%total(valence+core)' END IF - ic = lapw%nv(jsp) + atoms%nlotot - ic1 = ic*(ic+1)/2 - -#if( !defined CPP_INVERSION ) - v_x(:ic1) = conjg(v_x(:ic1)) -#endif - + + ! calculate exchange contribution of current k-point nk to total energy (te_hfex) ! in the case of a spin-unpolarized calculation the factor 2 is added in eigen_hf.F - + if (.not.v_x%l_real) v_x%data_c=conjg(v_x%data_c) exch = 0 + call v_x%multiply(z,tmp) DO iband = 1,hybdat%nbands(nk) - exch(iband,iband) = dotprod(z(:ic,iband),matvec(v_x(:ic1),z(:ic,iband))) + if (z%l_real) THEN + exch(iband,iband) = dot_product(z%data_r(:,iband),tmp%data_r(:,iband)) + else + exch(iband,iband) = dot_product(z%data_r(:,iband),tmp%data_r(:,iband)) + endif IF( iband .le. hybdat%nobd(nk) ) THEN - results%te_hfex%valence = results%te_hfex%valence -& - & a_ex*results%w_iks(iband,nk,jsp)*exch(iband,iband) + results%te_hfex%valence = results%te_hfex%valence -a_ex*results%w_iks(iband,nk,jsp)*exch(iband,iband) END IF IF(hybrid%l_calhf) THEN WRITE(6, '( '' ('',F5.3,'','',F5.3,'','',F5.3,'')'',I4,4X,3F10.5)')& @@ -476,8 +337,6 @@ MODULE m_hsfock END DO - DEALLOCATE( z,v_x ) - CALL timestop("total time hsfock") END SUBROUTINE hsfock diff --git a/hybrid/kp_perturbation.F90 b/hybrid/kp_perturbation.F90 index e7e32999..e460ce96 100644 --- a/hybrid/kp_perturbation.F90 +++ b/hybrid/kp_perturbation.F90 @@ -96,7 +96,7 @@ #endif - 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) @@ -811,6 +811,7 @@ USE m_dr2fdr USE m_constants USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat @@ -868,11 +869,9 @@ #endif ! read in cmt coefficients from direct access file cmt at kpoint nk - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct', recl=irecl_cmt) - READ(777,rec=nk) cmt - CLOSE(777) - + + call read_cmt(cmt,nk) + ! read in z coefficients from direct access file z at kpoint nk #ifdef CPP_INVERSION irecl_z = dimension%nbasfcn*dimension%neigd*8 @@ -884,7 +883,7 @@ CLOSE(778) - 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) gpt(1,1:lapw%nv(jsp)) = lapw%k1(1:lapw%nv(jsp),jsp) gpt(2,1:lapw%nv(jsp)) = lapw%k2(1:lapw%nv(jsp),jsp) diff --git a/hybrid/symm_hf.F90 b/hybrid/symm_hf.F90 index eff3dd7f..4149f676 100644 --- a/hybrid/symm_hf.F90 +++ b/hybrid/symm_hf.F90 @@ -31,6 +31,7 @@ USE m_olap ,ONLY: wfolap,wfolap1,wfolap_init USE m_trafo ,ONLY: waveftrafo_symm USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat @@ -308,12 +309,7 @@ ! P(R,T)\phi_n,k = \sum_{n'\in degenerat(n)} rep_v(n',n) *\phi_n',k ! read in cmt and z at current k-point (nk) - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct',& - & recl=irecl_cmt) - READ(777,rec=nk) cmt - CLOSE(777) - + CALL read_cmt(cmt,nk) #ifdef CPP_INVERSION irecl_z = dimension%nbasfcn*dimension%neigd*8 #else @@ -465,13 +461,9 @@ ELSE ! read in cmt and z at current k-point (nk) - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct',& - & recl=irecl_cmt) - READ(777,rec=nk) cmt - CLOSE(777) - - CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf) + + CALL read_cmt(cmt,nk) + !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf) IF( allocated(olapmt)) deallocate(olapmt) ALLOCATE( olapmt(hybrid%maxindx,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype),stat=ok) @@ -622,7 +614,7 @@ & exp( 2*pi*img*dot_product(g,atoms%taual(:,ratom))) rep_c(:,:,:,iop,iatom) = & - & sym%d_wgn(-hybdat%lmaxcd:hybdat%lmaxcd,-hybdat%lmaxcd:hybdat%lmaxcd,0:hybdat%lmaxcd,isym) * cdum + & hybrid%d_wgn2(-hybdat%lmaxcd:hybdat%lmaxcd,-hybdat%lmaxcd:hybdat%lmaxcd,0:hybdat%lmaxcd,isym) * cdum END DO END DO iatom0 = iatom0 + atoms%neq(itype) diff --git a/hybrid/symmetrizeh.F90 b/hybrid/symmetrizeh.F90 index 9a8659c2..87e9670f 100644 --- a/hybrid/symmetrizeh.F90 +++ b/hybrid/symmetrizeh.F90 @@ -9,7 +9,7 @@ & bk,dimension,jsp,lapw,gpt,& & sym,kveclo,cell,& & nsymop,psym,& - & a) + & hmat) USE m_constants @@ -30,12 +30,7 @@ INTEGER,INTENT(IN) :: psym(nsymop) REAL,INTENT(IN) :: bk(3) - -#if ( defined(CPP_INVERSION) && !defined(CPP_SOC) ) - REAL,INTENT(INOUT) :: a(dimension%nbasfcn*(dimension%nbasfcn+1)/2) -#else - COMPLEX,INTENT(INOUT) :: a(dimension%nbasfcn*(dimension%nbasfcn+1)/2) -#endif + TYPE(T_mat),INTENT(INOUT) :: hmat ! - local scalars - INTEGER :: ilotot,itype,itype1,ilo,ilo1 @@ -67,12 +62,6 @@ REAL :: rtaual(3),kghlp(3) REAL :: rotkpthlp(3),rotkpt(3) REAL :: trans(3,nsymop) - -#ifdef CPP_INVERSION - REAL :: v(dimension%nbasfcn,dimension%nbasfcn) -#else - COMPLEX :: v(dimension%nbasfcn,dimension%nbasfcn) -#endif COMPLEX,ALLOCATABLE :: c_lo(:,:,:,:),c_rot(:,:,:,:,:),y(:) COMPLEX,ALLOCATABLE :: cfac(:,:),chelp(:,:) @@ -169,22 +158,6 @@ END DO - ! unpack a into v and initialize indx - ic = 0 - DO j = 1,lapw%nv(jsp)+atoms%nlotot - DO i = 1,j - ic = ic + 1 - indx(j,i) = ic - indx(i,j) = ic - v(j,i) = a(ic) -#ifdef CPP_INVERSION - v(i,j) = v(j,i) -#else - v(i,j) = conjg(v(j,i)) -#endif - END DO - END DO - @@ -203,13 +176,21 @@ IF( iop .le. sym%nop) THEN IF( igpt .ne. 0 .and. igpt1 .ne. 0 ) THEN - ic = ic + 1 - cdum = cdum + conjg(cfac(i,isym))*v(igpt,igpt1)*cfac(j,isym) + ic = ic + 1 + if (hmat%l_real) THEN + cdum = cdum + conjg(cfac(i,isym))*hmat%data_r(igpt,igpt1)*cfac(j,isym) + else + cdum = cdum + conjg(cfac(i,isym))*hmat%data_c(igpt,igpt1)*cfac(j,isym) + endif END IF ELSE IF( igpt .ne. 0 .and. igpt1 .ne. 0 ) THEN - ic = ic + 1 - cdum = cdum + conjg(conjg(cfac(i,isym))*v(igpt,igpt1)*cfac(j,isym)) + ic = ic + 1 + if (hmat%l_real) THEN + cdum = cdum + conjg(conjg(cfac(i,isym))*hmat%data_r(igpt,igpt1)*cfac(j,isym)) + else + cdum = cdum + conjg(conjg(cfac(i,isym))*hmat%data_c(igpt,igpt1)*cfac(j,isym)) + end if END IF END IF @@ -224,13 +205,23 @@ IF( igpt .eq. 0 .or. igpt1 .eq. 0 ) CYCLE IF( igpt1 .gt. igpt ) CYCLE IF( ldum(igpt,igpt1) ) THEN - IF( iop .le. sym%nop) THEN - a(indx(igpt,igpt1)) = cdum/(conjg(cfac(i,isym))*cfac(j,isym)) - ldum(igpt,igpt1) = .false. - ELSE - a(indx(igpt,igpt1)) = conjg(cdum/(conjg(cfac(i,isym))*cfac(j,isym))) - ldum(igpt,igpt1) = .false. - END IF + if (hmat%l_real) THEN + IF( iop .le. sym%nop) THEN + hmat%data_r(igpt,igpt1) = cdum/(conjg(cfac(i,isym))*cfac(j,isym)) + ldum(igpt,igpt1) = .false. + ELSE + hmat%data_r(igpt,igpt1) = conjg(cdum/(conjg(cfac(i,isym))*cfac(j,isym))) + ldum(igpt,igpt1) = .false. + END IF + else + IF( iop .le. sym%nop) THEN + hmat%data_c(igpt,igpt1) = cdum/(conjg(cfac(i,isym))*cfac(j,isym)) + ldum(igpt,igpt1) = .false. + ELSE + hmat%data_c(igpt,igpt1) = conjg(cdum/(conjg(cfac(i,isym))*cfac(j,isym))) + ldum(igpt,igpt1) = .false. + END IF + end if END IF END DO END DO @@ -459,8 +450,12 @@ cdum2 = 0 ic1 = 0 DO igpt2 = igpt_lo1,igpt_lo2 - ic1 = ic1 + 1 - cdum2 = cdum2 + conjg(c_rot(ic1,igpt,ilo,ratom,isym)) * v(lapw%nv(jsp)+igpt2,igpt1) + ic1 = ic1 + 1 + if (hmat%l_real) THEN + cdum2 = cdum2 + conjg(c_rot(ic1,igpt,ilo,ratom,isym)) * hmat%data_r(lapw%nv(jsp)+igpt2,igpt1) + else + cdum2 = cdum2 + conjg(c_rot(ic1,igpt,ilo,ratom,isym)) * hmat%data_c(lapw%nv(jsp)+igpt2,igpt1) + end if END DO IF( iop .le. sym%nop) THEN @@ -469,8 +464,11 @@ cdum = cdum + conjg(cdum2*conjg(cfac(lapw%nv(jsp)+i,isym)) * cfac(j,isym)) END IF END DO - - a(indx(lapw%nv(jsp)+i,j)) = cdum/ic + if (hmat%l_real) THEN + hmat%data_r(lapw%nv(jsp)+i,j) = cdum/ic + else + hmat%data_c(lapw%nv(jsp)+i,j) = cdum/ic + end if END DO END DO END DO @@ -554,22 +552,29 @@ ic1 = ic1 + 1 ic2 = 0 DO igpt3 = igpt1_lo1,igpt1_lo2 - ic2 = ic2 + 1 - cdum2 = cdum2 + conjg(c_rot(ic1,igpt,ilo,ratom,isym))& - & * v(lapw%nv(jsp)+igpt2,lapw%nv(jsp)+igpt3) * c_rot(ic2,igpt1,ilo1,ratom1,isym) - END DO - END DO -! - ic = ic + 1 + ic2 = ic2 + 1 + if (hmat%l_real) THEN + cdum2 = cdum2 + conjg(c_rot(ic1,igpt,ilo,ratom,isym))& + * hmat%data_r(lapw%nv(jsp)+igpt2,lapw%nv(jsp)+igpt3) * c_rot(ic2,igpt1,ilo1,ratom1,isym) + else + cdum2 = cdum2 + conjg(c_rot(ic1,igpt,ilo,ratom,isym))& + * hmat%data_c(lapw%nv(jsp)+igpt2,lapw%nv(jsp)+igpt3) * c_rot(ic2,igpt1,ilo1,ratom1,isym) + end if + end DO + end DO + ic = ic + 1 IF( iop .le. sym%nop ) THEN - cdum = cdum + cdum2 *conjg(cfac(lapw%nv(jsp)+i,isym))*cfac(lapw%nv(jsp)+j,isym) - ELSE + cdum = cdum + cdum2*conjg(cfac(lapw%nv(jsp)+i,isym))*cfac(lapw%nv(jsp)+j,isym) + ELSE cdum = cdum + conjg(cdum2*conjg(cfac(lapw%nv(jsp)+i,isym))*cfac(lapw%nv(jsp)+j,isym)) - END IF - END DO ! isym - - a(indx(lapw%nv(jsp)+i,lapw%nv(jsp)+j)) = cdum/ic - + + END IF + END DO + if (hmat%l_real) THEN + hmat%data_r(lapw%nv(jsp)+i,lapw%nv(jsp)+j) = cdum!/ic + else + hmat%data_c(lapw%nv(jsp)+i,lapw%nv(jsp)+j) = cdum!/ic + end if END DO ! igpt_lo1 END DO ! ilo1 diff --git a/hybrid/wavefproducts.F90 b/hybrid/wavefproducts.F90 index dc6b029e..8168b506 100644 --- a/hybrid/wavefproducts.F90 +++ b/hybrid/wavefproducts.F90 @@ -15,6 +15,7 @@ USE m_util , ONLY : modulo1 USE m_wrapper USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat @@ -136,9 +137,7 @@ READ(778,rec=nk) z_nk ! read in cmt coefficients from direct access file cmt - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct', recl=irecl_cmt) - READ(777,rec=nk) cmt_nk(:,:,:) + call read_cmt(cmt_nk,nk) ! ! open file cprod ! irecl_cprod = maxbasm*mnobd*(bandf-bandi+1)*16 @@ -195,8 +194,7 @@ ! ! MT contribution ! - - READ(777,rec=nkqpt) cmt(:,:,:) + call read_cmt(cmt,nkqpt) lm_0 = 0 ic = 0 @@ -301,7 +299,6 @@ END DO - CLOSE(777) CLOSE(778) @@ -321,6 +318,7 @@ USE m_wrapper USE m_constants USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat TYPE(t_dimension),INTENT(IN) :: dimension @@ -476,15 +474,12 @@ END DO ! read in cmt coefficient at k-point nk - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct', recl=irecl_cmt) - READ(777,rec=nk) ccmt_nk + + CALL read_cmt(ccmt_nk,nk) !read in cmt coefficients at k+q point - READ(777,rec=nkqpt) ccmt(:,:,:) + call read_cmt(ccmt,nkqpt) - CLOSE(777) - iatom = 0 DO itype = 1,atoms%ntype @@ -1276,6 +1271,7 @@ USE m_constants USE m_apws USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat TYPE(t_dimension),INTENT(IN) :: dimension @@ -1505,15 +1501,10 @@ ALLOCATE(ccmt_nk(dimension%neigd,hybrid%maxlmindx,atoms%nat),ccmt(dimension%neigd,hybrid%maxlmindx,atoms%nat), stat=ok) IF( ok .ne. 0 ) STOP 'wavefproducts_inv5: error allocation ccmt_nk/ccmt' - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct', recl=irecl_cmt) - READ(777,rec=nk) ccmt_nk - + call read_cmt(ccmt_nk,nk) !read in cmt coefficients at k+q point - READ(777,rec=nkqpt) ccmt(:,:,:) + call read_cmt(ccmt,nkqpt) - CLOSE(777) - iatom = 0 DO itype = 1,atoms%ntype @@ -2317,6 +2308,7 @@ USE m_trafo USE m_wrapper USE m_types + USE m_io_hybrid IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat TYPE(t_dimension),INTENT(IN) :: dimension @@ -2545,12 +2537,9 @@ ! read in cmt coefficients from direct access file cmt - irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16 - OPEN(unit=777,file='cmt',form='unformatted',access='direct', recl=irecl_cmt) - READ(777,rec=nk) cmt_nk(:,:,:) - READ(777,rec=nkqpt) cmt(:,:,:) - CLOSE(777) - + call read_cmt(cmt_nk(:,:,:),nk) + call read_cmt(cmt(:,:,:),nkqpt) + lm_0 = 0 ic = 0 diff --git a/init/gen_bz.F90 b/init/gen_bz.F90 index 38301e52..176f9dcf 100644 --- a/init/gen_bz.F90 +++ b/init/gen_bz.F90 @@ -76,8 +76,17 @@ SUBROUTINE gen_bz( kpts,sym) ! add it if it is not yet in this set. kpts%bkf = 0 - ic = 0 - + + !Add existing vectors to list of full vectors + print *,"WARNING from gen_bz" + print *,"Assuming Identity to be fist symmetry op!" + DO ic=1,kpts%nkpt + kpts%bkf(:,ic) = kpts%bk(:,ic) + kpts%bkp(ic) = ic + kpts%bksym(ic) = 1 + ENDDO + ic=ic-1 + DO iop=1,nsym DO ikpt=1,kpts%nkpt l_found = .FALSE. diff --git a/init/init_hybrid.f90 b/init/init_hybrid.f90 new file mode 100644 index 00000000..dffee10b --- /dev/null +++ b/init/init_hybrid.f90 @@ -0,0 +1,37 @@ + IF ( (xcpot%icorr.EQ.icorr_hf ) .OR. (xcpot%icorr.EQ.icorr_pbe0)& + & .OR.(xcpot%icorr.EQ.icorr_exx) .OR. (xcpot%icorr.EQ.icorr_hse)& + & .OR.(xcpot%icorr.EQ.icorr_vhse) ) THEN + IF (input%film .OR. oneD%odi%d1)& + & CALL juDFT_error("2D film and 1D calculations not implemented"& + & //"for HF/EXX/PBE0/HSE", calledby ="fleur",& + & hint="Use a supercell or a different functional") + + IF( ANY( atoms%l_geo ) )& + & CALL juDFT_error("Forces not implemented for HF/PBE0/HSE ",& + & calledby ="fleur") + + IF (.NOT. obsolete%pot8) STOP 'Choose pot8=T' + !calculate whole Brilloun zone + CALL gen_bz(kpts,sym) + CALL gen_map(& + & atoms,sym,oneD,hybrid) + ! + ! calculate d_wgn + ! + ALLOCATE (hybrid%d_wgn2(-atoms%lmaxd:atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,sym%nsym)) + CALL d_wigner(sym%nop,sym%mrot,cell%bmat,atoms%lmaxd,hybrid%d_wgn2(:,:,1:,:sym%nop)) + hybrid%d_wgn2(:,:,0,:) = 1 + + DO isym = sym%nop+1,sym%nsym + iisym = isym - sym%nop + DO l = 0,atoms%lmaxd + DO m2 = -l,l + DO m1 = -l,-1 + cdum = hybrid%d_wgn2( m1,m2,l,iisym) + hybrid%d_wgn2( m1,m2,l,isym) = hybrid%d_wgn2(-m1,m2,l,iisym)*(-1)**m1 + hybrid%d_wgn2(-m1,m2,l,isym) = cdum *(-1)**m1 + END DO + hybrid%d_wgn2(0,m2,l,isym) = hybrid%d_wgn2(0,m2,l,iisym) + END DO + END DO + END DO diff --git a/init/kpoints.f90 b/init/kpoints.f90 new file mode 100644 index 00000000..02048995 --- /dev/null +++ b/init/kpoints.f90 @@ -0,0 +1,61 @@ +!-------------------------------------------------------------------------------- +! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany +! This file is part of FLEUR and available as free software under the conditions +! of the MIT license as expressed in the LICENSE file in more detail. +!-------------------------------------------------------------------------------- + +module m_kpoints +contains + subroutine kpoints(oneD,jij,sym,cell,input,noco,banddos,kpts,l_kpts) + use m_types + use m_julia + use m_kptgen_hybrid + use m_od_kptsgen + USE m_gen_bz + + implicit none + TYPE(t_input),INTENT(IN) :: input + TYPE(t_sym),INTENT(IN) :: sym + TYPE(t_oneD),INTENT(IN) :: oneD + TYPE(t_Jij),INTENT(IN) :: Jij + TYPE(t_cell),INTENT(IN) :: cell + TYPE(t_banddos),INTENT(IN) :: banddos + TYPE(t_noco),INTENT(IN) :: noco + TYPE(t_kpts),INTENT(INOUT) :: kpts + LOGICAL,INTENT(IN) :: l_kpts + + TYPE(t_sym):: sym_hlp + + if (.not.l_kpts) THEN + IF (.NOT.oneD%odd%d1) THEN + IF (jij%l_J) THEN + sym_hlp=sym + sym_hlp%nop=1 + sym_hlp%nop2=1 + CALL julia(sym_hlp,cell,input,noco,banddos,kpts,.FALSE.,.TRUE.) + ELSE IF(kpts%l_gamma .and. banddos%ndir .eq. 0) THEN + CALL kptgen_hybrid(kpts,sym%invs,noco%l_soc,sym%nop,sym%mrot,sym%tau) + ELSE + CALL julia(sym,cell,input,noco,banddos,kpts,.FALSE.,.TRUE.) + END IF + ELSE + STOP 'Error: No kpoint set generation for 1D systems yet!' + CALL od_kptsgen (kpts%nkpt) + END IF + endif + !Rescale weights and kpoints + + kpts%wtkpt(:) = kpts%wtkpt(:) / sum(kpts%wtkpt) + kpts%bk(:,:) = kpts%bk(:,:) / kpts%posScale + kpts%posScale = 1.0 + IF (kpts%nkpt3(3).EQ.0) kpts%nkpt3(3) = 1 + + !Calculate kpoint in the full BZ + IF (all(kpts%nkpt3>0)) THEN + CALL gen_bz(kpts,sym) + else + kpts%nkptf=0 + endif + +end subroutine kpoints +end module m_kpoints diff --git a/io/CMakeLists.txt b/io/CMakeLists.txt index cfb28ee3..1ea779e3 100644 --- a/io/CMakeLists.txt +++ b/io/CMakeLists.txt @@ -5,6 +5,8 @@ io/rw_symfile.f io/write_gw.F ) set(fleur_F90 ${fleur_F90} +io/io_matrix.F90 +io/io_hybrid.F90 io/eig66_da.F90 io/eig66_data.F90 io/eig66_hdf.F90 diff --git a/io/eig66_da.F90 b/io/eig66_da.F90 index 8da0b231..dcaa3647 100644 --- a/io/eig66_da.F90 +++ b/io/eig66_da.F90 @@ -197,9 +197,9 @@ CONTAINS IF (PRESENT(wk)) wk=wtkpt IF (PRESENT(nv)) nv=nv_s IF (PRESENT(nmat)) nmat=nmat_s - IF (PRESENT(k1)) k1=d%kvec_s(:,1) - IF (PRESENT(k2)) k2=d%kvec_s(:,2) - IF (PRESENT(k3)) k3=d%kvec_s(:,3) + IF (PRESENT(k1)) k1=d%kvec_s(:size(k1),1) + IF (PRESENT(k2)) k2=d%kvec_s(:size(k1),2) + IF (PRESENT(k3)) k3=d%kvec_s(:size(k1),3) IF (PRESENT(kveclo)) kveclo=d%kveclo_s ENDIF @@ -282,9 +282,9 @@ CONTAINS PRESENT(bk).AND.PRESENT(el).AND.PRESENT(ello).AND.PRESENT(evac))) THEN CALL juDFT_error("BUG:Direct access IO of eig-file only with all data") ENDIF - d%kvec_s(:,1)=k1 - d%kvec_s(:,2)=k2 - d%kvec_s(:,3)=k3 + d%kvec_s(:size(k1),1)=k1 + d%kvec_s(:size(k1),2)=k2 + d%kvec_s(:size(k1),3)=k3 IF ((SIZE(el).NE.SIZE(d%el_s)).OR.(SIZE(ello).NE.SIZE(d%ello_s).OR.(SIZE(evac).NE.SIZE(d%evac_s)))) THEN WRITE(*,*) SHAPE(el),SHAPE(d%el_s) WRITE(*,*) SHAPE(ello),SHAPE(d%ello_s) diff --git a/io/io_hybrid.F90 b/io/io_hybrid.F90 new file mode 100644 index 00000000..701cabe6 --- /dev/null +++ b/io/io_hybrid.F90 @@ -0,0 +1,98 @@ +module m_io_hybrid + use m_io_matrix + use m_judft + use m_types + implicit none + !private + integer,save :: id_olap,id_z,id_v_x + !public:: open_hybrid_io,read_cmt,write_cmt +contains + + subroutine open_hybrid_io(hybrid,dimension,atoms,l_real) + implicit none + TYPE(t_hybrid),INTENT(IN) :: hybrid + TYPE(t_dimension),INTENT(IN):: dimension + TYPE(t_atoms),INTENT(IN) :: atoms + LOGICAL,INTENT(IN) :: l_real + + + OPEN(unit=777,file='cmt',form='unformatted',access='direct',& + & recl=dimension%neigd*hybrid%maxlmindx*atoms%nat*16) + print *,"Open olap.mat" + id_olap=OPEN_MATRIX(l_real,dimension%nbasfcn,1,"olap.mat") + print *,"Open z.mat" + id_z=OPEN_MATRIX(l_real,dimension%nbasfcn,1,"z.mat") + print *,"Open v_x.mat" + id_v_x=OPEN_MATRIX(l_real,dimension%nbasfcn,1,"v_x.mat") + + end subroutine open_hybrid_io + + subroutine write_cmt(cmt,nk) + implicit none + complex,INTENT(IN):: cmt(:,:,:) + integer,INTENT(IN):: nk + + write(777,rec=nk) cmt + end subroutine write_cmt + + subroutine read_cmt(cmt,nk) + implicit none + complex,INTENT(OUT):: cmt(:,:,:) + integer,INTENT(IN):: nk + + read(777,rec=nk) cmt + end subroutine read_cmt + + + subroutine read_olap(mat,rec) + implicit none + TYPE(t_mat),INTENT(INOUT):: mat + INTEGER,INTENT(IN) :: rec + + CALL read_matrix(mat,rec,id_olap) + END subroutine read_olap + + subroutine write_olap(mat,rec) + implicit none + TYPE(t_mat),INTENT(IN) :: mat + INTEGER,INTENT(IN) :: rec + + CALL write_matrix(mat,rec,id_olap) + END subroutine write_olap + + subroutine read_z(mat,rec) + implicit none + TYPE(t_mat),INTENT(INOUT):: mat + INTEGER,INTENT(IN) :: rec + + CALL read_matrix(mat,rec,id_z) + END subroutine read_z + + subroutine write_z(mat,rec) + implicit none + TYPE(t_mat),INTENT(IN) :: mat + INTEGER,INTENT(IN) :: rec + + CALL write_matrix(mat,rec,id_z) + END subroutine write_z + + subroutine read_v_x(mat,rec) + implicit none + TYPE(t_mat),INTENT(INOUT):: mat + INTEGER,INTENT(IN) :: rec + + CALL read_matrix(mat,rec,id_v_x) + END subroutine read_v_x + + subroutine write_v_x(mat,rec) + implicit none + TYPE(t_mat),INTENT(IN) :: mat + INTEGER,INTENT(IN) :: rec + + CALL write_matrix(mat,rec,id_v_x) + END subroutine write_v_x + + + + +end module m_io_hybrid diff --git a/io/io_matrix.F90 b/io/io_matrix.F90 new file mode 100644 index 00000000..d5858597 --- /dev/null +++ b/io/io_matrix.F90 @@ -0,0 +1,101 @@ +MODULE m_io_matrix + USE m_types + IMPLICIT NONE + private + INTEGER:: mode=1 + public:: open_matrix,read_matrix,write_matrix +CONTAINS + INTEGER FUNCTION OPEN_matrix(l_real,matsize,no_rec,filename) + LOGICAL,INTENT(IN) :: l_real + INTEGER,INTENT(in) :: matsize,no_rec + CHARACTER(len=*),INTENT(in) :: filename + SELECT CASE (mode) + CASE (1) + OPEN_matrix=open_DA(l_real,matsize,no_rec,filename) + CASE default + CALL judft_error("BUG in io_matrix") + END SELECT + END FUNCTION OPEN_MATRIX + + SUBROUTINE read_matrix(mat,rec,id) + TYPE(t_Mat),INTENT(INOUT):: mat + INTEGER,INTENT(IN) :: rec,id + + CALL mat%alloc() + SELECT CASE (mode) + CASE (1) + CALL read_matrix_DA(mat,rec,id) + CASE default + CALL judft_error("BUG in io_matrix") + END SELECT + END SUBROUTINE read_matrix + + SUBROUTINE write_matrix(mat,rec,id) + TYPE(t_Mat),INTENT(IN):: mat + INTEGER,INTENT(IN) :: rec,id + + SELECT CASE (mode) + CASE (1) + CALL write_matrix_DA(mat,rec,id) + CASE default + CALL judft_error("BUG in io_matrix") + END SELECT + END SUBROUTINE write_matrix + + !Now the implementation in terms of fortran DA-files + INTEGER FUNCTION open_DA(l_real,matsize,no_rec,filename) + LOGICAL,INTENT(IN) :: l_real + INTEGER,INTENT(in) :: matsize,no_rec + CHARACTER(len=*),INTENT(in) :: filename + + LOGICAL :: used_unit + REAL :: r + COMPLEX :: c + INTEGER :: datasize + + !Determine size of data + IF (l_real) THEN + INQUIRE(IOLENGTH=datasize) r + ELSE + INQUIRE(IOLENGTH=datasize) c + END IF + + !find free unit starting at 901 + open_DA=901 + DO + INQUIRE(unit=open_DA,opened=used_unit) + IF (.NOT.used_unit) EXIT + open_DA=open_DA+1 + END DO + !openfile + OPEN(unit=open_DA,file=filename,access='direct',recl=datasize*matsize*matsize) + + + END FUNCTION open_DA + + SUBROUTINE read_matrix_DA(mat,rec,id) + TYPE(t_Mat),INTENT(INOUT):: mat + INTEGER,INTENT(IN) :: rec,id + + INTEGER:: err + IF (mat%l_real) THEN + READ(id,rec=rec,iostat=err) mat%data_r + ELSE + READ(id,rec=rec,iostat=err) mat%data_c + END IF + IF (err.NE.0) CALL judft_error("Failed in reading of matrix") + END SUBROUTINE read_matrix_DA + + SUBROUTINE write_matrix_DA(mat,rec,id) + TYPE(t_Mat),INTENT(IN):: mat + INTEGER,INTENT(IN) :: rec,id + INTEGER:: err + IF (mat%l_real) THEN + WRITE(id,rec=rec,iostat=err) mat%data_r + ELSE + WRITE(id,rec=rec,iostat=err) mat%data_c + END IF + IF (err.NE.0) CALL judft_error("Failed in writing of matrix") + END SUBROUTINE write_matrix_DA + +END MODULE m_io_matrix diff --git a/math/util.F b/math/util.F index ae1419eb..f42ecff3 100644 --- a/math/util.F +++ b/math/util.F @@ -217,7 +217,7 @@ c Initializes fast numerical integration intgrf (see below). INTEGER, INTENT(IN) :: ntype,jmtd INTEGER, INTENT(IN) :: jri(ntype) REAL, INTENT(IN) :: dx(ntype),rmsh(jmtd,ntype) - REAL, ALLOCATABLE :: gridf(:,:) + REAL, INTENT(OUT), ALLOCATABLE :: gridf(:,:) c - local - INTEGER :: i,j,itype diff --git a/xc-pot/trafo.F90 b/xc-pot/trafo.F90 index abc7b705..722b8d87 100644 --- a/xc-pot/trafo.F90 +++ b/xc-pot/trafo.F90 @@ -160,7 +160,7 @@ SUBROUTINE waveftrafo_genwavf( & cmt_out,z_rout,z_cout,cmt,l_real,z_r,z_c,nk,iop,atoms,& hybrid,kpts,sym,jsp,dimension,nbands,& - cell,gpt_nk,nv_nk,gpt_rkpt,nv_rkpt,phase) + cell,lapw_nk,lapw_rkpt,phase) use m_juDFT USE m_constants @@ -175,14 +175,12 @@ TYPE(t_cell),INTENT(IN) :: cell TYPE(t_kpts),INTENT(IN) :: kpts TYPE(t_atoms),INTENT(IN) :: atoms - + TYPE(t_lapw),INTENT(IN) :: lapw_nk,lapw_rkpt ! - scalars - INTEGER,INTENT(IN) :: nk,jsp ,nbands INTEGER,INTENT(IN) :: iop LOGICAL :: phase ! - arrays - - INTEGER,INTENT(IN) :: nv_nk (dimension%jspd),gpt_nk (3,nv_nk (jsp)) - INTEGER,INTENT(IN) :: nv_rkpt(dimension%jspd),gpt_rkpt(3,nv_rkpt(jsp)) COMPLEX,INTENT(IN) :: cmt(dimension%neigd,hybrid%maxlmindx,atoms%nat) LOGICAL,INTENT(IN) :: l_real REAL,INTENT(IN) :: z_r(dimension%nbasfcn,dimension%neigd) @@ -281,18 +279,18 @@ ! PW coefficients zhlp = 0 - DO igpt = 1,nv_rkpt(jsp) - g = matmul( invrrot,gpt_rkpt(:,igpt)+g1 ) + DO igpt = 1,lapw_rkpt%nv(jsp) + g = matmul( invrrot,(/lapw_rkpt%k1(igpt,jsp),lapw_rkpt%k2(igpt,jsp),lapw_rkpt%k3(igpt,jsp)/)+g1 ) !determine number of g igpt1 = 0 - DO i=1,nv_nk(jsp) - IF ( maxval( abs( g - gpt_nk(:,i) ) ) .le. 1E-06 ) THEN + DO i=1,lapw_nk%nv(jsp) + IF ( maxval( abs( g - (/lapw_nk%k1(i,jsp),lapw_nk%k2(i,jsp),lapw_nk%k3(i,jsp)/) ) ) .le. 1E-06 ) THEN igpt1 = i EXIT END IF END DO IF ( igpt1 .eq. 0 ) CYCLE - cdum = exp(tpiimg*dotprod(rkpt+gpt_rkpt(:,igpt),trans ) ) + cdum = exp(tpiimg*dotprod(rkpt+(/lapw_rkpt%k1(igpt,jsp),lapw_rkpt%k2(igpt,jsp),lapw_rkpt%k3(igpt,jsp)/),trans ) ) if (l_real) THEN zhlp(igpt,:nbands) = cdum*z_r(igpt1,:nbands) else -- 2.26.2