fertetra.f90 5.68 KB
Newer Older
1
2
3
MODULE m_fertetra

   USE m_types
4
   USE m_constants
5
6
   USE m_juDFT
   USE m_tetrahedronInit
7
   USE m_xmlOutput
8
9
10
11
12

   IMPLICIT NONE

   CONTAINS

Henning Janssen's avatar
Henning Janssen committed
13
   SUBROUTINE fertetra(input,noco,kpts,mpi,ne,eig,ef,w,seigv)
14
15

      TYPE(t_kpts),        INTENT(IN)     :: kpts
Henning Janssen's avatar
Henning Janssen committed
16
      TYPE(t_noco),        INTENT(IN)     :: noco
17
18
19
20
21
22
23
24
      TYPE(t_input),       INTENT(IN)     :: input
      TYPE(t_mpi),         INTENT(IN)     :: mpi
      INTEGER,             INTENT(IN)     :: ne(:,:)
      REAL,                INTENT(IN)     :: eig(:,:,:)
      REAL,                INTENT(INOUT)  :: seigv
      REAL,                INTENT(INOUT)  :: ef
      REAL,                INTENT(INOUT)  :: w(:,:,:)

Henning Janssen's avatar
Henning Janssen committed
25
      INTEGER :: jspin,jspins,ikpt,it,iBand
26
      REAL    :: dlow,dup,dfermi,s1,s,chmom,seigvTemp
27
      REAL    :: lowBound,upperBound,weightSum
28
      CHARACTER(LEN=20)    :: attributes(2)
29

Henning Janssen's avatar
Henning Janssen committed
30
31

      jspins = MERGE(1,input%jspins,noco%l_noco)
32
33
34
35
36
37
38
39
40
      !---------------------------------------------
      !Find the interval, where ef should be located
      !---------------------------------------------
      !Initial guess
      lowBound = MINVAL(eig)-0.01
      upperBound = ef + 0.2

      !First check the lower bound
      dlow = 0.0
Henning Janssen's avatar
Henning Janssen committed
41
      DO jspin = 1, jspins
42
43
         CALL tetrahedronInit(kpts,input,eig(:,:,jspin),MINVAL(ne(:,jspin)),&
                              lowBound,weightSum=weightSum)
44
         dlow = dlow + weightSum * 2.0/input%jspins
45
46
47
      ENDDO

      IF (dlow.GT.input%zelec) THEN
48
         WRITE(oUnit,9000) lowBound,dlow,input%zelec
49
50
51
52
53
54
55
56
57
         CALL juDFT_error("valence band too high ",calledby="fertetra")
      ENDIF
9000  FORMAT (' valence band too high ',/,&
              '  elow ',f10.5,' dlow ',f10.5,' nelec ',f10.5)

      it = 0
      DO
         !Now check the upper bound
         dup = 0.0
Henning Janssen's avatar
Henning Janssen committed
58
         DO jspin = 1, jspins
59
60
            CALL tetrahedronInit(kpts,input,eig(:,:,jspin),MINVAL(ne(:,jspin)),&
                                 upperBound,weightSum=weightSum)
61
            dup = dup + weightSum * 2.0/input%jspins
62
63
64
65
66
67
68
69
70
         ENDDO

         IF (dup.GT.input%zelec) THEN
            EXIT
         ELSE
            !Raise the upper bound
            upperBound = upperBound + 0.2
            it = it + 1
            IF(it.GT.10) THEN
71
               WRITE (oUnit,9100) upperBound,dup,input%zelec
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
9100           FORMAT (' valence band too low ',/,&
                       '  eup  ',f10.5,' dup  ',f10.5,' nelec ',f10.5)
               CALL juDFT_error("valence band too low ",calledby ="fertetra")
            ENDIF
         ENDIF
      ENDDO

      !-----------------------------------------------------------------------------------
      !Now that the fermi energy is guaranteed to be in the interval [lowBound,upperBound]
      !We use the bisection method to find it
      !-----------------------------------------------------------------------------------
      DO WHILE(upperBound-lowBound.GT.1e-10)

         ef = (lowBound+upperBound)/2.0
         dfermi = 0.0
