	use MSFLIB
c uncomment the above line for Microsoft Fortran compiler
c
c     This program was written by William S. "Bill" Bennett.
c     It is unpublished, but Bill should be cited if it is used.
c     The superposition code was written by S. James Remington
c
C     OVRLAP:  A PROGRAM TO CALCULATE EQUIVALENCE PROBABILITIES  FOR
C     TWO PROTEIN STRUCTURES USING THE PROCEDURES OF ROSSMANN &
C     ARGOS [J. BIOL. CHEM. 250, P. 7525 (1975)].
C
C     GENERAL DESCRIPTION
C
C     THIS PROGRAM SUPERIMPOSES ONE STRUCTURE (THE "TRIAL" STRUC-
C     TURE) ON THE OTHER (THE "REFERENCE" STRUCTURE) WITH A RIGID-
C     BODY TRANSFORMATION AND THEN CONSTRUCTS A TABLE FOR EACH TRIAL
C     RESIDUE OF THE MOST-PROBABLY EQUIVALENT REF RESIDUES (1<JSAVE<
C     10 OF THEM).  THE "PROBABILITY" OF EQUIVALENCE, PIJ, IS DEFINED
C     BELOW (SEE CARD 3).  "EQUIVALENT" RESIDUES ARE THEN CHOSEN IN
C     THREE WAYS:
C       1. SIMPLE PROGRESSION RULE.  THE "PROGRESSION RULE" (SEE
C     R&A REFERENCE ABOVE) IS APPLIED TO FORWARD RUNS OF MOST-PROB-
C     ABLY EQUIV RESIDUES; SURVIVING RUNS ARE THEN EXTENDED WITH LESS-
C     PROBABLY EQUIV RESIDUES IN PIJ TABLES OF ADJACENT RESIDUES.
C     THIS IS THE PROCEDURE OF ROSSMAN & ARGOS AND SHOULD GIVE SIMILAR
C     RESULTS.
C       2. EXTENDED PROGRESSION RULE.  THE PROGRESSION RULE IS
C     APPLIED TO EXTENDED FORWARD RUNS.  RUNS SMALLER THAN "IMINR"
C     ARE OMITTED AND A MORE RESTRICTIVE SELECTION ALGORITHM IS
C     USED (DESCRIBED IN "SELEQ" SOURCE).
C       3. WITHOUT PROGRESSION RULE.  THE SET OF NONOVERLAPPING
C     EXTENDED RUNS WITH THE LARGEST TOTAL PROBABILITY OF EQUIV
C     IS FOUND.  AGAIN RUNS SMALLER THAN IMINR ARE OMITTED AND
C     THE MORE RESTRICTIVE ALGORITHM IS USED.
C       ANY OF THE THREE EQUIVALENCE SETS CAN BE REFINED BY ITERA-
C     TIVELY SELECTING AN EQUIV SET AND THEN REFINING THE RIGID-
C     BODY TRANSFORM TO OBTAIN THE BEST FIT OF THE EQUIV CA ATOMS.
C     THE NEW TRANSFORM IS USED TO SELECT A NEW EQUIV SET IN THE
C     NEXT CYCLE.  ITERATIONS CONTINUE UNTIL (1) A PRESET LIMIT
C     "ITMAX" IS REACHED (SEE CARD 2), (2) THE EQUIVALENCE SET
C     "CONVERGES"--I.E., IS THE SAME IN TWO SUCCESSIVE ITERATIONS,
C     OR (3) THE CHANGES IN THE TRANSFORMATION PARAMETERS ARE LESS
C     THAN A PRESET MULTIPLE OF THEIR STANDARD ERRORS, "ESIG"
C     (SEE CARD 3).
C       THE USER MUST SPECIFY EITHER SOME EQUIVALENT ATOMS TO BE SUPER-
C     IMPOSED (POSSIBLY ONLY A SMALL SUBSET OF THOSE EXPECTED)
C
C     I/O ASSIGNMENTS:                             IBM      VAX
C                                                  ===      ===
C       TRIAL STRUCTURE PERMANENT BINARY FILE       1      TRIBIN 
C       REFERENCE STRUCTURE PERMANENT BINARY FILE   2      REFBIN 
C       TRIAL STRUCTURE BINARY SCRATCH FILE         3    "trial.tmp"
C       REFERENCE STRUCTURE BINARY SCRATCH FILE     4    "ref.tmp"
C       SYSTEM INPUT STREAM                         5    SYS$INPUT
C       SYSTEM PRINT STREAM                         6    SYS$OUTPUT
C       SYSTEM PUNCH (IBM) FOR MATRIX               7   "RTN.DAT"
C       BINARY SCRATCH FILE FOR PCALC               9    "pcl.tmp"
C       TRIAL STRUCTURE CARD INPUT        ICARD OR IIN     TRICRD
C       REFERENCE STRUCTURE CARD INPUT    ICARD OR JIN     REFCRD
C
C     PROGRAMMING NOTES:
C
C     --SEE ALSO INDIVIDUAL SUBROUTINES.
C     --ARRAY LIMITS:  THE PROGRAM IS CURRENTLY LIMITED TO 
C     500 ATOMS (EACH OF TRIAL & REF STRUCTURES); SEE
C     "GETCRD" SOURCE FOR DETAILS ON CHANGING THESE.
C     --IT IS POSSIBLE (BUT NOT IMPLEMENTED) TO OVERLAY
C     THE /COORDS/ ARRAYS EXCEPT IN "GETCRD" AND "LSCTE";
C     SEE "LSCE" SOURCE FOR DETAILS.
C     --THE DISPLAY OF RUNS & EQUIVS (SEE "RLIST") ALLOWS ROOM
C     TO LIST ALSO THE SEQUENCE; THIS IS NOT IMPLEMENTED, THO.
C
C
      LOGICAL NPL,PTMP
      REAL*8 TXT8,BIN,CARD,PLIT(3),AEL1,AML1
      COMMON /PARAMS/ ICARD,IDISK,ITRAN,JSAVE,NPPRT,IMINR,IPROG,ISTAR,
     1ITMAX,NCYCLE,NEQ,IECRD,IHOLD,IDRAW,IORNT,NPL,DCUT,PCUT,E1,E2,
     2e1m,e2m,ITER,NI,NJ,NCHNG
      COMMON /XFRMS/ A(3,3),OA(3,3),OR(3),AS,AT,AP,V(3),ESIG,TSIG
      COMMON /DIREC/ NRUN,NSRUN,ISFIN(60),JSFIN(60),ISST(60),JSST(60),
     1JOF(60),ISTRT(60),JSTRT(60),IFIN(60),JFIN(60),SSPIJ(60),SPIJ(60),
     2SD2(60),SS2(60),SDX(60),SDY(60),SDZ(60),JCNT(60),IPRG1(60),
     4IPRG2(60),NPROG(60),ISPR(60),IEPR(60),INPR(60),NEXT(60),IGAP(60)
      COMMON /T99COM/ STRNG(64)
      DATA BIN,CARD,BLANK,RNOT/'BINARY  ','CARD    ','    ','NOT '/
      DATA PLIT,AEL1,AEL2,AML1,AML2/' SIMPLE ','EXTENDED','   NO   ',
     1'EQUIVALE','NCE ','TRANSFOR','M   '/
C
C
c      READ(5,1001) ICARD,IDISK,ITRAN,JSAVE,NPPRT,IMINR,IPROG,ISTAR,
c     1ITMAX,NCYCLE,NEQ,IECRD,IHOLD,IDRAW,IORNT
c1001  FORMAT(16I5)
	icard=1
	idisk=0
	itran=0
	jsave=5
	npprt=0
	iminr=3
	iprog=1   ! progression rule see comments ***
	istar=1
	itmax=8   ! maximum equivalence iterations***
	ncycle=1  ! (NOT USED) 
	neq=0     ! set neq later *** !
	iecrd=0
	ihold=1
	idraw=0
	iornt=0
c
c open list file
c
	open(unit=8,name='ovrlap.lis',status='unknown')
	write(8,1009)
1009	format(' OVRLAP version 2.0:  A program to align two protein' 
     1 ,' structures using the method of'/ 
     2 ' Rossmann & Argos, ',
     3 ' J. Biol. Chem. 250, 7525 (1975)'/)
c
C - - - - - - - - -  INPUT PARAMETERS  - - - - - - - - - - - - - - - - -
C
C
C       ICARD = LOGICAL UNIT NUMBER FOR CARD INPUT (1 FOR VAX)
C             = 0 FOR BINARY INPUT
C             = -1 TO INVOKE EDITING OF INPUT CARDS (SEE CARD 14)
C                  (SEE LOGICAL UNIT ASSIGNMENTS BELOW)
C
C       IDISK = 1 TO WRITE BINARY FILES FOR FUTURE RUNS
C
C       ITRAN = 0 TO CALC THE INITIAL TRANSFORM FROM THE INPUT 
C                 EQUIVALENCE (SEE CARD 14)
C             = 1 TO USE THE INPUT TRANSFORM TO CALCULATE AN INITIAL
C                 EQUIVALENCE
C
C       JSAVE = NUMBER OF LARGEST P(I,J) TO BE KEPT FOR EACH I
C               JSAVE MUST BE IN THE RANGE (1,10). GENERALLY JSAVE=5
C               IS ADEQUATE, JSAVE=1 FOR STRUCTURES WITH A ONE-TO-ONE
C               CORRESPONDENCE BETWEEN RESIDUES.
C             = 0 TO SKIP PROBABILITY AND EQUIVALENCE CALCS
C               ...THIS ALLOWS YOU TO PERFORM ONLY:
C               --BINARY FILE CREATION (ITRAN=IDRAW=0)
C               --LEAST-SQUARES REFINEMENT OF TRANSFORMATION BETWEEN
C               ATOMS IN "INIT EQUIVALENCE" (ITRAN=1,NEQ>0)
C               --ORTEP CARD GENERATION (IDRAW=1)
C
C       NPPRT = 1 TO PRINT ALL (KEPT) P(I,J) IN ADDITION TO LINEAR
C                 DISPLAY OF RUNS AND EQUIVALENCES; THIS IS ONLY
C                 USEFUL FOR DIAGNOSTICS OR SEEING HOW THINGS WORK.
C             = 0 TO PRINT LINEAR DISPLAY OF RUNS & EQUIVS ALONE.
C               (NPPRT APPLIES ONLY TO THE LAST (OR ONLY) ITERATION)
C
C       IMINR = MINIMUM LENGTH ALLOWED FOR A RUN (APPLIED ONLY TO
C               EXTENDED-PROG-RULE AND NO-PROG-RULE EQUIVS);
C               --SUGGESTED: 3 (TO ELIMINATE ONLY 2-RES RUNS)
C
C       IPROG = 0 NO ITERATIVE REFINEMENT OF THE EQUIVALENCE
C             = 1 THE SIMPLE PROGRESSION RULE EQUIVALENCE IS REFINED
C             = 2 THE EXTENDED PROGRESSION RULE EQUIV IS REFINED
C             = 3 THE EQUIV WITHOUT PROGRESSION RULE IS REFINED
C
C       ISTAR = FIRST CYCLE OF ITERATION TO BE PERFORMED (COSMETIC)
C               NORMALLY 1, BUT YOU CAN CONTINUE ITERATIVE REFINEMENT
C               OF AN EQUIV WHERE IT LEFT OFF BY STARTING FROM THE
C               OUTPUT EQUIVALENCE.  THEN ISTAR CAN BE SET TO THE
C               LAST ITERATION +1 AS A REMINDER.
C
C       ITMAX = MAXIMUM NUMBER OF EQUIVALENCE ITERATIONS TO BE
C               PERFORMED.  THIS SERVES ONLY TO PREVENT THE PROGRAM
C               FROM USING TOO MUCH CPU TIME; THE TRANSFORMATION
C               NORMALLY ATTAINS STABLE VALUES AFTER THE FIRST FEW
C               ITERATIONS, BUT THE REFINEMENT OF THE EQUIV CAN GO
C               ON FOR MANY ITERATIONS MAKING SMALL CHANGES IN THE
C               EQUIV AND IN E1 AND E2.  NORMALLY THE ITERATIVE
C               REFINEMENT IS STOPPED WITH THE "ESIG" PARAMETER
C               (CARD 3); THE PROGRAM WILL ALSO STOP IF THE EQIV
C               CONVERGES.
C
C       NCYCLE= NUMBER OF CYCLES OF NON-LINEAR LEAST-SQUARES REFINEMENT
C               OF ORIENTATION TO PERFORM ON EACH ITERATION.  THIS IS
C               IS ALSO A SAFETY-VALVE LEFT OVER FROM EARLIER VERSIONS
C               OF THE PROGRAM; FOR THIS VERSION, SET NYCYCLE LARGE
C               (10 OR 15) AND STOP THE LEAST-SQUARES PROCEDURE WITH
C               THE PARAMETER "TSIG" (CARD 3).
C             < 0 TO PRINT A FULL LIST OF SUPERIMPOSED ATOMS AND
C                 AND DISTANCES BETWEEN THEM AFTER THE LEAST-
C                 SQUARES PROCEDURE; OTHERWISE THIS IS DONE ONLY
C                 IN THE FIRST ITERATION WHEN STARTING FROM AN
C                 INPUT EQUIVALENCE.  ABS(NCYCLE) IS STILL USED
C                 TO LIMIT THE NUMBER OF LEAST-SQUARES CYCLES.
C
C       NEQ   = NUMBER OF EQUIVALENT PAIRS OF RESIDUES TO BE READ
C               FROM CARDS.  YOU MAY READ IN AS MANY EQUIVALENT PAIRS
C               AS THERE ARE ATOMS AND AS FEW AS THREE; IF THERE ARE
C               FEWER THAN 6 PAIRS, NO LEAST-SQUARES REFINEMENT IS
C               PERFORMED, BUT A QUICK-AND-DIRTY ESTIMATE OF THE 
C               TRANSFORM CAN BE MADE (SEE CARD 5) AND EQUIV RESIDUES
C               CALCULATED WITH THAT.
C
C       IHOLD = NUMBER OF ITERATIONS TO HOLD E1 AND E2 AT INPUT VALUES
C             = 0 TO RESET E1 AND E2 WITH RMS D AND RMS S (RESPECTIVELY)
C                 FROM LSCTE IN EACH CYCLE
C               --IHOLD SHOULD BE SET TO 1 WHEN STARTING FROM A SMALL
C               INPUT EQUIVALENCE TO PREVENT E1 AND E2 FROM BEING
C               RESET TO UNREASONABLY SMALL VALUES, AND MAY BE SET TO
C               ITMAX TO PREVENT E1 AND E2 FROM BEING RESET AT ALL.
C
C
c      READ(5,1051) DCUT,PCUT,E1,E2,ESIG,TSIG
c1051  FORMAT(6F10.5)
c
c set reasonable defaults
c
	dcut=8.
	pcut=0.1
	e1=5.
	e1m=0.
	e2=5.
	e2m=0.
	esig=1.0
	tsig=0.5
C
C       DCUT = MAXIMUM D(I,J) TO BE INCLUDED IN PROBABILITY CALCS.
C              D(I,J) IS THE DISTANCE IN ANGSTROMS BETWEEN THE ITH TRIAL
C              CA AND THE JTH REF CA AFTER APPLICATION OF THE TRANSFORM.
C              THIS PARAMETER SPEEDS UP THE PROGRAM BY SKIPPING
C              ATOMS THAT ARE MUCH TOO FAR APART TO BE EQUIVALENT.
C              --SUGGESTED: 8.5 ANGSTROMS
C
C       PCUT = MIN P(I,J) TO BE INCLUDED IN EQUIVALENCE CALCS
C        E1  = WEIGHTING APPLIED TO DISTANCE COMPONENT OF PROBABILITY
C        E2  = WIEGHTING APPLIED TO SCATTER COMPONENT OF PROBABILITY
C              --THESE THREE PARAMETERS INTERACT TO DETERMINE HOW FAR
C              A RUN IS EXTENDED.  THE PROBABILITY THAT THE ITH TRIAL
C              RES IS EQUIVALENT TO THE JTH REF RES IS THAT DEFINED BY
C              ROSSMANN & ARGOS:
C                P(I,J) = EXP(-D(I,J)**2/E1**2)*EXP(-S(I,J)**2/E2**2)
C              WHERE D(I,J) IS AS ABOVE AND S(I,J) IS A "SCATTER" TERM
C              THAT IS THE SUM OF THE VECTOR DIFFERENCES BETWEEN THE
C              VECTORS CA(I)-CA(I-1) AND CA(I+1)-CA(I) IN THE TRIAL
C              AND REF STRUCTURES.  THIS QUANTITY ALSO HAS DIMENSIONS OF
C              LENGTH AND IS NUMERICALLY ABOUT THE SAME AS D(I,J).  THUS
C              E1 AND E2 ARE THE SAME AS USED BY ARGOS & ROSSMANN.
C              --SUGGESTED VALUES: SET PCUT=0.1; THEN E1 AND E2 CAN BE
C              IN THE RANGE 1.5-4.0 ANGSTROMS, DEPENDING ON THE ACCURACY
C              OF YOUR COORDS.
C
C       ESIG = SIGNIFICANCE LEVEL OF LSQ SHIFTS REQUIRED FOR CONTINUED
C              EQUIV ITERATION, IN MULTIPLES OF THE STANDARD DEVIATION
C              OF SHIFTS; ITERATIONS ARE STOPPED IF ALL SHIFTS OF
C              1ST LSQ CYCLE ARE LESS THAN ESIG*SIG(I).
C              --SUGGESTED: 2.0
C
C       TSIG = SIGNIFICANCE LEVEL OF LSQ SHIFTS REQUIRED FOR CONTINUED
C              LSQ REFINEMENT OF TRANSFORM, IN MULTIPLES OF ESD.
C              --SUGGESTED: 1.0
10    WRITE(*,1008) DCUT,PCUT,E1,E2
1008  FORMAT(' Program OVRLAP'/
	1' Input parameters: DCUT = ',F5.1,', PCUT = ',F6.3,
     2', E1 = ',F5.1,', E2 = ',F5.1)
      WRITE(*,1053) JSAVE,IMINR
