CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C                  PROTROT - (c) Sean Parkin 96.01.12                  C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  A program to take two protein models with paired atoms and find     C
C  the optimum translation and rotations necessary to give the best    C
C  superposition of user selected atom types.  The r.m.s. deviation    C
C  for those atoms after transformation is also given.  Output         C
C  consists of a file with rotation and superposition information      C
C  'file2_rot.out', a file 'file2_rot.pdb' with the 2nd model rotated  C
C  onto the 1st, and a file 'file1-file2_ave.pdb'.  Other stuff,       C
C  notably the info. in 'file2_rot.out' is written to the screen.      C
C  Modified input files with just the specified atoms are also made.   C
C                                                                      C
C  The program is based on the quaternion method outlined in an        C
C  article by S.K. Kearsley (Acta Crystallographica (1989) A45,        C
C  208-210.).  The code for the Jacobi rotation in subroutine JACOBI   C
C  was adapted from a routine given in the book Numerical Recipes in   C
C  FORTRAN by Press, Teukolsky, Vetterling and Flannery.               C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  History of changes:                                                 C
C                                                                      C
C  Original version ---------------------------------------- 96.01.12  C
C  Various bugs fixed -------------------------------------- 96.10.24  C
C  Notes added within the program -------------------------- 96.10.25  C
C  Small change made to mode of information input ---------- 97.02.05  C
C  Small change made to mode of information output --------- 97.02.07  C
C  Output nicely formatted --------------------------------- 97.02.24  C
C  Incorporated PDBStrip program as subroutine PDBSTR ------ 97.09.17  C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This program can be used free of charge by anyone PROVIDED they     C
C  make adequate reference to it.  The most suitable reference is:     C
C                                                                      C
C  Parkin, S. (1996) "PROTROT - A Fortran program for superimposing    C
C                     two protein structures in pdb format."           C
C                     Lawrence Livermore National Laboratory           C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   Although it works fine on all systems tested the author makes no   C
C    guarantee concerning efficacy, efficiency or freedom from bugs    C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      PROGRAM PROTROT

      CHARACTER*80 MODEL1, MODEL2, NEWMOD
      INTEGER FLAG, N1, NUMROT
      DOUBLE PRECISION QUTRNN, EIGVEC, EIGVAL, ROTMAT, TRNVEC, THETAX,
     &   THETAY, THETAZ
      DIMENSION QUTRNN(4,4), EIGVEC(4,4), EIGVAL(4), ROTMAT(3,3),
     &   TRNVEC(6)

C  Call a subroutine to print out some identifying information.

      CALL BANNER

C  Ask user for two files in Brookhaven PDB file format, and determine
C  whether the fit is to be made over 1) All atoms, 2) Backbone atoms
C  only, 3) C(alpha) carbons only, 4) Side chain atoms only or 5) Back
C  bone atoms plus C beta carbons.

      CALL GETNFO (MODEL1, MODEL2, FLAG)

C  Create new files with names based on the two input model file names.
C  The names are amended so that 'name.pdb' becomes 'name_XX.pdb',
C  where XX is either AA (all atoms), CA (C alpha only), BB (backbone
C  only), PA (polyalanine chain) or SC (sidechains only).

      CALL PDBSTR (MODEL1, MODEL2, FLAG)

C  Check that the two PDB files have the same number of atoms.  As
C  the program progresses, more checking criteria will be added

      CALL MODCHK (MODEL1, MODEL2, N1)

C  Find the pair of vectors needed to translate the centroids of 
C  both models to the origin (0,0,0).  The vector for model 1 is 
C  in TRNVEC(1-3) and that of model 2 is in TRNVEC(4-6).

      CALL MODTRN (MODEL1, MODEL2, N1, TRNVEC)

C  Then build up the quaternion according to the method given in:
C  Kearsley, S.K. (1989) Acta Cryst. A45, 208-210.

      CALL QUTBLD (MODEL1, MODEL2, QUTRNN, N1, TRNVEC)

C  Diagonalize the symmetric matrix to produce four orthogonal unit
C  quaternions.

      CALL JACOBI (QUTRNN, EIGVEC, EIGVAL, NUMROT)

C  The rotation matrix is built according to the method given in:
C  Kearsley, S.K. (1989) Acta Cryst. A45, 208-210.

      CALL ROTBLD (EIGVEC, ROTMAT)

C  Figure out how much the rotation matrix corresponds to rotations
C  about the x,y and z axes.  This is done by calculating the angle
C  between the projection of vector (1,1,1) on the yz, xz and xy
C  planes before and after rotation.

      CALL ROTXYZ (ROTMAT, THETAX, THETAY, THETAZ)

C  Model 2 is then shifted to the origin, rotated according to the
C  rotation matrix calculated in ROTBLD, and then shifted to sit on
C  top of model 1.  A new file with '.rot' appended to the name of
C  model 2 is written.

      CALL MODROT (MODEL1, MODEL2, NEWMOD, N1, ROTMAT, TRNVEC)

C  Write a file that is the average of the target (model 1) and the
C  rotated-translated probe (model 2).

      CALL MODAVE (MODEL1, MODEL2, NEWMOD, N1)

