eigen.F90 11.1 KB
Newer Older
1 2 3 4 5
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------
6 7 8
!<@brief
!<The eigen routine sets up and solves the generalized eigenvalue problem
!More description at end of file
9
MODULE m_eigen
10
  USE m_juDFT
11
  IMPLICIT NONE
12
CONTAINS
13 14 15 16 17 18 19 20
  !>The eigenvalue problem is constructed and solved in this routine. The following steps are performed:
  !> 1. Preparation: generate energy parameters, open eig-file
  !> 2. CALL to mt_setup() : this constructs the local Hamiltonian (i.e. the Hamiltonian in the \f$ u,\dot u, u_{lo} \f$ basis) LDA+U is also added here
  !> 3. within the (collinear)spin and k-point loop: CALL to eigen_hssetup() to generate the matrices, CALL to eigen_diag() to perform diagonalization 
  !> 4. writing (saving) of eigenvectors
  !>
  !> The matrices generated and diagonalized here are of type m_mat as defined in m_types_mat. 
  !>@author D. Wortmann
21
  SUBROUTINE eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,sym,kpts,DIMENSION,vacuum,input,&
22
                   cell,enpara,banddos,noco,oneD,hybrid,iter,eig_id,results,inden,v,vx)
23

24
#include"cpp_double.h"
25
    USE m_constants, ONLY : pi_const,sfp_const
26 27
    USE m_types
    USE m_apws
28 29
    USE m_eigen_hssetup
    USE m_pot_io
30
    USE m_eigen_diag
31
    USE m_add_vnonlocal
32
    USE m_subvxc
33
    !USE m_hsefunctional
34
    USE m_mt_setup
35
    USE m_util
36 37
    USE m_io_hybrid
    !USE m_icorrkeys
38
    USE m_eig66_io, ONLY : open_eig, write_eig, close_eig,read_eig
39
    USE m_xmlOutput
40
#ifdef CPP_MPI
41
    USE m_mpi_bc_potden
42
#endif
43
    USE m_symmetrize_matrix
44 45
    USE m_unfold_band_kpts !used for unfolding bands
    USE m_types_mpimat
46

47
    
48 49
    IMPLICIT NONE
    TYPE(t_results),INTENT(INOUT):: results
50
    CLASS(t_xcpot),INTENT(IN)    :: xcpot
51
    TYPE(t_mpi),INTENT(IN)       :: mpi
52
    TYPE(t_dimension),INTENT(IN) :: DIMENSION
53
    TYPE(t_oneD),INTENT(IN)      :: oneD
54
    TYPE(t_hybrid),INTENT(INOUT) :: hybrid
55
    TYPE(t_enpara),INTENT(INOUT) :: enpara
56 57 58 59
    TYPE(t_obsolete),INTENT(IN)  :: obsolete
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_vacuum),INTENT(IN)    :: vacuum
    TYPE(t_noco),INTENT(IN)      :: noco
60
    TYPE(t_banddos),INTENT(IN)   :: banddos
Daniel Wortmann's avatar
Daniel Wortmann committed
61
    TYPE(t_sym),INTENT(IN)       :: sym  
62 63
    TYPE(t_stars),INTENT(IN)     :: stars
    TYPE(t_cell),INTENT(IN)      :: cell
64
    TYPE(t_kpts),INTENT(INOUT)   :: kpts
65
    TYPE(t_sphhar),INTENT(IN)    :: sphhar
66 67 68 69
    TYPE(t_atoms),INTENT(IN)     :: atoms
    TYPE(t_potden),INTENT(IN)    :: inden,vx
    TYPE(t_potden),INTENT(INOUT) :: v    !u_setup will modify the potential matrix
   
70 71 72
#ifdef CPP_MPI
    INCLUDE 'mpif.h'
#endif
73

74 75
!    EXTERNAL MPI_BCAST    !only used by band_unfolding to broadcast the gvec

76
    ! Scalar Arguments
77
    INTEGER,INTENT(IN)    :: iter
