PROGRAM ALLIGN
      implicit none
      INTEGER*4 NTA,MENDELEV,MXB,i,ia,ix,iy,iz,ii,ity,itt,NAT,
     1igeo,isub,nd,idumm,ierr,ig,iter,j,MB,it,maxiter,itrial,
     2IX3
      PARAMETER (NTA=1100,MENDELEV=89,MXB=25,maxiter=1000)
      real*8 r,pi,x,y,z,ax,r1,ro,c,d,u,A,XS,XB,RTOL,dd,am,
     1oxx,amxx,ck,sk,dr1,dr2,fi,fk,pp,rat,wx,yy,z0,
     2tol
      real*4 amas
      DIMENSION X(NTA),ITY(NTA),r(NTA,3),r1(NTA,3),ro(NTA,3),
     1c(3),d(NTA,3),u(3),itt(NTA)
      CHARACTER*80 s80
      CHARACTER*1 OK,fp
      dimension amas(MENDELEV)
      data amas/1.007825,4.003,
     2  6.941, 9.012, 10.810,12.0  ,14.003074,15.994915,18.9984,20.179,
     3 22.990,24.305, 26.981,28.086,30.9737,31.972,34.968852,39.948,
     4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847,
     4 58.933,58.700,63.546,65.380,
     4                  69.720,72.590,74.922,78.960,79.904,83.800,
     5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070,
     5 102.906,106.400,107.868,112.410,
     5                 114.82,118.69,121.75,127.600,126.905,131.300,
     6 132.905,137.330,138.906,
     6                 140.120,140.908,144.240,145.000,150.400,
     6 151.960,157.250,158.925,162.500,164.930,167.260,168.934,
     6 173.040,174.970,
     6 178.490,180.948,183.850,186.207,190.200,192.220,195.090,
     6 196.967,207.590,204.370,207.200,208.980,210.000,210.001,
     6 222.02,
     7 223.000,226.025,227.028/     
      data z0/0.0d0/
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,MB,ITER
C     
      pi=4.0d0*atan(1.0d0)
      tol=1.0d-7

      WRITE(6,600)
600   format(/,' Alligns geometries in   FILE.X',/,
     1         ' so that mass center and rotation do not change',/)
c
      OPEN(4,FILE='FILE.X',STATUS='OLD')
      AM=z0
      READ(4,*)
      READ(4,*)NAT
      if(NAT.gt.NTA)then
       write(6,*)'too many atoms'
       close(4)
       stop
      endif
      DO 7 I=1,NAT
      READ(4,*)ITY(I),(r(I,ix),ix=1,3)
      do 71 ix=1,3
71    r1(I,ix)=r(I,ix)
      AM=AM+dble(AMAS(ITY(I)))
7     X(I)=dble(AMAS(ITY(I)))
      CLOSE(4)

      igeo=0
      open(4,file='FILE.X')
22    read(4,*,end=222,err=222)
      read(4,*,end=222,err=222)
      do 21 I=1,NAT
21    read(4,*,end=222,err=222)
      igeo=igeo+1
      goto 22
222   close(4)
c
2     WRITE(6,609)NAT,igeo,AM
609   format(i4,' ATOMS, ',i4,' GEOMETRIES, MOLECULAR MASS ',f10.3,/,/,
     1' HOW MANY ISOTOPIC SUBSTITUTION (negative for auto-D)? ',$)
      READ(5,*)ISUB
      write(6,*)
      if(ISUB.lt.0)then
       ISUB=0
       ND=0
c
c      Substitute all acidic Hs by Ds
       do 8 i=1,NAT
       ax=r(i,1)
       y=r(i,2)
       z=r(i,3)
       if(ITY(i).eq.1)then
        do 3 j=1,NAT
        it=ITY(j)
        dd=sqrt((ax-r(j,1))**2+(y-r(j,2))**2+(z-r(j,3))**2)
        if(it.eq.7.or.it.eq.8.or.it.eq.9.or.it.eq.16.or.it.eq.17.or.
     1  it.eq.35.or.it.eq.53)then
         if(dd.lt.1.23)then
          ND=ND+1
          X(i)=2.014000d0
         endif
        endif
3       continue
       endif
8      continue
       write(6,7009)ND
