program sortconf
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
c
c     looks at gaussian outputs, finds final structures
c     and sorts them out
c     ns0 ... number of structures
c     e(*) energies
c     ind(*) index mapping
c     re(*) relative energies
c     ig(*) number of geometries found
c     ifollow  .. number of parameters to follow (nt0 max)
c
      parameter (ns0=10000,nat0=200,nt0=20)
      character*80 ts,s80,tso
      character*10 name(ns0)
      character*3 o(ns0)
      dimension e(ns0),ind(ns0),re(ns0),ig(ns0),ie(ns0),r(3*nat0*ns0),
     1iz(nat0),i4(nt0,4),afj(nt0),af(nt0,ns0),ib1(nat0),ib2(nat0),
     2bo(nat0)
      logical lex


      ic=0
      io=0
      write(6,700)
700   format('Looking at conformers in the list NAMES.LST')
      open(3,file='NAMES.LST',status='old')
1     read(3,3000,end=3333,err=3333)ts
3000  format(a80)
      write(6,*)ts(1:40)
      ic=ic+1
      if(ic.gt.ns0)then
       write(6,*)'too many structures'
       stop
      endif
      name(ic)=ts(1:10)
      ie(ic)=0
      o(ic)=' No'
      ig(ic)=0
      ist=1

      do 87 i=1,80
87    tso(i:i)=' '
      open(31,file=ts)
11    read(31,3000,err=33331,end=33331)ts

      if(ts(1:9).eq.' SCF Done'.or.ts(1:9).eq.' Energy= ')then
       ie(ic)=ie(ic)+1
       do 3 is=1,len(ts)
3      if(ts(is:is).eq.'=')goto 4
4      open(32,file='scr')
       write(32,*)ts(is+1:len(ts))
       rewind 32
       read(32,*)e(ic)
       close(32)
      endif

      if(ts(28:50).eq.'!   Optimized Parameter'.or.
     1   ts(26:48).eq.'!   Optimized Parameter'.or.
     1   ts( 1:24).eq.' Optimization completed.')then
       io=io+1
       o(ic)='Yes'
      endif 

      IF(ts(19:39).EQ.'Z-Matrix orientation:'.OR.
     1   ts(26:46).EQ.'Z-Matrix orientation:'.OR.
     2   ts(20:37).EQ.'Input orientation:'.OR.
     3   ts(27:44).EQ.'Input orientation:')then
       ist=0
       call readgeox(ts,s80,r,iz,nat,ig,ic,nat0)
      ENDIF
      IF(ist.eq.1.and.
     1  (ts(20:40).EQ.'Standard orientation:'.OR.
     2   ts(26:46).EQ.'Standard orientation:'))
     3   call readgeox(ts,s80,r,iz,nat,ig,ic,nat0)

      do 3002 i=1,71
      if(ts(i:i+1).eq.'%%'.or.(ts(1:1).eq.'%'.and.tso(71:71).eq.'%'))
     1then
       iu=i
       call readgeoa(ts,r,nat,ig,ic,nat0,iu)
       goto 3003
      endif
3002  continue
3003  continue


      tso=ts
      goto 11 
33331 close(31)

      goto 1
3333  close(3)

      inquire(file='SORT.MOL',exist=lex)
      if(lex)then
c      mcm .mol format + accuracy in the second line
       open(3,file='SORT.MOL')
       read(3,*)
       read(3,*)ibonds,tolb
       read(3,*)
       read(3,*)
       if(ibonds.gt.nat0)then
        write(6,*)'Too many bonds'
	stop
       endif
       do 1013 i=1,ibonds
 1013  read(3,*)ib1(i),ib2(i)
       close(3)
       write(6,560)ibonds
  560  format(i4,' bonds will be monitored')
      else
       ibonds=0
      endif

      inquire(file='SORT.OPT',exist=lex)
      if(lex)then
       open(3,file='SORT.OPT')
       read(3,*)ifollow,tol
       if(ifollow.gt.nt0)then
        write(6,*)'Too many torsions'
	stop
       endif
       do 101 i=1,ifollow
 101   read(3,*)(i4(i,ii),ii=1,4)
       close(3)
      else
       ifollow=0
      endif

      nc=ic

      if(ibonds.gt.0)then
       write(6,6005)ibonds
