Branch data Line data Source code
1 : : !!****m* ABINIT/m_berryphase
2 : : !! NAME
3 : : !! m_berryphase
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2000-2026 ABINIT group (NSAI,XG,MKV)
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_berryphase
22 : :
23 : : use defs_basis
24 : : use m_errors
25 : : use m_abicore
26 : : use m_hdr
27 : : use m_dtset
28 : :
29 : : use m_geometry, only : xred2xcart
30 : : use m_hide_lapack, only : dzgedi, dzgefa
31 : : use m_matrix, only : matr3inv
32 : :
33 : : implicit none
34 : :
35 : : private
36 : : !!***
37 : :
38 : : public :: berryphase
39 : : !!***
40 : :
41 : : contains
42 : : !!***
43 : :
44 : : !!****f* ABINIT/berryphase
45 : : !! NAME
46 : : !! berryphase
47 : : !!
48 : : !! FUNCTION
49 : : !! This routine is called in scfcv.f to compute the electronic Berry Phase
50 : : !! polarization and the ionic contribution to the polarization
51 : : !! Work for nsppol=1 or 2 ,but only accept nspinor=1, and mkmem=nkpt
52 : : !! or 0, kptopt = 2 or 3
53 : : !!
54 : : !! INPUTS
55 : : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
56 : : !! bdberry(4)=band limits for Berry phase contributions,
57 : : !! spin up and spin down (bdberry(3:4) is irrelevant when nsppol=1)
58 : : !! cg(2,mcg)=planewave coefficients of wavefunctions
59 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
60 : : !! istwfk(nkpt_)=input option parameter that describes the storage of wfs
61 : : !! kberry(3,nberry)= different delta k for Berry phases, in unit of kptrlatt
62 : : !! only kberry(1:3,1:nberry) is relevant
63 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates
64 : : !! kpt_(3,nkpt_)=reduced coordinates of k points generated by ABINIT,
65 : : !! kpt_ sampels half the BZ if time-reversal symetrie is used
66 : : !! kptopt=2 when time-reversal symetrie is used
67 : : !! =3 when time-reversal symetrie is not used
68 : : !! kptrlatt(3,3)=k-point lattice specification
69 : : !! mband=maximum number of bands
70 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
71 : : !! mkmem=number of k points treated by this node.
72 : : !! mpw=maximum dimensioned size of npw
73 : : !! natom=number of atoms in cell
74 : : !! nattyp(ntypat)= # atoms of each type.
75 : : !! nband(nkpt*nsppol)=number of bands at each k point, for each polarization
76 : : !! nberry=number of Berry phases to be computed
77 : : !! nkpt=number of k points
78 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
79 : : !! nspinor=number of spinorial components (on current proc)
80 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
81 : : !! ntypat=number of types of atoms in unit cell
82 : : !! nkpt_=number of k points generated by ABINIT, (see kpt_)
83 : : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
84 : : !! ucvol=unit cell volume in bohr**3.
85 : : !! xred(3,natom)=reduced atomic coordinates
86 : : !! zion(ntypat)=valence charge of each type of atom
87 : : !!
88 : : !! OUTPUT
89 : : !! (the polarization is printed)
90 : : !!
91 : : !! SIDE EFFECTS
92 : : !!
93 : : !! TODO
94 : : !! Cleaning, checking for rules.
95 : : !! Should allow for time-reversal symmetry (istwfk)
96 : : !! Should use randac to scan rapidly the wf file
97 : : !!
98 : : !! NOTES
99 : : !! Local Variables:
100 : : !! cmatrix(:,:,:)= overlap matrix of size maxband*maxband
101 : : !! cg_index(:,:,:)= unpacked cg index array for specific band,
102 : : !! k point and polarization.
103 : : !! det(2,2)= intermediate output of Lapack routine zgedi.f
104 : : !! determinant(:,:)= determinant of cmatrix
105 : : !! det_average(2)= averaged det_string over all strings
106 : : !! det_string(:,:)= determinant product of cmatrices along each string
107 : : !! dk(3)= step taken to the next k mesh point along the kberry direction
108 : : !! dkptnext(3)= step between the next and current k point
109 : : !! dphase= phase angle computed from rel_string(2)
110 : : !! gpard(3)= dimensionalreciprocal lattice vector G along which the
111 : : !! polarization is computed
112 : : !! kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
113 : : !! planewave and k point
114 : : !! kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
115 : : !! kpt_flag(nkpt)=kpt_flag(ikpt)=0 when the wf was generated by the ABINIT code
116 : : !! kpt_flag(ikpt) gives the indices of the k-point related
117 : : !! to ikpt by time reversal symetrie
118 : : !! kpt_mark(nkpt)= 0, if k point is unmarked; 1, if k point has been marked
119 : : !! maxband/minband= control the minimum and maximum band calculated in the
120 : : !! overlap matrix
121 : : !! nkstr= number of k points per string
122 : : !! npw_k= npwarr(ikpt), number of planewaves in basis at this k point
123 : : !! nstr= number of k point strings
124 : : !! nkpt=number of k points in the whole BZ
125 : : !! phase0= phase angle computed from det_average(2)
126 : : !! polberry(:)= berry phase of each string (2/nsppol)*(phase0+dphase)/two_pi
127 : : !! polb(isppol) = total berry phase polarization for each spin
128 : : !! polbtot= total berry phase polarization
129 : : !! polion= ionic polarization for each ion
130 : : !! politot= total ionic polarization
131 : : !! poltot= total polarization = polbtot + politot
132 : : !! rel_string(2)= det_string(2)/det_average(2)
133 : : !! shift_g(nkpt)= .true. if the k point should be shifted by a G vector;
134 : : !! .false. if not
135 : : !! tr(2)=variable that changes k to -k
136 : : !! G to -G
137 : : !! $c_g$ to $c_g^*$
138 : : !! when time-reversal symetrie is used
139 : : !! xcart(3,natom)= cartesian coordinates of atoms (bohr)
140 : : !! xcart_reindex(:,:,:)= unpack xcart for each atomic species and number
141 : : !! of atoms for each species
142 : : !!
143 : : !! WARNING
144 : : !! This routine is not yet memory optimized
145 : : !! It might be also rather time-consuming, since there is a
146 : : !! double loop on the number of plane waves.
147 : : !!
148 : : !! SOURCE
149 : :
150 : 0 : subroutine berryphase(atindx1,bdberry,cg,gprimd,istwfk,kberry,kg,kpt_,&
151 : : & kptopt,kptrlatt,mband,mcg,&
152 : 0 : & mkmem,mpw,natom,nattyp,nband,nberry,npwarr,nspinor,nsppol,ntypat,&
153 : 0 : & nkpt_,rprimd,ucvol,xred,zion)
154 : :
155 : : !Arguments ------------------------------------
156 : : !scalars
157 : : integer,intent(in) :: kptopt,mband,mcg,mkmem,mpw,natom,nberry,nkpt_
158 : : integer,intent(in) :: nspinor,nsppol,ntypat
159 : : real(dp),intent(in) :: ucvol
160 : : !arrays
161 : : integer,intent(in) :: atindx1(natom),bdberry(4),istwfk(nkpt_),kberry(3,nberry)
162 : : integer,intent(in) :: kg(3,mpw*mkmem),kptrlatt(3,3),nattyp(ntypat)
163 : : integer,intent(in) :: nband(nkpt_*nsppol),npwarr(nkpt_)
164 : : real(dp),intent(in) :: cg(2,mcg),gprimd(1:3,1:3)
165 : : real(dp),intent(in) :: kpt_(3,nkpt_),rprimd(3,3),xred(3,natom),zion(ntypat)
166 : :
167 : : !Local variables -------------------------
168 : : !scalars
169 : : integer :: band_in,cg_index_iband,cg_index_jband,flag1,iatom
170 : : integer :: iattyp,iband,iberry,icg,ii,ikpt,ikpt2,index,index1,info
171 : : integer :: ipw,isppol,istr,itypat,iunmark,jband,jj,jkpt,jkstr
172 : : integer :: jkstr_ori,jpw,lkstr,lkstr_ori,lkstr_ori_,maxband
173 : : integer :: minband,nband_k,nkpt,nkstr,npw_k,nstr,read_k
174 : : real(dp) :: det_mod,dphase,fac,gmod,phase0,pol,polbtot,polion,politot
175 : : real(dp) :: poltot
176 : : character(len=500) :: message
177 : : !arrays
178 : 0 : integer :: dg(3),kpt_flag(2*nkpt_),kpt_mark(2*nkpt_)
179 : 0 : integer,allocatable :: cg_index(:,:,:),ikpt_dk(:),ikstr(:,:),ipvt(:)
180 : 0 : integer,allocatable :: kg_dum(:,:),kg_kpt(:,:,:)
181 : : real(dp) :: det(2,2),det_average(2),diffk(3),dk(3),gpard(3)
182 : 0 : real(dp) :: klattice(3,3),kptrlattr(3,3),polb(nsppol),rel_string(2),tr(2)
183 : 0 : real(dp) :: xcart(3,natom)
184 : 0 : real(dp),allocatable :: cmatrix(:,:,:),det_string(:,:)
185 : 0 : real(dp),allocatable :: det_tmp(:,:),determinant(:,:),kpt(:,:)
186 : 0 : real(dp),allocatable :: polberry(:),xcart_reindex(:,:,:)
187 : 0 : real(dp),allocatable :: zgwork(:,:)
188 : 0 : logical,allocatable :: shift_g(:)
189 : :
190 : : ! ***********************************************************************
191 : :
192 : : !DEBUG
193 : : !write(std_out,*)' berryphase : enter '
194 : : !ENDDEBUG
195 : :
196 [ # # ]: 0 : if(nspinor==2)then
197 : 0 : message = ' berryphase : does not yet work for nspinor=2'
198 : 0 : ABI_ERROR(message)
199 : : end if
200 : :
201 [ # # # # ]: 0 : if(maxval(istwfk(:))/=1)then
202 : : write(message, '(a,a,a)' )&
203 : 0 : & ' Sorry, this routine does not work yet with istwfk/=1.',ch10,&
204 : 0 : & ' This should have been tested previously ...'
205 : 0 : ABI_BUG(message)
206 : : end if
207 : :
208 : : !change8: set up the whole k point grid in the case where kptopt = 2
209 [ # # ]: 0 : if (kptopt==3) then
210 : 0 : nkpt = nkpt_
211 [ # # # # : 0 : ABI_MALLOC(kpt,(3,nkpt))
# # # # #
# # # ]
212 [ # # # # ]: 0 : kpt(:,:)=kpt_(:,:)
213 [ # # ]: 0 : else if (kptopt==2) then
214 : 0 : nkpt = nkpt_*2
215 [ # # # # : 0 : ABI_MALLOC(kpt,(3,nkpt))
# # # # #
# # # ]
216 [ # # ]: 0 : do ikpt = 1,nkpt/2
217 : 0 : kpt_flag(ikpt) = 0
218 [ # # ]: 0 : kpt(:,ikpt)=kpt_(:,ikpt)
219 : : end do
220 : 0 : index = 0
221 [ # # ]: 0 : do ikpt = (nkpt/2+1),nkpt
222 : 0 : flag1 = 0
223 [ # # ]: 0 : do jkpt = 1, nkpt/2
224 : : if (((abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt))<1.0d-8).or.&
225 : : & (abs(1-abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt)))<1.0d-8))&
226 : : & .and.((abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt))<1.0d-8).or.&
227 : : & (abs(1-abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt)))<1.0d-8))&
228 [ # # # # : 0 : & .and.((abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt))<1.0d-8).or.&
# # # # #
# ]
229 [ # # ]: 0 : & (abs(1-abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt)))<1.0d-8))) then
230 : 0 : flag1 = 1
231 : 0 : index = index + 1
232 : : exit
233 : : end if
234 : : end do
235 : 0 : if (flag1==0) then
236 : 0 : kpt_flag(ikpt-index)=ikpt-nkpt/2
237 [ # # ]: 0 : kpt(:,ikpt-index)=-kpt_(:,ikpt-nkpt/2)
238 : : end if
239 : : end do
240 : 0 : nkpt = nkpt - index
241 : : end if
242 : :
243 : : !change8
244 : :
245 [ # # # # : 0 : ABI_MALLOC(shift_g,(nkpt))
# # # # #
# # # ]
246 [ # # ]: 0 : ABI_MALLOC(kg_dum,(3,0))
247 : :
248 : : !Compute primitive vectors of the k point lattice
249 : : !Copy to real(dp)
250 [ # # # # ]: 0 : kptrlattr(:,:)=kptrlatt(:,:)
251 : : !Go to reciprocal space (in reduced coordinates)
252 : 0 : call matr3inv(kptrlattr,klattice)
253 : :
254 [ # # ]: 0 : do iberry=1,nberry
255 : :
256 : : ! Calculate dimensional recip lattice vector along which P is calculated
257 : : ! dk = step to the nearest k point along that direction
258 : : ! in reduced coordinates
259 : : dk(:)=kberry(1,iberry)*klattice(:,1)+&
260 : : & kberry(2,iberry)*klattice(:,2)+&
261 [ # # ]: 0 : & kberry(3,iberry)*klattice(:,3)
262 [ # # ]: 0 : gpard(:)=dk(1)*gprimd(:,1)+dk(2)*gprimd(:,2)+dk(3)*gprimd(:,3)
263 [ # # ]: 0 : gmod=sqrt(dot_product(gpard,gpard))
264 : :
265 : : ! *****************************************************************************
266 : : ! Select the k grid points along the kberry direction
267 : : ! dk = step to the nearest k point along that direction
268 : :
269 : : ! For each k point, find k_prim such that k_prim= k + dk mod(G)
270 : : ! where G is a vector of the reciprocal lattice
271 [ # # # # : 0 : ABI_MALLOC(ikpt_dk,(nkpt))
# # # # ]
272 [ # # ]: 0 : shift_g(:)= .false.
273 [ # # ]: 0 : do ikpt=1,nkpt
274 [ # # ]: 0 : do ikpt2=1,nkpt
275 [ # # ]: 0 : diffk(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)-dk(:))
276 [ # # # # ]: 0 : if(sum(abs(diffk(:)-nint(diffk(:))))<3*tol8)then
277 : 0 : ikpt_dk(ikpt)=ikpt2
278 [ # # # # ]: 0 : if(sum(diffk(:))>=3*tol8)shift_g(ikpt2) = .true.
279 : : exit
280 : : end if
281 : : end do
282 : : end do
283 : :
284 : : ! DEBUG
285 : : ! do ikpt = 1,nkpt
286 : : ! write(100,*)'ikpt_dk = ',ikpt_dk(ikpt)
287 : : ! if (shift_g(ikpt))then
288 : : ! write(100,*)'true'
289 : : ! else
290 : : ! write(100,*)'false'
291 : : ! end if
292 : : ! write(100,*)''
293 : : ! end do
294 : : ! ENDDEBUG
295 : :
296 : : ! Find the string length, starting from k point 1
297 : : ! (all strings must have the same number of points)
298 : 0 : nkstr=1
299 : 0 : ikpt2=1
300 [ # # ]: 0 : do ikpt=1,nkpt
301 : 0 : ikpt2=ikpt_dk(ikpt2)
302 [ # # ]: 0 : if(ikpt2==1)exit
303 : 0 : nkstr=nkstr+1
304 : : end do
305 : :
306 : : ! Check that the string length is a divisor of nkpt
307 [ # # ]: 0 : if(mod(nkpt,nkstr)/=0)then
308 : 0 : write(message,'(a,a,a,a,i5,a,i7)')ch10,&
309 : 0 : & ' berryphase: BUG -',ch10,&
310 : 0 : & ' The string length=',nkstr,', is not a divisor of nkpt=',nkpt
311 : 0 : call wrtout(std_out,message,'COLL')
312 : : end if
313 : 0 : nstr=nkpt/nkstr
314 : :
315 : 0 : write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
316 : 0 : & ' Computing the polarization (Berry phase) for reciprocal vector:',ch10,&
317 : 0 : & dk(:),' (in reduced coordinates)',ch10,&
318 : 0 : & gpard(1:3),' (in cartesian coordinates - atomic units)'
319 : 0 : call wrtout(ab_out,message,'COLL')
320 : 0 : call wrtout(std_out,message,'COLL')
321 : :
322 : : write(message,'(a,i5,a,a,i5)')&
323 : 0 : & ' Number of strings: ',nstr,ch10,&
324 : 0 : & ' Number of k points in string:', nkstr
325 : 0 : call wrtout(std_out,message,'COLL')
326 : :
327 [ # # ]: 0 : if(nsppol==1)then
328 : : write(message, '(a,i5,a,i5)')&
329 : 0 : & ' From band number',bdberry(1),' to band number',bdberry(2)
330 : : else
331 : : write(message, '(a,i5,a,i5,a,a,a,i5,a,i5,a)')&
332 : 0 : & ' From band number',bdberry(1),' to band number',bdberry(2),' for spin up,',&
333 : 0 : & ch10,&
334 : 0 : & ' from band number',bdberry(3),' to band number',bdberry(4),' for spin down.'
335 : : end if
336 : 0 : call wrtout(ab_out,message,'COLL')
337 : 0 : call wrtout(std_out,message,'COLL')
338 : :
339 : : ! DEBUG
340 : : ! write(std_out,*)' berryphase : find nkpt,nkstr,nstr=',nkpt,nkstr,nstr
341 : : ! stop
342 : : ! ENDDEBUG
343 : :
344 : : ! Build the different strings
345 [ # # # # : 0 : ABI_MALLOC(ikstr,(nkstr,nstr))
# # # # #
# # # # #
# # ]
346 : :
347 : 0 : iunmark=1
348 [ # # ]: 0 : kpt_mark(:)=0
349 [ # # ]: 0 : do istr = 1, nstr
350 [ # # ]: 0 : do while(kpt_mark(iunmark)/=0)
351 : 0 : iunmark = iunmark + 1
352 : : end do
353 : 0 : ikstr(1, istr) = iunmark
354 : 0 : kpt_mark(iunmark)=1
355 [ # # ]: 0 : do jkstr = 2, nkstr
356 : 0 : ikstr(jkstr,istr)=ikpt_dk(ikstr(jkstr-1,istr))
357 : 0 : kpt_mark(ikstr(jkstr,istr))=1
358 : : end do
359 : : end do ! istr
360 : :
361 : : ! DEBUG
362 : : ! do istr = 1,nstr
363 : : ! do jkstr = 1,nkstr
364 : : ! if (shift_g(ikstr(jkstr,istr))) then
365 : : ! write(99,*) ikstr(jkstr,istr),'true'
366 : : ! else
367 : : ! write(99,*) ikstr(jkstr,istr),'false'
368 : : ! end if
369 : : ! end do
370 : : ! end do
371 : : ! ENDDEBUG
372 : :
373 [ # # ]: 0 : ABI_FREE(ikpt_dk)
374 : : ! DEBUG!
375 : : ! write(100,*) 'list all the k points strings:'
376 : : ! do istr=1,nstr
377 : : ! write(100,*) (ikstr(jkstr,istr),jkstr=1,nkstr)
378 : : ! end do
379 : : ! ENDDEBUG!
380 : :
381 : : ! *****************************************************************************
382 : : ! Find the location of each wavefunction
383 [ # # # # : 0 : ABI_MALLOC(cg_index,(mband,nkpt,nsppol))
# # # # #
# # # # #
# # # # #
# ]
384 : :
385 : 0 : icg = 0
386 [ # # ]: 0 : do isppol=1,nsppol
387 [ # # ]: 0 : do ikpt=1,nkpt_
388 : 0 : nband_k=nband(ikpt+(isppol-1)*nkpt_)
389 : 0 : npw_k=npwarr(ikpt)
390 [ # # ]: 0 : do iband=1,nband_k
391 : 0 : cg_index(iband,ikpt,isppol)=(iband-1)*npw_k*nspinor+icg
392 : : end do
393 : 0 : icg=icg+npw_k*nspinor*nband(ikpt)
394 : : end do
395 : : end do
396 : :
397 : : ! change5
398 [ # # ]: 0 : if (mkmem/=0) then
399 : : ! Find the planewave vectors and their indexes for each k point
400 [ # # # # : 0 : ABI_MALLOC(kg_kpt,(3,mpw*nspinor,nkpt_))
# # # # #
# # # # #
# # # # ]
401 [ # # # # : 0 : kg_kpt(:,:,:) = 0
# # ]
402 : : index1 = 0
403 [ # # ]: 0 : do ikpt=1,nkpt_
404 : 0 : npw_k=npwarr(ikpt)
405 [ # # ]: 0 : do ipw=1,npw_k*nspinor
406 [ # # ]: 0 : kg_kpt(1:3,ipw,ikpt)=kg(1:3,ipw+index1)
407 : : end do
408 : 0 : index1=index1+npw_k*nspinor
409 : : end do
410 : : end if !change5
411 : : ! *****************************************************************************
412 [ # # # # : 0 : ABI_MALLOC(det_string,(2, nstr))
# # # # #
# # # ]
413 [ # # # # : 0 : ABI_MALLOC(det_tmp,(2, nstr))
# # # # ]
414 [ # # # # : 0 : ABI_MALLOC(polberry,(nstr))
# # # # #
# # # ]
415 : :
416 : : ! Initialize berry phase polarization for each spin and the total one
417 [ # # ]: 0 : polb(1:nsppol) = 0.0_dp
418 : 0 : polbtot=0.0_dp
419 : :
420 : : ! Loop over spins
421 [ # # ]: 0 : do isppol=1,nsppol
422 : :
423 : 0 : minband=bdberry(2*isppol-1)
424 : 0 : maxband=bdberry(2*isppol)
425 : :
426 [ # # ]: 0 : if(minband<1)then
427 : 0 : write(message,'(a,i0,a)')' The band limit minband=',minband,', is lower than 0.'
428 : 0 : ABI_BUG(message)
429 : : end if
430 : :
431 [ # # ]: 0 : if(maxband<1)then
432 : 0 : write(message,'(a,i0,a)')' The band limit maxband=',maxband,', is lower than 0.'
433 : 0 : ABI_BUG(message)
434 : : end if
435 : :
436 [ # # ]: 0 : if(maxband<minband)then
437 : 0 : write(message,'(a,i0,a,i0)')' maxband=',maxband,', is lower than minband=',minband
438 : 0 : ABI_BUG(message)
439 : : end if
440 : :
441 : : ! Initialize det_string and det_average
442 [ # # # # ]: 0 : det_string(1, 1:nstr) = 1.0_dp; det_string(2, 1:nstr) = 0.0_dp
443 : 0 : det_average(1:2)=0.0_dp; det_average(2)=0.0_dp
444 : :
445 : : ! Loop over strings
446 [ # # ]: 0 : do istr = 1, nstr
447 : :
448 : : ! change7
449 : 0 : read_k = 0
450 : :
451 : : ! DEBUG!
452 : : ! write(100,'(a,i4)') 'This is in string', istr
453 : : ! ENDDEBUG!
454 : :
455 : : ! Loop over k points per string
456 [ # # # # : 0 : ABI_MALLOC(determinant,(2, nkstr))
# # # # #
# # # ]
457 : :
458 [ # # ]: 0 : do jkstr=1,nkstr
459 : :
460 [ # # # # : 0 : ABI_MALLOC(cmatrix,(2,maxband,maxband))
# # # # #
# # # # #
# # ]
461 [ # # ]: 0 : if(jkstr < nkstr) then
462 : 0 : lkstr=jkstr+1
463 : : else
464 : 0 : lkstr= jkstr+1-nkstr
465 : : end if
466 : 0 : jkstr_ori=ikstr(jkstr,istr)
467 : 0 : lkstr_ori=ikstr(lkstr,istr)
468 : :
469 : : ! change9
470 : 0 : lkstr_ori_=lkstr_ori
471 : 0 : tr(1) = 1.0_dp
472 : 0 : tr(2) = 1.0_dp
473 [ # # ]: 0 : if (kptopt==2) then
474 [ # # ]: 0 : if (read_k == 0) then
475 [ # # ]: 0 : if (kpt_flag(jkstr_ori)/=0) then
476 : 0 : tr(1) = -1.0_dp
477 : 0 : jkstr_ori = kpt_flag(jkstr_ori)
478 : : end if
479 [ # # ]: 0 : if (kpt_flag(lkstr_ori)/=0) then
480 : 0 : tr(2) = -1.0_dp
481 : 0 : lkstr_ori = kpt_flag(lkstr_ori)
482 : : end if
483 : : else !read_k
484 [ # # ]: 0 : if (kpt_flag(jkstr_ori)/=0) then
485 : 0 : tr(-1*read_k+3) = -1.0_dp
486 : 0 : jkstr_ori = kpt_flag(jkstr_ori)
487 : : end if
488 [ # # ]: 0 : if (kpt_flag(lkstr_ori)/=0) then
489 : 0 : tr(read_k) = -1.0_dp
490 : 0 : lkstr_ori = kpt_flag(lkstr_ori)
491 : : end if
492 : : end if !read_k
493 : : end if !kptopt
494 : : ! change9
495 : :
496 : 0 : nband_k=nband(jkstr_ori+(isppol-1)*nkpt_)
497 [ # # ]: 0 : if(nband_k<maxband)then
498 : 0 : write(message,'(a,i0,a,i0)')' maxband=',maxband,', is larger than nband(j,isppol)=',nband_k
499 : 0 : ABI_BUG(message)
500 : : end if
501 : :
502 : 0 : nband_k=nband(lkstr_ori+(isppol-1)*nkpt_)
503 [ # # ]: 0 : if(nband_k<maxband)then
504 : : write(message,'(a,i0,a,i0)')&
505 : 0 : & ' maxband=',maxband,', is larger than nband(l,isppol)=',nband_k
506 : 0 : ABI_BUG(message)
507 : : end if
508 : :
509 [ # # ]: 0 : if (jkstr==1) read_k = 2
510 : : ! Compute the overlap matrix <u_k|u_k+b>
511 [ # # # # : 0 : cmatrix(1:2,1:maxband,1:maxband)=zero
# # ]
512 : 0 : jj = read_k
513 : 0 : ii = -1*read_k+3
514 [ # # ]: 0 : if(.not. shift_g(lkstr_ori_) ) then
515 : : ! Change3
516 [ # # ]: 0 : do ipw=1,npwarr(jkstr_ori)
517 [ # # ]: 0 : do jpw=1,npwarr(lkstr_ori)
518 : :
519 : : ! Check if Fourier components of jkstr and jkstr+1 matches
520 : :
521 : : if((tr(ii)*kg_kpt(1,ipw,jkstr_ori)==tr(jj)*kg_kpt(1,jpw,lkstr_ori))&
522 : : & .and.(tr(ii)*kg_kpt(2,ipw,jkstr_ori) == tr(jj)*kg_kpt(2,jpw,lkstr_ori))&
523 [ # # # # : 0 : & .and.(tr(ii)*kg_kpt(3,ipw,jkstr_ori) == tr(jj)*kg_kpt(3,jpw,lkstr_ori)))&
# # ]
524 : 0 : & then
525 : :
526 [ # # ]: 0 : do iband=minband,maxband
527 : 0 : cg_index_iband=cg_index(iband,jkstr_ori,isppol)
528 [ # # ]: 0 : do jband=minband,maxband
529 : 0 : cg_index_jband=cg_index(jband,lkstr_ori,isppol)
530 : :
531 : : cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
532 : : & cg(1,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)+&
533 : 0 : & tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)
534 : : cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
535 : : & cg(1,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)-&
536 : 0 : & tr(ii)*cg(2,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)
537 : :
538 : : end do !jband
539 : : end do !iband
540 : : exit !stop loop over jpw if Fourier components of jkstr and jkstr + 1 matches
541 : : end if
542 : :
543 : : end do ! jpw
544 : : end do ! ipw
545 : :
546 : : ! But there is a special pair of k points which involves the shift of a
547 : : ! G vector
548 : :
549 : : else
550 : :
551 [ # # ]: 0 : dg(:) = -1*nint(tr(jj)*kpt(:,lkstr_ori)-tr(ii)*kpt(:,jkstr_ori)-dk(:))
552 : :
553 : : ! DEBUG
554 : : ! write(100,*)dg
555 : : ! write(100,*)kberry(:,iberry)
556 : : ! write(100,*)''
557 : : ! ENDDEBUG
558 : :
559 : : ! change4
560 [ # # ]: 0 : do ipw=1,npwarr(jkstr_ori)
561 [ # # ]: 0 : do jpw=1,npwarr(lkstr_ori)
562 : :
563 : : ! Check if Fourier components of jkstr and jkstr+1
564 : : ! matches by comparing the G vectors
565 : :
566 : : if((tr(ii)*kg_kpt(1,ipw,jkstr_ori)==tr(jj)*kg_kpt(1,jpw,lkstr_ori)-dg(1))&
567 : : & .and.(tr(ii)*kg_kpt(2,ipw,jkstr_ori) == tr(jj)*kg_kpt(2,jpw,lkstr_ori)-dg(2))&
568 [ # # # # : 0 : & .and.(tr(ii)*kg_kpt(3,ipw,jkstr_ori) == tr(jj)*kg_kpt(3,jpw,lkstr_ori)-dg(3)))&
# # ]
569 : 0 : & then
570 : :
571 [ # # ]: 0 : do iband=minband,maxband
572 : 0 : cg_index_iband=cg_index(iband,jkstr_ori,isppol)
573 : :
574 [ # # ]: 0 : do jband=minband,maxband
575 : 0 : cg_index_jband=cg_index(jband,lkstr_ori,isppol)
576 : :
577 : : cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
578 : : & cg(1,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)+&
579 : 0 : & tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)
580 : : cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
581 : : & cg(1,ipw+cg_index_iband)*tr(jj)*cg(2,jpw+cg_index_jband)-&
582 : 0 : & tr(ii)*cg(2,ipw+cg_index_iband)*cg(1,jpw+cg_index_jband)
583 : :
584 : : end do ! jband
585 : : end do ! iband
586 : : exit !stop loop over jpw if Fourier components of jkstr and jkstr + 1 matches
587 : : end if
588 : : end do ! jpw
589 : : end do ! ipw
590 : : end if
591 : :
592 : : ! Compute the determinant of cmatrix(1:2,minband:maxband, minband:maxband)
593 : :
594 : 0 : band_in = maxband - minband + 1
595 : :
596 [ # # # # : 0 : ABI_MALLOC(ipvt,(maxband))
# # # # #
# # # ]
597 [ # # # # : 0 : ABI_MALLOC(zgwork,(2,1:maxband))
# # # # #
# # # ]
598 : :
599 : : ! Last argument of zgedi means calculate determinant only.
600 : 0 : call dzgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info)
601 : 0 : call dzgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,10)
602 : :
603 [ # # ]: 0 : ABI_FREE(zgwork)
604 [ # # ]: 0 : ABI_FREE(ipvt)
605 : :
606 : 0 : fac=exp(log(10._dp)*det(1,2))
607 : : determinant(1, jkstr) = fac*(det(1,1)*cos(log(10._dp)*det(2,2)) - &
608 : 0 : & det(2,1)*sin(log(10._dp)*det(2,2)))
609 : : determinant(2, jkstr) = fac*(det(1,1)*sin(log(10._dp)*det(2,2)) + &
610 : 0 : & det(2,1)*cos(log(10._dp)*det(2,2)))
611 : : ! DEBUG!
612 : : ! write(100,*) 'det',jkstr,lkstr,'=', determinant(1:2,jkstr)
613 : : ! ENDDEBUG!
614 : :
615 : : det_tmp(1,istr) = det_string(1,istr)*determinant(1,jkstr) - &
616 : 0 : & det_string(2,istr)*determinant(2,jkstr)
617 : : det_tmp(2,istr) = det_string(1,istr)*determinant(2,jkstr) + &
618 : 0 : & det_string(2,istr)*determinant(1,jkstr)
619 [ # # ]: 0 : det_string(1:2,istr) = det_tmp(1:2,istr)
620 : :
621 [ # # ]: 0 : ABI_FREE(cmatrix)
622 : :
623 : : ! Close loop over k points along string
624 : 0 : read_k = -1*read_k + 3 ! read_k=2 <-> read_k=1
625 : : end do
626 : :
627 : : ! DEBUG!
628 : : ! write(100,*) 'det_string =', det_string(1:2,istr)
629 : : ! write(100,*)
630 : : ! ENDDEBUG!
631 : :
632 : 0 : det_average(1) = det_average(1) + det_string(1,istr)/nstr
633 : 0 : det_average(2) = det_average(2) + det_string(2,istr)/nstr
634 : :
635 [ # # ]: 0 : ABI_FREE(determinant)
636 : :
637 : : ! Close loop over strings
638 : : end do
639 : :
640 : :
641 : : ! *****************************************************************************
642 : : ! Calculate the electronic contribution to the polarization
643 : :
644 : 0 : write(message,'(a,a)')ch10,&
645 : 0 : & ' Compute the electronic contribution to polarization'
646 : 0 : call wrtout(std_out,message,'COLL')
647 : :
648 : : ! First berry phase that corresponds to det_average
649 : 0 : phase0 = atan2(det_average(2),det_average(1))
650 : 0 : det_mod = det_average(1)**2+det_average(2)**2
651 : :
652 : : ! Then berry phase that corresponds to each string relative to the average
653 [ # # ]: 0 : do istr = 1, nstr
654 : : rel_string(1) = (det_string(1,istr)*det_average(1) + &
655 : 0 : det_string(2,istr)*det_average(2))/det_mod
656 : : rel_string(2) = (det_string(2,istr)*det_average(1) - &
657 : 0 : det_string(1,istr)*det_average(2))/det_mod
658 : 0 : dphase = atan2(rel_string(2),rel_string(1))
659 : 0 : polberry(istr) = (2.0_dp/nsppol)*(phase0+dphase)/two_pi
660 : 0 : polb(isppol) = polb(isppol) + polberry(istr)/nstr
661 : : end do
662 : :
663 : : ! Output berry phase polarization
664 : 0 : write(message,'(a,10x,a,10x,a)')ch10,&
665 : 0 : & 'istr','polberry(istr)'
666 : 0 : call wrtout(std_out,message,'COLL')
667 [ # # ]: 0 : do istr=1,nstr
668 : 0 : write(message,'(10x,i4,7x,e16.9)')istr,polberry(istr)
669 : 0 : call wrtout(std_out,message,'COLL')
670 : : end do
671 : :
672 : : write(message,'(9x,a,7x,e16.9,1x,a,i4,a,a)')&
673 : 0 : & 'total',polb(isppol),'(isppol=',isppol,')',ch10
674 : 0 : call wrtout(std_out,message,'COLL')
675 : :
676 : 0 : polbtot=polbtot+polb(isppol)
677 : :
678 : : end do ! isppol
679 : :
680 [ # # ]: 0 : ABI_FREE(polberry)
681 [ # # ]: 0 : ABI_FREE(det_tmp)
682 [ # # ]: 0 : ABI_FREE(det_string)
683 [ # # ]: 0 : ABI_FREE(ikstr)
684 [ # # ]: 0 : ABI_FREE(cg_index)
685 : : ! change6
686 [ # # ]: 0 : if (mkmem /=0) then
687 [ # # ]: 0 : ABI_FREE(kg_kpt)
688 : : end if
689 : : ! *****************************************************************************
690 : : ! Reindex xcart according to atom and type
691 : 0 : call xred2xcart(natom,rprimd,xcart,xred)
692 [ # # # # : 0 : ABI_MALLOC(xcart_reindex,(3,natom,ntypat))
# # # # #
# # # # #
# # ]
693 : 0 : index=1
694 [ # # ]: 0 : do itypat=1,ntypat
695 [ # # ]: 0 : do iattyp=1,nattyp(itypat)
696 : 0 : iatom=atindx1(index)
697 [ # # ]: 0 : xcart_reindex(1:3,iattyp,itypat) = xcart(1:3,iatom)
698 : 0 : index = index+1
699 : : end do
700 : : end do
701 : :
702 : : ! Compute the ionic contribution to the polarization
703 : 0 : politot = 0.0_dp
704 : 0 : write(message,'(a)')' Compute the ionic contributions'
705 : 0 : call wrtout(std_out,message,'COLL')
706 : :
707 : 0 : write(message,'(a,2x,a,2x,a,15x,a)')ch10,&
708 : 0 : & 'itypat', 'iattyp', 'polion'
709 : 0 : call wrtout(std_out,message,'COLL')
710 : :
711 [ # # ]: 0 : do itypat=1,ntypat
712 [ # # ]: 0 : do iattyp=1,nattyp(itypat)
713 : : polion=zion(itypat)*nkstr*&
714 [ # # ]: 0 : & dot_product(xcart_reindex(1:3,iattyp,itypat),gpard(1:3))
715 : : ! Fold into interval (-1,1)
716 : 0 : polion=polion-2._dp*nint(polion/2.0_dp)
717 : 0 : politot=politot+polion
718 : 0 : write(message,'(2x,i2,5x,i2,10x,e16.9)') itypat,iattyp,polion
719 : 0 : call wrtout(std_out,message,'COLL')
720 : : end do
721 : : end do
722 : :
723 : : ! Fold into interval (-1,1) again
724 : 0 : politot=politot-2.0_dp*nint(politot/2.0_dp)
725 : :
726 : 0 : write(message,'(9x,a,7x,es19.9)') 'total',politot
727 : 0 : call wrtout(std_out,message,'COLL')
728 : :
729 [ # # ]: 0 : ABI_FREE(xcart_reindex)
730 : :
731 : : ! Compute the total polarizations
732 : :
733 : 0 : poltot=politot+polbtot
734 : :
735 : 0 : write(message,'(a,a)')ch10,&
736 : 0 : & ' Summary of the results'
737 : 0 : call wrtout(std_out,message,'COLL')
738 : 0 : call wrtout(ab_out,message,'COLL')
739 : :
740 : : write(message,'(a,es19.9)')&
741 : 0 : & ' Electronic Berry phase ' ,polbtot
742 : 0 : call wrtout(std_out,message,'COLL')
743 : 0 : call wrtout(ab_out,message,'COLL')
744 : :
745 : : write(message,'(a,es19.9)') &
746 : 0 : & ' Ionic phase ', politot
747 : 0 : call wrtout(std_out,message,'COLL')
748 : 0 : call wrtout(ab_out,message,'COLL')
749 : :
750 : : write(message,'(a,es19.9)') &
751 : 0 : & ' Total phase ', poltot
752 : 0 : call wrtout(std_out,message,'COLL')
753 : 0 : call wrtout(ab_out,message,'COLL')
754 : :
755 : 0 : poltot=poltot-2.0_dp*nint(poltot/2._dp)
756 : : write(message,'(a,es19.9)') &
757 : 0 : & ' Remapping in [-1,1] ', poltot
758 : 0 : call wrtout(std_out,message,'COLL')
759 : 0 : call wrtout(ab_out,message,'COLL')
760 : :
761 : : ! Transform the phase into a polarization
762 : 0 : fac = 1._dp/(gmod*nkstr)
763 : 0 : fac = fac/ucvol
764 : 0 : pol = fac*poltot
765 : :
766 : 0 : write(message,'(a,a,es19.9,a,a,a,es19.9,a,a)')ch10,&
767 : 0 : & ' Polarization ', pol,' (a.u. of charge)/bohr^2',ch10,&
768 : 0 : & ' Polarization ', pol*(e_Cb)/(Bohr_Ang*1d-10)**2,&
769 : 0 : & ' C/m^2',ch10
770 : 0 : call wrtout(std_out,message,'COLL')
771 : 0 : call wrtout(ab_out,message,'COLL')
772 : :
773 : : end do ! iberry
774 : :
775 [ # # ]: 0 : ABI_FREE(shift_g)
776 [ # # ]: 0 : ABI_FREE(kpt)
777 [ # # ]: 0 : ABI_FREE(kg_dum)
778 : :
779 [ # # ]: 0 : end subroutine berryphase
780 : : !!***
781 : :
782 : : end module m_berryphase
783 : : !!***
|
