PROGRAM chkproj
C     error from translation noninvariance
c     works in cartesian coordinates, for one molecule only
c     TRED12 diagonalization replaced by JACOBI_EIGENVALUE
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      real*8,allocatable::FCAR(:,:),A(:,:)
C
      OPEN(17,FILE='FILE.UMA')
      READ(17,*)
      READ(17,*)
      READ(17,*)NOAT
      NA=3*NOAT
      allocate(FCAR(NA,NA),A(NA,NA))
C
C     Read-in cartesian FF
      OPEN(20,FILE='FILE.FC',STATUS='OLD')
      CALL READFF(NA,FCAR)
      CLOSE(20)
C
      CALL PROJECT(A,FCAR,NA)
      end
C
      SUBROUTINE READFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(4X,5D14.6)
C
      CONST=4.359828d0/0.5291772d0/0.5291772d0
      DO 3 I=1,N
      DO 3 J=I,N
3     FCAR(J,I)=FCAR(J,I)*CONST
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 FCAR(I,J)=FCAR(J,I)
      WRITE(*,*)' Cartesian FF read in ... '
      RETURN
      END

C
      SUBROUTINE TRERR(F,N)
      IMPLICIT none
      INTEGER*4 N,NAT,I,ia
      real*8 F(N,N),rr,s1,s2,s3
      rr=0.0d0
      NAT=N/3
      do 1 I=1,N
      s1=0.0d0
      s2=0.0d0
      s3=0.0d0
      do 2 ia=2,nat
      s1=s1-F(I,1+3*(ia-1))
      s2=s2-F(I,2+3*(ia-1))
2     s3=s3-F(I,3+3*(ia-1))
1     rr=rr+(s1-F(I,1))**2+(s2-F(I,2))**2+(s3-F(I,3))**2
      write(6,600)rr
600   format(' Translational FF error ',f20.10)
      return
      end

      SUBROUTINE PROJECT(A,F,N)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(N,N),F(N,N)
      real*8,allocatable::TEM(:,:)
C
      call TRERR(F,N)
      allocate(TEM(N,N))
      WRITE(*,*)'Projecting transl/rotations from the force field...'
      CALL DOMA(A,N)
      DO 1 I=1,N
      DO 1 J=1,N
      TEM(I,J)=0.0d0
      DO 1 II=1,N
1     TEM(I,J)=TEM(I,J)+A(I,II)*F(II,J)
      DO 2 I=1,N
      DO 2 J=1,N
      F(I,J)=0.0d0
      DO 2 II=1,N
2     F(I,J)=F(I,J)+TEM(I,II)*A(J,II)
      WRITE(*,*)'  ... done.'
      call TRERR(F,N)
      RETURN
      END

      SUBROUTINE DOMA(A,N)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(N,N)
      real*8,allocatable::C(:,:),TEM(:,:)
C
      N=N
      allocate(C(3,N/3),TEM(N,N))
      AMACH=0.00000000000001d0
      NAT=N/3
      OPEN(4,FILE='FILE.X',STATUS='OLD')
      READ(4,*)
      READ(4,*)NAT
      DO 11 I=1,NAT
11    READ(4,*)C(1,I),(C(II,I),II=1,3)
      CLOSE(4)
      DO 12 I=1,N
      DO 12 J=1,6
12    A(I,J)=0.0d0
      DO 1 I=1,NAT
      A(3*I-2,1)= 1.00d0
      A(3*I-1,2)= 1.00d0
      A(3*I  ,3)= 1.00d0
      A(3*I-2,4)=-C(2,I)
      A(3*I-1,4)= C(1,I)
      A(3*I-1,5)=-C(3,I)
      A(3*I  ,5)= C(2,I)
      A(3*I-2,6)= C(3,I)
1     A(3*I  ,6)=-C(1,I)
      N6=6
      CALL ORT(A,N6,NAT)
      ICOL=1
5     IF(ICOL.GT.N6)GOTO 8
      S=0.0d0
      DO 6 J=1,N
6     S=S+A(J,ICOL)**2
      IF(S.LT.1000.0d0*AMACH)THEN
       DO 7 J=1,N
7      A(J,ICOL)=A(J,N6)
       N6=N6-1
       GOTO 5
      ENDIF
      ICOL=ICOL+1
      GOTO 5
8     CONTINUE
      WRITE(*,*)N6,' coordinates projected'
      DO 2 I=1,N
      DO 2 J=1,N
      TEM(I,J)=0.0d0
      DO 2 II=1,N6
2     TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II)
      DO 3 I=1,N
      DO 4 J=1,N
4     A(I,J)=-TEM(I,J)
3     A(I,I)=A(I,I)+1.0d0
      RETURN
      END
C     ============================================================ 
      SUBROUTINE ORT(A,N6,NAT)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(3*NAT,N6)
      AMACH=0.00000000000001d0
      NAT3=3*NAT
      C=0.0d0
      DO 1 I=1,NAT3
1     C=C+A(I,1)**2
      IF(C.GT.AMACH)C=1.0d0/SQRT(C)
      DO 2 I=1,NAT3
2     A(I,1)=A(I,1)*C
      DO 3 IC=2,N6
C     Column IC, orthonormalize to previous columns:
      DO 3 ICP=1,IC-1
      SP=0.0d0
      DO 5 J=1,NAT3
5     SP=SP+A(J,IC)*A(J,ICP)
C     Subtract projection to ICP:
      DO 6 J=1,NAT3
6     A(J,IC)=A(J,IC)-SP*A(J,ICP)
C     Normalize :
      SP=0.0d0
      DO 7 J=1,NAT3
7     SP=SP+A(J,IC)**2
      IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP)
      DO 3 J=1,NAT3
3     A(J,IC)=A(J,IC)*SP
      RETURN
      END