TABLE OF CONTENTS


ABINIT/mlwfovlp_ylmfar [ Functions ]

[ Top ] [ 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