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(i)*a1(i)
         amat(i,2) = aa*scale(i)*a2(i)
         amat(i,3) = aa*scale(i)*a3(i)
53
54
55
      ENDDO

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

      DO i=1,3
62
63
64
         amatinv(1,i) = (1.0/(aa*scale(i))) * b1(i)
         amatinv(2,i) = (1.0/(aa*scale(i))) * b2(i)
         amatinv(3,i) = (1.0/(aa*scale(i))) * 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