1053  FORMAT(' The top ',I2,' probabilites will be saved',
     1'. Minimum run length ',I2)
14    NPL = .FALSE.
      PTMP = .FALSE.
      IF(NPPRT .NE. 1) then
      NPL = .TRUE.
      PTMP = .TRUE.
      WRITE(*,1052)
1052  FORMAT(' No listing of P(I,J)')
	end if
C
C     MAIN PROGRAM. Get coordinates and assign initial equivalences
C
      CALL GETCRD
	write(*,11)
11	format(' New E-values: initial, lower limit: (5,5,0,0)?',$)
	read(*,12) E1N,E2N,E1MN,E2MN
12	format(4f10.0)
	if(E1N.gt.0.) e1=e1n
	if(E2N.gt.0.) e2=e2n
	if(E1MN.gt.0.) e1m=e1mn
	if(E2MN.gt.0.) e2m=e2mn
	write (*,*)
	write(8,13) e1,e2,e1m,e2m
13	format(//5x,'E-values, initial, lower limit:'4f6.2)
c
C     SUPPRESS PIJ LISTING ON ALL BUT LAST ITERATION
c
      IF(IPROG .NE. 0) NPL = .TRUE.
      IF(IPROG .EQ. 0) ITMAX = 1
C
C     EQUIVALENCE ITERATION LOOP
C
	DO 25 IT=1,ITMAX
      ITER = ISTAR + IT - 1
      CALL LSCTE
      CALL PCALC
      CALL RUNS
      IF(NRUN .EQ. 0) stop 'No runs found'
      CALL SELEQ
C     SELEQ SETS ITER BACK TO 0 UNLESS THERE IS NO CHANGE IN EQUIVS
      IF(ITER .NE. 0) GO TO 26
25    CONTINUE
c
26    IF(PTMP) GO TO 27
      NPL = .FALSE.
      ITER = ISTAR + IT - 1
      CALL PCALC
27    CALL RLIST
C
28    CONTINUE
C 
C     DELETE TEMPORARY COORD FILES
505   OPEN (UNIT=3,NAME='trial.tmp',status='OLD',FORM='UNFORMATTED')
      CLOSE (UNIT=3,status='DELETE')
      OPEN (UNIT=4,NAME='ref.tmp',status='OLD',FORM='UNFORMATTED')
      CLOSE (UNIT=4,status='DELETE')
c
c report rotation and translation
c
	CALL MATCON(A,V)
c
	write(*,1001)
1001	format(/' Output written to OVRLAP.LIS'/
	1 ' Matrix written to RTN.DAT')
c
      close (unit=8)
	pause 'Hit Enter to quit'
	STOP
      END
      SUBROUTINE GETCRD
C
C     GET CA COORDS
C
C *** ARRAY LIMITS CURRENTLY SET FOR A MAXIMUM OF 500 CA ATOMS
C     TO CHANGE LIMTS TO NEW MAXIMUM, "NMX":
C       (1) SET ARRAY LIMITS IN VARIOUS OCCURENCES OF COMMON STATEMENTS
C           /COORDS/ XI(3,NMX+2),XJ(3,NMX+2)
C           /SEQ/ IRES(NMX+2),JRES(NMX+2)
C           /EQUIV/ IEQ(NMX+2)
C           /PIJ/ IJ(NMX+2,10)
C           /SAVE/ NSV(NMX+2)
C       (2) SET ARRAY LIMITS IN SUBROUTINES
C           RUNS:  JMP(NMX+2)
C           SELEQ:  IOLD(NMX+2)
C           DRAW:  IAT(NMX+2),JAT(NMX+2)
C       (3) CHANGE NMX AND NIJMX=NMX+1 IN DATA STATEMENT OF GETCRD
C       (4) CHANGE COMMENT IN GETCRD
C
      LOGICAL NPL,LOK
      COMMON /PARAMS/ ICARD,IDISK,ITRAN,JSAVE,NPPRT,IMINR,IPROG,ISTAR,
     1ITMAX,NCYCLE,NEQ,IECRD,IHOLD,IDRAW,IORNT,NPL,DCUT,PCUT,E1,E2,
     2e1m,e2m,ITER,NI,NJ,NCHNG
      COMMON /COORDS/ XI(3,702),XJ(3,702)
      COMMON /SEQ/ IRES(702),JRES(702),AIRES(702),AJRES(702)
      COMMON /EQUIV/ IEQ(702),IQND(2,3)
      COMMON /XFRMS/ A(3,3),OA(3,3),OR(3),AS,AT,AP,V(3),ESIG,TSIG
      DIMENSION VT(3),INEQ(16),IS(3),IF(3),JS(3),JF(3),initeq(2,10)
      EQUIVALENCE (X,VT(1)),(Y,VT(2)),(Z,VT(3))
	character*60 tricrd,refcrd
	character*80 line
	real card(20),x(3)
c
	integer ca
	DATA CA/'CA  '/
      DATA NMX,NIJMX/700,701/
      DATA IS,JS,IF,JF,NOMAT/13*')()('/
c
	data nb/4h    /,atom/4hATOM/,het/4hHETA/
C
      NI = 1
C
C     NOTE THAT NI, NJ = NUMBER OF ATOMS + 1
C	*******Coordinates run from 2 to NI,NJ*******
C
	write(*,*)
	write(*,1) 'Trial coordinate file name: '
1	format(1x,a,$) 
	read(*,'(a)') tricrd
	write(*,*)
	write(8,8) '    Trial coordinate filename: ',tricrd
8	format(1x,2a)
      OPEN (UNIT=1,NAME=TRICRD,status='old')
c
c
c loop over input file
cTOM      1  N   ALA     1     -18.102   5.182   4.182  1.00 20.35      5CPA 156
c
20    read(1,100,end=220) line
100   format(a)
c
	if(line(1:4).ne.'ATOM') goto 20
c                   atom res   resn
	read(line,fmt=101) nam,rtype,nr,x,w,b
101	format(13x,a4,a3,2x,a4,4x,3f8.3,2f6.2,f6.0)
cTOM   1115  OH  OXY   169      27.369  37.098  19.298  1.00100.00
c
	if (nam .eq. ca) then
	nr=ljust(nr)

c 19    THR  CA       22.935    42.432    23.611
c 20    ILE  CA       22.230    40.439    26.735
c 21    READ(1,1011,END=29) NRES,RES,ATOM,X,Y,Z
c1011  FORMAT(A4,4X,A4,1X,A4,1X,3F10.5)
      NI = NI + 1
      IRES(NI) = nr       !NRES
      AIRES(NI) = rtype   !RES
      XI(1,NI) = x(1)
      XI(2,NI) = x(2)
      XI(3,NI) = x(3)
c	if(ni.lt.7) write(*,101) nam,rtype,nr,x
	end if
      IF(NI .lt. NIJMX) GO TO 20
	if(NI  .eq. NIJMX) WRITE(*,1069) NMX,IRES(NIJMX)
1069  FORMAT(' INPUT FILE TRUNCATED TO ',I3,
     1' RESIDUES, last accepted:',A4)
c
220   CLOSE (UNIT=1)
      N = NI - 1
      WRITE(*,1012) N
1012  FORMAT(1x,i4,' atoms in trial structure')
c
      NJ = 1
 	write(*,1) 'Reference coordinate file name: '
	read(*,'(a)') refcrd
	write(*,*)
	write(8,8) 'Reference coordinate filename: ',refcrd
      OPEN (UNIT=1,NAME=refcrd,status='old')
c
c loop over input file
cTOM      1  N   ALA     1     -18.102   5.182   4.182  1.00 20.35      5CPA 156
c
30    read(1,100,end=320) line
c
	if(line(1:4).ne.'ATOM') goto 30
c                   atom res   resn
	read(line,101) nam,rtype,nr,x,w,b
c
	if (nam .eq. ca) then
	nr=ljust(nr)

      NJ = NJ + 1
      JRES(NJ) = nr
      AJRES(NJ) = rtype
      XJ(1,NJ) = x(1)
      XJ(2,NJ) = x(2)
      XJ(3,NJ) = x(3)
c	if(nj.lt.7) write(*,101) nam,rtype,nr,x
	end if
c
	IF(NJ .lt. NIJMX) GO TO 30
	if(NJ .eq. NIJMX) WRITE(*,1069) NMX,JRES(NIJMX)
320   CLOSE (UNIT=1)
      N = NJ - 1
      WRITE(*,1013) N
1013  FORMAT(1x,i4,' atoms in reference structure')
C
C     FILL IN ENDS OF DATA ARRAYS FOR SCATTER CALCS
C
      XI(1,1) = XI(1,2)
      XI(2,1) = XI(2,2)
      XI(3,1) = XI(3,2)
      N = NI + 1
      XI(1,N) = XI(1,NI)
      XI(2,N) = XI(2,NI)
      XI(3,N) = XI(3,NI)
      XJ(1,1) = XJ(1,2)
      XJ(2,1) = XJ(2,2)
      XJ(3,1) = XJ(3,2)
      N = NJ + 1
      XJ(1,N) = XJ(1,NJ)
      XJ(2,N) = XJ(2,NJ)
      XJ(3,N) = XJ(3,NJ)
C
C     KEEP CURRENT UNTRANSFORMED COORDS ON SCRATCH FILES
C
      OPEN (UNIT=3,NAME='trial.tmp',status='unknown',FORM='UNFORMATTED')
      WRITE(3) NI,IRES,XI
      CLOSE (UNIT=3)
      OPEN (UNIT=4,NAME='ref.tmp',status='unknown',FORM='UNFORMATTED')
      WRITE(4) NJ,JRES,XJ
      CLOSE (UNIT=4)
C
C     READ EQUIVALENCE FROM CARDS
C
	write(*,1) 'Number of initial equivalences? (5) '
	read(*,2) neq
2	format(i5)
	if(neq.le.0) neq=5
	if(neq.gt.10) neq=10
	write(*,*)
	write(*,3) neq
3	format(' Enter',i3,' residue IDs, one pair per line (trial, ref)')
	DO 43 I=2,NI
43    IEQ(I) = 0
c
	do 50 n=1,neq
	call get2a4(I1,J1)
C
C	*******Coordinates run from 2 to NI,NJ*******
C
	I2 = 0
      J2 = 0
      DO 46 I=2,NI
      IF(I1 .EQ. IRES(I) .and. i2.eq.0) I2=I
46    CONTINUE
      DO 48 J=2,NJ
      IF(J1 .EQ. JRES(J) .and. j2.eq.0) J2=J
48    CONTINUE
      IF(I2 .EQ. 0 .OR. J2 .EQ. 0) GO TO 45
      IEQ(I2) = J2
	initeq(1,n)=i1
	initeq(2,n)=j1
50	continue
	write(8,4) (initeq(1,i),initeq(2,i),i=1,neq)
4	format(//5x,'Starting equivalences: (trial, ref)'
	1//(5x,a4,'<-> ',a4))
c
      RETURN
C     ======
45    WRITE(*,1066) I1,I2,J1,J2
1066  FORMAT(' Residue number not found. Trial= ',A4,I4,
	1'  Reference= 'A4,I4)
      STOP
      END
	INTEGER FUNCTION ljust(NUM)
C
C LEFT JUSTIFY 4 CHARACTERS IN NUM
C
	BYTE S(4),BLANK
	INTEGER N
	EQUIVALENCE (N,S(1))
	DATA BLANK/1H /
	N=NUM
	DO I=1,3
	IF(S(1).eq.BLANK ) THEN
		DO J=2,4
		S(J-1)=S(J)
		END DO
		S(4)=BLANK
	    END IF
	END DO
	ljust=N
	RETURN
	END
	subroutine get2a4(I,J)
c
c	read a line from standard input and get two 4-character names into 4-byte integers I and J
c	result is left-justified
c
	integer line(5),m4
	byte n(20),m(4),blank,tab,comma
	equivalence (m4,m(1))
	equivalence (n(1),line(1))
	data blank,tab,comma/1h ,8,1h,/
	read(*,1) line
1	format(5a4)
c
c replace separators
c
	do l=1,20
	if(n(l).ne.blank) then
	    if(n(l).eq.comma .or. n(l).eq.tab) n(l)=blank
	end if
	end do
c
c find 1st group
c
	i1=0
	j1=0
	do l=1,20
	if(n(l).ne.blank .and. i1.eq.0) i1=l
	if(n(l).eq.blank .and. i1.ne.0) then
	    if(j1.eq.0) j1=l-1 !last nonblank
	    end if
	end do 
	if ((j1-i1).gt.3) j1=i1+3
	nc=0
	do l=1,4
	m(l)=blank
	end do
	do l=i1,j1
	nc=nc+1
	m(nc)=n(l)
	n(l)=blank !wipe out character extracted
	end do
	I=m4  !return I
c
c find 2nd group
c
	i1=0
	j1=0
	do l=1,20
	if(n(l).ne.blank .and. i1.eq.0) i1=l
	if(n(l).eq.blank .and. i1.ne.0) then
	    if(j1.eq.0) j1=l-1 !last nonblank
	    end if
	end do
	if ((j1-i1).gt.3) j1=i1+3
	nc=0
	do l=1,4
	m(l)=blank
	end do
	do l=i1,j1
	nc=nc+1
	m(nc)=n(l)
	n(l)=blank !wipe out character extracted
	end do
	J=m4 !return J
	return
	end

      SUBROUTINE LSCTE
C     LEAST SQUARES RIGID-BODY ORIENTATION REFINEMENT FOR EQUIVALENCE
C     CALCULATIONS.  THE PROGRAM REFINES THREE EULERIAN ANGLES
C     AND THREE TRANSLATIONAL SHIFTS.
C     THE EULER ANGLES IN THIS VERSION ARE COUNTER-CLOCKWISE ROTATIONS
C     OF THE COORDINATE SYSTEM:
C       PSI   AROUND THE Z-AXIS
C       THETA AROUND THE Y'-AXIS
C       PHI   AROUND THE Z"-AXIS
C NOTATION:  THE "TRIAL" STRUCTURE IS TRANSFORMED TO THE "REFERENCE"
C
C     NEWER VERSION: STOPS LSQ ITERATIONS IF ALL SHIFTS LESS THAN
C     TSIG*SIG AND STOPS EQUIV ITERATIONS IF 1ST LSQ CYCLE GIVES
C     ALL SHIFTS LESS THAN ESIG*SIG.
C *** TEMP KLUDGE FOR E-SCAN: SEND MOST LSCTE OUTPUT TO UNIT 7
C     EXCEPT E-RESETS; UNIT 7 MUST BE ASSIGNED BEFORE RUN
C
      LOGICAL NPL,LSKP
      COMMON /PARAMS/ ICARD,IDISK,ITRAN,JSAVE,NPPRT,IMINR,IPROG,ISTAR,
     1ITMAX,NCYCLE,NEQ,IECRD,IHOLD,IDRAW,IORNT,NPL,DCUT,PCUT,E1,E2,
     2e1m,e2m,ITER,NI,NJ,NCHNG
      COMMON /COORDS/ XI(3,702),XJ(3,702)
      COMMON /SEQ/ IRES(702),JRES(702),AIRES(702),AJRES(702)
	COMMON /EQUIV/ IEQ(702),IQND(2,3)
      COMMON /XFRMS/ PTSM(3,3),OA(3,3),OR(3),AS,AT,AP,XO(3),ESIG,TSIG
      COMMON /MATJNK/ PM(3,3),TM(3,3),SM(3,3),PTM(3,3),
     1                CS,SS,CT,ST,CP,SP
      DIMENSION UP(3),U(3),XS(3),XSP(3),XP(3),XOBS(3),XCAL(3),DX(3),
     1          V1(3),W(21),V( 6),DXDP(3, 6),SDX(3)
      DIMENSION DXSV(15,3),DSSV(15),PSV(15,6),SIGSV(15,6),DPSV(15,6)
     1,SXO(3),P(6),DIAG(6),DP(6)
      DIMENSION QP(3),QQ(3),QR(3),QU(3),QV(3),QW(3)
      EQUIVALENCE (P(1),AS),(NEQ,NDATA)
      EQUIVALENCE (X,XCAL(1)),(Y,XCAL(2)),(Z,XCAL(3))
      EQUIVALENCE (PM(1,1),QP(1)),(PM(1,2),QQ(1)),(PM(1,3),QR(1)),
     1            (TM(1,1),QU(1)),(TM(1,2),QV(1)),(TM(1,3),QW(1))
	dimension xa(3,702),xb(3,702)
      DATA LSPRT/0/
	data ca/4hCA  /,ibl/4h    /
  520 FORMAT(/,(10F10.5))
C
C     GET NEW COPY OF COORDS (HAVE TO GET NEW TRIAL; REREAD REF
C     TO PERMIT OVERLAYING OF /COORDS/.)
C
      OPEN (UNIT=3,NAME='trial.tmp',status='OLD',FORM='UNFORMATTED')
      READ(3) NI,IRES,XI
      CLOSE (UNIT=3)
      OPEN (UNIT=4,NAME='ref.tmp',status='OLD',FORM='UNFORMATTED')
      READ(4) NJ,JRES,XJ
      CLOSE (UNIT=4)
	T=float(neq)
	np=0
	do 100 iobs=2,NI
	J = IEQ(IOBS)
      IF(J .EQ. 0) GO TO 100
	np=np+1
      xb(1,np) = XI(1,iobs)
      xb(2,np) = XI(2,iobs)
      xb(3,np) = XI(3,iobs)
	xa(1,np) = XJ(1,J)
      xa(2,np) = XJ(2,J)
      xa(3,np) = XJ(3,J)
100	continue
	if(np .ne. neq) stop 'LSCTE: NP is not equal to NEQ'
	if(np.lt.3) stop 'Too few points to fit, try another set'
c      SUBROUTINE QKFIT (N, XA, XB, R, V, RMSE)
      call qkfit(np,xa,xb,ptsm,xo,rmse)
	write(*,1) np,' points for QKFIT, rmsD = ',rmse
1	format(1x,i3,a,f5.2)
C
C     APPLY TRANSFORMATION TO TRIAL COORDS
C     DISPLAY INDIVIDUAL XO, XC, DIST AND CORREL MTX...
C
c      WRITE(8,181)
181   FORMAT(//' Results: '/
     3,18X,'Trial',33X,'Reference'/
     1'  NR  RES  ATOM',4X,'XCAL',3X,'YCAL',3X,'ZCAL',5X,'NR  RES  ATOM'
     2,4X,'XOBS',3X,'YOBS',3X,'ZOBS',8X,'DX',6X,'DY',6X,'DZ',14X,'D')
C
3133  S2 = 0.0
      SDX(1) = 0.0
      SDX(2) = 0.0
      SDX(3) = 0.0
      NN = 0
      INXT = 0
      SCAT2 = 0.0
	nrmsd = 0
      DO 8020 IOBS=2,NI
      CALL MV(PTSM,XI(1,IOBS),XP)
      X = XP(1) + XO(1)
      Y = XP(2) + XO(2)
      Z = XP(3) + XO(3)
      XI(1,IOBS) = X
      XI(2,IOBS) = Y
      XI(3,IOBS) = Z
C     PRINT ONLY ATOMS IN THE EQUIVALENCE; SKIP PRINT IF NO CHANGE
      J = IEQ(IOBS)
      IF(J .EQ. 0) GO TO 8020
	nrmsd = nrmsd + 1
      XOBS(1) = XJ(1,J)
      XOBS(2) = XJ(2,J)
      XOBS(3) = XJ(3,J)
      DO 8030 I=1,3
      DX(I) = XOBS(I) - XCAL(I)  !XCAL(I) is equivalenced to X,Y,X
8030  SDX(I) = SDX(I) + DX(I)*DX(I)
      DS2 = VV(DX,DX)
      S2 = S2 + DS2
      DS = SQRT(DS2)
c     WRITE(8,190) IRES(IOBS),IBL,CA,XCAL,JRES(J),IBL,CA,XOBS
c    1,DX,DS
190   FORMAT(1X,2(A4,1X,A4,1X,A4,1X,3F7.2,4X),3F8.3,6X,F10.3)
8020  CONTINUE
	rmsd = sqrt(s2/float(nrmsd))
C
C      ESTIMATE SCATTER FOR OVRLAP
C
	nscat=0
	DO 8040 IOBS=2,NI
      J=IEQ(IOBS)
      IF(J .EQ. 0) GO TO 8040
      nscat=nscat+1
	X = XI(1,IOBS)
      Y = XI(2,IOBS)
      Z = XI(3,IOBS)
      I1 = IOBS - 1
      DXM = XI(1,I1) - X
      DYM = XI(2,I1) - Y
      DZM = XI(3,I1) - Z
      I1 = IOBS + 1
      DXP = XI(1,I1) - X
      DYP = XI(2,I1) - Y
      DZP = XI(3,I1) - Z
      US = XJ(1,J)
      VS = XJ(2,J)
      WS = XJ(3,J)
      J1 = J - 1
      SXM = XJ(1,J1) - US
      SYM = XJ(2,J1) - VS
      SZM = XJ(3,J1) - WS
      J1 = J + 1
      SXP = XJ(1,J1) - US
      SYP = XJ(2,J1) - VS
      SZP = XJ(3,J1) - WS
C     "SP" SCATTER WILL APPLY IF A RUN HS JOFF = + 1
      T1 = DXM - SXM
      T2 = DYM - SYM
      T3 = DZM - SZM
      T4 = DXP - SXP
      T5 = DYP - SYP
      T6 = DZP - SZP
      SP2 = T1*T1 + T2*T2 + T3*T3 + T4*T4 + T5*T5 + T6*T6
C     "SN" SCATTER WILL APPLY IF A RUN HAS JOFF = -1
      T1 = DXM - SXP
      T2 = DYM - SYP
      T3 = DZM - SZP
      T4 = DXP - SXM
      T5 = DYP - SYM
      T6 = DZP - SZM
      SN2 = T1*T1 + T2*T2 + T3*T3 + T4*T4 + T5*T5 + T6*T6
      SS2 = SP2
      IF(SN2 .GE. SP2) GO TO 109
      SS2 = SN2
109   SCAT2 = SCAT2 + SS2
8040  CONTINUE
C
C SCATTER
C
      SCAT = SQRT(SCAT2/float(nscat))
      WRITE(8,9104) rmsd,scat,nscat
9104  FORMAT(/5x,'RmsD = 'f6.2, ' and estimated scatter =',F6.2,
	1 ' for',i4,' equivalences')
c
      WRITE(8, 545) ((PTSM(K,J),J=1,3),K=1,3),XO
545   FORMAT(/5x,'Rotation matrix and vector',
     1 /3(/5x,3F10.4)//5x,3F10.2)
C
	open(unit=7,name='rtn.dat',status='unknown')
	WRITE(7,1545) ((PTSM(K,J),J=1,3),K=1,3),XO
1545  FORMAT(4(3f10.5/))
	close(unit=7)
C
C     UPDATE E1 & E2 TO ESTIMATES DERIVED FROM NEWLY-XFRMED COORDS
C
C
9060  IF(ITER - ISTAR + 1 .LE. IHOLD) RETURN
      E1 = amax1(rmsd,e1m)
      E2 = amax1(scat,e2m)
      WRITE(8,9107) E1,E2
9107  FORMAT(//' *** RESET E1 to',F6.2,' and E2 to',F6.2, ' ***')
      RETURN
      END
      SUBROUTINE QKFIT (N, XA, XB, R, V, RMSE)
      DIMENSION   XA(3,N), XB(3,N), R(3,3), V(3)
C
C     SUBROUTINE TO FIT THE COORD SET XA(3,N) TO THE SET XB(3,N)
C     IN THE SENSE OF:
C            XA= R*XB +V
C     R IS A UNITARY 3.3 RIGHT HANDED ROTATION MATRIX
C     AND V IS THE OFFSET VECTOR. THIS IS AN EXACT SOLUTION
C
C
C     THIS SUBROUTINE IS A COMBINATION OF MCLACHLAN'S AND KABSCH'S
C     TECHNIQUES. SEE
C     KABSCH, W. ACTA CRYST A34, 827,1978
C     MCLACHAN, A.D., J. MOL. BIO. NNN, NNNN 1978
C     WRITTEN BY S.J. REMINGTON 11/78.
C
C     THIS SUBROUTINE USES THE IBM SSP EIGENVALUE ROUTINE 'EIGEN'
C
      DIMENSION CMA(3),CMB(3),UMAT(3,3)
      XN=N
      XASQ=0.0
      XBSQ=0.0
      XNI=1.0/XN
C
C     ACCUMULATE UNCORRECTED (FOR C.M.) SUMS AND SQUARES
C
      DO 40 I=1,3
      CMA(I) = 0.0
      CMB(I) = 0.0
      DO 10 J=1,3
   10 UMAT(I,J) = 0.0
C
      DO 30 J=1,N
      DO 20 K=1,3
   20 UMAT(I,K) = UMAT(I,K) + XA(I,J)*XB(K,J)
C
      T = XA(I,J)
      CMA(I) = CMA(I) + T
      XASQ = XASQ + T*T
      T = XB(I,J)
      CMB(I) = CMB(I) + T
      XBSQ = XBSQ + T*T
   30 CONTINUE
   40 CONTINUE
C
C     SUBTRACT CM OFFSETS
C
      DO 50 I=1,3
      XASQ = XASQ - CMA(I)*CMA(I)*XNI
      XBSQ = XBSQ - CMB(I)*CMB(I)*XNI
      DO 50 J=1,3
      UMAT(I,J) = (UMAT(I,J) - CMA(I)*CMB(J)*XNI)*XNI
   50 CONTINUE
C
C     FIT IT
C
      CALL QKFIT2( UMAT, RTSUM, R)
      RMSE =(XASQ + XBSQ)*XNI - 2.0*RTSUM
      IF( RMSE .LT. 0.0 ) RMSE = 0.0
      RMSE = SQRT (RMSE)
C
C      CALCULATE OFFSET
C
      DO 70 I=1,3
      T = 0.0
      DO 60 J=1,3
   60 T = T + R(I,J)*CMB(J)
      V(I) = ( CMA(I) - T )*XNI
   70 CONTINUE
      RETURN
      END
      SUBROUTINE QKFIT2 (UMAT, RTSUM, R)
      DIMENSION UMAT(3,3), R(3,3)
      DIMENSION ROOT(3),A(3,3),B(3,3),UTR(6)
C
      ISIG=1
C
C     FORM USQ = (UT).U    (IN UPPER TRIANGULAR SYMMETRIC STORAGE MODE)
C
      DO 220 I=1,3
      DO 220 J=I,3
      T = 0.0
      DO 210 K=1,3
      T = T + UMAT(K,I)*UMAT(K,J)
  210 CONTINUE
      IA = I + (J*J-J)/2
      UTR(IA) = T
  220 CONTINUE
C
C     CALCULATE EIGENVALUES AND VECTORS
C
      CALL EIGEN (UTR, A, 3, 0)
C
      ROOT(1) = UTR(1)
      ROOT(2) = UTR(3)
      ROOT(3) = UTR(6)
C%    WRITE(6,999) ROOT
C%    WRITE(6,999) A
C
C     SET A3 = A1 CROSS A2
C     ROOTS ARE IN ORDER R(1) >= R(2) >= R(3) >= 0
C
      A(1,3) = A(2,1)*A(3,2) - A(3,1)*A(2,2)
      A(2,3) = A(3,1)*A(1,2) - A(1,1)*A(3,2)
      A(3,3) = A(1,1)*A(2,2) - A(2,1)*A(1,2)
C%    WRITE(6,999) A
C
C     VECTOR SET B=U.A
C
      DO 240 I=1,3
      DO 240 J=1,3
      T = 0.0
      DO 230 K=1,3
  230 T = T + UMAT(J,K)*A(K,I)
      B(J,I) = T
  240 CONTINUE
C
C      NORMALIZE B1 AND B2 AND CALCULATE B3 = B1 CROSS B2
C
      B1 = SQRT( B(1,1)*B(1,1) + B(2,1)*B(2,1) + B(3,1)*B(3,1) )
      B2 = SQRT( B(1,2)*B(1,2) + B(2,2)*B(2,2) + B(3,2)*B(3,2) )
      DO 250 I=1,3
      B(I,1) = B(I,1)/B1
  250 B(I,2) = B(I,2)/B2
C
C      CHECK FOR LEFT HANDED ROTATION
C
      B13 = B(2,1)*B(3,2) - B(3,1)*B(2,2)
      B23 = B(3,1)*B(1,2) - B(1,1)*B(3,2)
      B33 = B(1,1)*B(2,2) - B(2,1)*B(1,2)
C
      S = B13*B(1,3) + B23*B(2,3) + B33*B(3,3)
      IF (S .LT. 0.0) ISIG = -1
      B(1,3) = B13
      B(2,3) = B23
      B(3,3) = B33
C%    WRITE(6,999) B
C
C     CALCULATE ROTATION MATRIX R
C
      DO 270 I=1,3
      DO 270 J=1,3
      T = 0.0
      DO 260 K=1,3
  260 T = T + B(I,K)*A(J,K)
      R(I,J) = T
  270 CONTINUE
C
C     RMS ERROR
C
      DO 280 I=1,3
      IF (ROOT(I) .LT. 0.0) ROOT(I) = 0.0
      ROOT(I) = SQRT( ROOT(I) )
  280 CONTINUE
C
C     CHANGE SIGN OF EVAL #3 IF LEFT HANDED
C
      IF (ISIG .LT. 0) ROOT(3)=-ROOT(3)
      RTSUM = ROOT(3) + ROOT(2) + ROOT(1)
      RETURN
      END
      SUBROUTINE EIGEN(A,R,N,MV)
      DIMENSION A(1),R(1)
C     SUBROUTINE EIGEN
C
C        PURPOSE
C           COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC
C           MATRIX
C
C        USAGE
C           CALL EIGEN(A,R,N,MV)
C
C        DESCRIPTION OF PARAMETERS
C           A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION.
C               RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF
C               MATRIX A IN DESCENDING ORDER.
C  P.S. THE EQUATION FOR THE INDICES OF A IS A(I,J) -> A(I+(J*J-J)/2))
C           R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE,
C               IN SAME SEQUENCE AS EIGENVALUES)
C           N - ORDER OF MATRICES A AND R
C           MV- INPUT CODE
C                   0   COMPUTE EIGENVALUES AND EIGENVECTORS
C                   1   COMPUTE EIGENVALUES ONLY (R NEED NOT BE
C                       DIMENSIONED BUT MUST STILL APPEAR IN CALLING
C                       SEQUENCE)
C
C        REMARKS
C           ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1)
C           MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED
C           BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN 'MATHEMATICAL
C           METHODS FOR DIGITAL COMPUTERS', EDITED BY A. RALSTON AND
C           H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C        STATEMENT WHICH FOLLOWS.
C
C     DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,
C    1                 COSX2,SINCS,RANGE
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C
C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  SQRT IN STATEMENTS
C        40, 68, 75, AND 78 MUST BE CHANGED TO DSQRT.  ABS IN STATEMENT
C        62 MUST BE CHANGED TO DABS. THE CONSTANT IN STATEMENT 5 SHOULD
C        BE CHANGED TO 1.0D-12.
C
C        ...............................................................
C
C        GENERATE IDENTITY MATRIX
C
    5 RANGE=1.0E-6
      IF(MV-1) 10,25,10
   10 IQ=-N
      DO 20 J=1,N
      IQ=IQ+N
      DO 20 I=1,N
      IJ=IQ+I
      R(IJ)=0.0
      IF(I-J) 20,15,20
   15 R(IJ)=1.0
   20 CONTINUE
C
C        COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX)
C
   25 ANORM=0.0
      DO 35 I=1,N
      DO 35 J=I,N
      IF(I-J) 30,35,30
   30 IA=I+(J*J-J)/2
      ANORM=ANORM+A(IA)*A(IA)
   35 CONTINUE
      IF(ANORM) 165,165,40
