c read PDB file in orthogonal system (later read any pdb, convert to orth axes)
c     no, read RDI. convert pdb to rdi w pdb2rdi AFTER REMOVING nonatom records
c operate with specified rotation + translation
c put into arbitrary unit cell (given orth axis relative to cell vectors)
c (write out this file for use by sfclc)
c read MEP file
c check for overlap of protein between asu's:
c	for each atom
c	operate on atom with each MEP operation
c	check atom in new position for overlap with each atom of primary position
c	 (If primary position does not overlap with any other, that ensures
c	  there will be no overlap between an positions?)


	REAL*4 U(3),V(3),W(3),RM(3,3),A,B,ANGLE,COSA,SINA,SYMOP(3,4,12)
	REAL*4 D(3),CELLA,CELLB,CELLC,ALPHA,BETA,GAMMA,CELL(3,3)
	REAL*8 RCELL(4,4)
	INTEGER I,J,K,I1,I2,I3
	CHARACTER*15 RESI
	CHARACTER*40 FNAME
	CHARACTER*1 AXIS
	character*80 astring

	WRITE(6,*) 'PROGRAM TO generate sym equiv positions.        '
	WRITE(6,*) '1st position given in A or fractional units      '
	WRITE(6,*) 'Sym equiv listed in fractional and A units      '

	pi=4*atan(1.0)
	radius=5

	WRITE(6,*)'ENTER A B C ALPHA BETA GAMMA:'
	READ(5,*) CELLA,CELLB,CELLC,ALPHA,BETA,GAMMA
	alpha=alpha*pi/180.
	beta=beta*pi/180
	gamma=gamma*pi/180

	CELL(1,1)=CELLA
	CELL(2,1)=0
	CELL(3,1)=0
	CELL(1,2)=CELLB*COS(GAMMA)
	CELL(2,2)=CELLB*SIN(GAMMA)
	CELL(3,2)=0
	CELL(1,3)=CELLC*COS(BETA)
	CELL(2,3)=CELLC*COS(ALPHA)*SIN(BETA)
	CELL(3,3)=CELLC*SIN(ALPHA)*SIN(BETA)

	DO 405 I=1,3
	DO 405 J=1,3
405	RCELL(I,J)=CELL(I,J)
	N=3
	m=3
	call matprn(rcell,n,m)
	CALL MATINV(RCELL,N,DET)
	write(6,*)
	call matprn(rcell,n,m)

c	TEST MATRIX INVERSION
	do 407 i=1,3
	do 406 j=1,3
	v(j)=0.
	do 406 k=1,3
406	v(j)=v(j)+rcell(i,k)*cell(k,j)
407	write(6,*)(v(j),j=1,3)

C IF V IS COOORD IN FRACTION ALONG UNIT CELL => X = CELL * V
C  SO V CAN BE OBTAINED AS V = RCELL * X
C  MEP operators work on V, not on X.

	WRITE(6,*) 'ENTER NAME OF .MEP FILE WITH SYM OPS'
	READ(5,40)FNAME
 40    FORMAT (A40)
	OPEN (UNIT=4,FILE=FNAME,STATUS='OLD')
	READ(4,*)NSYM
	DO 41 K=1,NSYM
  	READ(4,*) ((SYMOP(I,J,K),I=1,3),J=1,4)
	WRITE(6,*)
  	DO 41 I=1,3
	SYMOP(I,4,K)=SYMOP(I,4,K)/12.		!TRANSLATION STORED AS 12 * TRANSL
41	WRITE(6,42) (SYMOP(I,J,K),J=1,4)
42	FORMAT(12F8.4)
	CLOSE (UNIT=4)


50	write(6,*)'Start with (1)fractional or (2)catesian/Angstrom units?'
	write(6,*) '(anything else to quit)'
	read(5,*) iflag
	if (iflag.eq.1) then
	  write (6,*)' Enter fractional coord of site:'
	  read (5,*) (u(i),i=1,3)
	  goto 552
	else if (iflag.eq.2) then
	  write (6,*)' Enter Angstrom coord of site:'
	  read (5,*) (D(i),i=1,3)
	else 
	  goto 700
  	endif



c  CONVERT TO FRACTIONAL COORDINATES: D( ,L) -> U( )
 	do 550 i=1,3
	U(i)=0
	do 550 j=1,3