78
    INTEGER,INTENT(IN)    :: eig_id
79 80

    ! Local Scalars
81
    INTEGER jsp,nk,nred,ne_all,ne_found
82
    INTEGER ne,lh0
83
    INTEGER isp,i,j,err
Daniel Wortmann's avatar
Daniel Wortmann committed
84
    LOGICAL l_wu,l_file,l_real,l_zref
85 86 87 88

    ! Local Arrays
    INTEGER              :: ierr(3)
    INTEGER              :: neigBuffer(kpts%nkpt,input%jspins)
89 90 91

    COMPLEX              :: unfoldingBuffer(SIZE(results%unfolding_weights,1),kpts%nkpt,input%jspins) ! needed for unfolding bandstructure mpi case

92
    INTEGER, PARAMETER   :: lmaxb = 3
93
    REAL,    ALLOCATABLE :: bkpt(:)
94
    REAL,    ALLOCATABLE :: eig(:)
95
    COMPLEX, ALLOCATABLE :: vs_mmp(:,:,:,:)
96

97 98
	INTEGER                   :: jsp_m, i_kpt_m, i_m

99 100 101 102 103
    TYPE(t_tlmplm)            :: td
    TYPE(t_usdus)             :: ud
    TYPE(t_lapw)              :: lapw
    CLASS(t_mat), ALLOCATABLE :: zMat
    CLASS(t_mat), ALLOCATABLE :: hmat,smat
104
    CLASS(t_mat), ALLOCATABLE :: smat_unfold !used for unfolding bandstructure
105 106

    ! Variables for HF or hybrid functional calculation
107 108
    INTEGER                   :: comm(kpts%nkpt),irank2(kpts%nkpt),isize2(kpts%nkpt), dealloc_stat
    character(len=300)        :: errmsg
109
    
Daniel Wortmann's avatar
Daniel Wortmann committed
110
    call ud%init(atoms,DIMENSION%jspd)
111
    ALLOCATE (eig(DIMENSION%neigd),bkpt(3))
Daniel Wortmann's avatar
Daniel Wortmann committed
112
    
113
    l_real=sym%invs.AND..NOT.noco%l_noco
114

115
    ! check if z-reflection trick can be used
116 117 118 119
    l_zref=(sym%zrfs.AND.(SUM(ABS(kpts%bk(3,:kpts%nkpt))).LT.1e-9).AND..NOT.noco%l_noco) 
    IF (mpi%n_size > 1) l_zref = .FALSE.
 
#ifdef CPP_MPI
120
    CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,v)
121
#endif
122

123
    !IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber      '/),(/iter,v%iter/),&
124
    !                                                RESHAPE((/19,13,5,5/),(/2,2/)))
125
    
126 127 128
    ! Set up and solve the eigenvalue problem
    !   loop over spins
    !     set up k-point independent t(l'm',lm) matrices
129
    CALL mt_setup(atoms,sym,sphhar,input,noco,enpara,inden,v,mpi,results,DIMENSION,td,ud)
130

131
    neigBuffer = 0
132
    results%neig = 0
133
    results%eig = 1.0e300
134
    unfoldingBuffer = CMPLX(0.0,0.0)
135

136
    DO jsp = 1,MERGE(1,input%jspins,noco%l_noco)
Daniel Wortmann's avatar
Daniel Wortmann committed
137
       k_loop:DO nk = mpi%n_start,kpts%nkpt,mpi%n_stride
138

139
          ! Set up lapw list        
140 141
          CALL lapw%init(input,noco, kpts,atoms,sym,nk,cell,l_zref, mpi)
          call timestart("Setup of H&S matrices")
142
          CALL eigen_hssetup(jsp,mpi,DIMENSION,hybrid,enpara,input,vacuum,noco,sym,&
143
                             stars,cell,sphhar,atoms,ud,td,v,lapw,l_real,smat,hmat)
144
          CALL timestop("Setup of H&S matrices")
145

146
          IF(hybrid%l_hybrid) THEN
