program clustersort
      implicit none
      integer*4 m,nat,i,ia,n,nc,ii,imin,ip,j,icy,
     1nt,mt,np,iu,mint,numin,iumin,nu,iustart,iuend
      real*8 da,dd,dmin,dist,ta
      real*8,allocatable::r(:,:),tau(:,:)
      integer*4,allocatable::lp(:),npp(:),
     1list(:,:),al(:),il(:)
      character*1 ,allocatable::let(:)
      character*80 ,allocatable::fn(:)
      write(6,600)
 600  format(/,' Sorts geometries into N groups',/,/,
     1         '   Input:   STRUCT.LST   list of geometries',/,
     1         '            CLUSTER.OPT  options',/)

      mint=1
      allocate(al(1),il(1))
      call readopt(0,nat,nc,al,mint,il)
      deallocate(al,il)
      allocate(al(nat),il(4*mint+1))
      call readopt(1,nat,nc,al,mint,il)
      n=3*nat

c     load structures and put them to geometrical centers
      allocate(r(1,1),let(1),fn(1))
      call ldstr(0,m,nat,n,r, 1, 1,let,al,fn)
      deallocate(r,let,fn)
      allocate(r(m,n),let(m),fn(m))
      mt=m
      nt=n
      call ldstr(1,m,nat,n,r,mt,nt,let,al,fn)

      if(mint.gt.0)then
c      calculate internal coordinates = torsion angles:
       allocate(tau(m,mint))
       do 35 i=1,m
       do 35 j=1,mint
35     tau(i,j)=ta(il(1+4*(j-1)),il(2+4*(j-1)),
     1             il(3+4*(j-1)),il(4+4*(j-1)),r,m,n,i)
      else
c      rotate to first geometry:
       call rgeo(m,n,nat,r)
      endif

      if(nc.lt.2)call report('too few points')
c     one by one building
      allocate(lp(m),npp(nc),list(nc,m))

      do 2000 icy=1,2
      numin=2**30
      if(icy.eq.1)then
       iustart=1
       iuend=m
       iumin=1
      else
       iustart=iumin
       iuend=iustart
      endif
      do 1000 iu=iustart,iuend
c     start from structure iu:
      np=1
      lp(1)=iu
      if(icy.eq.2)write(6,614)np,iu,let(iu)
77    da=0.0d0
      ia=0
c     find farthest structure not included:
      do 22 i=1,m
      do 25 j=1,np
25    if(lp(j).eq.i)goto 22
c     minimal distance to any point:
      dmin=1.0d10
      do 251 j=1,np
c     distance to this point
      dd=0.0d0
      if(mint.eq.0)then
       do 26 ii=1,n
26     dd=dd+(r(i,ii)-r(lp(j),ii))**2
      else
       do 261 ii=1,mint
261    dd=dd+(dist(tau(i,ii),tau(lp(j),ii)))**2
      endif
251   if(dd.lt.dmin)dmin=dd
      if(dmin.gt.da)then
       da=dmin
       ia=i
      endif
22    continue
      if(ia.eq.0)call report('Farthest structure not found')
      np=np+1
      lp(np)=ia
      if(icy.eq.2)write(6,614)np,ia,let(ia)
614   format(' Point ',i3,1x,' structure ;',i6,' ',A1,2f10.2)
      if(np.lt.nc)goto 77

c     npp ... number of geometries in each ref. point
      npp=0
c     loop over geometries:
      do 27 i=1,m
c     select point with the smallest distance:
      dmin=1.0d100
      imin=1
      do 28 ip=1,np
      dd=0.0d0
      if(mint.eq.0)then
       do 29 ii=1,n
 29    dd=dd+(r(i,ii)-r(lp(ip),ii))**2
      else
       do 291 ii=1,mint
 291   dd=dd+(dist(tau(i,ii),tau(lp(ip),ii)))**2
      endif
      if(dd.lt.dmin)then
       dmin=dd
       imin=ip
      endif
28    continue
c     assign to it:
      npp(imin)=npp(imin)+1
27    list(imin,npp(imin))=i    

      if(icy.eq.2)then
       do 30 i=1,np
       write(6,611)i
 611   format(i6,':',$)
       do 34 j=1,npp(i)
34     write(6,613)let(list(i,j))
 613   format(a1,$)
30     write(6,*)
       write(6,*)
       do 40 i=1,np
       write(6,611)i
       do 41 j=1,npp(i)
41     write(6,615)list(i,j)
 615   format(i5,$)
40     write(6,*)
       write(6,*)
       call wrfn(fn,nc,m,npp,list)
       write(6,*)'Best distrubution ',iumin
      endif

      nu=0
      do 39 i=1,np
