Branch data 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 : : !!***
|