7009   format(i5,' acidic hydrogens have been substituted by deuteria')
      endif
      DO 5 II=1,ISUB
      WRITE(*,*)' H     1.007825  D    2.014000  T     3.016050'
      WRITE(*,*)' C12  12.000000  C13 13.003355  N14  14.003074'
      WRITE(*,*)' N15  15.000108  O16 15.994915  O17  16.999131'
      WRITE(*,*)' O18  17.999160  F   18.998403  CL35 34.968852'
      WRITE(*,*)' CL37 36.965903  P   30.973762  S    31.972070'
      WRITE(*,*)' GIVE NUMBER OF THE ATOM AND THE NEW MASS'
      READ(*,*)I,X(I)
5     CONTINUE
      DO 6 I=1,NAT,3
      IF(NAT-I.GE.2)WRITE(*,6676)I,X(I),I+1,X(I+1),I+2,X(I+2)
      IF(NAT-I.EQ.1)WRITE(*,6676)I,X(I),I+1,X(I+1)
      IF(NAT-I.EQ.0)WRITE(*,6676)I,X(I)
6     CONTINUE
6676  FORMAT(3(5X,I3,F12.6))
      AM=z0
      do 6001 ia=1,nat
6001  AM=AM+X(ia)
      WRITE(*,6002)AM
6002  format(/,' Mass ',f10.3,' CORRECT (Y/N) ?',$)
      READ(*,'(A)')OK
      IF ((OK.EQ.'n').OR.(OK.EQ.'N')) GOTO 2

      write(6,6003)
6003  format(/,' Reference First, Previous or fIt geometry (F,P,I)? ',$)
      read(5,'(a)')fp
      if(fp.eq.'f')fp='F'
      if(fp.eq.'i')fp='I'
      if(fp.eq.'I')then
       write(6,6006)
6006   format(/,' How many atoms? ',$)
       read(5,*)MB
       if(MB.gt.MXB)then
        write(6,*)' Too many atoms.'
	stop
       endif
       write(6,6005)MB
6005   format(/,' List the ',i2,' atoms: ',$)
       read(5,*)(ITT(ii),ii=1,MB)
      endif

      open(4,file='FILE.X')
      open(41,file='FILEC.X')
      do 204 ig=1,igeo
      read(4,404)s80
      write(41,404)s80
      read(4,404)s80
      write(41,404)s80
404   format(a80)
      do 201 ix=1,3
201   c(ix)=z0
      AM=z0
      do 200 ia=1,NAT
      read(4,*)idumm,(r(ia,ix),ix=1,3)
      AM=AM+X(ia)
      do 200 ix=1,3
200   c(ix)=c(ix)+r(ia,ix)*X(ia)
      do 202 ix=1,3
202   c(ix)=c(ix)/AM
      do 203 ia=1,NAT
      do 203 ix=1,3
203   r(ia,ix)=r(ia,ix)-c(ix)

      if(ig.gt.1)then
       if(fp.eq.'F'.or.fp.eq.'I')then
        do 207 ia=1,NAT
        do 207 ix=1,3
207     ro(ia,ix)=r1(ia,ix)
       endif
       if(fp.eq.'I')then
        do 212 ii=1,MB
	do 212 ix=1,3
	XB(ii,ix)=ro(ITT(ii),ix)
212     XS(ii,ix)=r(ITT(ii),ix)
        call DOMATRIX(IERR)
        write(6,6004)A
6004    format(' Transformation matrix:',/,3(3F10.4,/))
        do 213 ia=1,NAT
	ax=r(ia,1)
	y=r(ia,2)
	z=r(ia,3)
	r(ia,1)=A(1,1)*ax+A(1,2)*y+A(1,3)*z
	r(ia,2)=A(2,1)*ax+A(2,2)*y+A(2,3)*z
213     r(ia,3)=A(3,1)*ax+A(3,2)*y+A(3,3)*z
       endif
c
       itrial=1
88881  iter=0
8888   iter=iter+1
       if(iter.gt.maxiter)then
        write(6,*)'Maximum number of iterations exceeded'
	itrial=itrial+1
	if(itrial.lt.4)then
	 write(6,*)'Trial ',itrial
	 goto 88881
	endif
	write(6,*)'Iteration did not converge'
	goto 881
       endif

       do 210 ia=1,NAT
       do 210 ix=1,3
