SUBROUTINE symop (irrep, j, i, gamma, n, isymop, n1, kase, fac)
USE FileUnits_Module
USE SymGroup_Module
IMPLICIT NONE
!==============================================================================
INTEGER :: n, irrep, j, i, gamma, gammaC6v (12, 8), isymop, &
 gammaC1 (1, 1), gammaC2 (2, 2), gammaC1v (2, 2), gammaC2v (4, 4), &
 gammaC3v (6, 4), gammaC3 (3, 3), gammaC6 (6, 6), n1 (8), &
 kase, fac
DATA gammaC1 / 1 /
DATA gammaC2 / 1, 1, 1, - 1 /
DATA gammaC1v / 1, 1, 1, - 1 /
DATA gammaC2v / 1, 1, 1, 1, 1, 1, - 1, - 1, 1, - 1, 1, - 1, 1, &
 - 1, - 1, 1 /
DATA gammaC3 / 1, 1, 1, 2, - 1, - 1, 2, - 1, - 1 /
DATA gammaC3v / 1, 1, 1, 1, 1, 1, 1, 1, 1, - 1, - 1, - 1, 2, &
 - 1, - 1, 0, 0, 0, 0, 1, - 1, 0, 0, 0 /
DATA gammaC6 / 1, 1, 1, 1, 1, 1, 1, - 1, 1, - 1, 1, - 1, 2, 1, &
 - 1, - 2, - 1, 1, 0, 1, 1, 0, - 1, - 1, 2, - 1, - 1, 2, - 1, &
 - 1, 0, 1, - 1, 0, 1, - 1 /

DATA gammaC6v / 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
 1, - 1, - 1, - 1, - 1, - 1, - 1, 1, - 1, - 1, 1, 1, - 1, 1, 1, 1, &
 - 1, - 1, - 1, 1, - 1, - 1, 1, 1, - 1, - 1, - 1, - 1, 1, 1, 1, 2, &
 1, 1, - 1, - 1, - 2, 0, 0, 0, 0, 0, 0, 0, 1, - 1, 1, - 1, 0, 0, 0, &
 0, 0, 0, 0, 2, - 1, - 1, - 1, - 1, 2, 0, 0, 0, 0, 0, 0, 0, 1, &
 - 1, - 1, 1, 0, 0, 0, 0, 0, 0, 0 /
!==============================================================================
fac = 0
IF (group.eq.'C1 ') THEN
   IF (n.lt.0) THEN
      WRITE(Output_Unit, * ) 'Error: n=', n
      STOP 'symop'
   ENDIF
   IF (kase.eq.0) THEN
WRITE(Output_Unit,  * ) 'Warning NOT exploiting full symmetry for ABC'
      nirrep = 1
      nsymop = 1
      n1 (1) = n
      par (1) = 0
      IF (system.eq.'ABC'.or.system.eq.'ABB'.or.system.eq.'AAA') THEN
         WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for', ' system=', system
      ENDIF
      RETURN
   ENDIF
   IF (irrep.lt.1.or.irrep.gt.1) THEN
      WRITE(Output_Unit, * ) 'irrep=', irrep
      STOP 'symop'
   ENDIF
   IF (isymop.lt.1.or.isymop.gt.1) THEN
      WRITE(Output_Unit, * ) 'isymop=', isymop
      STOP 'symop'
   ENDIF
   gamma = gammaC1 (isymop, irrep)
   IF (isymop.eq.1) THEN
      i = j
   ELSE
      WRITE(Output_Unit, * ) 'Error: isymop=', isymop
   ENDIF
ELSEIF (group.eq.'C2 ') THEN
   IF (n.lt.0.or. (n / 2) * 2.ne.n) THEN
      WRITE(Output_Unit, * ) 'Error: n=', n
      STOP 'symop'
   ENDIF
   IF (kase.eq.0) THEN
      WRITE(Output_Unit, * ) 'Exploiting full symmetry for ABC'
      nirrep = 2
      nsymop = 2
      n1 (1) = n / 2
      n1 (2) = n / 2
      par (1) = 0
      par (2) = 1
      IF (system.eq.'ABC'.or.system.eq.'ABB'.or.system.eq.'AAA') THEN
         WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for',  ' system=', system
      ENDIF
      RETURN
   ENDIF
   IF (irrep.lt.1.or.irrep.gt.2) THEN
      WRITE(Output_Unit, * ) 'irrep=', irrep
      STOP 'symop'
   ENDIF
   IF (isymop.lt.1.or.isymop.gt.2) THEN
      WRITE(Output_Unit, * ) 'isymop=', isymop
      STOP 'symop'
   ENDIF
   gamma = gammaC2 (isymop, irrep)
   IF (isymop.eq.1) THEN
      i = j
   ELSEIF (isymop.eq.2) THEN
      i = mod (j - 1 + 6 * n / 12, n) + 1
   ELSE
      WRITE(Output_Unit, * ) 'Error: isymop=', isymop
   ENDIF
