TABLE OF CONTENTS

ABINIT/mlwfovlp_radial [ Functions ]

[ Top ] [ Functions ]

NAME

 mlwfovlp_radial

FUNCTION

 Calculates the radial part of the initial functions given to Wannier90 
 as an starting point for the minimization.
 The trial functions are a set of solutions to the radial part of the hydrogenic
 Schrodinger equation as it is explained in Table 3.3 of the Wannier90 user guide.

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (trangel,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

  alpha= Z/a = zona
  lmax= maximum value of l
  rvalue= integer defining the choice for radial functions R(r).
   It can take values from 1-3. 
   It is associted to the radial part of the hydrogenic Schrodinger equation for l=0,
   See the manual of Wannier90 for more information. (www.wannier.org)
  xx= scalar number used to calculate the spherical bessel function. J_il(xx)

OUTPUT

  mlwfovlp_radial= radial part for initial projections used to construct MLWF

SIDE EFFECTS

  None

NOTES

  Calculates the radial part of the initial functions given as an initial
  guess by the user to construct the MLWF.

PARENTS

      mlwfovlp_proj

CHILDREN

      besjm,simpson_int

SOURCE

 45 #if defined HAVE_CONFIG_H
 46 #include "config.h"
 47 #endif
 48 
 49 #include "abi_common.h"
 50 
 51 subroutine mlwfovlp_radial(alpha,lmax,lmax2,radial,rvalue,xx)
 52 
 53  use defs_basis
 54  use m_errors
 55  use m_profiling_abi
 56 
 57  use m_numeric_tools,   only : simpson_int
 58  use m_special_funcs,   only : besjm
 59 
 60 !This section has been created automatically by the script Abilint (TD).
 61 !Do not modify the following lines by hand.
 62 #undef ABI_FUNC
 63 #define ABI_FUNC 'mlwfovlp_radial'
 64 !End of the abilint section
 65 
 66  implicit none
 67 
 68 !Arguments ------------------------------------
 69 !scalars
 70  integer,intent(in) :: lmax,lmax2,rvalue
 71  real(dp),intent(in) :: alpha,xx
 72 !arrays
 73  real(dp),intent(out) :: radial(lmax2)
 74 
 75 !Local variables
 76 !scalars
 77  integer :: ir,ll,lm,mesh,mm
 78  real(dp),parameter :: dx=0.015d0,rmax=10.d0,xmin=0.d0
 79  real(dp) :: aa,ftmp,gauss,rtmp,x
 80  character(len=500) :: message
 81 !arrays
 82  real(dp),save :: dblefact(4)=(/1_dp,3_dp,15_dp,105_dp/)
 83  real(dp),allocatable :: aux(:),bes(:),cosr(:),func_r(:),r(:),rad_int(:)
 84  real(dp),allocatable :: sinr(:)
 85 
 86 ! *************************************************************************
 87  
 88 !Radial functions in the form of hydrogenic orbitals as defined in the 
 89 !wannier90 manual.
 90  if(( rvalue > 0 ).and.(rvalue < 4)) then
 91 
 92 !  mesh
 93    mesh= nint((rmax - xmin ) / dx + 1)
 94    ABI_ALLOCATE( bes,(mesh))
 95    ABI_ALLOCATE(func_r,(mesh))
 96    ABI_ALLOCATE(r,(mesh))
 97    ABI_ALLOCATE(rad_int,(mesh))
 98    ABI_ALLOCATE( aux,(mesh))
 99    ABI_ALLOCATE(cosr,(mesh))
100    ABI_ALLOCATE(sinr,(mesh))
101    do ir=1, mesh
102      x=xmin+DBLE(ir-1)*dx
103      r(ir)=x
104    end do   !ir
105 
106 !  radial functions shown in table 3.3 of wannier90 manual
107    if (rvalue==1) func_r(:) = 2.d0 * alpha**(3.d0/2.d0) * exp(-alpha*r(:))
108    if (rvalue==2) func_r(:) = 1.d0/(2.d0*sqrt(2.d0))*alpha**(3.d0/2.d0) *&
109 &   (2.d0 - alpha*r(:))*exp(-alpha*r(:)/2.d0)
110    if (rvalue==3) func_r(:) = sqrt(4.d0/27.d0)*alpha**(3.d0/2.d0)&
111 &   * (1.d0 - 2.d0*alpha*r(:)/3.d0 + 2.d0*alpha**2*r(:)**2/27.d0)&
112 &   * exp(-alpha * r(:)/3.d0)
113 
114 !  compute spherical bessel functions
115    cosr(:)=cos(xx*r(:))
116    sinr(:)=sin(xx*r(:))
117    lm=0
118    do ll=0,lmax
119      call besjm(xx,bes,cosr,ll,mesh,sinr,r)
120      aux(:)=bes(:)*func_r(:)*r(:)
121 !    do ir=1,mesh
122 !    write(310,*) r(ir),bes(ir)
123 !    end do
124      call simpson_int(mesh,dx,aux,rad_int)
125      rtmp=rad_int(mesh)/mesh
126      do mm=-ll,ll
127        lm=lm+1
128        radial(lm)=rtmp
129      end do !mm
130    end do !ll
131    ABI_DEALLOCATE(bes)
132    ABI_DEALLOCATE(func_r)
133    ABI_DEALLOCATE(r)
134    ABI_DEALLOCATE(aux)
135    ABI_DEALLOCATE(rad_int)
136    ABI_DEALLOCATE(cosr)
137    ABI_DEALLOCATE(sinr)
138 
139 !  Radial part in the form of Gaussian functions of a given width 
140 !  Taken by code of made by drh.
141  elseif ( rvalue == 4) then
142    aa=1._dp/alpha
143    gauss=exp(-0.25_dp*(aa*xx)**2)
144    lm=0
145    do ll=0,lmax
146      ftmp=(0.5_dp*pi)**(0.25_dp)*aa*sqrt(aa/dblefact(ll+1))*(aa*xx)**ll*gauss
147      do mm=-ll,ll
148        lm=lm+1
149        radial(lm)=ftmp
150      end do
151    end do
152  else ! rvalue < 0 of rvalue > 4
153    write(message,'(a,i6,5a)')&
154 &   '  Radial function r=',rvalue,ch10,&
155 &   '  is not defined',ch10,&
156 &   '  Modify .win file',ch10
157    MSG_BUG(message)
158  end if !rvalue
159 
160 end subroutine mlwfovlp_radial