C? 40 ANORM=1.414D0*SQRT(ANORM)
   40 ANORM=1.41421356*SQRT(ANORM)
      ANRMX=ANORM*RANGE/FLOAT(N)
C
C        INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR
C
      IND=0
      THR=ANORM
   45 THR=THR/FLOAT(N)
   50 L=1
   55 M=L+1
C
C        COMPUTE SIN AND COS
C
   60 MQ=(M*M-M)/2
      LQ=(L*L-L)/2
      LM=L+MQ
   62 IF( ABS(A(LM))-THR) 130,65,65
   65 IND=1
      LL=L+LQ
      MM=M+MQ
      X=0.5*(A(LL)-A(MM))
   68 Y=-A(LM)/ SQRT(A(LM)*A(LM)+X*X)
      IF(X) 70,75,75
   70 Y=-Y
   75 SINX=Y/ SQRT(2.0*(1.0+( SQRT(1.0-Y*Y))))
      SINX2=SINX*SINX
   78 COSX= SQRT(1.0-SINX2)
      COSX2=COSX*COSX
      SINCS =SINX*COSX
C
C        ROTATE L AND M COLUMNS
C
      ILQ=N*(L-1)
      IMQ=N*(M-1)
      DO 125 I=1,N
      IQ=(I*I-I)/2
      IF(I-L) 80,115,80
   80 IF(I-M) 85,115,90
   85 IM=I+MQ
      GO TO 95
   90 IM=M+IQ
   95 IF(I-L) 100,105,105
  100 IL=I+LQ
      GO TO 110
  105 IL=L+IQ
  110 X=A(IL)*COSX-A(IM)*SINX
      A(IM)=A(IL)*SINX+A(IM)*COSX
      A(IL)=X
  115 IF(MV-1) 120,125,120
  120 ILR=ILQ+I
      IMR=IMQ+I
      X=R(ILR)*COSX-R(IMR)*SINX
      R(IMR)=R(ILR)*SINX+R(IMR)*COSX
      R(ILR)=X
  125 CONTINUE
      X=2.0*A(LM)*SINCS
      Y=A(LL)*COSX2+A(MM)*SINX2-X
      X=A(LL)*SINX2+A(MM)*COSX2+X
      A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
      A(LL)=Y
      A(MM)=X
C
C        TESTS FOR COMPLETION
C
C        TEST FOR M = LAST COLUMN
C
  130 IF(M-N) 135,140,135
  135 M=M+1
      GO TO 60
C
C        TEST FOR L = SECOND FROM LAST COLUMN
C
  140 IF(L-(N-1)) 145,150,145
  145 L=L+1
      GO TO 55
  150 IF(IND-1) 160,155,160
  155 IND=0
      GO TO 50
C
C        COMPARE THRESHOLD WITH FINAL NORM
C
  160 IF(THR-ANRMX) 165,165,45
C
C        SORT EIGENVALUES AND EIGENVECTORS
C
  165 IQ=-N
      DO 185 I=1,N
      IQ=IQ+N
      LL=I+(I*I-I)/2
      JQ=N*(I-2)
      DO 185 J=I,N
      JQ=JQ+N
      MM=J+(J*J-J)/2
      IF(A(LL)-A(MM)) 170,185,185
  170 X=A(LL)
      A(LL)=A(MM)
      A(MM)=X
      IF(MV-1) 175,185,175
  175 DO 180 K=1,N
      ILR=IQ+K
      IMR=JQ+K
      X=R(ILR)
      R(ILR)=R(IMR)
  180 R(IMR)=X
  185 CONTINUE
      RETURN
      END
	subroutine matcon(DM,DV)
