PROGRAM	limits

c	Calculate CM, radii, distances from local symmetry axes for
c	a PDB file

	CHARACTER*66	line
	REAL*8		s1, s2, s31, s32

c	First open file

	WRITE (6, 1000) 'Enter name of PDB file:'
1000	FORMAT (a)
	READ (5, 1000) line
	OPEN (unit=1, form='FORMATTED', status='OLD', readonly,
	1	file=line)

c	Get position to use

	WRITE (6, 1000) 'Enter center position to use for distances:'
	READ (5, 1001, end=1) xc, yc, zc
1001	FORMAT (3f10.4)

c	Scan through file

1	s1 = 0.0d0
	s2 = 0.0d0
	s31 = 0.0d0
	s32 = 0.0d0
	r1min = 999.
	r2min = 999.
	r31min = 999.
	r32min = 999.
	r1max = 0.
	r2max = 0.
	r31max = 0.
	r32max = 0.
	num = 0

	DO WHILE (.TRUE.)

5	  READ (1, 1000, end=9) line

	  IF (line(18:20) .EQ. 'WAT') GOTO 5
	  IF (line(14:17) .NE. 'CA  ') GOTO 5
	  num = num + 1
	  READ (line(31:54), 1002) x, y, z
1002	  FORMAT (3f8.4)

	  r1 = (x-xc)**2 + (y-yc)**2 + (z-zc)**2
	  r = x*x + y*y + z*z
	  r2 = x*x + y*y
	  r31 = r-(x+y+z)**2/3.
	  r32 = r-(x-y+z)**2/3.

	  s1 = s1 + r1
	  s2 = s2 + r2
	  s31 = s31 + r31
	  s32 = s32 + r32

	  r1min = MIN(r1,r1min)
	  r1max = MAX(r1,r1max)
	  r2min = MIN(r2,r2min)
	  r2max = MAX(r2,r2max)
	  r31min = MIN(r31,r31min)
	  r31max = MAX(r31,r31max)
	  r32min = MIN(r32,r32min)
	  r32max = MAX(r32,r32max)

	endDO

9	CLOSE (unit=1)

	WRITE (6, 1003) 'From  point:', SQRT(r1min), SQRT(r1max),
	1		SQRT(s1/num)
	WRITE (6, 1003) 'From  (010):', SQRT(r2min), SQRT(r2max),
	1		SQRT(s2/num)
	WRITE (6, 1003) 'From  (111):', SQRT(r31min), SQRT(r31max),
	1		SQRT(s31/num)
	WRITE (6, 1003) 'From (1-11):', SQRT(r32min), SQRT(r32max),
	1		SQRT(s32/num)
1003	FORMAT (a, 3f10.4)

	end