Line data Source code
1 : MODULE m_fergwt
2 : USE m_juDFT
3 : !****************************************************************
4 : ! determines the fermi energy and weights for the k-space
5 : ! integration using gaussing-smearing method.
6 : ! c.l.fu
7 : !*****************************************************************
8 : CONTAINS
9 0 : SUBROUTINE fergwt(kpts,input,fmpi, ne,eig, ef,w_iks,seigv,l_output)
10 :
11 : USE m_constants
12 : USE m_types
13 : IMPLICIT NONE
14 :
15 : TYPE(t_mpi),INTENT(IN) :: fmpi
16 : TYPE(t_input),INTENT(IN) :: input
17 : TYPE(t_kpts),INTENT(IN) :: kpts
18 : REAL, INTENT(INOUT) :: ef,seigv
19 : REAL,INTENT(INOUT) :: w_iks(:,:,:)
20 : ! ..
21 : ! ..
22 : ! .. Array Arguments ..
23 : INTEGER, INTENT (IN) :: ne(:,:) !(kpts%nkpt,dimension%jspd)
24 : REAL, INTENT (IN) :: eig(:,:,:) !input%neig,kpts%nkpt,dimension%jspd)
25 : ! .. Logical Arguments ..
26 : LOGICAL,INTENT(IN) :: l_output
27 : ! ..
28 : ! .. Local Scalars ..
29 : REAL chmom,de,ef0,ef1,elow,en,eps,eup,fac,fact1,s,s0,s1,s2,&
30 : workf,wt,wtk,zcdiff,zero,seigv1
31 : INTEGER i,ifl,it,jspin,k,nbnd
32 : ! ..
33 : ! .. External Functions ..
34 : REAL erf
35 : ! ..
36 : ! .. Data statements ..
37 : DATA zero/0.e0/,eps/1.e-5/,eup/3.0e0/,elow/-3.0e0/
38 : ! ..
39 0 : fact1 = input%tkb/SQRT(pi_const)
40 : ! ---> determines ef
41 0 : ifl = 0
42 : conv_loop:DO WHILE (.TRUE.)
43 0 : DO it = 1,50
44 0 : s = 0.
45 0 : DO jspin = 1,input%jspins
46 0 : DO k = 1,kpts%nkpt
47 0 : wtk = kpts%wtkpt(k)
48 0 : nbnd = ne(k,jspin)
49 0 : DO i = 1,nbnd
50 0 : en = eig(i,k,jspin)
51 0 : de = (en-ef)/input%tkb
52 0 : wt = 2.0
53 0 : IF (de.GT.eup) wt = 0.0
54 0 : IF (de.GE.elow .AND. de.LE.eup) THEN
55 0 : IF (de.LT.zero) THEN
56 0 : wt = 1. + ERF(-de)
57 : ELSE
58 0 : wt = 1. - ERF(de)
59 : END IF
60 : END IF
61 0 : s = s + wt*wtk
62 0 : w_iks(i,k,jspin) = wt/2.
63 : ENDDO
64 : ENDDO
65 : ENDDO
66 0 : s = s/REAL(input%jspins)
67 0 : zcdiff = input%zelec - s
68 0 : IF (ABS(zcdiff).LT.eps) EXIT conv_loop
69 0 : IF (ifl.EQ.0) THEN
70 0 : ifl = 1
71 0 : ef0 = ef
72 0 : ef = ef + 0.003
73 0 : s0 = s
74 : ELSE
75 0 : fac = (s0-s)/ (input%zelec-s)
76 0 : IF (ABS(fac).LT.1.0e-1) THEN
77 0 : ef0 = ef
78 0 : s0 = s
79 0 : IF (zcdiff.GE.zero) THEN
80 0 : ef = ef + 0.003
81 : ELSE
82 0 : ef = ef - 0.003
83 : END IF
84 : ELSE
85 0 : ef1 = ef
86 0 : ef = ef + (ef0-ef)/fac
87 0 : ef0 = ef1
88 0 : s0 = s
89 : END IF
90 : END IF
91 : ENDDO
92 0 : eps = 1.25*eps
93 0 : IF ( fmpi%irank == 0 .and. l_output) WRITE (oUnit,FMT=8000) eps
94 : 8000 FORMAT (10x,'warning: eps has been increased to',e12.5)
95 : ENDDO conv_loop
96 0 : workf = -hartree_to_ev_const*ef
97 0 : IF ( fmpi%irank == 0 .and. l_output) THEN
98 0 : WRITE (oUnit,FMT=8010) ef,workf,s
99 : END IF
100 : 8010 FORMAT (/,10x,'fermi energy=',f10.5,' har',3x,'work function=',&
101 : f10.5,' ev',/,10x,'number of valence electrons=',f10.5)
102 0 : IF (ABS(zcdiff).GT.5.0e-4) THEN
103 : CALL juDFT_error('Fermi-level determination did not converge',&
104 0 : hint ="change temperature or set input = F", calledby ="fergwt")
105 : ENDIF
106 0 : DO jspin = 1,input%jspins
107 0 : IF ( fmpi%irank == 0 .and. l_output) WRITE (oUnit,FMT=8020) jspin
108 : 8020 FORMAT (/,/,5x,'band-weighting factor for spin=',i5)
109 0 : DO k = 1,kpts%nkpt
110 0 : nbnd = ne(k,jspin)
111 0 : IF ( fmpi%irank == 0 .and. l_output) WRITE (oUnit,FMT=8030) k
112 : 8030 FORMAT (/,5x,'k-point=',i5,/)
113 0 : w_iks(:,k,jspin) = kpts%wtkpt(k)*w_iks(:,k,jspin)
114 0 : IF ( fmpi%irank == 0 .and. l_output) THEN
115 0 : WRITE (oUnit,FMT=8040) (w_iks(i,k,jspin),i=1,nbnd)
116 : END IF
117 : 8040 FORMAT (5x,16f6.3)
118 : ENDDO
119 : ENDDO
120 0 : s1 = 0.
121 0 : s2 = 0.
122 0 : DO jspin = 1,input%jspins
123 0 : s = 0.
124 0 : DO k = 1,kpts%nkpt
125 0 : DO i = 1,ne(k,jspin)
126 0 : s = s + w_iks(i,k,jspin)
127 0 : seigv = seigv + w_iks(i,k,jspin)*eig(i,k,jspin)
128 0 : en = eig(i,k,jspin)
129 0 : de = (en-ef)/input%tkb
130 : ! ---> correction term
131 0 : IF (ABS(de).LT.3.) THEN
132 0 : de = de*de
133 0 : s2 = s2 + EXP(-de)*kpts%wtkpt(k)
134 : END IF
135 : ENDDO
136 : ENDDO
137 0 : s1 = s1 + s
138 : ENDDO
139 0 : seigv = (2/input%jspins)*seigv
140 0 : seigv1 = (1/input%jspins)*fact1*s2
141 0 : chmom = s1 - input%jspins*s
142 0 : IF ( fmpi%irank == 0 .and. l_output) THEN
143 0 : WRITE (oUnit,FMT=8050) seigv - seigv1,s1,chmom
144 : END IF
145 : 8050 FORMAT (/,10x,'sum of eigenvalues-correction=',f12.5,/,10x,&
146 : 'sum of weight =',f12.5,/,10x,&
147 : 'total moment =',f12.5,/)
148 :
149 0 : END SUBROUTINE fergwt
150 : END MODULE m_fergwt
|