C=========================================================================
C     SETUP CONSTRAINTS WHICH ARE SUPPOSED TO BE LINEAR WITH
C     RESPECT TO COORDINATES. THEN INSTEAD OF THE PHYSICAL COORDINATES
C     Ri,i=1,3*N, WE INTRODUCE A SET OF AUXILIARY COORDINATES, 
C     Pi,i=1,M, M < 3*N. SUBROUTINE CONSTR SETS UP THE MATRIX OF
C     DERIVATIVES RCONSTR(i,j)=dRi/dPj.     
C
C     NCNSTR - NUMBER OF CONSTRAINTS 
C     NEWCOO = M = 3 * N - NCNSTR
C            __
C      dRi = >  RCNSTR(i,j) * dPj ,  i=1,3*NIONST
C            --
C          j=1,M
C
C     THIS RELATION IS WRITTEN IN CARTESIAN COORDINATES FOR 
C     DISPLACEMENTS (dR and dP).
C     FOR ARBITRARY BASIS VECTORS THE SYMMETRY WE ARE LOOKING FOR 
C     WILL BE BETTER SEEN IN THE CARTESIANS THAN IN UNITS OF 
C     BASIS VECTORS.
C=========================================================================
C
      SUBROUTINE CONSTR(NSPEC,NIONS,NIONST,RCNSTR,DIRC,POSIC,
     &                  RMOVE,NCNSTR,NEWCOO,WORK)
C
      IMPLICIT COMPLEX (C)
      DIMENSION RCNSTR(3*NSPEC*NIONS,3*NSPEC*NIONS)
      DIMENSION DIRC(3,3)
      DIMENSION POSIC(3,NSPEC*NIONS)
      DIMENSION WORK(3,*)
      DIMENSION RMOVE(NSPEC*NIONS)
C
      NDIM = 3 * NSPEC * NIONS
      NCOO = 3 * NIONST
C
C=========================================================================
C  EXAMPLE - NO CONSTRAINT WHATSOEVER (EVEN RMOVE)
C=========================================================================
C
C      NEWCOO = NCOO
C      NCNSTR = 0
C      DO 100 I = 1 , NDIM
C        DO 101 J = 1 , NDIM
C          RCNSTR(I,J) = 0.0
C 101    CONTINUE
C 100  CONTINUE
C      DO 102 I = 1 , NCOO
C        RCNSTR(I,I) = 1.0
C 102  CONTINUE
C
C=========================================================================
C  EXAMPLE - CONSTRAINT ONLY IN RMOVE
C=========================================================================
C
      DO 200 I = 1 , NDIM
        DO 201 J = 1 , NDIM
          RCNSTR(I,J) = 0.0
 201    CONTINUE
 200  CONTINUE
C
      NCNSTR = 0 
      NINDP = 1
      DO 204 NI = 1 , NIONST
        IF (RMOVE(NI).NE.0.0) THEN
          DO 202 K =  1 , 3      
            RCNSTR(K+3*(NI-1),NINDP) = 1.0
            NINDP = NINDP + 1
 202      CONTINUE
        ELSE
          NCNSTR = NCNSTR + 3
        END IF
 204  CONTINUE
      NEWCOO = NCOO - NCNSTR
      WRITE (*,*) ' CONSTR in action'
      WRITE (*,1000) NCOO,NEWCOO,NCNSTR,NIONST
      DO 210 J = 1 , NDIM
        DO 211 I = 1 , NDIM
          IF (RCNSTR(I,J).NE.0.0) WRITE (*,1001) I,J,RCNSTR(I,J)
 211    CONTINUE
 210  CONTINUE
 1000 FORMAT(' Old coord:   ',I3,' New coord: ',I3/
     1       ' Constraints: ',I3,' Ions     : ',I3/
     2       '   CONSTRAINTS MATRIX (Rij)')
 1001 FORMAT(1X,'RCNSTR(',I3,',',I3,')=',F6.2)
 1002 FORMAT(' DCNSTR(',I3,')=',3F12.5)
      RETURN 
      END
