LCOV - code coverage report
Current view: top level - inpgen - bravais_symm.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 72 76 94.7 %
Date: 2019-09-08 04:53:50 Functions: 3 3 100.0 %

          Line data    Source code
       1             :       MODULE m_bravaissymm
       2             :       use m_juDFT
       3             : !********************************************************************
       4             : !     determines the point group of the bravais lattice given the
       5             : !     lattice vectors. the idea is to determine all the lattice
       6             : !     vectors that have the same length as a_{1,2,3}, and then use
       7             : !     these to determine the possible rotation matrices.
       8             : !     these rotation matrices are in lattice coordinates.    mw 12-99
       9             : !********************************************************************
      10             :       CONTAINS
      11           3 :       SUBROUTINE bravais_symm(
      12             :      >                        as,bs,scale,nop48,neig12,
      13           3 :      <                        nops,mrot)
      14             : 
      15             :       IMPLICIT NONE
      16             : 
      17             : !==> Arguments
      18             :       INTEGER, INTENT (IN) :: nop48  ! max. number of operations (typically 48)
      19             :       INTEGER, INTENT (IN) :: neig12 ! max. number of lattice vectors with same length
      20             :                                      ! (max occurs for close-packed fcc: 12)
      21             :       REAL,    INTENT (IN) :: as(3,3),bs(3,3),scale(3) ! lattice information
      22             :       INTEGER, INTENT(OUT) :: nops, mrot(3,3,nop48)    ! point group operations
      23             : 
      24             : !==> Locals
      25             :       REAL    amet(3,3),b1,b2,b3,d1,d2,d3,dmax,dt
      26             :       INTEGER i,k,k1,k2,k3,m1,m2,m3,n1,n2,n3
      27             :       INTEGER irot(3,3)
      28           6 :       INTEGER lv1(3,neig12),lv2(3,neig12),lv3(3,neig12)
      29             : 
      30             :       REAL, PARAMETER :: eps=1.0e-9
      31             : 
      32             : !---> set up metric for distances
      33           3 :       amet = 0.0
      34          12 :       DO i = 1,3
      35           9 :          amet(1,1) = amet(1,1) + (scale(i)**2)*as(i,1)*as(i,1)
      36           9 :          amet(2,2) = amet(2,2) + (scale(i)**2)*as(i,2)*as(i,2)
      37           9 :          amet(3,3) = amet(3,3) + (scale(i)**2)*as(i,3)*as(i,3)
      38           9 :          amet(2,1) = amet(2,1) + (scale(i)**2)*as(i,1)*as(i,2)
      39           9 :          amet(3,2) = amet(3,2) + (scale(i)**2)*as(i,2)*as(i,3)
      40          12 :          amet(3,1) = amet(3,1) + (scale(i)**2)*as(i,3)*as(i,1)
      41             :       ENDDO
      42             : 
      43           3 :       amet(1,2) = amet(2,1)
      44           3 :       amet(2,3) = amet(3,2)
      45           3 :       amet(1,3) = amet(3,1)
      46             : 
      47             : !---> distances for the lattice vectors
      48           3 :       d1 = amet(1,1)
      49           3 :       d2 = amet(2,2)
      50           3 :       d3 = amet(3,3)
      51             :       b1 = ( bs(1,1)/scale(1) )**2 + ( bs(1,2)/scale(2) )**2             &
      52           3 :      &   + ( bs(1,3)/scale(3) )**2
      53             :       b2 = ( bs(2,1)/scale(1) )**2 + ( bs(2,2)/scale(2) )**2             &
      54           3 :      &   + ( bs(2,3)/scale(3) )**2
      55             :       b3 = ( bs(3,1)/scale(1) )**2 + ( bs(3,2)/scale(2) )**2             &
      56           3 :      &   + ( bs(3,3)/scale(3) )**2
      57             : 
      58             : !---> determine the cutoffs along each direction a_i:
      59           3 :       dmax = max( d1,d2,d3)
      60             : 
      61           3 :       m1 = nint( dmax * b1 )
      62           3 :       m2 = nint( dmax * b2 )
      63           3 :       m3 = nint( dmax * b3 )
      64             : 
      65             : !---->loop over all possible lattice vectors to find those with the
      66             : !---->length, i.e., ones that could be rotations
      67           3 :       n1 = 1
      68           3 :       n2 = 1
      69           3 :       n3 = 1
      70             : 
      71          12 :       lv1(1:3,1) = (/ 1,0,0 /)
      72          12 :       lv2(1:3,1) = (/ 0,1,0 /)
      73          12 :       lv3(1:3,1) = (/ 0,0,1 /)
      74             : 
      75          18 :       DO k3=-m3,m3
      76          93 :          DO k2=-m2,m2
      77         465 :             DO k1=-m1,m1
      78             : 
      79         375 :                dt = distance2(k1,k2,k3)
      80             : 
      81             : !---->    check if the same length
      82         375 :                IF ( abs( dt - d1 ) < eps ) THEN
      83          36 :                   IF (.not.( k1==1 .and. k2==0 .and. k3==0 ) ) THEN
      84          33 :                      n1 = n1+1
      85          33 :                      IF(n1>neig12)  CALL juDFT_error("n1>neig12",
      86           0 :      +                   calledby ="bravais_symm")
      87          33 :                      lv1(1,n1) = k1
      88          33 :                      lv1(2,n1) = k2
      89          33 :                      lv1(3,n1) = k3
      90             :                   ENDIF
      91             :                ENDIF
      92             : 
      93         375 :                IF ( abs( dt - d2 ) < eps ) THEN
      94          36 :                   IF (.not.( k1==0 .and. k2==1 .and. k3==0 ) ) THEN
      95          33 :                      n2 = n2+1
      96          33 :                      IF(n2>neig12)  CALL juDFT_error("n2>neig12",
      97           0 :      +                    calledby="bravais_symm")
      98          33 :                      lv2(1,n2) = k1
      99          33 :                      lv2(2,n2) = k2
     100          33 :                      lv2(3,n2) = k3
     101             :                   ENDIF
     102             :                ENDIF
     103             : 
     104         450 :                IF ( abs( dt - d3 ) < eps ) THEN
     105          36 :                   IF (.not.( k1==0 .and. k2==0 .and. k3==1 ) ) THEN
     106          33 :                      n3 = n3+1
     107          33 :                      IF(n3>neig12)  CALL juDFT_error("n3>neig12",
     108           0 :      +                    calledby="bravais_symm")
     109          33 :                      lv3(1,n3) = k1
     110          33 :                      lv3(2,n3) = k2
     111          33 :                      lv3(3,n3) = k3
     112             :                   ENDIF
     113             :                ENDIF
     114             : 
     115             :             ENDDO
     116             :          ENDDO
     117             :       ENDDO
     118             : 
     119             : !---> the possible rotation matrices are given by the matrix of
     120             : !---> column vectors of lv_{1,2,3}
     121           3 :       nops = 0
     122          39 :       DO k3 = 1,n3
     123         471 :          DO k2 = 1,n2
     124        5652 :             DO k1 = 1,n1
     125             : 
     126             : !--->          check whether determinant is +/-1 (needs to be for rotation)
     127        5184 :                IF ( abs(mdet(k1,k2,k3)) .NE. 1 ) CYCLE
     128             : 
     129             : !--->          check whether this maintains lengths correctly
     130             : !--->          if M is the metric, then must have R^T M R = M 
     131             :                irot = reshape( (/ lv1(:,k1),lv2(:,k2),lv3(:,k3) /) ,
     132        2304 :      &                         (/ 3,3 /) )
     133        7488 :                IF ( any( abs(
     134             :      &           matmul( transpose(irot), matmul(amet,irot) ) - amet 
     135        4608 :      &                     ) > eps ) ) CYCLE
     136             : 
     137         144 :                nops = nops + 1
     138         144 :                IF ( nops > nop48 )  CALL juDFT_error("nop > nop48",
     139           0 :      &                                     calledby="bravais_symm")
     140         576 :                mrot(:,:,nops) = irot
     141             : 
     142             :             ENDDO
     143             :          ENDDO
     144             :       ENDDO
     145             : 
     146             :       WRITE (6,'(//," Point group of the Bravais lattice has ",i2,
     147           3 :      &        " operations")') nops
     148             : 
     149           3 :       RETURN
     150             : 
     151             :       CONTAINS   ! INTERNAL routines
     152             : 
     153         375 :       REAL FUNCTION distance2(l1,l2,l3)
     154             : !*********************************************************************
     155             : !     calculates the magnitude square for a vector (l1,l2,l3) given in
     156             : !     lattice units
     157             : !*********************************************************************
     158             :       IMPLICIT NONE
     159             : 
     160             :       INTEGER, INTENT(IN) :: l1,l2,l3
     161             : 
     162             :       distance2 = l1*(l1*amet(1,1) + 2*l2*amet(2,1))
     163             :      &          + l2*(l2*amet(2,2) + 2*l3*amet(3,2))
     164         375 :      &          + l3*(l3*amet(3,3) + 2*l1*amet(1,3))
     165             : 
     166             :       RETURN
     167             :       END FUNCTION distance2
     168             : 
     169        5184 :       INTEGER FUNCTION mdet(k1,k2,k3)
     170             : !*********************************************************************
     171             : !     determines the determinant for possible rotation matrix
     172             : !     ( lv1(:,k1) ; lv2(:,k2) ; lv3(:,k3) )
     173             : !*********************************************************************
     174             :       IMPLICIT NONE
     175             : 
     176             :       INTEGER, INTENT(IN) :: k1,k2,k3
     177             : 
     178             :       mdet = lv1(1,k1)*( lv2(2,k2)*lv3(3,k3) - lv2(3,k2)*lv3(2,k3) )
     179             :      &     + lv1(2,k1)*( lv2(3,k2)*lv3(1,k3) - lv2(1,k2)*lv3(3,k3) )
     180        5184 :      &     + lv1(3,k1)*( lv2(1,k2)*lv3(2,k3) - lv2(2,k2)*lv3(1,k3) )
     181             : 
     182             :       RETURN
     183             :       END FUNCTION mdet
     184             : 
     185             :       END SUBROUTINE bravais_symm
     186             :       END MODULE m_bravaissymm

Generated by: LCOV version 1.13