6005   format(' Determining ',i4,' bonds from the first structure')
       do 501 i=1,ibonds
       i1=ib1(i)
       i2=ib2(i)
       x1=r(1+3*(i1-1))
       y1=r(2+3*(i1-1))
       z1=r(3+3*(i1-1))
       x2=r(1+3*(i2-1))
       y2=r(2+3*(i2-1))
       z2=r(3+3*(i2-1))
501    bo(i)=dsqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
  
       j=1
5029   jo=3*nat0*(j-1)
       do 502 i=1,ibonds
       bold=bo(i)
       i1=ib1(i)
       i2=ib2(i)
       x1=r(1+3*(i1-1)+jo)
       y1=r(2+3*(i1-1)+jo)
       z1=r(3+3*(i1-1)+jo)
       x2=r(1+3*(i2-1)+jo)
       y2=r(2+3*(i2-1)+jo)
       z2=r(3+3*(i2-1)+jo)
       bnew=dsqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
       if(abs(bnew-bold).gt.tolb)then
        write(6,399)j,name(j),i,bnew,bold
399     format(i5,1x,a10,' bond ',i3,' broken, new ',f5.2,', old ',
     1f5.2,/,' structure deleted')
        do 253 k=j,ic-1
        e(k)=e(k+1)
        o(k)=o(k+1)
        ig(k)=ig(k+1)
        name(k)=name(k+1)
        do 253 it=1,3*nat
253     r(it+3*nat0*(k-1))=r(it+3*nat0*k)
        ic=ic-1
	j=j+1
	if(j.le.ic)goto 5029
       endif
502    continue
       j=j+1
       if(j.le.ic)goto 5029
      endif

      if(ifollow.gt.0)then
       i=1
29     iw=3*nat0*(i-1)
       do 211 ii=1,ifollow
       i0=i4(ii,1)
       i1=i4(ii,2)
       i2=i4(ii,3)
       i3=i4(ii,4)
211    af(ii,i)=torsion(r(1+3*(i0-1)+iw),r(1+3*(i1-1)+iw),
     1 r(1+3*(i2-1)+iw),r(1+3*(i3-1)+iw),r(2+3*(i0-1)+iw),
     1 r(2+3*(i1-1)+iw),r(2+3*(i2-1)+iw),r(2+3*(i3-1)+iw),
     1 r(3+3*(i0-1)+iw),r(3+3*(i1-1)+iw),r(3+3*(i2-1)+iw),
     1 r(3+3*(i3-1)+iw),ad,bangle)

       j=i+1
22     if(j.le.ic)then
        jo=3*nat0*(j-1)
        do 212 jj=1,ifollow
        i0=i4(jj,1)
        i1=i4(jj,2)
        i2=i4(jj,3)
        i3=i4(jj,4)
212     afj(jj)=torsion(r(1+3*(i0-1)+jo),r(1+3*(i1-1)+jo),
     1  r(1+3*(i2-1)+jo),r(1+3*(i3-1)+jo),r(2+3*(i0-1)+jo),
     1  r(2+3*(i1-1)+jo),r(2+3*(i2-1)+jo),r(2+3*(i3-1)+jo),
     1  r(3+3*(i0-1)+jo),r(3+3*(i1-1)+jo),r(3+3*(i2-1)+jo),
     1  r(3+3*(i3-1)+jo),ad,bangle)
        sum=0.0d0
        do 24 kk=1,ifollow
24      sum=sum+(af(kk,i)-afj(kk))**2
        rms=dsqrt(sum/dble(ifollow))
        if(rms.lt.tol)then
         write(6,6007)name(i),name(j),rms,name(j)
6007     format(a10,' and ',a10,' same geom(',f7.4,') ',a10,' removed')
         do 25 k=j,ic-1
         e(k)=e(k+1)
         o(k)=o(k+1)
         ig(k)=ig(k+1)
         name(k)=name(k+1)
         do 25 it=1,3*nat
25       r(it+3*nat0*(k-1))=r(it+3*nat0*k)
         ic=ic-1
        else
         j=j+1
        endif
        if(j.le.ic)goto 22
       endif	
       i=i+1
       if(i.le.ic)goto 29
      endif

      emin=0.0d0
      emax=0.0d0
      do 6 i=1,ic
      if(ie(i).gt.0)then
       if(e(i).lt.emin.or.i.eq.1)emin=e(i)
       if(e(i).gt.emax.or.i.eq.1)emax=e(i)
      endif