39    nu=nu+npp(i)**2
      if(nu.lt.numin)then
       numin=nu
       iumin=iu
      endif

1000  continue
2000  continue


      end

      subroutine wrfn(fn,nc,m,npp,list)
      implicit none
      character*80 fn(*)
      character*20 s20
      integer*4 nc,m,npp(*),list(nc,m),i,is,j,l,k,i1,i2
      do 1 i=1,nc
      write(s20,20)i
 20   format(i20)
      do 2 is=1,len(s20)
2     if(s20(is:is).ne.' ')goto 3
3     open(9,file=s20(is:len(s20))//'.lst')
      do 4 j=1,npp(i)
      l=list(i,j)
      do 5 i1=1,len(fn(l))
5     if(fn(l)(i1:i1).ne.' ')goto 6
6     do 7 i2=len(fn(l)),1,-1
7     if(fn(l)(i2:i2).ne.' ')goto 8
8     continue
      do 9 k=i1,i2
9     write(9,90)fn(l)(k:k)
90    format(a1,$)
4     write(9,*)
      close(9)
1     write(6,*)s20(is:len(s20))//'.lst'
      return
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine wr(XB,i,r,m,nat)
      implicit none
      integer*4 m,nat,ia,i,ix
      real*8 XB(nat,3),r(m,3*nat)
      do 1 ia=1,nat
      do 1 ix=1,3
1     XB(ia,ix)=r(i,ix+3*(ia-1))
      return
      end

      SUBROUTINE DOMATRIX(IERR,A,XB,XS,nat)
C     This is a Fortran version of Petr Malon's subroutine Amoeba
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3),A(3,3),
     1XS(nat,3),XB(nat,3)
      A=0.0d0
      IERR=0
C     STARTING VALUES FOR THE ITERATION:
      ANG(1)=0.2d0
      ANG(2)=0.1d0
      ANG(3)=0.0d0
      DO 1 I=1,3
1     P(1,I)=ANG(I)
      Y(1)=FU(nat,ANG,A,XS,XB)
      ANG(2)=ANG(1)
      ANG(1)=0.0d0
      DO 2 I=1,3
2     P(2,I)=ANG(I)
      Y(2)=FU(nat,ANG,A,XS,XB)
      ANG(3)=ANG(2)
      ANG(2)=0.0d0
      DO 3 I=1,3
3     P(3,I)=ANG(I)
      Y(3)=FU(nat,ANG,A,XS,XB)
      ANG(1)=ANG(3)
      ANG(2)=ANG(3)
      ANG(3)=0.0d0
      DO 4 I=1,3
4     P(4,I)=ANG(I)
      Y(4)=FU(nat,ANG,A,XS,XB)
      FTOL=1.0d-7
      ITMAX=5000
      ITER=0
99999 ILO=1
      IF(Y(1).GT.Y(2))THEN
       IHI=1
       INHI=2
      ELSE
       IHI=2
       INHI=1
      ENDIF
      DO 5 I=1,4
      IF(Y(I).LT.Y(ILO))ILO=I
      IF(Y(I).GT.Y(IHI))THEN
       INHI=IHI
       IHI=I
      ELSE
       IF(Y(I).GT.Y(INHI))THEN
        IF(I.NE.IHI)INHI=I
       ENDIF
      ENDIF
5     CONTINUE
      IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.1.0d-10)THEN
       RTOL=0.0d0
      ELSE
       RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0
      ENDIF
      IF(RTOL.LT.FTOL)then
       RETURN
      endif
      IF(ITER.EQ.ITMAX)THEN
       write(6,*)'nat',nat
       write(6,*)1,XS(1,1),XS(1,2),XS(1,3),XB(1,1),XB(1,2),XB(1,3)
       write(6,*)nat,XS(nat,1),XS(nat,2),XS(nat,3),
     1 XB(nat,1),XB(nat,2),XB(nat,3)
       write(6,*)'A',A
       WRITE(6,*)' Rotation has not converged !',ITER
       IERR=1
       RETURN
      ENDIF
      ITER=ITER+1
      DO 55 I=1,3
55    PBAR(I)=0.0d0
      DO 6 I=1,4
      IF(I.NE.IHI)THEN
       DO 7 J=1,3
7      PBAR(J)=PBAR(J)+P(I,J)
      ENDIF
6     CONTINUE
      DO 8 J=1,3
      PBAR(J)=PBAR(J)/3.0d0
8     PR(J)=2.0d0*PBAR(J)-P(IHI,J)
      YPR=FU(nat,PR,A,XS,XB)
      IF(YPR.LE.Y(ILO))THEN
       DO 9 J=1,3