ELSEIF (group.eq.'C1v') THEN
   IF (system.eq.'ABC') THEN
      WRITE(Output_Unit, * ) 'system=', system
      STOP 'symop'
   ENDIF
   IF (n.lt.0.or. (n / 2) * 2.ne.n) THEN
      WRITE(Output_Unit, * ) 'Error: n=', n
      STOP 'symop'
   ENDIF
   IF (kase.eq.0) THEN
WRITE(Output_Unit,  * ) 'Warning NOT exploiting full symmetry for AB2'
      nirrep = 2
      nsymop = 2
      n1 (1) = n / 2 + 1
      n1 (2) = n / 2 - 1
      par (1) = 0
      par (1) = 0
      IF (system.eq.'ABB'.or.system.eq.'AAA') THEN
         WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for',' system=', system
      ENDIF
      RETURN
   ENDIF
   IF (irrep.lt.1.or.irrep.gt.2) THEN
      WRITE(Output_Unit, * ) 'irrep=', irrep
      STOP 'symop'
   ENDIF
   IF (isymop.lt.1.or.isymop.gt.2) THEN
      WRITE(Output_Unit, * ) 'isymop=', isymop
      STOP 'symop'
   ENDIF
   gamma = gammaC1v (isymop, irrep)
   IF (isymop.eq.1) THEN
      i = j
   ELSEIF (isymop.eq.2) THEN
      i = mod (n + 1 - j, n) + 1
   ELSE
      WRITE(Output_Unit, * ) 'Error: isymop=', isymop
   ENDIF
ELSEIF (group.eq.'C2v') THEN
   IF (system.eq.'ABC') THEN
      WRITE(Output_Unit, * ) 'system=', system
      STOP 'symop'
   ENDIF
   IF (n.lt.0.or. (n / 4) * 4.ne.n) THEN
      WRITE(Output_Unit, * ) 'Error: n=', n
      STOP 'symop'
   ENDIF
   IF (kase.eq.0) THEN
      IF (system.eq.'AB2') THEN
         WRITE(Output_Unit, * ) 'Exploiting full symmetry for AB2'
      ENDIF
      nirrep = 4
      nsymop = 4
      n1 (1) = n / 4 + 1
      n1 (2) = n / 4 - 1
      n1 (3) = n / 4
      n1 (4) = n / 4
      par (1) = 0
      par (2) = 1
      par (3) = 0
      par (4) = 1
      chishift_ir(1)=0
      chishift_ir(2)=1
      chishift_ir(3)=0
      chishift_ir(4)=1
      IF (system.eq.'AAA') THEN
         WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for',  ' system=', system
      ENDIF
      RETURN
   ENDIF
   IF (irrep.lt.1.or.irrep.gt.4) THEN
      WRITE(Output_Unit, * ) 'irrep=', irrep
      STOP 'symop'
   ENDIF
   IF (isymop.lt.1.or.isymop.gt.4) THEN
      WRITE(Output_Unit, * ) 'isymop=', isymop
      STOP 'symop'
   ENDIF
   gamma = gammaC2v (isymop, irrep)
   IF (isymop.eq.1) THEN
      i = j
   ELSEIF (isymop.eq.2) THEN
      i = mod (j - 1 + 6 * n / 12, n) + 1
   ELSEIF (isymop.eq.3) THEN
      i = mod (n + 1 - j, n) + 1
   ELSEIF (isymop.eq.4) THEN
      i = mod (n + 1 - j + 6 * n / 12, n) + 1
   ELSE
      WRITE(Output_Unit, * ) 'Error: isymop=', isymop
   ENDIF
ELSEIF (group.eq.'C3 ') THEN
   IF (system.eq.'ABC'.or.system.eq.'ABB') THEN
      WRITE(Output_Unit, * ) 'Error: system=', system
      STOP 'symop'
   ENDIF
   IF (n.lt.0.or. (n / 6) * 6.ne.n) THEN
      WRITE(Output_Unit, * ) 'Error: n=', n
      STOP 'symop'
   ENDIF
   IF (kase.eq.0) THEN
      WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for A3'
      nirrep = 3
      nsymop = 3
      n1 (1) = n / 3
      n1 (2) = n / 3
      n1 (3) = n / 3
      par (1) = 0
      par (1) = 0
      par (1) = 0
      IF (system.eq.'AAA') THEN
         WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for', ' system=', system
      ENDIF
      RETURN
   ENDIF
   IF (irrep.lt.1.or.irrep.gt.3) THEN
      WRITE(Output_Unit, * ) 'irrep=', irrep
      STOP 'symop'
   ENDIF
   IF (isymop.lt.1.or.isymop.gt.3) THEN
      WRITE(Output_Unit, * ) 'isymop=', isymop
      STOP 'symop'
   ENDIF
   IF (irrep.eq.1) THEN
      fac = 0
   ELSEIF (irrep.eq.2) THEN
      fac = 1
   ELSE
      fac = - 1
   ENDIF
   gamma = gammaC3 (isymop, irrep)
   IF (isymop.eq.1) THEN
      i = j
   ELSEIF (isymop.eq.2) THEN
      i = mod (j - 1 + 4 * n / 12, n) + 1
   ELSEIF (isymop.eq.3) THEN
      i = mod (j - 1 + 8 * n / 12, n) + 1
   ELSE
      WRITE(Output_Unit, * ) 'Error: isymop=', isymop
   ENDIF