C	MATCON: PROGRAM TO DECOMPOSE GENERAL LINEAR TRANSFORM INTO
C	ORTHOGONAL AND NON-ORTHOGONAL PARTS, THEN CONVERT
C	FORMER INTO EULER ROTATIONS, POLAR ROTATIONS, AND
C	AN EQUIVALENT SCREW OPERATION.
C
C	REFER TO:
C		M. ROSSMANN AND D. BLOW, ACTA CRYST. 15, P. 24 (1962)
C		R. DIAMOND, ACTA CRYST. A32, P. 1 (1976)
C		J. COX, J. MOL. BIOL. 28, P. 152 (1967)
C
C
C	THE PROGRAM WILL ALSO GENERATE A ROTATION MATRIX TO MAKE
C	THE EQUIVALENT SCREW AXIS COINCIDENT WITH THE Z-COORDINATE
C	AXIS ON REQUEST.
C
C
	DIMENSION DM(3,3),DV(3)
	DIMENSION DT(3,3), DTD(3,3), UI(3,3), UM(3,3),
     . RM(3,3), TM(3,3), R(3), TV(3), O(3)
	character*30 matfile
C
	dpr = 180./3.14159265
c
C	write(*,*) 'Matrix file name? (rtn.dat)'
C	read(*,2) matfile
C2	format(a)
C	if (matfile(1:1) .eq. ' ') matfile='rtn.dat'
C	open (unit=1,file=matfile)
C
C	read(1,1) ((DM(I,J), J=1,3),I=1,3)
C	read(1,1) (dv(j),j=1,3)
1	FORMAT(3F10.0)
	CALL DET (DM, D)
	write(*,*) '    Matrix and vector...'
	write(*,1012) ((DM(J,I),J=1,3),I=1,3),(DV(I),I=1,3)
1012	format(/(5x,3f10.4))
c	write(*,1001) D
1001	FORMAT(' Determinant is ',f5.2)
C
C	GET EULER (ROSSMAN & BLOW) ANGLES FROM DM
C       NOTE THIS IS THE TRANSPOSED DEFINITION per their errata 
C
	TEMP = SQRT(DM(3,1)*DM(3,1) + DM(3,2)*DM(3,2))
	THETA = ATAN2(TEMP,DM(3,3))
	PHI = ATAN2(DM(3,1),-DM(3,2))
	PSI = ATAN2(DM(1,3),DM(2,3))
	PHA = PHI*DPR
	THA = THETA*DPR
	PSA = PSI*DPR
	write(8,1009) PHA,THA,PSA
1009	FORMAT(/' R&B Euler angles (Z,X'',Z") ',3F8.1)
C
C	GET SCREW ROTATION AND DIRECTION COSINES, VECTOR
	ALFA = (DM(1,3)*DM(2,1) - (DM(1,1)-1.0)*DM(2,3))/
	1      (DM(2,3)*DM(1,2) - (DM(2,2)-1.0)*DM(1,3))
	BETA = (DM(1,2)*DM(3,1) - (DM(1,1)-1.0)*DM(3,2))/
	1      (DM(3,2)*DM(1,3) - (DM(3,3)-1.0)*DM(1,2))
	R(1) = SQRT(1.0/(1.0 + ALFA*ALFA + BETA*BETA))
	R(2) = R(1)*ALFA
	R(3) = R(1)*BETA
	COSD = (DM(1,1) + DM(2,2) + DM(3,3) - 1.0)/2.0
	SIND = (R(1)*(DM(2,3) - DM(3,2)) + R(2)*(DM(3,1) - DM(1,3))
	1       + R(3)*(DM(1,2) - DM(2,1)))/2.0
	DELTA = ATAN2(SIND,COSD)
	DA = DELTA*DPR
	TAND2 = 1.0E6
	IF(COSD .GT. -0.9999) TAND2 = SIND/(1.0 + COSD)
C	COMPONENT OF DV ALONG R (IN ANGSTROMS) AND POINT
C	ON ROTATION AXIS
	S = 0.0
	DO 35 I=1,3
	J = MOD(I,3) + 1
	K = MOD(I+1,3) + 1
	O(I) = 0.5*(DV(I)*(1.0 - R(I)*R(I)) - DV(J)*R(I)*R(J)
	1    - DV(K)*R(I)*R(K) + (R(K)*DV(J) - R(J)*DV(K))/TAND2)
35	S = S + R(I)*DV(I)
	write(*,1011) DA,R,O,S
	WRITE(8,1011) DA,R,O,S
1011	FORMAT(/' This corresponds to a rotation of',F8.2,
	1 ' degrees'/' around an axis with direction cosines:'
	2 3F8.4/' and passing through the point:'
	3 3F6.1/' The translation of the screw axis is',
	4 F6.1,' Ang')
C	COMBINE SCREW TRANSLATION WITH ORIGINAL TRANSLATION
	TEMP = 0.0
	DO 36 I=1,3
	TV(I) = R(I)*S+ DV(I)
36	TEMP = TEMP + TV(I)*TV(I)
	TEMP = SQRT(TEMP)
	write(8,1019) TEMP
1019	FORMAT(/' Net translation of original transform:',F8.1,
	1 ' Ang')
c	write(8,1012) TV
Cx
C	GET POLAR ANGLES (ROSSMAN & BLOW DEF.)
C
	TEMP = SQRT(R(1)*R(1) + R(3)*R(3))
	PSA = ATAN2(TEMP,R(2))*DPR
	PHA = ATAN2(-R(3),R(1))*DPR
	write(8,1020) PSA,PHA,DA
1020	FORMAT(/' R&B Polar angles (Psi,Phi,Chi) of rotation:'
	1 3F8.1)
C
C	GET ROTATIONS TO PLACE ROTATION AXIS ON Z AXIS
C	OF NEW COORD SYSTEM ... FOR PICTURES
C
	DO 40 I=1,3
	DO 39 J=1,3
	RM(I,J) = 0.0
39	TM(I,J) = 0.0
	RM(I,I) = 1.0
40	TM(I,I) = 1.0
	TEMP = SQRT(R(1)*R(1) + R(2)*R(2))
C	ROTATE AROUND Z AXIS TO BRING R INTO X-Z PLANE
	TM(1,1) = R(1)/TEMP
	TM(2,2) = R(1)/TEMP
	TM(1,2) = R(2)/TEMP
	TM(2,1) = -R(2)/TEMP
C	THEN ROTATE AROUND Y AXIS TO BRING R ONTO Z AXIS
	RM(1,1) = R(3)
	RM(3,3) = R(3)
	RM(1,3) = -TEMP
	RM(3,1) = TEMP
	CALL MM(RM,TM,DT)
	write(8,1015)
1015	FORMAT(/' Matrix to align screw axis with Z axis')
	write(8,1012) ((DT(I,J),J=1,3),I=1,3)
	END
	SUBROUTINE JACOBI(A,U,MAXIT)
C
C	BRUTE-FORCE JACOBI DIAGONALIZATION OF REAL,
C	SYMMETRIC MATRIX
C		INPUT MATRIX IS A
C		EIGENVALUES ARE RETURNED AS DIAGONALS OF A
C		EIGENVECTORS ARE RETURNED AS COLUMNS OF U
C Dimensions of A, U, UK and TM below are NxN. 
C N should be in call arguments
C
	DIMENSION A(3,3),U(3,3),UK(3,3),TM(3,3)
	N=3
	IF(MAXIT .LE. 0) MAXIT = 1000
C	INIT U, GET SUM OF ELEMENTS **2
	SUMA = 0.0
	DO 10 I=1,N
	DO 9 J=1,N
	TEMP = A(I,J)
	SUMA = SUMA + TEMP*TEMP
9	U(I,J) = 0.0
10	U(I,I) = 1.0
	EPS = SUMA/1.0E6
	NPASS = 0
C	CHECK SUM OF DIAGONALS **2 ... IS DIAGONAL?
20	NPASS = NPASS + 1
	SUMD = 0.0
	DO 25 I=1,N
	TEMP = A(I,I)
25	SUMD = SUMD + TEMP*TEMP
	IF(SUMA - SUMD .LT. EPS) GO TO 100	!DONE
	IF(NPASS .LE. MAXIT) GO TO 30
	write(*,1000) MAXIT
1000	FORMAT(' ***',I10,' ITERATIONS EXCEEDED IN'
	1,' JACOBI DIAGONALIZATION ROUTINE')
	GO TO 100
C	FIND LARGEST REMAINING NON-DIAGONAL ELEMENT
30	AMAX = 0.0
	DO 35 I=1,N
	DO 34 J=1,N
	IF(J .EQ. I) GO TO 34
	TEMP = ABS(A(I,J))
	IF(TEMP .LE. AMAX) GO TO 34
	AMAX = TEMP
	IM = I
	JM = J
34	CONTINUE
35	CONTINUE
	IF(AMAX .EQ. 0.0) GO TO 100
C	FIND ROTATION ANGLE
	TOP = 2.0*A(IM,JM)
	BOT = A(IM,IM) - A(JM,JM)
	TEMP = SIGN(1.0,BOT)
	P = TOP*TEMP
	Q = ABS(BOT)
	TEMP = SQRT(P*P + Q*Q)
	COSA = SQRT(0.5*(1.0 + Q/TEMP))
	SINA = P/(2.0*COSA*TEMP)
C	FORM UK
	DO 40 I=1,N
	DO 39 J=1,N
39	UK(I,J) = 0.0
40	UK(I,I) = 1.0
	UK(IM,IM) = COSA
	UK(IM,JM) = SINA
	UK(JM,JM) = -COSA
	UK(JM,IM) = SINA
C	APPLY TO U (NET TRANSFORMATION)
	DO 45 I=1,N
	DO 44 J=1,N
	TM(I,J) = 0.0
	DO 43 K=1,N
43	TM(I,J) = TM(I,J) + U(I,K)*UK(K,J)
44	CONTINUE
45	CONTINUE
	DO 50 I=1,N
	DO 50 J=1,N
50	U(I,J) = TM(I,J)
C	APPLY TO A
	DO 55 I=1,N
	DO 54 J=1,N
	TM(I,J) = 0.0
	DO 53 K=1,N
53	TM(I,J) = TM(I,J) + UK(I,K)*A(K,J)
54	CONTINUE
55	CONTINUE
	DO 60 I=1,N
	DO 59 J=1,N
	A(I,J) = 0.0
	DO 58 K=1,N
58	A(I,J) = A(I,J) + TM(I,K)*UK(K,J)
59	CONTINUE
60	CONTINUE
	GO TO 20
C
100	MAXIT = NPASS - 1
	RETURN
	END
C
	SUBROUTINE MT(A,B)
	DIMENSION A(3,3),B(3,3)
	DO 10 I=1,3
	DO 10 J=1,3
10	B(J,I) = A(I,J)
	RETURN
	END
C
	SUBROUTINE DET(A,D)
	DIMENSION A(3,3)
	D = A(1,1)*A(2,2)*A(3,3) + A(1,2)*A(2,3)*A(3,1)
	1 + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)
	2 - A(1,2)*A(2,1)*A(3,3) - A(1,1)*A(2,3)*A(3,2)
	RETURN
	END
	SUBROUTINE INVERT(A,B)
	DIMENSION IP(3),A(3,3),B(3,3)
	N = 3
	DO 102 K=1,N
	DO 101 J=1,N
101	B(J,K) = 0.0
	B(K,K) = 1.0
102	IP(K) = 0
	DO 150 I=1,N
	AMAX = 0.0
	DO 203 J=1,N
	IF(IP(J)) 203,104,203
104	DO 103 K=1,N
	IF(IP(K)) 103,106,103
106	IF(ABS(A(J,K)) - ABS(AMAX)) 103,103,107
107	AMAX = A(J,K)
	IRO = J
	ICO = K
103	CONTINUE
203	CONTINUE
	IP(ICO) = 1
	DO 108 L=1,N
	EX = A(IRO,L)
	A(IRO,L) = A(ICO,L)
	A(ICO,L) = EX/AMAX
	EX = B(IRO,L)
	B(IRO,L) = B(ICO,L)
108	B(ICO,L) = EX/AMAX
	A(ICO,ICO) = 1.0
	DO 109 LR=1,N
	IF(LR - ICO) 110,109,110
110	RED = A(LR,ICO)
	DO 111 LC=1,N
	A(LR,LC) = A(LR,LC) - A(ICO,LC)*RED
111	B(LR,LC) = B(LR,LC) - B(ICO,LC)*RED
	A(LR,ICO) = 0.0
109	CONTINUE
150	CONTINUE
	RETURN
	END
C
      SUBROUTINE MM(A,B,C)
      DIMENSION A(3,3),B(3,3),C(3,3)
      DO 1010 I=1,3
      DO 1010 J=1,3
      C(I,J)=A(I,1)*B(1,J)+A(I,2)*B(2,J)+A(I,3)*B(3,J)
 1010 CONTINUE
      RETURN
      END
C
      SUBROUTINE  MV(U,A,B)
      DIMENSION U(3,3),A(3),B(3)
       DO 1010 I=1,3
 1010 B(I)=U(I,1)*A(1)+U(I,2)*A(2)+U(I,3)*A(3)
      RETURN
      END
C
      FUNCTION VV(A,B)
      DIMENSION A(3),B(3)
      VV=A(1)*B(1)+A(2)*B(2)+A(3)*B(3)
      RETURN
      END
C
      SUBROUTINE UV(X)
      DIMENSION X(3)
      S = 0.0
      DO 10 I=1,3
10    S = S + X(I)*X(I)
      S = SQRT(S)
	IF(S.LE.0.) RETURN
      DO 20 I=1,3
20    X(I) = X(I)/S
      RETURN
      END
      SUBROUTINE VXV(A,B,C)
      DIMENSION A(3),B(3),C(3)
      C(1)=A(2)*B(3)-A(3)*B(2)
      C(2)=A(3)*B(1)-A(1)*B(3)
      C(3)=A(1)*B(2)-A(2)*B(1)
      RETURN
      END
      SUBROUTINE MATINV(AM,N,NFAIL)
      DIMENSION AM(21)
C
C                  SEGMENT 1 OF CHOLESKI INVERSION *****
C                  ***** FACTOR MATRIX INTO LOWER TRIANGLE X TRANSPOSE
      K=1
      IF (N-1) 10,8,9
    8 AM(1)=1.0/AM(1)
      GO TO 204
C                  ***** LOOP M OF A(L,M) *****
    9 DO 7 M=1,N
      IMAX=M-1
C                  ***** LOOP L OF A(L,M) *****
      DO 6 L=M,N
      SUMA=0.0
      KLI=L
      KMI=M
      IF (IMAX) 23,23,1
C                  ***** SUM OVER I=1,M-1 A(L,I)*A(M,I) *****
    1 DO 2 I=1,IMAX
      SUMA=SUMA+AM(KLI)*AM(KMI)
      J=N-I
      KLI=KLI+J
    2 KMI=KMI+J
C                  ***** TERM=C(L,M)-SUM *****
   23 TERM=AM(K)-SUMA
      IF (L-M) 3,3,5
    3 IF (TERM) 10,10,4
C                  ***** A(M,M)=SQRT(TERM) *****
    4 DENOM=SQRT(TERM)
      AM(K)=DENOM
      GO TO 6
C                  ***** A(L,M)=TERM/A(M,M) *****
    5 AM(K)=TERM/DENOM
    6 K=K+1
    7 CONTINUE
C
C                  SEGMENT 2 OF CHOLESKI INVERSION *****
C                  ***** INVERSION OF TRIANGULAR MATRIX *****
  100 AM(1)=1.0/AM(1)
      KDM=1
C                  ***** STEP L OF B(L,M) *****
      DO 104 L=2,N
      KDM=KDM+N-L+2
C                  ***** RECIPROCAL OF DIAGONAL TERM *****
      TERM=1.0/AM(KDM)
      AM(KDM)=TERM
      KMI=0
      KLI=L
      IMAX=L-1
C                  ***** STEP M OF B(L,M) *****
      DO 103 M=1,IMAX
      K=KLI
C                  ***** SUM TERMS *****
      SUMA=0.0
      DO 102 I=M,IMAX
      II=KMI+I
      SUMA=SUMA-AM(KLI)*AM(II)
  102 KLI=KLI+N-I
C                  ***** MULT SUM * RECIP OF DIAGONAL *****
      AM(K)=SUMA*TERM
      J=N-M
      KLI=K+J
  103 KMI=KMI+J
  104 CONTINUE
C
C                  SEGMENT 3 OF CHOLESKI INVERSION *****
C                  ***** PREMULTIPLY LOWER TRIANGLE BY TRANSPOSE *****
  200 K=1
      DO 203 M=1,N
      KLI=K
      DO 202 L=M,N
      KMI=K
      IMAX=N-L+1
      SUMA=0.0
      DO 201 I=1,IMAX
      SUMA=SUMA+AM(KLI)*AM(KMI)
      KLI=KLI+1
  201 KMI=KMI+1
      AM(K)=SUMA
  202 K=K+1
  203 CONTINUE
  204 NFAIL=0
  300 RETURN
   10 NFAIL=K
      RETURN
      END
      SUBROUTINE PCALC
C
C     EQUIVALENCE PROBABILITY CALCS
C
      LOGICAL NPL,NSL
      COMMON /PARAMS/ ICARD,IDISK,ITRAN,JSAVE,NPPRT,IMINR,IPROG,ISTAR,
     1ITMAX,NCYCLE,NEQ,IECRD,IHOLD,IDRAW,IORNT,NPL,DCUT,PCUT,E1,E2,
     2e1m,e2m,ITER,NI,NJ,NCHNG
      COMMON /COORDS/ XI(3,702),XJ(3,702)
      COMMON /SEQ/ IRES(702),JRES(702),AIRES(702),AJRES(702)
	COMMON /PIJ/ IJ(702,10)
      COMMON /SAVE/ NSV(702)
      COMMON /SCRAT/ JSV(10),PSV(10),D2SV(10),S2SV(10),DXSV(3,10)
      COMMON /XFRMS/ A(3,3),OA(3,3),OR(3),AS,AT,AP,V(3),ESIG,TSIG
      DIMENSION SUMSV(10),DOESV(10),SOESV(10),DSV(10),SSV(10),PNSV(10)
      DIMENSION IO(70),VT(3)
      EQUIVALENCE (X,VT(1)),(Y,VT(2)),(Z,VT(3))
      EQUIVALENCE (IO(1),JSV(1))