6     continue
      q=0.0d0
      rt=0.592d0
      io=0
      do 7 i=1,ic
      if(o(i).eq.'Yes')io=io+1
      if(ie(i).gt.0)then
       ind(i)=i
       re(i)=(e(i)-emin)*627.5d0
       q=q+exp(-re(i)/rt)
      else
       re(i)=100000.0d0
       e(i)=10000.0d0
      endif
7     continue

      write(6,6000)nc,ic,io
6000  format(i8,' conformers ',i8,' unique ',i8,' optimized')

      do 8 i=1,ic
      ii=ind(i)
      do 8 j=i+1,ic
      jj=ind(j)
      if(e(jj).lt.e(ii))then
       it=ind(i)
       ind(i)=ind(j)
       ind(j)=it
       ii=ind(i)
      endif
8     continue

10    write(6,*)'How many conformers to list:'
      read(5,*)no
      if(no.gt.ic)then
       write(6,*)'Not so many left!'
       goto 10
      endif
      if(no.eq.0)goto 888
      open(61,file='SORT.OUT')
      write(6,6004)
      write(61,6004)
6004  format(' conformer rel. en. ratio  opt. #ens #geos',/,
     1       '           kcal/mol   %                        ')
      do 9 j=1,no
      i=ind(j)
      p=exp(-re(i)/rt)/q
      if(ifollow.gt.0)then
       write(61,6009)j,i,re(i),p*100.0d0,o(i),ie(i),ig(i),name(i),
     1 (af(io,i),io=1,ifollow)
       write(6,6009)j,i,re(i),p*100.0d0,o(i),ie(i),ig(i),name(i),
     1 (af(io,i),io=1,ifollow)
6009   format(i4,i6,f9.2,f6.2,3x,a3,i4,i4,2x,a10,20f7.2)
      else
       write(6,4004)j,i,re(i),p*100.0d0,o(i),ie(i),ig(i),name(i)
       write(61,4004)j,i,re(i),p*100.0d0,o(i),ie(i),ig(i),name(i)
4004   format(i4,i6,f9.2,f6.2,3x,a3,i4,i4,2x,a10)
      endif
      do 201 ir=1,len(name(i))