C  A file with '_rot.out' appended to the stem name of model 2 is
C  written.  It contains contains some statistics and the various
C  eigenvectors, eigenvalues, rotation matrix etc. generated in the 
C  run.

      CALL OUTPUT (MODEL1, MODEL2, N1, EIGVAL, EIGVEC, ROTMAT, TRNVEC,
     &    THETAX, THETAY, THETAZ)

      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE BANNER

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This simple subroutine prints out some identifying information to   C
C  the screen.                                                         C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      WRITE (*,'(1X,A78)') '+------------------------------------------'
     &   //'----------------------------------+'
      WRITE (*,'(A2,76X,A1)') '|', '|'
      WRITE (*,'(1X,A78)') '|                 Program PROTROT - (c) Sea'
     &   //'n Parkin 96.01.12                 |'
      WRITE (*,'(1X,A78)') '|      A program to calculate the optimal s'
     &   //'uperposition for two protein      |'
      WRITE (*,'(1X,A78)') '|                  structure files in Brook'
     &   //'haven PDB format                  |'
      WRITE (*,'(A2,76X,A1)') '|', '|'
      WRITE (*,'(1X,A78)') '+------------------------------------------'
     &   //'----------------------------------+'
      WRITE (*,*)

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE GETNFO (MODEL1, MODEL2, FLAG)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine asks for the name of two files in Brookhaven PDB    C
C  file format.  It also asks whether the fit is to be made based on   C
C  all atoms, backbone atoms only, C alpha carbons, side chain only    C
C  or a polyalanine chain.                                             C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 MODEL1, MODEL2
      INTEGER FLAG

100   FORMAT (A80)
200   FORMAT (I1)

      FLAG = 0

      WRITE (*,'(1X,A$)') 'Input filename for model 1: '
      READ (*,100) MODEL1
      WRITE (*,'(1X,A$)') 'Input filename for model 2: '
      READ (*,100) MODEL2

      MODEL1 = MODEL1(1:INDEX(MODEL1,' ')-1)
      MODEL2 = MODEL2(1:INDEX(MODEL2,' ')-1)

      WRITE (*,*)  
      WRITE (*,*) 'Which atoms would you like to base the comparison on'
     &   //' ?'
      WRITE (*,*)
      WRITE (*,*) '  All atoms (default): [1] '
      WRITE (*,*) '  Alpha carbons only:  [2] '
      WRITE (*,*) '  Backbone atoms only: [3] '
      WRITE (*,*) '  Polyalanine chain:   [4] '
      WRITE (*,*) '  Side chains only:    [5] '
      WRITE (*,*)

      READ 200, FLAG
      IF ((FLAG.LT.1).OR.(FLAG.GT.5)) THEN
         WRITE (*,*) 'Crazy input ignored, fitting all atoms ... '
      ENDIF  

c      WRITE (*,*) 'The fit will be made based on all atoms.  If you wan'
c     &   //'t to do the fit on ' 
c      WRITE (*,*) 'some subset of atoms, the accompanying program PDBST'
c     &   //'RIP can be used to '
c      WRITE (*,*) 'prepare input files. '

      WRITE (*,*)
      WRITE (*,*) 'Running ... '
      WRITE (*,*)

      RETURN
      END

C=======================================================================


C-----------------------------------------------------------------------

      SUBROUTINE PDBSTR (MODEL1, MODEL2, FLAG)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine takes each of the input protein model files and     C
C  writes out new files containing just the specified atoms.  The new  C
C  files have either AA, CA, BB, PA or SC added to their names so as   C
C  to distinguish them from the original files and also to indicate    C
C  their contents.                                                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 MODEL1, MODEL2, LINE
      CHARACTER*4 WORD, ATOM
      CHARACTER*7 ATOMS
      INTEGER FLAG, CTRL

 100  FORMAT (A80)

      OPEN (10, FILE = MODEL1, STATUS = 'OLD', ERR = 9990)
      OPEN (11, FILE = MODEL2, STATUS = 'OLD', ERR = 9991)

      IF (FLAG.EQ.1) ATOMS = '_AA.pdb'
      IF (FLAG.EQ.2) ATOMS = '_CA.pdb'
      IF (FLAG.EQ.3) ATOMS = '_BB.pdb'
      IF (FLAG.EQ.4) ATOMS = '_PA.pdb'
      IF (FLAG.EQ.5) ATOMS = '_SC.pdb'

      MODEL1 = MODEL1(1:INDEX(MODEL1,' ')-5)//ATOMS
      MODEL2 = MODEL2(1:INDEX(MODEL2,' ')-5)//ATOMS

      OPEN (12, FILE = MODEL1, STATUS = 'UNKNOWN')
      OPEN (13, FILE = MODEL2, STATUS = 'UNKNOWN')

      DO 1001 N = 10, 11
