      SUBROUTINE ROSYM4 (A1,A2,A3,B1,B2,B3,ISY,NC,IB,V,R,RB,
     +  NPLWV,NG1,NG2,NG3,ROINP,ROSYMM,LPCTX,LPCTY,LPCTZ,
     +  LPCTXI,LPCTYI,LPCTZI,IPRINT)
      IMPLICIT COMPLEX (C)
C
C     WRITTEN ON JANUARY 2ND, 1980 - FROM ROSYM1.
C     LOOKUP OF G-VECTORS MODIFIED 19-MAY-82 BY OHN.
C     ROTATION MATRICES MODIFIED 02-MAR-84 BY OHN.
C
C SYMMETRIZATION OF (FOURIER COMPONENTS OF) CHARGE DENSITY
C BY POINT-GROUP SYMMETRY OPERATIONS
C SUBROUTINES NEEDED: LOOKUP
C ROSYMM(G)=(1/NC)*SUM( EXP(I*G*V(R(J)))*ROINP(R(J)**(-1)*G) )
C WHERE THE SUM RUNS OVER ALL NC POINT-GROUP ROTATIONS R.
C
C INPUT DATA:
C      A1....B3 ... ELEMENTARY TRANSLATIONS OF THE DIRECT
C               AND RECIPROCAL LATTICES
C      NPLWV ...NG1*NG2*NG3
C      NG1,...IGLIST ... SEE THE SUBROUTINE RLV4
C               INPUT DATA WHICH MAY BE OBTAINED FROM SBRT. GROUP1:
C      ISY .... CODE INDICATING WHETHER THE SPACE GROUP OF CRYSTAL
C               IS SYMMORPHIC OR NONSYMMORPHIC
C               ISY=0 MEANS NONSYMMORPHIC GROUP
C               ISY=1 MEANS SYMMORPHIC GROUP
C      NC ..... TOTAL NUMBER OF ELEMENTS IN THE POINT GROUP OF THE
C               CRYSTAL
C      IB ..... LIST OF THE ROTATIONS CONSTITUTING THE POINT GROUP
C               OF THE CRYSTAL. THE NUMBERING IS THAT DEFINED IN
C               WORLTON AND WARREN, I.E. THE ONE MATERIALIZED IN THE
C               ARRAY R (SEE BELOW)
C               ONLY THE FIRST NC ELEMENTS OF THE ARRAY IB ARE
C               MEANINGFUL
C      V ...... NONPRIMITIVE TRANSLATIONS (IN THE CASE OF NONSYMMOR-
C               PHIC GROUPS). V(I,N) IS THE I-TH COMPONENT
C               OF THE TRANSLATION CONNECTED WITH THE N-TH ELEMENT
C               OF THE POINT GROUP (I.E. WITH THE ROTATION
C               NUMBER IB(N) ).
C               ATTENTION: V(I) HAVE TO BE GIVEN IN THE
C               BASIS A1,A2,A3 AND NOT AS CARTESIAN COMPONENTS.
C               (I.E. ONE CAN USE THE OUTPUT OF THE
C               SUBROUTINE GROUP1 WITHOUT ANY MODIFICATION.)
C      R ...... LIST OF THE 3 X 3 ROTATION MATRICES
C               (XYZ REPRESENTATION OF THE O(H) OR D(6)H GROUPS)
C               ALL 48 OR 24 MATRICES ARE LISTED.
C      ROINP .. ARRAY WITH ELEMENTS RO(G) CORRESPONDING TO
C               G-VECTORS LISTED IN THE IGLIST.
C               THIS IS THE ARRAY TO BE SYMMETRIZED. (IN NORMAL
C               CASES, IT WILL BE OBTAINED AS WEIGHTED SUM OF
C               RO(WVK,G) (SBRT. RO5) OVER A SET OF SPECIAL
C               POINTS WVK.)
C      ROSYMM . WORK ARRAY
C
C OUTPUT DATA:
C      ROINP .. THE SYMMETRIZED ARRAY
C      ROSYMM . EQUAL TO ROINP ON EXIT
C
      DIMENSION A1(3),A2(3),A3(3),B1(3),B2(3),B3(3),IB(48),V(3,48)
      DIMENSION RB(48,3,3),SUM(3,3),IGL(3,48),R(49,3,3)
      DIMENSION LPCTX(NG1),LPCTY(NG2),LPCTZ(NG3)
      DIMENSION LPCTXI(-NG1/2:NG1/2),LPCTYI(-NG2/2:NG2/2)
      DIMENSION LPCTZI(-NG3/2:NG3/2)
      COMPLEX ROINP(NPLWV),ROSYMM(NPLWV),CROSUM
