C     INTRAMOLECULAR AND INTERMOLECULAR DISTANCES
C                  CRYSTALLOGRAPHY LABORATORY, UNIV. OF PITTSBURGH
C                            R.SHIONO AND SHIRLEY S.C.CHU
C                                MAY 1964 AND FEBRUARY 1965
C              IBM 1130 SYSTEM WITH 16K CORE AND 1442 CARD READ/PUNCH
C                                       R.SHIONO
C                          APRIL 1968, FEB. 1969
C
	DIMENSION  C1(192),I1(192),J1(192),C2(192),I2(192),J2(192)
	DIMENSION  C3(192),I3(192),J3(192)
	DIMENSION  U(175),V(175),W(175),P(7)
	DIMENSION  NP(3,175),AIDA(175),AIDB(175)
	DIMENSION  TITLE(20),MSY(7),ORCD(3)
	DIMENSION  TST(4,250),KST(3,250)
	COMMON  ICDI,ICDO,ITYP,LDIS,A1,A2,A3,B1,B2,B3,TITLE
	COMMON  U,V,W,C1,C2,C3
	DATA MSY/'-Z','-Y','-X','  ','+X','+Y','+Z'/
C     CONSTANTS FOR ALPHABETIC +X,-X,+Y      ETC.  IN EBCDIC
	ICDI=5
	ICDO=7
	ITYP=6
	WRITE ( ITYP,91 )
	READ ( ICDI, 71 )  (TITLE(K),K=1,20 )
	WRITE  ( ITYP,71 )  TITLE
	WRITE ( ITYP,91 )
C     READ THE CELL CONSTANT CARD
	READ ( ICDI, 73) A1,A2,A3,B1,B2,B3
	WRITE ( ITYP,93 ) A1,A2,A3,B1,B2,B3
	JTRN=27
	READ (ICDI,78) N,M,IAGL,ITRN,ISOR,NCARD,LFMT,LSYM,DMAX,DMIN,DANG
	IF  ( ITRN )  87,  87,  86
86	JTRN=1
87	IAGL=IAGL+1
	NCARD=NCARD+1
C     READ SYMMETRY CARDS
	READ ( ICDI,75) (C1(K),I1(K),J1(K),C2(K),I2(K),J2(K),C3(K),
	1   I3(K),J3(K), K=1,M)
	IF  ( LSYM )  89,  89,  90
89	WRITE ( ITYP,79 )
	DO  61  L=1,M
	K1=I1(L)+4
	K2=J1(L)+4
	K3=I2(L)+4
	K4=J2(L)+4
	K5=I3(L)+4
	K6=J3(L)+4
61	WRITE (ITYP,80) L,C1(L),MSY(K1),MSY(K2),C2(L),MSY(K3),MSY(K4),
	1      C3(L),MSY(K5),MSY(K6)
90	IF ( DANG )  88,  88,  97
88	DANG=DMAX
97	DO  146  K=1,N
	IF ( LFMT )  150,  150,  151
150	READ ( ICDI,74 ) AIDA(K),AIDB(K),U(K),V(K),W(K)
	GO TO 152
C     X-RAY 67 FORMAT
151	READ ( ICDI,70 ) AIDA(K),AIDB(K),U(K),V(K),W(K)
152	GO TO ( 146, 147 ), NCARD
147	READ  ( ICDI, 74 )
146	CONTINUE
	WRITE ( ITYP, 8) DMAX,DMIN,DANG
	WRITE  ( ITYP,72 )
C     ADJUST THE COORDINATES TO WITHIN ONE HALF
	DO  39  L=1,N
	NP(1,L)=0
	NP(2,L)=0
	NP(3,L)=0
	IF ( ITRN )  180,  180,  39
180	IF ( ABS(U(L))-0.5 )  64,  64,  63
63	NP(1,L)=SIGN(1.0,U(L))
	IF   (U(L))   65,  64,  66
65	U(L)=U(L)+1.
	GO TO 64
66	U(L)=U(L)-1.
64	IF  ( ABS(V(L) ) -0.5 )  67,  67,  68
68	NP(2,L)=SIGN(1.0,V(L))
	IF   (V(L))   69,  67,  51
69	V(L)=V(L)+1.
	GO TO 67
51	V(L)=V(L)-1.
67	IF  ( ABS(W(L) ) -0.5 )  39,  39,  52
52	NP(3,L)=SIGN(1.0,W(L))
	IF   (W(L))   53,  39,  54
53	W(L)=W(L)+1.
	GO TO 39