201   if(name(i)(ir:ir).eq.'.'.or.name(i)(ir:ir).eq.' ')goto 202
202   open(33,file=name(i)(1:ir-1)//'o.x')
      write(33,*)e(i)
      write(33,*)nat
      do 203 ia=1,nat
203   write(33,33339)iz(ia),(r(ix+3*(ia-1)+3*nat0*(i-1)),ix=1,3)
33339 format(i6,3f12.6)
      close(33)
9     continue
      close(61)

      if(no.gt.0)goto 10

888   stop
      end
     
      subroutine readgeox(ts,s80,r,iz,nat,ig,ic,nat0)
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      character*(*) ts,s80
      dimension r(*),iz(*),ig(*)
      ig(ic)=ig(ic)+1
      ig98=0
      if(ts(26:46).EQ.'Z-Matrix orientation:'.OR.
     1   ts(27:44).EQ.'Input orientation:'.OR.
     1   ts(26:46).EQ.'Standard orientation:')ig98=1
      DO 4 I=1,4
4     READ(31,*)
      l=0
5     READ(31,3000)s80
3000  FORMAT(a80)
      IF(s80(2:4).NE.'---')THEN
       l=l+1
       if(l.gt.nat0)then
        write(6,*)'Too many atoms'
	stop
       endif
       ll=3*(l-1)+3*nat0*(ic-1)
       BACKSPACE 31
       if(ig98.eq.0)then
        READ(31,*)idumm,iz(l),(r(i+ll),i=1,3)
       else
        READ(31,*)idumm,iz(l),idumm,(r(i+ll),i=1,3)
       endif
       GOTO 5
      ENDIF
      nat=l
      return
      end

       
      subroutine readgeoa(ts,r,nat,ig,ic,nat0,I)
c     reads geometry from the archive
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      character*1 L,LS,LS1
      character*80 ts
      character*1 NUMBER(20)
      dimension r(*),ig(*)
C
      LS='\\'
      LS1='|'
      ig(ic)=ig(ic)+1
      I=I+1
C     I .. index of the letter to be gotten
c     archive line has format (1X,A70)
      IF(I.GT.71)THEN
       READ(31,3000)ts
3000   format(a80)
       I=2
      ENDIF

      DO 73 IA=1,NAT
C     Get element and forget it right away:
      I=I+1
      CALL GETLET(L,ts,I)
      I=I+1
      CALL GETLET(L,ts,I)
      if(L.ne.',')then
       I=I+1
       CALL GETLET(L,ts,I)
      endif
      DO 73 IX=1,3
      DO 75 J=1,20
75    NUMBER(J)=' '
      IN=0
7558  IN=IN+1
      I=I+1
      CALL GETLET(L,ts,I)
      NUMBER(IN)=L
      IF(L.NE.','.AND.(L.NE.LS.and.L.NE.LS1))GOTO 7558
      NUMBER(IN)=' '
      open(8,file='scr')
      WRITE(8,*)(NUMBER(J),J=1,20)
      REWIND 8
      READ(8,*)r(IX+3*(IA-1)+3*nat0*(ic-1))
73    close(8)
      I=79
      return
      end
c
      SUBROUTINE GETLET(L,ts,I)
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*(*) ts
      CHARACTER*1 L
      IF(I.GT.71)THEN
       READ(31,7000)ts
7000   FORMAT(a80)
       I=2
      ENDIF
      L=ts(I:I)
      RETURN
      END

      function torsion(x0,x1,x2,x3,y0,y1,y2,y3,z0,z1,z2,z3,ad,bangle)
c     torsion, distance 0 1:ad , band angle:bangle
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      DIMENSION A(3),P(3),E0(3),XX(3),YY(3)
      PI=3.141592653589D0
      A(1)=X2-X3
      A(2)=Y2-Y3
      A(3)=Z2-Z3
      P(1)=X0-X1
      P(2)=Y0-Y1
      P(3)=Z0-Z1
      E0(1)=X1-X2
      E0(2)=Y1-Y2
      E0(3)=Z1-Z2
      SSE0E0=E0(1)*E0(1)+E0(2)*E0(2)+E0(3)*E0(3)
      DO 8 KK=1,3
8     E0(KK)=E0(KK)/DSQRT(SSE0E0)
      AD=DSQRT(P(1)*P(1)+P(2)*P(2)+P(3)*P(3))
      SSAE0=A(1)*E0(1)+A(2)*E0(2)+A(3)*E0(3)
      SSPE0=P(1)*E0(1)+P(2)*E0(2)+P(3)*E0(3)
      DO 9 KK=1,3
      XX(KK)=A(KK)-E0(KK)*SSAE0
9     YY(KK)=P(KK)-E0(KK)*SSPE0
      SSXXXX=XX(1)*XX(1)+XX(2)*XX(2)+XX(3)*XX(3)
      SSYYYY=YY(1)*YY(1)+YY(2)*YY(2)+YY(3)*YY(3)
      SSXXYY=XX(1)*YY(1)+XX(2)*YY(2)+XX(3)*YY(3)
      APX=XX(2)*YY(3)-XX(3)*YY(2)
      APY=XX(3)*YY(1)-XX(1)*YY(3)
      APZ=XX(1)*YY(2)-XX(2)*YY(1)
      SSE0AP=APX*E0(1)+APY*E0(2)+APZ*E0(3)
      SIGN=1.0D0
      IF (SSE0AP.NE.0.0D0)SIGN=-SSE0AP/ABS(SSE0AP)
      CB=0.0D0
      IF (AD.NE.0.0D0)CB=-SSPE0/AD
      CA=0.0D0
      IF (SSXXXX.NE.0.0D0.AND.SSYYYY.NE.0.0D0)
     1CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY)
      BANGLE=0.0D0
      ANGLE=0.0D0
      IF (CB.EQ.-1.0D0)BANGLE=PI
      IF (ABS(CB).LT.1.0D0)BANGLE=ACOS(CB)
      IF (CA.LE.-1.0D0) ANGLE=PI
      IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA)
      ANGLE=180.0D0* ANGLE/PI
      BANGLE=180.0D0*BANGLE/PI
      torsion=ANGLE*SIGN
      return 
      end