147 148

      DO i = 1, hmat%matsize1
149 150 151 152 153 154 155 156
	 DO j = 1, i
	    IF (hmat%l_real) THEN
	       IF ((i.LE.5).AND.(j.LE.5)) THEN
		  WRITE(1233,'(2i7,2f15.8)') i, j, hmat%data_r(i,j), hmat%data_r(j,i)
	       END IF
	    ELSE
	    ENDIF
	 END DO
157 158
      END DO

159
             ! Write overlap matrix smat to direct access file olap
Daniel Wortmann's avatar
Daniel Wortmann committed
160
             print *,"Wrong overlap matrix used, fix this later"
161
             CALL write_olap(smat,(jsp-1)*kpts%nkpt+nk) ! Note: At this moment this only works without MPI parallelization
162 163
             PRINT *,"TODO"
!             STOP "TODO"
164
             PRINT *,"BASIS:", lapw%nv(jsp), atoms%nlotot
165
             IF (hybrid%l_addhf) CALL add_Vnonlocal(nk,lapw,atoms,hybrid,dimension,kpts,jsp,results,xcpot,noco,hmat)
166 167

             IF(hybrid%l_subvxc) THEN
168 169
                CALL subvxc(lapw,kpts%bk(:,nk),DIMENSION,input,jsp,v%mt(:,0,:,:),atoms,ud,hybrid,enpara%el0,enpara%ello0,&
                            sym,cell,sphhar,stars,xcpot,mpi,oneD,hmat,vx)
170
             END IF
171
          END IF ! hybrid%l_hybrid
172

173 174
          l_wu=.FALSE.
          ne_all=DIMENSION%neigd
175 176 177 178 179 180
          if (allocated(zmat)) then
             deallocate(zmat, stat=dealloc_stat, errmsg=errmsg)
             if(dealloc_stat /= 0) call juDFT_error("deallocate failed for zmat",&
                                                hint=errmsg, calledby="eigen.F90")
          endif

181 182 183
          !Try to symmetrize matrix
          CALL symmetrize_matrix(mpi,noco,kpts,nk,hmat,smat)
          
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
          IF (banddos%unfoldband) THEN
		select type(smat)
		type is (t_mat)
			allocate(t_mat::smat_unfold)
			select type(smat_unfold)
		             type is (t_mat)
			     smat_unfold=smat
			end select
		type is (t_mpimat)
			allocate(t_mpimat::smat_unfold)
			select type(smat_unfold)
		             type is (t_mpimat)
			     smat_unfold=smat
			end select
		end select
          END IF

201
          CALL eigen_diag(mpi,hmat,smat,nk,jsp,iter,ne_all,eig,zMat)
202
          CALL smat%free()
203 204 205
          DEALLOCATE(hmat,smat, stat=dealloc_stat, errmsg=errmsg)
          if(dealloc_stat /= 0) call juDFT_error("deallocate failed for hmat or smat",&
                                             hint=errmsg, calledby="eigen.F90")
206 207

          ! Output results
208 209
          CALL timestart("EV output")
#if defined(CPP_MPI)
210
          ! Collect number of all eigenvalues
211
          ne_found=ne_all
Daniel Wortmann's avatar
Daniel Wortmann committed
212
          CALL MPI_ALLREDUCE(ne_found,ne_all,1,MPI_INTEGER,MPI_SUM, mpi%sub_comm,ierr)
213
          ne_all=MIN(DIMENSION%neigd,ne_all)
214 215 216
#else
          ne_found=ne_all
#endif          
Daniel Wortmann's avatar
Daniel Wortmann committed
217
          IF (.NOT.zmat%l_real) THEN
218
             zMat%data_c(:lapw%nmat,:ne_found) = CONJG(zMat%data_c(:lapw%nmat,:ne_found))
219
          END IF
220 221
          CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,&
                         eig(:ne_found),n_start=mpi%n_size,n_end=mpi%n_rank,zmat=zMat)
222
          neigBuffer(nk,jsp) = ne_found