54	W(L)=W(L)-1.
39	CONTINUE
	DMIN=DMIN*DMIN
	DMTWO=DMAX+DMAX
	DMA=DMTWO
	DMB=DMTWO
	DMC=DMTWO
	JA=1
	JB=1
	JC=1
	DMSQ=DMAX*DMAX
	IF  (B1)  43, 40, 43
40	JA=2
	IF  (B2)  48, 41, 48
41	JB=2
	DMC=DMAX
	IF  (B3)  9, 42, 9
42	DMB=DMAX
	DMA=DMAX
	JA=3
	GO TO 9
43	IF  (B2)  44, 46, 44
44	IF  (B3)  9, 50, 9
46	JB=2
	IF  (B3)  9, 47, 9
47	DMA=DMAX
	GO TO 50
48	IF  (B3)  9, 49, 9
49	DMB=DMAX
50	JC=2
C     LOOP FOR THE STARTING ATOM  1-N
9	DO  60  L=1,N
C     RECREATE THE ORIGINAL COORDINATE    FOR OUTPUT
	ORCD(1)=U(L)
	ORCD(2)=V(L)
	ORCD(3)=W(L)
	DO  145  II=1,3
	IF  ( NP(II,L) )     143,  145,  144
143	ORCD(II)=ORCD(II)-1.0
	GO TO 145
144	ORCD(II)=ORCD(II)+1.0
145	CONTINUE
	WRITE (ITYP,76) L,AIDA(L),AIDB(L),(ORCD(II),II=1,3)
	KDIS=0
C     FOR THE ENDING ATOMS   1-N
	DO  62  K=1,N
	P(4)=0.0
	P(5)=U(K)
	P(3)=-P(5)
	P(6)=V(K)
	P(2)=-P(6)
	P(7)=W(K)
	P(1)=-P(7)
C     GO THROUGH ALL SYMMETRY RELATIONS   1-M
	DO 62 I=1,M
	M1=I1(I)+4
	N1=J1(I)+4
	M2=I2(I)+4
	N2=J2(I)+4
	M3=I3(I)+4
	N3=J3(I)+4
	NPX=0
	NPY=0
	NPZ=0
	IJ=IABS(I1(I))
	IF  ( IJ )  111,  111,  110
110	NPX=NPX+NP(IJ,K)*ISIGN(1,I1(I))
111	IJ=IABS(J1(I))
	IF  ( IJ )  113,  113,  112
112	NPX=NPX+NP(IJ,K)*ISIGN(1,J1(I))
113	IJ=IABS(I2(I))
	IF  (  IJ )  115,  115,  114
114	NPY=NPY+NP(IJ,K)*ISIGN(1,I2(I))
115	IJ=IABS(J2(I))
	IF  ( IJ )  117,  117,  116
116	NPY=NPY+NP(IJ,K)*ISIGN(1,J2(I))
117	IJ=IABS(I3(I))
	IF  ( IJ )  119,  119,  118
118	NPZ=NPZ+NP(IJ,K)*ISIGN(1,I3(I))
119	IJ=IABS(J3(I))
	IF  ( IJ )  121,  121,  120
120	NPZ=NPZ+NP(IJ,K)*ISIGN(1,J3(I))
121	X=C1(I)+P(M1)+P(N1)
	Y=C2(I)+P(M2)+P(N2)
	Z=C3(I)+P(M3)+P(N3)
	DO  62  J=1,JTRN
	GO TO (11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30
	1,31,32,33,34,35,36,37),J
11	E=X
	R=Y
	G=Z
	NTRR=555
	GO TO 45
12	E=1.+X
	NTRR=655
	GO TO 45
13	E=X
	R=1.+Y
	NTRR=565
	GO TO 45
14	R=Y
	G=1.+Z
	NTRR=556
	GO TO 45
15	E=-1.+X
	G=Z
	NTRR=455
	GO TO 45
16	E=X
	R=-1.+Y
	NTRR=545
	GO TO 45
17	R=Y
	G=-1.+Z
	NTRR=554
	GO TO 45
18	E=1.+X
	R=1.+Y
	G=Z
	NTRR=665
	GO TO 45
19	E=X
	G=1.+Z
	NTRR=566
	GO TO 45
20	E=1.+X
	R=Y
	NTRR=656
	GO TO 45
21	E=-1.+X
	R=-1.+Y
	G=Z
	NTRR=445
	GO TO 45
22	E=X
	G=-1.+Z
	NTRR=544
	GO TO 45