C
C     APPLY TRANSFORM TO TRIAL STRUCTURE -- ONLY IF LSCTE DID NOT
C
      IF(ITER .NE. ISTAR .OR. ITRAN .EQ. 0) GO TO 50
      N = NI + 1
      DO 45 I=1,N
      CALL MV(A,XI(1,I),VT)
      XI(1,I) = X + V(1)
      XI(2,I) = Y + V(2)
45    XI(3,I) = Z + V(3)
50    DC2 = DCUT*DCUT
      E12 = 2*E1*E1
      E22 = 2*E2*E2
      KK = JSAVE - 1
      OPEN (UNIT=9,NAME='pcl.tmp',status='unknown',FORM='UNFORMATTED')
      IF(NPL) GO TO 52
      WRITE(8,1015) JSAVE
1015  FORMAT(//' Equivalence probabilities:  ',I3,
	1' Most probable equivalences with reference structure',
     2' listed by atom in trial structure')
C
C --- OUTER LOOP:  TRIAL STRUCTURE ATOMS
C
52    REWIND 9
      DO 500 I=2,NI
      IJ(I,1) = 0
      DO 51 K=1,10
51    SUMSV(K) = 1.0E10
      NSAVE = 0
      X = XI(1,I)
      Y = XI(2,I)
      Z = XI(3,I)
      I1 = I - 1
      DXM = XI(1,I1) - X
      DYM = XI(2,I1) - Y
      DZM = XI(3,I1) - Z
      I1 = I + 1
      DXP = XI(1,I1) - X
      DYP = XI(2,I1) - Y
      DZP = XI(3,I1) - Z
C
C ---  INNER LOOP:  REFERENCE STRUCTURE ATOMS
C
      DO 400 J=2,NJ
      U = XJ(1,J)
      V1 = XJ(2,J)
      W = XJ(3,J)
      DX = U - X
      IF(DX .GE. DCUT) GO TO 400
      DY = V1 - Y
      DZ = W - Z
      D2 = DX*DX + DY*DY + DZ*DZ
      IF(D2 .GE. DC2) GO TO 400
C
C     SCATTER CALCS
C
      J1 = J - 1
      SXM = XJ(1,J1) - U
      SYM = XJ(2,J1) - V1
      SZM = XJ(3,J1) - W
      J1 = J + 1
      SXP = XJ(1,J1) - U
      SYP = XJ(2,J1) - V1
      SZP = XJ(3,J1) - W
C     "SP" SCATTER WILL APPLY IF A RUN HS JOFF = + 1
      T1 = DXM - SXM
      T2 = DYM - SYM
      T3 = DZM - SZM
      T4 = DXP - SXP
      T5 = DYP - SYP
      T6 = DZP - SZP
      SP2 = T1*T1 + T2*T2 + T3*T3 + T4*T4 + T5*T5 + T6*T6
C     "SN" SCATTER WILL APPLY IF A RUN HAS JOFF = -1
      T1 = DXM - SXP
      T2 = DYM - SYP
      T3 = DZM - SZP
      T4 = DXP - SXM
      T5 = DYP - SYM
      T6 = DZP - SZM
      SN2 = T1*T1 + T2*T2 + T3*T3 + T4*T4 + T5*T5 + T6*T6
      S2 = SP2
      PN = 1.0
      IF(SN2 .GE. SP2) GO TO 109
      S2 = SN2
      PN = -1.0
109   DOE = D2/E12
      SOE = S2/E22
      SUM = SOE + DOE
C     SAVE ONLY TOP JSAVE PROBABILITIES ... SMALLEST SUMS
C     SMALLER THAN LARGEST?  (SUMS IN SUMSV ARE IN ASCENDING ORDER)
      IF(SUM .GE. SUMSV(JSAVE)) GO TO 400
C     IF SO, FIND PLACE IN LIST
      K1 = JSAVE
      DO 110 K=1,KK
      IF(SUM .GE. SUMSV(K)) GO TO 110
      K1 = K
      GO TO 111
110   CONTINUE
111   IF(K .EQ. JSAVE) GO TO 105
C     IF PLACE IS IN MIDDLE OF LIST,EXPAND LIST FOR INSERTION
      L1 = JSAVE + 1
      L2 = JSAVE
      DO 100 L=K1,KK
      L1 = L1 - 1
      L2 = L2 - 1
      SUMSV(L1) = SUMSV(L2)
      DXSV(1,L1) = DXSV(1,L2)
      DXSV(2,L1) = DXSV(2,L2)
      DXSV(3,L1) = DXSV(3,L2)
      D2SV(L1) = D2SV(L2)
      S2SV(L1) = S2SV(L2)
      DOESV(L1) = DOESV(L2)
      SOESV(L1) = SOESV(L2)
      PNSV(L1) = PNSV(L2)
100   JSV(L1) = JSV(L2)
C     INSERT NEW SUM IN LIST
105   SUMSV(K1) = SUM
      DXSV(1,K1) = DX
      DXSV(2,K1) = DY
      DXSV(3,K1) = DZ
      D2SV(K1) = D2
      S2SV(K1) = S2
      DOESV(K1) = DOE
      SOESV(K1) = SOE
      PNSV(K1) = PN
      JSV(K1) = J
      NSAVE = NSAVE + 1
400   CONTINUE
      IF(NSAVE .GT. JSAVE) NSAVE = JSAVE
C
C      COMPUTE P(I,J)'S
C
      NSL = .TRUE.
      IF(NSAVE .EQ. 0) GO TO 451
      NSL = .FALSE.
      DO 450 K=1,NSAVE
C     PROTECT AGAINST EXP OVERFOLW ON LARGE DISTANCES
      PD = 0.0
      ARG = DOESV(K)
      IF(ARG .LT. 50.0) PD = EXP(-ARG)
C     PROTECT AGAINST EXPONENT OVERFLOW ON LARGE SCATTERS
      PS = 0.0
      ARG = SOESV(K)
      IF(ARG .LE. 50.0) PS = EXP(-ARG)
      P = PD*PS
      PSV(K) = P
450   CONTINUE
C
C     EITHER STORE GOODIES IN SCRATCH FILE & VARIOUS COMMON BLOCKS,
C     OR PRINT P(I,J)'S... BUT DON'T CHANGE ANYTHING ON PRINT RUN
C
451   IF(NPL) GO TO 490
C     PRINT ONLY
      WRITE(8,1018) IRES(I)
1018  FORMAT(/,7X,A4)
      IF(NSL) GO TO 460
      DO 455 K=1,NSAVE
      DSV(K) = SQRT(D2SV(K))
455   SSV(K) = SQRT(S2SV(K))
      WRITE(8,1019) (JRES(JSV(K)),K=1,NSAVE)
1019  FORMAT(13X,'JRES',4X,10(6X,A4))
      WRITE(8,1119) (PSV(K),K=1,NSAVE)
1119  FORMAT(12X,'P(I,J)',4X,10F10.5)
      WRITE(8,1219) (DSV(K),K=1,NSAVE)
1219  FORMAT(15X,'D',5X,10F10.2)
      WRITE(8,1319) (SSV(K),K=1,NSAVE)
1319  FORMAT(15X,'S',5X,10F10.2)
      GO TO 500
460   WRITE(8,1020)
1020  FORMAT(31X,'NO EQUIVALENT ATOMS')
      GO TO 500
C     STORE RESULTS
490   NSV(I) = NSAVE
      IF(NSAVE .EQ. 0) GO TO 500
      DO 495 K=1,NSAVE
C     NOTE THAT -J IS STORED IN IJ IF THE "SN" SCATTER WAS USED
      J = JSV(K)*PNSV(K)
      IF(PSV(K) .LT. PCUT) J = 0
      IJ(I,K) = J
495   CONTINUE
      WRITE(9) IO
500   CONTINUE
      CLOSE (UNIT=9)
      RETURN
      END
      SUBROUTINE RLIST
      LOGICAL NPL
      COMMON /PARAMS/ ICARD,IDISK,ITRAN,JSAVE,NPPRT,IMINR,IPROG,ISTAR,
     1ITMAX,NCYCLE,NEQ,IECRD,IHOLD,IDRAW,IORNT,NPL,DCUT,PCUT,E1,E2,
     2e1m,e2m,ITER,NI,NJ,NCHNG
      COMMON /SEQ/ IRES(702),JRES(702),AIRES(702),AJRES(702)
      COMMON /SAVE/ NSV(702)
      COMMON /DIREC/ NRUN,NSRUN,ISFIN(60),JSFIN(60),ISST(60),JSST(60),
     1JOF(60),ISTRT(60),JSTRT(60),IFIN(60),JFIN(60),SSPIJ(60),SPIJ(60),
     2SD2(60),SS2(60),SDX(60),SDY(60),SDZ(60),JCNT(60),IPRG1(60),
     4IPRG2(60),NPROG(60),ISPR(60),IEPR(60),INPR(60),NEXT(60),IGAP(60)
      COMMON /T99COM/ STRNG(64)
      DIMENSION LINE(20),KEY(10)
      EQUIVALENCE (STRNG(1),LIT1)
      EQUIVALENCE (LITR,XLITR)
      DATA IBLNK,ISTR /'    ','  * '/
C
C     SIDE-BY-SIDE LISTING OF RUNS AND EQUIVALENCES
C
      write(8,1044)
