File: genapwfr.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (131 lines) | stat: -rw-r--r-- 4,103 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

!BOP
! !ROUTINE: genapwfr
! !INTERFACE:
subroutine genapwfr
! !USES:
use modmain
! !DESCRIPTION:
!   Generates the APW radial functions. This is done by integrating the scalar
!   relativistic Schr\"{o}dinger equation (or its energy deriatives) at the
!   current linearisation energies using the spherical part of the Kohn-Sham
!   potential. The number of radial functions at each $l$-value is given by the
!   variable {\tt apword} (at the muffin-tin boundary, the APW functions have
!   continuous derivatives up to order ${\tt apword}-1$). Within each $l$, these
!   functions are orthonormalised with the Gram-Schmidt method. The radial
!   Hamiltonian is applied to the orthonormalised functions and the results are
!   stored in the global array {\tt apwfr}.
!
! !REVISION HISTORY:
!   Created March 2003 (JKD)
!   Copied to equivalent atoms, February 2010 (A. Kozhevnikov and JKD)
!EOP
!BOC
implicit none
! local variables
integer is,ia,ja,ias,jas
integer nr,nri,ir,i
integer nn,l,io,jo
real(8) e,t1
! automatic arrays
logical done(natmmax)
real(8) vr(nrmtmax),fr(nrmtmax)
real(8) p0(nrmtmax,apwordmax),p1(nrmtmax),p1s(apwordmax)
real(8) q0(nrmtmax),q1(nrmtmax),ep0(nrmtmax,apwordmax)
! external functions
real(8) splint
external splint
do is=1,nspecies
  nr=nrmt(is)
  nri=nrmti(is)
  done(:)=.false.
  do ia=1,natoms(is)
    if (done(ia)) cycle
    ias=idxas(ia,is)
! use spherical part of potential
    i=1
    do ir=1,nri
      vr(ir)=vsmt(i,ias)*y00
      i=i+lmmaxi
    end do
    do ir=nri+1,nr
      vr(ir)=vsmt(i,ias)*y00
      i=i+lmmaxo
    end do
    do l=0,lmaxapw
      do io=1,apword(l,is)
! linearisation energy accounting for energy derivative
        e=apwe(io,l,ias)+dble(apwdm(io,l,is))*deapwlo
! integrate the radial Schrodinger equation
        call rschrodint(solsc,l,e,nr,rsp(:,is),vr,nn,p0(:,io),p1,q0,q1)
! multiply by the linearisation energy
        ep0(1:nr,io)=e*p0(1:nr,io)
! normalise radial functions
        fr(1:nr)=p0(1:nr,io)**2
        t1=splint(nr,rsp(:,is),fr)
        t1=1.d0/sqrt(abs(t1))
        call dscal(nr,t1,p0(:,io),1)
        p1s(io)=t1*p1(nr)
        call dscal(nr,t1,ep0(:,io),1)
! subtract linear combination of previous vectors
        do jo=1,io-1
          fr(1:nr)=p0(1:nr,io)*p0(1:nr,jo)
          t1=-splint(nr,rsp(:,is),fr)
          call daxpy(nr,t1,p0(:,jo),1,p0(:,io),1)
          p1s(io)=p1s(io)+t1*p1s(jo)
          call daxpy(nr,t1,ep0(:,jo),1,ep0(:,io),1)
        end do
! normalise radial functions again
        fr(1:nr)=p0(1:nr,io)**2
        t1=splint(nr,rsp(:,is),fr)
        t1=abs(t1)
        if (t1.lt.1.d-25) then
          write(*,*)
          write(*,'("Error(genapwfr): degenerate APW radial functions")')
          write(*,'(" for species ",I4)') is
          write(*,'(" atom ",I4)') ia
          write(*,'(" angular momentum ",I4)') l
          write(*,'(" and order ",I4)') io
          write(*,*)
          stop
        end if
        t1=1.d0/sqrt(t1)
        call dscal(nr,t1,p0(:,io),1)
        p1s(io)=t1*p1s(io)
        call dscal(nr,t1,ep0(:,io),1)
! divide by r and store in global array
        do ir=1,nr
          t1=1.d0/rsp(ir,is)
          apwfr(ir,1,io,l,ias)=t1*p0(ir,io)
          apwfr(ir,2,io,l,ias)=t1*ep0(ir,io)
        end do
! derivative at the muffin-tin surface
        apwdfr(io,l,ias)=(p1s(io)-p0(nr,io)*t1)*t1
      end do
    end do
    done(ia)=.true.
! copy to equivalent atoms
    do ja=1,natoms(is)
      if ((.not.done(ja)).and.(eqatoms(ia,ja,is))) then
        jas=idxas(ja,is)
        do l=0,lmaxapw
          do io=1,apword(l,is)
            call dcopy(nr,apwfr(:,1,io,l,ias),1,apwfr(:,1,io,l,jas),1)
            call dcopy(nr,apwfr(:,2,io,l,ias),1,apwfr(:,2,io,l,jas),1)
            apwdfr(io,l,jas)=apwdfr(io,l,ias)
          end do
        end do
        done(ja)=.true.
      end if
    end do
! end loop over atoms and species
  end do
end do
return
end subroutine
!EOC