23	E=-1.+X
	R=Y
	NTRR=454
	GO TO 45
24	E=1.+X
	R=-1.+Y
	G=Z
	NTRR=645
	GO TO 45
25	E=-1.+X
	R=1.+Y
	NTRR=465
	GO TO 45
26	E=X
	G=-1.+Z
	NTRR=564
	GO TO 45
27	R=-1.+Y
	G=1.+Z
	NTRR=546
	GO TO 45
28	E=-1.+X
	R=Y
	NTRR=456
	GO TO 45
29	E=1.+X
	G=-1.+Z
	NTRR=654
	GO TO 45
30	R=1.+Y
	G=1.+Z
	NTRR=666
	GO TO 45
31	G=-1.+Z
	NTRR=664
	GO TO 45
32	R=-1.+Y
	G=1.+Z
	NTRR=646
	GO TO 45
33	E=-1.+X
	R=1.+Y
	NTRR=466
	GO TO 45
34	E=1.+X
	R=-1.+Y
	G=-1.+Z
	NTRR=644
	GO TO 45
35	E=-1.+X
	R=1.+Y
	NTRR=464
	GO TO 45
36	R=-1.+Y
	G=1.+Z
	NTRR=446
	GO TO 45
37	G=-1.+Z
	NTRR=444
45	D4=0.
	D1=A1*(U(L)-E)
	IF  (D1-DMA)  55, 55, 62
55	D2=A2*(V(L)-R)
	IF  (D2-DMB)  56, 56, 62
56	D3=A3*(W(L)-G)
	IF  (D3-DMC)  57, 57, 62
57	GO TO (58, 59, 84), JA
58	D4=D4+D2*D3*B1
59	GO TO ( 81, 82), JB
81	D4=D4+D3*D1*B2
82	GO TO ( 83, 84), JC
83	D4=D4+D1*D2*B3
84	DIST=D1*D1+D2*D2+D3*D3+D4+D4
	IF  (DIST-DMSQ)  85, 85, 62
85	IF  ( DIST )  62,  122, 140
140	IF  ( DIST-DMIN )  62,  122,  122
122	DIST=SQRT(DIST)
	IF  ( NP(1,L) ) 130,  132,  131
130	E=E-1.0
	GO TO 132
131	E=E+1.0
132	IF  ( NP(2,L) )  133,  135,  134
133	R=R-1.0
	GO TO 135
134	R=R+1.0
135	IF  ( NP(3,L) )  136,  138,  137
136	G=G-1.0
	GO TO 138
137	G=G+1.0
138	NTRR=NTRR+(NP(1,L)-NPX)*100+(NP(2,L)-NPY)*10+NP(3,L)-NPZ
	GO TO ( 148, 148, 141 ), IAGL
C     IF STORAGE SPACE IS TIGHT, THIS SORTING PART MAY BE OMMITTED OR
C         LINKED
C     OMIT THE DISTANCE TO ITSELF
148	IF ( DIST )  173,  173,  174
173	IF ( NTRR-555 )  174,  175,  174
175	IF ( L-K )  174,  176,  174
176	IF ( I-1 )  174,  141,  174
174	IF ( ISOR )  142,  142,  149
142	MTX=C1(I)*100.
	MTY=C2(I)*100.
	MTZ=C3(I)*100.
	WRITE (ITYP,77) L,K,AIDA(K),AIDB(K),DIST,E,R,G,MTX,MSY(M1),MSY(N1)
	1,MTY,MSY(M2),MSY(N2),MTZ,MSY(M3),MSY(N3),I,NTRR
	GO TO 141
C     PREPARE FOR SORTING THE DISTANCES
149	KDIS=KDIS+1
	KST(1,KDIS)=K
	KST(2,KDIS)=I
	KST(3,KDIS)=NTRR
	TST(1,KDIS)=DIST
	TST(2,KDIS)=E
	TST(3,KDIS)=R
	TST(4,KDIS)=G
141	GO TO ( 62, 92, 92 ), IAGL
92	IF  (DIST-DANG)  98,  98,  62
98	WRITE (1)      L,K,DIST,E,R,G,I,NTRR,AIDA(K),AIDB(K)
62	CONTINUE
	IF ( ISOR )  60,  60,  160
C     START SORTING DATA
160	NDATA=KDIS
	IF  ( NDATA )  60,  60,  161
161	M1=NDATA
163	M1=M1/2
	IF ( M1 )  164,  170,  164
164	M2=NDATA-M1
	M3=1
