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