Line data Source code
1 : !!****m* ABINIT/m_berrytk
2 : !! NAME
3 : !! m_berrytk
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2000-2026 ABINIT group (MVeithen)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_berrytk
22 :
23 : use defs_basis
24 : use m_abicore
25 : use m_errors
26 :
27 : use m_cgtools, only : overlap_g
28 : use m_hide_lapack, only : dzgedi, dzgefa
29 :
30 : implicit none
31 :
32 : private
33 : !!***
34 :
35 : public :: smatrix
36 : public :: polcart
37 : !!***
38 :
39 : contains
40 : !!***
41 :
42 : !!****f* ABINIT/smatrix
43 : !! NAME
44 : !! smatrix
45 : !!
46 : !! FUNCTION
47 : !! Compute the overlap matrix between the k-points k and k + dk.
48 : !! Depending on the value of job and ddkflag, compute also its determinant,
49 : !! its inverse and the product of its inverse with the wavefunctions at k.
50 : !!
51 : !! INPUTS
52 : !! cg(2,mcg_k) = planewave coefficients of wavefunctions at k
53 : !! cgq(2,mcg_q) = planewave coefficients of wavefunctions at q = k + dk
54 : !! ddkflag = 1 : compute product of the inverse overlap matrix
55 : !! with the wavefunction at k (job = 1 or 11)
56 : !! 0 : do not compute the product of the inverse overlap matrix
57 : !! with the wavefunction at k
58 : !! icg = shift applied to the wavefunctions in the array cg
59 : !! icg1 = shift applied to the wavefunctions in the array cgq
60 : !! itrs = variable that governs the use of time-reversal symmetry
61 : !! when converting the wavefunctions from the iBZ to the fBZ
62 : !! job = type of calculation
63 : !! 0 : update overlap matrix only
64 : !! 1 : like 0 but also compute inverse of the overlap matrix
65 : !! 10 : like 0 but also compute determinant of the overlap matrix
66 : !! 11 : like 0 but also compute determinant and inverse of the overlap matrix
67 : !! 20 : like 0 but also transfer cgq to cg1_k without multiplication by S^{-1}
68 : !! 21 : like 1 but also transfer cgq to cg1_k without multiplication by S^{-1}
69 : !! maxbd = used in case ddkflag = 1, defines the highest band for
70 : !! which the ddk will be computed
71 : !! mcg_k = second dimension of cg
72 : !! mcg_q = second dimension of cg_q
73 : !! mcg1_k = second dimension of cg1_k, should be equal to
74 : !! mpw*nsppol*nspinor*(maxbd - minbd + 1)
75 : !! minbd = used in case ddkflag = 1, defines the lowest band for
76 : !! which the ddk will be computed
77 : !! mpw = maximum dimensioned size of npw
78 : !! mband_occ = max number of occupied valence bands for both spins
79 : !! nband_occ = number of (occupied) valence bands
80 : !! npw_k1 = number of plane waves at k
81 : !! npw_k2 = number of plane waves at k + dk
82 : !! nspinor = number of spinorial components of the wavefunctions
83 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
84 : !! pwind_k = array used to compute the overlap matrix
85 : !! pwnsfac = phase factors for non-symmorphic translations
86 : !! shifrbd = shift applied to the location of the WF in the cg-array
87 : !! after each loop over bands
88 : !! 0 : apply no shift, this is allowed in case cg
89 : !! contains only the wf of one single band
90 : !! 1 : apply a shift of npw_k1*nspinor, this is the usual option
91 : !! when cg contains the wf for all occupied bands
92 : !! smat_k_paw : overlap matrix due to on-site PAW terms between different bands
93 : !! at k and k+b. Only relevant when usepaw = 1, and is to be computed
94 : !! previously by smatrix_k_paw.F90
95 : !! usepaw = flag governing use of PAW: 0 means no PAW, 1 means PAW is used
96 : !!
97 : !! OUTPUT
98 : !! cg1_k(2,mcg1_k) = product of the inverse overlap matrix with the
99 : !! wavefunctions at k; computed in case job = 1 or 11;
100 : !! or just cgq in case of job = 20 or 21
101 : !! dtm_k(2) = determinant of the overlap matrix between k and k + dk;
102 : !! computed in case job = 10 or 11
103 : !! smat_inv = inverse of the overlap matrix
104 : !!
105 : !! SIDE EFFECTS
106 : !! Input/Output
107 : !! sflag_k(iband) = 1 if the elements smat_k(:,iband,:) are up to date
108 : !! -> they will not be recomputed
109 : !! 0 the elements smat_k(:,iband,:) will be recomputed
110 : !! at the end of the routine, sflag_k(1:mband_occ) = 1
111 : !! (the whole overlap matrix is up to date)
112 : !! smat_k = overlap matrix between k, k + dk
113 : !! only the lines for which sflag_k = 0 are computed
114 : !! smat_k(:,n,m) = < u_{n,k} | u_{m,k+dk} >
115 : !!
116 : !! NOTES
117 : !! This routine is quite flexible in the way it deals with the wavefunctions:
118 : !! - cg (WF at k) can contain either the whole WF array (all k-points
119 : !! and bands), in which case the location of the WF at k is specified
120 : !! by icg and shiftbd = 1, or the WF of a single k-point/band, in which case
121 : !! shiftbd = 0 and icg = 0.
122 : !! - cgq (WF at k + dk) can contain either the whole WF array (all k-points
123 : !! and bands), in which case the location of the WF at k is specified
124 : !! by icg1, or the WF of a single k-point, in which case
125 : !! icg1 = 0. cgq must contain the WF of ALL occupied bands.
126 : !! - cg1_k can either be computed for all valence bands or
127 : !! for a group of valence bands defined by minbd and maxbd.
128 : !!
129 : !! SOURCE
130 :
131 2369386 : subroutine smatrix(cg,cgq,cg1_k,ddkflag,dtm_k,icg,icg1,itrs,job,maxbd,&
132 : & mcg_k,mcg_q,mcg1_k,minbd,mpw,mband_occ,nband_occ,npw_k1,npw_k2,nspinor,&
133 2369386 : & pwind_k,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_k,smat_k_paw,usepaw)
134 :
135 : !Arguments ------------------------------------
136 : !scalars
137 : integer,intent(in) :: ddkflag,icg,icg1,itrs,job,maxbd,mcg1_k,mcg_k,mcg_q
138 : integer,intent(in) :: minbd,mpw,mband_occ,npw_k1,npw_k2,nspinor,shiftbd
139 : integer,intent(in) :: nband_occ
140 : integer,intent(in) :: usepaw
141 : !arrays
142 : integer,intent(in) :: pwind_k(mpw)
143 : integer,intent(inout) :: sflag_k(mband_occ)
144 : real(dp),intent(in) :: cg(2,mcg_k),cgq(2,mcg_q),pwnsfac_k(4,mpw)
145 : real(dp),intent(in) :: smat_k_paw(2,usepaw*mband_occ,usepaw*mband_occ)
146 : real(dp),intent(inout) :: smat_k(2,mband_occ,mband_occ)
147 : real(dp),intent(out) :: cg1_k(2,mcg1_k),dtm_k(2)
148 : real(dp),intent(out) :: smat_inv(2,mband_occ,mband_occ)
149 :
150 : !Local variables -------------------------
151 : !scalars
152 : integer :: count,dzgedi_job,iband,info,ipw,ispinor,jband,jband1,jpw,pwmax,pwmin
153 : integer :: spnshft_k1, spnshft_k2
154 : real(dp) :: doti,dotr,fac,wfi,wfr
155 : character(len=500) :: message
156 : !arrays
157 2369386 : integer,allocatable :: ipvt(:)
158 : ! integer,allocatable :: my_pwind_k(:) ! used in debugging below
159 : real(dp) :: det(2,2)
160 2369386 : real(dp),allocatable :: vect1(:,:),vect2(:,:),zgwork(:,:)
161 :
162 : ! ***********************************************************************
163 :
164 : !DEBUG
165 : !write(std_out,*)'smatrix : enter'
166 : !write(std_out,*)'sflag_k = ',sflag_k
167 : !write(std_out,'(a,4i4)' )'job, ddkflag, shiftbd, itrs = ',job,ddkflag,shiftbd,itrs
168 : !write(std_out,'(a,2i6)')' JWZ smatrix.F90 debug : npw_k1, npw_k2 ',npw_k1,npw_k2
169 : !stop
170 : !ENDDEBUG
171 :
172 7108158 : ABI_MALLOC(ipvt,(nband_occ))
173 7108158 : ABI_MALLOC(zgwork,(2,nband_occ))
174 7108158 : ABI_MALLOC(vect1,(2,0:mpw*nspinor))
175 4738772 : ABI_MALLOC(vect2,(2,0:mpw*nspinor))
176 11846930 : vect1(:,0) = zero ; vect2(:,0) = zero
177 :
178 : !Check if the values of ddkflag and job are compatible
179 :
180 2369386 : if ((job /= 0).and.(job /= 1).and.(job /= 10).and.(job /= 11).and.(job/=20).and.(job/=21)) then
181 : write(message,'(a,i3,a,a)')&
182 0 : & ' job is equal to ',job,ch10,&
183 0 : & ' while only the values job = 0, 1, 10, 11, 20, or 21 are allowed.'
184 0 : ABI_ERROR(message)
185 : end if
186 :
187 2369386 : if (ddkflag == 1) then
188 1055328 : if ((job/=1).and.(job/=11)) then
189 : write(message,'(a,i0,a,a)')&
190 0 : & ' job is equal to ',job,ch10,&
191 0 : & ' while ddkflag = 1. This is not allowed.'
192 0 : ABI_ERROR(message)
193 : end if
194 : end if
195 :
196 : !Check the values of sflag_k
197 31767322 : do iband=1,nband_occ
198 31767322 : if (sflag_k(iband)/=0 .and. sflag_k(iband)/=1)then
199 : write(message,'(3a,i4,a,i4)')&
200 0 : & ' The content of sflag_k must be 0 or 1.',ch10,&
201 0 : & ' However, for iband=',iband,', sflag_k(iband)=',sflag_k(iband)
202 0 : ABI_ERROR(message)
203 : end if
204 : end do
205 :
206 : !Check if shiftbd is consistent with sflag_k
207 2369386 : if (shiftbd == 0) then
208 1049248 : count = 0
209 15128228 : do iband = 1, nband_occ
210 15128228 : if (sflag_k(iband) == 0) count = count + 1
211 : end do
212 1049248 : if (count > 1) then
213 0 : message = 'in case shiftbd = 0, only 1 element of sflag can be 0'
214 0 : ABI_ERROR(message)
215 : end if
216 : end if
217 :
218 : !Update the lines of the overlap matrix for which sflag = 0
219 : !MVeithen: because of sflag, it is more efficient to perform
220 : !the loop over jband inside the loop over iband
221 :
222 : !DEBUG
223 : !write(std_out,*)' smatrix : smat_k(1,1,1)=',smat_k(1,1,1)
224 : !write(std_out,*)' smatrix : sflag_k=',sflag_k
225 : !ENDDEBUG
226 :
227 : !!
228 : !! debugging based on norm of k1 vector
229 : !
230 : !ABI_MALLOC(my_pwind_k,(mpw))
231 : !!
232 : !do iband = 1, nband_occ
233 : !!
234 : !pwmin = (iband-1)*npw_k1*nspinor*shiftbd
235 : !pwmax = pwmin + npw_k1*nspinor
236 : !!
237 : !!! Multiply the bra wave function by the phase factor
238 : !if (itrs==1.or.itrs==11) then ! take complex conjugate of bra
239 : !do ispinor = 1, nspinor
240 : !spnshft_k1=(ispinor-1)*npw_k1
241 : !do ipw = 1,npw_k1
242 : !vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
243 : !& +cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw)
244 : !vect1(2,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw) &
245 : !& -cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
246 : !end do
247 : !end do
248 : !else
249 : !do ispinor=1, nspinor
250 : !spnshft_k1=(ispinor-1)*npw_k1
251 : !do ipw = 1,npw_k1
252 : !vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
253 : !& -cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw)
254 : !vect1(2,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw) &
255 : !& +cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
256 : !end do
257 : !end do
258 : !end if
259 : !
260 : !!
261 : !if (npw_k1*nspinor < mpw*nspinor) vect1(:,npw_k1*nspinor+1:mpw*nspinor) = zero
262 : !
263 : !do jband = 1, nband_occ
264 : !
265 : !pwmin = (jband-1)*npw_k1*nspinor
266 : !pwmax = pwmin + npw_k1*nspinor
267 : !
268 : !if (itrs==10.or.itrs==11) then ! take complex conjugate of ket
269 : !do ispinor=1, nspinor
270 : !spnshft_k2=(ispinor-1)*npw_k1
271 : !do ipw = 1, npw_k1
272 : !my_pwind_k(ipw) = ipw
273 : !vect2(1,spnshft_k2+ipw) = cg(1,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(1,ipw) &
274 : !& +cg(2,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(2,ipw)
275 : !vect2(2,spnshft_k2+ipw) = cg(1,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(2,ipw) &
276 : !& -cg(2,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(1,ipw)
277 : !end do
278 : !end do
279 : !else
280 : !do ispinor=1, nspinor
281 : !spnshft_k2=(ispinor-1)*npw_k1
282 : !do ipw = 1, npw_k1
283 : !my_pwind_k(ipw) = ipw
284 : !vect2(1,spnshft_k2+ipw) = cg(1,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(1,ipw) &
285 : !& -cg(2,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(2,ipw)
286 : !vect2(2,spnshft_k2+ipw) = cg(1,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(2,ipw) &
287 : !& +cg(2,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(1,ipw)
288 : !end do
289 : !end do
290 : !end if
291 : !!
292 : !if (npw_k1*nspinor < mpw*nspinor) vect2(:,npw_k1*nspinor+1:mpw*nspinor) = zero
293 : !!
294 : !call overlap_g(doti,dotr,mpw,npw_k1,npw_k1,nspinor,my_pwind_k,vect1,vect2)
295 : !!
296 : !smat_k(1,iband,jband) = dotr
297 : !smat_k(2,iband,jband) = doti
298 : !!
299 : !if (usepaw == 1) then
300 : !smat_k(1,iband,jband) = smat_k(1,iband,jband)+smat_k_paw(1,iband,jband)
301 : !smat_k(2,iband,jband) = smat_k(2,iband,jband)+smat_k_paw(2,iband,jband)
302 : !end if
303 : !!
304 : !if( (smat_k(1,iband,jband)**2+smat_k(2,iband,jband)**2)>tol8) then
305 : !write(std_out,'(a,i2,a,i2,a,es16.8,a,es16.8)')' JWZ Debug: <',iband,'|',jband,'> = ',&
306 : !& smat_k(1,iband,jband),' + i ',smat_k(2,iband,jband)
307 : !end if
308 : !!
309 : !end do ! jband
310 : !!
311 : !end do ! iband
312 : !!
313 : !ABI_FREE(my_pwind_k)
314 : !!
315 31767322 : do iband = 1, nband_occ
316 :
317 31767322 : if (sflag_k(iband) == 0) then
318 :
319 6796856 : pwmin = (iband-1)*npw_k1*nspinor*shiftbd
320 6796856 : pwmax = pwmin + npw_k1*nspinor
321 : !
322 : ! old version (*** multiply by nspinor missing??? ***)
323 : ! vect1(:,1:npw_k1) = cg(:,icg + 1 + pwmin:icg + pwmax)
324 : !
325 :
326 : ! Multiply the bra wave function by the phase factor
327 6796856 : if (itrs==1.or.itrs==11) then ! take complex conjugate of bra
328 1166608 : do ispinor = 1, nspinor
329 583304 : spnshft_k1=(ispinor-1)*npw_k1
330 70407968 : do ipw = 1,npw_k1
331 : vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
332 69241360 : & +cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw)
333 : vect1(2,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw) &
334 69824664 : & -cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
335 : end do
336 : end do
337 : else
338 12720604 : do ispinor=1, nspinor
339 6507052 : spnshft_k1=(ispinor-1)*npw_k1
340 960082526 : do ipw = 1,npw_k1
341 : vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
342 947361922 : & -cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw)
343 : vect1(2,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw) &
344 953868974 : & +cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
345 : end do
346 : end do
347 : end if
348 :
349 : !
350 112427678 : if (npw_k1*nspinor < mpw*nspinor) vect1(:,npw_k1*nspinor+1:mpw*nspinor) = zero
351 :
352 143501904 : do jband = 1, nband_occ
353 :
354 136705048 : pwmin = (jband-1)*npw_k2*nspinor
355 136705048 : pwmax = pwmin + npw_k2*nspinor
356 :
357 136705048 : if (itrs==10.or.itrs==11) then ! take complex conjugate of ket
358 74295704 : do ispinor=1, nspinor
359 37147852 : spnshft_k2=(ispinor-1)*npw_k2
360 6185481456 : do ipw = 1, npw_k2
361 : vect2(1,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw) &
362 6111185752 : & +cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(4,ipw)
363 : vect2(2,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(4,ipw) &
364 6148333604 : & -cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw)
365 : end do
366 : end do
367 : else
368 201462392 : do ispinor=1, nspinor
369 101905196 : spnshft_k2=(ispinor-1)*npw_k2
370 16955831628 : do ipw = 1, npw_k2
371 : vect2(1,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw) &
372 16754369236 : & -cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(4,ipw)
373 : vect2(2,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(4,ipw) &
374 16856274432 : & +cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw)
375 : end do
376 : end do
377 : end if
378 :
379 2207792428 : if (npw_k2*nspinor < mpw*nspinor) vect2(:,npw_k2*nspinor+1:mpw*nspinor) = zero
380 :
381 : ! DEBUG
382 : ! if(iband==1 .and. jband==1 .and. itrs==0 .and. npw_k1==68 .and. npw_k2==74)then
383 : ! write(std_out,'(a)' )' smatrix : ii,vect1,cg,pwnsfac='
384 : ! do ii=1,npw_k1
385 : ! write(std_out,'(i4,6es16.6)' )ii,vect1(:,ii),cg(:,icg+ii+pwmin),pwnsfac_k(1:2,ii)
386 : ! end do
387 : ! do ii=1,npw_k2
388 : ! write(std_out,'(i4,6es16.6)' )ii,vect2(:,ii),cg(:,icg1+ii+pwmin),pwnsfac_k(3:4,ii)
389 : ! end do
390 : ! end if
391 : ! ENDDEBUG
392 :
393 136705048 : call overlap_g(doti,dotr,mpw,npw_k1,npw_k2,nspinor,pwind_k,vect1,vect2)
394 :
395 136705048 : smat_k(1,iband,jband) = dotr
396 136705048 : smat_k(2,iband,jband) = doti
397 :
398 143501904 : if (usepaw == 1) then
399 8172096 : smat_k(1,iband,jband) = smat_k(1,iband,jband)+smat_k_paw(1,iband,jband)
400 8172096 : smat_k(2,iband,jband) = smat_k(2,iband,jband)+smat_k_paw(2,iband,jband)
401 : end if
402 :
403 : ! DEBUG
404 : ! if(iband==1 .and. jband==1)then
405 : ! write(std_out,'(a,2es16.6,3i4)' )' smatrix : dotr,smat_k(1,iband,jband),mpw,npw_k1,npw_k2',dotr,smat_k(1,iband,jband),mpw,npw_k1,npw_k2
406 : ! end if
407 : ! ENDDEBUG
408 :
409 : end do ! jband
410 :
411 : end if ! sflag_k(iband) == 0
412 :
413 : end do ! iband
414 :
415 : !DEBUG
416 : !do iband=1,nband_occ
417 : !do jband=1,nband_occ
418 : !write(std_out,'(a,2i4,2e20.10)') 'smat',iband,jband,smat_k(1,iband,jband),smat_k(2,iband,jband)
419 : !end do
420 : !end do
421 : !write(std_out,*)' smatrix : smat_k(1,1,1)=',smat_k(1,1,1)
422 : !ENDDEBUG
423 :
424 : !Update sflag_k
425 31767322 : sflag_k(:) = 1
426 :
427 : !Depending on the value of job, compute the determinant of the
428 : !overlap matrix, its inverse or the product of the inverse
429 : !overlap matrix with the WF at k.
430 :
431 2369386 : if ((job==1).or.(job==10).or.(job==11).or.(job==21)) then
432 :
433 2373532522 : smat_inv(:,:,:) = smat_k(:,:,:)
434 :
435 : ! DEBUG
436 : ! write(std_out,*)' smatrix : smat_inv=',smat_inv
437 : ! ENDDEBUG
438 :
439 2369386 : dzgedi_job=job; if(job==21) dzgedi_job=1
440 : ! TODO: should this be over nband_occ(isppol)?
441 2369386 : call dzgefa(smat_inv,mband_occ,nband_occ,ipvt,info)
442 2369386 : call dzgedi(smat_inv,mband_occ,nband_occ,ipvt,det,zgwork,dzgedi_job)
443 :
444 : ! DEBUG
445 : ! write(std_out,*)' smatrix : det=',det
446 : ! ENDDEBUG
447 :
448 : ! Compute the determinant of the overlap matrix
449 2369386 : dtm_k(:) = zero
450 2369386 : if (job==10 .or. job==11) then
451 2363402 : fac = exp(log(10._dp)*det(1,2))
452 : dtm_k(1) = fac*(det(1,1)*cos(log(10._dp)*det(2,2)) - &
453 2363402 : & det(2,1)*sin(log(10._dp)*det(2,2)))
454 : dtm_k(2) = fac*(det(1,1)*sin(log(10._dp)*det(2,2)) + &
455 2363402 : & det(2,1)*cos(log(10._dp)*det(2,2)))
456 : end if
457 :
458 : ! Compute the product of the inverse overlap matrix with the WF
459 :
460 2369386 : if (ddkflag == 1) then
461 :
462 446370900 : cg1_k(:,:) = zero
463 1055328 : jband1 = 0
464 :
465 1055328 : if (itrs == 10 .or. itrs == 11) then
466 :
467 377368 : do jband = minbd, maxbd
468 188684 : jband1 = jband1 + 1
469 3696414 : do iband = 1, nband_occ
470 :
471 6826776 : do ispinor = 1, nspinor
472 3319046 : spnshft_k1 = (ispinor-1)*npw_k1
473 3319046 : spnshft_k2 = (ispinor-1)*npw_k2
474 544516846 : do ipw = 1, npw_k1
475 :
476 537878754 : jpw = pwind_k(ipw)
477 :
478 541197800 : if (jpw > 0) then
479 :
480 : wfr = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
481 503615708 : & -cgq(2,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)
482 : wfi = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)&
483 503615708 : & +cgq(2,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
484 :
485 : cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = &
486 : & cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) + &
487 503615708 : & smat_inv(1,iband,jband)*wfr + smat_inv(2,iband,jband)*wfi
488 :
489 : cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = &
490 : & cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) - &
491 503615708 : & smat_inv(1,iband,jband)*wfi + smat_inv(2,iband,jband)*wfr
492 :
493 : end if
494 :
495 : end do ! end loop over npw_k1
496 : end do ! end loop over nspinor
497 :
498 : end do
499 : end do
500 :
501 : else
502 :
503 1752040 : do jband = minbd, maxbd
504 885396 : jband1 = jband1 + 1
505 12615398 : do iband = 1, nband_occ
506 :
507 23143808 : do ispinor = 1, nspinor
508 11395054 : spnshft_k1 = (ispinor-1)*npw_k1
509 11395054 : spnshft_k2 = (ispinor-1)*npw_k2
510 1799568542 : do ipw = 1, npw_k1
511 :
512 1777310130 : jpw = pwind_k(ipw)
513 :
514 1788705184 : if (jpw > 0) then
515 :
516 : wfr = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
517 1666755170 : & -cgq(2,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)
518 : wfi = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)&
519 1666755170 : & +cgq(2,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
520 :
521 : cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = &
522 : & cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) + &
523 1666755170 : & smat_inv(1,iband,jband)*wfr - smat_inv(2,iband,jband)*wfi
524 :
525 : cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = &
526 : & cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) + &
527 1666755170 : & smat_inv(1,iband,jband)*wfi + smat_inv(2,iband,jband)*wfr
528 :
529 : end if
530 :
531 : end do ! end loop over npw_k1
532 : end do ! end loop over nspinor
533 :
534 : end do
535 : end do
536 :
537 : end if ! itrs
538 :
539 : end if
540 :
541 : end if !
542 :
543 2369386 : if(job == 20 .or. job == 21) then ! special case transfering cgq to cg1_k without use of S^{-1}, used in
544 : ! magnetic field case
545 :
546 0 : cg1_k(:,:) = zero
547 0 : jband1 = 0
548 :
549 0 : if (itrs == 10 .or. itrs == 11) then
550 :
551 0 : do jband = minbd, maxbd
552 0 : jband1 = jband1 + 1
553 0 : do ispinor = 1, nspinor
554 0 : spnshft_k1 = (ispinor-1)*npw_k1
555 0 : spnshft_k2 = (ispinor-1)*npw_k2
556 0 : do ipw = 1, npw_k1
557 0 : jpw = pwind_k(ipw)
558 :
559 0 : if (jpw > 0) then
560 : wfr = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
561 0 : & -cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)
562 : wfi = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)&
563 0 : & +cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
564 :
565 0 : cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfr
566 0 : cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfi
567 :
568 : end if
569 :
570 : end do ! end loop over npw_k1
571 : end do ! end loop over nspinor
572 :
573 : end do
574 :
575 : else
576 :
577 0 : do jband = minbd, maxbd
578 0 : jband1 = jband1 + 1
579 0 : do ispinor = 1, nspinor
580 0 : spnshft_k1 = (ispinor-1)*npw_k1
581 0 : spnshft_k2 = (ispinor-1)*npw_k2
582 0 : do ipw = 1, npw_k1
583 0 : jpw = pwind_k(ipw)
584 :
585 0 : if (jpw > 0) then
586 :
587 : wfr = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
588 0 : & -cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)
589 : wfi = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)&
590 0 : +cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
591 :
592 0 : cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfr
593 0 : cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfi
594 :
595 : end if
596 :
597 : end do ! end loop over npw_k1
598 : end do ! end loop over nspinor
599 :
600 : end do
601 :
602 : end if ! itrs
603 :
604 : end if ! end job == 20 .or. job == 21 case
605 :
606 2369386 : ABI_FREE(ipvt)
607 2369386 : ABI_FREE(zgwork)
608 2369386 : ABI_FREE(vect1)
609 2369386 : ABI_FREE(vect2)
610 :
611 : !DEBUG
612 : !write(std_out,*)' dtm_k=',dtm_k(:)
613 : !write(std_out,*)' smatrix : exit '
614 : !ENDDEBUG
615 :
616 2369386 : end subroutine smatrix
617 : !!***
618 :
619 : !!****f* ABINIT/polcart
620 : !! NAME
621 : !! polcart
622 : !!
623 : !! FUNCTION
624 : !! Transform polarization from reduced to cartesian coordinates,
625 : !! divide by ucvol and write the result to an output file
626 : !!
627 : !! INPUTS
628 : !! pel(3) = reduced coordinates of the electronic polarization
629 : !! pion(3) = reduced coordinates of the ionic polarization
630 : !! polunit = units used for the output of the polarization
631 : !! 1 : use atomic units
632 : !! 2 : use MKS units
633 : !! 3 : use both atomic and MKS units
634 : !! rprimd(3,3) = dimensional primitive translations (bohr)
635 : !! ucvol = volume of the primitive unit cell
636 : !! unit_out = unit for output of the results
637 : !! usepaw = 1 for PAW computation, zero else
638 : !!
639 : !! OUTPUT
640 : !! pel_cart(3) = cartesian coords of the electronic polarization
641 : !! in atomic units
642 : !! pelev(3)= expectation value polarization term (PAW only) in cartesian coordinates
643 : !! pion_cart(3)= cartesian coords of the ionic polarization
644 : !! in atomic units
645 : !! ptot_cart(3)= cartesian coords of the total polarization
646 : !! in atomic units
647 : !!
648 : !! NOTES
649 : !! - The sum of the electronic and ionic Berry phase is folded into
650 : !! [-1,1] before it is transformed to cartesian coordinates.
651 : !! This means that in some cases, ptot_cart /= pel_cart + pion_cart
652 : !!
653 : !! - pel and pion do not take into account the factor 1/ucvol.
654 : !! At the opposite, this factor is taken into account in
655 : !! pel_cart and pion_cart
656 : !! - unit_out = 0 is allowed, in this case, there will be no
657 : !! output of the results
658 : !!
659 : !! SOURCE
660 :
661 270 : subroutine polcart(red_ptot,pel,pel_cart,pelev,pion,pion_cart,polunit,&
662 : & ptot_cart,rprimd,ucvol,unit_out,usepaw)
663 :
664 : !Arguments ------------------------------------
665 : !scalars
666 : integer,intent(in) :: polunit,unit_out,usepaw
667 : real(dp),intent(in) :: ucvol
668 : !arrays
669 : real(dp),intent(in) :: red_ptot(3) !!REC
670 : real(dp),intent(in) :: pel(3),pelev(3),pion(3),rprimd(3,3)
671 : real(dp),intent(out) :: pel_cart(3),pion_cart(3),ptot_cart(3)
672 :
673 : !Local variables -------------------------
674 : !scalars
675 : integer :: idir
676 : character(len=500) :: message
677 : !arrays
678 : real(dp) :: pel_mks(3),pelev_mks(3),pion_mks(3),ptot(3),ptot_mks(3)
679 :
680 : ! GA: - The input/ouput description above is wrong.
681 : ! - Why do we care whether we used PAW?
682 :
683 : ! ***********************************************************************
684 : !!REC Note ptot has already been folded and kept onto same branch
685 : !unless ptot=0d0, in which case ptot has not been computed yet
686 1080 : if( sum(abs(red_ptot(:))) < tol8 )then
687 24 : ptot(:) = pel(:) + pion(:)
688 : ! Fold ptot into [-1, 1]
689 24 : do idir = 1, 3
690 24 : ptot(idir) = ptot(idir) - 2_dp*nint(ptot(idir)/2_dp)
691 : end do
692 : else !!REC
693 264 : ptot=red_ptot !!REC
694 : end if !!REC
695 :
696 : !Transform pel, pion and ptot to cartesian coordinates
697 270 : pel_cart(:) = zero ; pion_cart(:) = zero ; ptot_cart(:) = zero
698 1080 : do idir = 1, 3
699 : pel_cart(idir) = rprimd(idir,1)*pel(1) + rprimd(idir,2)*pel(2) + &
700 810 : & rprimd(idir,3)*pel(3)
701 : pion_cart(idir) = rprimd(idir,1)*pion(1) + rprimd(idir,2)*pion(2) + &
702 810 : & rprimd(idir,3)*pion(3)
703 : ptot_cart(idir) = rprimd(idir,1)*ptot(1) + rprimd(idir,2)*ptot(2) + &
704 1080 : & rprimd(idir,3)*ptot(3)
705 : end do
706 :
707 : !Divide by the unit cell volume
708 1080 : pel_cart(:) = pel_cart(:)/ucvol
709 1080 : pion_cart(:) = pion_cart(:)/ucvol
710 1080 : ptot_cart(:) = ptot_cart(:)/ucvol
711 : !pelev is either zero (in NCPP case), or possibly non-zero (in PAW case)
712 : !however, in the PAW case, it is already implicitly included in the computation
713 : !of pel and should not be added in here. Below in the PAW case we report its
714 : !value to the user, which is helpful for comparison to the USPP theory.
715 : !13 June 2012 J Zwanziger
716 : !note that pelev was AlREADY in cartesian frame.
717 : !ptot_cart(:) = (ptot_cart(:)+pelev(:))/ucvol
718 :
719 : !Write the results to unit_out (if /= 0)
720 : !Use the coordinates specified by the value polunit
721 :
722 : !Atomic units
723 :
724 270 : if (((polunit == 1).or.(polunit == 3)).and.(unit_out /= 0)) then
725 :
726 : ! write(message,'(7(a),3(e16.9,2x),a,a,3(e16.9,2x),a,a,3(e16.9,2x),a,a,3(e16.9,2x))')ch10,&
727 : ! & ' Polarization in cartesian coordinates (a.u.):',ch10,&
728 : ! & ' (the sum of the electronic and ionic Berry phase',&
729 : ! & ' has been fold into [-1, 1])',ch10,&
730 : ! & ' Electronic berry phase: ', (pel_cart(idir), idir = 1, 3), ch10,&
731 : ! & ' Expectation value (PAW only): ', (pelev(idir)/ucvol, idir = 1, 3), ch10,&
732 : ! & ' Ionic: ', (pion_cart(idir), idir = 1, 3), ch10, &
733 : ! & ' Total: ', (ptot_cart(idir), idir = 1, 3)
734 : ! call wrtout(unit_out,message,'COLL')
735 79 : write(message,'(7(a),3(e16.9,2x),a,a,3(e16.9,2x))')ch10,&
736 79 : & ' Polarization in cartesian coordinates (a.u.):',ch10,&
737 79 : & ' (the sum of the electronic and ionic Berry phase',&
738 79 : & ' has been folded into [-1, 1])',ch10,&
739 158 : & ' Electronic berry phase: ', (pel_cart(idir), idir = 1, 3)
740 79 : call wrtout(unit_out,message,'COLL')
741 79 : if(usepaw==1) then
742 : write(message,'(a,3(e16.9,2x))')&
743 120 : & ' ...includes PAW on-site term: ', (pelev(idir)/ucvol, idir = 1, 3)
744 30 : call wrtout(unit_out,message,'COLL')
745 : end if
746 : write(message,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')&
747 79 : & ' Ionic: ', (pion_cart(idir), idir = 1, 3), ch10, &
748 158 : & ' Total: ', (ptot_cart(idir), idir = 1, 3)
749 79 : call wrtout(unit_out,message,'COLL')
750 :
751 : end if
752 :
753 : !MKS units
754 :
755 270 : if (((polunit == 2).or.(polunit == 3)).and.(unit_out /= 0)) then
756 :
757 316 : pel_mks(:) = pel_cart(:)*(e_Cb)/(Bohr_Ang*1d-10)**2
758 316 : pion_mks(:) = pion_cart(:)*(e_Cb)/(Bohr_Ang*1d-10)**2
759 316 : pelev_mks(:) = pelev(:)/ucvol*(e_Cb)/(Bohr_Ang*1d-10)**2
760 316 : ptot_mks(:) = (ptot_cart(:))*(e_Cb)/(Bohr_Ang*1d-10)**2
761 :
762 79 : write(message,'(7(a),3(e16.9,2x),a,a,3(e16.9,2x),a,a,3(e16.9,2x),a,a,3(e16.9,2x))')ch10,&
763 79 : & ' Polarization in cartesian coordinates (C/m^2):',ch10,&
764 79 : & ' (the sum of the electronic and ionic Berry phase',&
765 79 : & ' has been folded into [-1, 1])',ch10,&
766 158 : & ' Electronic berry phase: ', (pel_mks(idir), idir = 1, 3)
767 79 : call wrtout(unit_out,message,'COLL')
768 79 : if(usepaw==1) then
769 : write(message,'(a,3(e16.9,2x))')&
770 30 : & ' ...includes PAW on-site term: ', (pelev_mks(idir), idir = 1, 3)
771 30 : call wrtout(unit_out,message,'COLL')
772 : end if
773 : write(message,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')&
774 79 : & ' Ionic: ', (pion_mks(idir), idir = 1, 3), ch10, &
775 158 : & ' Total: ', (ptot_mks(idir), idir = 1, 3)
776 79 : call wrtout(unit_out,message,'COLL')
777 :
778 : end if
779 :
780 270 : end subroutine polcart
781 : !!***
782 :
783 : end module m_berrytk
784 : !!***
|
