setab.f 3 KB
Newer Older
1 2 3 4 5 6 7
      MODULE m_setab
      use m_juDFT
!*********************************************************************
!     set up lattice quantities and matrices 
!*********************************************************************
      CONTAINS
      SUBROUTINE setab(
8
     >                 a1,a2,a3,aa,scale,
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
     <                 amat,bmat,aamat,bbmat,amatinv,omtil)

      USE m_constants

      IMPLICIT NONE

!==>  Arguments
      REAL, INTENT (IN)  :: aa
      REAL, INTENT (IN)  :: a1(3),a2(3),a3(3),scale(3)
      REAL, INTENT (OUT) :: amat(3,3),bmat(3,3),amatinv(3,3)
      REAL, INTENT (OUT) :: aamat(3,3),bbmat(3,3)
      REAL, INTENT (OUT) :: omtil

!==>  Locals
      INTEGER i, j
      REAL    volume
      LOGICAL lerr
      REAL    tmat(3,3),b1(3),b2(3),b3(3)

!  volume in scaled Cartesian units
      volume  = a1(1)*a2(2)*a3(3) + a2(1)*a3(2)*a1(3) +
     &          a3(1)*a1(2)*a2(3) - a1(3)*a2(2)*a3(1) -
     &          a2(3)*a3(2)*a1(1) - a3(3)*a1(2)*a2(1)

!  reciprocal lattice vectors in scaled Cartesian units
      b1(1) = (a2(2)*a3(3) - a2(3)*a3(2))/volume
      b1(2) = (a2(3)*a3(1) - a2(1)*a3(3))/volume
      b1(3) = (a2(1)*a3(2) - a2(2)*a3(1))/volume
      b2(1) = (a3(2)*a1(3) - a3(3)*a1(2))/volume
      b2(2) = (a3(3)*a1(1) - a3(1)*a1(3))/volume
      b2(3) = (a3(1)*a1(2) - a3(2)*a1(1))/volume
      b3(1) = (a1(2)*a2(3) - a1(3)*a2(2))/volume
      b3(2) = (a1(3)*a2(1) - a1(1)*a2(3))/volume
      b3(3) = (a1(1)*a2(2) - a1(2)*a2(1))/volume

!  volume and area (assuming a1 and a2 define surface periodicity)
      omtil = (aa**3)*scale(1)*scale(2)*scale(3)*volume

!  matrices of lattice vectors in full Cartesian units
48 49

      DO i=1,3
50 51 52
         amat(i,1) = aa*scale(1)*a1(i)
         amat(i,2) = aa*scale(2)*a2(i)
         amat(i,3) = aa*scale(3)*a3(i)
53 54 55
      ENDDO

      DO i=1,3
56 57 58
         bmat(1,i) = (pi_const/(aa*scale(1))) * b1(i)
         bmat(2,i) = (pi_const/(aa*scale(2))) * b2(i)
         bmat(3,i) = (pi_const/(aa*scale(3))) * b3(i)
59 60 61
      ENDDO

      DO i=1,3
62 63 64
         amatinv(1,i) = (1.0/(aa*scale(1))) * b1(i)
         amatinv(2,i) = (1.0/(aa*scale(2))) * b2(i)
         amatinv(3,i) = (1.0/(aa*scale(3))) * b3(i)
65 66
      ENDDO

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
!--->  check that amat and amatinv consistent 
!      (amat*amatinv should be identity)

      tmat = matmul( amat, amatinv )
      lerr = .false.
      DO j=1,3
         if( abs( tmat(j,j) - 1.000 ) .gt. 1.e-10 ) lerr = .true.
         DO i=1,3
            if(i.eq.j) cycle
            if( abs( tmat(i,j) ) .gt. 1.e-10 ) lerr = .true.
         ENDDO
      ENDDO
      IF (lerr) THEN
         WRITE(6,'(" error in set-up of amat and amatinv matrices")')
         WRITE(6,'(" (",3f12.8," )")') tmat(1,1),tmat(1,2),tmat(1,3)
         WRITE(6,'(" (",3f12.8," )")') tmat(2,1),tmat(2,2),tmat(2,3)
         WRITE(6,'(" (",3f12.8," )")') tmat(3,1),tmat(3,2),tmat(3,3)
         CALL juDFT_error("ERROR in amat,amatinv matrices",calledby
     +        ="setab")
      ENDIF

      aamat=matmul(transpose(amat),amat)
89
      bbmat=matmul(bmat,transpose(bmat))
90

91 92
      END SUBROUTINE setab
      END MODULE m_setab