PROGRAM chksolv

c	Applies symmetry to potential solvent list and looks for redundancy

	PARAMETER	NDIM = 720

	CHARACTER*80	title
	CHARACTER*15	label(NDIM/24)
	REAL		mat(3,3,5), dummy(3,3), x(NDIM), y(NDIM), z(NDIM),
	1		pk(NDIM/24), dh(NDIM/24)

	aa = 197.17
	bb = 127.03
	cc = 134.18
	beta = 97.64
	conv = 45./atan(1.)
	sinb = sin(beta/conv)
	cosb = cos(beta/conv)

	WRITE (6, 1000) 'Enter starting solvent number:  '
1000	FORMAT (a)

	READ (5, 1001) nstart
1001	FORMAT (i5)
	WRITE (6, 1001) nstart

c	Read local symmetry operations

	OPEN (unit=10, status='OLD', form='FORMATTED', readonly,
	1	file='AVER.MEP')

	READ (10, 1002)
	DO n = 1,5
	  READ (10, 1002) ((mat(i,j,n), i=1,3), j=1,3)
1002	  FORMAT (9f10.5)
	  mat(1,2,n) = mat(1,2,n)*aa/bb
	  mat(1,3,n) = mat(1,3,n)*aa/cc
	  mat(2,1,n) = mat(2,1,n)*bb/aa
	  mat(2,3,n) = mat(2,3,n)*bb/cc
	  mat(3,1,n) = mat(3,1,n)*cc/aa
	  mat(3,2,n) = mat(3,2,n)*cc/bb
	endDO

	CLOSE unit=10

c	Get file names

	WRITE (6, 1000) 'Enter filename for input solvents:'
	READ (5, 1000) title
	WRITE (6, 1000) title
	OPEN (unit=1, status='OLD', file=title, readonly, form='FORMATTED')

	WRITE (6, 1000) 'Enter filename for output solvents:'
	READ (5, 1000) title
	WRITE (6, 1000) title
	OPEN (unit=2, status='UNKNOWN', form='FORMATTED', file=title,
	1	carriagecontrol='LIST')

c	Store solvent list

	DO n = 1, 100000
	  READ (1, 1003, end=5) pk(n), label(n), x(n), y(n), z(n), dh(n)
1003	  FORMAT (f15.4, a15, 4f10.3)
	endDO

5	num = n - 1
	CLOSE (unit=1)

c	Apply local symmetry to atoms

	natom = num
	DO m = 1,5
	  DO n = 1,num
	    natom = natom + 1
	    x(natom) = x(n)*mat(1,1,m) + y(n)*mat(1,2,m) + z(n)*mat(1,3,m)
	    y(natom) = x(n)*mat(2,1,m) + y(n)*mat(2,2,m) + z(n)*mat(2,3,m)
	    z(natom) = x(n)*mat(3,1,m) + y(n)*mat(3,2,m) + z(n)*mat(3,3,m)
	   endDO
	endDO

c	Apply I2 symmetry to atoms

	DO n = 1,num*6
	  natom = natom + 1
	  x(natom) = -x(n)
	  y(natom) = y(n)
	  z(natom) = -z(n)
	endDO
	DO n = 1,num*12
	  natom = natom + 1
	  x(natom) = x(n)+aa/2.
	  y(natom) = y(n)+bb/2.
	  z(natom) = z(n)+cc/2.
	endDO

c	Check atoms for redundancy

	DO n = 1,num
	  DO m = n+num,num*24,num
	    xx = MOD(x(n)-x(m)+5*aa,aa)
	    IF (xx .GT. aa/2.) xx = xx - aa
	    yy = MOD(y(n)-y(m)+5*bb,bb)
	    IF (yy .GT. bb/2.) yy = yy - bb
	    zz = MOD(z(n)-z(m)+5*cc,cc)
	    IF (zz .GT. cc/2.) zz = zz - cc
	    zz = zz + xx*cosb
	    xx = xx*sinb
	    d = xx*xx + yy*yy + zz*zz
	    IF (d .LE. 4.) pk(n) = 0.0
	  endDO
100	  continue
	endDO

	DO n = 1,num
	  IF (pk(n) .EQ. 0) GOTO 200
	  DO m = n+1,num*24
	    xx = MOD(x(n)-x(m)+5*aa,aa)
	    IF (xx .GT. aa/2.) xx = xx - aa
	    yy = MOD(y(n)-y(m)+5*bb,bb)
	    IF (yy .GT. bb/2.) yy = yy - bb
	    zz = MOD(z(n)-z(m)+5*cc,cc)
	    IF (zz .GT. cc/2.) zz = zz - cc
	    zz = zz + xx*cosb
	    xx = xx*sinb
	    d = xx*xx + yy*yy + zz*zz
	    IF (d .LE. 4.) pk(MOD(m,num)) = 0.0
	  endDO
200	  continue
	endDO

	DO n = 1, num
	  IF (pk(n) .NE. 0.) THEN
	    WRITE (2, 1004) 'WAT ', nstart, 'O   ', x(n), y(n), z(n), 
	1	10., 1., pk(n), label(n), dh(n)
1004	    FORMAT (3x, a4, 1x, i3, a4, 6f10.5, a15, f10.5)
	    nstart = nstart + 1
	  endIF
	endDO

c	Always add a some solvents around origin in case I want to
c	move them into density

	DO n = 1,10
	  WRITE (2, 1004) 'WAT ', nstart, 'O   ', FLOAT(n), FLOAT(n),
	1	FLOAT(n), 10., 1.
	  nstart = nstart + 1
	endDO

	CLOSE (unit=2)

	end