165	M4=M3
166	M5=M4+M1
	IF ( TST(1,M4)-TST(1,M5) )  168,  168,  167
167	DO  169  IJ=1,4
	AHOLD=TST(IJ,M4)
	TST(IJ,M4)=TST(IJ,M5)
169	TST(IJ,M5)=AHOLD
	DO  172  IJ=1,3
	IHOLD=KST(IJ,M4)
	KST(IJ,M4)=KST(IJ,M5)
172	KST(IJ,M5)=IHOLD
	M4=M4-M1
	IF ( M4-1 )  168,  166,  166
168	M3=M3+1
	IF ( M3-M2 )  165,  165,  163
C     END OF SORTING
170	DO  171  IJ=1,NDATA
	I=KST(2,IJ)
	MTX=C1(I)*100.0
	MTY=C2(I)*100.0
	MTZ=C3(I)*100.0
	M1=I1(I)+4
	N1=J1(I)+4
	M2=I2(I)+4
	N2=J2(I)+4
	M3=I3(I)+4
	N3=J3(I)+4
	IA=KST(1,IJ)
171	WRITE (ITYP,77) L,IA,AIDA(IA),AIDB(IA),(TST(IK,IJ),IK=1,4),MTX,
	1 MSY(M1),MSY(N1),MTY,MSY(M2),MSY(N2),MTZ,MSY(M3),MSY(N3),I,
	2  KST(3,IJ)
60	CONTINUE
	GO TO ( 38, 99, 99), IAGL
99	L=0
	K=0
	DIST=0.
	WRITE (1      ) L,K,DIST,E,R,G,I,NTRR,AIDA(1),AIDB(1)
C     GO TO CALCULATE ANGLES
C              THIS SUBROUTINE MAY BE SEPARATED IF NOT ENOUGH STORAGE
	CALL ANGLE
38	WRITE  ( ITYP,94 )
	CALL EXIT
8	FORMAT ( / 5X5HDMAX=,F6.2,2X5HDMIN=,F6.2,3X5HDANG=,F6.2/)
70	FORMAT ( 7X,A4,A2,3F8.4 )
71	FORMAT (1X, 20A4 )
72	FORMAT( 2X5HATOMS,9X8HDISTANCE,6X1HX,7X1HY,7X1HZ,9X8HSYMMETRY,
	1  11X5HTRLN. /)