C
C.....PHYSICAL AND MATHEMATICAL CONSTANTS
      DOUBLE PRECISION PI
C
      DATA INITLZ /0/
      SAVE INITLZ
           PI   = 3.141592654D0
C-----------------------------------------------------------------------
C
      INIPRN = 0
      IF (INITLZ .NE. 0) GOTO 50
      INIPRN = 1
C
C     INITIALIZATION:
C
C     CHECK THAT THE FIRST SYMMETRY OPERATION IN THE LIST IS IDENTITY
      IF (IB(1) .NE. 1) THEN
        WRITE (*,*) '*** ERROR *** ROSYM4'
        WRITE (*,99 ) (IB(I),I=1,NC)
        CALL EXIT
      ENDIF
C
C     EXPRESS THE ROTATION MATRICES IN THE BASIS B1,B2,B3:
C     (ACTUALLY, RB CONTAINS R**(-1) IN THIS BASIS;
C     R**(-1) IS JUST R(TRANSPOSE), BUT NOT SO FOR RB !!!).
C
      DO 40 N = 1 , NC
        IC = IB(N)
        DO 10 I = 1 , 3
          DO 10 J = 1 , 3
            SUM(I,J) = 0.0
10      CONTINUE
        DO 20 I = 1 , 3
          DO 20 J = 1 , 3
            SUM(1,1) = SUM(1,1) + A1(I) * R(IC,J,I) * B1(J)
            SUM(2,1) = SUM(2,1) + A2(I) * R(IC,J,I) * B1(J)
            SUM(3,1) = SUM(3,1) + A3(I) * R(IC,J,I) * B1(J)
            SUM(1,2) = SUM(1,2) + A1(I) * R(IC,J,I) * B2(J)
            SUM(2,2) = SUM(2,2) + A2(I) * R(IC,J,I) * B2(J)
            SUM(3,2) = SUM(3,2) + A3(I) * R(IC,J,I) * B2(J)
            SUM(1,3) = SUM(1,3) + A1(I) * R(IC,J,I) * B3(J)
            SUM(2,3) = SUM(2,3) + A2(I) * R(IC,J,I) * B3(J)
            SUM(3,3) = SUM(3,3) + A3(I) * R(IC,J,I) * B3(J)
20      CONTINUE
        DO 30 I = 1 , 3
          DO 30 J = 1 , 3
            RB(N,I,J) = SUM(I,J)
30      CONTINUE
        IF (IPRINT.GE.3) WRITE (*,550) N,IC,
     &                  ((SUM(I,J),J=1,3),I=1,3)
40    CONTINUE
C
      DO 90 I = 1 , NG1
         IG1 = LPCTX(I)
         LPCTXI(IG1) = I
 90   CONTINUE
C
      DO 91 I = 1 , NG2
         IG2 = LPCTY(I)
         LPCTYI(IG2) = I
 91   CONTINUE
C
      DO 92 I = 1 , NG3
         IG3 = LPCTZ(I)
         LPCTZI(IG3) = I
 92   CONTINUE
C
      INITLZ = 1
C-----------------------------------------------------------------------
      IF (IPRINT.GE.1) THEN
        WRITE (*,*) '  A VECTORS IN ROSYM...'
        WRITE (*,1112) (A1(I),I=1,3)
        WRITE (*,1112) (A2(I),I=1,3)
        WRITE (*,1112) (A3(I),I=1,3)
        WRITE (*,*) '  B VECTORS IN ROSYM...'
        WRITE (*,1112) (B1(I),I=1,3)
        WRITE (*,1112) (B2(I),I=1,3)
        WRITE (*,1112) (B3(I),I=1,3)
      END IF
C
50    IF (NC .EQ. 1) THEN
        RETURN
      ELSE IF (NC .LT. 1 .OR. NC .GT. 48) THEN
        WRITE(*,*) '***ERROR*** ROSYM4'
        WRITE(*,*) 'NC = ',NC,' IS NOT ALLOWED'
        CALL EXIT
      ENDIF
