program breduce
      implicit none
      integer*4 nat,n,nob,m,i,j,k,no,noa,iargc
      character*10 s10
      real*8 sij,sp,vn,an
      real*8 ,allocatable::B(:,:),BT(:,:)
      integer*4 ,allocatable::IB(:)
      write(6,600)
600   format(/,' Reduces number of internal coordinates',/,/,
     1'      Input: B.MAT',/,
     1'             FILE.UMA',/,
     1'     Output: BN.MAT',/,
     1'             FILEN.UMA',/,/,
     1'     Usage: "breduce <m>", m .. number of conserved coordinates')
      write(6,*)
      if(iargc().ne.1)stop
      call getarg(1,s10)
      read(s10,*)m
      open(9,file='FILE.UMA',status='old')
      read(9,*)
      read(9,*)
      read(9,*)nat,nob
      close(9)
      n=3*nat
      allocate(B(nob,n),BT(nob,n),IB(nob))
      call rb(B,nob,n)
c     orthogonalize coordinates, mark non-redundant:
      noa=0
      no=0
      BT=B
      IB=0
      do 1 i=1,nob
      do 2 j=1,i-1
      sij=sp(i,j,BT,BT,nob,n)
      do 2 k=1,n
2     BT(i,k)=BT(i,k)-sij*BT(j,k)
      an=vn(BT,i,nob,n)
      if(an.gt.1.0d-3)then
       do 4 k=1,n
4      BT(i,k)=BT(i,k)/an
c      mark i-coordinate as useful:
       IB(i)=1
       no=no+1
       if(i.gt.m)noa=noa+1
      else
       do 12 k=1,n
12     BT(i,k)=0.0d0
      endif
1     continue 
      write(6,601)no,nob,no-noa,m,noa,3*nat-6,m+noa
601   format(/,i6,' of',i6,' non-redundant',/,i6,' within 1 ..',i6,/,
     1i6,' in the rest',/,i6,' = 3 x Nat -6',/,
     1i6,' new redundant recorded',/)
      deallocate(BT)
      allocate(BT(m+noa,n))
      k=0
      do 9 i=1,nob
      if(IB(i).eq.1.or.i.le.m)then
       k=k+1
       do 10 j=1,n
10     BT(k,j)=B(i,j)
      endif
9     continue
      call wrb(BT,m+noa,n,1.0d-8,'BN.MAT')
      call rwuma(nat,m,noa,nob,IB)
      end

      subroutine rwuma(nat,m,noa,nob,IB)
      implicit none
      integer*4 nat,m,noa,nob,i,j,nn,no,k,n1,n2,IB(nat)
      integer*4 ,allocatable::li(:)
      character*120 s120
      allocate(li(nat))
      open(9,file='FILE.UMA')
      open(8,file='FILEN.UMA')
      read(9,80)s120
      call wst(s120,8)
      read(9,80)s120
      call wst(s120,8)
 80   format(a120)
      read(9,*)
      write(8,801)nat,m+noa,((m+noa)*(m+noa))/2
 801  format(3i12)
      do 3 i=1,nat
      read(9,80)s120
  3   call wst(s120,8)
      no=0
      do 5 i=1,nob
      read(9,80)s120
      read(s120,*)k,n1,n2
      nn=0
      if(k.eq.4)nn=2
c     inter-molecular:
      if(k.ge.8.and.k.le.13)nn=2
c     monomolecular:
      if(k.ge.24.and.k.le.29)nn=1
      if(IB(i).eq.1.or.i.le.m)then
       no=no+1
       call ws(s120,8)
       write(8,802)no
 802   format(i6)
       if(k.ge.24.and.k.le.29)then
        read(9,*)(li(j),j=1,n1)
        do 8 j=1,n1
  8     write(8,803)li(j)
 803    format(i6,$)
        write(8,*)
       elseif(k.ge.8.and.k.le.13)then
        read(9,*)(li(j),j=1,n1)
        do 1 j=1,n1
  1     write(8,803)li(j)
        write(8,*)
        read(9,*)(li(j),j=1,n2)
        do 2 j=1,n2
  2     write(8,803)li(j)
        write(8,*)
       elseif(k.eq.14)then