ELSEIF (group.eq.'C3v') THEN
   IF (system.eq.'ABC'.or.system.eq.'ABB') THEN
      WRITE(Output_Unit, * ) 'Error: system=', system
      STOP 'symop'
   ENDIF
   IF (n.lt.0.or. (n / 6) * 6.ne.n) THEN
      WRITE(Output_Unit, * ) 'Error: n=', n
      STOP 'symop'
   ENDIF
   IF (kase.eq.0) THEN
      WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for A3'
      nirrep = 4
      nsymop = 6
      n1 (1) = n / 6 + 1
      n1 (2) = n / 6 - 1
      n1 (3) = 2 * n / 6
      n1 (4) = 2 * n / 6
      par (1) = 0
      par (2) = 0
      par (3) = 0
      par (4) = 0
      IF (system.eq.'AAA') THEN
         WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for', ' system=', system
      ENDIF
      RETURN
   ENDIF
   IF (irrep.lt.1.or.irrep.gt.4) THEN
      WRITE(Output_Unit, * ) 'irrep=', irrep
      STOP 'symop'
   ENDIF
   IF (isymop.lt.1.or.isymop.gt.6) THEN
      WRITE(Output_Unit, * ) 'isymop=', isymop
      STOP 'symop'
   ENDIF
   IF (irrep.le.2) THEN
      fac = 0
   ELSEIF (irrep.eq.3) THEN
      fac = 1
   ELSE
      fac = - 1
   ENDIF
   gamma = gammaC3v (isymop, irrep)
   IF (isymop.eq.1) THEN
      i = j
   ELSEIF (isymop.eq.2) THEN
      i = mod (j - 1 + 4 * n / 12, n) + 1
   ELSEIF (isymop.eq.3) THEN
      i = mod (j - 1 + 8 * n / 12, n) + 1
   ELSEIF (isymop.eq.4) THEN
      i = mod (n + 1 - j, n) + 1
   ELSEIF (isymop.eq.5) THEN
      i = mod (n + 1 - j + 4 * n / 12, n) + 1
   ELSEIF (isymop.eq.6) THEN
      i = mod (n + 1 - j + 8 * n / 12, n) + 1
   ELSE
      WRITE(Output_Unit, * ) 'Error: isymop=', isymop
   ENDIF
ELSEIF (group.eq.'C6 ') THEN
   IF (system.eq.'ABC'.or.system.eq.'ABB') THEN
      WRITE(Output_Unit, * ) 'Error: system=', system
      STOP 'symop'
   ENDIF
   IF (n.lt.0.or. (n / 6) * 6.ne.n) THEN
      WRITE(Output_Unit, * ) 'Error: n=', n
      STOP 'symop'
   ENDIF
   IF (kase.eq.0) THEN
      WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for A3'
      nirrep = 6
      nsymop = 6
      n1 (1) = n / 6
      n1 (2) = n / 6
      n1 (3) = n / 6
      n1 (4) = n / 6
      n1 (5) = n / 6
      n1 (6) = n / 6
      par (1) = 0
      par (2) = 0
      par (3) = 0
      par (4) = 0
      par (5) = 0
      par (6) = 0
      IF (system.eq.'AAA') THEN
         WRITE(Output_Unit, * ) 'Warning NOT exploiting full symmetry for', ' system=', system
      ENDIF
      RETURN
   ENDIF
   IF (irrep.lt.1.or.irrep.gt.6) THEN
      WRITE(Output_Unit, * ) 'irrep=', irrep
      STOP 'symop'
   ENDIF
   IF (isymop.lt.1.or.isymop.gt.6) THEN
      WRITE(Output_Unit, * ) 'isymop=', isymop
      STOP 'symop'
   ENDIF
   IF (irrep.le.2) THEN
      fac = 0
   ELSEIF (irrep.eq.3.or.irrep.eq.5) THEN
      fac = 1
   ELSE
      fac = - 1
   ENDIF
   gamma = gammaC6 (isymop, irrep)
   IF (isymop.eq.1) THEN
      i = j
   ELSEIF (isymop.eq.2) THEN
      i = mod (j - 1 + 2 * n / 12, n) + 1
   ELSEIF (isymop.eq.3) THEN
      i = mod (j - 1 + 4 * n / 12, n) + 1
   ELSEIF (isymop.eq.4) THEN
      i = mod (j - 1 + 6 * n / 12, n) + 1
   ELSEIF (isymop.eq.5) THEN
      i = mod (j - 1 + 8 * n / 12, n) + 1
   ELSEIF (isymop.eq.6) THEN
      i = mod (j - 1 + 10 * n / 12, n) + 1
   ELSE
      WRITE(Output_Unit, * ) 'Error: isymop=', isymop
   ENDIF