C
C==================================================================
C  LOOP OVER G-VECTORS
C==================================================================
      INDEX = 0
C
      DO 100 K = 1 , NG3
        IG3 = LPCTZ(K)
        DO 100 J = 1 , NG2
          IG2 = LPCTY(J)
          DO 100 I = 1 , NG1
            IG1 = LPCTX(I)
            INDEX = INDEX + 1
            CROSUM = (0.0,0.0)
C
            DO 70 N = 2 , NC
              NN = N
C==================================================================
C   APPLICATION OF R**(-1) TO G (SEE NOTE ABOVE)
C==================================================================
              IGL1 = NINT( RB(N,1,1) * IG1 + RB(N,1,2) * IG2 + 
     &                     RB(N,1,3) * IG3 )
              IGL2 = NINT( RB(N,2,1) * IG1 + RB(N,2,2) * IG2 + 
     &                     RB(N,2,3) * IG3 )
              IGL3 = NINT( RB(N,3,1) * IG1 + RB(N,3,2) * IG2 + 
     &                     RB(N,3,3) * IG3 )
              IF (IGL1.GT.NG1/2 .OR. IGL1.LT.-NG1/2 .OR.
     &            IGL2.GT.NG2/2 .OR. IGL2.LT.-NG2/2 .OR.
     &            IGL3.GT.NG3/2 .OR. IGL3.LT.-NG3/2) GO TO 199
C
              IF (IGL1.EQ.-NG1/2) IGL1 = NG1 / 2
              IF (IGL2.EQ.-NG2/2) IGL2 = NG2 / 2
              IF (IGL3.EQ.-NG3/2) IGL3 = NG3 / 2
C
              IPLACE =  LPCTXI(IGL1) + NG1 * ( LPCTYI(IGL2) - 1 ) 
     +                  + NG1 * NG2 * ( LPCTZI(IGL3) - 1 )
              IF (ISY .EQ. 0) THEN
C==================================================================
C  NONSYMMORPHIC GROUP - SIGN OF ARG IS DIFFERENT FROM THAT IN
C  RICHARD NEEDS'S PROGRAM!
C==================================================================     
C        ARG = -2.0 * PI * ( IG1*V(1,N) + IG2*V(2,N) + IG3*V(3,N) )
C
                ARG = 2.0 * PI * ( IG1 * V(1,N) + IG2 * V(2,N) + 
     &                             IG3 * V(3,N) )
                CROSUM = CROSUM + 
     &                   ROINP(IPLACE) * CMPLX ( COS(ARG),SIN(ARG) )
              ELSE
C==================================================================
C  SYMMORPHIC GROUP
C==================================================================
                CROSUM = CROSUM + ROINP(IPLACE)
              ENDIF
C
 70         CONTINUE
C==================================================================
            ROSYMM(INDEX) = ROINP(INDEX) + CROSUM
            IF (IPRINT.GE.2) THEN
              CDUM = ROSYMM(INDEX) / FLOAT(NC)
              IF (ABS(CDUM).GT.0.1 .OR. ABS(ROINP(INDEX)).GT.0.1) 
     &           WRITE (*,661) INDEX,IG1,IG2,IG3,ROINP(INDEX),CDUM
            END IF
            GO TO 100
C
199         ROSYMM(INDEX) = (0.0,0.0)
            IF (IPRINT.GE.2.AND.INIPRN.EQ.1) 
     &      WRITE (*,1113) IG1,IG2,IG3,IGL1,IGL2,IGL3,NN
100   CONTINUE
C
      DO 110 I = 1 , NPLWV
        ROSYMM(I) = ROSYMM(I) / FLOAT(NC)
        ROINP(I)  = ROSYMM(I)
110   CONTINUE
C
1111  FORMAT(5X,2(3I5,5X))
1112  FORMAT(5X,3F15.6)
1113  FORMAT(' IG:',3I4,' IMAGE OUT OF GRID:',3I4,' N=',I2)
661   FORMAT(1X,I5,1X,3I2,1X,2E11.4,2X,2E11.4)
 550  FORMAT(1X,'N,IC=',2I3/3(3F8.4))
99    FORMAT (' THERE IS SOMETHING WRONG WITH THE LIST OF '/
     +  'SYMMETRY OPERATIONS OF THE POINT GROUP, IB:'/4(/1X,12I3))
      RETURN
      END