210    d(ia,ix)=r(ia,ix)-ro(ia,ix)
c
c      loop over the three axes,(try three different orderings):

       do 208 IX3=1,3
       IX=itrial+IX3-1
       if(IX.eq.4)IX=1
       if(IX.eq.5)IX=2
       IY=IX+1
       if(IY.eq.4)IY=1
       IZ=IY+1
       if(IZ.eq.4)IZ=1
       oxx=z0
       amxx=z0
       do 209 ia=1,NAT
       dr1=sqrt(ro(ia,IY)**2+ro(ia,IZ)**2)
       dr2=sqrt( r(ia,IY)**2+ r(ia,IZ)**2)
       pp=ro(ia,IY)*r(ia,IY)+ro(ia,IZ)*r(ia,IZ)
       if(abs(dr1*dr2).gt.1.0e-8)then
        rat=pp/dr1/dr2
	if(rat.gt. 1.0d0-1.0d-13)then
	 fi=z0
	else
	 if(rat.lt.-1.0d0+1.0d-13)then
	  fi=pi
	 else
          fi=acos(rat)
	 endif
	endif
	wx=ro(ia,IY)*d(ia,IZ)-ro(ia,IZ)*d(ia,IY)
	if(wx.lt.0.0)fi=-fi
c       amxx=amxx+X(ia)*fi*(5.0d0*dr1**2+2.0d0*dr2**2-dr1*dr2)/6.0d0
        amxx=amxx+X(ia)*fi*(dr1**2+dr2**2+dr1*dr2)/3.0d0
        oxx=oxx+X(ia)*dr2**2
       endif
209    continue

       if(abs(oxx).gt.1.0e-8)then
        fk=amxx/oxx
       else
        fk=z0
       endif
       u(IX)=fk
c
c      rotate back by this angle:
       ck=cos(-fk)
       sk=sin(-fk)
       do 208 ia=1,NAT
       YY=r(ia,IY)
       r(ia,IY)=YY*ck-r(ia,IZ)*sk
208    r(ia,IZ)=YY*sk+r(ia,IZ)*ck

c      write(6,*)iter,u(1)**2+u(2)**2+u(3)**2
       if(u(1)**2+u(2)**2+u(3)**2.gt.tol)goto 8888

      endif

881   write(6,608)ig,c,iter
608   format(i4,' center: ',3f10.6,', ',i7,' iterations')
      write(6,*)

      do 206 ia=1,NAT
      do 206 ix=1,3
206   ro(ia,ix)=r(ia,ix)

      if(ig.eq.1)then
       do 205 ia=1,NAT
       do 205 ix=1,3
205    r1(ia,ix)=r(ia,ix)
      endif

      do 204 ia=1,NAT
204   write(41,4141)ITY(ia),(r(ia,ix),ix=1,3)
4141  format(i4,3f12.6)

      close(4)
      close(41)
      write(6,*)'FILEC.X written'
      stop
      END

      
  

      SUBROUTINE DOMATRIX(IERR)
C     This is a Fortran version of Petr Malon's subroutine Amoeba
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MXB=25)
      DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3)
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER
      IERR=0
C     STARTING VALUES FOR THE ITERATION:
      ANG(1)=0.2d0
      ANG(2)=0.0d0
      ANG(3)=0.0d0
      DO 1 I=1,3
1     P(1,I)=ANG(I)
      Y(1)=FU(ANG)
      ANG(2)=ANG(1)
      ANG(1)=0.0d0
      DO 2 I=1,3
2     P(2,I)=ANG(I)
      Y(2)=FU(ANG)
      ANG(3)=ANG(2)
      ANG(2)=0.0d0
      DO 3 I=1,3
3     P(3,I)=ANG(I)
      Y(3)=FU(ANG)
      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(ANG)
      FTOL=1.0d-7
      ITMAX=1500
      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-11)THEN
       RTOL=0.0d0
      ELSE
       RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0
      ENDIF
      IF(RTOL.LT.FTOL)RETURN
      IF(ITER.EQ.ITMAX)THEN
       WRITE(3,*)' Rotation has not converged !'
       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(PR)
      IF(YPR.LE.Y(ILO))THEN
       DO 9 J=1,3
9      PRR(J)=2.0d0*PR(J)-PBAR(J)
       YPRR=FU(PRR)
       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(PRR)
        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(PR)
         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(ANG)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MXB=25)
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER
      DIMENSION ANG(3)
      S1=SIN(ANG(1))
      S2=SIN(ANG(2))
      S3=SIN(ANG(3))
      C1=COS(ANG(1))
      C2=COS(ANG(2))
      C3=COS(ANG(3))
C     ang1 = theta
C     ang2 = phi
C     ang3 = 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
      R=0.0d0
      DO 1 I=1,NAT
      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     R=R+(TX-XB(I,1))**2+(TY-XB(I,2))**2+(TZ-XB(I,3))**2
      FU=R
      RETURN
      END