1        READ (N,100, ERR = 1001), LINE
         ATOM = LINE(1:4)
         WORD = LINE(14:17)
 
         IF (FLAG.EQ.1) THEN
            IF (ATOM.EQ.'ATOM') THEN
               WRITE (N+2,100) LINE
            ENDIF
         ELSE IF (FLAG.EQ.2) THEN
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'CA')) THEN
               WRITE (N+2,100) LINE
            ENDIF
         ELSE IF (FLAG.EQ.3) THEN
            CTRL = 0
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'N ')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'CA')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'C ')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'O ')) CTRL = 1
            IF (CTRL.EQ.1) THEN
               WRITE (N+2,100) LINE
            ENDIF
         ELSE IF (FLAG.EQ.4) THEN
            CTRL = 0
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'N ')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'CA')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'C ')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'O ')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'CB')) CTRL = 1
            IF (CTRL.EQ.1) THEN
               WRITE (N+2,100) LINE
            ENDIF
         ELSE IF (FLAG.EQ.5) THEN
            CTRL = 0
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'N ')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'CA')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'C ')) CTRL = 1
            IF ((ATOM.EQ.'ATOM').AND.(WORD(1:2).EQ.'O ')) CTRL = 1
            IF (CTRL.EQ.0) THEN
               WRITE (N+2,100) LINE
            ENDIF
         ENDIF
         GOTO 1
1001  CONTINUE

      CLOSE (10)
      CLOSE (11)
      CLOSE (12)
      CLOSE (13)

      RETURN

9990  WRITE (*,*) 'There is no such file as ', MODEL1
      STOP
9991  WRITE (*,*) 'There is no such file as ', MODEL2
      STOP

      END      
C=======================================================================