73	FORMAT ( 3F7.3, 3F7.4  )
74	FORMAT ( A4,A2,21X, 3F9.6  )
75	FORMAT  ( F11.6, 2I2,F11.6,2I2,F11.6, 2I2  )
76	FORMAT  ( /1X, I4, 2X A4,A2,3F10.5 / )
77	FORMAT (1X, I3,I4,1XA4,A2,F9.3,2X3F8.4,1X3(I3,2A2),1XI3,2XI3 )
78	FORMAT ( 8I4,3F6.2 )
79	FORMAT ( / 2X 3HNO.,24X8HSYMMETRY /)
80	FORMAT (1X, I4,F12.6,2XA3,A3,F12.6,2XA3,A3,F12.6,2XA3,A3 )
91	FORMAT ( / 1H0 )
93	FORMAT ( 1H0 3X 2HA=,F8.3,4H, B=,F8.3,4H, C=,F8.3,7H, COSA=,F7.4,
	1        7H, COSB=,F7.4,7H, COSC=,F7.4   ///    )
94	FORMAT  ( //1X,' END   OF DISTANCE CALCULATION' // )
	END
	SUBROUTINE ANGLE
C      CALCULATION OF ANGLES FROM DISTANCE PROGRAM RESULTS
C           CALCULATES ALL ANGLES, REJECTS REDUNDANT POSITIONS
C                                      BY R.SHIONO   APRIL 1963
C                            MODIFIED   MARCH 1965
C                  REVISED FOR NEW FORMAT       DEC. 1966
C             MODIFIED AS AN 1130 SUBROUTINE     APRIL 1968
C
	DIMENSION M(100),X(100),Y(100),Z(100),D(100),NS(100),NT(100)
	DIMENSION  TITLE(20)
	DIMENSION  AA(100),AB(100)
	COMMON  ICDI,ICDO,ITYP,LDIS,A,B,C,COSA,COSB,COSC,TITLE
C     THESE VALUES ARE IN COMMON TO SAVE STORAGE SPACE
	COMMON  X,Y,Z,D,M,NS,NT,AA,AB
C     SET THE STARTING DISK ADDRESS
	REWIND 1
	DATA KLP/'('/, KRP/')'/
	WRITE  (ITYP,4 )  TITLE
C     CONSTANTS FOR PUNCHING PARENTHESES TO BE PRINTED ON 407
C       EBCDIC FOR CHARACTERS TANSLATED INTO BINARY NUMBERS.
	READ (1)      N1,M(1),D(1),X(1),Y(1),Z(1),NS(1),NT(1),AA(1),AB(1)
41	N0=N1
	I=1
	GO TO 30
17	I=I+1
38	READ (1)      N1,M(I),D(I),X(I),Y(I),Z(I),NS(I),NT(I),AA(I),AB(I)
39	IF  (N1-N0)   11,  30,  11
30	IF  (N1-M(I))  10, 31,  10
31	IF  (NS(I)-1)  10, 32, 10
32	IF  ( NT(I)-555 )  10,  46,  10
C     CENTRAL ATOM WITH ZERO DISTANCE
46	WRITE (ITYP,5 ) N1,AA(I),AB(I),X(I),Y(I),Z(I)
	GO TO 38
10	IF  ( D(I) )  38,  38,  23
23	IF  ( I-1 )  17,  17,  34
34	II=I-1
	DO  35  N=1,II
	IF  (X(I)-X(N))  35, 36, 35
36	IF  (Y(I)-Y(N))  35, 37, 35
37	IF  (Z(I)-Z(N))  35, 38, 35
35	CONTINUE
	GO TO 17
11	K=I-1
	IF  (K-2)  18, 12, 12
12	K1=K-1
	WRITE ( ITYP,7 ) K
	DO  47  J=1,K
47	WRITE ( ITYP,6) M(J),AA(J),AB(J),KLP,NS(J),NT(J),KRP,X(J),Y(J),
	1  Z(J)
	WRITE ( ITYP,8 )  KLP,KRP,KLP,KRP
	DO  43  J=1,K1
	K2=J+1
	DO  14  L=K2,K
	A1 = (X(J) - X(L)) * A
	B1=(Y(J)-Y(L))*B
	C1=(Z(J)-Z(L))*C
	IF  (A1)  27, 25, 27
25	IF  (B1)  27, 26, 27
26	IF  (C1)  27, 14, 27
27	DS=A1*A1+B1*B1+C1*C1+2.*(A1*B1*COSC+B1*C1*COSA+C1*A1*COSB)
	AG=.5*(-DS+D(J)*D(J)+D(L)*D(L))/(D(J)*D(L))
	IF  (AG)  15,16,15
15	AGSQ=AG*AG
	IF  (1.-AGSQ)  28,  28,  29
29	AGL=ATAN(SQRT(1.0-AGSQ)/AG)*57.29578
	IF  (AGL)  13, 42, 22
13	AGL=180.+AGL
	GO TO 22
28	AGL=0.
42	DS=SQRT(DS)+0.001
	IF  (D(J)+D(L)-DS)  13, 13, 22
16	AGL=90.
22	WRITE ( ITYP,3) M(J),KLP,NS(J),NT(J),KRP,N0,M(L),KLP,NS(L),NT(L),
	1     KRP,D(J),D(L),AGL
14	CONTINUE
43	CONTINUE
18	IF  (N1)  19, 19, 40
40	M(1)=M(I)
	X(1)=X(I)
	Y(1)=Y(I)
	Z(1)=Z(I)
	D(1)=D(I)
	NS(1)=NS(I)
	NT(1)=NT(I)
	AA(1)=AA(I)
	AB(1)=AB(I)
	GO TO 41
19	RETURN
3	FORMAT (1X I4,1XA1,I3,I4,A1,I5,I6,1XA1,I3,I4,A1,F9.3,F8.3,F10.2 )
4	FORMAT (///1H1, 20A4, // 5X 23HANGLES AROUND EACH ATOM // )
5	FORMAT (/1X,16HCENTRAL ATOM NO.,I4,2XA4,A2,2X10HCOORDINATE,3F9.4/)
6	FORMAT ( 4XI4,2XA4,A2,1XA1,I3,I4,A1,1X3F9.5 )
7	FORMAT ( 3X23HATOMS AROUND THE CENTER, 10X12HTOTAL NUMBER, I4 / )
8	FORMAT (/2X6HATOM B,10X1HA,4X6HATOM C,10X1HD,A1,2HBA,A1,3X1HD,
	1  A1,2HAC,A1,5X5HANGLE   /     )
	END
