TABLE OF CONTENTS
ABINIT/smatrix [ 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