PROGRAM ligands

c	PROGRAM looks for S,N,O atoms near 6 Fe's and calculates stats

	CHARACTER*76	line, filename
	CHARACTER*13	lab(6,99)/594*'             '/
	REAL		fe(3,6), x(3), lig(6,99)/594*0./, d(6)
	INTEGER		num(6)/6*0/

1	WRITE (6, 1000) 'Enter input PDB filename:'
1000	FORMAT (a)
	READ (5, 1000) filename
	OPEN (unit=1, form='FORMATTED', status='OLD', file=filename,
	1	err=1)

c	First get Fe coordinates

	nA = ICHAR('A')-1
	DO WHILE (.TRUE.)
5	  READ (1, 1000, end=9) line
	  IF (line(1:4) .NE. 'ATOM' .AND. line(1:6) .NE. 'HETATM') GOTO 5
	  IF (line(14:15) .EQ. 'FE') THEN
	    n = MOD(ICHAR(line(73:73))-nA-1,6)+1
	    READ (line(31:54), 1001) (fe(i,n), i=1,3)
1001	    FORMAT (3f8.3)
	  endIF
	endDO

c	Close and reopen PDB file

9	CLOSE (unit=1)
	OPEN (unit=1, form='FORMATTED', status='OLD', file=filename)

c	Get radius for ligation sphere

	WRITE (6, 1000) 'Enter radius of ligation sphere:'
	READ (5, 1001) rad

c	Compile ligand atoms

	rad = rad*rad
	DO WHILE (.TRUE.)
15	  READ (1, 1000, end=19) line
	  IF (line(1:4) .NE. 'ATOM' .AND. line(1:6) .NE. 'HETATM') GOTO 15
	  IF (line(14:14) .EQ. 'O' .OR. line(14:14) .EQ. 'N' .OR.
	1	line(14:14) .EQ. 'S') THEN
	    n = MOD(ICHAR(line(73:73))-nA-1,6)+1
	    READ (line(31:54), 1001) x
	    dd = (x(1)-fe(1,n))**2+(x(2)-fe(2,n))**2+(x(3)-fe(3,n))**2
	    IF (dd .LE. rad) THEN
	      num(n) = num(n)+1
	      lig(n,num(n)) = SQRT(dd)
	      lab(n,num(n)) = line(14:26)
	    endIF
	  endIF
	endDO

19	CLOSE (unit=1)

	DO i = 1, num(1)
	  d(1) = lig(1,i)
	  DO j = 2, 6
	    DO k = 1, num(j)
	      IF (lab(1,i) .EQ. lab(j,k)) THEN
	        d(j) = lig(j,k)
	        GOTO 25
	      endIF
	    endDO
	  GOTO 30
25	  continue
	  endDO
	  dd = (d(1) + d(2) + d(3) + d(4) + d(5) + d(6))/6	
	  dd2 = (d(1)**2 + d(2)**2 + d(3)**2 + d(4)**2 + d(5)**2 + d(6)**2)/6	
	  sig = SQRT(dd2-dd**2)
	  WRITE (6, 1005) lab(1,i), d, dd, sig
1005	  FORMAT (a, 3x, 6f6.3, 2x, 2f6.3)
30	  continue
	endDO

	end