Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2022 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 :
7 : MODULE m_dfpt_eigen
8 :
9 : #ifdef CPP_MPI
10 : USE mpi
11 : #endif
12 : USE m_juDFT
13 :
14 : #ifdef _OPENACC_later
15 : USE cublas
16 : #define CPP_zgemv cublaszgemv
17 : #else
18 : #define CPP_zgemv zgemv
19 : #endif
20 :
21 : IMPLICIT NONE
22 :
23 : CONTAINS
24 :
25 0 : SUBROUTINE dfpt_eigen(fi, sphhar, results, resultsq, results1, fmpi, enpara, nococonv, starsq, v1real, v1imag, vTot, inden, bqpt, &
26 : eig_id, q_eig_id, dfpt_eig_id, iDir, iDtype, killcont, l_real, sh_den, dfpt_eig_id2)
27 :
28 : USE m_types
29 : USE m_constants
30 : USE m_dfpt_eigen_hssetup
31 : USE m_pot_io
32 : USE m_util
33 : USE m_eig66_io, ONLY : write_eig, read_eig
34 : USE m_xmlOutput
35 : USE m_types_mpimat
36 : USE m_dfpt_tlmplm
37 : USE m_local_hamiltonian
38 :
39 : IMPLICIT NONE
40 :
41 : type(t_fleurinput), intent(in) :: fi
42 : TYPE(t_sphhar), INTENT(IN) :: sphhar
43 : TYPE(t_results),INTENT(INOUT):: results, resultsq, results1
44 : TYPE(t_mpi),INTENT(IN) :: fmpi
45 : TYPE(t_enpara),INTENT(IN) :: enpara
46 : TYPE(t_nococonv),INTENT(IN) :: nococonv
47 : TYPE(t_stars),INTENT(IN) :: starsq
48 : TYPE(t_potden),INTENT(IN) :: inden, v1real, v1imag, vTot
49 : REAL, INTENT(IN) :: bqpt(3)
50 : INTEGER, INTENT(IN) :: eig_id, q_eig_id, dfpt_eig_id, iDir, iDtype, killcont(6)
51 : LOGICAL, INTENT(IN) :: l_real, sh_den
52 : INTEGER, OPTIONAL, INTENT(IN) :: dfpt_eig_id2
53 :
54 : INTEGER n_size,n_rank
55 : INTEGER i,err,nk,jsp,nk_i,neigd2
56 :
57 : INTEGER :: ierr, iNupr
58 :
59 : REAL :: bkpt(3), q_loop(3)
60 :
61 : INTEGER :: nu
62 :
63 : LOGICAL :: old_and_wrong
64 :
65 : COMPLEX :: wtfq
66 :
67 0 : TYPE(t_tlmplm) :: td, tdV1
68 0 : TYPE(t_potden) :: vx
69 0 : TYPE(t_hub1data) :: hub1data
70 0 : TYPE(t_usdus) :: ud
71 0 : TYPE(t_lapw) :: lapw, lapwq
72 0 : CLASS(t_mat), ALLOCATABLE :: zMatk, zMatq, zMat1, zMat2
73 0 : CLASS(t_mat), ALLOCATABLE :: hmat,smat
74 :
75 : INTEGER :: dealloc_stat, nbasfcnq, nbasfcn, neigk, neigq, noccbd, noccbdq, noccbdmin
76 : character(len=300) :: errmsg
77 0 : INTEGER, ALLOCATABLE :: ev_list(:), q_ev_list(:)
78 0 : COMPLEX, ALLOCATABLE :: tempVec(:), tempMat1(:), tempMat2(:), z1H(:,:), z1S(:,:), tempMat3(:), z1H2(:,:), z1S2(:,:)
79 0 : REAL, ALLOCATABLE :: eigk(:), eigq(:), eigs1(:), eigBuffer(:,:,:)
80 :
81 : COMPLEX zdotc
82 : EXTERNAL zdotc
83 :
84 :
85 0 : old_and_wrong = .FALSE.
86 :
87 0 : CALL vx%copyPotDen(vTot)
88 0 : ALLOCATE(vx%pw_w, mold=vx%pw)
89 0 : vx%pw_w = vTot%pw_w
90 :
91 : ! Get the (lm) matrix elements for V1 and H0
92 0 : CALL ud%init(fi%atoms,fi%input%jspins)
93 0 : CALL dfpt_tlmplm(fi%atoms,fi%sym,sphhar,fi%input,fi%noco,enpara,fi%hub1inp,hub1data,vTot,fmpi,tdV1,v1real,v1imag,.FALSE.)
94 0 : CALL local_ham(sphhar,fi%atoms,fi%sym,fi%noco,nococonv,enpara,fmpi,vTot,vx,inden,fi%input,fi%hub1inp,hub1data,td,ud,0.0,.TRUE.)
95 :
96 0 : ALLOCATE(eigBuffer(fi%input%neig,fi%kpts%nkpt,fi%input%jspins))
97 0 : eigBuffer = 0.0
98 0 : results1%eig = 1.0e300
99 :
100 0 : DO jsp = 1, MERGE(1,fi%input%jspins,fi%noco%l_noco)
101 0 : k_loop:DO nk_i = 1,size(fmpi%k_list)
102 0 : nk=fmpi%k_list(nk_i)
103 :
104 : ! Get the required eigenvectors and values at k for occupied bands:
105 0 : bkpt = fi%kpts%bk(:, nk)
106 :
107 0 : q_loop = bqpt
108 :
109 0 : CALL lapw%init(fi%input, fi%noco, nococonv, fi%kpts, fi%atoms, fi%sym, nk, fi%cell, fmpi)
110 0 : CALL lapwq%init(fi%input, fi%noco, nococonv, fi%kpts, fi%atoms, fi%sym, nk, fi%cell, fmpi, q_loop)
111 :
112 0 : noccbd = COUNT(results%w_iks(:,nk,jsp)*2.0/fi%input%jspins>1.e-8)
113 0 : noccbdq = COUNT(resultsq%w_iks(:,nk,jsp)*2.0/fi%input%jspins>1.e-8)
114 :
115 0 : noccbdmin = MIN(noccbdq,noccbd)
116 :
117 0 : nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*fi%atoms%nlotot,lapw%nv(1)+fi%atoms%nlotot,fi%noco%l_noco)
118 0 : nbasfcnq = MERGE(lapwq%nv(1)+lapwq%nv(2)+2*fi%atoms%nlotot,lapwq%nv(1)+fi%atoms%nlotot,fi%noco%l_noco)
119 :
120 0 : IF (fmpi%n_size == 1) THEN
121 0 : ALLOCATE (t_mat::zMatk)
122 0 : ALLOCATE (t_mat::zMatq)
123 : ELSE
124 0 : ALLOCATE (t_mpimat::zMatk)
125 0 : ALLOCATE (t_mpimat::zMatq)
126 : END IF
127 :
128 : ! Initialize the expansion coefficient matrices at k and k+q
129 : ! Then read all the stuff into it
130 0 : CALL zMatk%init(l_real,nbasfcn,noccbd)
131 0 : CALL zMatq%init(l_real,nbasfcnq,nbasfcnq)
132 :
133 0 : ALLOCATE(ev_list(noccbd))
134 0 : ev_list = (/(i, i=1,noccbd, 1)/)
135 0 : ALLOCATE(q_ev_list(nbasfcnq))
136 0 : q_ev_list = (/(i, i=1,nbasfcnq, 1)/)
137 :
138 0 : ALLOCATE(eigk(noccbd))
139 0 : ALLOCATE(eigq(nbasfcnq))
140 0 : ALLOCATE(eigs1(noccbd))
141 :
142 0 : CALL timestart("Read eigenstuff at k/k+q")
143 0 : CALL read_eig(eig_id, nk, jsp, list=ev_list, neig=neigk, eig=eigk, zmat=zMatk)
144 0 : CALL read_eig(q_eig_id, nk, jsp, list=q_ev_list, neig=neigq, eig=eigq, zmat=zMatq)
145 0 : CALL timestop("Read eigenstuff at k/k+q")
146 :
147 : ! Construct the perturbed Hamiltonian and Overlap matrix perturbations:
148 0 : CALL timestart("Setup of matrix perturbations")
149 0 : CALL dfpt_eigen_hssetup(jsp,fmpi,fi,enpara,nococonv,starsq,ud,td,tdV1,vTot,v1real,lapw,lapwq,iDir,iDtype,hmat,smat,nk,killcont)
150 0 : CALL timestop("Setup of matrix perturbations")
151 :
152 0 : IF (fmpi%n_size == 1) THEN
153 0 : ALLOCATE (t_mat::zMat1)
154 : ELSE
155 0 : ALLOCATE (t_mpimat::zMat1)
156 : END IF
157 :
158 : ! Initialize the expansion coefficient perturbation matrix
159 0 : CALL zMat1%init(.FALSE.,nbasfcnq,noccbd)
160 :
161 0 : ALLOCATE(z1H,mold=zMat1%data_c)
162 0 : ALLOCATE(z1S,mold=zMat1%data_c)
163 0 : z1H = CMPLX(0.0,0.0)
164 0 : z1S = CMPLX(0.0,0.0)
165 :
166 : ! For the dynmat calculation: initialize the auxiliary coefficients
167 0 : IF (.NOT.sh_den.AND..NOT.old_and_wrong) THEN
168 0 : IF (fmpi%n_size == 1) THEN
169 0 : ALLOCATE (t_mat::zMat2)
170 : ELSE
171 0 : ALLOCATE (t_mpimat::zMat2)
172 : END IF
173 :
174 0 : CALL zMat2%init(.FALSE.,nbasfcnq,noccbd)
175 :
176 0 : ALLOCATE(z1H2,mold=zMat2%data_c)
177 0 : ALLOCATE(z1S2,mold=zMat2%data_c)
178 0 : z1H2 = CMPLX(0.0,0.0)
179 0 : z1S2 = CMPLX(0.0,0.0)
180 : END IF
181 :
182 : ! Allocate auxiliary quantities
183 0 : ALLOCATE(tempVec(nbasfcnq))
184 0 : ALLOCATE(tempMat1(nbasfcnq))
185 0 : ALLOCATE(tempMat2(nbasfcnq))
186 0 : IF (.NOT.sh_den.AND..NOT.old_and_wrong) ALLOCATE(tempMat3(nbasfcnq))
187 :
188 0 : CALL timestart("Matrix multiplications")
189 0 : DO nu = 1, noccbd
190 0 : eigs1(nu) = 0.0
191 :
192 : ! TODO: At the moment H1 and S1 are handled seperately; this was
193 : ! for debugging purposes; possibly revert for efficiency.
194 0 : IF (l_real) THEN ! l_real for zMatk
195 0 : tempVec(:nbasfcnq) = MATMUL(hmat%data_c,zMatk%data_r(:nbasfcn,nu))
196 : ELSE
197 0 : CALL CPP_zgemv('N',nbasfcnq,nbasfcn,CMPLX(1.0,0.0),hmat%data_c,nbasfcnq,zMatk%data_c(:nbasfcn,nu),1,CMPLX(0.0,0.0),tempVec,1)
198 : !tempVec(:nbasfcnq) = MATMUL(hmat%data_c,zMatk%data_c(:nbasfcn,nu))
199 : END IF
200 :
201 0 : IF (norm2(q_loop).LT.1e-8) THEN
202 0 : IF (nbasfcnq.NE.nbasfcn) CALL juDFT_error("nbasfcnq/=nbasfcn for q=0", calledby="dfpt_eigen.F90")
203 0 : IF (l_real) THEN
204 0 : eigs1(nu) = REAL(DOT_PRODUCT(zMatk%data_r(:nbasfcn,nu),tempVec))
205 : ELSE
206 0 : eigs1(nu) = REAL(zdotc(nbasfcn,zMatk%data_c(:nbasfcn,nu),1,tempVec,1))
207 : !eigs1(nu) = REAL(DOT_PRODUCT(zMatk%data_c(:nbasfcn,nu),tempVec))
208 : END IF
209 : ELSE
210 : eigs1(nu) = 0
211 : END IF
212 :
213 0 : IF (zMatq%l_real) THEN ! l_real for zMatq
214 0 : tempMat1(:nbasfcnq) = MATMUL(TRANSPOSE(zMatq%data_r),tempvec)
215 : ELSE
216 0 : CALL CPP_zgemv('C',nbasfcnq,nbasfcnq,CMPLX(1.0,0.0),zmatq%data_c,nbasfcnq,tempvec,1,CMPLX(0.0,0.0),tempMat1,1)
217 : !tempMat1(:nbasfcnq) = MATMUL(CONJG(TRANSPOSE(zMatq%data_c)),tempvec)
218 : END IF
219 :
220 : ! tempMat1 = H^{(1}_{\nu'\nu}
221 0 : DO iNupr = 1, nbasfcnq
222 0 : IF (.NOT.sh_den.AND.old_and_wrong) THEN
223 : IF (norm2(bqpt)<1e-8.AND.iNupr==nu) THEN
224 : tempMat2(iNupr) = 0.0
225 : ELSE IF (ABS(eigq(iNupr)-eigk(nu))<fi%juPhon%eDiffCut) THEN
226 : tempMat2(iNupr) = 0.0
227 : ELSE
228 : tempMat2(iNupr) = 1.0/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr)
229 : END IF
230 0 : ELSE IF (.NOT.sh_den.AND..NOT.old_and_wrong) THEN
231 0 : IF (norm2(bqpt)<1e-8.AND.iNupr==nu) THEN
232 0 : tempMat2(iNupr) = 0.0
233 0 : tempMat3(iNupr) = 0.0
234 0 : ELSE IF (ABS(eigq(iNupr)-eigk(nu))<fi%juPhon%eDiffCut) THEN
235 0 : tempMat2(iNupr) = 0.0
236 : ! Additional correction term that constitutes new
237 : ! coefficients:
238 0 : tempMat3(iNupr) = 0.5 * tempMat1(iNupr)
239 : ! TODO: This part of the correction had no effect whatsoever yet.
240 : ! Reactivate for misbehaving materials and see if there are
241 : ! changes.
242 0 : ELSE IF (iNuPr<=noccbdmin.AND.nu<=noccbdmin) THEN
243 0 : wtfq = resultsq%w_iks(iNupr,nk,jsp)/fi%kpts%wtkpt(nk)
244 : tempMat2(iNupr) = 1.0/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr) &
245 0 : & *(1.0-wtfq)
246 : ! Additional correction term that constitutes new
247 : ! coefficients:
248 0 : tempMat3(iNupr) = 0.5 * tempMat1(iNupr) * wtfq
249 : ELSE
250 0 : tempMat2(iNupr) = 1.0/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr)
251 0 : tempMat3(iNupr) = 0.0
252 : END IF
253 : ELSE
254 0 : IF (norm2(bqpt)<1e-8.AND.iNupr==nu) THEN
255 0 : tempMat2(iNupr) = 0.0
256 0 : ELSE IF (ABS(eigq(iNupr)-eigk(nu))<fi%juPhon%eDiffCut) THEN
257 0 : tempMat2(iNupr) = 0.0
258 0 : ELSE IF (iNuPr<=noccbdmin.AND.nu<=noccbdmin) THEN
259 0 : wtfq = resultsq%w_iks(iNupr,nk,jsp)/fi%kpts%wtkpt(nk)
260 : tempMat2(iNupr) = 1.0/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr) &
261 0 : & *(1.0-wtfq)
262 : ELSE
263 0 : tempMat2(iNupr) = 1.0/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr)
264 : END IF
265 : END IF
266 : END DO
267 :
268 0 : IF (zMatq%l_real) THEN
269 0 : z1H(:nbasfcnq,nu) = -MATMUL(zMatq%data_r,tempMat2(:nbasfcnq))
270 0 : IF (.NOT.sh_den.AND..NOT.old_and_wrong) z1H2(:nbasfcnq,nu) = -MATMUL(zMatq%data_r,tempMat3(:nbasfcnq))
271 : ELSE
272 0 : CALL CPP_zgemv('N',nbasfcnq,nbasfcnq,CMPLX(-1.0,0.0),zMatq%data_c,nbasfcnq,tempMat2,1,CMPLX(0.0,0.0),z1H(:nbasfcnq,nu),1)
273 : !z1H(:nbasfcnq,nu) = -MATMUL(zMatq%data_c,tempMat2(:nbasfcnq))
274 0 : IF (.NOT.sh_den.AND..NOT.old_and_wrong) CALL CPP_zgemv('N',nbasfcnq,nbasfcnq,CMPLX(-1.0,0.0),zmatq%data_c,nbasfcnq,tempMat3,1,CMPLX(0.0,0.0),z1H2(:nbasfcnq,nu),1)
275 : !IF (.NOT.sh_den.AND..NOT.old_and_wrong) z1H2(:nbasfcnq,nu) = -MATMUL(zMatq%data_c,tempMat3(:nbasfcnq))
276 : END IF
277 :
278 0 : IF (l_real) THEN ! l_real for zMatk
279 0 : tempVec(:nbasfcnq) = MATMUL(smat%data_c,zMatk%data_r(:nbasfcn,nu))
280 : ELSE
281 0 : CALL CPP_zgemv('N',nbasfcnq,nbasfcn,CMPLX(1.0,0.0),smat%data_c,nbasfcnq,zMatk%data_c(:nbasfcn,nu),1,CMPLX(0.0,0.0),tempVec,1)
282 : !tempVec(:nbasfcnq) = MATMUL(smat%data_c,zMatk%data_c(:nbasfcn,nu))
283 : END IF
284 :
285 0 : IF (norm2(q_loop).LT.1e-8) THEN
286 0 : IF (nbasfcnq.NE.nbasfcn) CALL juDFT_error("nbasfcnq/=nbasfcn for q=0", calledby="dfpt_eigen.F90")
287 0 : IF (l_real) THEN
288 0 : eigs1(nu) = eigs1(nu) - eigk(nu)*REAL(DOT_PRODUCT(zMatk%data_r(:nbasfcn,nu),tempVec))
289 : ELSE
290 0 : eigs1(nu) = eigs1(nu) - eigk(nu)*REAL(zdotc(nbasfcn,zMatk%data_c(:nbasfcn,nu),1,tempVec,1))
291 : !eigs1(nu) = eigs1(nu) - eigk(nu)*REAL(DOT_PRODUCT(zMatk%data_c(:nbasfcn,nu),tempVec))
292 : END IF
293 : ELSE
294 0 : eigs1(nu) = 0
295 : END IF
296 :
297 0 : IF (zMatq%l_real) THEN ! l_real for zMatq
298 0 : tempMat1(:nbasfcnq) = MATMUL(TRANSPOSE(zMatq%data_r),tempvec)
299 : ELSE
300 0 : CALL CPP_zgemv('C',nbasfcnq,nbasfcnq,CMPLX(1.0,0.0),zmatq%data_c,nbasfcnq,tempvec,1,CMPLX(0.0,0.0),tempMat1,1)
301 : !tempMat1(:nbasfcnq) = MATMUL(CONJG(TRANSPOSE(zMatq%data_c)),tempvec)
302 : END IF
303 :
304 : ! tempMat1 = S^{(1}_{\nu'\nu}
305 0 : DO iNupr = 1, nbasfcnq
306 0 : IF (.NOT.sh_den.AND.old_and_wrong) THEN
307 : IF (norm2(bqpt)<1e-8.AND.iNupr==nu) THEN
308 : tempMat2(iNupr) = 0.0
309 : ELSE IF (ABS(eigq(iNupr)-eigk(nu))<fi%juPhon%eDiffCut) THEN
310 : tempMat2(iNupr) = 0.5*tempMat1(iNupr)
311 : ELSE
312 : tempMat2(iNupr) = -eigk(nu)/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr)
313 : END IF
314 0 : ELSE IF (.NOT.sh_den.AND..NOT.old_and_wrong) THEN
315 0 : IF (norm2(bqpt)<1e-8.AND.iNupr==nu) THEN
316 0 : tempMat2(iNupr) = 0.0
317 0 : tempMat3(iNupr) = 0.0
318 0 : ELSE IF (ABS(eigq(iNupr)-eigk(nu))<fi%juPhon%eDiffCut) THEN
319 0 : tempMat2(iNupr) = 0.5 * tempMat1(iNupr)
320 : ! Additional correction term that constitutes new
321 : ! coefficients:
322 0 : tempMat3(iNupr) = -0.5 * eigk(nu) * tempMat1(iNupr)
323 : ! TODO: This part of the correction had no effect whatsoever yet.
324 : ! Reactivate for misbehaving materials and see if there are
325 : ! changes.
326 0 : ELSE IF (iNuPr<=noccbdmin.AND.nu<=noccbdmin) THEN
327 0 : wtfq = resultsq%w_iks(iNupr,nk,jsp)/fi%kpts%wtkpt(nk)
328 : tempMat2(iNupr) = -eigk(nu)/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr) &
329 : & *(1.0-wtfq) &
330 0 : & +0.5*tempMat1(iNupr)*wtfq
331 : ! Additional correction term that constitutes new
332 : ! coefficients:
333 0 : tempMat3(iNupr) = -0.5 * eigq(iNupr) * tempMat1(iNupr) * wtfq
334 : ELSE
335 0 : tempMat2(iNupr) = -eigk(nu)/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr)
336 0 : tempMat3(iNupr) = 0.0
337 : END IF
338 : ELSE
339 0 : IF (norm2(bqpt)<1e-8.AND.iNupr==nu) THEN
340 0 : tempMat2(iNupr) = 0.0
341 0 : ELSE IF (ABS(eigq(iNupr)-eigk(nu))<fi%juPhon%eDiffCut) THEN
342 0 : tempMat2(iNupr) = 0.5*tempMat1(iNupr)
343 0 : ELSE IF (iNuPr<=noccbdmin.AND.nu<=noccbdmin) THEN
344 0 : wtfq = resultsq%w_iks(iNupr,nk,jsp)/fi%kpts%wtkpt(nk)
345 :
346 : tempMat2(iNupr) = -eigk(nu)/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr) &
347 : & *(1.0-wtfq) &
348 0 : & +0.5*tempMat1(iNupr)*wtfq
349 : ELSE
350 0 : tempMat2(iNupr) = -eigk(nu)/(eigq(iNupr)-eigk(nu))*tempMat1(iNupr)
351 : END IF
352 : END IF
353 : END DO
354 :
355 0 : IF (zMatq%l_real) THEN
356 0 : z1S(:nbasfcnq,nu) = -MATMUL(zMatq%data_r,tempMat2(:nbasfcnq))
357 0 : IF (.NOT.sh_den.AND..NOT.old_and_wrong) z1S2(:nbasfcnq,nu) = -MATMUL(zMatq%data_r,tempMat3(:nbasfcnq))
358 : ELSE
359 0 : CALL CPP_zgemv('N',nbasfcnq,nbasfcnq,CMPLX(-1.0,0.0),zmatq%data_c,nbasfcnq,tempMat2,1,CMPLX(0.0,0.0),z1S(:nbasfcnq,nu),1)
360 : !z1S(:nbasfcnq,nu) = -MATMUL(zMatq%data_c,tempMat2(:nbasfcnq))
361 0 : IF (.NOT.sh_den.AND..NOT.old_and_wrong) CALL CPP_zgemv('N',nbasfcnq,nbasfcnq,CMPLX(-1.0,0.0),zmatq%data_c,nbasfcnq,tempMat3,1,CMPLX(0.0,0.0),z1S2(:nbasfcnq,nu),1)
362 : !IF (.NOT.sh_den.AND..NOT.old_and_wrong) z1S2(:nbasfcnq,nu) = -MATMUL(zMatq%data_c,tempMat3(:nbasfcnq))
363 : END IF
364 :
365 0 : zMat1%data_c(:nbasfcnq,nu) = z1H(:nbasfcnq,nu) + z1S(:nbasfcnq,nu)
366 0 : IF (.NOT.sh_den.AND..NOT.old_and_wrong) zMat2%data_c(:nbasfcnq,nu) = z1H2(:nbasfcnq,nu) + z1S2(:nbasfcnq,nu)
367 : END DO
368 :
369 0 : results1%neig = results%neig
370 :
371 0 : CALL timestop("Matrix multiplications")
372 :
373 0 : CALL smat%free()
374 0 : CALL hmat%free()
375 0 : DEALLOCATE(hmat,smat, stat=dealloc_stat, errmsg=errmsg)
376 0 : IF(dealloc_stat /= 0) CALL juDFT_error("Deallocation failed for hmat or smat", hint=errmsg, calledby="dfpt_eigen.F90")
377 :
378 : !#ifdef CPP_MPI
379 : !CALL MPI_BARRIER(fmpi%mpi_comm,iErr) ! Synchronizes the RMA operations
380 : !#endif
381 :
382 : ! Output results
383 0 : CALL timestart("EV1 output")
384 :
385 0 : IF (fmpi%n_rank == 0) THEN
386 : #ifdef CPP_MPI
387 0 : CALL MPI_COMM_RANK(fmpi%diag_sub_comm,n_rank,err)
388 0 : CALL MPI_COMM_SIZE(fmpi%diag_sub_comm,n_size,err)
389 : #else
390 : n_rank = 0; n_size=1;
391 : #endif
392 : CALL write_eig(dfpt_eig_id, nk, jsp, noccbd, noccbd, &
393 0 : eigs1(:noccbd), n_start=n_size,n_end=n_rank,zMat=zMat1)
394 0 : eigBuffer(:noccbd,nk,jsp) = eigs1(:noccbd)
395 0 : IF (.NOT.sh_den.AND..NOT.old_and_wrong) CALL write_eig(dfpt_eig_id2, nk, jsp, noccbd, noccbd, &
396 0 : eigs1(:noccbd), n_start=n_size,n_end=n_rank,zMat=zMat2)
397 : ELSE
398 0 : IF (fmpi%pe_diag) CALL write_eig(dfpt_eig_id, nk, jsp, noccbd, &
399 0 : n_start=fmpi%n_size,n_end=fmpi%n_rank,zMat=zMat1)
400 0 : IF ((.NOT.sh_den).AND.(.NOT.old_and_wrong).AND.fmpi%pe_diag) CALL write_eig(dfpt_eig_id2, nk, jsp, noccbd, &
401 0 : n_start=fmpi%n_size,n_end=fmpi%n_rank,zMat=zMat2)
402 : END IF
403 :
404 : #if defined(CPP_MPI)
405 : ! RMA synchronization
406 0 : CALL MPI_BARRIER(fmpi%MPI_COMM,ierr)
407 : #endif
408 0 : CALL timestop("EV1 output")
409 :
410 0 : IF (ALLOCATED(ev_list)) DEALLOCATE(ev_list)
411 0 : IF (ALLOCATED(q_ev_list)) DEALLOCATE(q_ev_list)
412 0 : IF (ALLOCATED(eigk)) DEALLOCATE(eigk)
413 0 : IF (ALLOCATED(eigq)) DEALLOCATE(eigq)
414 0 : IF (ALLOCATED(eigs1)) DEALLOCATE(eigs1)
415 0 : IF (ALLOCATED(tempVec)) DEALLOCATE(tempVec)
416 0 : IF (ALLOCATED(tempMat1)) DEALLOCATE(tempMat1)
417 0 : IF (ALLOCATED(tempMat2)) DEALLOCATE(tempMat2)
418 0 : IF (ALLOCATED(tempMat3)) DEALLOCATE(tempMat3)
419 0 : IF (ALLOCATED(z1H)) DEALLOCATE(z1H)
420 0 : IF (ALLOCATED(z1S)) DEALLOCATE(z1S)
421 0 : IF (ALLOCATED(z1H2)) DEALLOCATE(z1H2)
422 0 : IF (ALLOCATED(z1S2)) DEALLOCATE(z1S2)
423 0 : IF (ALLOCATED(zmatk)) THEN
424 0 : CALL zMatk%free()
425 0 : DEALLOCATE(zMatk)
426 : END IF
427 0 : IF (ALLOCATED(zmatq)) THEN
428 0 : CALL zMatq%free()
429 0 : DEALLOCATE(zMatq)
430 : END IF
431 : IF (ALLOCATED(zmat1)) THEN
432 0 : CALL zMat1%free()
433 0 : DEALLOCATE(zMat1)
434 : END IF
435 0 : IF (ALLOCATED(zmat2)) THEN
436 0 : CALL zMat2%free()
437 0 : DEALLOCATE(zMat2)
438 : END IF
439 : END DO k_loop
440 : END DO ! spin loop ends
441 0 : neigd2 = MIN(fi%input%neig,lapw%dim_nbasfcn())
442 : #ifdef CPP_MPI
443 0 : CALL MPI_ALLREDUCE(eigBuffer(:neigd2,:,:),results1%eig(:neigd2,:,:),neigd2*fi%kpts%nkpt*fi%input%jspins,MPI_DOUBLE_PRECISION,MPI_SUM,fmpi%mpi_comm,ierr)
444 0 : CALL MPI_BARRIER(fmpi%MPI_COMM,ierr)
445 : #else
446 : results1%eig(:neigd2,:,:) = eigBuffer(:neigd2,:,:)
447 : #endif
448 :
449 0 : END SUBROUTINE dfpt_eigen
450 0 : END MODULE m_dfpt_eigen
|