Henning Janssen's avatar
Henning Janssen committed
87
         DO jspin = 1, jspins
88
            !-------------------------------------------------------
89
            ! Compute the current occupation
90
            !-------------------------------------------------------
91
92
            CALL tetrahedronInit(kpts,input,eig(:,:,jspin),MINVAL(ne(:,jspin)),&
                                 ef,weightSum=weightSum)
93
            dfermi = dfermi + weightSum * 2.0/input%jspins
94
         ENDDO
95
96
97
98
         IF(ABS(dfermi-input%zelec).LT.1e-12) THEN
            EXIT
         ELSE IF(dfermi-input%zelec.GT.0.0) THEN
            !Occupation to large -> search in the left interval
99
            upperBound = ef
100
         ELSE
101
102
103
104
            !Occupation to small -> search in the right interval
            lowBound = ef
         ENDIF
      ENDDO
105
106
107
108

      !---------------------------------------------------------------------
      !Calculate final occupation and weights for individual kpts and bands
      !---------------------------------------------------------------------
109
      ef = (lowBound+upperBound)/2.0
110
111
112
113
114
      dfermi = 0.0
      DO jspin = 1, jspins
         !-------------------------------------------------------
         ! Compute the weights for charge density integration
         !-------------------------------------------------------
115
116
         CALL tetrahedronInit(kpts,input,eig(:,:,jspin),MINVAL(ne(:,jspin)),&
                              ef,weightSum=weightSum,weights=w(:,:,jspin))
117
118
119
         dfermi = dfermi + weightSum * 2.0/input%jspins
      ENDDO

120
      WRITE (oUnit,9200) ef,dfermi,input%zelec
121
9200  FORMAT (//,'Tetrahedron method: ',//,'   fermi energy =',f10.5,&
122
123
124
125
126
                 ' dtot ',f10.5,' nelec ',f10.5)

      !----------------------------------------------
      ! Obtain sum of weights and valence eigenvalues
      !----------------------------------------------
127
128
      s1 = 0.0
      seigv = 0.0
Henning Janssen's avatar
Henning Janssen committed
129
      DO jspin = 1,jspins
130
         s = 0.0
131
132
         DO ikpt = 1,kpts%nkpt
            DO iBand = 1,ne(ikpt,jspin)
133
134
135
136
137
138
139
               s = s + w(iBand,ikpt,jspin)
               seigv = seigv + w(iBand,ikpt,jspin)*eig(iBand,ikpt,jspin)
            ENDDO
         ENDDO
         s1 = s1 + s
      ENDDO
      seigv = 2.0/input%jspins*seigv
Henning Janssen's avatar
Henning Janssen committed
140
      chmom = s1 - jspins*s
141
142
143
144
145
146
      
      seigvTemp = seigv
      IF (noco%l_soc .AND. (.NOT. noco%l_noco)) THEN
         seigvTemp = seigvTemp / 2.0
      END IF
      
147
      IF ( mpi%irank == 0 ) THEN
148
149
150
151
152
         attributes = ''
         WRITE(attributes(1),'(f20.10)') seigvTemp
         WRITE(attributes(2),'(a)') 'Htr'
         CALL writeXMLElement('sumValenceSingleParticleEnergies',(/'value','units'/),attributes)
         WRITE (oUnit,FMT=9300) seigvTemp,s1,chmom
153
      END IF
154
9300  FORMAT (/,10x,'sum of valence eigenvalues=',f20.10,5x,&
155
156
157
             'sum of weights=',f10.6,/,10x,'moment=',f12.6)

   END SUBROUTINE fertetra
158
END MODULE m_fertetra