TABLE OF CONTENTS
ABINIT/mlwfovlp_ylmfar [ Functions ]
NAME
mlwfovlp_ylmfar
FUNCTION
Routine that produces a fator 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] This function is similar to mlwfovlp_ylmfac, but the factors it uses real spherical harmonics instead of complex spherical harmonics. Remember that real spherical harmonics are linear combinations of complex spherical harmonics
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (T. Rangel) 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
NOTES
PARENTS
mlwfovlp_projpaw
CHILDREN
initylmr,matrginv,rotmat
SOURCE
54 #if defined HAVE_CONFIG_H 55 #include "config.h" 56 #endif 57 58 #include "abi_common.h" 59 60 subroutine mlwfovlp_ylmfar(ylmr_fac,lmax,lmax2,mband,nwan,proj_l,proj_m,proj_x,proj_z) 61 62 use defs_basis 63 use m_errors 64 use m_profiling_abi 65 66 use m_geometry, only : rotmat 67 use m_numeric_tools, only : uniformrandom 68 use m_abilasi, only : matrginv 69 use m_paw_sphharm, only : initylmr 70 71 !This section has been created automatically by the script Abilint (TD). 72 !Do not modify the following lines by hand. 73 #undef ABI_FUNC 74 #define ABI_FUNC 'mlwfovlp_ylmfar' 75 !End of the abilint section 76 77 implicit none 78 79 !Arguments ------------------------------------ 80 integer, intent(in):: lmax,lmax2,nwan,mband 81 ! arrays 82 integer,intent(in) :: proj_l(mband),proj_m(mband) 83 real(dp),intent(in) :: proj_x(3,mband),proj_z(3,mband) 84 real(dp),intent(out)::ylmr_fac(lmax2,nwan) 85 ! 86 !Local variables------------------------------- 87 ! 88 integer :: idum,ii,inversion_flag 89 integer :: ir,iwan,jj,ll,lm,mm,mr 90 real(dp):: onem,test 91 ! arrays 92 real(dp),allocatable::dummy(:,:),nrm(:) 93 real(dp)::r(3,lmax2),rp(3,lmax2) 94 real(dp)::rs2,rs3,rs6,rs12,umat(3,3) 95 real(dp)::rot(lmax2,lmax2),tor(lmax2,lmax2),orb_lm(lmax2,-5:3,7) 96 real(dp):: ylmrp(lmax2) 97 real(dp):: ylmr_rr(lmax2,lmax2),ylmr_rr_save(lmax2,lmax2) 98 real(dp):: ylmr_rrinv(lmax2,lmax2),ylmr_rp(lmax2,lmax2) 99 character(len=500) :: message ! to be uncommented, if needed 100 !no_abirules 101 !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 102 103 ! ************************************************************************* 104 105 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 106 !DEBUG 107 !write(std_out,*)'lmax ',lmax,'lmax2 ',lmax2 108 !write(std_out,*)'mband ',mband,'nwan ',nwan 109 ! 110 !do iwan=1,nwan 111 !write(std_out,*)'iwan,proj_l, proj_m',proj_l(iwan),proj_m(iwan) 112 !write(std_out,*)'iwan,proj_x, proj_z',iwan,proj_x(:,iwan),proj_z(:,iwan) 113 !end do 114 !!END DEBUG 115 116 !constants for linear combinations of ylm's 117 rs2=1._dp/sqrt(2._dp) 118 rs3=1._dp/sqrt(3._dp) 119 rs6=1._dp/sqrt(6._dp) 120 rs12=1._dp/sqrt(12._dp) 121 122 ! 123 !mapping lm coefficients for real spherical harmonics 124 !table 3.1 of Wannier90 user guide with real spherical harmonics in routine initylmr 125 !s, py,pz,px, dxy,dyz,dz2,dxz,dx2-y2, fy(3x2-y2),fxyz,fyz2,fz3,fxz2, 126 !fz(x2-y2),fx(x2-3y2) 127 !note: check ordering of f orbitals, it might be wrong 128 129 tor(:,:)=0.d0 130 lm=0 131 do ll=0,lmax 132 do mm=-ll,ll 133 onem=(-1.d0)**mm 134 lm=lm+1 135 if(ll == 0) then 136 tor(lm,lm)=1.d0 137 else 138 tor(lm,lm)=onem*1.d0 139 end if 140 end do !mm 141 end do !ll 142 !do lm=1,16 143 !write(std_out,*)'tor lm=',lm,tor(:,lm) 144 !end do 145 146 !coefficients for basic wannier orbitals in Table 3.1 order 147 orb_lm(:,:,:)=0.d0 148 ii=0 149 do ll=0,lmax 150 do mr=1,2*ll+1 151 ii=ii+1 152 orb_lm(:,ll,mr)= tor(:,ii) 153 ! write(std_out,*)'ii',ii,'orb_lm',orb_lm(:,ll,mr) 154 end do 155 end do 156 157 158 159 !coefficients for linear combinations in table 3.2 order 160 if(lmax>=1) then 161 ! s px 162 orb_lm(:,-1,1)=rs2*tor(:,1)+rs2*tor(:,4) 163 orb_lm(:,-1,2)=rs2*tor(:,1)-rs2*tor(:,4) 164 ! s px py 165 orb_lm(:,-2,1)=rs3*tor(:,1)-rs6*tor(:,4)+rs2*tor(:,2) 166 orb_lm(:,-2,2)=rs3*tor(:,1)-rs6*tor(:,4)-rs2*tor(:,2) 167 orb_lm(:,-2,3)=rs3*tor(:,1)+2._dp*rs6*tor(:,4) 168 ! s px py pz 169 orb_lm(:,-3,1)=half*(tor(:,1)+tor(:,4)+tor(:,2)+tor(:,3)) 170 orb_lm(:,-3,2)=half*(tor(:,1)+tor(:,4)-tor(:,2)-tor(:,3)) 171 orb_lm(:,-3,3)=half*(tor(:,1)-tor(:,4)+tor(:,2)-tor(:,3)) 172 orb_lm(:,-3,4)=half*(tor(:,1)-tor(:,4)-tor(:,2)+tor(:,3)) 173 end if 174 if(lmax>=2) then 175 ! s px py 176 orb_lm(:,-4,1)=rs3*tor(:,1)-rs6*tor(:,4)+rs2*tor(:,2) 177 orb_lm(:,-4,2)=rs3*tor(:,1)-rs6*tor(:,4)-rs2*tor(:,2) 178 orb_lm(:,-4,3)=rs3*tor(:,1)+2._dp*rs6*tor(:,4) 179 ! pz dz2 180 orb_lm(:,-4,4)= rs2*tor(:,3)+rs2*tor(:,7) 181 orb_lm(:,-4,5)=-rs2*tor(:,3)+rs2*tor(:,7) 182 ! s px dz2 dx2-y2 183 orb_lm(:,-5,1)=rs6*tor(:,1)-rs2*tor(:,4)-rs12*tor(:,7)+half*tor(:,9) 184 orb_lm(:,-5,2)=rs6*tor(:,1)+rs2*tor(:,4)-rs12*tor(:,7)+half*tor(:,9) 185 ! s py dz2 dx2-y2 186 orb_lm(:,-5,3)=rs6*tor(:,1)-rs2*tor(:,2)-rs12*tor(:,7)-half*tor(:,9) 187 orb_lm(:,-5,4)=rs6*tor(:,1)+rs2*tor(:,2)-rs12*tor(:,7)-half*tor(:,9) 188 ! s pz dz2 189 orb_lm(:,-5,5)=rs6*tor(:,1)-rs2*tor(:,3)+rs3*tor(:,7) 190 orb_lm(:,-5,6)=rs6*tor(:,1)+rs2*tor(:,3)+rs3*tor(:,7) 191 end if 192 193 !real wannier orbital coefficient array 194 do iwan=1,nwan 195 ylmr_fac(:,iwan)=orb_lm(:,proj_l(iwan),proj_m(iwan)) 196 end do 197 198 199 !setup to rotate ylmr_fac to new axes if called for 200 !skip if only s projetors are used 201 if ( lmax>0 ) then 202 ! generate a set of nr=lmax2 random vetors 203 idum=123456 204 do ir=1,lmax2 205 do ii=1,3 206 r(ii,ir) = uniformrandom(idum)-0.5d0 207 end do !ii 208 end do !ir 209 ABI_ALLOCATE(nrm,(lmax2)) 210 nrm(:)=sqrt(r(1,:)**2+r(2,:)**2+r(3,:)**2)**0.5 211 call initylmr(lmax+1,1,lmax2,nrm,1,r(:,:),ylmr_rr_save(:,:),dummy) 212 ylmr_rr(:,:)=ylmr_rr_save(:,:) 213 do ir=1,lmax2 214 ylmr_rr_save(ir,:)=ylmr_rr(:,ir) 215 end do 216 ABI_DEALLOCATE(nrm) 217 218 ylmr_rrinv(:,:)=0.d0 219 do ii=1,lmax2 220 ylmr_rrinv(ii,ii)=1.d0 221 end do !ii 222 ! calculate inverse of ylmr(ir,lm) matrix 223 ylmr_rrinv(:,:)=ylmr_rr_save(:,:) 224 call matrginv(ylmr_rrinv,lmax2,lmax2) 225 226 ! check that r points are independent (ie., that matrix inversion wasn't 227 ! too close to singular) 228 ylmr_rr=matmul(ylmr_rrinv,ylmr_rr_save) 229 test=0.d0 230 do ii=1,lmax2 231 ylmr_rr(ii,ii)=ylmr_rr(ii,ii)-1.d0 232 do jj=1,lmax2 233 test=max(abs(ylmr_rr(ii,jj)),test) 234 end do !ii 235 end do !jj 236 if(test>tol8) then 237 write(message, '(5a)' )& 238 & ' matrix inversion error for wannier rotations',ch10,& 239 & ' random vetors r(j,1:nr) are not all independent !! ',ch10,& 240 & ' Action : re-seed uniformrandom or maybe just try again' 241 MSG_ERROR(message) 242 end if !test>tol8 243 244 ! end of the preliminaries, now to the rotations of the wannier orbitals 245 do iwan=1,nwan 246 ! don't bother for s orbitals 247 if(proj_l(iwan)==0) cycle 248 ! check for default axes and cycle if found 249 if(proj_z(1,iwan)==0.d0 .and. proj_z(2,iwan)==0.d0 .and.& 250 & proj_z(3,iwan)== 1.d0 .and. proj_x(1,iwan)==1.d0 .and.& 251 & proj_x(2,iwan)==0.d0 .and. proj_x(3,iwan)==0.d0) cycle 252 253 ! get the u matrix that rotates the reference frame 254 call rotmat(proj_x(:,iwan),proj_z(:,iwan),inversion_flag,umat) 255 ! 256 ! find rotated r-vetors. Optional inversion 257 ! operation is an extension of the wannier90 axis-setting options 258 ! which only allow for proper axis rotations 259 if(inversion_flag==1) then 260 rp(:,:)= -matmul ( umat(:,:), r(:,:) ) 261 else 262 rp(:,:) = matmul ( umat(:,:) , r(:,:) ) 263 end if !inversion_flag 264 265 ! get the ylm representation of the rotated vetors 266 ABI_ALLOCATE(nrm,(lmax2)) 267 nrm(:)=sqrt(rp(1,:)**2+rp(2,:)**2+rp(3,:)**2)**0.5 268 call initylmr(lmax+1,1,lmax2,nrm,1,rp(:,:),ylmr_rp(:,:),dummy) 269 ylmr_rr(:,:)=ylmr_rp(:,:) 270 do ir=1,lmax2 271 ylmr_rp(ir,:)=ylmr_rr(:,ir) 272 end do 273 ABI_DEALLOCATE(nrm) 274 ! the matrix product sum(ir) ylmr_rrinv(lm,ir)*ylmr_rp(ir,lm') gives the 275 ! the lmXlm matrix representation of the coordinate rotation 276 277 rot(:,:)=matmul(ylmr_rrinv(:,:),ylmr_rp(:,:)) 278 ! 279 ! now rotate the current wannier orbital 280 ylmrp(:)=matmul(rot(:,:),ylmr_fac(:,iwan)) 281 ylmr_fac(:,iwan)=ylmrp(:) 282 end do !iwan 283 end if !lmax>0 284 285 !DEBUG 286 !write (std_out,*) ' mlwfovlp_ylmfar : exit' 287 !stop 288 !ENDDEBUG 289 290 end subroutine mlwfovlp_ylmfar