1044  FORMAT('1'//' Display of runs and equivalences ')
      WRITE(8,1045) 
1045  FORMAT(/'  Trial,     Reference, Reference, ...'/
     1,13X,'SIMPLE',5X,'EXTENDED',7X,'NO'/8x
     1,3(3X,'PROG RULE'),6X,37('-'),' RUNS ',37('-'))
      DO 954 J=1,10
      KEY(J) = 0
954   CONTINUE
      DO 955 K=1,NRUN
955   JCNT(K) = JSTRT(K) - JOF(K)
      NOTE = 0
      DO 970 I=2,NI
      IF(NSV(I) .EQ. 0) GO TO 969
      IF(NOTE .EQ. 1) WRITE(8,1040)
1040  FORMAT(//)
      NOTE = 0
      LSPR = IBLNK
      LSPRR = IBLNK
      LEPR = IBLNK
      LEPRR = IBLNK
      LNPR = IBLNK
      LNPRR = IBLNK
      DO 956 J=1,20
956   LINE(J) = IBLNK
      DO 965 K=1,NRUN
      IS = ISTRT(K)
      IF = IFIN(K)
      JS = JSTRT(K)
      JS1 = JS + JOF(K)
      IF((I .LT. IS .AND. I .LT. IF) .OR.
     1   (I .GT. IS .AND. I .GT. IF)) GO TO 965
      JCNT(K) = JCNT(K) + JOF(K)
      KK = JCNT(K)
      LIT = JRES(KK)
      XLITR=AJRES(KK)
C     IF LISTING OF RUN IS ALREADY STARTED, PLACE J IN SAME COL
      DO 958 J=1,10
      IF(KEY(J) .NE. K) GO TO 958
      J1 = 2*J
      LINE(J1) = LIT
      IF(KK .EQ. JS) GO TO 961
C     KEY IS SET TO -2 WHEN A RUN IS OVER ... THIS INSURES THAT RUNS
C     ARE SEPARATED BY SPACERS
      IF(KK .NE. JFIN(K)) GO TO 957
      KEY(J) = -2
C     OR SAVE COLUMN FOR FOLLOWING RUN IN REGISTER
      IF(IGAP(K) .NE. 0) KEY(J) = IGAP(K)
      GO TO 957
958   CONTINUE
C     IF STARTING TO LIST NEW RUN, FIND 1ST UNUSED COL ... IF NO FREE
C     SLOTS  (MAX 10), FORGET THIS RUN
      DO 960 J=1,10
      IF(KEY(J) .NE. 0) GO TO 960
      KEY(J) = K
      J1 = 2*J
      LINE(J1) = LIT
      GO TO 961
960   CONTINUE
      GO TO 965
C     LABEL RUN IF THIS IS FIRST RESIDUE
	j1=j1-1
961   continue
c961   ENCODE(4,1046,LIT1) K
c1046  FORMAT(I3,'-')
C     NOW SEE IF THIS RUN IS IN AN EQUIV
957   DO 964 K1=1,NRUN
      IF(K .NE. ISPR(K1)) GO TO 963
      IF(I .GE. ISST(K) .AND. I .LE. ISFIN(K)) THEN
         LSPR = LIT
         LSPRR = LITR
         END IF
963   IF(K .EQ. IEPR(K1)) THEN
         LEPR = LIT
         LEPRR = LITR
         END IF
      IF(K .EQ. INPR(K1)) THEN
         LNPR = LIT
         LNPRR = LITR
         END IF
964   CONTINUE
      IF(KK .NE. JS1) GO TO 965
C     LIST PROB IF THIS IS SECOND RESIDUE OF RUN
	j1 = j1-1
c      ENCODE(4,1047,LIT1) SPIJ(K)
c1047  FORMAT(F4.1)
965   CONTINUE
      DO 966 J=1,10
      IF(KEY(J)) 967,966,968
967   KEY(J) = KEY(J) + 1
      GO TO 966
968   J1 = 2*J
      IF(LINE(J1) .EQ. IBLNK) LINE(J1) = ISTR
966   CONTINUE
      WRITE(8,1043)
     .  AIRES(I),IRES(I),LSPRR,LSPR,LEPRR,LEPR,LNPRR,LNPR,LINE
1043  FORMAT(2X,A4,1X,A4,2X,3(A4,1X,A4,2X),20A4)
      GO TO 970
969   NOTE = 1
970   CONTINUE
      RETURN
      END
      SUBROUTINE RUNS
C *** TEMP KLUDGE FOR E-SCANS: PRINTED SUMMARY OF RUNS TO UNIT 7
      LOGICAL NPL
      COMMON /PARAMS/ ICARD,IDISK,ITRAN,JSAVE,NPPRT,IMINR,IPROG,ISTAR,
     1ITMAX,NCYCLE,NEQ,IECRD,IHOLD,IDRAW,IORNT,NPL,DCUT,PCUT,E1,E2,
     2e1m,e2m,ITER,NI,NJ,NCHNG
      COMMON /SEQ/ IRES(702),JRES(702),AIRES(702),AJRES(702)
	COMMON /PIJ/ IJ(702,10)
      COMMON /SAVE/ NSV(702)
      COMMON /SCRAT/ JSV(10),PSV(10),D2SV(10),S2SV(10),DXSV(3,10)
      COMMON /DIREC/ NRUN,NSRUN,ISFIN(60),JSFIN(60),ISST(60),JSST(60),
     1JOF(60),ISTRT(60),JSTRT(60),IFIN(60),JFIN(60),SSPIJ(60),SPIJ(60),
     2SD2(60),SS2(60),SDX(60),SDY(60),SDZ(60),JCNT(60),IPRG1(60),
     4IPRG2(60),NPROG(60),ISPR(60),IEPR(60),INPR(60),NEXT(60),IGAP(60)
      DIMENSION IO(70),JMP(702)
      EQUIVALENCE (IO(1),JSV(1)),(JMP(1),IJ(1,1))
C
C     MAXIMUM ALLOWED GAP (FOR RUNS IN REGISTER) ... POSS SHOULD BE READ
      IMAXG = 3
      OPEN (UNIT=9,NAME='pcl.tmp',status='OLD',FORM='UNFORMATTED')
C
C     LOCATE SIMPLE RUNS (I.E., RUNS IN MOST PROB EQUIVS)
C
      JOFF = 0
      NRUN = 0
      DO 750 I=2,NI
      I1 = I + 1
      J = JMP(I)
      J1 = JMP(I1)
      IF(JOFF .EQ. 0) GO TO 720
      IF(J1 - J .EQ. 1) GO TO 750
      ISFIN(NRUN) = I
      JSFIN(NRUN) = J
      GO TO 740
720   IF(J .EQ. 0 .OR. J1 .EQ. 0) GO TO 750
      JOFF= J1 - J
C     NOTE THAT INVERSE RUNS ARE ONLY ALLOWED FOR RESIDUES THAT USED
C     'REVERSE' SCATTER, WHICH ARE NEGATED IN IJ ... THUS JOFF MUST = +1
      IF(JOFF .NE. 1) GO TO 740
      NRUN = NRUN + 1
      ISST(NRUN) = I
      JSST(NRUN) = J
      JOF(NRUN) = ISIGN(1,J)
C     THE EXHAUSTIVE RUN ROUTINE WILL BOMB IF IT CAN'T PUT IN THE FIRST
C     RUN IT FINDS ... SO STOP THE SIMPLE RUNS 1 EARLY
      IF(NRUN .LT. 60) GO TO 750
      NRUN = NRUN - 1
      WRITE(8,1034)
      GO TO 745
740   JOFF = 0
750   CONTINUE
      IF(NRUN .EQ. 0) GO TO 796
745   NSRUN = NRUN
C
C     SKIP EXTENSIONS AND EXHAUSTIVE RUN SEARCH IF JSAVE = 1
C
      IF(JSAVE .GT. 1) GO TO 751
      DO 758 K=1,NRUN
      ISTRT(K) = ISST(K)
      IFIN(K) = ISFIN(K)
      JSTRT(K) = JSST(K)
758   JFIN(K) = JSFIN(K)
      GO TO 797
C
C     EXTEND RUNS
C
751   DO 790 K=1,NRUN
C     EXTEND BELOW
      I = ISST(K)
      J = JSST(K)
760   I = I - 1
      IF(I .LT. 2) GO TO 775
      NSAVE = NSV(I)
      IF(NSAVE .EQ. 0) GO TO 775
      DO 770 L=1,NSAVE
      IF(J - IJ(I,L) .NE. 1) GO TO 770
      J = J - 1
      GO TO 760
770   CONTINUE
775   ISTRT(K) = I + 1
      JSTRT(K) = J
C     EXTEND ABOVE
      I = ISFIN(K)
      J = JSFIN(K)
780   I = I + 1
      IF(I .GT. NI) GO TO 785
      NSAVE = NSV(I)
      IF(NSAVE .EQ. 0) GO TO 785
      DO 781 L=1,NSAVE
      IF(IJ(I,L) - J .NE. 1) GO TO 781
      J = J + 1
      GO TO 780
781   CONTINUE
785   IFIN(K) = I - 1
      JFIN(K) = J
790   CONTINUE
C
C     ELIMINATE DUPLICATE RUNS
C
C     (WHICH OCCUR WHEN THE SAME EXTENDED RUN CAN BE FORMED FROM TWO
C     UNCONNECTED SIMPLE RUNS)
792   NR1 = NRUN - 1
      NR2 = NRUN
      DO 795 K=1,NR1
      IS = ISTRT(K)
      IF= IFIN(K)
      JS = JSTRT(K)
      JF = JFIN(K)
      K1 = K + 1
      DO 794 L=K1,NR2
      IF(ISTRT(L) .NE. IS .OR. IFIN(L) .NE. IF .OR.
     1   JSTRT(L) .NE. JS .OR. JFIN(L) .NE. JF) GO TO 794
      DO 793 L1=L,NR1
      L2 = L1 + 1
      ISST(L1) = ISST(L2)
      ISFIN(L1) = ISFIN(L2)
      JSST(L1) = JSST(L2)
      JSFIN(L1) = JSFIN(L2)
      ISTRT(L1) = ISTRT(L2)
      IFIN(L1) = IFIN(L2)
      JSTRT(L1) = JSTRT(L2)
      JFIN(L1) = JFIN(L2)
793   JOF(L1) = JOF(L2)
      NRUN = NRUN - 1
      IF(K1 .LE. NSRUN) NSRUN = NSRUN - 1
      GO TO 792
794   CONTINUE
795   CONTINUE
C
C     LOCATE RUNS EXHAUSTIVELY
C
796   N = NI - 1
      DO 950 I=2,N
      I1 = I + 1
      NUM1 = NSV(I)
      NUM2 = NSV(I1)
      IF(NUM1 .EQ. 0 .OR. NUM2 .EQ. 0) GO TO 950
      DO 949 J=1,NUM1
      NUM3 = IJ(I,J)
      DO 948 J1=1,NUM2
      JOFF = IJ(I1,J1) - NUM3
      IF(JOFF .NE. 1) GO TO 948
C     THIS NUCLEATES A 'GENERAL' RUN ... NOW EXTEND BELOW
      JSIGN = ISIGN(1,NUM3)
      L1 = I
      L2 = NUM3
910   L1 = L1 - 1
      IF(L1 .LT. 2) GO TO 915
      NSAVE = NSV(L1)
      IF(NSAVE .EQ. 0) GO TO 915
      DO 914 L=1,NSAVE
      IF(L2 - IJ(L1,L) .NE. 1) GO TO 914
      L2 = L2 - 1
      GO TO 910
914   CONTINUE
915   IS = L1 + 1
      JS = L2
C     AND ABOVE
      L1 = I1
      L2 = NUM3 + 1
920   L1 = L1 + 1
      IF(L1 .GT. NI) GO TO 925
      NSAVE= NSV(L1)
      IF(NSAVE .EQ. 0) GO TO 925
      DO 924 L=1,NSAVE
      IF(IJ(L1,L) - L2 .NE. 1) GO TO 924
      L2 = L2 + 1
      GO TO 920
924   CONTINUE
925   IF = L1 - 1
      JF = L2
C     NOW TRY TO ELIMINATE IT
      DO 947 K=1,NRUN
      IF(IS .EQ. ISTRT(K) .AND. IF .EQ. IFIN(K) .AND.
     1   JS .EQ. JSTRT(K) .AND. JF .EQ. JFIN(K)) GO TO 948
947   CONTINUE
      NUM = IF - IS + 1
      IF(NUM .LT. IMINR) GO TO 933
      DO 930 K=1,NRUN
      ISK = ISTRT(K)
      IFK = IFIN(K)
      JSK = JSTRT(K)
      JFK = JFIN(K)
      IF((IS .LT. ISK .AND. IF .LT. ISK) .OR.
     1   (IS .GT. IFK .AND. IF .GT. IFK)) GO TO 930
      IF((JS .LT. JSK .AND. JF .LT. JSK) .OR.
     1   (JS .GT. JFK .AND. JF .GT. JFK)) GO TO 930
      IF(NUM .GT. IFK - ISK + 1) GO TO 930
      GO TO 933
930   CONTINUE
932   NRUN = NRUN + 1
      ISTRT(NRUN) = IS
      IFIN(NRUN) = IF
      JSTRT(NRUN) = JS
      JFIN(NRUN) = JF
      JOF(NRUN) = JSIGN
      ISST(NRUN) = 0
      ISFIN(NRUN) = 0
      IF(NRUN .LT. 60) GO TO 948
      WRITE(8,1034)
1034  FORMAT(//' *** TOO MANY RUNS')
      GO TO 797
C     FINALLY, TRY TO CONCATENATE TO AN EXISTING RUN
933   DO 929 K1=1,NRUN
      IF(JSIGN .NE. JOF(K1)) GO TO 929
      IF(IS .GT. ISTRT(K1)) GO TO 922
      NGAP = ISTRT(K1) - IF
      IF(NGAP .NE. JSTRT(K1) - JF) GO TO 929
      GO TO 923
922   NGAP = IS - IFIN(K1)
      IF(NGAP .NE. JS - JFIN(K1)) GO TO 929
923   IF(NGAP .GT. 1 .AND. NGAP .LE. IMAXG) GO TO 932
929   CONTINUE
948   CONTINUE
949   CONTINUE
950   CONTINUE
      IF(NRUN .NE. 0) GO TO 797
      WRITE(8,1037) ITER
1037  FORMAT(//' *** NO RUNS IN ITERATION ',I2,5X,'-- I QUIT')
      STOP
C
C     MIX AND MATCH; ELIM SMALL END-OVERLAPS
C
C     1. MIX AND MATCH
C     FOR MODERATELY SIMILAR STRUCTURES, WE TREAT RUNS AS UNITS
C     OF SECONDARY STRUCTURE, RATHER THAN "CROSSING OVER" FREELY
C     AMONG RUNS AS HOMOMGR DOES; HOWEVER, FOR REGIONS OF HIGH
C     SIMILARITY, SMALL INSERTIONS OR DELETIONS CAUSE ARTIFACTS SINCE
C     THE UNMATCHED ENDS CAN GENERALLY BE EXTENDED TO LOWER PIJ'S,
C     ESPECIALLY WHEN E1 AND E2 ARE HIGH.  HERE WE ATTEMPT TO
C     ELIMINATE INAPPROPRIATE EXTENSIONS FOR RUNS WHOSE REF SEQUENCES
C     ARE IN REGISTER WITH ONLY SMALL GAPS OR OMISSIONS...  SINCE
C     WE ONLY WISH TO APPLY THIS CORRECTION TO REGIONS THAT ARE VERY
C     SIMILAR, WE CONSIDER ONLY ABUTTING SIMPLE RUNS; NOTE THOUGH THAT
C     WE TRUNCATE RUNS *AFTER* THE EXHAUSTIVE SEARCH TO ALLOW THE FULLY
C     EXTENDED RUNS TO SCREEN OUT SECONDARY RUNS.
C     2. ELIMINATE IRRITATING 1- OR 2-RES OVERLAPS
C     A SIMILAR PROBLEM OCCURS FOR GENERAL RUNS (I.E., NOT
C     NECESSARLY ABUTTING SIMPLE RUNS), BUT IN PRACTICE
C     INVOLVES ONLY ONE (OR AT HIGH E'S) TWO OVERLAPPING
C     RESIDUES AT THE ENDS OF ADJACENT RUNS... HERE WE ALSO
C     TRIM THESE OFF.  TO INSURE AGAIN THAT WE ONLY SPLICE IN SIMILAR
C     REGIONS, WE REQUIRE THAT THE RUNS GO IN THE SAME DIRECTION.
C     (BESIDES, SPLICE CAN ONLY DEAL WITH SUCH RUNS.)
C
C     COUNT THE NUMBER OF SPLICES IN SIMPLE RUNS (NSS) AND
C     IN GEN. RUNS (NGS) SEPARATELY; SET AD HOC CONSTANTS
C     FOR THE MIN LENGTH OF RUN TO FOOL WITH AND THE MAX
C     OVERLAP ALLOWED AT ENDS OF GEN RUNS (THIS IS ONE LARGER
C     THAN THE LARGEST NUICANCE OVERLAP I'VE SEEN SO FAR).
797   NSS = 0
      NGS = 0
      MINSPL = 5
      MAXOVR = 2
      NR1 = NRUN - 1
C     POINTER TO NEXT RECORD IN SEARCH FILE
      I1 = 2
      DO 975 K=1,NR1
C     SKIP LITTLE ONES
      IF(IFIN(K)-ISTRT(K) .LT. MINSPL) GO TO 975
      L = K + 1
      DO 974 K1=L,NRUN
      IF(JOF(K) .NE. JOF(K1)) GO TO 974
      IF(IFIN(K1)-ISTRT(K1) .LT. MINSPL) GO TO 974
C     SPLICE ABUTTING SIMPLE RUNS (WHICH WILL BE IN ORDER ON I)
      IF(K1 .GT. NSRUN) GO TO 960
      IF(ISFIN(K)+1 .NE. ISST(K1)) GO TO 960
      CALL SPLICE(K,K1,I1,2,KBAD,NSS)
      GO TO 974
C     ELIM DANGLING ENDS: EITHER K OR K1 MAY BE FIRST
960   IF(ISTRT(K1) .LT. ISTRT(K)) GO TO 965
      ITMP = IFIN(K) - ISTRT(K1)
      IF(ITMP .LT. 0) GO TO 961
      IF(ITMP .LE. MAXOVR) GO TO 884
      GO TO 974
961   JTMP = JFIN(K) - JSTRT(K1)
      IF(JTMP .LT. 0) GO TO 974
      IF(JTMP .GT. MAXOVR) GO TO 974
884   CALL SPLICE(K,K1,I1,2,KBAD,NGS)
      GO TO 974
965   ITMP = IFIN(K1) - ISTRT(K)
      IF(ITMP .LT. 0) GO TO 966
      IF(ITMP .LE. MAXOVR) GO TO 885
      GO TO 974
966   JTMP = JFIN(K1) - JSTRT(K)
      IF(JTMP .LT. 0) GO TO 974
      IF(JTMP .GT. MAXOVR) GO TO 974
885   CALL SPLICE(K1,K,I1,2,KBAD,NGS)
974   CONTINUE
975   CONTINUE
C
C     CORRECT RESIDUE NUMBERS OF REVERSE RUNS
C
      DO 958 K=1,NRUN
      IF(JOF(K) .EQ. 1) GO TO 958
      JSST(K) = -JSST(K)
      JSTRT(K) = -JSTRT(K)
      JSFIN(K) = -JSFIN(K)
      JFIN(K) = -JFIN(K)
958   CONTINUE
C
C     ACCUMULATE RUN STATISTICS
C
951   DO 791 K=1,NRUN
      IGAP(K) = 0
      NEXT(K) = 0
      SPIJ(K) = 0.0
      SD2(K) = 0.0
      SS2(K) = 0.0
      SDX(K) = 0.0
      SDY(K) = 0.0
      SDZ(K) = 0.0
      SSPIJ(K) = 0.0
791   JCNT(K) = JSTRT(K) - JOF(K)
      REWIND 9
      DO 800 I=2,NI
      NSAVE = NSV(I)
      IF(NSAVE .EQ. 0) GO TO 800
      READ(9) IO
      DO 799 K=1,NRUN
      IF(I .LT. ISTRT(K) .OR. I .GT. IFIN(K)) GO TO 799
      J1 = JCNT(K) + JOF(K)
      JCNT(K) = J1
      DO 798 J=1,NSAVE
      IF(JSV(J) .NE. J1) GO TO 798
      SPIJ(K) = SPIJ(K) + PSV(J)
      SD2(K) = SD2(K) + D2SV(J)
      SS2(K) = SS2(K) + S2SV(J)
      SDX(K) = SDX(K) + DXSV(1,J)
      SDY(K) = SDY(K) + DXSV(2,J)
      SDZ(K) = SDZ(K) + DXSV(3,J)
      IF(I .LT. ISST(K) .OR. I .GT. ISFIN(K)) GO TO 798
      SSPIJ(K) = SSPIJ(K) + PSV(J)
798   CONTINUE
799   CONTINUE
800   CONTINUE
      CLOSE (UNIT=9)
C
C     LIST RUN DIRECTORY AND STATS
C
      WRITE(8,1021) ITER
1021  FORMAT(//5x,'Summary of runs for iteration',I3)
      WRITE(8,1026)
1026  FORMAT(//5X,'RUN',3X,'Trial limits',5X,'Ref limits',4X,
     1'Residues',3X,'rms D',5X,'rms S',4X,'sum P(I,J)  Ssum P(I,J)',
     24X,'<DX>  <DY>  <DZ>   <D>')
      DO 810 K=1,NRUN
      IF = IFIN(K)
      IS = ISTRT(K)
      JF = JFIN(K)
      JS = JSTRT(K)
      NUM = IF - IS + 1
      D = SQRT(SD2(K)/NUM)
      S = SQRT(SS2(K)/NUM)
      DX = SDX(K)/NUM
      DY = SDY(K)/NUM
      DZ = SDZ(K)/NUM
      AD = SQRT(DX*DX + DY*DY + DZ*DZ)
      WRITE(8,1022) K,IRES(IS),IRES(IF),JRES(JS),JRES(JF),NUM,D,S,
     1SPIJ(K),SSPIJ(K),DX,DY,DZ,AD
1022  FORMAT(3X,I5,2(4X,A4,'-- ',A4),I9,1X,2F10.2,3X,F10.5,2X,F10.5,
     14X,4F6.2)
810   CONTINUE
C
C     LOOK FOR RUNS IN REGISTER AND NOTE FOR USER
C
      WRITE(8,1036)
1036  FORMAT(/5X,'Runs in register',4X,'Gap',4X,
     1'Combined prob',4X,'Combined residues')
      NR1 = NRUN - 1
      DO 815 K=1,NR1
      JOFF = JOF(K)
      IS = ISTRT(K)
      IF = IFIN(K)
      JS = JSTRT(K)
      JF = JFIN(K)
      L = K + 1
      DO 814 K1=L,NRUN
      IF(JOFF .NE. JOF(K1)) GO TO 814
      MOD = 1
      IF(IS .GT. ISTRT(K1)) MOD = MOD + 1
      IF(JOFF .EQ. -1) MOD = MOD + 2
      GO TO (808,809,808,809),MOD
808   NGAP = ISTRT(K1) - IF
      IF(NGAP .NE. (JSTRT(K1) - JF)*JOFF) GO TO 814
      GO TO 811
809   NGAP = IS - IFIN(K1)
      IF(NGAP .NE. (JS - JFIN(K1))*JOFF) GO TO 814
811   IF(NGAP .LE. 1 .OR. NGAP .GT. IMAXG) GO TO 814
      NUM = NGAP - 1
      P = SPIJ(K) + SPIJ(K1)
      NUM1 = IFIN(K) - ISTRT(K) + IFIN(K1) - ISTRT(K1) + 2
      WRITE(8,1035) K,K1,NUM,P,NUM1
1035  FORMAT(5X,I4,' and',I4,5X,I5,8X,F6.3,12X,I5)
      IF(NEXT(K) .NE. 0) GO TO 814
      NEXT(K) = K1
      IF(IS .GT. IFIN(K1)) GO TO 817
      IGAP(K) = K1
      GO TO 814
817   IGAP(K1) = K
814   CONTINUE
815   CONTINUE
D     WRITE(8,2000) NSS,NGS
D2000 FORMAT(' NSS =',I5,3X,'NGS =',I5)
      RETURN
      END
      SUBROUTINE SELEQ
C
C     SELECT BEST EQUIVALENCES -- REVISED 4.80 (WSB)
C
C     NEW SELECTION ALGORITHMS:
C       IN EARLIER VERSIONS OF SELEQ, THE SELECTION ALGORITHMS WERE
C     ORDER-DEPENDENT.  IN PARTICULAR, A GIVEN RUN (E.G., I) MIGHT
C     BE USED TO SCREEN OUT SEVERAL RUNS AND THEN BE ELIMINATED
C     ITSELF BY ANOTHER RUN (E.G., J); THEN RUNS FOLLOWING J IN
C     THE DIRECTORY WOULD NEVER BE COMPARED TO RUN I.
C       IN THE PRESENT VERSION, WE USE TWO LIMIT-CASE ALGORITHMS THAT
C     ARE NOT ORDER-DEPENDENT:
C       (1) A "MORE-PERMISSIVE ALGORITHM" IN WHICH A RUN IS ELIMINATED
C     ONLY IF IT CONFLICTS WITH ANOTHER RUN OF HIGHER P(I,J) THAT IS
C     ALREADY IN THE EQUIVALENCE.  THIS IS INSURED BY CONSIDERING THE
C     RUNS IN DECREASING ORDER OF P(I,J); A RUN IS NOT USED FOR SCREEN-
C     ING UNLESS IT IS COMPATIBLE WITH ALL RUNS OF HIGHER P(I,J) IN THE
C     EQUIVALENCE.
C       THIS ALGORITHM IS USED FOR THE SIMPLE-PROG-RULE SELECTION AND
C     SHOULD GIVE RESULTS SIMILAR TO THOSE OBTAINED BY ROSSMANN &
C     ARGOS.  THEY USE A SELECTION PROCEDURE INVOLVING SUCCESSIVE PAIR-
C     WISE COMPARISONS OF RUNS IN SEQUENCE; THUS THE PROCEDURE IS
C     PERMISSIVE IN THE SENSE GIVEN ABOVE, BUT ORDER-DEPENDENT.
C       (2) A "LESS-PERMISSIVE ALGORITHM" IN WHICH A RUN IS ELIMINATED
C     IF IT CONFLICTS WITH *ANY* RUN OF HIGHER P(I,J), WHETHER OR NOT
C     THE RUN IS IN THE EQUIVALENCE; I.E., A RUN MUST BE COMPATIBLE
C     NOT ONLY WITH OTHER RUNS IN THE EQUIV, BUT ALSO WITH RUNS OF
C     HIGHER P(I,J) THAT WERE PREVIOUSLY EXCLUDED FROM THE EQUIV.
C       THIS ALGORITHM IS USED FOR THE EXTENDED-PROG-RULE AND NO-PROG-
C     RULE EQUIVS.
C       --NOTE THAT WE MUST NOW DISTINGUISH BETWEEN RUNS THAT *CANNOT*
C     BE IN THE EQUIV (I.E., REVERSE RUNS CANNOT, BY DEFINITION, APPEAR
C     IN A TOPOLOGICAL EQUIV), AND RUNS THAT COULD BE IN THE EQUIV BUT
C     CONFLICT WITH OTHER SUCH RUNS.
C
C     --ALSO NEW IN THIS VERSION:  THE TRIMMING OF OVERLAPPING ENDS
C     OF ADJACENT RUNS (USED HERE TO EMULATE THE ROSSMANN-ARGOS
C     PROCEDURE AND IN "RUNS" TO MAKE AD-HOC CORRECTIONS) HAVE BEEN
C     CONSOLIDATED IN A SINGLE SUBROUTINE "SPLICE", WHICH DOES ALL
C     THE RE-SCANNING OF THE P(I,J) TABLES.  NOTE THAT SPLICE
C     REQUIRES THE CALLING MODULE TO OPEN AND CLOSE THE SCRATCH FILE
C     CONTAINING THE P(I,J) TABLES; ELSEWHERE, ALL FILES EXCEPT
C     SYS$INPUT AND SYS$OUTPUT (IN VAX-ESE) ARE OPENED AND CLOSED
C     BY EACH SUBROUTINE THAT USES THEM.
C
C
      LOGICAL NPL
      COMMON /PARAMS/ ICARD,IDISK,ITRAN,JSAVE,NPPRT,IMINR,IPROG,ISTAR,
     1ITMAX,NCYCLE,NEQ,IECRD,IHOLD,IDRAW,IORNT,NPL,DCUT,PCUT,E1,E2,
     2e1m,e2m,ITER,NI,NJ,NCHNG
      COMMON /SEQ/ IRES(702),JRES(702),AIRES(702),AJRES(702)
	COMMON /EQUIV/ IEQ(702),IQND(2,3)
      COMMON /SAVE/ NSV(702)
      COMMON /DIREC/ NRUN,NSRUN,ISFIN(60),JSFIN(60),ISST(60),JSST(60),
     1JOF(60),ISTRT(60),JSTRT(60),IFIN(60),JFIN(60),SSPIJ(60),SPIJ(60),
     2SD2(60),SS2(60),SDX(60),SDY(60),SDZ(60),JCNT(60),IPRG1(60),
     4IPRG2(60),NPROG(60),ISPR(60),IEPR(60),INPR(60),NEXT(60),IGAP(60)
      DIMENSION INEQ(16),IOLD(702),IRANK(60)
      DATA IBL/'    '/
C
C     START WITH ALL POSSIBLE RUN IN EACH EQUIVALENCE.  FIRST SCREEN
C     OUT OBVIOUSLY INELIGIBLE RUNS (DIRECTION AND LENGTH CRITERIA),
C     SETTING THESE RUNS TO -1 IN THE EQUIVALENCE ARRAY.
C     THEN SCREEN BY PROGRESSION RULE, AS NEEDED, AND END-OVERLAP
C     CRITERIA FOR EACH TYPE OF RUN BELOW.
C
      DO 821 K=1,NRUN
      IRANK(K) = 0
      ISPR(K) = 0
      IEPR(K) = 0
      INPR(K) = 0
      IPRG1(K) = 1
      IF(JOF(K) .EQ. -1) IPRG1(K) = -1
      IPRG2(K) = 1
      IF(JOF(K) .EQ. -1) IPRG2(K) = -1
      NPROG(K) = 1
      IF(IFIN(K)-ISTRT(K)+1 .GE. IMINR) GO TO 821
      IPRG2(K) = -2
      NPROG(K) = -2
821   CONTINUE
      NR1 = NRUN - 1
      NS1 = NSRUN - 1
      WRITE(8,1027) ITER
1027  FORMAT(//5x,'Equivalences for interation',I3)
C
C     WITH PROGRESSION RULE APPLIED TO SIMPLE RUNS (ROSSMAN'S WAY)
C     ...NOTE THAT A LOT OF THE CODE HERE IS INTENDED TO EMULATE
C     M. G. ROSSMANN'S PROGRAM "HOMOMGR" (HIS NAME, NO SLUR INTENDED),
C     WHICH FILLS UP A SINGLE ARRAY OF EQUIV RESIDUES WITH THE
C     MOST PROBABLE RESIDUE AVAILABLE; THE SIMPLE PROG RULE RESULTS
C     SHOULD GIVE SOMETHING VERY SIMILAR TO THE EQUIV HOMOMGR WOULD FIND
C
C     FIRST APPLY PROG RULE TO SIMPLE RUNS
      WRITE(8,1028)
1028  FORMAT(//5X,'Best equivalence using progression rule',
     1' on simple runs'/)
C     SORT SIMPLE RUNS INTO DECREASING PIJ IN ARRAY IRANK(I)
C     NOTE IRANK INIT ALL 0 ABOVE
      IRANK(1) = 1
      K1 = 1
      DO 50 K=2,NSRUN
      SSPIJK = SSPIJ(K)
      DO 40 L=1,NSRUN
      IF(L .GT. K1) GO TO 30
      M = IRANK(L)
      IF(SSPIJK .LE. SSPIJ(M)) GO TO 40
C     HAVE FOUND POSITION OF RUN K AND IT'S NOT ON END; SO
C     EXPAND ARRAY TO INSERT AT L
      M = K1
20    IRANK(M+1) = IRANK(M)
      M = M - 1
      IF(M .LT. L) GO TO 30
      GO TO 20
30    IRANK(L) = K
      K1 = K1 + 1
      GO TO 50
40    CONTINUE
50    CONTINUE
D     write(8,2000)
D2000 FORMAT(' IRANK:')
D     write(8,2001) (IRANK(K),K=1,NSRUN)
D2001 FORMAT(1X,44I3)
D     write(8,2002)
D2002 FORMAT(' PASS    R1 ELIMS R2    IST,FI R1    IST,FI R2',
D    1 4X,'JST,FI R1    JST,FI R2')
C
C     NOW SCREEN SIMPLE RUNS EXHAUSTIVELY FOR VIOLATIONS OF
C     PROGRESSION RULE; NOTE THAT ONLY RUNS IN CONFLICT WITH
C     A RUN OF HIGHER PIJ ARE OMITTED HERE, WHILE RUNS THAT
C     MERELY OVERLAP (I.E., EXTEND INTO BUT NOT BEYOND THE
C     SECOND RUN) ARE KEPT TO BE TRIMMED BY THE CROSSOVER
C     ROUTINE BELOW.
      DO 830 K=1,NSRUN
      KR = IRANK(K)
      IF(IPRG1(KR) .LE. 0) GO TO 830
C     RUNS MARKED -2 IN IPRG1 ARE EITHER THE CURRENT TEST RUN (KR)
C     OR A TEST RUN THAT HAS SURVIVED THE INNER LOOP; SINCE ITS
C     PIJ IS LARGER THAN SUBSEQUENT TEST RUNS, NO NEED TO RECHECK IT.
      IPRG1(KR) = -2
      DO 829 K1=1,NSRUN
      IF(IPRG1(K1) .LE. 0) GO TO 829
      IF(ISFIN(KR) .GE. ISFIN(K1)) GO TO 823
      IF(ISST(KR) .GE. ISST(K1)) GO TO 824
      IF(JSFIN(KR) .GE. JSFIN(K1)) GO TO 824
      IF(JSST(KR) .GE. JSST(K1)) GO TO 824
      GO TO 829
823   IF(ISST(K1) .GE. ISST(KR)) GO TO 824
      IF(JSFIN(K1) .GE. JSFIN(KR)) GO TO 824
      IF(JSST(K1) .LT. JSST(KR)) GO TO 829
824   IF(SSPIJ(KR) - SSPIJ(K1)) 826,827,828
826   IPRG1(KR) = 0
D     write(8,2003) K,K1,KR,ISST(K1),ISFIN(K1),ISST(KR),ISFIN(KR),
D    1 JSST(K1),JSFIN(K1),JSST(KR),JSFIN(KR)
D2003 FORMAT(2X,I2,5X,I2,7X,I2,4(5X,I3,'-',I3,1X))
      GO TO 830
827   WRITE(8,1025) KR,K1
828   IPRG1(K1) = 0
D     write(8,2003) K,KR,K1,ISST(KR),ISFIN(KR),ISST(K1),ISFIN(K1),
D    1 JSST(KR),JSFIN(KR),JSST(K1),JSFIN(K1)
829   CONTINUE
830   CONTINUE
C     CORRECT SIGNS IN IPRG1
      DO 831 K=1,NSRUN
      IF(IPRG1(K) .EQ. -2) IPRG1(K) = 1
831   CONTINUE
C
C     CHECK FOR OVERLAPS IN THE EXTENDED RUN LIMITS
C     --PREPARE TO MODIFY EXTENDED RUN LIMITS TO EMULATE HOMOMGR: 
C     SAVE EXTENDED RUN LIMITS IN SIMPLE RUN LIMITS (WHICH ARE NO
C     LONGER NEEDED), SINCE SPLICE WILL MODIFY THE EXTENDED RUN
C     LIMITS.
      DO 869 K=1,NSRUN
      ISST(K) = ISTRT(K)
      JSST(K) = JSTRT(K)
      ISFIN(K) = IFIN(K)
      JSFIN(K) = JFIN(K)
869   SSPIJ(K) = SPIJ(K)
C     COUNT CROSS-OVERS
      NCO = 0
C
      OPEN (UNIT=9,NAME='pcl.tmp',status='OLD',FORM='UNFORMATTED')
      IS = 2
C
      DO 876 K=1,NS1
      IF(IPRG1(K) .LE. 0) GO TO 876
      KK = K + 1
      DO 877 M=KK,NSRUN
      IF(IPRG1(M) .LE. 0) GO TO 877
C     TENTATIVELY ASSIGN K1 AND K2 (FIRST AND SECOND RUNS IN
C     LINEAR SEQUENCE) WITH EXTENDED RUN LIMITS
      K1 = K
      K2 = M
      IF(ISTRT(K) .LE. ISTRT(M)) GO TO 870
      K1 = M
      K2 = K
C     CHECK FOR OVERLAP
870   IF(IFIN(K1).LT.ISTRT(K2) .AND. JFIN(K1).LT.JSTRT(K2)) GO TO 877
C
C     ALAS.  NOW WE MUST SPLIT THE CONFLICTING RUNS AT THE RESIDUE
C     WHERE THE DOWNWARD EXTENSION OF THE SECOND SIMPLE RUN (K2)
C     IS LESS PROBABLE THAN THE UPWARD EXTENSION OF THE FIRST
C     SIMPLE RUN (K1).
      LIM = 2
      CALL SPLICE(K1,K2,IS,LIM,KBAD,NCO)
C
C     ELIM ANY RUN THAT IS TRIMMED AWAY BY CROSS-OVER
      IF(KBAD .EQ. 0) GO TO 965
      IPRG1(KBAD) = 0
      GO TO 877
C
C     REPORT ACTION TAKEN (RESULTS STILL IN EXT RUN LIMS)
965   I = IFIN(K1)
      J = ISTRT(K2)
      WRITE(8,1033) K1,K2,IRES(I),IRES(J),K1,SPIJ(K1),K2,SPIJ(K2)
1033  FORMAT(5x,'Runs ',I2,' and ',I2,' cross over between trial',
     1' residues ',A4,'and ',A4,' PSUM(',I2,') =',F6.3,
     2'   PSUM(',I2,') =',F6.3/)
C
877   CONTINUE
876   CONTINUE
      CLOSE (UNIT=9,status='DELETE')
C
C     IF ANY CHANGES, SWITCH BACK SIMPLE AND EXTENDED LIMS
      IF(NCO .EQ. 0) GO TO 885
      DO 880 K=1,NSRUN
      ITMP = ISTRT(K)
      ISTRT(K) = ISST(K)
      ISST(K) = ITMP
      ITMP = JSTRT(K)
      JSTRT(K) = JSST(K)
      JSST(K) = ITMP
      ITMP = IFIN(K)
      IFIN(K) = ISFIN(K)
      ISFIN(K) = ITMP
      ITMP = JFIN(K)
      JFIN(K) = JSFIN(K)
      JSFIN(K) = ITMP
      TEMP = SPIJ(K)
      SPIJ(K) = SSPIJ(K)
880   SSPIJ(K) = TEMP
C
C     COLLECT EQUIV STATS
885   I = 0
      PS1 = 0.0
      D21 = 0.0
      S21 = 0.0
      NUM1 = 0
      DO 851 K=1,NSRUN
      IF(IPRG1(K) .LE. 0) GO TO 851
      I = I + 1
      ISPR(I) = K
      PS1 = PS1 + SSPIJ(K)
      D21 = D21 + SD2(K)
      S21 = S21 + SS2(K)
      NUM1 = NUM1 + ISFIN(K) - ISST(K) + 1
851   CONTINUE
      IF(NUM1 .EQ. 0) GO TO 891
      D = SQRT(D21/NUM1)
      S = SQRT(S21/NUM1)
      WRITE(8,1029) (ISPR(L),L=1,I)
1029  FORMAT(5X,10I5)
891   WRITE(8,1030) PS1,D,S,NUM1
1030  FORMAT(/5x,'Total P(I,J) =',F10.5,5X,'rms D =',F6.2,5X,
     1'rms S =',F6.2,5X,'for',I5,' equivalent residues')
C
C     WITH PROGRESSION RULE APPLIED TO EXTENDED RUNS ("ALTERNATE" WAY)
C
      WRITE(8,1031)
1031  FORMAT(/5x,'Best equivalence using progression rule',
     1' on extended runs'/)
D     write(8,2002)
      DO 840 K=1,NR1
      IF(IPRG2(K) .LT. 0) GO TO 840
      L = K + 1
      DO 839 K1=L,NRUN
      IF(IPRG2(K1) .LE. 0) GO TO 839
      IF(IFIN(K) .GE. ISTRT(K1)) GO TO 833
      IF(JFIN(K) .LT. JSTRT(K1)) GO TO 839
      GO TO 834
833   IF(IFIN(K1) .GE. ISTRT(K)) GO TO 834
      IF(JFIN(K1) .LT. JSTRT(K)) GO TO 839
834   IF(SPIJ(K) - SPIJ(K1)) 836,837,838
836   IPRG2(K) = 0
D     write(8,2003) K,K1,K,ISTRT(K1),IFIN(K1),ISTRT(K),IFIN(K),
D    1 JSTRT(K1),JFIN(K1),JSTRT(K),JFIN(K)
      GO TO 839
837   WRITE(8,1025) K,K1
1025  FORMAT(//' *** AMBIGUOUS CONFLICT BETWEEN RUNS ',I2,' AND '
     1,I2,' ***'//)
838   IPRG2(K1) = 0
D     write(8,2003) K,K,K1,ISTRT(K),IFIN(K),ISTRT(K1),IFIN(K1),
D    1 JSTRT(K),JFIN(K),JSTRT(K1),JFIN(K1)
839   CONTINUE
840   CONTINUE
      I1 = 0
      PS2 = 0.0
      D22 = 0.0
      S22 = 0.0
      NUM2 = 0
      DO 852 K=1,NRUN
      IF(IPRG2(K) .LE. 0) GO TO 852
      I1 = I1 + 1
      IEPR(I1) = K
      PS2 = PS2 + SPIJ(K)
      D22 = D22 + SD2(K)
      S22 = S22 + SS2(K)
      NUM2 = NUM2 + IFIN(K) - ISTRT(K) + 1
852   CONTINUE
      IF(NUM2 .EQ. 0) GO TO 893
      D = SQRT(D22/NUM2)
      S = SQRT(S22/NUM2)
      WRITE(8,1029) (IEPR(L),L=1,I1)
893   WRITE(8,1030) PS2,D,S,NUM2
C
C     WITHOUT PROGRESSION RULE
C
      WRITE(8,1032)
1032  FORMAT(/5x,'Best equivalence without using progression rule'/)
D     write(8,2002)
      DO 850 K=1,NR1
      IF(NPROG(K) .LT. 0) GO TO 850
      L = K + 1
      DO 849 K1 = L,NRUN
      IF(NPROG(K1) .LE. 0) GO TO 849
      JF = JFIN(K)
      JFK = JFIN(K1)
      JS = JSTRT(K)
      JSK = JSTRT(K1)
      JOFF = JOF(K)
      IF(JOFF .NE. JOF(K1)) GO TO 843
      IF(JOFF .EQ. -1) GO TO 842
      IF(JF .LT. JSK .OR. JFK .LT. JS) GO TO 848
      GO TO 844
842   IF(JF .GT. JSK .OR. JFK .GT. JS) GO TO 848
      GO TO 844
843   MOD = 1
      IF(JOFF .EQ. -1) MOD = MOD + 1
      IF(JFK .LT. JS) MOD = MOD + 2
      GO TO (931,932,933,934),MOD
931   IF(JF .LT. JFK) GO TO 848
      GO TO 844
932   IF(JS .LT. JSK) GO TO 848
      GO TO 844
933   IF(JSK .LT. JS) GO TO 848
      GO TO 844
934   IF(JFK .LT. JF) GO TO 848
844   IF(SPIJ(K) - SPIJ(K1)) 845,846,847
845   NPROG(K) = 0
D     write(8,2003) K,K1,K,ISTRT(K1),IFIN(K1),ISTRT(K),IFIN(K),
D    1 JSTRT(K1),JFIN(K1),JSTRT(K),JFIN(K)
      GO TO 849
848   IF(IFIN(K) .LT. ISTRT(K1) .OR. IFIN(K1) .LT. ISTRT(K)) GO TO 849
      GO TO 844
846   WRITE(8,1025) K,K1
847   NPROG(K1) = 0
D     write(8,2003) K,K,K1,ISTRT(K),IFIN(K),ISTRT(K1),IFIN(K1),
D    1 JSTRT(K),JFIN(K),JSTRT(K1),JFIN(K1)
849   CONTINUE
850   CONTINUE
      J = 0
      PS3 = 0.0
      D23 = 0.0
      S23 = 0.0
      NUM3 = 0
      DO 855 K=1,NRUN
      IF(NPROG(K) .LE. 0) GO TO 855
      J = J + 1
      INPR(J) = K
      PS3 = PS3 + SPIJ(K)
      D23 = D23 + SD2(K)
      S23 = S23 + SS2(K)
      NUM3 = NUM3 + IFIN(K) - ISTRT(K) + 1
855   CONTINUE
      IF(NUM3 .EQ. 0) GO TO 895
      D = SQRT(D23/NUM3)
      S = SQRT(S23/NUM3)
      WRITE(8,1029) (INPR(L),L=1,J)
895   WRITE(8,1030) PS3,D,S,NUM3
C
C     STORE EQUIV FOR NEXT ITERATION
C
      IF(IPROG .EQ. 0) return
      IF(ITER .NE. ISTAR .OR. NEQ .NE. 0) GO TO 857
      DO 858 I=2,NI
      IOLD(I) = 0
858   IEQ(I) = 0
      GO TO 859
857   DO 856 I=2,NI
      IOLD(I) = IEQ(I)
856   IEQ(I) = 0
859   DO 910 K=1,NRUN
      GO TO (901,902,903),IPROG
901   IF(IPRG1(K) .LE. 0) GO TO 910
      JCNT(K) = JSST(K) - JOF(K)
      I1 = ISST(K)
      I2 = ISFIN(K)
904   DO 905 I=I1,I2
      J = JCNT(K) + JOF(K)
      JCNT(K) = J
      IEQ(I) = J
905   CONTINUE
      GO TO 910
902   IF(IPRG2(K) .LE. 0) GO TO 910
906   JCNT(K) = JSTRT(K) - JOF(K)
      I1 = ISTRT(K)
      I2 = IFIN(K)
      GO TO 904
903   IF(NPROG(K) .LE. 0) GO TO 910
      GO TO 906
910   CONTINUE
C
C     ANY CHANGES FROM LAST CYCLE?
C
      N = 0
      NA = 0
      ND = 0
      NC = 0
      DO 915 I=2,NI
      I1 = 1
      IF(IOLD(I) .NE. 0) I1 = I1 + 1
      IF(IEQ(I) .NE. 0) I1 = I1 + 2
      GO TO (915,911,912,913),I1
911   ND = ND + 1
      GO TO 915
912   NA = NA + 1
      N = N + 1
      GO TO 915
913   N = N + 1
      IF(IOLD(I) .NE. IEQ(I)) NC = NC + 1
915   CONTINUE
      NEQ = N
      WRITE(8,1049) N,NA,ND,NC
1049  FORMAT(//' ***',I4,' EQUIVALENT RESIDUES for next cycle: ',I3,
     1' added, ',I3,' deleted and ',I3,' changed from last cycle')
      IF(NA+ND+NC .EQ. 0 .OR. ITER .EQ. ISTAR+ITMAX-1) return
      ITER = 0
      RETURN
      END
      SUBROUTINE SPLICE(K1,K2,IS,LIM,KBAD,NUM)
C     SPLICE IS CALLED BY RUNS OR SELEQ TO TRIM OVERLAPPING
C     ENDS OF RUNS K1 AND K2 ... THE CALLING PROGRAM DECIDES WHEN
C     SUCH MIXING & MATCHING IS APPRORIATE.
C
C     NOTE:
C       1) SCRATCH FILE 9 SHOULD BE OPENED BY THE CALLING PROGRAM;
C     IS SHOULD POINT TO THE NEXT RECORD TO BE READ.
C       2) IF EITHER K1 OR K2 IS SHORTER THAN LIM RESIDUES AFTER
C     TRIMMING, THE SPLICE IS RETRACTED: THE INITIAL RUN LIMITS ARE
C     RESTORED AND THE OVER-TRIMMED RUN NUMBER IS RETURNED IN KBAD.
C     THE CALLING PROGRAM MUST DECIDE IF THE RUN IS TO BE ELIMINATED.
C       3) OTHERWISE, THE *EXTENDED* RUN LIMITS AND SPIJ ARE SET
C     TO THE NEW VALUES, KBAD=0, AND NUM IS INCREMENTED.
C       4) SPLICING ONLY WORKS FOR RUN THAT GO IN THE SAME DIRECTION.
C
      COMMON /SEQ/ IRES(702),JRES(702),AIRES(702),AJRES(702)
	COMMON /SAVE/ NSV(702)
      COMMON /SCRAT/ JSV(10),PSV(10),D2SV(10),S2SV(10),DXSV(3,10)
      COMMON /DIREC/ NRUN,NSRUN,ISFIN(60),JSFIN(60),ISST(60),JSST(60),
     1JOF(60),ISTRT(60),JSTRT(60),IFIN(60),JFIN(60),SSPIJ(60),SPIJ(60),
     2SD2(60),SS2(60),SDX(60),SDY(60),SDZ(60),JCNT(60),IPRG1(60),
     4IPRG2(60),NPROG(60),ISPR(60),IEPR(60),INPR(60),NEXT(60),IGAP(60)
      DIMENSION IO(70),PT1(10),PT2(10)
      EQUIVALENCE (IO(1),JSV(1))
C
D     write(8,2010) K1,K2
D2010 FORMAT(' *FIND CROSS-OVER BETWEEN RUNS ',I2,' AND ',I2/
D    1 '  I1 IRES  JRES1   PIJ1  JRES2   PIJ2  KSIMPLE',5X,
D    2 'IST-FI R1',3X,'JST-FI R1',3X,'IST-FI R2',3X,
D    3 'JST-FI R2',6X,'J1   J2')
C
C     SAVE INITIAL LIMITS & SET UP NORMAL RETURN
C
      IS1 = ISTRT(K1)
      JS1 = JSTRT(K1)
      IF1 = IFIN(K1)
      JF1 = JFIN(K1)
      SP1 = SPIJ(K1)
      IS2 = ISTRT(K2)
      JS2 = JSTRT(K2)
      IF2 = IFIN(K2)
      JF2 = JFIN(K2)
      SP2 = SPIJ(K2)
      JO1 = JOF(K1)
      JO2 = JOF(K2)
      NUM = NUM + 1
      KBAD = 0
      LM1 = LIM - 1
C
C     POSITION UNIT 9 AT START OF K1; NOTE THAT "IS" ALWAYS
C     POINTS TO THE NEXT RECORD TO BE READ.
      I = ISTRT(K1) - 1
      IF(I .EQ. IS) GO TO 883
      IF = I - 1
      IF(IF .GE. IS) GO TO 871
      IS = 2
      REWIND 9
871   DO 875 I1=IS,IF
      IF(NSV(I1) .NE. 0) READ(9)
875   CONTINUE
883   IF(NSV(I) .NE. 0) READ(9) IO
C     SET SCAN LIMITS TO INCLUDE ALL OF BOTH RUNS
      IS = I + 1
      IF = IFIN(K1)
      IF(IFIN(K2) .GT. IF) IF = IFIN(K2)
C     SET J1 AND J2 TO CORRESPOND TO THE START OF RUN K1 (-1)
      J1 = JSTRT(K1) - JO1
      J2 = JSTRT(K2) + (IS-ISTRT(K2)) - JO2
C
C     SCAN THRU PIJ'S TO FIND CROSS-OVER POINT;
C     KEEP TRACK OF THE TOTAL PROBABILITY BEFORE THE C-O POINT
C     IN PS1 AND PS2 FOR RUNS K1 & K2, RESP.
      PS1 = 0.0
      PS2 = 0.0
      DO 940 L=1,10
      PT1(L) = 0.0
940   PT2(L) = 0.0
      N = 0
      KS = 0
      DO 879 I1=IS,IF
      READ(9) IO
C     COUNT RECORDS READ FROM UNIT 9 IN N
      N = N + 1
      J1 = J1 + JO1
      J2 = J2 + JO2
      NSAVE = NSV(I1)
      P1 = 0.0
      P2 = 0.0
      DO 878 L=1,NSAVE
      IF(JSV(L) .EQ. J1) P1 = PSV(L)
      IF(JSV(L) .EQ. J2) P2 = PSV(L)
878   CONTINUE
C     INSURE THAT CORRESPONDENCE IS NOT FORTUITOUS
      IF(I1 .LT. ISTRT(K1) .OR. I1 .GT. IFIN(K1)) P1 = 0.0
      IF(I1 .LT. ISTRT(K2) .OR. I1 .GT. IFIN(K2)) P2 = 0.0
      IF(KS .NE. 0) GO TO 882
C     NOW, IF KS=0, THE START OF A SIMPLE RUN HASN'T BEEN FOUND;
C     THIS MUST BE REQUIRED TO INSURE THAT THE CROSS-OVER IS
C     BETWEEN THE SIMPLE RUNS.
      IF(KS .NE. 0) GO TO 882
      IF(JSV(1) .NE. J1 .OR. P1 .EQ. 0.0) GO TO 874
C     OK, HAVE FOUND SIMPLE RUN, AND ASSIGNMENT OF K1 & K2 IS OK
      KS = 1
      GO TO 882
874   IF(JSV(1) .NE. J2 .OR. P2 .EQ. 0.0) GO TO 873
C     HAVE FOUND SIMPLE RUN, BUT K1 & K2 ARE REVERSED ... SHOULD
C     BE FIRST & SECOND SIMPLE RUNS IN SEQUENCE; APPARENTLY K2
C     CAN BE EXTENDED DOWN FARTHER THAN K1.  SWITCH K1 & K2.
D     write(8,2005)
D2005 FORMAT(' SWAP')
      ITMP = K1
      K1   = K2
      K2   = ITMP
      ITMP = J1
      J1   = J2
      J2   = ITMP
      ITMP = JO1
      JO1  = JO2
      JO2  = ITMP
      P  = P1
      P1 = P2
      P2 = P
      P   = PS1
      PS1 = PS2
      PS2 = P
      DO 955 L=1,10
      P = PT1(L)
      PT1(L) = PT2(L)
955   PT2(L) = P
      KS = 1
C
C     CHECK FOR CROSSOVER
882   IF(P1 .GT. P2) GO TO 873
C     ADJUST RUN LIMITS ON BASIS OF C-O ON I
      IF(I1-1 .LT. IFIN(K1)) IFIN(K1) = I1 - 1
      IF(J1-JO1 .LT. JFIN(K1)) JFIN(K1) = J1 - JO1
      IF(I1 .GT. ISTRT(K2)) ISTRT(K2) = I1
      IF(JO2 .EQ.  1 .AND. J2 .GT. JSTRT(K2)) JSTRT(K2) = J2
      IF(JO2 .EQ. -1 .AND. J2 .LT. JSTRT(K2)) JSTRT(K2) = J2
      SPIJ(K1) = PS1
      SPIJ(K2) = SPIJ(K2) - PS2
D     write(8,2008)
D2008 FORMAT(' *CROSS-OVER POINT IN TRIAL SEQ ("ON I"):')
D     write(8,2006) I1,IRES(I1),JRES(J1),P1,JRES(J2),P2,KS
D     write(8,2007) IRES(ISTRT(K1)),IRES(IFIN(K1)),
D    1 JRES(JSTRT(K1)),JRES(JFIN(K1)),IRES(ISTRT(K2)),
D    2 IRES(IFIN(K2)),JRES(JSTRT(K2)),JRES(JFIN(K2)),J1,J2
D2007 FORMAT(52X,4(A4,'-',A4,3X),2I5)
      GO TO 880
C
873   PS1 = PS1 + P1
      PS2 = PS2 + P2
D     write(8,2006) I1,IRES(I1),JRES(J1),P1,JRES(J2),P2,KS
D2006 FORMAT(1X,I3,1X,A4,3X,A4,F7.3,3X,A4,F7.3,3X,I3)
C     SAVE THE LAST 10 PIJ'S OF EACH RUN IN PT1 OR PT2
C     SHIFT EXCESS OF END
      IF(P1 .EQ. 0.0) GO TO 945
      DO 941 L=9,1,-1
941   PT1(L+1) = PT1(L)
      PT1(1) = P1
C     NO NEED TO KEEP TRACK OF PIJ FOR K2, ONCE WE ARE SURE
C     IT IS SECOND SIMPLE RUN.
945   IF(KS .NE. 0) GO TO 879
      IF(P2 .EQ. 0.0) GO TO 879
      DO 946 L=9,1,-1
946   PT2(L+1) = PT2(L)
      PT2(1) = P2
879   CONTINUE
C
C     NOW CHECK FOR OVERLAP ON J; NOTE THAT EVERYBODY RETAINS
C     VALUES SET AT C-O POINT ON I.
880   IF(JO1 .EQ.  1 .AND. J1 .LE. J2) GO TO 890
      IF(JO1 .EQ. -1 .AND. J2 .LE. J1) GO TO 890
D     write(8,2012)
D2012 FORMAT(' *FIND CROSS-OVER IN REF SEQUENCE ("ON J")'/
D    1 '  I1 IRES J1-JT  PT1(JT) JRES2   PIJ2    JT')
C     FIND INDEX OF PT1 CORRESPONDING TO J2
c
c problem (sjr): jt can be > 10. hack: where it violates array bounds, return 0
c
      JT = J1 - J2
      IF(JO1 .EQ. -1) JT = J2 - J1
C     CAN WE JUST TRUNCATE K1?
	xpt=0.
	if(jt.le.10) xpt=pt1(jt)
      IF(xpt .LE. P2) GO TO 886
      PS2 = P2
      IF(JT .EQ. 1) GO TO 888
C     NO; UPDATE IS AND SCAN FOR C-O ON J
      IS = IS + N
      N = 0
      IF = IS + JT - 1
      DO 885 I1=IS,IF
      JT = JT - 1
      READ(9) IO
      N = N + 1
      J2 = J2 + JO2
      NSAVE = NSV(I1)
      P2 = 0.0
      DO 881 L=1,NSAVE
      IF(JSV(L) .NE. J2) GO TO 881
      P2 = PSV(L)
      GO TO 884
881   CONTINUE
884   IF(I1 .LT. ISTRT(K2) .OR. I1 .GT. IFIN(K2)) P2 = 0.0
c
c hack
c
      xpt=0.
	if(jt.le.10) xpt=pt1(jt)
	IF(xpt .LE. P2) GO TO 886
      PS2 = PS2 + P2
D     write(8,2006) I1,IRES(I1),JRES(J1-JT),PT1(JT),JRES(J2),P2,JT
885   CONTINUE
C     FINAL LIMIT CORRECTIONS WHEN CROSS-OVER ON J IS AFTER
C     THE END OF K1, OR WHEN (JT.EQ.1 .AND. P2.LT.PT1(JT)).
888   ISTRT(K2) = I1 + 1
      JSTRT(K2) = J2 + JO2
      GO TO 889
C
C     APPLY FINAL CORRECTIONS TO RUN LIMITS
C     NOTE THAT YOUR FORTRAN MUST PRESERVE I1 OUT OF THE LOOP
886   IFIN(K1) = IFIN(K1) - JT
      JFIN(K1) = JFIN(K1) - JO1*JT
      IF(I1 .GT. ISTRT(K2)) ISTRT(K2) = I1
      IF(JO2 .EQ.  1 .AND. J2 .GT. JSTRT(K2)) JSTRT(K2) = J2
      IF(JO2 .EQ. -1 .AND. J2 .LT. JSTRT(K2)) JSTRT(K2) = J2
D     write(8,2013)
D2013 FORMAT(' *CROSS-0VER POINT ON J:')
D     write(8,2006) I1,IRES(I1),JRES(J1-JT),PT1(JT),JRES(J2),P2,JT
      PS1 = 0.0
c
c hack
c
	xpt=0.
	if(jt.le.10) xpt=pt1(jt)
	DO 887 I1=1,JT
887   PS1 = PS1 + xpt
      SPIJ(K1) = SPIJ(K1) - PS1
889   SPIJ(K2) = SPIJ(K2) - PS2
D     write(8,2007) IRES(ISTRT(K1)),IRES(IFIN(K1)),
D    1 JRES(JSTRT(K1)),JRES(JFIN(K1)),IRES(ISTRT(K2)),
D    2 IRES(IFIN(K2)),JRES(JSTRT(K2)),JRES(JFIN(K2)),J1,J2
C
C     CHECK THAT NEITHER RUN WAS TRUNCATED AWAY
890   IF(IFIN(K1)-ISTRT(K1) .GT. LM1) GO TO 960
      KBAD = K1
D     write(8,2014) K1
D2014 FORMAT(' *RUN ',I2,' TOO SHORT AFTER TRUNCATION')
      GO TO 962
960   IF(IFIN(K2)-ISTRT(K2) .GT. LM1) GO TO 965
      KBAD = K2
D     write(8,2014) K2
C     RESTORE ORIGINAL RUN LIMITS
962   ISTRT(K1) = IS1
      JSTRT(K1) = JS1
      IFIN(K1) = IF1
      JFIN(K1) = JF1
      SPIJ(K1) = SP1
      ISTRT(K2) = IS2
      JSTRT(K2) = JS2
      IFIN(K2) = IF2
      JFIN(K2) = JF2
      SPIJ(K2) = SP2
      NUM = NUM - 1
C
C     UPDATE IS
C
965   IS = IS + N
      RETURN
      END
