TABLE OF CONTENTS


ABINIT/m_berrytk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_berrytk

FUNCTION

COPYRIGHT

  Copyright (C) 2000-2022 ABINIT  group (MVeithen)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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

ABINIT/polcart [ Functions ]

[ Top ] [ Functions ]

NAME

 polcart

FUNCTION

 Transform polarization from reduced to cartesian coordinates,
 divide by ucvol and write the result to an output file

INPUTS

 pel(3)   = reduced coordinates of the electronic polarization
 pion(3)  = reduced coordinates of the ionic polarization
 polunit  = units used for the output of the polarization
             1 : use atomic units
             2 : use MKS units
             3 : use both atomic and MKS units
 rprimd(3,3) = dimensional primitive translations (bohr)
 ucvol    = volume of the primitive unit cell
 unit_out = unit for output of the results
 usepaw = 1 for PAW computation, zero else

OUTPUT

 pel_cart(3) = cartesian coords of the electronic polarization
               in atomic units
 pelev(3)= expectation value polarization term (PAW only) in cartesian coordinates
 pion_cart(3)= cartesian coords of the ionic polarization
               in atomic units
 ptot_cart(3)= cartesian coords of the total polarization
               in atomic units

NOTES

 - The sum of the electronic and ionic Berry phase is folded into
   [-1,1] before it is transformed to cartesian coordinates.
   This means that in some cases, ptot_cart /= pel_cart + pion_cart

 - pel and pion do not take into account the factor 1/ucvol.
   At the opposite, this factor is taken into account in
   pel_cart and pion_cart
 - unit_out = 0 is allowed, in this case, there will be no
   output of the results

SOURCE