c       first helical coordinate (from 14..23):
        read(9,*)(li(j),j=1,n1)
        do 4 j=1,n1
  4     write(8,803)li(j)
        write(8,*)
       else
        do 6 j=1,nn
        read(9,80)s120
 6      call wst(s120,8)
       endif
      else
       do 7 j=1,nn
 7     read(9,*)
      endif
  5   continue
      CLOSE(9)
      WRITE(8,*)' 0'
      DO 771 I=1,3*nat-6
771   WRITE(8,777)I,1.0d0
777   FORMAT(' 1',/,I5,F5.1)
      CLOSE(8)
      write(6,*)'FILEN.UMA done'
      return
      end


      subroutine ws(s,io)
      implicit none
      integer*4 ie,io,i
      character*(*) s
      do 1 ie=len(s),1,-1
 1    if(s(ie:ie).ne.' ')goto 2
 2    do 3 i=1,ie
 3    write(io,80)s(i:i)
 80   format(a1,$)
      return
      end

      subroutine wst(s,io)
      implicit none
      integer*4 ie,io,i
      character*(*) s
      do 1 ie=len(s),1,-1
 1    if(s(ie:ie).ne.' ')goto 2
 2    do 3 i=1,ie
 3    write(io,80)s(i:i)
 80   format(a1,$)
      write(io,*)
      return
      end


      subroutine wrb(B0,M,NA,tol,s)
      implicit none
      integer*4 M,NA,I,J,nn,ia,ix
      real*8 B0(M,NA),tol
      character*(*) s
      DO 155 I=1,M
      DO 155 J=1,NA                                                   
155   IF(DABS(B0(I,J)).LT.tol) B0(I,J) = 0.0D0                      
      OPEN(9,FILE=s)
      WRITE(9,901)M,NA
901   format(I6,' Internal coordinates',/,
     1I6,' Cartesian coordinates',/,
     1' Transformation matrix, non-zero elements:')
      do 130  I=1,M
      WRITE(9,9011)I
9011  format(I6,' internal:')
      nn=0
      do 140 J=1,NA
140   if(DABS(B0(I,J)).gt.tol)nn=nn+1
      write(9,9012)nn
9012  format(I6)
      nn=0
      do 130 J=1,NA
      if(DABS(B0(I,J)).gt.tol)then
       nn=nn+1
       ia=(J+2)/3
       ix=J-3*(ia-1)
       WRITE(9,9100)ia,ix,B0(I,J)
9100   format(2I5,F15.8)
      endif
130   continue
      WRITE(*,*)s,' written'
      CLOSE(9)
      return
      end

      function vn(t,i,nob,n)
      integer*4 n,i,nob,j
      real*8 t(nob,n),vn,s
      s=0.0d0
      do 1 j=1,n
 1    s=s+t(i,j)**2
      vn=dsqrt(s)
      return
      end

      function sp(i,j,B,O,nob,n)
      implicit none
      integer*4 i,j,k,nob,n
      real*8 B(nob,n),s,O(nob,n),sp
      s=0.0d0
      do 1 k=1,n
  1   s=s+B(i,k)*O(j,k)
      sp=s
      return
      end



      subroutine rb(b,nob0,na)
      implicit none
      integer*4 nob0,na,i,nn,ii,ix,ia,ic
      real*8 b(nob0,na)
      b=0.0d0
      open(9,file='B.MAT')
      read(9,*)ic
      if(ic.ne.nob0)then
       write(6,*)ic,nob0,' inconsistent number of internals'
       close(9)
       stop
      endif
      read(9,*)na
      read(9,*)
      do 13  i=1,nob0
      read(9,*)
      read(9,*)nn
      do 13 ii=1,nn
13    read(9,*)ia,ix,b(i,ix+3*(ia-1))
      write(6,*)'B.MAT read'
      close(9)
      return
      end