LCOV - code coverage report
Current view: top level - global - triang.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 101 126 80.2 %
Date: 2019-09-08 04:53:50 Functions: 1 3 33.3 %

          Line data    Source code
       1             :       MODULE m_triang
       2             :       use m_juDFT
       3             : !-------------------------------------------------------------------
       4             : c     find a triangular decomposition of the irreducible wedge of
       5             : c     the first brillouin zone for a given k-mesh. k-points at all
       6             : c     vertices are assumed.
       7             : c     erich wimmer     july 1981
       8             : !-------------------------------------------------------------------
       9             : 
      10             :       IMPLICIT NONE
      11             : 
      12             :       CONTAINS
      13           7 :       SUBROUTINE triang(
      14           7 :      >                  v,nkpt,
      15           7 :      <                  it,ntria,at,att,l_f_t)
      16             : 
      17             : c     Arguments
      18             :       INTEGER, INTENT(IN)  :: nkpt
      19             :       REAL,    INTENT(IN)  :: v(3,nkpt)
      20             :       REAL,    INTENT(OUT) :: att , at(2*nkpt)
      21             :       INTEGER, INTENT(OUT) :: ntria , it(3,2*nkpt)
      22             :       LOGICAL, INTENT(INOUT) :: l_f_t
      23             : c     locals
      24             :       REAL    :: a, a1, d, dm, s0, s1, x1, x2, x3, y1, y2, y3
      25             :       INTEGER :: j , j1 , j2 , jc , jj , k , k1 , k2 , k3 , kk , l1 , 
      26             :      &           l2 , l3 , n , nt0 , nkpts
      27             :       LOGICAL :: new
      28             : c     constants
      29             :       REAL , PARAMETER :: zero = 0.0 , big = 1.e8 , tol = 1.e-5
      30             : 
      31           7 :       ntria = 0
      32           7 :       IF ( nkpt.LT.3 ) RETURN
      33             : c
      34             : c l_f_t=.true. means that we call from fertri and on output gives 'film'
      35             : c
      36           7 :       IF (l_f_t) THEN ! check for bulk
      37         133 :          DO j= 2,nkpt
      38          70 :             IF (abs(v(3,j)-v(3,1)).GT.1.0e-12) THEN
      39           4 :               l_f_t = .false.
      40           4 :               nkpts = j - 1
      41           4 :               IF ((mod(nkpt,nkpts).NE.0).OR.(j.LT.3)) THEN
      42           3 :                 WRITE (6,*) 'tria=T & film=F requires k-point planes'
      43           3 :                 WRITE (6,*) 'with equally distributed k-points !'
      44             : c                 CALL juDFT_error("not a k-point set for tria=T & film=F",calledby="triang")
      45             :               ENDIF
      46             : !              RETURN
      47             :               GOTO 10
      48             :             ENDIF
      49             :          ENDDO
      50           3 :          l_f_t = .true.
      51           3 :          nkpts = nkpt
      52             :  10      CONTINUE
      53             :       ELSE
      54             : c
      55             : c l_f_t=.false. means that we call from evaldos
      56             : c
      57             :          nkpts = nkpt
      58           0 :          DO j= 2,nkpt
      59           0 :             IF (abs(v(3,j)-v(3,1)).GT.1.0e-12) THEN
      60           0 :               nkpts = j - 1
      61           0 :               GOTO 11
      62             :            ENDIF
      63             :         ENDDO
      64             :  11     CONTINUE
      65           0 :         IF ((mod(nkpt,nkpts).NE.0).OR.(j.LT.3)) THEN
      66           0 :           l_f_t = .false.
      67           0 :           RETURN
      68             :         ELSE
      69           0 :           l_f_t = .true.
      70             :         ENDIF
      71             :       ENDIF
      72             : c
      73             : c     create first triangle
      74             : c
      75           7 :       att = zero
      76           7 :       k1 = 1
      77           7 :       dm = big
      78          70 :       DO k = 2 , nkpts
      79         126 :          d = vd2(v(1,K1),v(2,k1),v(1,k),v(2,k))
      80          70 :          IF ( d.LT.dm ) THEN
      81           5 :             dm = d
      82           5 :             k2 = k
      83             :          ENDIF
      84             :       ENDDO
      85             : c
      86             :       dm = big
      87             :       new = .false.
      88         147 :       DO k = 1 , nkpts
      89          77 :          IF ( k.NE.k1 .AND. k.NE.k2 ) THEN
      90             :             d = VD2(V(1,k1),V(2,k1),V(1,k),V(2,k))
      91         174 :      &          + VD2(V(1,k2),V(2,k2),V(1,k),V(2,k))
      92         116 :             a = AREA(V(1,k1),V(2,k1),V(1,k2),V(2,k2),V(1,k),V(2,k))
      93          58 :             IF ( ABS(a).GE.TOL ) THEN
      94           3 :                IF ( d.LT.dm ) THEN
      95           2 :                   new = .TRUE.
      96           2 :                   dm = d
      97           2 :                   a1 = a
      98           2 :                   k3 = k
      99             :                ENDIF
     100             :             ENDIF
     101             :          ENDIF
     102             :       ENDDO
     103           7 :       IF ( new ) THEN
     104           2 :          Ntria = 1
     105           2 :          It(1,Ntria) = k1
     106           2 :          IF ( k3.LE.k2 ) THEN
     107           0 :             kk = k2
     108           0 :             k2 = k3
     109           0 :             k3 = kk
     110             :          ENDIF
     111           2 :          It(2,Ntria) = k2
     112           2 :          It(3,Ntria) = k3
     113           2 :          At(Ntria) = ABS(a1)/2
     114           2 :          Att = Att + At(Ntria)
     115             : c
     116             : c     create, if possible, a new triangle from each side of the mother
     117             : c     triangle nt0 with minimal sum of sides
     118             : c
     119           2 :          nt0 = 0
     120             :          DO WHILE ( .TRUE. )
     121           7 :             nt0 = nt0 + 1
     122           7 :             IF ( nt0.GT.Ntria ) RETURN
     123          35 :             DO l1 = 1 , 3
     124          15 :                l2 = MOD(l1,3) + 1
     125          15 :                l3 = MOD(l2,3) + 1
     126          15 :                k1 = It(l1,nt0)
     127          15 :                k2 = It(l2,nt0)
     128          15 :                IF ( k2.LE.k1 ) THEN
     129           5 :                   kk = k1
     130           5 :                   k1 = k2
     131           5 :                   k2 = kk
     132             :                ENDIF
     133             : c--->>     check if side k1-k2 belongs already to another triangle
     134          15 :                new = .TRUE.
     135          58 :                DO n = 1 , Ntria
     136          58 :                   IF ( n.NE.nt0 ) THEN
     137             :                      IF ( (k1.EQ.It(1,n) .AND. k2.EQ.It(2,n)) .OR. 
     138          28 :      &                    (k1.EQ.It(1,n) .AND. k2.EQ.It(3,n)) .OR. 
     139             :      &                    (k1.EQ.It(2,n) .AND. k2.EQ.It(3,n)) )
     140             :      &                    new = .FALSE.
     141             :                   ENDIF
     142             :                ENDDO
     143          20 :                IF ( new ) THEN
     144          10 :                   k3 = It(l3,nt0)
     145             :                   a = AREA(V(1,k1),V(2,k1),V(1,k2),V(2,k2),V(1,k3),
     146          20 :      &                V(2,k3))
     147          10 :                   s0 = SIGN(1.,a)
     148             : c--->>     a new triangle sharing the side k1-k2 with the mother
     149             : c--->>     triangle nt0 has to ly on the other side, i.e. the cross
     150             : c--->>     products (k2-k1)x(k3(old)-k1) and (k2-k1)x(k3(new)-k1)
     151             : c--->>     have to have opposite signs
     152          10 :                   dm = BIG
     153          10 :                   new = .FALSE.
     154          54 :                   DO k = 1 , Nkpts
     155          44 :                      IF ( k.NE.k1 .AND. k.NE.k2 ) THEN
     156             : c--->>     check if a new side, (k1,k) or (k2,k), belongs
     157             : c--->>     already to an older triangle
     158             :                         j1 = k1
     159             :                         j2 = k
     160          92 :                         DO j = 1 , 2
     161          41 :                            IF ( j2.LE.j1 ) THEN
     162          15 :                               jj = j1
     163          15 :                               j1 = j2
     164          15 :                               j2 = jj
     165             :                            ENDIF
     166          41 :                            jc = 0
     167         150 :                            DO n = 1 , Ntria
     168             :                               IF ( j1.EQ.It(1,n) .AND. j2.EQ.It(2,n)
     169             :      &                             .OR. j1.EQ.It(1,n) .AND. 
     170         109 :      &                             j2.EQ.It(3,n) .OR. j1.EQ.It(2,n)
     171          73 :      &                             .AND. j2.EQ.It(3,n) ) jc = jc + 1
     172             :                            ENDDO
     173          41 :                            IF ( jc.EQ.2 ) GOTO 5
     174          34 :                            j1 = k2
     175          51 :                            j2 = k
     176             :                         ENDDO
     177             :                         a = AREA(V(1,k1),V(2,k1),V(1,k2),V(2,k2),V(1,k),
     178          34 :      &                      V(2,k))
     179          17 :                         IF ( ABS(a).GE.TOL ) THEN
     180          16 :                            s1 = SIGN(1.,a)
     181          16 :                            IF ( ABS(s1-s0).GE.TOL ) THEN
     182             :                               d = VD2(V(1,k1),V(2,k1),V(1,k),V(2,k))
     183           9 :      &                            + VD2(V(1,k2),V(2,k2),V(1,k),V(2,k))
     184           3 :                               IF ( d.LT.dm ) THEN
     185           3 :                                  new = .TRUE.
     186           3 :                                  dm = d
     187           3 :                                  a1 = a
     188           3 :                                  k3 = k
     189             :                               ENDIF
     190             :                            ENDIF
     191             :                         ENDIF
     192             :                      ENDIF
     193          10 :  5                ENDDO
     194          10 :                   IF ( new ) THEN
     195           3 :                      Ntria = Ntria + 1
     196           3 :                      IF ( Ntria>2*nkpt )  CALL juDFT_error("ntriad"
     197           0 :      +                    ,calledby ="triang")
     198             : c
     199           3 :                      IF ( k2.LE.k1 ) THEN
     200           0 :                         kk = k1
     201           0 :                         k1 = k2
     202           0 :                         k2 = kk
     203             :                      ENDIF
     204           3 :                      IF ( k3.LE.k1 ) THEN
     205           0 :                         kk = k1
     206           0 :                         k1 = k3
     207           0 :                         k3 = kk
     208             :                      ENDIF
     209           3 :                      IF ( k3.LE.k2 ) THEN
     210           0 :                         kk = k2
     211           0 :                         k2 = k3
     212           0 :                         k3 = kk
     213             :                      ENDIF
     214           3 :                      It(1,Ntria) = k1
     215           3 :                      It(2,Ntria) = k2
     216           3 :                      It(3,Ntria) = k3
     217           3 :                      At(Ntria) = ABS(a1)/2
     218           3 :                      Att = Att + At(Ntria)
     219             :                   ENDIF
     220             :                ENDIF
     221             :             ENDDO
     222             :          ENDDO
     223             :       ELSE
     224             :       ENDIF
     225             : c
     226             : 99001 FORMAT (' $$$ error in triang: collinear k-points'/(5x,2F12.6))
     227             : c
     228             :       END SUBROUTINE triang
     229             : !-----------------------------------------------------
     230           0 :       REAL FUNCTION vd2(x1,y1,x2,y2)         ! distance between (x1,y1) and (x2,y2)
     231             :          REAL x1,y1,x2,y2
     232         185 :          vd2 = (x2-x1)**2 + (y2-y1)**2
     233           0 :       END FUNCTION vd2
     234           0 :       REAL FUNCTION area(x1,y1,x2,y2,x3,y3)  ! area of triangle (x1,y1),(x2,y2),(x3,y3)
     235             :          REAL x1,y1,x2,y2,x3,y3
     236          85 :          area = (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)
     237           0 :       END FUNCTION area
     238             : !-----------------------------------------------------
     239             :       END MODULE m_triang

Generated by: LCOV version 1.13