Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : !<@brief
7 : !<The eigen routine sets up and solves the generalized eigenvalue problem
8 : !More description at end of file
9 : MODULE m_eigen
10 : #ifdef CPP_MPI
11 : use mpi
12 : #endif
13 : USE m_juDFT
14 : IMPLICIT NONE
15 : CONTAINS
16 : !>The eigenvalue problem is constructed and solved in this routine. The following steps are performed:
17 : !> 1. Preparation: generate energy parameters, open eig-file
18 : !> 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
19 : !> 3. within the (collinear)spin and k-point loop: CALL to eigen_hssetup() to generate the matrices, CALL to eigen_diag() to perform diagonalization
20 : !> 4. writing (saving) of eigenvectors
21 : !>
22 : !>@author D. Wortmann
23 : !
24 : ! Modifications done to use this with DFPT phonons:
25 : ! a) We need additional MT-integrals from mt_setup that cover the potential variation V1.
26 : ! b) The eigenvalues are to be evaluated for k+q, not k.
27 : ! c) Additionally, load in the occupied states for k without q.
28 : ! d) The work isn't done once the eigenvectors and eigenvalues are found. There is post-
29 : ! processing to gain the perturbed eigenvalues and eigenvectors.
30 : ! e) The latter are the actual output for the routine when used for DFPT. They are saved
31 : ! the same way as the eigenvalues before, but for a shifted eig_id.
32 686 : SUBROUTINE eigen(fi,fmpi,stars,sphhar,xcpot,forcetheo,enpara,nococonv,&
33 : mpdata,hybdat,iter,eig_id,results,inden,v,vx,hub1data,&
34 : bqpt, hmat_out, smat_out)
35 :
36 : USE m_types
37 : USE m_constants
38 : USE m_eigen_hssetup
39 : USE m_pot_io
40 : USE m_eigen_diag
41 : !USE m_hsefunctional
42 : USE m_local_Hamiltonian
43 : USE m_util
44 : !USE m_icorrkeys
45 : USE m_eig66_io, ONLY : write_eig, read_eig
46 : USE m_xmlOutput
47 :
48 : USE m_symmetrize_matrix
49 : USE m_unfold_band_kpts !used for unfolding bands
50 : USE m_types_mpimat
51 : use m_store_load_hybrid
52 :
53 : IMPLICIT NONE
54 :
55 : type(t_fleurinput), intent(in) :: fi
56 : TYPE(t_results),INTENT(INOUT):: results
57 : CLASS(t_xcpot),INTENT(IN) :: xcpot
58 : TYPE(t_mpi),INTENT(IN) :: fmpi
59 : CLASS(t_forcetheo),INTENT(IN):: forcetheo
60 : TYPE(t_mpdata), intent(inout):: mpdata
61 : TYPE(t_hybdat), INTENT(INOUT):: hybdat
62 : TYPE(t_enpara),INTENT(INOUT) :: enpara
63 : TYPE(t_nococonv),INTENT(IN) :: nococonv
64 : TYPE(t_stars),INTENT(IN) :: stars
65 : TYPE(t_sphhar),INTENT(IN) :: sphhar
66 : TYPE(t_potden),INTENT(IN) :: inden !
67 : TYPE(t_hub1data),INTENT(INOUT):: hub1data
68 : TYPE(t_potden), INTENT(IN) :: vx
69 : TYPE(t_potden),INTENT(IN) :: v
70 :
71 : ! EXTERNAL MPI_BCAST !only used by band_unfolding to broadcast the gvec
72 :
73 : ! Scalar Arguments
74 : INTEGER,INTENT(IN) :: iter
75 : INTEGER,INTENT(IN) :: eig_id
76 :
77 : REAL, OPTIONAL, INTENT(IN) :: bqpt(3)
78 : CLASS(t_mat), OPTIONAL, INTENT(INOUT) :: hmat_out, smat_out
79 :
80 : ! Local Scalars
81 : INTEGER jsp,nk,ne_all,ne_found,neigd2,dim_mat
82 : INTEGER nk_i,n_size,n_rank
83 : LOGICAL l_needs_vectors
84 : INTEGER :: solver=0
85 : ! Local Arrays
86 : INTEGER :: ierr
87 686 : INTEGER :: neigBuffer(fi%kpts%nkpt,fi%input%jspins)
88 :
89 686 : COMPLEX :: unfoldingBuffer(SIZE(results%unfolding_weights,1),fi%kpts%nkpt,fi%input%jspins) ! needed for unfolding bandstructure fmpi case
90 :
91 1372 : INTEGER, ALLOCATABLE :: nvBuffer(:,:), nvBufferTemp(:,:), k_selection(:)
92 686 : REAL, ALLOCATABLE :: bkpt(:)
93 686 : REAL, ALLOCATABLE :: eig(:), eigBuffer(:,:,:)
94 :
95 686 : TYPE(t_tlmplm) :: td
96 686 : TYPE(t_usdus) :: ud
97 686 : TYPE(t_lapw) :: lapw
98 1372 : CLASS(t_mat), ALLOCATABLE :: zMat
99 2744 : CLASS(t_mat), ALLOCATABLE :: hmat,smat
100 686 : CLASS(t_mat), ALLOCATABLE :: smat_unfold !used for unfolding bandstructure
101 2744 : TYPE(t_kpts) :: kpts_mod!kqpts ! basically kpts, but with q added onto each one.
102 :
103 : ! Variables for HF or fi%hybinp functional calculation
104 : INTEGER :: dealloc_stat, iqdir
105 : character(len=300) :: errmsg
106 : real :: alpha_hybrid
107 :
108 : REAL :: qphon(3)
109 :
110 : !ALLOCATE(k_selection(16))
111 : !k_selection = [25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40]
112 686 : ALLOCATE(k_selection(1))
113 : !k_selection = [40,61,453]
114 2058 : k_selection = [-1]
115 :
116 :
117 : !kqpts = fi%kpts
118 1244 : kpts_mod = fi%kpts
119 :
120 686 : l_needs_vectors = .true.
121 686 : if (forcetheo%l_in_forcetheo_loop) then
122 0 : l_needs_vectors = forcetheo%l_needs_vectors
123 : endif
124 :
125 : ! Modify this from kpts only in DFPT case.
126 686 : ALLOCATE(bkpt(3))
127 686 : IF (PRESENT(bqpt)) THEN
128 0 : DO nk_i = 1, size(fmpi%k_list)
129 : !kqpts%bk(:, nk_i) = kqpts%bk(:, nk_i) + bqpt
130 0 : nk=fmpi%k_list(nk_i)
131 0 : bkpt = fi%kpts%bk(:, nk)
132 : DO iqdir = 1, 3
133 : !IF (bkpt(iqdir)+bqpt(iqdir)>=0.5) bkpt(iqdir) = bkpt(iqdir) - 1.0
134 : !IF (bkpt(iqdir)+bqpt(iqdir)<-0.5) bkpt(iqdir) = bkpt(iqdir) + 1.0
135 : !IF (bkpt(iqdir)+bqpt(iqdir)>=0.5.AND.ABS(bqpt(iqdir))>1e-8) bkpt(iqdir) = bkpt(iqdir) - 1.0
136 : !IF (bkpt(iqdir)+bqpt(iqdir)<-0.5.AND.ABS(bqpt(iqdir))>1e-8) bkpt(iqdir) = bkpt(iqdir) + 1.0
137 : END DO
138 0 : kpts_mod%bk(:, nk) = bkpt
139 : END DO
140 : END IF
141 :
142 686 : call ud%init(fi%atoms,fi%input%jspins)
143 :
144 2058 : ALLOCATE(eig(fi%input%neig))
145 3430 : ALLOCATE(eigBuffer(fi%input%neig,fi%kpts%nkpt,fi%input%jspins))
146 4610 : ALLOCATE(nvBuffer(fi%kpts%nkpt,MERGE(1,fi%input%jspins,fi%noco%l_noco)),nvBufferTemp(fi%kpts%nkpt,MERGE(1,fi%input%jspins,fi%noco%l_noco)))
147 :
148 : !IF (fmpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber '/),(/iter,v%iter/),&
149 : ! RESHAPE((/19,13,5,5/),(/2,2/)))
150 :
151 : ! Set up and solve the eigenvalue problem
152 : ! loop over spins
153 : ! set up k-point independent t(l'm',lm) matrices
154 :
155 686 : alpha_hybrid = MERGE(xcpot%get_exchange_weight(),0.0,hybdat%l_subvxc)
156 686 : CALL local_ham(sphhar,fi%atoms,fi%sym,fi%noco,nococonv,enpara,fmpi,v,vx,inden,fi%input,fi%hub1inp,hub1data,td,ud,alpha_hybrid)
157 12364 : neigBuffer = 0
158 12364 : results%neig = 0
159 406248 : results%eig = 1.0e300
160 352188 : eigBuffer = 1.0e300
161 406248 : unfoldingBuffer = CMPLX(0.0,0.0)
162 11388 : nvBuffer = 0
163 11388 : nvBufferTemp = 0
164 :
165 2070 : DO jsp = 1, MERGE(1,fi%input%jspins,fi%noco%l_noco)
166 8232 : k_loop:DO nk_i = 1,size(fmpi%k_list)
167 7342 : nk=fmpi%k_list(nk_i)
168 :
169 36710 : bkpt = kpts_mod%bk(:, nk)
170 : ! Set up lapw list
171 : !CALL lapw%init(fi%input,fi%noco,nococonv, kqpts, fi%atoms, fi%sym, nk, fi%cell, fmpi)
172 7342 : CALL lapw%init(fi%input,fi%noco,nococonv, kpts_mod, fi%atoms, fi%sym, nk, fi%cell, fmpi, bqpt)
173 :
174 7342 : call timestart("Setup of H&S matrices")
175 7342 : CALL eigen_hssetup(jsp,fmpi,fi,mpdata,results,inDen,vx,xcpot,enpara,nococonv,stars,sphhar,hybdat,ud,td,v,lapw,nk,smat,hmat)
176 7342 : CALL timestop("Setup of H&S matrices")
177 :
178 7342 : IF (PRESENT(hmat_out)) THEN
179 0 : IF (hmat_out%l_real) THEN
180 0 : dim_mat = SIZE(smat%data_r(:,1))
181 :
182 0 : CALL smat_out%init(.TRUE., dim_mat, dim_mat, fmpi%sub_comm, .false.)
183 0 : CALL hmat_out%init(smat_out)
184 :
185 0 : hmat_out%data_r(:dim_mat,:dim_mat) = hmat%data_r
186 0 : smat_out%data_r(:dim_mat,:dim_mat) = smat%data_r
187 : ELSE
188 0 : dim_mat = SIZE(smat%data_c(:,1))
189 :
190 0 : CALL smat_out%init(.FALSE., dim_mat, dim_mat, fmpi%sub_comm, .false.)
191 0 : CALL hmat_out%init(smat_out)
192 :
193 0 : hmat_out%data_c(:dim_mat,:dim_mat) = hmat%data_c
194 0 : smat_out%data_c(:dim_mat,:dim_mat) = smat%data_c
195 : END IF
196 : END IF
197 :
198 14684 : IF (ANY(nk==k_selection)) THEN
199 0 : CALL save_npy(int2str(eig_id)//"_"//int2str(nk)//"_h0.npy",hmat%data_r)
200 0 : CALL save_npy(int2str(eig_id)//"_"//int2str(nk)//"_s0.npy",smat%data_r)
201 0 : CALL save_npy(int2str(eig_id)//"_"//int2str(nk)//"_vk.npy",lapw%vk)
202 0 : CALL save_npy(int2str(eig_id)//"_"//int2str(nk)//"_gvec.npy",lapw%gvec)
203 0 : write(9530,*) nk, lapw%bkpt, qphon
204 : END IF
205 7342 : nvBuffer(nk,jsp) = lapw%nv(jsp)
206 :
207 7342 : ne_all=fi%input%neig
208 7342 : IF(ne_all < 0) ne_all = lapw%nmat
209 7342 : IF(ne_all > lapw%nmat) ne_all = lapw%nmat
210 :
211 : !Try to symmetrize matrix
212 : !CALL symmetrize_matrix(fmpi,fi%noco,kqpts,nk,hmat,smat,.FALSE.)
213 7342 : IF (.NOT.PRESENT(bqpt)) THEN
214 7342 : CALL symmetrize_matrix(fmpi,fi%noco,kpts_mod,nk,hmat,smat,.FALSE.)
215 : END IF
216 :
217 7342 : IF (fi%banddos%unfoldband .AND. (.NOT. fi%noco%l_soc)) THEN
218 : select type(smat)
219 : type is (t_mat)
220 24 : allocate(t_mat::smat_unfold)
221 24 : select type(smat_unfold)
222 : type is (t_mat)
223 24 : smat_unfold=smat
224 2082944 : smat_unfold%data_c=CONJG(smat%data_c)
225 : end select
226 : type is (t_mpimat)
227 0 : allocate(t_mpimat::smat_unfold)
228 0 : select type(smat_unfold)
229 : type is (t_mpimat)
230 0 : smat_unfold=smat
231 0 : smat_unfold%data_c=CONJG(smat%data_c)
232 : end select
233 : end select
234 : END IF
235 :
236 : ! Solve generalized eigenvalue problem.
237 : ! ne_all ... number of eigenpairs searched (and found) on this node
238 : ! on input, overall number of eigenpairs searched,
239 : ! on output, local number of eigenpairs found
240 : ! eig ...... all eigenvalues, output
241 : ! zMat ..... local eigenvectors, output
242 7342 : if (fmpi%pe_diag) THEN
243 7342 : CALL eigen_diag(solver,hmat,smat,ne_all,eig,zMat,nk,jsp,iter)
244 : ELSE
245 0 : ne_all=0
246 : endif
247 7342 : CALL smat%free()
248 7342 : CALL hmat%free()
249 22026 : DEALLOCATE(hmat,smat, stat=dealloc_stat, errmsg=errmsg)
250 7342 : if(dealloc_stat /= 0) call juDFT_error("deallocate failed for hmat or smat",&
251 0 : hint=errmsg, calledby="eigen.F90")
252 :
253 : ! Output results
254 7342 : CALL timestart("EV output")
255 7342 : ne_found=ne_all
256 7342 : if (fmpi%pe_diag) THEN
257 : #if defined(CPP_MPI)
258 : ! Collect number of all eigenvalues
259 7342 : CALL MPI_ALLREDUCE(ne_found,ne_all,1,MPI_INTEGER,MPI_SUM, fmpi%diag_sub_comm,ierr)
260 7342 : ne_all=MIN(fi%input%neig,ne_all)
261 : #endif
262 : endif
263 :
264 7342 : IF (fmpi%n_rank == 0) THEN
265 : ! Only process 0 writes out the value of ne_all and the
266 : ! eigenvalues.
267 : #ifdef CPP_MPI
268 4906 : call MPI_COMM_RANK(fmpi%diag_sub_comm,n_rank,ierr)
269 4906 : call MPI_COMM_SIZE(fmpi%diag_sub_comm,n_size,ierr)
270 : #else
271 : n_rank = 0; n_size=1;
272 : #endif
273 :
274 4906 : if (l_needs_vectors) then
275 4906 : call write_eig(eig_id, nk,jsp,ne_found,ne_all,eig(:ne_all),n_start=n_size,n_end=n_rank,zMat=zMat)
276 4906 : IF (fi%noco%l_soc .AND. fi%hybinp%l_hybrid) &
277 0 : CALL write_eig(hybdat%eig_id, nk,jsp,ne_found,ne_all,eig(:ne_all),n_start=n_size,n_end=n_rank,zMat=zMat)
278 : else
279 0 : CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,eig(:ne_all))
280 0 : IF (fi%noco%l_soc .AND. fi%hybinp%l_hybrid) &
281 0 : CALL write_eig(hybdat%eig_id, nk,jsp,ne_found,ne_all,eig(:ne_all))
282 : endif
283 158806 : eigBuffer(:ne_all,nk,jsp) = eig(:ne_all)
284 : ELSE
285 2436 : if (fmpi%pe_diag.and.l_needs_vectors) THEN
286 2436 : CALL write_eig(eig_id, nk,jsp,ne_found,n_start=fmpi%n_size,n_end=fmpi%n_rank,zMat=zMat)
287 2436 : IF (fi%noco%l_soc .AND. fi%hybinp%l_hybrid) &
288 0 : CALL write_eig(hybdat%eig_id, nk,jsp,ne_found,n_start=fmpi%n_size,n_end=fmpi%n_rank,zMat=zMat)
289 : ENDIF
290 : ENDIF
291 7342 : neigBuffer(nk,jsp) = ne_found
292 : #if defined(CPP_MPI)
293 : ! RMA synchronization
294 7342 : CALL MPI_BARRIER(fmpi%MPI_COMM,ierr)
295 : #endif
296 7342 : CALL timestop("EV output")
297 :
298 7342 : IF (fi%banddos%unfoldband .AND. (.NOT. fi%noco%l_soc)) THEN
299 24 : IF(modulo (fi%kpts%nkpt,fmpi%n_size).NE.0) call juDFT_error("number fi%kpts needs to be multiple of number fmpi threads",&
300 0 : hint=errmsg, calledby="eigen.F90")
301 24 : CALL calculate_plot_w_n(fi%banddos,fi%cell,kpts_mod,zMat,lapw,nk,jsp,eig,results,fi%input,fi%atoms,unfoldingBuffer,fmpi,fi%noco%l_soc,smat_unfold=smat_unfold)
302 24 : CALL smat_unfold%free()
303 48 : DEALLOCATE(smat_unfold, stat=dealloc_stat, errmsg=errmsg)
304 24 : if(dealloc_stat /= 0) call juDFT_error("deallocate failed for smat_unfold",&
305 0 : hint=errmsg, calledby="eigen.F90")
306 : END IF
307 :
308 15574 : IF (allocated(zmat)) THEN
309 7342 : call zMat%free()
310 14684 : deallocate(zMat)
311 : ENDIF
312 : END DO k_loop
313 : #ifdef CPP_MPI
314 : !print *,"Remaining Barriers:",size(fmpi%k_list)+1,fmpi%max_length_k_list
315 1576 : DO nk=size(fmpi%k_list)+1,fmpi%max_length_k_list
316 890 : CALL MPI_BARRIER(fmpi%MPI_COMM,ierr)
317 : ENDDO
318 : #endif
319 : END DO ! spin loop ends
320 :
321 686 : neigd2 = MIN(fi%input%neig,lapw%dim_nbasfcn())
322 : #ifdef CPP_MPI
323 686 : IF (fi%banddos%unfoldband .AND. (.NOT. fi%noco%l_soc)) THEN
324 3894 : results%unfolding_weights = CMPLX(0.0,0.0)
325 2 : CALL MPI_ALLREDUCE(unfoldingBuffer,results%unfolding_weights,SIZE(results%unfolding_weights,1)*SIZE(results%unfolding_weights,2)*SIZE(results%unfolding_weights,3),MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
326 : END IF
327 686 : CALL MPI_ALLREDUCE(neigBuffer,results%neig,fi%kpts%nkpt*fi%input%jspins,MPI_INTEGER,MPI_SUM,fmpi%mpi_comm,ierr)
328 686 : CALL MPI_ALLREDUCE(eigBuffer(:neigd2,:,:),results%eig(:neigd2,:,:),neigd2*fi%kpts%nkpt*fi%input%jspins,MPI_DOUBLE_PRECISION,MPI_MIN,fmpi%mpi_comm,ierr)
329 2058 : CALL MPI_ALLREDUCE(nvBuffer(:,:),nvBufferTemp(:,:),size(nvbuffer),MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,ierr)
330 686 : CALL MPI_BARRIER(fmpi%MPI_COMM,ierr)
331 : #else
332 : results%neig(:,:) = neigBuffer(:,:)
333 : results%eig(:neigd2,:,:) = eigBuffer(:neigd2,:,:)
334 : results%unfolding_weights(:,:,:) = unfoldingBuffer(:,:,:)
335 : nvBufferTemp(:,:) = nvBuffer(:,:)
336 : #endif
337 : !CALL save_npy(int2str(eig_id)//"_eig0.npy",results%eig(:neigd2,:,1))
338 :
339 686 : IF(fmpi%irank.EQ.0) THEN
340 343 : WRITE(oUnit,'(a)') ''
341 343 : WRITE(oUnit,'(a)') ' basis set size:'
342 343 : WRITE(oUnit,'(a)') ' jsp ikpt nv LOs overall'
343 788 : DO jsp = 1, MERGE(1,fi%input%jspins,fi%noco%l_noco)
344 5694 : DO nk = 1, fi%kpts%nkpt
345 5351 : WRITE(oUnit,'(5i8)') jsp, nk, nvBufferTemp(nk,jsp), fi%atoms%nlotot, nvBufferTemp(nk,jsp) + fi%atoms%nlotot
346 : END DO
347 : END DO
348 : END IF
349 :
350 21396 : enpara%epara_min = MINVAL(enpara%el0)
351 5616 : enpara%epara_min = MIN(MINVAL(enpara%ello0),enpara%epara_min)
352 1930 : END SUBROUTINE eigen
353 24 : END MODULE m_eigen
|