661 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 ! ***********************************************************************
681 !!REC Note ptot has already been folded and kept onto same branch
682 !unless ptot=0d0, in which case ptot has not been computed yet
683  if( sum(abs(red_ptot(:))) < tol8 )then
684    ptot(:) = pel(:) + pion(:)
685 !  Fold ptot into [-1, 1]
686    do idir = 1, 3
687      ptot(idir) = ptot(idir) - 2_dp*nint(ptot(idir)/2_dp)
688    end do
689  else !!REC
690    ptot=red_ptot !!REC
691  end if !!REC
692 
693 !Transform pel, pion and ptot to cartesian coordinates
694  pel_cart(:) = zero ; pion_cart(:) = zero ; ptot_cart(:) = zero
695  do idir = 1, 3
696    pel_cart(idir) =  rprimd(idir,1)*pel(1) + rprimd(idir,2)*pel(2) + &
697 &   rprimd(idir,3)*pel(3)
698    pion_cart(idir) = rprimd(idir,1)*pion(1) + rprimd(idir,2)*pion(2) + &
699 &   rprimd(idir,3)*pion(3)
700    ptot_cart(idir) = rprimd(idir,1)*ptot(1) + rprimd(idir,2)*ptot(2) + &
701 &   rprimd(idir,3)*ptot(3)
702  end do
703 
704 !Divide by the unit cell volume
705  pel_cart(:)  = pel_cart(:)/ucvol
706  pion_cart(:) = pion_cart(:)/ucvol
707  ptot_cart(:)  = ptot_cart(:)/ucvol
708 !pelev is either zero (in NCPP case), or possibly non-zero (in PAW case)
709 !however, in the PAW case, it is already implicitly included in the computation
710 !of pel and should not be added in here. Below in the PAW case we report its
711 !value to the user, which is helpful for comparison to the USPP theory.
712 !13 June 2012 J Zwanziger
713 !note that pelev was AlREADY in cartesian frame.
714 !ptot_cart(:)  = (ptot_cart(:)+pelev(:))/ucvol
715 
716 !Write the results to unit_out (if /= 0)
717 !Use the coordinates specified by the value polunit
718 
719 !Atomic units
720 
721  if (((polunit == 1).or.(polunit == 3)).and.(unit_out /= 0)) then
722 
723 !  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,&
724 !  &   ' Polarization in cartesian coordinates (a.u.):',ch10,&
725 !  &   ' (the sum of the electronic and ionic Berry phase',&
726 !  &   ' has been fold into [-1, 1])',ch10,&
727 !  &   '     Electronic berry phase:       ', (pel_cart(idir), idir = 1, 3), ch10,&
728 !  &   '     Expectation value (PAW only): ', (pelev(idir)/ucvol, idir = 1, 3), ch10,&
729 !  &   '     Ionic:                        ', (pion_cart(idir), idir = 1, 3), ch10, &
730 !  &   '     Total:                        ', (ptot_cart(idir), idir = 1, 3)
731 !  call wrtout(unit_out,message,'COLL')
732    write(message,'(7(a),3(e16.9,2x),a,a,3(e16.9,2x))')ch10,&
733 &   ' Polarization in cartesian coordinates (a.u.):',ch10,&
734 &   ' (the sum of the electronic and ionic Berry phase',&
735 &   ' has been folded into [-1, 1])',ch10,&
736 &   '     Electronic berry phase:       ', (pel_cart(idir), idir = 1, 3)
737    call wrtout(unit_out,message,'COLL')
738    if(usepaw==1) then
739      write(message,'(a,3(e16.9,2x))')&
740 &     '     ...includes PAW on-site term: ', (pelev(idir)/ucvol, idir = 1, 3)
741      call wrtout(unit_out,message,'COLL')
742    end if
743    write(message,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')&
744 &   '     Ionic:                        ', (pion_cart(idir), idir = 1, 3), ch10, &
745 &   '     Total:                        ', (ptot_cart(idir), idir = 1, 3)
746    call wrtout(unit_out,message,'COLL')
747 
748  end if
749 
750 !MKS units
751 
752  if (((polunit == 2).or.(polunit == 3)).and.(unit_out /= 0)) then
753 
754    pel_mks(:)  = pel_cart(:)*(e_Cb)/(Bohr_Ang*1d-10)**2
755    pion_mks(:) = pion_cart(:)*(e_Cb)/(Bohr_Ang*1d-10)**2
756    pelev_mks(:) = pelev(:)/ucvol*(e_Cb)/(Bohr_Ang*1d-10)**2
757    ptot_mks(:) = (ptot_cart(:))*(e_Cb)/(Bohr_Ang*1d-10)**2
758 
759    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,&
760 &   ' Polarization in cartesian coordinates (C/m^2):',ch10,&
761 &   ' (the sum of the electronic and ionic Berry phase',&
762 &   ' has been folded into [-1, 1])',ch10,&
763 &   '     Electronic berry phase:       ', (pel_mks(idir), idir = 1, 3)
764    call wrtout(unit_out,message,'COLL')
765    if(usepaw==1) then
766      write(message,'(a,3(e16.9,2x))')&
767 &     '     ...includes PAW on-site term: ', (pelev_mks(idir), idir = 1, 3)
768      call wrtout(unit_out,message,'COLL')
769    end if
770    write(message,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')&
771 &   '     Ionic:                        ', (pion_mks(idir), idir = 1, 3), ch10, &
772 &   '     Total:                        ', (ptot_mks(idir), idir = 1, 3)
773    call wrtout(unit_out,message,'COLL')
774 
775  end if
776 
777 end subroutine polcart

ABINIT/smatrix [ Functions ]

[ Top ] [ Functions ]

NAME

 smatrix

FUNCTION

 Compute the overlap matrix between the k-points k and k + dk.
 Depending on the value of job and ddkflag, compute also its determinant,
 its inverse and the product of its inverse with the wavefunctions at k.

INPUTS

 cg(2,mcg_k) = planewave coefficients of wavefunctions at k
 cgq(2,mcg_q) = planewave coefficients of wavefunctions at q = k + dk
 ddkflag = 1 : compute product of the inverse overlap matrix
               with the wavefunction at k (job = 1 or 11)
           0 : do not compute the product of the inverse overlap matrix
               with the wavefunction at k
 icg = shift applied to the wavefunctions in the array cg
 icg1 = shift applied to the wavefunctions in the array cgq
 itrs = variable that governs the use of time-reversal symmetry
        when converting the wavefunctions from the iBZ to the fBZ
 job = type of calculation
        0 : update overlap matrix only
        1 : like 0 but also compute inverse of the overlap matrix
       10 : like 0 but also compute determinant of the overlap matrix
       11 : like 0 but also compute determinant and inverse of the overlap matrix
       20 : like 0 but also transfer cgq to cg1_k without multiplication by S^{-1}
       21 : like 1 but also transfer cgq to cg1_k without multiplication by S^{-1}
 maxbd = used in case ddkflag = 1, defines the highest band for
         which the ddk will be computed
 mcg_k = second dimension of cg
 mcg_q = second dimension of cg_q
 mcg1_k = second dimension of cg1_k, should be equal to
          mpw*nsppol*nspinor*(maxbd - minbd + 1)
 minbd = used in case ddkflag = 1, defines the lowest band for
         which the ddk will be computed
 mpw = maximum dimensioned size of npw
 mband_occ = max number of occupied valence bands for both spins
 nband_occ = number of (occupied) valence bands
 npw_k1 = number of plane waves at k
 npw_k2 = number of plane waves at k + dk
 nspinor = number of spinorial components of the wavefunctions
 nsppol = 1 for unpolarized, 2 for spin-polarized
 pwind_k = array used to compute the overlap matrix
 pwnsfac = phase factors for non-symmorphic translations
 shifrbd = shift applied to the location of the WF in the cg-array
           after each loop over bands
           0 : apply no shift, this is allowed in case cg
               contains only the wf of one single band
           1 : apply a shift of npw_k1*nspinor, this is the usual option
               when cg contains the wf for all occupied bands
 smat_k_paw : overlap matrix due to on-site PAW terms between different bands
              at k and k+b. Only relevant when usepaw = 1, and is to be computed
              previously by smatrix_k_paw.F90
 usepaw = flag governing use of PAW: 0 means no PAW, 1 means PAW is used

OUTPUT

 cg1_k(2,mcg1_k) = product of the inverse overlap matrix with the
                   wavefunctions at k; computed in case job = 1 or 11;
                   or just cgq in case of job = 20 or 21
 dtm_k(2) = determinant of the overlap matrix between k and k + dk;
            computed in case job = 10 or 11
 smat_inv = inverse of the overlap matrix

SIDE EFFECTS

 Input/Output
 sflag_k(iband) = 1 if the elements smat_k(:,iband,:) are up to date
                    -> they will not be recomputed
                  0 the elements smat_k(:,iband,:) will be recomputed
      at the end of the routine, sflag_k(1:mband_occ) = 1
      (the whole overlap matrix is up to date)
 smat_k = overlap matrix between k, k + dk
          only the lines for which sflag_k = 0 are computed
          smat_k(:,n,m) = < u_{n,k} | u_{m,k+dk} >

NOTES

 This routine is quite flexible in the way it deals with the wavefunctions:
  - cg (WF at k) can contain either the whole WF array (all k-points
    and bands), in which case the location of the WF at k is specified
    by icg and shiftbd = 1, or the WF of a single k-point/band, in which case
    shiftbd = 0 and icg = 0.
  - cgq (WF at k + dk) can contain either the whole WF array (all k-points
    and bands), in which case the location of the WF at k is specified
    by icg1, or the WF of a single k-point, in which case
    icg1 = 0. cgq must contain the WF of ALL occupied bands.
  - cg1_k can either be computed for all valence bands or
    for a group of valence bands defined by minbd and maxbd.

SOURCE

131 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 &  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  integer,allocatable :: ipvt(:)
158 ! integer,allocatable :: my_pwind_k(:) ! used in debugging below
159  real(dp) :: det(2,2)
160  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  ABI_MALLOC(ipvt,(nband_occ))
173  ABI_MALLOC(zgwork,(2,nband_occ))
174  ABI_MALLOC(vect1,(2,0:mpw*nspinor))
175  ABI_MALLOC(vect2,(2,0:mpw*nspinor))
176  vect1(:,0) = zero ; vect2(:,0) = zero
177 
178 !Check if the values of ddkflag and job are compatible
179 
180  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 &   ' job is equal to ',job,ch10,&
183 &   ' while only the values job = 0, 1, 10, 11, 20, or 21 are allowed.'
184    ABI_ERROR(message)
185  end if
186 
187  if (ddkflag == 1) then
188    if ((job/=1).and.(job/=11)) then
189      write(message,'(a,i0,a,a)')&
190 &     ' job is equal to ',job,ch10,&
191 &     ' while ddkflag = 1. This is not allowed.'
192      ABI_ERROR(message)
193    end if
194  end if
195 
196 !Check the values of sflag_k
197  do iband=1,nband_occ
198    if (sflag_k(iband)/=0 .and. sflag_k(iband)/=1)then
199      write(message,'(3a,i4,a,i4)')&
200 &     '  The content of sflag_k must be 0 or 1.',ch10,&
201 &     '  However, for iband=',iband,', sflag_k(iband)=',sflag_k(iband)
202      ABI_ERROR(message)
203    end if
204  end do
205 
206 !Check if shiftbd is consistent with sflag_k
207  if (shiftbd == 0) then
208    count = 0
209    do iband = 1, nband_occ
210      if (sflag_k(iband) == 0) count = count + 1
211    end do
212    if (count > 1) then
213      message = 'in case shiftbd = 0, only 1 element of sflag can be 0'
214      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  do iband = 1, nband_occ
316 
317    if (sflag_k(iband) == 0) then
318 
319      pwmin = (iband-1)*npw_k1*nspinor*shiftbd
320      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      if (itrs==1.or.itrs==11) then  ! take complex conjugate of bra
328        do ispinor = 1, nspinor
329          spnshft_k1=(ispinor-1)*npw_k1
330          do ipw = 1,npw_k1
331            vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
332 &           +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 &           -cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
335          end do
336        end do
337      else
338        do ispinor=1, nspinor
339          spnshft_k1=(ispinor-1)*npw_k1
340          do ipw = 1,npw_k1
341            vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
342 &           -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 &           +cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
345          end do
346        end do
347      end if
348 
349 !
350      if (npw_k1*nspinor < mpw*nspinor) vect1(:,npw_k1*nspinor+1:mpw*nspinor) = zero
351 
352      do jband = 1, nband_occ
353 
354        pwmin = (jband-1)*npw_k2*nspinor
355        pwmax = pwmin + npw_k2*nspinor
356 
357        if (itrs==10.or.itrs==11) then ! take complex conjugate of ket
358          do ispinor=1, nspinor
359            spnshft_k2=(ispinor-1)*npw_k2
360            do ipw = 1, npw_k2
361              vect2(1,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw) &
362 &             +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 &             -cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw)
365            end do
366          end do
367        else
368          do ispinor=1, nspinor
369            spnshft_k2=(ispinor-1)*npw_k2
370            do ipw = 1, npw_k2
371              vect2(1,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw) &
372 &             -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 &             +cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw)
375            end do
376          end do
377        end if
378 
379        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        call overlap_g(doti,dotr,mpw,npw_k1,npw_k2,nspinor,pwind_k,vect1,vect2)
394 
395        smat_k(1,iband,jband) = dotr
396        smat_k(2,iband,jband) = doti
397 
398        if (usepaw == 1) then
399          smat_k(1,iband,jband) = smat_k(1,iband,jband)+smat_k_paw(1,iband,jband)
400          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  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  if ((job==1).or.(job==10).or.(job==11).or.(job==21)) then
432 
433    smat_inv(:,:,:) = smat_k(:,:,:)
434 
435 !  DEBUG
436 !  write(std_out,*)' smatrix : smat_inv=',smat_inv
437 !  ENDDEBUG
438 
439    dzgedi_job=job; if(job==21) dzgedi_job=1
440 ! TODO: should this be over nband_occ(isppol)?
441    call dzgefa(smat_inv,mband_occ,nband_occ,ipvt,info)
442    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    dtm_k(:) = zero
450    if (job==10 .or. job==11) then
451      fac = exp(log(10._dp)*det(1,2))
452      dtm_k(1) = fac*(det(1,1)*cos(log(10._dp)*det(2,2)) - &
453 &     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 &     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    if (ddkflag == 1) then
461 
462      cg1_k(:,:) = zero
463      jband1 = 0
464 
465      if (itrs == 10 .or. itrs == 11) then
466 
467        do jband = minbd, maxbd
468          jband1 = jband1 + 1
469          do iband = 1, nband_occ
470 
471            do ispinor = 1, nspinor
472              spnshft_k1 = (ispinor-1)*npw_k1
473              spnshft_k2 = (ispinor-1)*npw_k2
474              do ipw = 1, npw_k1
475 
476                jpw = pwind_k(ipw)
477 
478                if (jpw > 0) then
479 
480                  wfr = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
481 &                 -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 &                 +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 &                 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 &                 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        do jband = minbd, maxbd
504          jband1 = jband1 + 1
505          do iband = 1, nband_occ
506 
507            do ispinor = 1, nspinor
508              spnshft_k1 = (ispinor-1)*npw_k1
509              spnshft_k2 = (ispinor-1)*npw_k2
510              do ipw = 1, npw_k1
511 
512                jpw = pwind_k(ipw)
513 
514                if (jpw > 0) then
515 
516                  wfr = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
517 &                 -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 &                 +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 &                 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 &                 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  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    cg1_k(:,:) = zero
547    jband1 = 0
548 
549    if (itrs == 10 .or. itrs == 11) then
550 
551      do jband = minbd, maxbd
552        jband1 = jband1 + 1
553        do ispinor = 1, nspinor
554          spnshft_k1 = (ispinor-1)*npw_k1
555          spnshft_k2 = (ispinor-1)*npw_k2
556          do ipw = 1, npw_k1
557            jpw = pwind_k(ipw)
558 
559            if (jpw > 0) then
560              wfr = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
561 &             -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 &             +cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
564 
565              cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfr
566              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      do jband = minbd, maxbd
578        jband1 = jband1 + 1
579        do ispinor = 1, nspinor
580          spnshft_k1 = (ispinor-1)*npw_k1
581          spnshft_k2 = (ispinor-1)*npw_k2
582          do ipw = 1, npw_k1
583            jpw = pwind_k(ipw)
584 
585            if (jpw > 0) then
586 
587              wfr = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
588 &             -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              +cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
591 
592              cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfr
593              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  ABI_FREE(ipvt)
607  ABI_FREE(zgwork)
608  ABI_FREE(vect1)
609  ABI_FREE(vect2)
610 
611 !DEBUG
612 !write(std_out,*)' dtm_k=',dtm_k(:)
613 !write(std_out,*)' smatrix : exit '
614 !ENDDEBUG
615 
616 end subroutine smatrix