File: genlofr.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 (156 lines) | stat: -rw-r--r-- 4,759 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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156

! 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: genlofr
! !INTERFACE:
subroutine genlofr
! !USES:
use modmain
! !DESCRIPTION:
!   Generates the local-orbital 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. For each local-orbital, a linear combination of
!   {\tt lorbord} radial functions is constructed such that its radial
!   derivatives up to order ${\tt lorbord}-1$ are zero at the muffin-tin radius.
!   This function is normalised and the radial Hamiltonian applied to it. The
!   results are stored in the global array {\tt lofr}.
!
! !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 ilo,jlo,io,jo
integer nn,l,info
real(8) e,t1
! automatic arrays
logical done(natmmax)
integer ipiv(nplorb)
real(8) vr(nrmtmax),fr(nrmtmax)
real(8) p0(nrmtmax,lorbordmax),p1(nrmtmax)
real(8) q0(nrmtmax),q1(nrmtmax),ep0(nrmtmax,lorbordmax)
real(8) p0s(nrmtmax,nlomax),ep0s(nrmtmax,nlomax)
real(8) xa(nplorb),ya(nplorb)
real(8) a(nplorb,nplorb),b(nplorb),c(nplorb)
! external functions
real(8) splint,polynom
external splint,polynom
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 ilo=1,nlorb(is)
      l=lorbl(ilo,is)
      do jo=1,lorbord(ilo,is)
! linearisation energy accounting for energy derivative
        e=lorbe(jo,ilo,ias)+dble(lorbdm(jo,ilo,is))*deapwlo
! integrate the radial Schrodinger equation
        call rschrodint(solsc,l,e,nr,rsp(:,is),vr,nn,p0(:,jo),p1,q0,q1)
        ep0(1:nr,jo)=e*p0(1:nr,jo)
! normalise radial functions
        fr(1:nr)=p0(1:nr,jo)**2
        t1=splint(nr,rsp(:,is),fr)
        t1=1.d0/sqrt(abs(t1))
        call dscal(nr,t1,p0(:,jo),1)
        call dscal(nr,t1,ep0(:,jo),1)
! set up the matrix of radial derivatives
        do i=1,nplorb
          ir=nr-nplorb+i
          xa(i)=rsp(ir,is)
          ya(i)=p0(ir,jo)/rsp(ir,is)
        end do
        do io=1,lorbord(ilo,is)
          a(io,jo)=polynom(io-1,nplorb,xa,ya,c,rmt(is))
        end do
      end do
! set up the target vector
      b(:)=0.d0
      b(lorbord(ilo,is))=1.d0
      call dgesv(lorbord(ilo,is),1,a,nplorb,ipiv,b,nplorb,info)
      if (info.ne.0) goto 10
! generate linear superposition of radial functions
      p0s(:,ilo)=0.d0
      ep0s(:,ilo)=0.d0
      do io=1,lorbord(ilo,is)
        t1=b(io)
        call daxpy(nr,t1,p0(:,io),1,p0s(:,ilo),1)
        call daxpy(nr,t1,ep0(:,io),1,ep0s(:,ilo),1)
      end do
! normalise radial functions
      fr(1:nr)=p0s(1:nr,ilo)**2
      t1=splint(nr,rsp(:,is),fr)
      t1=1.d0/sqrt(abs(t1))
      call dscal(nr,t1,p0s(:,ilo),1)
      call dscal(nr,t1,ep0s(:,ilo),1)
! subtract linear combination of previous local-orbitals with same l
      do jlo=1,ilo-1
        if (lorbl(jlo,is).eq.l) then
          fr(1:nr)=p0s(1:nr,ilo)*p0s(1:nr,jlo)
          t1=-splint(nr,rsp(:,is),fr)
          call daxpy(nr,t1,p0s(:,jlo),1,p0s(:,ilo),1)
          call daxpy(nr,t1,ep0s(:,jlo),1,ep0s(:,ilo),1)
        end if
      end do
! normalise radial functions again
      fr(1:nr)=p0s(1:nr,ilo)**2
      t1=splint(nr,rsp(:,is),fr)
      t1=abs(t1)
      if (t1.lt.1.d-25) goto 10
      t1=1.d0/sqrt(t1)
      call dscal(nr,t1,p0s(:,ilo),1)
      call dscal(nr,t1,ep0s(:,ilo),1)
! divide by r and store in global array
      do ir=1,nr
        t1=1.d0/rsp(ir,is)
        lofr(ir,1,ilo,ias)=t1*p0s(ir,ilo)
        lofr(ir,2,ilo,ias)=t1*ep0s(ir,ilo)
      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 ilo=1,nlorb(is)
          call dcopy(nr,lofr(:,1,ilo,ias),1,lofr(:,1,ilo,jas),1)
          call dcopy(nr,lofr(:,2,ilo,ias),1,lofr(:,2,ilo,jas),1)
        end do
        done(ja)=.true.
      end if
    end do
! end loop over atoms and species
  end do
end do
return
10 continue
write(*,*)
write(*,'("Error(genlofr): degenerate local-orbital radial functions")')
write(*,'(" for species ",I4)') is
write(*,'(" atom ",I4)') ia
write(*,'(" and local-orbital ",I4)') ilo
write(*,*)
stop
end subroutine
!EOC