make_tetra.f 2.92 KB
 Markus Betzinger committed Apr 26, 2016 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 `````` MODULE m_maketetra CONTAINS SUBROUTINE make_tetra( > nkpt,bk,ntria,itria,atr, < ntetra,itetra,voltet) C---------------------------------------------------------------- c c Make tetrahedrons out of triangles; assume that layers of c k-points exist, that have been grouped to triangles in c subroutine fertri. Build prisms from them and cut each into c three tetrahedrons. c c bk(1-3,nk) ... coordinates of k-point nk , nk = 1,nkpt c ntria ... number of triangles per layer c itria(1-3,nt) ... index of k-points forming triangle nt c atr(nt) ... area of triangle nt c c ntetra) ... number of tetrahedrons c itetra(1-4,nt)... index of k-points forming tetrahedron nt c voltet(nt) ... volume of tetrahedron nt c omega_bz ... volume of irreducible part of BZ c c Note that the following assumes layers of equally distributed k-points ! c C---------------------------------------------------------------- IMPLICIT NONE INTEGER, INTENT (IN) :: nkpt,ntria INTEGER, INTENT (OUT) :: ntetra INTEGER, INTENT (IN) :: itria(:,:) !(3,2*nkptd) REAL, INTENT (IN) :: atr(:) !(2*nkptd) REAL, INTENT (IN) :: bk(:,:) !(3,nkptd) INTEGER, INTENT (OUT) :: itetra(:,:) !(4,6*nkptd) REAL, INTENT (OUT) :: voltet(:) !(6*nkpt) INTEGER ikpt,nkpp,itri,itet,ip1,ip2,ip3,ip4,ip5,ip6,i,ilay REAL h,h_thrd,tol,omega_bz c c determine distance between planes (h) and number of k-points per plane (nkpp) c tol = 1.0e-15 DO ikpt = 2,nkpt h = abs(bk(3,ikpt)-bk(3,1)) IF (h.GT.tol) EXIT ENDDO nkpp = ikpt - 1 h_thrd = h / 3.0 c c make tetrahedrons c ntetra = 0 DO ilay = 0, (nkpt/nkpp)-2 DO itri = 1,ntria ip1 = itria(1,itri) + nkpp*ilay ; ip4 = ip1 + nkpp ip2 = itria(2,itri) + nkpp*ilay ; ip5 = ip2 + nkpp ip3 = itria(3,itri) + nkpp*ilay ; ip6 = ip3 + nkpp c ntetra = ntetra + 1 itetra(1,ntetra) = ip1 ; itetra(2,ntetra) = ip2 itetra(3,ntetra) = ip3 ; itetra(4,ntetra) = ip4 voltet(ntetra) = h_thrd * atr(itri) c ntetra = ntetra + 1 itetra(1,ntetra) = ip4 ; itetra(2,ntetra) = ip5 itetra(3,ntetra) = ip6 ; itetra(4,ntetra) = ip2 voltet(ntetra) = h_thrd * atr(itri) c ntetra = ntetra + 1 itetra(1,ntetra) = ip2 ; itetra(2,ntetra) = ip3 itetra(3,ntetra) = ip4 ; itetra(4,ntetra) = ip6 voltet(ntetra) = h_thrd * atr(itri) c ENDDO ENDDO omega_bz = 0.0 DO itet = 1,ntetra omega_bz = omega_bz + voltet(itet) ENDDO DO itet = 1,ntetra voltet(itet) = voltet(itet) /omega_bz ENDDO RETURN END SUBROUTINE make_tetra END MODULE m_maketetra``````