TABLE OF CONTENTS
ABINIT/mlwfovlp_ylmfac [ Functions ]
NAME
mlwfovlp_ylmfac
FUNCTION
Routine that produces a factor by which the initial guess of functions will be multiplied for the Wannier90 interface. It is just used if there are rotations, or if the functions required are linear combinations of the ylm real functions. Example, For a function G(r)= 1/2 s + 1/3 px - 1/2 pz it would produce a matrix of the following form: [1/2,-1/2,1/3,0,0...0] The real spherical harmonics are given as factors of complex spherical harmonics The real spherical harmonics are given in table 3.1 of Wannier90 user guide.
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (T. Rangel, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
lmax= maximum l value for spherical harmonics lmax2=number of ylm functions mband=maximum number of bands nwan = number of wannier functions proj_l(mband)= angular part of the projection function (quantum number l) proj_m(mband)= angular part of the projection function (quantum number m) proj_x(3,mband)= x axis for the projection. proj_z(3,mband)= z axis for the projection.
OUTPUT
ylmc_fac(lmax2,nwan)=matrix containig a factor for ylm hybrid orbitals
SIDE EFFECTS
(only writing, printing)
NOTES
PARENTS
mlwfovlp_proj
CHILDREN
rotmat,ylm_cmplx,zgesv
SOURCE
52 #if defined HAVE_CONFIG_H 53 #include "config.h" 54 #endif 55 56 #include "abi_common.h" 57 58 59 subroutine mlwfovlp_ylmfac(ylmc_fac,lmax,lmax2,mband,nwan,proj_l,proj_m,proj_x,proj_z) 60 61 use defs_basis 62 use m_profiling_abi 63 use m_errors 64 65 use m_geometry, only : rotmat 66 use m_numeric_tools, only : uniformrandom 67 use m_paw_sphharm, only : ylm_cmplx 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'mlwfovlp_ylmfac' 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 integer, intent(in):: lmax,lmax2,nwan,mband 79 ! arrays 80 integer,intent(in) :: proj_l(mband),proj_m(mband) 81 real(dp),intent(in) :: proj_x(3,mband),proj_z(3,mband) 82 complex(dp),intent(out)::ylmc_fac(lmax2,nwan) 83 ! 84 !Local variables------------------------------- 85 ! 86 integer :: orb_idx(16)=(/1,3,4,2,7,8,6,9,5,13,14,12,15,11,16,10/) !Tab3.1 Wannier90 user guide 87 integer :: idum,ii,info,inversion_flag 88 integer :: ir,iwan,jj,ll,lm,lmc,mm,mr 89 real(dp):: onem,test 90 ! arrays 91 integer:: ipiv(lmax2) 92 real(dp)::r(3,lmax2),rp(3,lmax2) 93 real(dp)::rs2,rs3,rs6,rs12,umat(3,3) 94 complex(dp)::crot(lmax2,lmax2),ctor(lmax2,lmax2),orb_lm(lmax2,-5:3,7) 95 complex(dp):: ylmcp(lmax2) 96 complex(dp):: ylmc_rr(lmax2,lmax2),ylmc_rr_save(lmax2,lmax2) 97 complex(dp):: ylmc_rrinv(lmax2,lmax2),ylmc_rp(lmax2,lmax2) 98 complex(dp),parameter :: c0=(0._dp,0._dp),c1=(1._dp,0._dp),ci=(0._dp,1._dp) 99 character(len=500) :: message ! to be uncommented, if needed 100 101 ! ************************************************************************* 102 103 104 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 105 !DEBUG 106 !write(std_out,*)'lmax ',lmax,'lmax2 ',lmax2 107 !write(std_out,*)'mband ',mband,'nwan ',nwan 108 ! 109 !do iwan=1,nwan 110 !write(std_out,*)'iwan,proj_l, proj_m',proj_l(iwan),proj_m(iwan) 111 !write(std_out,*)'iwan,proj_x, proj_z',iwan,proj_x(:,iwan),proj_z(:,iwan) 112 !end do 113 !!END DEBUG 114 115 !constants for linear combinations of ylm's 116 rs2=1._dp/sqrt(2._dp) 117 rs3=1._dp/sqrt(3._dp) 118 rs6=1._dp/sqrt(6._dp) 119 rs12=1._dp/sqrt(12._dp) 120 121 !complex lm coefficients for real spherical harmonics in conventional order 122 !s, py,pz,px, dxy,dyz,dz2,dxz,dx2-y2, fy(3x2-y2),fxyz,fyz2,fz3,fxz2, 123 !fz(x2-y2),fx(x2-3y2) 124 ctor(:,:)=c0 125 do ll=0,lmax 126 mm=0 127 lm= ll**2+ll+mm+1 128 ctor(lm,lm)=c1 129 if(ll>0) then 130 onem=one 131 do mm=1,ll 132 onem=-onem !(-1^mm) 133 lm= ll**2+ll+mm+1 134 lmc=ll**2+ll-mm+1 135 ctor(lm ,lm )=rs2*c1 136 ctor(lmc,lm )=onem*rs2*c1 137 ctor(lm ,lmc)=rs2*ci 138 ctor(lmc,lmc)=-onem*rs2*ci 139 end do 140 end if 141 end do 142 143 lm=0 144 do ll=0,lmax 145 do mm=-ll,ll 146 lm=lm+1 147 ctor(:,lm)=ctor(:,lm)*conjg(ci)**ll 148 end do !mm 149 end do !ll 150 151 152 !coefficients for basic wannier orbitals in Table 3.1 order 153 orb_lm(:,:,:)=c0 154 ii=0 155 do ll=0,lmax 156 do mr=1,2*ll+1 157 ii=ii+1 158 orb_lm(:,ll,mr)=ctor(:,orb_idx(ii)) 159 end do 160 end do 161 162 163 164 !coefficients for linear combinations in table 3.2 order 165 if(lmax>=1) then 166 ! s px 167 orb_lm(:,-1,1)=rs2*ctor(:,1)+rs2*ctor(:,4) 168 orb_lm(:,-1,2)=rs2*ctor(:,1)-rs2*ctor(:,4) 169 ! s px py 170 orb_lm(:,-2,1)=rs3*ctor(:,1)-rs6*ctor(:,4)+rs2*ctor(:,2) 171 orb_lm(:,-2,2)=rs3*ctor(:,1)-rs6*ctor(:,4)-rs2*ctor(:,2) 172 orb_lm(:,-2,3)=rs3*ctor(:,1)+2._dp*rs6*ctor(:,4) 173 ! s px py pz 174 orb_lm(:,-3,1)=half*(ctor(:,1)+ctor(:,4)+ctor(:,2)+ctor(:,3)) 175 orb_lm(:,-3,2)=half*(ctor(:,1)+ctor(:,4)-ctor(:,2)-ctor(:,3)) 176 orb_lm(:,-3,3)=half*(ctor(:,1)-ctor(:,4)+ctor(:,2)-ctor(:,3)) 177 orb_lm(:,-3,4)=half*(ctor(:,1)-ctor(:,4)-ctor(:,2)+ctor(:,3)) 178 end if 179 if(lmax>=2) then 180 ! s px py 181 orb_lm(:,-4,1)=rs3*ctor(:,1)-rs6*ctor(:,4)+rs2*ctor(:,2) 182 orb_lm(:,-4,2)=rs3*ctor(:,1)-rs6*ctor(:,4)-rs2*ctor(:,2) 183 orb_lm(:,-4,3)=rs3*ctor(:,1)+2._dp*rs6*ctor(:,4) 184 ! pz dz2 185 orb_lm(:,-4,4)= rs2*ctor(:,3)+rs2*ctor(:,7) 186 orb_lm(:,-4,5)=-rs2*ctor(:,3)+rs2*ctor(:,7) 187 ! s px dz2 dx2-y2 188 orb_lm(:,-5,1)=rs6*ctor(:,1)-rs2*ctor(:,4)-rs12*ctor(:,7)+half*ctor(:,9) 189 orb_lm(:,-5,2)=rs6*ctor(:,1)+rs2*ctor(:,4)-rs12*ctor(:,7)+half*ctor(:,9) 190 ! s py dz2 dx2-y2 191 orb_lm(:,-5,3)=rs6*ctor(:,1)-rs2*ctor(:,2)-rs12*ctor(:,7)-half*ctor(:,9) 192 orb_lm(:,-5,4)=rs6*ctor(:,1)+rs2*ctor(:,2)-rs12*ctor(:,7)-half*ctor(:,9) 193 ! s pz dz2 194 orb_lm(:,-5,5)=rs6*ctor(:,1)-rs2*ctor(:,3)+rs3*ctor(:,7) 195 orb_lm(:,-5,6)=rs6*ctor(:,1)+rs2*ctor(:,3)+rs3*ctor(:,7) 196 end if 197 198 !stuff complex wannier orbital coefficient array 199 do iwan=1,nwan 200 ylmc_fac(:,iwan)=orb_lm(:,proj_l(iwan),proj_m(iwan)) 201 end do 202 203 204 !setup to rotate ylmc_fac to new axes if called for 205 !skip if only s projectors are used 206 if ( lmax>0 ) then 207 ! generate a set of nr=lmax2 random vectors 208 ! idum=123456 209 do ir=1,lmax2 210 do ii=1,3 211 r(ii,ir) = uniformrandom(idum)-0.5d0 212 end do !ii 213 call ylm_cmplx(lmax,ylmcp,r(1,ir),r(2,ir),r(3,ir)) 214 ylmc_rr(ir,:)=conjg(ylmcp(:)) 215 ylmc_rr_save(ir,:)=conjg(ylmcp(:)) 216 end do !ir 217 218 ylmc_rrinv(:,:)=c0 219 do ii=1,lmax2 220 ylmc_rrinv(ii,ii)=c1 221 end do !ii 222 ! calculate inverse of ylmc(ir,lm) matrix 223 call ZGESV(lmax2,lmax2,ylmc_rr,lmax2,ipiv,ylmc_rrinv,lmax2,info) 224 225 ! check that r points are independent (ie., that matrix inversion wasn't 226 ! too close to singular) 227 ylmc_rr=matmul(ylmc_rrinv,ylmc_rr_save) 228 test=zero 229 do ii=1,lmax2 230 ylmc_rr(ii,ii)=ylmc_rr(ii,ii)-c1 231 do jj=1,lmax2 232 test=max(abs(ylmc_rr(ii,jj)),test) 233 end do !ii 234 end do !jj 235 if(test>tol8) then 236 write(message, '(5a)' )& 237 & ' matrix inversion error for wannier rotations',ch10,& 238 & ' random vectors r(j,1:nr) are not all independent !! ',ch10,& 239 & ' Action : re-seed uniformrandom or maybe just try again' 240 MSG_ERROR(message) 241 end if !test>tol8 242 243 ! end of the preliminaries, now to the rotations of the wannier orbitals 244 do iwan=1,nwan 245 ! don't bother for s orbitals 246 if(proj_l(iwan)==0) cycle 247 ! check for default axes and cycle if found 248 if(proj_z(1,iwan)==zero .and. proj_z(2,iwan)==zero .and.& 249 & proj_z(3,iwan)== one .and. proj_x(1,iwan)==one .and.& 250 & proj_x(2,iwan)==zero .and. proj_x(3,iwan)==zero) cycle 251 252 ! get the u matrix that rotates the reference frame 253 call rotmat(proj_x(:,iwan),proj_z(:,iwan),inversion_flag,umat) 254 255 ! find rotated r-vectors. Optional inversion 256 ! operation is an extension of the wannier90 axis-setting options 257 ! which only allow for proper axis rotations 258 if(inversion_flag==1) then 259 rp(:,:)= -matmul ( umat(:,:), r(:,:) ) 260 else 261 rp(:,:) = matmul ( umat(:,:) , r(:,:) ) 262 end if !inversion_flag 263 264 do ir=1,lmax2 265 ! get the ylm representation of the rotated vectors 266 call ylm_cmplx(lmax,ylmcp,rp(1,ir),rp(2,ir),rp(3,ir)) 267 ylmc_rp(ir,:)=conjg(ylmcp(:)) 268 end do !ir 269 ! the matrix product sum(ir) ylmc_rrinv(lm,ir)*ylmc_rp(ir,lm') gives the 270 ! the complex lmXlm matrix representation of the coordinate rotation 271 crot(:,:)=matmul(ylmc_rrinv(:,:),ylmc_rp(:,:)) 272 273 ! now rotate the current wannier orbital 274 ylmcp(:)=matmul(crot(:,:),ylmc_fac(:,iwan)) 275 ylmc_fac(:,iwan)=ylmcp(:) 276 277 ! write(std_out,*)'ylmc_fac',ylmc_fac(:,iwan) 278 end do !iwan 279 end if !lmax>0 280 281 end subroutine mlwfovlp_ylmfac