1!
2! Copyright (C) 2004 PWSCF group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!--------------------------------------------------------------------------
9subroutine compute_phi_tm(lam,ik,chir,phi_out,iflag,xc,e,els_in)
10  !--------------------------------------------------------------------------
11  !
12  !     This routine computes the phi functions by pseudizing the
13  !     all_electron chi functions. In input it receives, the point
14  !     ik where the cut is done, the angular momentum lam,
15  !     and the all electron wavefunction
16  !
17  !
18  !
19  use io_global, only : stdout
20  use kinds, only : DP
21  use constants, only: pi
22  use radial_grids, only: ndmx
23  use ld1inc, only: grid, vpot
24  implicit none
25
26  integer ::    &
27       lam,  &   ! input: the angular momentum
28       ik,   &   ! input: the point corresponding to rc
29       iflag     ! input: if 1 print the message
30
31  real(DP) :: &
32       e         ! input: the energy of the level
33
34  real(DP) ::       &
35       xc(8),       &    ! output: parameters of the fit
36       chir(ndmx),   &    ! input: the all-electron function
37       phi_out(ndmx)      ! output: the phi function
38
39  character(len=2) :: els_in  ! input: the label of the state
40  !
41  real(DP) :: &
42       faenor    ! the norm of the function
43
44  integer ::    &
45       n, nst
46
47  real(DP) :: &
48       gi(ndmx), &
49       cn(6), c2
50
51  real(DP), external :: deriv_7pts, deriv2_7pts, int_0_inf_dr, pr
52  !
53  !
54  !   compute the norm of the all-electron function
55  !
56  nst=(lam+1)*2
57  do n=1,ik+1
58     gi(n)=chir(n)**2
59  enddo
60  faenor=int_0_inf_dr(gi,grid,ik,nst)
61  !
62  !
63  ! TM: the pseudo-wavefunction is written as polynomial times exponential
64  !
65  call find_coefficients &
66          (ik, chir, e, grid%r, grid%dx, faenor, vpot, cn, c2, lam, grid%mesh)
67  do n=1,ik
68     phi_out(n)=sign(grid%r(n)**(lam+1)*exp(pr(cn,c2,grid%r(n))),chir(ik))
69  end do
70  xc(1:6) = cn(:)
71  xc(7) = c2
72  !
73  !      for r > r(ik) the pseudo and all-electron psi(r) coincide
74  !
75  do n=ik+1,grid%mesh
76     phi_out(n)= chir(n)
77  enddo
78  !
79  !      check for the norm of the pseudo wavefunction
80  !
81!  do n=1,ik
82!     gi(n)=phi_out(n)**2
83!  enddo
84!  psnor=int_0_inf_dr(gi,grid,ik,nst)
85! write(stdout,'(5x," AE norm = ",f6.3,"  PS norm = ",f6.3)') faenor, psnor
86
87  if (iflag == 1) then
88     write(stdout,120) els_in, grid%r(ik)
89120  format (/ /5x, ' Wfc  ',a3,'  rcut=',f6.3, &
90       '  Using Troullier-Martins method ')
91  end if
92
93  return
94end subroutine compute_phi_tm
95