223
#if defined(CPP_MPI)
224
          ! RMA synchronization
225 226
          CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
#endif
227
          CALL timestop("EV output")
228 229

          IF (banddos%unfoldband) THEN
230
               CALL calculate_plot_w_n(banddos,cell,kpts,smat_unfold,zMat,lapw,nk,jsp,eig,results,input,atoms,unfoldingBuffer,mpi)
231 232 233
	       DEALLOCATE(smat_unfold, stat=dealloc_stat, errmsg=errmsg)
          if(dealloc_stat /= 0) call juDFT_error("deallocate failed for smat_unfold",&
                                             hint=errmsg, calledby="eigen.F90")
234 235
          END IF

236 237 238 239
       END DO  k_loop
    END DO ! spin loop ends

#ifdef CPP_MPI
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
    IF (banddos%unfoldband) THEN
       write(1230+mpi%irank,*) SHAPE(results%unfolding_weights)
       write(1230+mpi%irank,*) SHAPE(unfoldingBuffer)
       DO jsp_m = 1, SIZE(unfoldingBuffer,3)
		DO i_kpt_m = 1, SIZE(unfoldingBuffer,2)
			DO i_m = 1, SIZE(unfoldingBuffer,1)
             write(1230+mpi%irank,'(2f15.8)') unfoldingBuffer(i_m, i_kpt_m,jsp_m)
			END DO
		END DO
	END DO
        write(*,*) SHAPE(results%unfolding_weights)
!       FLUSH(1230+mpi%irank)
       results%unfolding_weights = CMPLX(0.0,0.0)
       CALL MPI_ALLREDUCE(unfoldingBuffer,results%unfolding_weights,SIZE(results%unfolding_weights,1)*SIZE(results%unfolding_weights,2)*SIZE(results%unfolding_weights,3),CPP_MPI_COMPLEX,MPI_SUM,mpi%mpi_comm,ierr)
	       write(1240+mpi%irank,*) SHAPE(results%unfolding_weights)       
	DO jsp_m = 1, SIZE(results%unfolding_weights,3)
		DO i_kpt_m = 1, SIZE(results%unfolding_weights,2)
			DO i_m = 1, SIZE(results%unfolding_weights,1)
             write(1240+mpi%irank,'(2f15.8)') results%unfolding_weights(i_m, i_kpt_m,jsp_m)
			END DO
		END DO
	END DO
    END IF
263
    CALL MPI_ALLREDUCE(neigBuffer,results%neig,kpts%nkpt*input%jspins,MPI_INTEGER,MPI_SUM,mpi%sub_comm,ierr)
264
    CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
265
#else
266
    results%neig(:,:) = neigBuffer(:,:)
267
    results%unfolding_weights(:,:,:) = unfoldingBuffer(:,:,:)
268
#endif
269 270 271 272 273 274 275 276 277

    ! Sorry for the following strange workaround to fill the results%eig array.
    ! At some point someone should have a closer look at how the eigenvalues are
    ! distributed and fill the array without using the eigenvalue-IO.
    DO jsp = 1,MERGE(1,input%jspins,noco%l_noco)
       DO nk = 1,kpts%nkpt
          CALL read_eig(eig_id,nk,jsp,results%neig(nk,jsp),results%eig(:,nk,jsp))
#ifdef CPP_MPI
          CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
278
#endif
279 280 281
       END DO
    END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
282
    !IF (hybrid%l_hybrid.OR.hybrid%l_calhf) CALL close_eig(eig_id)
283

284
    IF( input%jspins .EQ. 1 .AND. hybrid%l_hybrid ) THEN
285
       results%te_hfex%valence = 2*results%te_hfex%valence
286
       IF(hybrid%l_calhf) results%te_hfex%core = 2*results%te_hfex%core
287
    END IF
288 289
    enpara%epara_min = MINVAL(enpara%el0)
    enpara%epara_min = MIN(MINVAL(enpara%ello0),enpara%epara_min)
290 291
  END SUBROUTINE eigen
END MODULE m_eigen
292