C-----------------------------------------------------------------------

      SUBROUTINE MODCHK (MODEL1, MODEL2, N1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine checks the two input files to make sure they        C
C  correspond.  Currently, the only check done is to ensure that       C
C  the number of atoms in each is the same.  If it is not, then        C
C  the execution of MODCHK and hence ProtRot is stopped.               C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 MODEL1, MODEL2
      CHARACTER*6 WORD
      INTEGER N1, N2

100   FORMAT (A6)

      N1 = 0
      N2 = 0

      OPEN (10, FILE = MODEL1, STATUS = 'OLD', ERR = 9990)
      OPEN (11, FILE = MODEL2, STATUS = 'OLD', ERR = 9991)

1     READ (10,100, ERR = 2) WORD
      IF ((WORD.EQ.'ATOM  ').OR.(WORD.EQ.'HETATM')) N1 = N1 + 1
      GOTO 1

2     READ (11,100, ERR = 3) WORD
      IF ((WORD.EQ.'ATOM  ').OR.(WORD.EQ.'HETATM')) N2 = N2 + 1
      GOTO 2

3     CONTINUE
      IF (N1.NE.N2) THEN 
         WRITE (*,*) 'The models have different numbers of atoms '
         CLOSE (10)
         CLOSE (11)
         STOP
      ENDIF

      CLOSE (10)
      CLOSE (11)

      RETURN

9990  WRITE (*,*) 'There is no such file as ', MODEL1
      STOP
9991  WRITE (*,*) 'There is no such file as ', MODEL2
      STOP

      END      

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE MODTRN (MODEL1, MODEL2, N1, TRNVEC)
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine obtains the centroid of both MODEL1 and MODEL2.     C
C  The difference is needed to superimpose the second model on the     C
C  first so that the subsequent rotation is valid.                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 MODEL1, MODEL2
      CHARACTER*66 LINEA, LINEB
      CHARACTER*6 WORD1, WORD2
      INTEGER N1, N
      DOUBLE PRECISION X1, Y1, Z1, X2, Y2, Z2, TRNVEC, X1SUM, Y1SUM, 
     &   Z1SUM, X2SUM, Y2SUM, Z2SUM
      DIMENSION TRNVEC(6)

100   FORMAT (30X, 3F8.3)
200   FORMAT (A66)

      X1SUM = 0.0
      Y1SUM = 0.0
      Z1SUM = 0.0
      X2SUM = 0.0
      Y2SUM = 0.0
      Z2SUM = 0.0

      OPEN (10, FILE = MODEL1, STATUS = 'OLD')
      OPEN (11, FILE = MODEL2, STATUS = 'OLD')

C  Sum up all the x, y and z coordinates of both models.  These 
C  numbers can get pretty big, hence they are double precision.
C  The strange internal reads are used to skip lines that are
C  not ATOM or HETATM records.

      DO 1 N = 1, N1

2        READ (10,200) LINEA
         WORD1 = LINEA(1:6)
         IF ((WORD1.EQ.'ATOM  ').OR.(WORD1.EQ.'HETATM')) THEN
            READ (LINEA,100) X1, Y1, Z1
            GOTO 3
         ELSE
            GOTO 2
         ENDIF
3        READ (11,200) LINEB
         WORD2 = LINEB(1:6)
         IF ((WORD2.EQ.'ATOM  ').OR.(WORD2.EQ.'HETATM')) THEN
            READ (LINEB,100) X2, Y2, Z2
            GOTO 4
         ELSE
            GOTO 3
         ENDIF
4        CONTINUE

         X1SUM = X1SUM + X1
         Y1SUM = Y1SUM + Y1
         Z1SUM = Z1SUM + Z1
         X2SUM = X2SUM + X2
         Y2SUM = Y2SUM + Y2
         Z2SUM = Z2SUM + Z2
1     CONTINUE

C  Divide these sums by the number of atoms to give the displacement
C  from the origin of each model.

      TRNVEC(1) = X1SUM / N1
      TRNVEC(2) = Y1SUM / N1
      TRNVEC(3) = Z1SUM / N1
      TRNVEC(4) = X2SUM / N1
      TRNVEC(5) = Y2SUM / N1
      TRNVEC(6) = Z2SUM / N1

      CLOSE (10)
      CLOSE (11)

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE QUTBLD (MODEL1, MODEL2, QUTRNN, N1, TRNVEC)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine builds the quaternion from coordinates of the       C
C  paired atoms of models 1 and 2.                                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 MODEL1, MODEL2
      CHARACTER*66 LINEA, LINEB
      CHARACTER*6 WORD1, WORD2
      INTEGER I, J, N, N1
      DOUBLE PRECISION QUTRNN, X1, Y1, Z1, X2, Y2, Z2, XM, YM, ZM, XP, 
     &   YP, ZP, TRNVEC
      DIMENSION QUTRNN(4,4), TRNVEC(6)

100   FORMAT (30X, 3F8.3)
200   FORMAT (A66)

C  Initialize the quaternion matrix elements to zero.

      DO 1 I = 1, 4
         DO 2 J = I, 4
            QUTRNN(I,J) = 0.0
            QUTRNN(J,I) = 0.0
2        CONTINUE
         QUTRNN(I,I) = 0.0
1     CONTINUE

      OPEN (10, FILE = MODEL1, STATUS = 'OLD')
      OPEN (11, FILE = MODEL2, STATUS = 'OLD')

C  Build the (4x4) symmetric matrix that gives the four unit
C  quaternions after diagonalization.  The internal READs
C  ensure that only ATOM or HETATM records are read.

      DO 3 N = 1, N1

4        READ (10,200) LINEA
         WORD1 = LINEA(1:6)
         IF ((WORD1.EQ.'ATOM  ').OR.(WORD1.EQ.'HETATM')) THEN
            READ (LINEA,100) X1, Y1, Z1
            GOTO 5
         ELSE
            GOTO 4
         ENDIF
5        READ (11,200) LINEB
         WORD2 = LINEB(1:6)
         IF ((WORD2.EQ.'ATOM  ').OR.(WORD2.EQ.'HETATM')) THEN
            READ (LINEB,100) X2, Y2, Z2
            GOTO 6
         ELSE
            GOTO 5
         ENDIF
6        CONTINUE

         XM = (X1 - TRNVEC(1)) - (X2 - TRNVEC(4))
         XP = (X1 - TRNVEC(1)) + (X2 - TRNVEC(4))
         YM = (Y1 - TRNVEC(2)) - (Y2 - TRNVEC(5))
         YP = (Y1 - TRNVEC(2)) + (Y2 - TRNVEC(5))
         ZM = (Z1 - TRNVEC(3)) - (Z2 - TRNVEC(6))
         ZP = (Z1 - TRNVEC(3)) + (Z2 - TRNVEC(6))
         QUTRNN(1,1) = QUTRNN(1,1) + (XM**2 + YM**2 + ZM**2)
         QUTRNN(2,2) = QUTRNN(2,2) + (YP**2 + ZP**2 + XM**2)
         QUTRNN(3,3) = QUTRNN(3,3) + (XP**2 + ZP**2 + YM**2)
         QUTRNN(4,4) = QUTRNN(4,4) + (XP**2 + YP**2 + ZM**2)
         QUTRNN(1,2) = QUTRNN(1,2) + (YP*ZM - YM*ZP)
         QUTRNN(1,3) = QUTRNN(1,3) + (XM*ZP - XP*ZM)
         QUTRNN(1,4) = QUTRNN(1,4) + (XP*YM - XM*YP)
         QUTRNN(2,3) = QUTRNN(2,3) + (XM*YM - XP*YP)
         QUTRNN(2,4) = QUTRNN(2,4) + (XM*ZM - XP*ZP)
         QUTRNN(3,4) = QUTRNN(3,4) + (YM*ZM - YP*ZP)
         QUTRNN(2,1) = QUTRNN(1,2)
         QUTRNN(3,1) = QUTRNN(1,3)
         QUTRNN(4,1) = QUTRNN(1,4)
         QUTRNN(3,2) = QUTRNN(2,3)
         QUTRNN(4,2) = QUTRNN(2,4)
         QUTRNN(4,3) = QUTRNN(3,4)
3     CONTINUE

      CLOSE (10)
      CLOSE (11)

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE JACOBI (QUTRNN, EIGVEC, EIGVAL, NUMROT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine performs the Jacobi rotation and sorts the          C
C  eigenvectors and eigenvalues.  Modified from a routine given        C
C  in Numerical Recipes in Fortran 2nd Ed. by Press, Teukolsky,        C
C  Vetterling and Flannery.                                            C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER NUMROT, I, J, K, P, Q
      DOUBLE PRECISION QUTRNN, EIGVEC, EIGVAL, BLD, ACC, SINPHI, 
     &   COSPHI, ARP, ARQ, SUM, T, TAU, THETA, THRESH, SWAP
      DIMENSION QUTRNN(4,4), EIGVEC(4,4), EIGVAL(4), BLD(4), ACC(4)

C  Set NUMROT to zero, initialize EIGVEC to the identity matrix, 
C  initialize EIGVAL and BLD to the diagonal of QUTRNN, and fill 
C  ACC with zeros.

      NUMROT = 0
      DO 1 P = 1, 4
         DO 2 Q = 1, 4
            EIGVEC(P,Q) = 0.0
2        CONTINUE
         EIGVEC(P,P) = 1.0
         EIGVAL(P) = QUTRNN(P,P)
         BLD(P) = QUTRNN(P,P)
         ACC(P) = 0.0
1     CONTINUE

C  JACOBI proper starts here.

      DO 3 I = 1, 50
         SUM = 0.0
         DO 4 P = 1, 3
            DO 5 Q = P+1, 4
               SUM = SUM + ABS(QUTRNN(P,Q))
5           CONTINUE
4        CONTINUE
         IF (SUM.EQ.0.0) GO TO 1000
         IF (I.LT.4) THEN
            THRESH = 0.0125 * SUM
         ELSE
            THRESH = 0.0
         ENDIF
         DO 6 P = 1, 3
            DO 7 Q = P+1, 4
               ARP = 100.0 * ABS(QUTRNN(P,Q))
               IF ((I.GT.4).AND.(ABS(EIGVAL(P))+ARP.EQ.ABS(EIGVAL(P)))
     &            .AND.(ABS(EIGVAL(Q))+ARP.EQ.ABS(EIGVAL(Q)))) THEN
                  QUTRNN(P,Q) = 0.0
               ELSE IF (ABS(QUTRNN(P,Q)).GT.THRESH) THEN
                  ARQ = EIGVAL(Q) - EIGVAL(P)
                  IF (ABS(ARQ)+ARP.EQ.ABS(ARQ)) THEN
                     T = QUTRNN(P,Q) / ARQ
                  ELSE
                     THETA = 0.5 * ARQ / QUTRNN(P,Q)
                     T = 1.0 / (ABS(THETA) + SQRT(1.0 + THETA ** 2))
                     IF (THETA.LT.0.0) T = 0.0 - T
                  ENDIF
                  COSPHI = 1.0 / SQRT(1.0 + T ** 2)
                  SINPHI = T * COSPHI
                  TAU = SINPHI / (1.0 + COSPHI)
                  ARQ = T * QUTRNN(P,Q)
                  ACC(P) = ACC(P) - ARQ
                  ACC(Q) = ACC(Q) + ARQ
                  EIGVAL(P) = EIGVAL(P) - ARQ
                  EIGVAL(Q) = EIGVAL(Q) + ARQ
                  QUTRNN(P,Q) = 0.0
                  DO 8 J = 1, P-1
                     ARP = QUTRNN(J,P)
                     ARQ = QUTRNN(J,Q)
                     QUTRNN(J,P) = ARP - SINPHI * (ARQ + ARP * TAU)
                     QUTRNN(J,Q) = ARQ + SINPHI * (ARP - ARQ * TAU)
8                 CONTINUE
                  DO 9 J = P+1, Q-1
                     ARP = QUTRNN(P,J)
                     ARQ = QUTRNN(J,Q)
                     QUTRNN(P,J) = ARP - SINPHI * (ARQ + ARP * TAU)
                     QUTRNN(J,Q) = ARQ + SINPHI * (ARP - ARQ * TAU)
9                 CONTINUE
                  DO 10 J = Q+1, 4
                     ARP = QUTRNN(P,J)
                     ARQ = QUTRNN(Q,J)
                     QUTRNN(P,J) = ARP - SINPHI * (ARQ + ARP * TAU)
                     QUTRNN(Q,J) = ARQ + SINPHI * (ARP - ARQ * TAU)
10                CONTINUE
                  DO 11 J = 1, 4
                     ARP = EIGVEC(J,P)
                     ARQ = EIGVEC(J,Q)
                     EIGVEC(J,P) = ARP - SINPHI * (ARQ + ARP * TAU)
                     EIGVEC(J,Q) = ARQ + SINPHI * (ARP - ARQ * TAU)
11                CONTINUE
                  NUMROT = NUMROT + 1
               ENDIF
7           CONTINUE
6        CONTINUE
         DO 12 P = 1, 4
            BLD(P) = BLD(P) + ACC(P)
            EIGVAL(P) = BLD(P)
            ACC(P) = 0.0
12       CONTINUE
3     CONTINUE

1000  CONTINUE

C  JACOBI is done, now sort the eigenvalues and eigenvectors into
C  the order of descending eigenvalues.

      DO 13 I = 1, 3
         K = I
         SWAP = EIGVAL(I)
         DO 14 J = I + 1, 4
            IF (EIGVAL(J).GE.SWAP) THEN 
               K = J
               SWAP = EIGVAL(J)
            ENDIF
14       CONTINUE
         IF (K.NE.I) THEN
            EIGVAL(K) = EIGVAL(I)
            EIGVAL(I) = SWAP
            DO 15 J = 1, 4
               SWAP = EIGVEC(J,I)
               EIGVEC(J,I) = EIGVEC(J,K)
               EIGVEC(J,K) = SWAP
15          CONTINUE
         ENDIF
13    CONTINUE

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE ROTBLD (EIGVEC, ROTMAT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine constructs the rotation matrix that transforms      C
C  model 2 onto model 1 from the elements of the eigenvector with      C
C  the lowest corresponding eigenvalue.                                C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DOUBLE PRECISION EIGVEC, ROTMAT, Q1, Q2, Q3, Q4
      DIMENSION EIGVEC(4,4), ROTMAT(3,3)

C  The smallest eigenvalue corresponds to the last eigenvector (sorting
C  is done within subroutine JACOBI).  This eigenvector contains the  
C  four quaternion components used to build the rotation matrix.

      Q1 = EIGVEC(1,4)
      Q2 = EIGVEC(2,4)
      Q3 = EIGVEC(3,4)
      Q4 = EIGVEC(4,4)

C  Build the rotation matrix.

      ROTMAT(1,1) = Q1**2 + Q2**2 - Q3**2 - Q4**2
      ROTMAT(1,2) = 2 * (Q2*Q3 + Q1*Q4)
      ROTMAT(1,3) = 2 * (Q2*Q4 - Q1*Q3)
      ROTMAT(2,1) = 2 * (Q2*Q3 - Q1*Q4)
      ROTMAT(2,2) = Q1**2 + Q3**2 - Q2**2 - Q4**2
      ROTMAT(2,3) = 2 * (Q3*Q4 + Q1*Q2)
      ROTMAT(3,1) = 2 * (Q2*Q4 + Q1*Q3)
      ROTMAT(3,2) = 2 * (Q3*Q4 - Q1*Q2)
      ROTMAT(3,3) = Q1**2 + Q4**2 - Q2**2 - Q3**2

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE ROTXYZ (ROTMAT, THETAX, THETAY, THETAZ)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine figures out the components of the rotation matrix   C
C  about the x,y and z axes - angles of rotation in the yz, xz and xy  C
C  planes.                                                             C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DOUBLE PRECISION ROTMAT, THETAX, THETAY, THETAZ, SUM11, SUM12, 
     &   SUM13, SUM21, SUM22, SUM23, SUM31, SUM32, SUM33, SQRT2, SQRTX, 
     &   SQRTY, SQRTZ
      DIMENSION ROTMAT(3,3)

      SUM11 = ROTMAT(1,2) + ROTMAT(1,3)
      SUM12 = ROTMAT(2,2) + ROTMAT(2,3)
      SUM13 = ROTMAT(3,2) + ROTMAT(3,3)
      SUM21 = ROTMAT(1,1) + ROTMAT(1,3)
      SUM22 = ROTMAT(2,1) + ROTMAT(2,3)
      SUM23 = ROTMAT(3,1) + ROTMAT(3,3)
      SUM31 = ROTMAT(1,1) + ROTMAT(1,2)
      SUM32 = ROTMAT(2,1) + ROTMAT(2,2)
      SUM33 = ROTMAT(3,1) + ROTMAT(3,2)

      SQRT2 = SQRT(2.0)
      SQRTX = SQRT((SUM11 ** 2) + (SUM12 ** 2) + (SUM13 ** 2))
      SQRTY = SQRT((SUM21 ** 2) + (SUM22 ** 2) + (SUM23 ** 2))
      SQRTZ = SQRT((SUM31 ** 2) + (SUM32 ** 2) + (SUM33 ** 2))

      THETAX = ACOS((SUM12 + SUM13) / (SQRT2 * SQRTX))
      THETAY = ACOS((SUM21 + SUM23) / (SQRT2 * SQRTY))
      THETAZ = ACOS((SUM31 + SUM32) / (SQRT2 * SQRTZ))

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE MODROT (MODEL1, MODEL2, NEWMOD, N1, ROTMAT, TRNVEC)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine applies the translation and rotation found earlier  C
C  to model 2 so that it superimposes on model 1.  A new file, with    C
C  '.rot' appended to the name of model 2 is created.                  C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 MODEL1, MODEL2, NEWMOD
      CHARACTER*66 LINEA
      CHARACTER*30 WORD1
      CHARACTER*12 WORD2
      CHARACTER*6  WORDA
      INTEGER N, N1
      DOUBLE PRECISION ROTMAT, X2, Y2, Z2, XNEW, YNEW, ZNEW, TRNVEC
      DIMENSION ROTMAT(3,3), TRNVEC(6)

100   FORMAT (A30, 3F8.3, A12)
200   FORMAT (A)
300   FORMAT (A66)

      NEWMOD = MODEL2(1:INDEX(MODEL2,' ')-5)//'_rot.pdb'

      OPEN (11, FILE = MODEL2, STATUS = 'OLD')
      OPEN (13, FILE = NEWMOD, STATUS = 'UNKNOWN')

      WRITE (13,200) 'HEADER  '//NEWMOD
      WRITE (13,200) 'REMARK  '
      WRITE (13,200) 'REMARK  FILE WAS WRITTEN BY PROGRAM PROTROT AFTER'
     &   //' SUPERPOSITION'
      WRITE (13,200) 'REMARK  ACCORDING TO THE QUATERNION METHOD OF S.K'
     &   //'.KEARSLEY'
      WRITE (13,200) 'REMARK  (1989) ACTA CRYST. A45, 208-210.'
      WRITE (13,200) 'REMARK  FILE CONTAINS THE STRUCTURE AFTER SUPERPO'
     &   //'SITION OF'
      WRITE (13,200) 'REMARK  '//MODEL2(1:INDEX(MODEL2,' '))//'ONTO '
     &   //MODEL1(1:INDEX(MODEL1,' '))
      WRITE (13,200) 'REMARK  '

      DO 1 N = 1, N1

C  The strange internal read ensures that only ATOM and HETATM records
C  are read.

2        READ (11,300) LINEA
         WORDA = LINEA(1:6)
         IF ((WORDA.EQ.'ATOM  ').OR.(WORDA.EQ.'HETATM')) THEN
            READ (LINEA,100) WORD1, X2, Y2, Z2, WORD2
            GOTO 3
         ELSE
            GOTO 2
         ENDIF
3        CONTINUE 

C  Move model 2 so that its centroid is at the origin.

         X2 = X2 - TRNVEC(4)
         Y2 = Y2 - TRNVEC(5)
         Z2 = Z2 - TRNVEC(6)

C  Perform the rotation.

         XNEW = (ROTMAT(1,1)*X2) + (ROTMAT(1,2)*Y2) + (ROTMAT(1,3)*Z2)
         YNEW = (ROTMAT(2,1)*X2) + (ROTMAT(2,2)*Y2) + (ROTMAT(2,3)*Z2)
         ZNEW = (ROTMAT(3,1)*X2) + (ROTMAT(3,2)*Y2) + (ROTMAT(3,3)*Z2)

C  Move the centroid of the rotated model to the centroid of model 1.

         XNEW = XNEW + TRNVEC(1)
         YNEW = YNEW + TRNVEC(2)
         ZNEW = ZNEW + TRNVEC(3)

         WRITE (13,100) WORD1, XNEW, YNEW, ZNEW, WORD2

1     CONTINUE

      CLOSE (11)
      CLOSE (13)

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

      SUBROUTINE MODAVE (MODEL1, MODEL2, NEWMOD, N1)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine writes a file with the average of the target        C
C  (model 1) and the rotated-translated probe (model 2).  This is      C
C  useful as an improved model for investigations of missed symmetry   C
C  in models with two (etc.) molecules in the asymmetric unit.         C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 MODEL1, MODEL2, NEWMOD, AVEMOD, WORD1A, WORD1B,
     &   WORD2A, WORD2B
      CHARACTER*66 LINEA, LINEB
      INTEGER N, N1
      REAL X1, Y1, Z1, B1, X2TRAN, Y2TRAN, Z2TRAN, B2, XAVE, YAVE, ZAVE,
     &   BAVE

100   FORMAT (A30, 3F8.3, A6, F6.2)
200   FORMAT (A)
300   FORMAT (A66)

      AVEMOD = MODEL1(1:INDEX(MODEL1,' ')-5)//'-'//
     &         MODEL2(1:INDEX(MODEL2,' ')-5)//'_ave.pdb'

      OPEN (11, FILE = MODEL1, STATUS = 'OLD')
      OPEN (12, FILE = NEWMOD, STATUS = 'OLD')
      OPEN (13, FILE = AVEMOD, STATUS = 'UNKNOWN')

      WRITE (13,200) 'HEADER  '//AVEMOD
      WRITE (13,200) 'REMARK  '
      WRITE (13,200) 'REMARK  FILE WAS WRITTEN BY PROGRAM PROTROT AFTER'
     &   //' SUPERPOSITION'
      WRITE (13,200) 'REMARK  ACCORDING TO THE QUATERNION METHOD OF S.K'
     &   //'.KEARSLEY'
      WRITE (13,200) 'REMARK  (1989) ACTA CRYST. A45, 208-210.'
      WRITE (13,200) 'REMARK  FILE CONTAINS THE AVERAGE STRUCTURE AFTER'
     &   //' SUPERPOSITION OF'
      WRITE (13,200) 'REMARK  '//MODEL2(1:INDEX(MODEL2,' '))//'ONTO '
     &   //MODEL1(1:INDEX(MODEL1,' '))
      WRITE (13,200) 'REMARK  '

      DO 1 N = 1, N1

C The strange internal reads here are to enable the non ATOM lines in 
C the pdb files to be skipped.

2        READ (11,300) LINEA
         WORD1A = LINEA(1:6)
         IF ((WORD1A.EQ.'ATOM  ').OR.(WORD1A.EQ.'HETATM')) THEN
            READ (LINEA,100) WORD1A, X1, Y1, Z1, WORD1B, B1
            GOTO 3
         ELSE
            GOTO 2
         ENDIF

3        READ (12,300) LINEB
         WORD2A = LINEB(1:6)
         IF ((WORD2A.EQ.'ATOM  ').OR.(WORD2A.EQ.'HETATM')) THEN
            READ (LINEB,100) WORD2A, X2TRAN, Y2TRAN, Z2TRAN, WORD2B, B2
            XAVE = (X1 + X2TRAN) / 2.0
            YAVE = (Y1 + Y2TRAN) / 2.0
            ZAVE = (Z1 + Z2TRAN) / 2.0
            BAVE = (B1 + B2) / 2.0
            WRITE (13,100) WORD1A, XAVE, YAVE, ZAVE, WORD1B, BAVE
         ELSE
            GOTO 3
         ENDIF
1     CONTINUE

      CLOSE (11)
      CLOSE (12)
      CLOSE (13)

      RETURN
      END

C=======================================================================

C-----------------------------------------------------------------------

       SUBROUTINE OUTPUT (MODEL1, MODEL2, N1, EIGVAL, EIGVEC, ROTMAT,
     &   TRNVEC, THETAX, THETAY, THETAZ)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine prints some of the relevant run information to a    C
C  file *_rot.out and also to the screen.                              C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 MODEL1, MODEL2, OUTFIL, LINE
      INTEGER N1, N
      DOUBLE PRECISION EIGVAL, EIGVEC, ROTMAT, TRNVEC, THETAX, THETAY,
     &   THETAZ
      DIMENSION EIGVAL(4), EIGVEC(4,4), ROTMAT(3,3), TRNVEC(6)

      OUTFIL = MODEL2(1:INDEX(MODEL2,' ')-5)//'_rot'//'.out'

      OPEN (12, FILE = OUTFIL, STATUS = 'UNKNOWN')

      WRITE (12,'(1X,A78)') '+-----------------------------------------'
     &   //'-----------------------------------+'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(1X,A78)') '|           Output from program PROTROT - '
     &   //'(c) Sean Parkin 96.01.12           |'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(1X,A78)') '+========================================='
     &   //'===================================+'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(A2,30X,3(11X,A1),10X,A1)') '|', 'X', 'Y', 'Z', '|'
      WRITE (12,'(A2,6X,A27,3(F12.5),7X,A1))')
     &   '|','Centroid of model 1 is at: ', TRNVEC(1), TRNVEC(2), 
     &   TRNVEC(3), '|'
      WRITE (12,'(A2,6X,A27,3(F12.5),7X,A1))')
     &   '|','Centroid of model 2 is at: ', TRNVEC(4), TRNVEC(5), 
     &   TRNVEC(6), '|'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(1X,A78)') '+-----------------------------------------'
     &   //'-----------------------------------+'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(A2,5X,A66,5X,A1)') '|', 'Eigenvalues and their corres'
     &   //'ponding eigenvectors are (by columns) ','|'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(A2,3X,4(F16.6),9X,A1)') '|', EIGVAL, '|'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(A2,3X,4(F16.6),9X,A1)') '|', EIGVEC(1,1), EIGVEC(1,2),
     &   EIGVEC(1,3), EIGVEC(1,4), '|'
      WRITE (12,'(A2,3X,4(F16.6),9X,A1)') '|', EIGVEC(2,1), EIGVEC(2,2), 
     &   EIGVEC(2,3), EIGVEC(2,4), '|'
      WRITE (12,'(A2,3X,4(F16.6),9X,A1)') '|', EIGVEC(3,1), EIGVEC(3,2),
     &   EIGVEC(3,3), EIGVEC(3,4), '|'
      WRITE (12,'(A2,3X,4(F16.6),9X,A1)') '|', EIGVEC(4,1), EIGVEC(4,2),
     &   EIGVEC(4,3), EIGVEC(4,4), '|'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(1X,A78)') '+-----------------------------------------'
     &   //'-----------------------------------+'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(A2,14X,A47,15X,A1)') '|', 'The rotation matrix and tr'
     &   //'anslation vectors are ', '|'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(A2,31X,A1,31X,A1,12X,A1)') '|', 'R', 'T', '|'
      WRITE (12,'(A2,3X,4(F16.6),9X,A1)') '|', ROTMAT(1,1), ROTMAT(1,2), 
     &   ROTMAT(1,3), TRNVEC(4) - TRNVEC(1), '|'
      WRITE (12,'(A2,3X,4(F16.6),9X,A1)') '|', ROTMAT(2,1), ROTMAT(2,2), 
     &   ROTMAT(2,3), TRNVEC(5) - TRNVEC(2), '|'
      WRITE (12,'(A2,3X,4(F16.6),9X,A1)') '|', ROTMAT(3,1), ROTMAT(3,2), 
     &   ROTMAT(3,3), TRNVEC(6) - TRNVEC(3), '|'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(1X,A78)') '+-----------------------------------------'
     &   //'-----------------------------------+'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(A2,3X,A69,4X,A1)') '|', 'Rotation about each axis in '
     &   //'degrees and rms deviation in angstroms is', '|'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(A2,14X,3(A1,15X),A3,11X,A1)') '|', 'X', 'Y', 'Z',
     &   'rms', '|'
      WRITE (12,'(A2,3X,4(F16.6),9X,A1)') '|', THETAX, THETAY, THETAZ,
     &   SQRT(EIGVAL(4) / N1), '|'
      WRITE (12,'(A2,76X,A1)') '|', '|'
      WRITE (12,'(1X,A78)') '+-----------------------------------------'
     &   //'-----------------------------------+'

      REWIND (12)

      DO 1 N = 1, 38
         READ (12,'(A80)') LINE
         WRITE (*,'(A80)') LINE
1     CONTINUE

      CLOSE (12)

      RETURN
      END

C=======================================================================