TABLE OF CONTENTS
ABINIT/subdiago [ Functions ]
NAME
subdiago
FUNCTION
This routine diagonalizes the Hamiltonian in the eigenfunction subspace
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MT) 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
icg=shift to be applied on the location of data in the array cg igsc=shift to be applied on the location of data in the array gsc istwf_k=input parameter that describes the storage of wfs mcg=second dimension of the cg array mgsc=second dimension of the gsc array nband_k=number of bands at this k point for that spin polarization npw_k=number of plane waves at this k point nspinor=number of spinorial components of the wavefunctions (on current proc) subham(nband_k*(nband_k+1))=Hamiltonian expressed in the WFs subspace subovl(nband_k*(nband_k+1)*use_subovl)=overlap matrix expressed in the WFs subspace use_subovl=1 if the overlap matrix is not identity in WFs subspace usepaw= 0 for non paw calculation; =1 for paw calculation me_g0=1 if this processors has G=0, 0 otherwise.
OUTPUT
eig_k(nband_k)=array for holding eigenvalues (hartree) evec(2*nband_k,nband_k)=array for holding eigenvectors
SIDE EFFECTS
cg(2,mcg)=wavefunctions gsc(2,mgsc)=<g|S|c> matrix elements (S=overlap)
PARENTS
rayleigh_ritz,vtowfk
CHILDREN
abi_xcopy,abi_xgemm,abi_xhpev,abi_xhpgv,cg_normev,hermit
SOURCE
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 54 subroutine subdiago(cg,eig_k,evec,gsc,icg,igsc,istwf_k,& 55 & mcg,mgsc,nband_k,npw_k,nspinor,paral_kgb,& 56 & subham,subovl,use_subovl,usepaw,me_g0) 57 58 use defs_basis 59 use m_profiling_abi 60 use m_errors 61 use m_cgtools 62 use m_linalg_interfaces 63 use m_abi_linalg 64 use m_xmpi 65 66 use m_numeric_tools, only : hermit 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'subdiago' 72 !End of the abilint section 73 74 implicit none 75 76 !Arguments ------------------------------------ 77 integer,intent(in) :: icg,igsc,istwf_k,mcg,mgsc,nband_k,npw_k,me_g0 78 integer,intent(in) :: nspinor,paral_kgb,use_subovl,usepaw 79 real(dp),intent(inout) :: subham(nband_k*(nband_k+1)),subovl(nband_k*(nband_k+1)*use_subovl) 80 real(dp),intent(out) :: eig_k(nband_k),evec(2*nband_k,nband_k) 81 real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc) 82 83 !Local variables------------------------------- 84 integer :: iband,ii,ierr,rvectsize,vectsize,use_slk 85 character(len=500) :: message 86 ! real(dp) :: tsec(2) 87 real(dp),allocatable :: evec_tmp(:,:),subovl_tmp(:),subham_tmp(:) 88 real(dp),allocatable :: work(:,:) 89 real(dp),allocatable :: blockvectora(:,:),blockvectorb(:,:),blockvectorc(:,:) 90 91 ! ********************************************************************* 92 93 if (paral_kgb<0) then 94 MSG_BUG('paral_kgb should be positive ') 95 end if 96 97 ! 1 if Scalapack version is used. 98 use_slk = paral_kgb 99 100 rvectsize=npw_k*nspinor 101 vectsize=2*rvectsize;if (me_g0==1) vectsize=vectsize-1 102 103 !Impose Hermiticity on diagonal elements of subham (and subovl, if needed) 104 ! MG FIXME: In these two calls we are aliasing the args 105 call hermit(subham,subham,ierr,nband_k) 106 if (use_subovl==1) then 107 call hermit(subovl,subovl,ierr,nband_k) 108 end if 109 110 !Diagonalize the Hamitonian matrix 111 if(istwf_k==2) then 112 ABI_ALLOCATE(evec_tmp,(nband_k,nband_k)) 113 ABI_ALLOCATE(subham_tmp,(nband_k*(nband_k+1)/2)) 114 subham_tmp=subham(1:nband_k*(nband_k+1):2) 115 evec_tmp=zero 116 if (use_subovl==1) then 117 ABI_ALLOCATE(subovl_tmp,(nband_k*(nband_k+1)/2)) 118 subovl_tmp=subovl(1:nband_k*(nband_k+1):2) 119 ! TO DO: Not sure this one has been fully tested 120 call abi_xhpgv(1,'V','U',nband_k, subham_tmp,subovl_tmp, eig_k,evec_tmp,istwf_k=istwf_k,use_slk=use_slk) 121 ABI_DEALLOCATE(subovl_tmp) 122 else 123 call abi_xhpev('V','U',nband_k,subham_tmp,eig_k,evec_tmp,istwf_k=istwf_k,use_slk=use_slk) 124 end if 125 evec(:,:)=zero;evec(1:2*nband_k:2,:) =evec_tmp 126 ABI_DEALLOCATE(evec_tmp) 127 ABI_DEALLOCATE(subham_tmp) 128 else 129 if (use_subovl==1) then 130 call abi_xhpgv(1,'V','U',nband_k,subham,subovl,eig_k,evec,istwf_k=istwf_k,use_slk=use_slk) 131 else 132 call abi_xhpev('V','U',nband_k,subham,eig_k,evec,istwf_k=istwf_k,use_slk=use_slk) 133 end if 134 end if 135 136 !Normalize each eigenvector and set phase: 137 !The problem with minus/plus signs might be present also if .not. use_subovl 138 !if(use_subovl == 0) then 139 call cg_normev(evec,nband_k,nband_k) 140 !end if 141 142 if(istwf_k==2)then 143 do iband=1,nband_k 144 do ii=1,nband_k 145 if(abs(evec(2*ii,iband))>1.0d-10)then 146 write(message,'(3a,2i0,2es16.6,a,a)')ch10,& 147 & ' subdiago: For istwf_k=2, observed the following element of evec :',ch10,& 148 & iband,ii,evec(2*ii-1,iband),evec(2*ii,iband),ch10,' with a non-negligible imaginary part.' 149 MSG_BUG(message) 150 end if 151 end do 152 end do 153 end if 154 155 !===================================================== 156 !Carry out rotation of bands C(G,n) according to evecs 157 ! ZGEMM if istwfk==1, DGEMM if istwfk==2 158 !===================================================== 159 if (istwf_k==2) then 160 161 ABI_STAT_ALLOCATE(blockvectora,(vectsize,nband_k), ierr) 162 ABI_CHECK(ierr==0, "out-of-memory in blockvectora") 163 ABI_STAT_ALLOCATE(blockvectorb,(nband_k,nband_k), ierr) 164 ABI_CHECK(ierr==0, "out-of-memory in blockvectorb") 165 ABI_STAT_ALLOCATE(blockvectorc,(vectsize,nband_k), ierr) 166 ABI_CHECK(ierr==0, "out-of-memory in blockvectorc") 167 168 do iband=1,nband_k 169 if (me_g0 == 1) then 170 call abi_xcopy(1,cg(1,cgindex_subd(iband)),1,blockvectora(1,iband),1) 171 call abi_xcopy(rvectsize-1,cg(1,cgindex_subd(iband)+1),2,blockvectora(2,iband),1) 172 call abi_xcopy(rvectsize-1,cg(2,cgindex_subd(iband)+1),2,blockvectora(rvectsize+1,iband),1) 173 else 174 call abi_xcopy(rvectsize,cg(1,cgindex_subd(iband)),2,blockvectora(1,iband),1) 175 call abi_xcopy(rvectsize,cg(2,cgindex_subd(iband)),2,blockvectora(rvectsize+1,iband),1) 176 end if 177 call abi_xcopy(nband_k,evec(2*iband-1,1),2*nband_k,blockvectorb(iband,1),nband_k) 178 end do 179 180 call abi_xgemm('N','N',vectsize,nband_k,nband_k,& 181 & cone,blockvectora,vectsize,blockvectorb,nband_k,czero,blockvectorc,vectsize) 182 183 do iband=1,nband_k 184 if (me_g0 == 1) then 185 call abi_xcopy(1,blockvectorc(1,iband),1,cg(1,cgindex_subd(iband)),1) 186 call abi_xcopy(rvectsize-1,blockvectorc(2,iband),1,cg(1,cgindex_subd(iband)+1),2) 187 call abi_xcopy(rvectsize-1,blockvectorc(rvectsize+1,iband),1,cg(2,cgindex_subd(iband)+1),2) 188 else 189 call abi_xcopy(rvectsize,blockvectorc(1,iband),1,cg(1,cgindex_subd(iband)),2) 190 call abi_xcopy(rvectsize,blockvectorc(rvectsize+1,iband),1,cg(2,cgindex_subd(iband)),2) 191 end if 192 end do 193 194 ! If paw, musb also rotate S.C(G,n): 195 if (usepaw==1) then 196 197 do iband=1,nband_k 198 if (me_g0 == 1) then 199 call abi_xcopy(1,gsc(1,gscindex_subd(iband)),1,blockvectora(1,iband),1) 200 call abi_xcopy(rvectsize-1,gsc(1,gscindex_subd(iband)+1),2,blockvectora(2,iband),1) 201 call abi_xcopy(rvectsize-1,gsc(2,gscindex_subd(iband)+1),2,blockvectora(rvectsize+1,iband),1) 202 else 203 call abi_xcopy(rvectsize ,gsc(1,gscindex_subd(iband)),2,blockvectora(1,iband),1) 204 call abi_xcopy(rvectsize ,gsc(2,gscindex_subd(iband)),2,blockvectora(rvectsize+1,iband),1) 205 end if 206 call abi_xcopy(nband_k,evec(2*iband-1,1),2*nband_k,blockvectorb(iband,1),nband_k) 207 end do 208 209 call abi_xgemm('N','N',vectsize,nband_k,nband_k,& 210 & cone,blockvectora,vectsize,blockvectorb,nband_k,czero,blockvectorc,vectsize) 211 212 do iband=1,nband_k 213 if (me_g0 == 1) then 214 call abi_xcopy(1,blockvectorc(1,iband),1,gsc(1,gscindex_subd(iband)),1) 215 call abi_xcopy(rvectsize-1,blockvectorc(2,iband),1,gsc(1,gscindex_subd(iband)+1),2) 216 call abi_xcopy(rvectsize-1,blockvectorc(rvectsize+1,iband),1,gsc(2,gscindex_subd(iband)+1),2) 217 else 218 call abi_xcopy(rvectsize,blockvectorc(1,iband),1,gsc(1,gscindex_subd(iband)),2) 219 call abi_xcopy(rvectsize,blockvectorc(rvectsize+1,iband),1,gsc(2,gscindex_subd(iband)),2) 220 end if 221 end do 222 223 end if 224 225 ABI_DEALLOCATE(blockvectora) 226 ABI_DEALLOCATE(blockvectorb) 227 ABI_DEALLOCATE(blockvectorc) 228 229 else 230 231 ABI_STAT_ALLOCATE(work,(2,npw_k*nspinor*nband_k), ierr) 232 ABI_CHECK(ierr==0, "out-of-memory in work") 233 234 ! MG: Do not remove this initialization. 235 ! telast_06 stops in fxphase on inca_debug and little_buda (very very strange, due to atlas?) 236 work=zero 237 238 call abi_xgemm('N','N',npw_k*nspinor,nband_k,nband_k,cone, & 239 & cg(:,icg+1:npw_k*nspinor*nband_k+icg),npw_k*nspinor, & 240 & evec,nband_k,czero,work,npw_k*nspinor,x_cplx=2) 241 242 call abi_xcopy(npw_k*nspinor*nband_k,work(1,1),1,cg(1,1+icg),1,x_cplx=2) 243 244 ! If paw, must also rotate S.C(G,n): 245 if (usepaw==1) then 246 call abi_xgemm('N','N',npw_k*nspinor,nband_k,nband_k,cone, & 247 & gsc(:,1+igsc:npw_k*nspinor*nband_k+igsc),npw_k*nspinor, & 248 & evec,nband_k,czero,work,npw_k*nspinor,x_cplx=2) 249 call abi_xcopy(npw_k*nspinor*nband_k, work(1,1),1,gsc(1,1+igsc),1,x_cplx=2) 250 end if 251 252 ABI_DEALLOCATE(work) 253 end if 254 255 contains 256 257 function cgindex_subd(iband) 258 259 260 !This section has been created automatically by the script Abilint (TD). 261 !Do not modify the following lines by hand. 262 #undef ABI_FUNC 263 #define ABI_FUNC 'cgindex_subd' 264 !End of the abilint section 265 266 integer :: iband,cgindex_subd 267 cgindex_subd=npw_k*nspinor*(iband-1)+icg+1 268 end function cgindex_subd 269 function gscindex_subd(iband) 270 271 272 !This section has been created automatically by the script Abilint (TD). 273 !Do not modify the following lines by hand. 274 #undef ABI_FUNC 275 #define ABI_FUNC 'gscindex_subd' 276 !End of the abilint section 277 278 integer :: iband,gscindex_subd 279 gscindex_subd=npw_k*nspinor*(iband-1)+igsc+1 280 end function gscindex_subd 281 282 end subroutine subdiago