Commit d3c897e3 authored by Daniel Wortmann's avatar Daniel Wortmann

Changed parallelization of cdnval

parent 3d075764
......@@ -94,8 +94,8 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
LOGICAL :: l_orbcomprot, l_real, l_dosNdir, l_corespec
! Local Arrays
REAL :: we(MAXVAL(cdnvalJob%noccbd(:)))
REAL :: eig(MAXVAL(cdnvalJob%noccbd(:)))
REAL,ALLOCATABLE :: we(:),eig(:)
INTEGER,ALLOCATABLE :: ev_list(:)
REAL, ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) ! radial functions
TYPE (t_lapw) :: lapw
......@@ -179,19 +179,20 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
CALL lapw%init(input,noco, kpts,atoms,sym,ikpt,cell,.false., mpi)
skip_t = skip_tt
noccbd = cdnvalJob%noccbd(ikpt)
we(1:noccbd) = cdnvalJob%weights(1:noccbd,ikpt)
eig(1:noccbd) = results%eig(cdnvalJob%ev_list(1:noccbd),ikpt,jsp)
ev_list=cdnvaljob%compact_ev_list(ikpt_i,.not.banddos%dos)
noccbd = SIZE(ev_list)
we = cdnvalJob%weights(ev_list,ikpt)
eig = results%eig(ev_list,ikpt,jsp)
IF (cdnvalJob%l_evp) THEN
IF (minval(cdnvalJob%ev_list(1:noccbd)) > skip_tt) skip_t = 0
IF (maxval(cdnvalJob%ev_list(1:noccbd)) <= skip_tt) skip_t = noccbd
IF ((minval(cdnvalJob%ev_list(1:noccbd)) <= skip_tt).AND.(maxval(cdnvalJob%ev_list(1:noccbd)) > skip_tt)) skip_t = mod(skip_tt,noccbd)
IF (minval(ev_list) > skip_tt) skip_t = 0
IF (maxval(ev_list) <= skip_tt) skip_t = noccbd
IF ((minval(ev_list) <= skip_tt).AND.(maxval(ev_list) > skip_tt)) skip_t = mod(skip_tt,noccbd)
END IF
nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
CALL zMat%init(l_real,nbasfcn,noccbd)
CALL read_eig(eig_id,ikpt,jsp,list=cdnvaljob%ev_list(1:noccbd),neig=nbands,zmat=zMat)
CALL read_eig(eig_id,ikpt,jsp,list=ev_list,neig=nbands,zmat=zMat)
#ifdef CPP_MPI
CALL MPI_BARRIER(mpi%mpi_comm,iErr) ! Synchronizes the RMA operations
#endif
......
......@@ -133,7 +133,7 @@ CONTAINS
DO jsp = 1,MERGE(1,input%jspins,noco%l_noco)
k_loop:DO nk_i = 1,size(mpi%k_list)
nk=mpi%k_list(i)
nk=mpi%k_list(nk_i)
! Set up lapw list
CALL lapw%init(input,noco, kpts,atoms,sym,nk,cell,l_zref, mpi)
call timestart("Setup of H&S matrices")
......
......@@ -97,7 +97,7 @@ CONTAINS
WRITE(6,*) "Non-SOC ev for nk,jsp:",nk,jsp
WRITE(6,"(6(f10.6,1x))") eig(:ne,jsp)
ENDIF
CALL read_eig(eig_id,nk,jsp,n_start=1,n_end=ne,zmat=zmat(jsp))
CALL read_eig(eig_id,nk,jsp,list=[(i,i=1,ne)],zmat=zmat(jsp))
! write(*,*) 'process',irank,' reads ',nk
......
......@@ -73,7 +73,7 @@ c n_bands(i) the number of ev's processed on n_rank=0...i-1
TYPE(t_mat), INTENT (INOUT) :: zmat !z(nbasfcn,noccbd) !can be real/complex
! Local variables ...
INTEGER :: isp
INTEGER :: isp,i
#ifdef CPP_MPI
INCLUDE 'mpif.h'
INTEGER mpiierr
......@@ -83,7 +83,7 @@ c n_bands(i) the number of ev's processed on n_rank=0...i-1
CALL timestart("cdn_read")
CALL read_eig(eig_id,ikpt,isp,n_start=n_start,n_end=n_end,
CALL read_eig(eig_id,ikpt,isp,list=[(i,i=n_start,n_end)],
< neig=nbands,eig=eig,zmat=zmat)
CALL timestop("cdn_read")
......
......@@ -115,12 +115,12 @@ CONTAINS
d%fname="eig"
CALL eig66_remove_data(id)
END SUBROUTINE close_eig
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zmat)
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,list,zmat)
IMPLICIT NONE
INTEGER, INTENT(IN) :: id,nk,jspin
INTEGER, INTENT(OUT),OPTIONAL :: neig
REAL, INTENT(OUT),OPTIONAL :: eig(:),w_iks(:)
INTEGER, INTENT(IN),OPTIONAL :: n_start,n_end
INTEGER, INTENT(IN),OPTIONAL :: list(:)
TYPE(t_mat),OPTIONAL :: zmat
!Local variables
......@@ -134,8 +134,8 @@ CONTAINS
CALL priv_find_data(id,d)
! check if io is performed correctly
IF (PRESENT(n_start)) THEN
IF (n_start/=1) &
IF (PRESENT(list)) THEN
IF (list(1)/=1) &
CALL juDFT_error("In direct access mode only all eigenstates can be read")
ENDIF
......
......@@ -192,15 +192,21 @@ CONTAINS
REAL, INTENT(OUT) :: z(:,:)
INTEGER :: nmat
INTEGER i,j,neig_l
INTEGER i
neig_l = n_end - n_start + 1
nmat=SIZE(z,1)
!read eigenvectors
CALL io_read_real2(d%evsetid,(/1,1,n_start,nk,jspin/),&
& (/1,nmat,neig_l,1,1/),&
& z(:nmat,:neig_l) )
IF (.NOT.PRESENT(list)) THEN
! read all eigenvectors
CALL io_read_real2(d%evsetid,(/1,1,1,nk,jspin/),&
(/1,nmat,SIZE(z,2),1,1/),z(:nmat,:) )
ELSE
DO i=1,SIZE(list)
CALL io_read_real2(d%evsetid,(/1,1,list(i),nk,jspin/),&
& (/1,nmat,1,1,1/),z(:nmat,i))
ENDDO
END IF
END SUBROUTINE priv_r_vec
......@@ -283,47 +289,49 @@ CONTAINS
!----------------------------------------------------------------------
SUBROUTINE priv_r_vecc(&
& d,nk,jspin,n_start,n_end,z)
& d,nk,jspin,list,z)
USE m_hdf_tools
IMPLICIT NONE
TYPE(t_data_HDF),INTENT(IN)::d
INTEGER, INTENT(IN) :: nk,jspin
INTEGER, INTENT(IN) :: n_start,n_end
INTEGER,OPTIONAL, INTENT(IN) :: list(:)
COMPLEX, INTENT(OUT) :: z(:,:)
REAL, ALLOCATABLE :: z1(:,:,:)
INTEGER i,j,neig_l
INTEGER i,j
INTEGER :: nmat
neig_l = n_end - n_start + 1
nmat=SIZE(z,1)
! read eigenvectors
ALLOCATE (z1(2,nmat,neig_l))
CALL io_read_real3(d%evsetid,(/1,1,n_start,nk,jspin/),&
& (/2,nmat,neig_l,1,1/),z1)
DO i=1,neig_l
DO j=1,nmat
z(j,i) = CMPLX( z1(1,j,i) ,z1(2,j,i) )
IF (.NOT.PRESENT(list)) THEN
! read all eigenvectors
ALLOCATE (z1(2,nmat,SIZE(z,2)))
CALL io_read_real3(d%evsetid,(/1,1,1,nk,jspin/),&
& (/2,nmat,SIZE(z,2),1,1/),z1)
DO i=1,SIZE(z,2)
z(:,i) = CMPLX( z1(1,:,i) ,z1(2,:,i) )
ENDDO
ENDDO
DEALLOCATE (z1)
ELSE
ALLOCATE (z1(2,nmat,1))
DO i=1,SIZE(list)
CALL io_read_real3(d%evsetid,(/1,1,list(i),nk,jspin/),&
& (/2,nmat,1,1,1/),z1
z(:,i) = CMPLX( z1(1,:,i) ,z1(2,:,i) )
ENDDO
END IF
END SUBROUTINE priv_r_vecc
!-----------------------------------------------------------------------
#endif
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zMat)
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,list,zMat)
IMPLICIT NONE
INTEGER, INTENT(IN) :: id,nk,jspin
INTEGER, INTENT(OUT),OPTIONAL :: neig
REAL, INTENT(OUT),OPTIONAL :: eig(:),w_iks(:)
INTEGER, INTENT(IN),OPTIONAL :: n_start,n_end
INTEGER, INTENT(IN),OPTIONAL :: list(:)
TYPE(t_mat),OPTIONAL :: zmat
#ifdef CPP_HDF
......@@ -348,14 +356,11 @@ CONTAINS
ENDIF
ENDIF
IF (PRESENT(n_start)) THEN
IF (.NOT.PRESENT(n_end)) CALL juDFT_error("BUG3 in read_eig")
IF (PRESENT(zMat)) THEN
IF (zmat%l_real) THEN
CALL priv_r_vec(d,nk,jspin,n_start,n_end,zmat%data_r)
ELSE
CALL priv_r_vecc(d,nk,jspin,n_start,n_end,zmat%data_c)
ENDIF
IF (PRESENT(zMat)) THEN
IF (zmat%l_real) THEN
CALL priv_r_vec(d,nk,jspin,list,zmat%data_r)
ELSE
CALL priv_r_vecc(d,nk,jspin,list,zmat%data_c)
ENDIF
ENDIF
#endif
......
......@@ -109,7 +109,7 @@ CONTAINS
END SUBROUTINE close_eig
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zmat)
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,list,zmat)
USE m_eig66_hdf,ONLY:read_eig_hdf=>read_eig
USE m_eig66_DA ,ONLY:read_eig_DA=>read_eig
USE m_eig66_mem,ONLY:read_eig_mem=>read_eig
......@@ -118,19 +118,19 @@ CONTAINS
INTEGER, INTENT(IN) :: id,nk,jspin
INTEGER, INTENT(OUT),OPTIONAL :: neig
REAL, INTENT(OUT),OPTIONAL :: eig(:),w_iks(:)
INTEGER, INTENT(IN),OPTIONAL :: n_start,n_end
INTEGER, INTENT(IN),OPTIONAL :: list(:)
TYPE(t_mat),INTENT(INOUT),OPTIONAL :: zmat
INTEGER::n
CALL timestart("IO (read)")
SELECT CASE (eig66_data_mode(id))
CASE (DA_mode)
CALL read_eig_DA(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zmat)
CALL read_eig_DA(id,nk,jspin,neig,eig,w_iks,list,zmat)
CASE (hdf_mode)
CALL read_eig_hdf(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zmat)
CALL read_eig_hdf(id,nk,jspin,neig,eig,w_iks,list,zmat)
CASE (mem_mode)
CALL read_eig_mem(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zmat)
CALL read_eig_mem(id,nk,jspin,neig,eig,w_iks,list,zmat)
CASE (mpi_mode)
CALL read_eig_mpi(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zmat)
CALL read_eig_mpi(id,nk,jspin,neig,eig,w_iks,list,zmat)
CASE (-1)
CALL juDFT_error("Could not read eig-file before opening", calledby = "eig66_io")
END SELECT
......
......@@ -131,15 +131,16 @@ CONTAINS
END SUBROUTINE priv_writetofile
END SUBROUTINE close_eig
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zmat)
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,list,zmat)
IMPLICIT NONE
INTEGER, INTENT(IN) :: id,nk,jspin
INTEGER, INTENT(OUT),OPTIONAL :: neig
REAL, INTENT(OUT),OPTIONAL :: eig(:),w_iks(:)
INTEGER, INTENT(IN),OPTIONAL :: n_start,n_end
INTEGER, INTENT(IN),OPTIONAL :: list(:)
TYPE(t_mat),OPTIONAL :: zmat
INTEGER::nrec, arrayStart
INTEGER::nrec, arrayStart,arrayStop,i
INTEGER,ALLOCATABLE :: ind(:)
TYPE(t_data_mem),POINTER:: d
CALL priv_find_data(id,d)
......@@ -161,23 +162,37 @@ CONTAINS
!data from d%eig_vec
arrayStart = 1
IF(PRESENT(n_start)) THEN
arrayStart = (n_start-1)*zMat%matsize1+1
END IF
IF (PRESENT(zmat)) THEN
IF(PRESENT(list)) THEN
ind=list
ELSE
ALLOCATE(ind(SIZE(zmat%data_r,2)))
ind=[(i,i=1,SIZE(ind))]
END IF
IF (zmat%l_real) THEN
IF (.NOT.ALLOCATED(d%eig_vecr)) THEN
IF (.NOT.ALLOCATED(d%eig_vecc)) CALL juDFT_error("BUG: can not read real/complex vectors from memory")
zmat%data_r=REAL(RESHAPE(d%eig_vecc(arrayStart:arrayStart+SIZE(zmat%data_r)-1,nrec),SHAPE(zmat%data_r)))
DO i=1,SIZE(ind)
arrayStart=(ind(i)-1)*zMat%matsize1+1
arrayStop=ind(i)*zMat%matsize1+1
zmat%data_r(:,i)=REAL(d%eig_vecc(arrayStart:arrayStop,nrec))
ENDDO
ELSE
zmat%data_r=RESHAPE(d%eig_vecr(arrayStart:arrayStart+SIZE(zmat%data_r)-1,nrec),SHAPE(zmat%data_r))
DO i=1,SIZE(ind)
arrayStart=(ind(i)-1)*zMat%matsize1+1
arrayStop=ind(i)*zMat%matsize1+1
zmat%data_r(:,i)=d%eig_vecr(arrayStart:arrayStop,nrec)
ENDDO
ENDIF
ELSE !TYPE is (COMPLEX)
IF (.NOT.ALLOCATED(d%eig_vecc)) CALL juDFT_error("BUG: can not read complex vectors from memory", calledby = "eig66_mem")
zmat%data_c=RESHAPE(d%eig_vecc(arrayStart:arrayStart+SIZE(zmat%data_c)-1,nrec),SHAPE(zmat%data_c))
DO i=1,SIZE(ind)
arrayStart=(ind(i)-1)*zMat%matsize1+1
arrayStop=ind(i)*zMat%matsize1+1
zmat%data_c(:,i)=d%eig_vecc(arrayStart:arrayStop,nrec)
END DO
END IF
ENDIF
END SUBROUTINE read_eig
......
......@@ -156,12 +156,12 @@ CONTAINS
IF (PRESENT(filename)) CALL judft_error("Storing of data not implemented for MPI case",calledby="eig66_mpi.F")
END SUBROUTINE close_eig
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zmat)
SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,list,zmat)
IMPLICIT NONE
INTEGER, INTENT(IN) :: id,nk,jspin
INTEGER, INTENT(OUT),OPTIONAL :: neig
REAL, INTENT(OUT),OPTIONAL :: eig(:),w_iks(:)
INTEGER, INTENT(IN),OPTIONAL :: n_start,n_end
INTEGER, INTENT(IN),OPTIONAL :: list(:)
TYPE(t_mat),OPTIONAL :: zmat
#ifdef CPP_MPI
......@@ -183,21 +183,18 @@ CONTAINS
ENDIF
IF (PRESENT(eig).or.PRESENT(w_iks)) THEN
ALLOCATE(tmp_real(d%size_eig))
ALLOCATE(tmp_real(MIN(SIZE(eig),d%size_eig)))
IF (PRESENT(eig)) THEN
CALL MPI_WIN_LOCK(MPI_LOCK_SHARED,pe,0,d%eig_handle,e)
CALL MPI_GET(tmp_real,d%size_eig,MPI_DOUBLE_PRECISION,pe,slot,d%size_eig,MPI_DOUBLE_PRECISION,d%eig_handle,e)
CALL MPI_GET(tmp_real,SIZE(tmp_real),MPI_DOUBLE_PRECISION,pe,slot,d%size_eig,MPI_DOUBLE_PRECISION,d%eig_handle,e)
CALL MPI_WIN_UNLOCK(pe,d%eig_handle,e)
n1=1;n3=1;n2=SIZE(eig)
IF (PRESENT(n_start)) n1=n_start
IF (PRESENT(n_end)) n2=n_end
eig(:n2-n1+1)=tmp_real(n1:n2)
eig(:size(tmp_real))=tmp_real
END IF
IF (PRESENT(w_iks)) THEN
CALL MPI_WIN_LOCK(MPI_LOCK_SHARED,pe,0,d%w_iks_handle,e)
CALL MPI_GET(tmp_real,d%size_eig,MPI_DOUBLE_PRECISION,pe,slot,d%size_eig,MPI_DOUBLE_PRECISION,d%w_iks_handle,e)
CALL MPI_GET(tmp_real,size(tmp_real),MPI_DOUBLE_PRECISION,pe,slot,d%size_eig,MPI_DOUBLE_PRECISION,d%w_iks_handle,e)
CALL MPI_WIN_UNLOCK(pe,d%w_iks_handle,e)
w_iks=tmp_real(:size(w_iks))
w_iks(:SIZE(tmp_real))=tmp_real
END IF
DEALLOCATE(tmp_real)
ENDIF
......@@ -207,11 +204,10 @@ CONTAINS
ALLOCATE(tmp_real(tmp_size))
ALLOCATE(tmp_cmplx(tmp_size))
DO n=1,zmat%matsize2
n1=n
IF (PRESENT(n_start)) n1=n_start+n-1
IF (PRESENT(n_end)) THEN
IF (n1>n_end) CYCLE
ENDIF
IF (PRESENT(list)) THEN
IF (n>SIZE(list)) CYCLE
n1=list(n)
END IF
slot=d%slot_ev(nk,jspin,n1)
pe=d%pe_ev(nk,jspin,n1)
......
......@@ -114,7 +114,8 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
jspmax = input%jspins
IF (noco%l_mperp) jspmax = 1
DO jspin = 1,jspmax
CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin,sliceplot,banddos)
CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin)
IF (sliceplot%slice) CALL cdnvalJob%select_slice(sliceplot,results,input,kpts,noco,jspin)
CALL cdnval(eig_id,mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,stars,vacuum,dimension,&
sphhar,sym,vTot,oneD,cdnvalJob,outDen,regCharges,dos,results,moments,coreSpecInput,mcd,slab,orbcomp)
END DO
......
......@@ -42,10 +42,7 @@ CONTAINS
mpi%n_rank=0
mpi%n_size=1
mpi%sub_comm=mpi%mpi_comm
IF (ALLOCATED(mpi%k_list)) DEALLOCATE(mpi%k_List,mpi%ev_list)
ALLOCATE(mpi%k_list(nkpt))
mpi%k_list=[(i,i=1,nkpt)]
ALLOCATE(mpi%ev_list(neigd))
mpi%ev_list=[(i,i=1,neigd)]
WRITE(*,*) "--------------------------------------------------------"
RETURN
......@@ -58,9 +55,10 @@ CONTAINS
IF (mpi%n_size>1.AND..NOT.parallel_solver_available()) &
CALL juDFT_error("MPI parallelization failed",hint="You have to either compile FLEUR with a parallel diagonalization library (ELPA,SCALAPACK...) or you have to run such that the No of kpoints can be distributed on the PEs")
#endif
!generate the MPI communicators
CALL priv_create_comm(nkpt,neigd,mpi)
mpi%k_list=[(i, i=INT(mpi%irank/mpi%n_size)+1,nkpt,mpi%isize/mpi%n_size )]
if (mpi%irank==0) WRITE(*,*) "--------------------------------------------------------"
END SUBROUTINE setupMPI
......@@ -111,9 +109,7 @@ CONTAINS
IF ((MOD(mpi%isize,n_members) == 0).AND.(MOD(nkpt,n_members) == 0) ) EXIT
n_members = n_members - 1
ENDDO
ALLOCATE(mpi%k_list(nkpt/n_members))
mpi%k_list=[(nk, nk=nkpt/n_members,nkpt,n_members )]
!mpi%n_groups = nkpt/n_members
mpi%n_size = mpi%isize/n_members
!mpi%n_stride = n_members
......@@ -196,7 +192,7 @@ CONTAINS
CALL MPI_COMM_RANK (mpi%SUB_COMM,mpi%n_rank,ierr)
ALLOCATE(mpi%ev_list(neigd/mpi%n_size+1))
mpi%ev_list=[(i,i=mpi%irank+1,neigd,mpi%n_size)]
mpi%ev_list=[(i,i=mpi%n_rank+1,neigd,mpi%n_size)]
#endif
END SUBROUTINE priv_create_comm
......
......@@ -80,7 +80,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
TYPE(t_mat) :: exMat, zMat, olap, trafo, invtrafo, tmpMat, exMatLAPW
TYPE(t_lapw) :: lapw
TYPE(t_hybdat) :: hybdat
INTEGER :: ikpt, iBand, jkpt, jBand, iAtom, i, na, itype, lh, j
INTEGER :: ikpt,ikpt_i,iband_i, iBand, jkpt, jBand, iAtom, i, na, itype, lh, j
INTEGER :: jspin, jspmax, jsp, isp, ispin, nbasfcn, nbands
INTEGER :: nsymop, nkpt_EIBZ, ikptf, iterHF, mnobd
INTEGER :: iState, iStep, numStates, maxHistoryLength, numRelevantStates
......@@ -279,35 +279,23 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
vTotSSDen = 0.0
! Calculate all single state densities
cdnvalJob%l_evp = .FALSE.
cdnvalJob%nkptExtended = kpts%nkpt
ALLOCATE(cdnvalJob%noccbd(kpts%nkpt))
ALLOCATE(cdnvalJob%nStart(kpts%nkpt))
ALLOCATE(cdnvalJob%nEnd(kpts%nkpt))
ALLOCATE(cdnvalJob%weights(1,kpts%nkpt))
CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin)
numStates = 0
DO jspin = 1, input%jspins
jsp = MERGE(1,jspin,noco%l_noco)
DO ikpt = 1, kpts%nkpt
DO iBand = 1, highestState(ikpt,jsp)
DO ikpt_i = 1, SIZE(mpi%k_list)
ikpt= mpi%k_list(ikpt_i)
DO iBand_i = 1,size(cdnvalJOB%ev_list)
iband=mpi%ev_list(iband_i)
IF (iband>highestState(ikpt,jsp)) CYCLE
numStates = numStates + 1
! Construct cdnvalJob object for this state
! (Reasonable parallelization is not yet done - should be placed over the loops enclosing this section)
cdnvalJob%ikptStart = ikpt
cdnvalJob%ikptIncrement = kpts%nkpt
IF(mpi%irank.EQ.0) THEN
cdnvalJob%noccbd = 1
cdnvalJob%nStart = iBand
cdnvalJob%nEnd = iBand
cdnvalJob%weights = 0.0
cdnvalJob%weights(1,ikpt) = spinDegenFac ! Note the doubling of the weight for non-spinpolarized calculations.
ELSE
cdnvalJob%noccbd = 0
cdnvalJob%nStart = 1
cdnvalJob%nEnd = 0
cdnvalJob%weights = 0.0
END IF
cdnvalJob%k_list=[ikpt]
cdnvalJob%ev_list=[iband]
cdnvalJob%weights(iBand,ikpt) = spinDegenFac
! Call cdnval to construct density
WRITE(*,*) 'Note: some optional flags may have to be reset in rdmft before the cdnval call'
WRITE(*,*) 'This is not yet implemented!'
......@@ -394,7 +382,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
jspmax = input%jspins
IF (noco%l_mperp) jspmax = 1
DO jspin = 1,jspmax
CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin,banddos=banddos)
CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin)
CALL cdnval(eig_id,mpi,kpts,jsp,noco,input,banddos,cell,atoms,enpara,stars,vacuum,dimension,&
sphhar,sym,vTot,oneD,cdnvalJob,overallDen,regCharges,dos,results,moments)
END DO
......@@ -515,7 +503,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
CALL zMat%init(olap%l_real,nbasfcn,dimension%neigd)
CALL read_eig(eig_id,ikpt,jspin,n_start=1,n_end=hybrid%nbands(ikpt),neig=nbands,zmat=zMat)
CALL read_eig(eig_id,ikpt,jspin,list=[(i,i=1,hybrid%nbands(ikpt))],neig=nbands,zmat=zMat)
! CALL read_z(zMat,kpts%nkpt*(jspin-1)+ikpt)
zMat%matsize2 = hybrid%nbands(ikpt) ! reduce "visible matsize" for the following computations
......
......@@ -122,11 +122,14 @@ PRIVATE
LOGICAL :: l_evp
INTEGER, ALLOCATABLE :: k_list(:)
INTEGER, ALLOCATABLE :: ev_list(:)
INTEGER, ALLOCATABLE :: noccbd(:)
REAL, ALLOCATABLE :: weights(:,:) ! weights(band_idx, kpt_idx)
INTEGER, ALLOCATABLE :: noccbd(:) ! Attention, these are for all k-points and all states
REAL, ALLOCATABLE :: weights(:,:) !
CONTAINS
PROCEDURE,PASS :: init => cdnvalJob_init
PROCEDURE :: select_slice
PROCEDURE :: compact_ev_list
END TYPE t_cdnvalJob
TYPE t_gVacMap
......@@ -441,7 +444,7 @@ SUBROUTINE orbcomp_init(thisOrbcomp,input,banddos,dimension,atoms,kpts)
END SUBROUTINE orbcomp_init
SUBROUTINE cdnvalJob_init(thisCdnvalJob,mpi,input,kpts,noco,results,jspin,sliceplot,banddos)
SUBROUTINE cdnvalJob_init(thisCdnvalJob,mpi,input,kpts,noco,results,jspin)
USE m_types_mpi
USE m_types_setup
......@@ -450,83 +453,103 @@ SUBROUTINE cdnvalJob_init(thisCdnvalJob,mpi,input,kpts,noco,results,jspin,slicep
IMPLICIT NONE
CLASS(t_cdnvalJob), INTENT(INOUT) :: thisCdnvalJob
CLASS(t_cdnvalJob), INTENT(OUT) :: thisCdnvalJob
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_input), INTENT(IN) :: input
TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_results), INTENT(IN) :: results
TYPE(t_sliceplot), OPTIONAL, INTENT(IN) :: sliceplot
TYPE(t_banddos), OPTIONAL, INTENT(IN) :: banddos
INTEGER, INTENT(IN) :: jspin
INTEGER :: jsp, iBand, ikpt, nslibd, noccbd_l, noccbd, ikpt_i
thisCdnvalJob%l_evp = mpi%n_size>1
IF (ALLOCATED(thisCdnvalJob%k_list))DEALLOCATE (thisCdnvalJob%k_list)
IF (ALLOCATED(thisCdnvalJob%ev_list)) DEALLOCATE (thisCdnvalJob%ev_list)
jsp = MERGE(1,jspin,noco%l_noco)
thisCdnvalJob%l_evp=mpi%n_size>1
thisCdnvalJob%k_list=mpi%k_list !includes allocate
thisCdnvalJob%ev_list=mpi%ev_list
IF (ALLOCATED(thisCdnvalJob%noccbd)) DEALLOCATE (thisCdnvalJob%noccbd)
IF (ALLOCATED(thisCdnvalJob%weights)) DEALLOCATE (thisCdnvalJob%weights)
ALLOCATE(thisCdnvalJob%noccbd(size(thiscdnvaljob%k_list)))
ALLOCATE(thisCdnvalJob%weights(size(thiscdnvaljob%ev_list),size(thiscdnvalJob%k_list)))
thisCdnvalJob%weights(:,:) = results%w_iks(thiscdnvaljob%ev_list,:,jsp)*2.0/input%jspins
thisCdnvalJob%weights = results%w_iks(:,:,jsp)*2.0/input%jspins
ALLOCATE(thisCdnvalJob%noccbd(kpts%nkpt))
thisCdnvalJob%noccbd = 0
jsp = MERGE(1,jspin,noco%l_noco)
! determine bands to be used for each k point, MPI process
DO ikpt_i = 1,size(thisCdnvalJob%k_list)
DO ikpt_i = 1,SIZE(thisCdnvalJob%k_list)
ikpt=thisCdnvalJob%k_list(ikpt_i)
IF(PRESENT(sliceplot)) THEN
IF (sliceplot%slice.AND.input%pallst) thisCdnvalJob%weights(:,ikpt) = kpts%wtkpt(ikpt)
END IF
!Max number of bands
thisCdnvalJob%noccbd(ikpt)= count(thiscdnvaljob%ev_list<=results%neig(ikpt,jsp))
thisCdnvalJob%noccbd(ikpt)= COUNT(thiscdnvaljob%ev_list<=results%neig(ikpt,jsp))
ENDDO
END SUBROUTINE cdnvalJob_init
SUBROUTINE select_slice(thiscdnvalJob,sliceplot,results,input,kpts,noco,jspin)
USE m_types_setup
USE m_types_misc
IMPLICIT NONE
CLASS(t_cdnvalJob),INTENT(INOUT) :: thisCdnvalJob
TYPE(t_sliceplot), INTENT(IN) :: sliceplot
TYPE(t_results), INTENT(IN) :: results
TYPE(t_input),INTENT(IN) :: input
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_noco),INTENT(IN) :: noco
INTEGER,INTENT(IN) :: jspin
INTEGER :: iband,iband_i,ikpt,ikpt_i,jsp
jsp = MERGE(1,jspin,noco%l_noco)
DO ikpt_i=1,SIZE(thiscdnvalJob%k_list)
ikpt=thiscdnvalJob%k_list(ikpt_i)
!---> if slice, only certain bands are taken into account
IF(PRESENT(sliceplot)) THEN
IF (sliceplot%slice.AND.thisCdnvalJob%noccbd(ikpt).GT.0) THEN
IF (sliceplot%kk.EQ.0) THEN
do iband=1,thisCdnvalJob%noccbd(ikpt)
if (results%eig(thiscdnvaljob%ev_list(iBand),ikpt,jsp).LT.sliceplot%e1s) thisCdnvalJob%weights(iband,ikpt)=0.0
if (results%eig(thiscdnvaljob%ev_list(iBand),ikpt,jsp).GT.sliceplot%e2s) thisCdnvalJob%weights(iband,ikpt)=0.0
end do
ELSE IF (sliceplot%kk.EQ.ikpt) THEN
IF ((sliceplot%e1s.EQ.0.0) .AND. (sliceplot%e2s.EQ.0.0)) THEN
do iband=1,thisCdnvalJob%noccbd(ikpt)
if (thiscdnvaljob%ev_list(iBand).ne.sliceplot%nnne) thisCdnvalJob%weights(iband,ikpt)=0.0
enddo
ELSE
do iband=1,thisCdnvalJob%noccbd(ikpt)
if (results%eig(thiscdnvaljob%ev_list(iBand),ikpt,jsp).LT.sliceplot%e1s) thisCdnvalJob%weights(iband,ikpt)=0.0
if (results%eig(thiscdnvaljob%ev_list(iBand),ikpt,jsp).GT.sliceplot%e2s) thisCdnvalJob%weights(iband,ikpt)=0.0
end do
END IF
else
thisCdnvalJob%weights(:,ikpt)=0.0
END IF
END IF ! sliceplot%slice
END IF
!remove unoccupied states
iband=thisCdnvalJob%noccbd(ikpt)
IF(PRESENT(banddos)) THEN
IF (.NOT.banddos%dos)THEN
DO WHILE (thisCdnvalJob%weights(iBand,ikpt).LT.1.e-8)
iBand = iBand - 1
IF (sliceplot%slice.AND.input%pallst) thisCdnvalJob%weights(:,ikpt) = kpts%wtkpt(ikpt)
IF (sliceplot%slice.AND.thisCdnvalJob%noccbd(ikpt).GT.0) THEN
IF (sliceplot%kk.EQ.0) THEN
DO iband_i=1,thisCdnvalJob%noccbd(ikpt)
iband=thiscdnvaljob%ev_list(iband_i)
IF (results%eig(iBand,ikpt,jsp).LT.sliceplot%e1s) thisCdnvalJob%weights(iband,ikpt)=0.0
IF (results%eig(iBand,ikpt,jsp).GT.sliceplot%e2s) thisCdnvalJob%weights(iband,ikpt)=0.0
END DO
ENDIF
ENDIF
thisCdnvalJob%noccbd(ikpt) = iBand
ELSE IF (sliceplot%kk.EQ.ikpt) THEN
IF ((sliceplot%e1s.EQ.0.0) .AND. (sliceplot%e2s.EQ.0.0)) THEN
DO iband_i=1,thisCdnvalJob%noccbd(ikpt)
iband=thiscdnvaljob%ev_list(iband_i)
IF (iBand.NE.sliceplot%nnne) thisCdnvalJob%weights(iband,ikpt)=0.0
ENDDO
ELSE
DO iband_i=1,thisCdnvalJob%noccbd(ikpt)
iband=thiscdnvaljob%ev_list(iband_i)
IF (results%eig(iBand,ikpt,jsp).LT.sliceplot%e1s) thisCdnvalJob%weights(iband,ikpt)=0.0
IF (results%eig(iBand,ikpt,jsp).GT.sliceplot%e2s) thisCdnvalJob%weights(iband,ikpt)=0.0
END DO
END IF
ELSE
thisCdnvalJob%weights(:,ikpt)=0.0
END IF
END IF ! sliceplot%slice
END DO
END SUBROUTINE select_slice
FUNCTION compact_ev_list(thiscdnvaljob,ikpt,l_empty)
IMPLICIT NONE
CLASS(t_cdnvalJob),INTENT(IN) :: thisCdnvalJob
INTEGER,INTENT(IN) :: ikpt
LOGICAL,INTENT(IN) :: l_empty
INTEGER,ALLOCATABLE :: compact_ev_list(:)