src/wavefmt.f90

   1
   2 ! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
   3 ! This file is distributed under the terms of the GNU General Public License.
   4 ! See the file COPYING for license details.
   5
   6 !BOP
   7 ! !ROUTINE: wavefmt
   8 ! !INTERFACE:
   9 subroutine wavefmt(lrstp,lmax,is,ia,ngp,apwalm,evecfv,ld,wfmt)
  10 ! !USES:
  11 use modmain
  12 ! !INPUT/OUTPUT PARAMETERS:
  13 !   lrstp  : radial step length (in,integer)
  14 !   lmax   : maximum angular momentum required (in,integer)
  15 !   is     : species number (in,integer)
  16 !   ia     : atom number (in,integer)
  17 !   ngp    : number of G+p-vectors (in,integer)
  18 !   apwalm : APW matching coefficients
  19 !            (in,complex(ngkmax,apwordmax,lmmaxapw,natmtot))
  20 !   evecfv : first-variational eigenvector (in,complex(nmatmax))
  21 !   ld     : leading dimension (in,integer)
  22 !   wfmt   : muffin-tin wavefunction (out,complex(ld,*))
  23 ! !DESCRIPTION:
  24 !   Calculates the first-variational wavefunction in the muffin-tin in terms of
  25 !   a spherical harmonic expansion. For atom $\alpha$ and a particular
  26 !   $p$-point, the $r$-dependent $(l,m)$-coefficients of the wavefunction for
  27 !   the $i$th state are given by
  28 !   $$ \Psi^{i{\bf p}\alpha}_{lm}(r)=\sum_{\bf G}\Phi^{i{\bf p}}_{\bf G}
  29 !    \sum_{j=1}^{M^{\alpha}_l}A^{\alpha}_{jlm}({\bf G+p})u^{\alpha}_{jl}(r)
  30 !    +\sum_{j=1}^{N^{\alpha}}\Phi^{i{\bf p}}_{(\alpha,j,m)}v^{\alpha}_j(r)
  31 !    \delta_{l,l_j}, $$
  32 !   where $\Phi^{i{\bf p}}$ is the $i$th eigenvector returned from routine
  33 !   {\tt seceqn}; $A^{\alpha}_{jlm}({\bf G+p})$ is the matching coefficient;
  34 !   $M^{\alpha}_l$ is the order of the APW; $u^{\alpha}_{jl}$ is the APW radial
  35 !   function; $N^{\alpha}$ is the number of local-orbitals; $v^{\alpha}_j$ is
  36 !   the $j$th local-orbital radial function; and $(\alpha,j,m)$ is a compound
  37 !   index for the location of the local-orbital in the eigenvector. See routines
  38 !   {\tt genapwfr}, {\tt genlofr}, {\tt match} and {\tt seceqn}.
  39 !
  40 ! !REVISION HISTORY:
  41 !   Created April 2003 (JKD)
  42 !   Fixed description, October 2004 (C. Brouder)
  43 !   Removed argument ist, November 2006 (JKD)
  44 !EOP
  45 !BOC
  46 implicit none
  47 ! arguments
  48 integer, intent(in) :: lrstp
  49 integer, intent(in) :: lmax
  50 integer, intent(in) :: is
  51 integer, intent(in) :: ia
  52 integer, intent(in) :: ngp
  53 complex(8), intent(in) :: apwalm(ngkmax,apwordmax,lmmaxapw,natmtot)
  54 complex(8), intent(in) :: evecfv(nmatmax)
  55 integer, intent(in) :: ld
  56 complex(8), intent(out) :: wfmt(ld,*)
  57 ! local variables
  58 integer ias,l,m,lm,i
  59 integer ir,nr,io,ilo
  60 real(8) a,b
  61 complex(8) zt1
  62 ! external functions
  63 complex(8) zdotu
  64 external zdotu
  65 if (lmax.gt.lmaxapw) then
  66   write(*,*)
  67   write(*,'("Error(wavefmt): lmax > lmaxapw : ",I8)') lmax
  68   write(*,*)
  69   stop
  70 end if
  71 ias=idxas(ia,is)
  72 ! zero the wavefunction
  73 nr=0
  74 do ir=1,nrmt(is),lrstp
  75   nr=nr+1
  76   wfmt(:,nr)=0.d0
  77 end do
  78 ! APW functions
  79 do l=0,lmax
  80   do m=-l,l
  81     lm=idxlm(l,m)
  82     do io=1,apword(l,is)
  83       zt1=zdotu(ngp,evecfv,1,apwalm(1,io,lm,ias),1)
  84       a=dble(zt1)
  85       b=aimag(zt1)
  86       call wavefmt_add(nr,ld,wfmt(lm,1),a,b,lrstp,apwfr(1,1,io,l,ias))
  87     end do
  88   end do
  89 end do
  90 ! local-orbital functions
  91 do ilo=1,nlorb(is)
  92   l=lorbl(ilo,is)
  93   if (l.le.lmax) then
  94     do m=-l,l
  95       lm=idxlm(l,m)
  96       i=ngp+idxlo(lm,ilo,ias)
  97       a=dble(evecfv(i))
  98       b=aimag(evecfv(i))
  99       call wavefmt_add(nr,ld,wfmt(lm,1),a,b,lrstp,lofr(1,1,ilo,ias))
 100     end do
 101   end if
 102 end do
 103 return
 104 end subroutine
 105 !EOC
 106
 107 !BOP
 108 ! !ROUTINE: wavefmt_add
 109 ! !INTERFACE:
 110 subroutine wavefmt_add(nr,ld,wfmt,a,b,lrstp,fr)
 111 ! !INPUT/OUTPUT PARAMETERS:
 112 !   nr     : number of radial mesh points (in,integer)
 113 !   ld     : leading dimension (in,integer)
 114 !   wfmt   : complex muffin-tin wavefunction passed in as a real array
 115 !            (inout,real(2*ld,*))
 116 !   a      : real part of complex constant (in,real)
 117 !   b      : imaginary part of complex constant (in,real)
 118 !   lrstp  : radial step length (in,integer)
 119 !   fr     : real radial function (in,real(lrstp,*))
 120 ! !DESCRIPTION:
 121 !   Adds a real function times a complex constant to a complex muffin-tin
 122 !   wavefunction as efficiently as possible. See routine {\tt wavefmt}.
 123 !
 124 ! !REVISION HISTORY:
 125 !   Created December 2006 (JKD)
 126 !EOP
 127 !BOC
 128 implicit none
 129 ! arguments
 130 integer, intent(in) :: nr
 131 integer, intent(in) :: ld
 132 real(8), intent(inout) :: wfmt(2*ld,*)
 133 real(8), intent(in) :: a
 134 real(8), intent(in) :: b
 135 integer, intent(in) :: lrstp
 136 real(8), intent(in) :: fr(lrstp,*)
 137 ! local variables
 138 integer ir
 139 ! values smaller than eps are taken to be zero
 140 real(8), parameter :: eps=1.d-14
 141 if (abs(b).lt.eps) then
 142 ! zero constant
 143   if (abs(a).lt.eps) return
 144 ! pure real constant
 145   do ir=1,nr
 146     wfmt(1,ir)=wfmt(1,ir)+a*fr(1,ir)
 147   end do
 148 else if (abs(a).lt.eps) then
 149 ! pure imaginary constant
 150   do ir=1,nr
 151     wfmt(2,ir)=wfmt(2,ir)+b*fr(1,ir)
 152   end do
 153 else
 154 ! general complex constant
 155   do ir=1,nr
 156     wfmt(1,ir)=wfmt(1,ir)+a*fr(1,ir)
 157     wfmt(2,ir)=wfmt(2,ir)+b*fr(1,ir)
 158   end do
 159 end if
 160 return
 161 end subroutine
 162 !EOC
 163