closure.f90 6.11 KB
Newer Older
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
MODULE m_closure

use m_juDFT

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Contains 3 subroutines that more or less check the closure:
!     closure :    checks whether the space group operations close
!     close_pt:    checks that the point group of the bravais 
!                  lattice closes
!     check_close: additionally calculate the multiplication table,
!                  inverse operations and also determines the type 
!                  of every operation                    mw99,gs00
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

CONTAINS

SUBROUTINE closure(mops,mrot,tau,nops,index_op,lclose)

   IMPLICIT NONE

   INTEGER, INTENT (IN)  :: mops           ! number of operations of the bravais lattice
   INTEGER, INTENT (IN)  :: nops           ! number of operations in space group
   INTEGER, INTENT (IN)  :: mrot(3,3,mops) ! refer to the operations of the 
   REAL,    INTENT (IN)  :: tau(3,mops)    ! bravais lattice
   INTEGER, INTENT (IN)  :: index_op(nops) ! mapping function between space group 
                                              ! op's and those of the bravais lattice
   LOGICAL, INTENT (OUT) :: lclose

   REAL    ttau(3),eps7
   INTEGER i,ii,j,jj,k,kk,mp(3,3),map(nops)

   eps7 = 1.0e-7

   ! loop over all operations
   DO jj = 1, nops
      j = index_op(jj)

      map(1:nops) = 0

      ! multiply {R_j|t_j}{R_i|t_i}
      DO ii = 1, nops
         i = index_op(ii)
         mp = matmul( mrot(:,:,j) , mrot(:,:,i) )
         ttau = tau(:,j) + matmul( mrot(:,:,j) , tau(:,i) )
         ttau = ttau - anint( ttau - eps7 )

         ! determine which operation this is
         DO kk=1,nops
            k = index_op(kk)
            IF ( all( mp(:,:) == mrot(:,:,k) ) .AND. all( abs( ttau(:)-tau(:,k) ) < eps7 ) ) THEN
               IF ( map(ii) .eq. 0 ) THEN
                  map(ii) = kk
               ELSE
                  write(6,*)'ERROR Closure: Multiplying ', jj,' with ',kk, ' and with ',map(ii)
                  write(6,*) 'yields the same matrix'
                  lclose = .false.
                  RETURN
               END IF
            END IF
         END DO

         IF (map(ii).eq.0) THEN
            write(6,*)'ERROR Closure:',ii,' times',jj,' leaves group'
            lclose = .false.
            RETURN
         END IF
      END DO
   END DO

   lclose = .true.

END SUBROUTINE closure

!*********************************************************************

SUBROUTINE close_pt(nops,mrot,mtable)

   IMPLICIT NONE

   INTEGER, INTENT (IN)  :: nops,mrot(3,3,nops)
   INTEGER, INTENT (OUT) :: mtable(nops,nops)   ! table(i,j) = {R_i|0}{R_j|0}

   INTEGER              :: i,j,k,mp(3,3),map(nops)

   ! loop over all operations
   DO j = 1, nops

      map(1:nops) = 0

      ! multiply {R_j|0}{R_i|0}
      DO i = 1, nops
         mp = matmul( mrot(:,:,j) , mrot(:,:,i) )

         ! determine which operation this is
         DO k = 1, nops
            IF ( all( mp(:,:)==mrot(:,:,k) ) ) THEN
               IF ( map(i) .eq. 0 ) THEN
                  map(i) = k
               ELSE
                  WRITE (6,'(" Symmetry error : multiple ops")')
                  CALL juDFT_error("close_pt: Multiple ops (Bravais)",calledby ="closure")
               END IF
            END IF
         END DO

         IF (map(i).eq.0) THEN
            WRITE(6,*) 'Symmetry operations:'
            DO k = 1, nops
               WRITE(6,*) 'Matrix ', k, ':'
               WRITE(6,'(3i7)') mrot(:,1,k)
               WRITE(6,'(3i7)') mrot(:,2,k)
               WRITE(6,'(3i7)') mrot(:,3,k)
               WRITE(6,*) ''
            END DO
            WRITE (6,'(" Group not closed (Bravais lattice)")')
            WRITE (6,'(" operation j=",i2,"  map=",12i4,:/,(21x,12i4))')  j, map(1:nops)
            WRITE(6,*) ''
            WRITE(6,*) 'Expected product of operations ', j, ' and ', i, ':'
            WRITE(6,'(3i7)') mp(:,1)
            WRITE(6,'(3i7)') mp(:,2)
            WRITE(6,'(3i7)') mp(:,3)
            WRITE(6,*) ''
            CALL juDFT_error("close_pt:Not closed",calledby="closure")
         END IF
      END DO
      mtable(j,1:nops) = map(1:nops)
   END DO

END SUBROUTINE close_pt

!*********************************************************************

SUBROUTINE check_close(nops,mrot,tau,multtab,inv_op,optype)

   IMPLICIT NONE

   INTEGER, INTENT (IN)  :: nops
   INTEGER, INTENT (IN)  :: mrot(3,3,nops)
   REAL,    INTENT (IN)  :: tau(3,nops)
   INTEGER, INTENT (OUT) :: inv_op(nops)
   INTEGER, INTENT (OUT) :: multtab(nops,nops)
   INTEGER, INTENT (OUT) :: optype(nops)

   REAL    ttau(3)
   INTEGER i,j,n,k,mp(3,3),mdet,mtr

   REAL,    PARAMETER :: eps=1.0e-7
   INTEGER, PARAMETER :: cops(-1:3)=(/ 2, 3, 4, 6, 1 /)

   inv_op(1:nops) = 0

   multtab = 0

   ! loop over all operations
   DO j = 1, nops

      ! multiply {R_j|t_j}{R_i|t_i}
      DO i = 1, nops
         mp = matmul( mrot(:,:,j) , mrot(:,:,i) )
         ttau = tau(:,j) + matmul( mrot(:,:,j) , tau(:,i) )
         ttau = ttau - anint( ttau - eps )

         ! determine which operation this is
         DO k=1,nops
            IF ( all( mp(:,:) == mrot(:,:,k) ) .and. all( abs( ttau(:)-tau(:,k) ) < eps ) ) THEN
               IF ( multtab(j,i) .eq. 0 ) THEN
                  multtab(j,i) = k
                  IF (k .eq. 1) inv_op(j)=i
               ELSE
                  WRITE(6,'(" Symmetry error: multiple ops")')
                  CALL juDFT_error("check_close: Multiple ops",calledby ="closure")
               END IF
            END IF
         END DO

         IF (multtab(j,i).eq.0) THEN
            WRITE (6,'(" Group not closed")')
            WRITE (6,'("  j , i =",2i4)') j,i
            CALL juDFT_error("check_close: Not closed",calledby="closure")
         END IF
      END DO
   END DO

   ! determine the type of each operation
   DO n = 1, nops
      mtr = mrot(1,1,n) + mrot(2,2,n) + mrot(3,3,n)
      mdet = mrot(1,1,n)*(mrot(2,2,n)*mrot(3,3,n)-mrot(3,2,n)*mrot(2,3,n)) +&
             mrot(1,2,n)*(mrot(3,1,n)*mrot(2,3,n)-mrot(2,1,n)*mrot(3,3,n)) +&
             mrot(1,3,n)*(mrot(2,1,n)*mrot(3,2,n)-mrot(3,1,n)*mrot(2,2,n))

      optype(n) = mdet*cops(mdet*mtr)

   END DO

END SUBROUTINE check_close

END MODULE m_closure