TABLE OF CONTENTS


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.

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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

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