C This is for self distances-  neglect distances within residue
C Don't need to read the file twice!
C If the purpose is for studying sec structure, better to make
C 2D array- residue number vs Ca-C-O-N
C for each N measure distance to O in 4th prvious residue
C or  distance O to all N's to look for beta sheet.
C note N is closer to O in prvious residue (same peptide bond, 2.2A)
C than to any H-bonded O.

	CHARACTER*40 INFILE1,INFILE2,FILENAME
	LOGICAL EFLG
	real*8 M1(3,3000),M2(3),d
	INTEGER*4 I,J,k
	character*14 nres(3000),nres2

1200	FORMAT(A)
	write(6,*)'Read 2 .PDB files and list all pairs of atoms '
	write(6,*) 'in one file within threshold distance'
	write(6,*) 'of atom in other file'


c 	read first file into m1
6	WRITE(6,*) 'ENTER NAME OF 1st pdb file (atom recds only):'
	READ(5,1200) INFILE1
	IF (INFILE1(:5).EQ.'     ') STOP
	  INQUIRE(FILE=INFILE1, NAME=FILENAME,EXIST=EFLG)
	  IF (.NOT.EFLG) GOTO 6
c        OPEN (UNIT=2,FILE=INFILE1,STATUS='OLD',READONLY)
        OPEN (UNIT=2,FILE=INFILE1,STATUS='OLD')
8	WRITE(6,*) 'ENTER NAME OF 2nd pdb file (atom recds only):'
	READ(5,1200) INFILE2
	IF (INFILE2(:5).EQ.'     ') STOP
	  INQUIRE(FILE=INFILE2, NAME=FILENAME,EXIST=EFLG)
	  IF (.NOT.EFLG) GOTO 8
c        OPEN (UNIT=3,FILE=INFILE2,STATUS='OLD',READONLY)
        OPEN (UNIT=3,FILE=INFILE2,STATUS='OLD')

	write(6,*)'List atoms closer than R. R=?'
	read (5,*) thresh
	write(6,*)'Listing intermolecular contacts closer than ',thresh
	write(6,*) '   between ',INFILE1(:15),' and ',INFILE2(:15)
	k=1
c 	read coord         
50	read(2,51,end=85) nres(k),(m1(i,k),I=1,3)
c51 	format (30x,3d8.3,2f6.2)
51      format (12x,A14,4x,3d8.3,2f6.2)
c	write(6,*) k,nres(k),(m1(i,k),I=1,3)
	k=k+1
	if (k.gt.3000) write(6,*) 'dimensioned for <3001 atoms in first file!'
	if (k.gt.3000) goto 85
80	goto 50


c123456789012345678901234567890
c            123456789012345678901234567890
cATOM    450  CA  ALA J  61      28.433 161.994  41.553  1.00 99.97          

	
85	read(3,51,end=1160) nres2,(m2(i),I=1,3)
c	write(6,*) nres2,(m2(i),I=1,3)


c	do 110 j=1,2
	do 110 j=1,k-1
	if (nres(j)(10:14).ne.nres2(10:14)) then
88	  x=0.
	  do 90 i=1,3
90	  x=x+(m1(i,j)-m2(i))**2
	  d=sqrt(x)
c	  write(6,*) j,nres(j),(m1(i,j),I=1,3), nres2,(m2(i),I=1,3),d
	  if (d.le.thresh) write (6,111) nres(j),nres2,d
		end if
110	continue
111	format (' ',2A20,16f6.1)

	goto 85

1160	CLOSE(UNIT=2)
	END