TABLE OF CONTENTS


ABINIT/m_berrytk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_berrytk

FUNCTION

COPYRIGHT

  Copyright (C) 2000-2018 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 .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_berrytk
27 
28  use defs_basis
29  use m_abicore
30  use m_errors
31 
32  use m_cgtools,   only : overlap_g
33  use m_hide_lapack,   only : dzgedi, dzgefa
34 
35  implicit none
36 
37  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

PARENTS

      berryphase_new,relaxpol

CHILDREN

      wrtout

SOURCE

687 subroutine polcart(red_ptot,pel,pel_cart,pelev,pion,pion_cart,polunit,&
688 &  ptot_cart,rprimd,ucvol,unit_out,usepaw)
689 
690 
691 !This section has been created automatically by the script Abilint (TD).
692 !Do not modify the following lines by hand.
693 #undef ABI_FUNC
694 #define ABI_FUNC 'polcart'
695 !End of the abilint section
696 
697  implicit none
698 
699 !Arguments ------------------------------------
700 !scalars
701  integer,intent(in) :: polunit,unit_out,usepaw
702  real(dp),intent(in) :: ucvol
703 !arrays
704  real(dp),intent(in) :: red_ptot(3) !!REC
705  real(dp),intent(in) :: pel(3),pelev(3),pion(3),rprimd(3,3)
706  real(dp),intent(out) :: pel_cart(3),pion_cart(3),ptot_cart(3)
707 
708 !Local variables -------------------------
709 !scalars
710  integer :: idir
711  character(len=500) :: message
712 !arrays
713  real(dp) :: pel_mks(3),pelev_mks(3),pion_mks(3),ptot(3),ptot_mks(3)
714 
715 ! ***********************************************************************
716 !!REC Note ptot has already been folded and kept onto same branch
717 !unless ptot=0d0, in which case ptot has not been computed yet
718  if( sum(abs(red_ptot(:))) < tol8 )then
719    ptot(:) = pel(:) + pion(:)
720 !  Fold ptot into [-1, 1]
721    do idir = 1, 3
722      ptot(idir) = ptot(idir) - 2_dp*nint(ptot(idir)/2_dp)
723    end do
724  else !!REC
725    ptot=red_ptot !!REC
726  end if !!REC
727 
728 !Transform pel, pion and ptot to cartesian coordinates
729  pel_cart(:) = zero ; pion_cart(:) = zero ; ptot_cart(:) = zero
730  do idir = 1, 3
731    pel_cart(idir) =  rprimd(idir,1)*pel(1) + rprimd(idir,2)*pel(2) + &
732 &   rprimd(idir,3)*pel(3)
733    pion_cart(idir) = rprimd(idir,1)*pion(1) + rprimd(idir,2)*pion(2) + &
734 &   rprimd(idir,3)*pion(3)
735    ptot_cart(idir) = rprimd(idir,1)*ptot(1) + rprimd(idir,2)*ptot(2) + &
736 &   rprimd(idir,3)*ptot(3)
737  end do
738 
739 !Divide by the unit cell volume
740  pel_cart(:)  = pel_cart(:)/ucvol
741  pion_cart(:) = pion_cart(:)/ucvol
742  ptot_cart(:)  = ptot_cart(:)/ucvol
743 !pelev is either zero (in NCPP case), or possibly non-zero (in PAW case)
744 !however, in the PAW case, it is already implicitly included in the computation
745 !of pel and should not be added in here. Below in the PAW case we report its
746 !value to the user, which is helpful for comparison to the USPP theory.
747 !13 June 2012 J Zwanziger
748 !note that pelev was AlREADY in cartesian frame.
749 !ptot_cart(:)  = (ptot_cart(:)+pelev(:))/ucvol
750 
751 !Write the results to unit_out (if /= 0)
752 !Use the coordinates specified by the value polunit
753 
754 !Atomic units
755 
756  if (((polunit == 1).or.(polunit == 3)).and.(unit_out /= 0)) then
757 
758 !  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,&
759 !  &   ' Polarization in cartesian coordinates (a.u.):',ch10,&
760 !  &   ' (the sum of the electronic and ionic Berry phase',&
761 !  &   ' has been fold into [-1, 1])',ch10,&
762 !  &   '     Electronic berry phase:       ', (pel_cart(idir), idir = 1, 3), ch10,&
763 !  &   '     Expectation value (PAW only): ', (pelev(idir)/ucvol, idir = 1, 3), ch10,&
764 !  &   '     Ionic:                        ', (pion_cart(idir), idir = 1, 3), ch10, &
765 !  &   '     Total:                        ', (ptot_cart(idir), idir = 1, 3)
766 !  call wrtout(unit_out,message,'COLL')
767    write(message,'(7(a),3(e16.9,2x),a,a,3(e16.9,2x))')ch10,&
768 &   ' Polarization in cartesian coordinates (a.u.):',ch10,&
769 &   ' (the sum of the electronic and ionic Berry phase',&
770 &   ' has been folded into [-1, 1])',ch10,&
771 &   '     Electronic berry phase:       ', (pel_cart(idir), idir = 1, 3)
772    call wrtout(unit_out,message,'COLL')
773    if(usepaw==1) then
774      write(message,'(a,3(e16.9,2x))')&
775 &     '     ...includes PAW on-site term: ', (pelev(idir)/ucvol, idir = 1, 3)
776      call wrtout(unit_out,message,'COLL')
777    end if
778    write(message,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')&
779 &   '     Ionic:                        ', (pion_cart(idir), idir = 1, 3), ch10, &
780 &   '     Total:                        ', (ptot_cart(idir), idir = 1, 3)
781    call wrtout(unit_out,message,'COLL')
782 
783  end if
784 
785 !MKS units
786 
787  if (((polunit == 2).or.(polunit == 3)).and.(unit_out /= 0)) then
788 
789    pel_mks(:)  = pel_cart(:)*(e_Cb)/(Bohr_Ang*1d-10)**2
790    pion_mks(:) = pion_cart(:)*(e_Cb)/(Bohr_Ang*1d-10)**2
791    pelev_mks(:) = pelev(:)/ucvol*(e_Cb)/(Bohr_Ang*1d-10)**2
792    ptot_mks(:) = (ptot_cart(:))*(e_Cb)/(Bohr_Ang*1d-10)**2
793 
794    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,&
795 &   ' Polarization in cartesian coordinates (C/m^2):',ch10,&
796 &   ' (the sum of the electronic and ionic Berry phase',&
797 &   ' has been folded into [-1, 1])',ch10,&
798 &   '     Electronic berry phase:       ', (pel_mks(idir), idir = 1, 3)
799    call wrtout(unit_out,message,'COLL')
800    if(usepaw==1) then
801      write(message,'(a,3(e16.9,2x))')&
802 &     '     ...includes PAW on-site term: ', (pelev_mks(idir), idir = 1, 3)
803      call wrtout(unit_out,message,'COLL')
804    end if
805    write(message,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')&
806 &   '     Ionic:                        ', (pion_mks(idir), idir = 1, 3), ch10, &
807 &   '     Total:                        ', (ptot_mks(idir), idir = 1, 3)
808    call wrtout(unit_out,message,'COLL')
809 
810  end if
811 
812 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.

PARENTS

      berryphase_new,cgwf,chern_number,getcgqphase,make_grad_berry

CHILDREN

      dzgedi,dzgefa,overlap_g

SOURCE

142 subroutine smatrix(cg,cgq,cg1_k,ddkflag,dtm_k,icg,icg1,itrs,job,maxbd,&
143 &  mcg_k,mcg_q,mcg1_k,minbd,mpw,mband_occ,nband_occ,npw_k1,npw_k2,nspinor,&
144 &  pwind_k,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_k,smat_k_paw,usepaw)
145 
146 
147 !This section has been created automatically by the script Abilint (TD).
148 !Do not modify the following lines by hand.
149 #undef ABI_FUNC
150 #define ABI_FUNC 'smatrix'
151 !End of the abilint section
152 
153  implicit none
154 
155 !Arguments ------------------------------------
156 !scalars
157  integer,intent(in) :: ddkflag,icg,icg1,itrs,job,maxbd,mcg1_k,mcg_k,mcg_q
158  integer,intent(in) :: minbd,mpw,mband_occ,npw_k1,npw_k2,nspinor,shiftbd
159  integer,intent(in) :: nband_occ
160  integer,intent(in) :: usepaw
161 !arrays
162  integer,intent(in) :: pwind_k(mpw)
163  integer,intent(inout) :: sflag_k(mband_occ)
164  real(dp),intent(in) :: cg(2,mcg_k),cgq(2,mcg_q),pwnsfac_k(4,mpw)
165  real(dp),intent(in) :: smat_k_paw(2,usepaw*mband_occ,usepaw*mband_occ)
166  real(dp),intent(inout) :: smat_k(2,mband_occ,mband_occ)
167  real(dp),intent(out) :: cg1_k(2,mcg1_k),dtm_k(2)
168  real(dp),intent(out) :: smat_inv(2,mband_occ,mband_occ)
169 
170 !Local variables -------------------------
171 !scalars
172  integer :: count,dzgedi_job,iband,info,ipw,ispinor,jband,jband1,jpw,pwmax,pwmin
173  integer :: spnshft_k1, spnshft_k2
174  real(dp) :: doti,dotr,fac,wfi,wfr
175  character(len=500) :: message
176 !arrays
177  integer,allocatable :: ipvt(:)
178 ! integer,allocatable :: my_pwind_k(:) ! used in debugging below
179  real(dp) :: det(2,2)
180  real(dp),allocatable :: vect1(:,:),vect2(:,:),zgwork(:,:)
181 
182 ! ***********************************************************************
183 
184 !DEBUG
185 !write(std_out,*)'smatrix : enter'
186 !write(std_out,*)'sflag_k = ',sflag_k
187 !write(std_out,'(a,4i4)' )'job, ddkflag, shiftbd, itrs = ',job,ddkflag,shiftbd,itrs
188 !write(std_out,'(a,2i6)')' JWZ smatrix.F90 debug : npw_k1, npw_k2 ',npw_k1,npw_k2
189 !stop
190 !ENDDEBUG
191 
192  ABI_ALLOCATE(ipvt,(nband_occ))
193  ABI_ALLOCATE(zgwork,(2,nband_occ))
194  ABI_ALLOCATE(vect1,(2,0:mpw*nspinor))
195  ABI_ALLOCATE(vect2,(2,0:mpw*nspinor))
196  vect1(:,0) = zero ; vect2(:,0) = zero
197 
198 !Check if the values of ddkflag and job are compatible
199 
200  if ((job /= 0).and.(job /= 1).and.(job /= 10).and.(job /= 11).and.(job/=20).and.(job/=21)) then
201    write(message,'(a,i3,a,a)')&
202 &   ' job is equal to ',job,ch10,&
203 &   ' while only the values job = 0, 1, 10, 11, 20, or 21 are allowed.'
204    MSG_ERROR(message)
205  end if
206 
207  if (ddkflag == 1) then
208    if ((job/=1).and.(job/=11)) then
209      write(message,'(a,i0,a,a)')&
210 &     ' job is equal to ',job,ch10,&
211 &     ' while ddkflag = 1. This is not allowed.'
212      MSG_ERROR(message)
213    end if
214  end if
215 
216 !Check the values of sflag_k
217  do iband=1,nband_occ
218    if (sflag_k(iband)/=0 .and. sflag_k(iband)/=1)then
219      write(message,'(3a,i4,a,i4)')&
220 &     '  The content of sflag_k must be 0 or 1.',ch10,&
221 &     '  However, for iband=',iband,', sflag_k(iband)=',sflag_k(iband)
222      MSG_ERROR(message)
223    end if
224  end do
225 
226 !Check if shiftbd is consistent with sflag_k
227  if (shiftbd == 0) then
228    count = 0
229    do iband = 1, nband_occ
230      if (sflag_k(iband) == 0) count = count + 1
231    end do
232    if (count > 1) then
233      message = 'in case shiftbd = 0, only 1 element of sflag can be 0'
234      MSG_ERROR(message)
235    end if
236  end if
237 
238 !Update the lines of the overlap matrix for which sflag = 0
239 !MVeithen: because of sflag, it is more efficient to perform
240 !the loop over jband inside the loop over iband
241 
242 !DEBUG
243 !write(std_out,*)' smatrix : smat_k(1,1,1)=',smat_k(1,1,1)
244 !write(std_out,*)' smatrix : sflag_k=',sflag_k
245 !ENDDEBUG
246 
247 !!
248 !! debugging based on norm of k1 vector
249 !
250 !ABI_ALLOCATE(my_pwind_k,(mpw))
251 !!
252 !do iband = 1, nband_occ
253 !!
254 !pwmin = (iband-1)*npw_k1*nspinor*shiftbd
255 !pwmax = pwmin + npw_k1*nspinor
256 !!
257 !!!    Multiply the bra wave function by the phase factor
258 !if (itrs==1.or.itrs==11) then  ! take complex conjugate of bra
259 !do ispinor = 1, nspinor
260 !spnshft_k1=(ispinor-1)*npw_k1
261 !do ipw = 1,npw_k1
262 !vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
263 !&           +cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw)
264 !vect1(2,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw) &
265 !&           -cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
266 !end do
267 !end do
268 !else
269 !do ispinor=1, nspinor
270 !spnshft_k1=(ispinor-1)*npw_k1
271 !do ipw = 1,npw_k1
272 !vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
273 !&           -cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw)
274 !vect1(2,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw) &
275 !&           +cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
276 !end do
277 !end do
278 !end if
279 !
280 !!
281 !if (npw_k1*nspinor < mpw*nspinor) vect1(:,npw_k1*nspinor+1:mpw*nspinor) = zero
282 !
283 !do jband = 1, nband_occ
284 !
285 !pwmin = (jband-1)*npw_k1*nspinor
286 !pwmax = pwmin + npw_k1*nspinor
287 !
288 !if (itrs==10.or.itrs==11) then ! take complex conjugate of ket
289 !do ispinor=1, nspinor
290 !spnshft_k2=(ispinor-1)*npw_k1
291 !do ipw = 1, npw_k1
292 !my_pwind_k(ipw) = ipw
293 !vect2(1,spnshft_k2+ipw) = cg(1,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(1,ipw) &
294 !&             +cg(2,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(2,ipw)
295 !vect2(2,spnshft_k2+ipw) = cg(1,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(2,ipw) &
296 !&             -cg(2,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(1,ipw)
297 !end do
298 !end do
299 !else
300 !do ispinor=1, nspinor
301 !spnshft_k2=(ispinor-1)*npw_k1
302 !do ipw = 1, npw_k1
303 !my_pwind_k(ipw) = ipw
304 !vect2(1,spnshft_k2+ipw) = cg(1,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(1,ipw) &
305 !&             -cg(2,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(2,ipw)
306 !vect2(2,spnshft_k2+ipw) = cg(1,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(2,ipw) &
307 !&             +cg(2,icg+spnshft_k2+ipw+pwmin)*pwnsfac_k(1,ipw)
308 !end do
309 !end do
310 !end if
311 !!
312 !if (npw_k1*nspinor < mpw*nspinor) vect2(:,npw_k1*nspinor+1:mpw*nspinor) = zero
313 !!
314 !call overlap_g(doti,dotr,mpw,npw_k1,npw_k1,nspinor,my_pwind_k,vect1,vect2)
315 !!
316 !smat_k(1,iband,jband) = dotr
317 !smat_k(2,iband,jband) = doti
318 !!
319 !if (usepaw == 1) then
320 !smat_k(1,iband,jband) = smat_k(1,iband,jband)+smat_k_paw(1,iband,jband)
321 !smat_k(2,iband,jband) = smat_k(2,iband,jband)+smat_k_paw(2,iband,jband)
322 !end if
323 !!
324 !if( (smat_k(1,iband,jband)**2+smat_k(2,iband,jband)**2)>tol8) then
325 !write(std_out,'(a,i2,a,i2,a,es16.8,a,es16.8)')' JWZ Debug: <',iband,'|',jband,'> = ',&
326 !&                     smat_k(1,iband,jband),' + i ',smat_k(2,iband,jband)
327 !end if
328 !!
329 !end do   ! jband
330 !!
331 !end do   ! iband
332 !!
333 !ABI_DEALLOCATE(my_pwind_k)
334 !!
335  do iband = 1, nband_occ
336 
337    if (sflag_k(iband) == 0) then
338 
339      pwmin = (iband-1)*npw_k1*nspinor*shiftbd
340      pwmax = pwmin + npw_k1*nspinor
341 !
342 !    old version  (*** multiply by nspinor missing??? ***)
343 !    vect1(:,1:npw_k1) = cg(:,icg + 1 + pwmin:icg + pwmax)
344 !
345 
346 !    Multiply the bra wave function by the phase factor
347      if (itrs==1.or.itrs==11) then  ! take complex conjugate of bra
348        do ispinor = 1, nspinor
349          spnshft_k1=(ispinor-1)*npw_k1
350          do ipw = 1,npw_k1
351            vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
352 &           +cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw)
353            vect1(2,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw) &
354 &           -cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
355          end do
356        end do
357      else
358        do ispinor=1, nspinor
359          spnshft_k1=(ispinor-1)*npw_k1
360          do ipw = 1,npw_k1
361            vect1(1,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw) &
362 &           -cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw)
363            vect1(2,spnshft_k1+ipw) = cg(1,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(2,ipw) &
364 &           +cg(2,icg+spnshft_k1+ipw+pwmin)*pwnsfac_k(1,ipw)
365          end do
366        end do
367      end if
368 
369 !
370      if (npw_k1*nspinor < mpw*nspinor) vect1(:,npw_k1*nspinor+1:mpw*nspinor) = zero
371 
372      do jband = 1, nband_occ
373 
374        pwmin = (jband-1)*npw_k2*nspinor
375        pwmax = pwmin + npw_k2*nspinor
376 
377        if (itrs==10.or.itrs==11) then ! take complex conjugate of ket
378          do ispinor=1, nspinor
379            spnshft_k2=(ispinor-1)*npw_k2
380            do ipw = 1, npw_k2
381              vect2(1,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw) &
382 &             +cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(4,ipw)
383              vect2(2,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(4,ipw) &
384 &             -cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw)
385            end do
386          end do
387        else
388          do ispinor=1, nspinor
389            spnshft_k2=(ispinor-1)*npw_k2
390            do ipw = 1, npw_k2
391              vect2(1,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw) &
392 &             -cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(4,ipw)
393              vect2(2,spnshft_k2+ipw) = cgq(1,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(4,ipw) &
394 &             +cgq(2,icg1+spnshft_k2+ipw+pwmin)*pwnsfac_k(3,ipw)
395            end do
396          end do
397        end if
398 
399        if (npw_k2*nspinor < mpw*nspinor) vect2(:,npw_k2*nspinor+1:mpw*nspinor) = zero
400 
401 !      DEBUG
402 !      if(iband==1 .and. jband==1 .and. itrs==0 .and. npw_k1==68 .and. npw_k2==74)then
403 !      write(std_out,'(a)' )' smatrix : ii,vect1,cg,pwnsfac='
404 !      do ii=1,npw_k1
405 !      write(std_out,'(i4,6es16.6)' )ii,vect1(:,ii),cg(:,icg+ii+pwmin),pwnsfac_k(1:2,ii)
406 !      end do
407 !      do ii=1,npw_k2
408 !      write(std_out,'(i4,6es16.6)' )ii,vect2(:,ii),cg(:,icg1+ii+pwmin),pwnsfac_k(3:4,ii)
409 !      end do
410 !      end if
411 !      ENDDEBUG
412 
413        call overlap_g(doti,dotr,mpw,npw_k1,npw_k2,nspinor,pwind_k,vect1,vect2)
414 
415        smat_k(1,iband,jband) = dotr
416        smat_k(2,iband,jband) = doti
417 
418        if (usepaw == 1) then
419          smat_k(1,iband,jband) = smat_k(1,iband,jband)+smat_k_paw(1,iband,jband)
420          smat_k(2,iband,jband) = smat_k(2,iband,jband)+smat_k_paw(2,iband,jband)
421        end if
422 
423 !      DEBUG
424 !      if(iband==1 .and. jband==1)then
425 !      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
426 !      end if
427 !      ENDDEBUG
428 
429      end do   ! jband
430 
431    end if    ! sflag_k(iband) == 0
432 
433  end do   ! iband
434 
435 !DEBUG
436 !do iband=1,nband_occ
437 !do jband=1,nband_occ
438 !write(std_out,'(a,2i4,2e20.10)') 'smat',iband,jband,smat_k(1,iband,jband),smat_k(2,iband,jband)
439 !end do
440 !end do
441 !write(std_out,*)' smatrix : smat_k(1,1,1)=',smat_k(1,1,1)
442 !ENDDEBUG
443 
444 !Update sflag_k
445  sflag_k(:) = 1
446 
447 !Depending on the value of job, compute the determinant of the
448 !overlap matrix, its inverse or the product of the inverse
449 !overlap matrix with the WF at k.
450 
451  if ((job==1).or.(job==10).or.(job==11).or.(job==21)) then
452 
453    smat_inv(:,:,:) = smat_k(:,:,:)
454 
455 !  DEBUG
456 !  write(std_out,*)' smatrix : smat_inv=',smat_inv
457 !  ENDDEBUG
458 
459    dzgedi_job=job; if(job==21) dzgedi_job=1
460 ! TODO: should this be over nband_occ(isppol)?
461    call dzgefa(smat_inv,mband_occ,nband_occ,ipvt,info)
462    call dzgedi(smat_inv,mband_occ,nband_occ,ipvt,det,zgwork,dzgedi_job)
463 
464 !  DEBUG
465 !  write(std_out,*)' smatrix : det=',det
466 !  ENDDEBUG
467 
468 !  Compute the determinant of the overlap matrix
469    dtm_k(:) = zero
470    if (job==10 .or. job==11) then
471      fac = exp(log(10._dp)*det(1,2))
472      dtm_k(1) = fac*(det(1,1)*cos(log(10._dp)*det(2,2)) - &
473 &     det(2,1)*sin(log(10._dp)*det(2,2)))
474      dtm_k(2) = fac*(det(1,1)*sin(log(10._dp)*det(2,2)) + &
475 &     det(2,1)*cos(log(10._dp)*det(2,2)))
476    end if
477 
478 !  Compute the product of the inverse overlap matrix with the WF
479 
480    if (ddkflag == 1) then
481 
482      cg1_k(:,:) = zero
483      jband1 = 0
484 
485      if (itrs == 10 .or. itrs == 11) then
486 
487        do jband = minbd, maxbd
488          jband1 = jband1 + 1
489          do iband = 1, nband_occ
490 
491            do ispinor = 1, nspinor
492              spnshft_k1 = (ispinor-1)*npw_k1
493              spnshft_k2 = (ispinor-1)*npw_k2
494              do ipw = 1, npw_k1
495 
496                jpw = pwind_k(ipw)
497 
498                if (jpw > 0) then
499 
500                  wfr = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
501 &                 -cgq(2,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)
502                  wfi = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)&
503 &                 +cgq(2,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
504 
505                  cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = &
506 &                 cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) + &
507 &                 smat_inv(1,iband,jband)*wfr + smat_inv(2,iband,jband)*wfi
508 
509                  cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = &
510 &                 cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) - &
511 &                 smat_inv(1,iband,jband)*wfi + smat_inv(2,iband,jband)*wfr
512 
513                end if
514 
515              end do ! end loop over npw_k1
516            end do ! end loop over nspinor
517 
518          end do
519        end do
520 
521      else
522 
523        do jband = minbd, maxbd
524          jband1 = jband1 + 1
525          do iband = 1, nband_occ
526 
527            do ispinor = 1, nspinor
528              spnshft_k1 = (ispinor-1)*npw_k1
529              spnshft_k2 = (ispinor-1)*npw_k2
530              do ipw = 1, npw_k1
531 
532                jpw = pwind_k(ipw)
533 
534                if (jpw > 0) then
535 
536                  wfr = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
537 &                 -cgq(2,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)
538                  wfi = cgq(1,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)&
539 &                 +cgq(2,icg1+(iband-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
540 
541                  cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = &
542 &                 cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) + &
543 &                 smat_inv(1,iband,jband)*wfr - smat_inv(2,iband,jband)*wfi
544 
545                  cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = &
546 &                 cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) + &
547 &                 smat_inv(1,iband,jband)*wfi + smat_inv(2,iband,jband)*wfr
548 
549                end if
550 
551              end do ! end loop over npw_k1
552            end do ! end loop over nspinor
553 
554          end do
555        end do
556 
557      end if     ! itrs
558 
559    end if
560 
561  end if         !
562 
563  if(job == 20 .or. job == 21) then ! special case transfering cgq to cg1_k without use of S^{-1}, used in
564 !  magnetic field case
565 
566    cg1_k(:,:) = zero
567    jband1 = 0
568 
569    if (itrs == 10 .or. itrs == 11) then
570 
571      do jband = minbd, maxbd
572        jband1 = jband1 + 1
573        do ispinor = 1, nspinor
574          spnshft_k1 = (ispinor-1)*npw_k1
575          spnshft_k2 = (ispinor-1)*npw_k2
576          do ipw = 1, npw_k1
577            jpw = pwind_k(ipw)
578 
579            if (jpw > 0) then
580              wfr = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
581 &             -cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)
582              wfi = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)&
583 &             +cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
584 
585              cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfr
586              cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfi
587 
588            end if
589 
590          end do ! end loop over npw_k1
591        end do ! end loop over nspinor
592 
593      end do
594 
595    else
596 
597      do jband = minbd, maxbd
598        jband1 = jband1 + 1
599        do ispinor = 1, nspinor
600          spnshft_k1 = (ispinor-1)*npw_k1
601          spnshft_k2 = (ispinor-1)*npw_k2
602          do ipw = 1, npw_k1
603            jpw = pwind_k(ipw)
604 
605            if (jpw > 0) then
606 
607              wfr = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)&
608 &             -cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)
609              wfi = cgq(1,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(4,jpw)&
610              +cgq(2,icg1+(jband1-1)*npw_k2*nspinor+spnshft_k2+jpw)*pwnsfac_k(3,jpw)
611 
612              cg1_k(1,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfr
613              cg1_k(2,(jband1-1)*npw_k1*nspinor + spnshft_k1 + ipw) = wfi
614 
615            end if
616 
617          end do ! end loop over npw_k1
618        end do ! end loop over nspinor
619 
620      end do
621 
622    end if     ! itrs
623 
624  end if ! end job == 20 .or. job == 21 case
625 
626  ABI_DEALLOCATE(ipvt)
627  ABI_DEALLOCATE(zgwork)
628  ABI_DEALLOCATE(vect1)
629  ABI_DEALLOCATE(vect2)
630 
631 !DEBUG
632 !write(std_out,*)' dtm_k=',dtm_k(:)
633 !write(std_out,*)' smatrix : exit '
634 !ENDDEBUG
635 
636 end subroutine smatrix