550	U(i)=U(i)+RCELL(i,j)*D(j)

552	continue

	DO 600 ksym=1,nsym
c  operate on atom with symmetry operator:  U( ) -> V( )
 	do 555 i=1,3
	V(i)=symop(i,4,ksym)
	do 554 j=1,3
554	V(i)=V(i)+symop(i,j,ksym)*U(j)
	IF (V(I).LT.-.5) V(I)=V(I)+1.0		!PUT IT BACK near origin (+/- .5)
555	IF (V(I).GT.0.5) V(I)=V(I)-1.0
c	IF (V(I).LT.0.) V(I)=V(I)+1.0		!PUT IT BACK IN THE CELL!
c555	IF (V(I).GT.1.) V(I)=V(I)-1.0

c  CONVERT TO ANGSTROM COORDINATES:         V( ) -> W() WHICH WILL BE COMPARED WITH ALL UNROTATED ATOMS
 	do 590 i=1,3
	W(i)=0
	do 590 j=1,3
590	W(i)=W(i)+CELL(i,j)*V(j)



600	write(6,601) (v(i),i=1,3),(w(i),i=1,3)
601	format(' ',3f10.6,'        ',3f10.3)
	goto 50

700	continue
	end
           
	SUBROUTINE MATPRN(ARRAY,M,N)
C	DOUBLE PRECISION ARRAY
	REAL*8 ARRAY(4,4)
	DO 10 I=1,M
10	WRITE (6,*) (ARRAY(I,J),J=1,N)
  	RETURN
        END


	SUBROUTINE MATINV (ARRAY, NORDER, DET)
	DOUBLE PRECISION ARRAY, AMAX, SAVE
2	DIMENSION ARRAY(4,4), IK (10), JK (10)
10	DET = 1
11 	DO 100  K=1, NORDER
C
C	FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX
C
	AMAX=0
21 	DO 30 I=K, NORDER
     	DO 30 J=K, NORDER
23 	IF (DABS (AMAX) - DABS (ARRAY(I,J) ) ) 24, 24, 30
24 	AMAX=ARRAY(I,J)
     	IK(K) =I
     	JK(K) =J
30 	CONTINUE
C
C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K)
C
31 	IF (AMAX) 41, 32, 41
32 	DET=0.
     	GO TO 140
41 	I=IK (K)
     	IF (I-K) 21, 51, 43
43 	DO 50 J=1, NORDER
     	SAVE= ARRAY(K,J)
      	ARRAY(K,J)= ARRAY(I,J)
50 	ARRAY(I,J)= -SAVE
51 	J=JK (K)
     	IF (J-K) 21, 61, 53
53 	DO 60 I=1, NORDER
     	SAVE= ARRAY(I,K)
     	ARRAY(I,K)= ARRAY(I,J)
60 	ARRAY(I,J)= -SAVE
C
C	ACCUMULATE ELEMENTS OF INVERSE MATRIX
C
61 	DO 70 I=1, NORDER
     	IF(I-K) 63, 70, 63
63 	ARRAY(I,K) = -ARRAY(I,K)  / AMAX
70 	CONTINUE
71 	DO 80 I=1, NORDER
     	DO 80 J=1, NORDER
     	IF (I-K) 74, 80, 74
74 	IF (J-K) 75, 80, 75
75 	ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K) *ARRAY(K,J)
80 	CONTINUE
81 	DO 90 J=1, NORDER
     	IF (J-K) 83, 90, 83
83 	ARRAY(K,J) = ARRAY(K,J) / AMAX
90 	CONTINUE
     	ARRAY(K,K) = 1. / AMAX
100 	DET=DET *AMAX
C
C 	RESTORE ORDERING OF MATRIX
C 
101 	DO 130 L=1, NORDER
       	K= NORDER  - L +1
       	J= IK (K)
       	IF (J-K) 111, 111, 105
105 	DO 110  I=1, NORDER
       	SAVE= ARRAY(I,K)
        ARRAY(I,K)= -ARRAY(I,J)
110 	ARRAY(I,J)= SAVE
111 	I= JK(K)
       	IF (I-K) 130, 130, 113
113 	DO 120  J=1, NORDER
       	SAVE= ARRAY(K,J)
       	ARRAY(K,J) = -ARRAY(I,J)
120 	ARRAY(I,J)= SAVE
130 	CONTINUE
140  	RETURN
	END