Line data Source code
1 : MODULE m_fertetra
2 :
3 : USE m_types
4 : USE m_constants
5 : USE m_juDFT
6 : USE m_tetrahedronInit
7 : USE m_xmlOutput
8 :
9 : IMPLICIT NONE
10 :
11 : CONTAINS
12 :
13 23 : SUBROUTINE fertetra(input,noco,kpts,mpi,ne,eig,ef,w,seigv,l_output)
14 :
15 : TYPE(t_kpts), INTENT(IN) :: kpts
16 : TYPE(t_noco), INTENT(IN) :: noco
17 : TYPE(t_input), INTENT(IN) :: input
18 : TYPE(t_mpi), INTENT(IN) :: mpi
19 : INTEGER, INTENT(IN) :: ne(:,:)
20 : REAL, INTENT(IN) :: eig(:,:,:)
21 : REAL, INTENT(INOUT) :: seigv
22 : REAL, INTENT(INOUT) :: ef
23 : REAL, INTENT(INOUT) :: w(:,:,:)
24 : LOGICAL, INTENT(IN) :: l_output
25 :
26 : INTEGER :: jspin,jspins,ikpt,it,iBand
27 : REAL :: dlow,dup,dfermi,s1,s,chmom,seigvTemp
28 : REAL :: lowBound,upperBound,weightSum
29 : CHARACTER(LEN=20) :: attributes(2)
30 :
31 :
32 23 : jspins = MERGE(1,input%jspins,noco%l_noco)
33 : !---------------------------------------------
34 : !Find the interval, where ef should be located
35 : !---------------------------------------------
36 : !Initial guess
37 35833 : lowBound = MINVAL(eig)-0.01
38 23 : upperBound = ef + 0.2
39 :
40 : !First check the lower bound
41 23 : dlow = 0.0
42 69 : DO jspin = 1, jspins
43 : CALL tetrahedronInit(kpts,input,eig(:,:,jspin),MINVAL(ne(:,jspin)),&
44 2110 : lowBound,weightSum=weightSum)
45 69 : dlow = dlow + weightSum * 2.0/input%jspins
46 : ENDDO
47 :
48 23 : IF (dlow.GT.input%zelec) THEN
49 0 : WRITE(oUnit,9000) lowBound,dlow,input%zelec
50 0 : CALL juDFT_error("valence band too high ",calledby="fertetra")
51 : ENDIF
52 : 9000 FORMAT (' valence band too high ',/,&
53 : ' elow ',f10.5,' dlow ',f10.5,' nelec ',f10.5)
54 :
55 : it = 0
56 : DO
57 : !Now check the upper bound
58 23 : dup = 0.0
59 69 : DO jspin = 1, jspins
60 : CALL tetrahedronInit(kpts,input,eig(:,:,jspin),MINVAL(ne(:,jspin)),&
61 2110 : upperBound,weightSum=weightSum)
62 69 : dup = dup + weightSum * 2.0/input%jspins
63 : ENDDO
64 :
65 23 : IF (dup.GT.input%zelec) THEN
66 : EXIT
67 : ELSE
68 : !Raise the upper bound
69 0 : upperBound = upperBound + 0.2
70 0 : it = it + 1
71 0 : IF(it.GT.10) THEN
72 0 : WRITE (oUnit,9100) upperBound,dup,input%zelec
73 : 9100 FORMAT (' valence band too low ',/,&
74 : ' eup ',f10.5,' dup ',f10.5,' nelec ',f10.5)
75 0 : CALL juDFT_error("valence band too low ",calledby ="fertetra")
76 : ENDIF
77 : ENDIF
78 : ENDDO
79 :
80 : !-----------------------------------------------------------------------------------
81 : !Now that the fermi energy is guaranteed to be in the interval [lowBound,upperBound]
82 : !We use the bisection method to find it
83 : !-----------------------------------------------------------------------------------
84 851 : DO WHILE(upperBound-lowBound.GT.1e-10)
85 :
86 828 : ef = (lowBound+upperBound)/2.0
87 828 : dfermi = 0.0
88 2484 : DO jspin = 1, jspins
89 : !-------------------------------------------------------
90 : ! Compute the current occupation
91 : !-------------------------------------------------------
92 : CALL tetrahedronInit(kpts,input,eig(:,:,jspin),MINVAL(ne(:,jspin)),&
93 75960 : ef,weightSum=weightSum)
94 2484 : dfermi = dfermi + weightSum * 2.0/input%jspins
95 : ENDDO
96 851 : IF(ABS(dfermi-input%zelec).LT.1e-12) THEN
97 : EXIT
98 828 : ELSE IF(dfermi-input%zelec.GT.0.0) THEN
99 : !Occupation to large -> search in the left interval
100 423 : upperBound = ef
101 : ELSE
102 : !Occupation to small -> search in the right interval
103 405 : lowBound = ef
104 : ENDIF
105 : ENDDO
106 :
107 : !---------------------------------------------------------------------
108 : !Calculate final occupation and weights for individual kpts and bands
109 : !---------------------------------------------------------------------
110 23 : ef = (lowBound+upperBound)/2.0
111 23 : dfermi = 0.0
112 69 : DO jspin = 1, jspins
113 : !-------------------------------------------------------
114 : ! Compute the weights for charge density integration
115 : !-------------------------------------------------------
116 : CALL tetrahedronInit(kpts,input,eig(:,:,jspin),MINVAL(ne(:,jspin)),&
117 2110 : ef,weightSum=weightSum,weights=w(:,:,jspin))
118 69 : dfermi = dfermi + weightSum * 2.0/input%jspins
119 : ENDDO
120 :
121 23 : IF (l_output) THEN
122 23 : WRITE (oUnit,9200) ef,dfermi,input%zelec
123 : 9200 FORMAT (//,'Tetrahedron method: ',//,' fermi energy =',f10.5,&
124 : ' dtot ',f10.5,' nelec ',f10.5)
125 : ENDIF
126 :
127 : !----------------------------------------------
128 : ! Obtain sum of weights and valence eigenvalues
129 : !----------------------------------------------
130 23 : s1 = 0.0
131 23 : seigv = 0.0
132 69 : DO jspin = 1,jspins
133 46 : s = 0.0
134 2110 : DO ikpt = 1,kpts%nkpt
135 34662 : DO iBand = 1,ne(ikpt,jspin)
136 32552 : s = s + w(iBand,ikpt,jspin)
137 34616 : seigv = seigv + w(iBand,ikpt,jspin)*eig(iBand,ikpt,jspin)
138 : ENDDO
139 : ENDDO
140 69 : s1 = s1 + s
141 : ENDDO
142 23 : seigv = 2.0/input%jspins*seigv
143 23 : chmom = s1 - jspins*s
144 :
145 23 : seigvTemp = seigv
146 23 : IF (noco%l_soc .AND. (.NOT. noco%l_noco)) THEN
147 0 : seigvTemp = seigvTemp / 2.0
148 : END IF
149 :
150 23 : IF ( mpi%irank == 0 .AND. l_output) THEN
151 69 : attributes = ''
152 23 : WRITE(attributes(1),'(f20.10)') seigvTemp
153 23 : WRITE(attributes(2),'(a)') 'Htr'
154 69 : CALL writeXMLElement('sumValenceSingleParticleEnergies',(/'value','units'/),attributes)
155 23 : WRITE (oUnit,FMT=9300) seigvTemp,s1,chmom
156 : END IF
157 : 9300 FORMAT (/,10x,'sum of valence eigenvalues=',f20.10,5x,&
158 : 'sum of weights=',f10.6,/,10x,'moment=',f12.6)
159 :
160 23 : END SUBROUTINE fertetra
161 : END MODULE m_fertetra
|