ELSEIF (group.eq.'C6v') THEN
   IF (system.eq.'ABC'.or.system.eq.'ABB') THEN
      WRITE(Output_Unit, * ) 'Error: system=', system
      STOP 'symop'
   ENDIF
   IF (n.lt.0.or. (n / 12) * 12.ne.n) THEN
      WRITE(Output_Unit, * ) 'Error: n=', n
      STOP 'symop'
   ENDIF
   IF (kase.eq.0) THEN
      WRITE(Output_Unit, * ) 'Exploiting full symmetry for A3'
      nirrep = 8
      nsymop = 12
      n1 (1) = n / 12 + 1
      n1 (2) = n / 12 - 1
      n1 (3) = n / 12
      n1 (4) = n / 12
      n1 (5) = 2 * n / 12
      n1 (6) = 2 * n / 12
      n1 (7) = 2 * n / 12
      n1 (8) = 2 * n / 12
      par (1) = 0
      par (2) = 0
      par (3) = 1
      par (4) = 1
      par (5) = 1
      par (6) = 1
      par (7) = 0
      par (8) = 0
      ! Added by jeff
      chishift_ir(1)=0
      chishift_ir(2)=1
      chishift_ir(3)=0
      chishift_ir(4)=1
      chishift_ir(5)=0
      chishift_ir(6)=0
      chishift_ir(7)=0
      chishift_ir(8)=0
      RETURN
   ENDIF
   IF (irrep.lt.1.or.irrep.gt.8) THEN
      WRITE(Output_Unit, * ) 'irrep=', irrep
      STOP 'symop'
   ENDIF
   IF (isymop.lt.1.or.isymop.gt.12) THEN
      WRITE(Output_Unit, * ) 'isymop=', isymop
      STOP 'symop'
   ENDIF
   IF (irrep.le.4) THEN
      fac = 0
   ELSEIF (irrep.eq.5.or.irrep.eq.7) THEN
      fac = 1
   ELSE
      fac = - 1
   ENDIF
   gamma = gammaC6v (isymop, irrep)
   IF (isymop.eq.1) THEN
      i = j
   ELSEIF (isymop.eq.2) THEN
      i = mod (j - 1 + 2 * n / 12, n) + 1
   ELSEIF (isymop.eq.3) THEN
      i = mod (j - 1 + 10 * n / 12, n) + 1
   ELSEIF (isymop.eq.4) THEN
      i = mod (j - 1 + 4 * n / 12, n) + 1
   ELSEIF (isymop.eq.5) THEN
      i = mod (j - 1 + 8 * n / 12, n) + 1
   ELSEIF (isymop.eq.6) THEN
      i = mod (j - 1 + 6 * n / 12, n) + 1
   ELSEIF (isymop.eq.7) THEN
      i = mod (n + 1 - j, n) + 1
   ELSEIF (isymop.eq.8) THEN
      i = mod (n + 1 - j + 4 * n / 12, n) + 1
   ELSEIF (isymop.eq.9) THEN
      i = mod (n + 1 - j + 8 * n / 12, n) + 1
   ELSEIF (isymop.eq.10) THEN
      i = mod (n + 1 - j + 2 * n / 12, n) + 1
   ELSEIF (isymop.eq.11) THEN
      i = mod (n + 1 - j + 6 * n / 12, n) + 1
   ELSEIF (isymop.eq.12) THEN
      i = mod (n + 1 - j + 10 * n / 12, n) + 1
   ELSE
      WRITE(Output_Unit, * ) 'isymop=', isymop
      STOP 'symop'
   ENDIF
ELSE
   WRITE(Output_Unit, * ) 'Incorrect group: group=', group
   STOP 'symop'
ENDIF
IF (i.lt.1.or.i.gt.n) THEN
WRITE(Output_Unit,  * ) 'i=', i, '   n=', n
   STOP 'symop'
ENDIF
RETURN
END SUBROUTINE symop