9      PRR(J)=2.0d0*PR(J)-PBAR(J)
       YPRR=FU(nat,PRR,A,XS,XB)
       IF(YPRR.LT.Y(ILO))THEN
        DO 10 J=1,3
10      P(IHI,J)=PRR(J)
        Y(IHI)=YPRR
       ELSE
        DO 11 J=1,3
11      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ELSE
       IF(YPR.GE.Y(INHI))THEN
        IF(YPR.LT.Y(IHI))THEN
         DO 12 J=1,3
12       P(IHI,J)=PR(J)
         Y(IHI)=YPR
        ENDIF
        DO 13 J=1,3
13      PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J)
        YPRR=FU(nat,PRR,A,XS,XB)
        IF(YPRR.LT.Y(IHI))THEN
         DO 14 J=1,3
14       P(IHI,J)=PRR(J)
         Y(IHI)=YPRR
        ELSE
         DO 15 I=1,4
         IF(I.NE.ILO)THEN
          DO 16 J=1,3
          PR(J)=0.5d0*(P(I,J)+P(ILO,J))
16        P(I,J)=PR(J)
          Y(I)=FU(nat,PR,A,XS,XB)
         ENDIF
15       CONTINUE
        ENDIF
       ELSE
        DO 17 J=1,3
17      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ENDIF
      GOTO 99999
      END

      FUNCTION FU(n,U,A,XS,XB)
      implicit none
      INTEGER*4 n,I
      REAL*8 S1,S2,S3,C1,C2,C3,TX,TY,TZ,U(3),FU,
     1A(3,3),XS(n,3),XB(n,3),RMS
      S1=SIN(U(1))
      S2=SIN(U(2))
      S3=SIN(U(3))
      C1=COS(U(1))
      C2=COS(U(2))
      C3=COS(U(3))
C     u1 = theta
C     u2 = phi
C     u3 = cappa
      A(1,1)=C1*C2*C3-S2*S3
      A(1,2)=C1*S2*C3+C2*S3
      A(1,3)=-S1*C3
      A(2,1)=-C1*C2*S3-S2*C3
      A(2,2)=-C1*S2*S3+C2*C3
      A(2,3)=S1*S3
      A(3,1)=S1*C2
      A(3,2)=S1*S2
      A(3,3)=C1
      RMS=0.0d0
      DO 1 I=1,n
      TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3)
      TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3)
      TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3)
1     RMS=RMS+(TX-XB(I,1))*(TX-XB(I,1))+
     1        (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3))
      FU=RMS
      RETURN
      END


      subroutine ldstr(ic,m,nat,n,r,m1,n1,let,al,fn)
c     loading structures
      implicit none
      integer*4 ic,m,nat1,nat,m1,n1,ia,i,ix,n,nata,al(*),il
      real*8 r(m1,n1),c
      character*80 fn(*)
      character*1 let(*)
      open(9,file='STRUCT.LST',status='old')
      if(ic.eq.0)then
       m=0
 1     read(9,90,end=99,err=99)fn(1)
90     format(a80)
       m=m+1
       if(m.eq.1)then
        open(90,file=fn(1),status='old')
        read(90,*)
        read(90,*)nat1
        close(90)
       endif
       goto 1
 99    write(6,60)m,nat1,nat
       n=3*nat
60     format(i6,' geometries,',i6,' atoms in first,',i6,' atoms read')
      else
       do 2 i=1,m
       read(9,90)fn(i)
       do 101 ia=1,len(fn(i))
       if(fn(i)(ia:ia).ne.' ')then
        let(i)=fn(i)(ia:ia)
        goto 102
       endif
101    continue
102    open(90,file=fn(i),status='old')
       read(90,*)
       read(90,*)nata
       do 3 ia=1,nata
       do 31 il=1,nat
       if(al(il).eq.ia)then
        read(90,*)(r(i,ix+3*(il-1)),ix=1,1),(r(i,ix+3*(il-1)),ix=1,3)
        goto 3
       endif
 31    continue
       read(90,*)
3      continue
c      put it to the geometrical center:
       do 4 ix=1,3
       c=0.0d0
       do 5 ia=1,nat
5      c=c+r(i,ix+3*(ia-1))
       c=c/dble(nat)
       do 4 ia=1,nat
4      r(i,ix+3*(ia-1))=r(i,ix+3*(ia-1))-c
2      close(90)
      endif
      close(9)
      return
      end

      subroutine rgeo(m,n,nat,r)
c     rotate all to first:
      implicit none
      integer*4 m,n,nat,i,is,ix,IERR
      real*8 r(m,n),A(3,3)
      real*8 ,allocatable::XB(:,:),XS(:,:)
      character*20 s20
      allocate(XB(nat,3),XS(nat,3))
c     first geometry in XB:
      do 1 i=1,nat
      do 1 ix=1,3
