-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathrelnodes.f90
135 lines (116 loc) · 3.9 KB
/
relnodes.f90
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
! The code was developed at the Fritz Haber Institute, and
! the intellectual properties and copyright of this file
! are with the Max Planck Society. When you use it, please
! cite R. Gomez-Abal, X. Li, C. Ambrosch-Draxl, M. Scheffler,
! Extended linear tetrahedron method for the calculation of q-dependent
! dynamical response functions, to be published in Comp. Phys. Commun. (2010)
!BOP
!
! !ROUTINE: relnodes
!
! !INTERFACE:
subroutine relnodes
!
! !DESCRIPTION:
!
! This subroutine select the relevant intnodes, that is, only those defining
! the polyhedron that has to be integrated. The intnodes selected are those for
! which $\varepsilon(\vec{k})\le \varepsilon_F$ and
! $\varepsilon'(\vec{k})\ge \varepsilon_F$ if ndabc=1. $\varepsilon(\vec{k})
! \le \varepsilon_F$ when ndabc=2. $\varepsilon(\vec{k})\le \varepsilon_F$ and
! $\varepsilon'(\vec{k})\le \varepsilon_F$ when ndabc=3. So the first region equals
! the second region minused by the third region.
!
! !USES:
use polyhedron
implicit none
! !LOCAL VARIABLES:
integer(4) :: inod, jnod, maxnod
real(8) :: e1, e2, zerotol
real(8), dimension(3) :: ncor
real(8), external :: tlinap
! !REVISION HISTORY:
!
! Created 22nd. April 2004 by RGA
! Last revised 14th. Dec 2004 by XZL
!
!EOP
!BOC
zerotol=1.00d-15
inod=1
maxnod=nnod
select case(ndabc)
case(1)
!-----------------------------------------------------------------------------------------------
! in this case, we select the nodes to get the region when e1<efer and e2>efer
!-----------------------------------------------------------------------------------------------
inod=1
maxnod=nnod
do while (inod.le.maxnod)
ncor(1:3)=intnodes(1:3,inod)
e1=tlinap(ncor,e)
e2=tlinap(ncor,f)
if(((e1.gt.ef).and.(dabs(e1-ef).gt.zerotol)).or.((e2.lt.ef).and.(dabs(e2-ef).gt.zerotol)))then
do jnod=inod+1,maxnod
intnodes(1:3,jnod-1)=intnodes(1:3,jnod)
ntype(jnod-1)=ntype(jnod)
enddo
maxnod=maxnod-1
else
inod=inod+1
endif
enddo
do inod=maxnod+1,nnod
intnodes(1:3,inod)=0.0d0
ntype(inod)=0
enddo
case(2)
!-----------------------------------------------------------------------------------------------
! in this case, we select the nodes to get the region when e1<efer
!-----------------------------------------------------------------------------------------------
do while (inod.le.maxnod)
ncor(1:3)=intnodes(1:3,inod)
e1=tlinap(ncor,e)
e2=tlinap(ncor,f)
if((e1.gt.ef).and.(dabs(e1-ef).gt.zerotol))then
do jnod=inod+1,maxnod
intnodes(1:3,jnod-1)=intnodes(1:3,jnod)
ntype(jnod-1)=ntype(jnod)
enddo
maxnod=maxnod-1
else
inod=inod+1
endif
enddo
do inod=maxnod+1,nnod
intnodes(1:3,inod)=0.0d0
ntype(inod)=0
enddo
case(3)
!-----------------------------------------------------------------------------------------------
! In this case, we get the region when e1<efer and e1<efer
!-----------------------------------------------------------------------------------------------
inod=1
maxnod=nnod
do while (inod.le.maxnod)
ncor(1:3)=intnodes(1:3,inod)
e1=tlinap(ncor,e)
e2=tlinap(ncor,f)
if(((e1.gt.ef).and.(dabs(e1-ef).gt.zerotol)).or.((e2.gt.ef).and.(dabs(e2-ef).gt.zerotol)))then
do jnod=inod+1,maxnod
intnodes(1:3,jnod-1)=intnodes(1:3,jnod)
ntype(jnod-1)=ntype(jnod)
enddo
maxnod=maxnod-1
else
inod=inod+1
endif
enddo
do inod=maxnod+1,nnod
intnodes(1:3,inod)=0.0d0
ntype(inod)=0
enddo
end select
nnod=maxnod
end subroutine relnodes
!EOC