Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_tetra_dos
7 : !----------------------------------------------------------------------
8 : !
9 : ! This subroutine evaluates the density of states (g) by the linear
10 : ! tetrahedron method on a energy grid (e) of 'ned' points.
11 : !
12 : ! eig() ... eigenvalues
13 : ! qal() ... partial charges
14 : ! ntet ... number of tetrahedrons
15 : ! ntetra(1-4,ntet)... index of k-points forming tetrahedron nt
16 : ! voltet(ntet) ... volume of tetrahedron nt
17 : !
18 : ! gb 2000
19 : !----------------------------------------------------------------------
20 : USE m_types_kpts
21 :
22 : IMPLICIT NONE
23 :
24 : CONTAINS
25 :
26 68 : SUBROUTINE tetra_dos(jspins,kpts,eGrid,neig,eig,qal,g,energyShift)
27 :
28 : INTEGER, INTENT(IN) :: jspins
29 : TYPE(t_kpts), INTENT(IN) :: kpts
30 : INTEGER, INTENT(IN) :: neig(:,:) !(nkpt,jspins)
31 : REAL, INTENT(IN) :: eGrid(:) !(ned)
32 : REAL, INTENT(IN) :: qal(:,:,:) !(neigd,nkpt,jspins)
33 : REAL, INTENT(IN) :: eig(:,:,:) !(neigd,nkpt,jspins)
34 : REAL, INTENT(OUT) :: g(:,:) !(ned,jspins)
35 : REAL, OPTIONAL, INTENT(IN) :: energyShift
36 :
37 : INTEGER :: i,j,iBand,ikpt,ie,idim,itet,icorn,jcorn,ispin
38 : REAL :: ener,w,sfac, shift, term
39 68 : REAL :: weight(4),eval(4),ecmax(SIZE(eig,1),jspins)
40 68 : REAL, ALLOCATABLE :: wpar(:,:,:)
41 68 : REAL, ALLOCATABLE :: eig_nondeg(:,:,:)
42 :
43 68 : shift = 0.0
44 68 : IF(PRESENT(energyShift)) shift = energyShift
45 :
46 25976 : eig_nondeg = eig - shift
47 :
48 340 : ALLOCATE(wpar,mold=qal)
49 25976 : wpar = 0.0
50 :
51 68 : sfac = 2.0/jspins
52 :
53 136 : DO ispin = 1, SIZE(qal,3)
54 1360 : DO iBand = 1,SIZE(eig,1)
55 1224 : ecmax(iBand,ispin) = -1.0e25
56 25772 : DO ikpt = 1,kpts%nkpt
57 25704 : IF(eig_nondeg(iBand,ikpt,ispin).GT.ecmax(iBand,ispin)) ecmax(iBand,ispin) = eig_nondeg(iBand,ikpt,ispin)
58 : ENDDO
59 : ENDDO
60 : ENDDO
61 : !
62 : ! check for energy degeneracies in tetrahedrons
63 : !
64 136 : DO ispin = 1, SIZE(qal,3)
65 3128 : DO itet = 1,kpts%ntet
66 56916 : DO iBand = 1,SIZE(eig_nondeg,1)
67 218416 : DO i = 1,3
68 161568 : icorn = kpts%ntetra(i,itet)
69 538560 : DO j = i+1,4
70 323136 : jcorn = kpts%ntetra(j,itet)
71 484704 : IF (abs(eig_nondeg(iBand,icorn,ispin)-eig_nondeg(iBand,jcorn,ispin)).LT.1.0e-7) THEN
72 0 : eig_nondeg(iBand,icorn,ispin) = eig_nondeg(iBand,icorn,ispin) + i*1.0e-7*itet
73 0 : eig_nondeg(iBand,jcorn,ispin) = eig_nondeg(iBand,jcorn,ispin) - i*1.0e-7*itet
74 : ENDIF
75 : ENDDO
76 : ENDDO
77 : ENDDO
78 : ENDDO
79 : ENDDO
80 :
81 : !WRITE (*,*) 'in tetra_dos' ! do not remove this statement
82 :
83 : !
84 : ! calculate partial weights
85 : !
86 136 : DO ispin = 1, SIZE(qal,3)
87 1496 : DO ikpt=1,kpts%nkpt
88 25908 : DO iBand = 1,neig(ikpt,ispin)
89 1102960 : DO itet = 1,kpts%ntet
90 4847040 : IF (ALL(kpts%ntetra(:,itet).ne.ikpt)) CYCLE
91 :
92 2154240 : eval = eig_nondeg(iBand,kpts%ntetra(:,itet),ispin)
93 :
94 1077120 : IF(ANY(eval.GE.9999.9)) CYCLE
95 :
96 1101600 : DO i=1,4
97 861696 : icorn = kpts%ntetra(i,itet)
98 :
99 861696 : weight(i)=1.0
100 4308480 : DO j=1,4
101 4308480 : IF (i.NE.j) weight(i)=weight(i)*(eval(j)-eval(i))
102 : ENDDO
103 861696 : weight(i)=6.0*kpts%voltet(itet)/(weight(i)*kpts%ntet)
104 :
105 : wpar(iBand,icorn,ispin) = wpar(iBand,icorn,ispin) &
106 1938816 : + sfac*0.25*weight(i)*qal(iBand,ikpt,ispin)
107 : ENDDO
108 :
109 : ENDDO
110 : ENDDO
111 : ENDDO
112 : ENDDO
113 : !
114 : !---------------------------------------------------
115 : ! evaluate d.o.s.
116 : !---------------------------------------------------
117 : !
118 89964 : g = 0.0
119 :
120 136 : DO ispin = 1, SIZE(qal,3)
121 1496 : DO ikpt = 1,kpts%nkpt
122 25908 : DO iBand = 1,neig(ikpt,ispin)
123 :
124 32362560 : IF(MINVAL(eGrid).GT.ecmax(iband,ispin)) cycle
125 :
126 24480 : ener = eig_nondeg(iBand,ikpt,ispin)
127 24480 : w = 0.5*wpar(iBand,ikpt,ispin)
128 32363920 : DO ie = 1,SIZE(eGrid)
129 32338080 : term = eGrid(ie) - ener
130 32338080 : IF(eGrid(ie).GT.ecmax(iBand,ispin)) cycle
131 29456240 : IF(term.LT.0.0e0) cycle
132 32362560 : g(ie,ispin) = g(ie,ispin) + w * term**2
133 : ENDDO
134 :
135 : ENDDO
136 : ENDDO
137 : ENDDO
138 :
139 68 : END SUBROUTINE tetra_dos
140 : END MODULE m_tetra_dos
|
