Line data Source code
1 : MODULE m_alineso
2 : USE m_juDFT
3 : !----------------------------------------------------------------------
4 : ! Set up SO-hamiltonian for 2nd variation (hsohelp and hsoham) and
5 : ! diagonalize by lapack routines.
6 : ! Eigenvalues and vectors (eig_so and zso) are returned
7 : !----------------------------------------------------------------------
8 : CONTAINS
9 548 : SUBROUTINE alineso(eig_id,lapw,fmpi,atoms,sym,kpts,input,noco,nococonv,&
10 548 : cell , nk, usdus,rsoc,nsize,nmat, eig_so,zso)
11 :
12 : USE m_types
13 : USE m_constants
14 : USE m_hsohelp
15 : USE m_hsoham
16 : USE m_eig66_io, ONLY : read_eig
17 : IMPLICIT NONE
18 : TYPE(t_mpi),INTENT(IN) :: fmpi
19 : TYPE(t_lapw),INTENT(IN) :: lapw
20 :
21 :
22 : TYPE(t_input),INTENT(IN) :: input
23 : TYPE(t_noco),INTENT(IN) :: noco
24 : TYPE(t_nococonv),INTENT(IN) :: nococonv
25 : TYPE(t_sym),INTENT(IN) :: sym
26 : TYPE(t_cell),INTENT(IN) :: cell
27 : TYPE(t_atoms),INTENT(IN) :: atoms
28 : TYPE(t_kpts),INTENT(IN) :: kpts
29 : TYPE(t_usdus),INTENT(IN) :: usdus
30 : TYPE(t_rsoc),INTENT(IN) :: rsoc
31 : ! ..
32 : ! .. Scalar Arguments ..
33 : INTEGER, INTENT (IN) :: eig_id
34 : INTEGER, INTENT (IN) :: nk
35 : INTEGER, INTENT (OUT):: nsize,nmat
36 : ! ..
37 : ! .. Array Arguments ..
38 : COMPLEX, INTENT (OUT) :: zso(:,:,:)!(lapw%dim_nbasfcn(),2*input%neig,wannierspin)
39 : REAL, INTENT (OUT) :: eig_so(2*input%neig)
40 : !-odim
41 : !+odim
42 : ! ..
43 : ! .. Local Scalars ..
44 : REAL r2
45 : INTEGER i,i1 ,j,jsp,jsp1,k,ne,nn,nn1,nrec,info
46 : INTEGER idim_c,idim_r,jsp2,nbas,j1,ierr
47 : CHARACTER vectors
48 : LOGICAL l_socvec,l_qsgw,l_open
49 : INTEGER irec,irecl_qsgw
50 : INTEGER nat_l, extra, nat_start, nat_stop
51 : COMPLEX cdum
52 : ! ..
53 : ! .. Local Arrays ..
54 : INTEGER :: nsz(2)
55 548 : REAL :: eig(input%neig,input%jspins),s(3)
56 548 : REAL, ALLOCATABLE :: rwork(:)
57 548 : COMPLEX,ALLOCATABLE :: cwork(:),chelp(:,:,:,:,:)
58 548 : COMPLEX,ALLOCATABLE :: ahelp(:,:,:,:),bhelp(:,:,:,:)
59 548 : COMPLEX,ALLOCATABLE :: zhelp1(:,:),zhelp2(:,:)
60 548 : COMPLEX,ALLOCATABLE :: hso(:,:),hsomtx(:,:,:,:)
61 548 : COMPLEX,ALLOCATABLE :: sigma_xc_apw(:,:),sigma_xc(:,:)
62 :
63 6406 : TYPE(t_mat)::zmat(input%jspins)
64 : ! ..
65 : ! .. External Subroutines ..
66 : EXTERNAL zheev
67 :
68 : ! read from eigenvalue and -vector file
69 : !
70 :
71 1610 : zmat%l_real=input%l_real
72 1610 : zMat(1:input%jspins)%matsize1=lapw%nv(1:input%jspins)+atoms%nlotot
73 1610 : zmat%matsize2=input%neig
74 :
75 548 : INQUIRE (4649,opened=l_socvec)
76 548 : INQUIRE (file='fleur.qsgw',exist=l_qsgw)
77 548 : if (input%l_real) THEN
78 16 : ALLOCATE (zmat(1)%data_r(zmat(1)%matsize1,input%neig) )
79 38516 : zmat(1)%data_r(:,:)= 0.
80 4 : if (size(zmat)==2)THEN
81 8 : ALLOCATE(zmat(2)%data_r(zmat(2)%matsize1,input%neig) )
82 8322 : zmat(2)%data_r=0.0
83 : ENDIF
84 : else
85 2176 : ALLOCATE (zmat(1)%data_c(zmat(1)%matsize1,input%neig) )
86 4044412 : zmat(1)%data_c(:,:)= 0.
87 544 : if (size(zmat)==2)THEN
88 2048 : ALLOCATE(zmat(2)%data_c(zmat(2)%matsize1,input%neig) )
89 3566912 : zmat(2)%data_c=0.0
90 : ENDIF
91 : endif
92 15315810 : zso(:,:,:)= CMPLX(0.,0.)
93 :
94 1610 : DO jsp = 1,input%jspins
95 1062 : CALL read_eig(eig_id,nk,jsp, neig=ne,eig=eig(:,jsp))
96 1062 : IF (judft_was_argument("-debugtime")) THEN
97 0 : WRITE(oUnit,*) "Non-SOC ev for nk,jsp:",nk,jsp
98 0 : WRITE(oUnit,"(6(f10.6,1x))") eig(:ne,jsp)
99 : ENDIF
100 110914 : CALL read_eig(eig_id,nk,jsp,list=[(i,i=1,ne)],zmat=zmat(jsp))
101 :
102 : ! write(*,*) 'process',irank,' reads ',nk
103 :
104 : !!$ DO i = 1, lapw%nv(1)
105 : !!$ s = lapw%bkpt +lapw%gvec(:,i,1)
106 : !!$ r2 = DOT_PRODUCT(s,MATMUL(s,cell%bbmat))
107 : !!$ lapw%rk(i,1) = SQRT(r2)
108 : !!$ ENDDO
109 :
110 1062 : IF (ne.GT.input%neig) THEN
111 0 : WRITE (oUnit,'(a13,i4,a8,i4)') 'alineso: ne=',ne,' > input%neig=',input%neig
112 0 : CALL juDFT_error("alineso: ne > neigd",calledby="alineso")
113 : ENDIF
114 1610 : nsz(jsp) = ne
115 : ENDDO
116 : #ifdef CPP_MPI
117 548 : CALL MPI_BARRIER(fmpi%sub_comm,ierr) ! Synchronizes the RMA operations
118 : #endif
119 : !
120 : ! set up size of e.v. problem in second variation: nsize
121 : !
122 548 : nsize = 0
123 1610 : DO jsp = 1,input%jspins
124 1610 : IF (input%jspins.EQ.1) THEN
125 34 : nsize = 2*nsz(jsp)
126 34 : nsz(2) = nsz(1)
127 : ELSE
128 1028 : nsize = nsize + nsz(jsp)
129 : ENDIF
130 : ENDDO
131 : !
132 : ! distribution of (abc)cof over atoms
133 : !
134 : !
135 : ! in case of ev-parallelization, now distribute the atoms:
136 : !
137 548 : IF (fmpi%n_size > 1) THEN
138 544 : nat_l=ceiling(real(atoms%ntype)/fmpi%n_size) !stride in types
139 544 : nat_start=fmpi%n_rank*nat_l+1 !now start and stop for types
140 544 : nat_stop=(fmpi%n_rank+1)*nat_l
141 544 : if (nat_start>atoms%ntype.or.nat_stop>atoms%ntype) THEN
142 19 : nat_start=1
143 19 : nat_stop=0
144 : ELSE
145 778 : nat_start=sum(atoms%neq(:nat_start-1))+1 !convert to na
146 1303 : nat_stop=sum(atoms%neq(:nat_stop))
147 : ENDIF
148 : ELSE
149 4 : nat_start = 1
150 4 : nat_stop = atoms%nat
151 : ENDIF
152 548 : nat_l = nat_stop - nat_start + 1
153 :
154 : ! set up A and B coefficients
155 3288 : ALLOCATE (ahelp(atoms%lmaxd*(atoms%lmaxd+2),nat_l,input%neig,input%jspins))
156 2740 : ALLOCATE (bhelp(atoms%lmaxd*(atoms%lmaxd+2),nat_l,input%neig,input%jspins))
157 3836 : ALLOCATE (chelp(-atoms%llod :atoms%llod, input%neig,atoms%nlod,nat_l,input%jspins))
158 :
159 548 : CALL timestart("alineso SOC: -help")
160 : CALL hsohelp(atoms,sym,input,lapw,nsz,cell,zmat,usdus,&
161 548 : zso,noco,nococonv,nat_start,nat_stop,nat_l,ahelp,bhelp,chelp)
162 548 : CALL timestop("alineso SOC: -help")
163 :
164 : ! set up hamilton matrix
165 548 : CALL timestart("alineso SOC: -ham")
166 : #ifdef CPP_MPI
167 548 : CALL MPI_BARRIER(fmpi%MPI_COMM,ierr)
168 : #endif
169 3288 : ALLOCATE (hsomtx(input%neig,input%neig,2,2))
170 : CALL hsoham(atoms,noco,input,nsz,input%neig,chelp,rsoc,ahelp,bhelp,&
171 548 : nat_start,nat_stop,fmpi%n_rank,fmpi%n_size,fmpi%SUB_COMM,hsomtx)
172 548 : DEALLOCATE (ahelp,bhelp,chelp)
173 548 : CALL timestop("alineso SOC: -ham")
174 548 : IF (fmpi%n_rank==0) THEN
175 :
176 : ! add e.v. on diagonal
177 : ! write(*,*) '!!!!!!!!!!! remove SOC !!!!!!!!!!!!!!'
178 : ! hsomtx = 0 !!!!!!!!!!!!
179 810 : DO jsp = 1,input%jspins
180 27840 : DO i = 1,nsz(jsp)
181 27030 : hsomtx(i,i,jsp,jsp) = hsomtx(i,i,jsp,jsp) + CMPLX(eig(i,jsp),0.)
182 27564 : IF (input%jspins.EQ.1) THEN
183 566 : hsomtx(i,i,2,2) = hsomtx(i,i,2,2) + CMPLX(eig(i,jsp),0.)
184 : ENDIF
185 : ENDDO
186 : ENDDO
187 :
188 : !
189 : ! resort H-matrix
190 : !
191 1380 : ALLOCATE (hso(2*input%neig,2*input%neig))
192 828 : DO jsp = 1,2
193 1932 : DO jsp1 = 1,2
194 1104 : IF (jsp.EQ.1) nn = 0
195 1104 : IF (jsp1.EQ.1) nn1 = 0
196 1104 : IF (jsp.EQ.2) nn = nsz(1)
197 1104 : IF (jsp1.EQ.2) nn1 = nsz(1)
198 :
199 : !write(3333,'(2i3,4e15.8)') jsp,jsp1,hsomtx(jsp,jsp1,8,8),hsomtx(jsp,jsp1,32,109)
200 56848 : DO i = 1,nsz(jsp)
201 3029656 : DO j = 1,nsz(jsp1)
202 3028552 : hso(i+nn,j+nn1) = hsomtx(i,j,jsp,jsp1)
203 : ENDDO
204 : ENDDO
205 : ENDDO
206 : ENDDO
207 276 : DEALLOCATE (hsomtx)
208 :
209 : !
210 : ! add Sigma-vxc (QSGW)
211 : !
212 276 : IF(l_qsgw) THEN
213 0 : nbas = lapw%nv(1) + atoms%nlotot
214 0 : WRITE(*,'(A,I3,A,I5,A)') 'Read fleur.qsgw (',nk,',',nbas,')'
215 0 : IF( input%jspins .EQ. 2 ) STOP 'alineso: GW+noco not implemented.'
216 0 : ALLOCATE (sigma_xc(2*nsz(1),2*nsz(1)))
217 0 : ALLOCATE (sigma_xc_apw(nbas,nbas))
218 0 : INQUIRE(667,opened=l_open)
219 0 : IF( .NOT.l_open ) THEN
220 0 : IF( nk.NE.1 ) STOP 'unit 667 not opened but not at 1st k'
221 0 : OPEN(667,file='fleur.qsgw',form='unformatted')
222 0 : ELSE IF( nk.EQ.1) THEN
223 0 : REWIND(667)
224 : ENDIF
225 0 : DO jsp1 = 1,2
226 0 : DO jsp2 = 1,jsp1
227 0 : IF(jsp1.EQ.jsp2) THEN
228 0 : READ(667) ((sigma_xc_apw(i,j),j=1,i),i=1,nbas)
229 0 : DO i = 1,nbas
230 0 : DO j = 1,i-1
231 0 : sigma_xc_apw(j,i) = CONJG(sigma_xc_apw(i,j))
232 : ENDDO
233 : ENDDO
234 : ELSE
235 0 : READ(667) sigma_xc_apw
236 : ENDIF
237 : ! write(*,*) 'lo part set to zero!'
238 : ! sigma_xc_apw(nv+1:,nv+1:) = 0
239 0 : i = nsz(1) * (jsp1-1) + 1 ; i1 = nsz(1) * jsp1
240 0 : j = nsz(1) * (jsp2-1) + 1 ; j1 = nsz(1) * jsp2
241 0 : if (input%l_real) THEN
242 : sigma_xc(i:i1,j:j1) = &
243 0 : MATMUL ( TRANSPOSE(zmat(1)%data_r(:nbas,:)) ,&
244 0 : MATMUL ( sigma_xc_apw, zmat(1)%data_r(:nbas,:) ) )
245 : else
246 : sigma_xc(i:i1,j:j1) = &
247 0 : MATMUL ( CONJG(TRANSPOSE(zmat(1)%data_c(:nbas,:))) ,&
248 0 : MATMUL ( sigma_xc_apw, zmat(1)%data_c(:nbas,:) ) )
249 : endif
250 0 : hso(i:i1,j:j1) = hso(i:i1,j:j1) + CONJG(sigma_xc(i:i1,j:j1))
251 0 : IF(jsp1.NE.jsp2) THEN
252 0 : sigma_xc(j:j1,i:i1) = TRANSPOSE(CONJG(sigma_xc(i:i1,j:j1)))
253 0 : hso(j:j1,i:i1) = hso(j:j1,i:i1)+CONJG(sigma_xc(j:j1,i:i1))
254 : ENDIF
255 : ENDDO
256 : ENDDO
257 0 : DEALLOCATE (sigma_xc_apw)
258 : ENDIF
259 :
260 : !
261 : ! diagonalize the hamiltonian using library-routines
262 : !
263 276 : idim_c = 4*input%neig
264 276 : idim_r = 6*input%neig
265 :
266 276 : CALL timestart("alineso SOC: -diag")
267 :
268 1380 : ALLOCATE (cwork(idim_c),rwork(idim_r))
269 :
270 276 : IF (input%eonly) THEN
271 0 : vectors= 'N'
272 : ELSE
273 276 : vectors= 'V'
274 : ENDIF
275 : CALL zheev(vectors,'U',nsize,hso,2*input%neig,eig_so,&
276 276 : cwork, idim_c, rwork, info)
277 276 : IF (info.NE.0) WRITE (oUnit,FMT=8000) info
278 : 8000 FORMAT (' AFTER zheev: info=',i4)
279 276 : CALL timestop("alineso SOC: -diag")
280 :
281 276 : DEALLOCATE (cwork,rwork)
282 :
283 276 : IF (input%eonly) THEN
284 0 : IF(l_socvec) CALL juDFT_error("EONLY set. Vectors not calculated.",calledby ="alineso")
285 : ELSE
286 1104 : ALLOCATE (zhelp2(input%neig,2*input%neig))
287 : !
288 : ! proj. back to G - space: old eigenvector 'z' to new one 'Z'
289 : ! +
290 : ! s --- s | z(G,j,s) ... z(ig,i,jsp)
291 : ! Z (G) = > z (G) * C (i,j) | Z(G,j,s) ... zso(ig,j,jsp)
292 : ! j --- i | C(i,j) ... hso(i ,j)
293 : ! i +
294 : ! reorder new e.w. in 2x2 spin space : zhelp(,1),zhelp(,2)
295 : !
296 :
297 828 : DO i1 = 1,2
298 :
299 552 : jsp = i1
300 552 : jsp2= i1
301 552 : IF (input%jspins.EQ.1) jsp = 1
302 552 : IF (input%jspins.EQ.1 .AND..NOT.(input%l_wann.OR.l_socvec)) jsp2=1
303 552 : IF (i1.EQ.1) nn = 0
304 276 : IF (i1.EQ.2) nn = nsz(1)
305 :
306 3029104 : zhelp2(:,:) = 0.0
307 55744 : DO j = 1,nsize
308 3029104 : DO i = 1,nsz(jsp)
309 3028552 : zhelp2(i,j) = CONJG(hso(i+nn,j))
310 : ENDDO
311 : ENDDO ! j
312 :
313 828 : if (input%l_real) THEN
314 : CALL zgemm("N","N",zmat(1)%matsize1,2*input%neig,input%neig,CMPLX(1.0,0.0),CMPLX(zmat(jsp)%data_r(:,:)),&
315 46838 : zmat(1)%matsize1, zhelp2,input%neig,CMPLX(0.0,0.0), zso(1,1,jsp2),zmat(1)%matsize1)
316 : else
317 : CALL zgemm("N","N",zmat(1)%matsize1,2*input%neig,input%neig, CMPLX(1.0,0.0),zmat(jsp)%data_c(:,:),&
318 546 : zmat(1)%matsize1, zhelp2,input%neig,CMPLX(0.0,0.0), zso(1,1,jsp2),zmat(1)%matsize1)
319 : endif
320 :
321 : ENDDO !isp
322 :
323 276 : IF(l_socvec) THEN
324 : !RS: write SOC vectors to SOCVEC
325 0 : WRITE(4649) lapw%nmat,nsize,input%jspins,nsz,2*input%neig,CONJG(hso)
326 : !CF: write qsgw
327 0 : IF(l_qsgw) THEN
328 0 : nn = 2*nsz(1)
329 0 : sigma_xc = MATMUL ( TRANSPOSE(hso(:nn,:nn)) ,&
330 0 : & MATMUL ( sigma_xc , CONJG(hso(:nn,:nn)) ) )
331 0 : WRITE(1014) nn
332 0 : WRITE(1014) ((sigma_xc(i,j),i=1,j),j=1,nn)
333 0 : DEALLOCATE (sigma_xc)
334 : ENDIF
335 : ENDIF
336 :
337 276 : DEALLOCATE (zhelp2)
338 : ENDIF ! (.NOT.input%eonly)
339 :
340 276 : DEALLOCATE ( hso )
341 : ENDIF ! (n_rank==0)
342 : !
343 548 : nmat=lapw%nmat
344 548 : RETURN
345 :
346 4830 : END SUBROUTINE alineso
347 0 : END MODULE m_alineso
|
