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