1     XB(i,ix)=r(1,ix+3*(i-1))

c     rotate each geometry to XB:
      do 4 is=1,m
      call wr(XS,is,r,m,nat)
      CALL DOMATRIX(IERR,A,XB,XS,nat)
      IF(IERR.EQ.1)CALL REPORT('CRMs DoMatrix error')
      write(s20,200)is
200   format(i20)
      do 4 I=1,nat
      do 4 ix=1,3
4     r(is,ix+3*(I-1))=A(ix,1)*XS(I,1)+A(ix,2)*XS(I,2)+A(ix,3)*XS(I,3)

      write(6,*)'geometries rotated to first'
      return
      end

      subroutine readopt(ic,nat,nc,al,mint,il)
      implicit none
      integer*4 ic,nc,nat,al(*),i,j,mint,il(*)
      character*3 key
c     number of atoms read (not all):
      if(ic.eq.0)nat=0
c     numeber of conformer classesL
      nc=0
c     use mint internal coordinates=torsion angles:
      mint=0
      open(9,file='CLUSTER.OPT',status='old')
1     read(9,90,end=99,err=99)key
90    format(a3)
      if(key.eq.'NAT')read(9,*)nat
      if(key.eq.'CLA')read(9,*)nc
      if(key.eq.'LIS'.and.ic.eq.1)read(9,*)(al(i),i=1,nat)
      if(key.eq.'INT'.and.ic.eq.0)read(9,*)mint
      if(key.eq.'INT'.and.ic.eq.1)
     1read(9,*)mint,((il(j+4*(i-1)),j=1,4),i=1,mint)
      goto 1
99    close(9)
      return
      end

      function ta(i1,i2,i3,i4,r,m,n,i)
      implicit none
      integer*4 i1,i2,i3,i4,m,n,i,ix,KK
      real*8 r(m,n),A(3),P(3),E0(3),B0(3),SSB0B0,AP(3),
     1SSE0E0,SSAE0,SSPE0,XX(3),YY(3),SSXXXX,sp,SSE0AP,SIGN,CA,
     1ANGLE,PI,SSYYYY,SSXXYY,ta
      PI=4.0d0*datan(1.0d0)
      do 23 ix=1,3
      A(ix) =r(i,ix+3*(i3-1))-r(i,ix+3*(i4-1))
      P(ix) =r(i,ix+3*(i1-1))-r(i,ix+3*(i2-1))
23    E0(ix)=r(i,ix+3*(i2-1))-r(i,ix+3*(i3-1))
      call vp(B0,E0,A)
      SSB0B0=sp(B0,B0)
      SSE0E0=sp(E0,E0)
      DO 8 KK=1,3
      B0(KK)=B0(KK)/DSQRT(SSB0B0)
8     E0(KK)=E0(KK)/DSQRT(SSE0E0)
      SSAE0=sp(A,E0)
      SSPE0=sp(P,E0)
      DO 9 KK=1,3
      XX(KK)=A(KK)-E0(KK)*SSAE0
9     YY(KK)=P(KK)-E0(KK)*SSPE0
      SSXXXX=sp(XX,XX)
      SSYYYY=sp(YY,YY)
      SSXXYY=sp(XX,YY)
      call vp(AP,XX,YY)
      SSE0AP=sp(AP,E0)
      SIGN=1.0D0
      IF (SSE0AP.NE.0.0D0)SIGN=-SSE0AP/ABS(SSE0AP)
      CA=0.0D0
      IF (SSXXXX.NE.0.0D0.AND.SSYYYY.NE.0.0D0)
     1CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY)
      ANGLE=0.0D0
      IF (CA.LE.-1.0D0) ANGLE=PI
      IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA)
      ANGLE=SIGN*ANGLE
      if(ANGLE.lt.0.0d0)ANGLE=2.0d0*PI+ANGLE
      ta=ANGLE
      return
      end

      function sp(a,b)
      IMPLICIT none
      real*8 a(*),b(*),sp
      sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
      return
      end

      subroutine vp(A,X,Y)
      implicit none
      real*8 A(*),X(*),Y(*)
      A(1)=X(2)*Y(3)-X(3)*Y(2)
      A(2)=X(3)*Y(1)-X(1)*Y(3)
      A(3)=X(1)*Y(2)-X(2)*Y(1)
      return
      end

      function dist(a,b)
      implicit none
      real*8 a,b,pi,d,dist
      pi=4.0d0*datan(1.0d0)
      d=a-b
      if(dabs(d+2.0d0*pi).lt.dabs(d))d=d+2.0d0*pi
      if(dabs(d-2.0d0*pi).lt.dabs(d))d=d-2.0d0*pi
      dist=d
      return
      end