File: gencore.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 (148 lines) | stat: -rw-r--r-- 4,484 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

! 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: gencore
! !INTERFACE:
subroutine gencore
! !USES:
use modmain
use modomp
! !DESCRIPTION:
!   Computes the core radial wavefunctions, eigenvalues and densities. The
!   radial Dirac equation is solved in the spherical part of the Kohn-Sham
!   potential to which the atomic potential has been appended for
!   $r>R_{\rm MT}$. In the case of spin-polarised calculations, and when
!   {\tt spincore} is set to {\tt .true.}, the Dirac equation is solved in the
!   spin-up and -down potentials created from the Kohn-Sham scalar potential and
!   magnetic field magnitude, with the occupancy divided equally between up and
!   down. The up and down densities determined in this way are added to both the
!   scalar density and the magnetisation in the routine {\tt rhocore}. Note
!   that this procedure is a simple, but inexact, approach to solving the radial
!   Dirac equation in a magnetic field.
!
! !REVISION HISTORY:
!   Created April 2003 (JKD)
!   Added polarised cores, November 2009 (JKD)
!EOP
!BOC
implicit none
! local variables
integer ist,ispn,idm
integer is,ia,ja,ias,jas
integer nr,nri,nrs,ir,nthd
real(8) t1
! automatic arrays
logical done(natmmax)
real(8) vr(nrspmax),eval(nstspmax)
! allocatable arrays
real(8), allocatable :: br(:),fr(:,:)
if (spincore) allocate(br(nrmtmax),fr(nrmtmax,ndmag))
! loop over species and atoms
do is=1,nspecies
  nr=nrmt(is)
  nri=nrmti(is)
  nrs=nrsp(is)
  done(:)=.false.
  do ia=1,natoms(is)
    if (done(ia)) cycle
    ias=idxas(ia,is)
! Kohn-Sham magnetic field for spin-polarised core
    if (spincore) then
      do idm=1,ndmag
        call rfmtlm(1,nr,nri,bxcmt(:,ias,idm),fr(:,idm))
      end do
      if (ncmag) then
        do ir=1,nr
          br(ir)=sqrt(fr(ir,1)**2+fr(ir,2)**2+fr(ir,3)**2)*y00
        end do
      else
        do ir=1,nr
          br(ir)=abs(fr(ir,1))*y00
        end do
      end if
    end if
! loop over spin channels
    do ispn=1,nspncr
! use the spherical part of the crystal Kohn-Sham potential
      call rfmtlm(1,nr,nri,vsmt(:,ias),vr)
      vr(1:nr)=vr(1:nr)*y00
! spin-up and -down potentials for polarised core
      if (spincore) then
        if (ispn.eq.1) then
          vr(1:nr)=vr(1:nr)+br(1:nr)
        else
          vr(1:nr)=vr(1:nr)-br(1:nr)
        end if
      end if
! append the Kohn-Sham potential from the atomic calculation for r > R_MT
      t1=vr(nr)-vrsp(nr,is)
      do ir=nr+1,nrs
        vr(ir)=vrsp(ir,is)+t1
      end do
      rhocr(:,ias,ispn)=0.d0
      call omp_hold(nstsp(is),nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(t1,ir) SHARED(is) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
      do ist=1,nstsp(is)
        if (spcore(ist,is)) then
! solve the Dirac equation
          eval(ist)=evalcr(ist,ias)
          call rdirac(solsc,nsp(ist,is),lsp(ist,is),ksp(ist,is),nrs,rsp(:,is), &
           vr,eval(ist),rwfcr(:,1,ist,ias),rwfcr(:,2,ist,ias))
          if (spincore) then
! use the spin-averaged eigenvalue for the polarised core
            if (ispn.eq.1) then
              evalcr(ist,ias)=eval(ist)
            else
              evalcr(ist,ias)=0.5d0*(evalcr(ist,ias)+eval(ist))
            end if
            t1=0.5d0*occcr(ist,ias)
          else
            evalcr(ist,ias)=eval(ist)
            t1=occcr(ist,ias)
          end if
! add to the core density
!$OMP CRITICAL(gencore_)
          do ir=1,nrs
            rhocr(ir,ias,ispn)=rhocr(ir,ias,ispn) &
             +t1*(rwfcr(ir,1,ist,ias)**2+rwfcr(ir,2,ist,ias)**2)
          end do
!$OMP END CRITICAL(gencore_)
        end if
      end do
!$OMP END DO
!$OMP END PARALLEL
      call omp_free(nthd)
      do ir=1,nrs
        rhocr(ir,ias,ispn)=rhocr(ir,ias,ispn)*y00/r2sp(ir,is)
      end do
! end loop over spin channels
    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 ist=1,nstsp(is)
          if (spcore(ist,is)) then
            evalcr(ist,jas)=evalcr(ist,ias)
            rwfcr(1:nrs,:,ist,jas)=rwfcr(1:nrs,:,ist,ias)
          end if
        end do
        rhocr(1:nrs,jas,:)=rhocr(1:nrs,ias,:)
        done(ja)=.true.
      end if
    end do
! end loop over species and atoms
  end do
end do
if (spincore) deallocate(br,fr)
return
end subroutine
!EOC