      SUBROUTINE BODAN1(MBODAN)
C
C BOND DISTANCE AND ANGLE SUBPROGRAM FOR U M REFINE
C  OUTPUT DELETED IF ARGUMENT NEGATIVE
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'A.INC'
      INCLUDE 'C.INC'
      INCLUDE 'E.INC'
      COMMON/D/B(9),DXX(3),DXY(3),DXZ(3),  F(6),XP(3),XXP(3),XYP(3),YP(3
     1),ATEMP(65),NUMB(300),XPR(3,300),LIST,SN,SX,NXX(6),ZTEMP(390)
      DIMENSION IDPNT(100),DLIST(100)
      EQUIVALENCE (NXX(1),NXL),(NXX(2),NYL),(NXX(3),NZL),(NXX(4),NXU),(N
     1 XX(5),NYU),(NXX(6),NZU)
C
   10 FORMAT(2X,2A3,3F8.5,43X,I5)
   12 FORMAT(1H1,18A4/15X,23HBOND DISTANCES BETWEEN F5.2,4H AND,F5.2,5H 
     1ANG.//2X,6HATOM 1,6X,1HX,7X,1HY,7X,1HZ,3X,6HATOM 2,6X,1HX,7X,1HY,
     2 7X,1HZ,4X,8HDISTANCE,2X,6HNUMBER/)
   14 FORMAT(34X,2A3,3F8.5,F11.5,I5)
   16 FORMAT(20H0ATOM ARRAY OVERFLOW)
   18 FORMAT(1H1,10X,18A4//10X,11HBOND ANGLES)
   20 FORMAT(1H0,5X,12HCENTRAL ATOM,5X,1HX,9X,1HY,9X,1HZ//9X,2A3,3F10.5/
     1/15X,6HATOM 1,2X,6HNUMBER,10X,6HATOM 2,2X,6HNUMBER,10X,5HANGLE/)
   22 FORMAT(15X,2A3,I6,12X,2A3,I6,F17.3)
C
      NUM=1
      IF(ICENT)310,300,310
  300 NUM=2
  310 IF(MBODAN)330,330,320
  320 WRITE(IOUT,12)TITLE,DMIN,DMAX
  330 SX=DMAX**2
      SN=DMIN**2
      LIST=0
			     DO 340 I=1,3
  340 DXX(I)=DMAX*RECELL(I)
C FIRST LOOP THROUGH ATOM LIST
			     DO 700 II=1,NATOM
			     DO 350 J=1,3
      DXY(J)=XYZ(J,II)+DXX(J)
  350 DXZ(J)=XYZ(J,II)-DXX(J)
      IF(II-1)360,400,360
C CHECK IF ATOM II IN GENERATED LIST
  360			     DO 390 I=1,LIST
      IF(NUMB(I)/64-II)390,370,390
  370			     DO 380 J=1,3
      IF(ABS(XPR(J,I)-XYZ(J,II))-0.00001)380,390,390
  380 CONTINUE
      GO TO 420
  390 CONTINUE
C ADD ATOM II TO LIST
  400 LIST=LIST+1
      I=LIST
      NUMB(LIST)=64*II+1
			     DO 410 J=1,3
  410 XPR(J,LIST)=XYZ(J,II)
  420 IF(MBODAN)440,430,430
  430 WRITE(IOUT,10)(TAG(J,II),J=1,2),(XYZ(J,II),J=1,3),I
C LOOP THROUGH GENERATED LIST
  440 NLST=0
			     DO 500 J=1,LIST
      S=0.0
			     DO 450 K=1,3
  450 XP(K)=XYZ(K,II)-XPR(K,J)
			     DO 470 N=1,3
      XZP=0.0
			     DO 460 K=1,3
  460 XZP=XZP+G(K,N)*XP(K)
  470 S=S+XZP*XP(N)
      IF((S-SX)*(SN-S))500,500,480
  480 S=SQRT(S)
      NLST=NLST+1
      IDPNT(NLST)=J
      DLIST(NLST)=S
  500 CONTINUE
C START SECOND LOOP THROUGH ATOM LIST
			     DO 690 JJ=1,NATOM
C START LOOP THROUGH SYMMETRY OPERATIONS
			     DO 690 J=1,NSYM
			     DO 510 L=1,3
      XP(L)=TRANS(L,J)
			     DO 510 M=1,3
  510 XP(L)=XP(L)+RMAT(L,M,J)*XYZ(M,JJ)
C LOOP FOR CENTER OF SYMMETRY
			     DO 690 ITY=1,NUM
      GO TO (540,520),ITY
  520			     DO 530 L=1,3
  530 XP(L)=-XP(L)
C LOOP THROUGH LATTICE POINTS
  540			     DO 690 KK=1,JCELL
			     DO 550 L=1,3
  550 XXP(L)=XP(L)+CELLP(L,KK)
C SET TRANSLATION LIMITS
			     DO 560 I=1,3
      NXX(I)=5-IFIX(XXP(I)-DXZ(I))
      NXX(I+3)=5.0+DXY(I)-XXP(I)
      IF(NXX(I)-NXX(I+3))560,560,690
  560 CONTINUE
      NXU=NXX(4)
      NYU=NXX(5)
      NZU=NXX(6)
C STEP THROUGH UNIT CELLS
			     DO 680 IA=NXL,NXU
      XYP(1)=XXP(1)+FLOAT(IA-5)
			     DO 680 IB=NYL,NYU
      XYP(2)=XXP(2)+FLOAT(IB-5)
			     DO 680 IC=NZL,NZU
      XYP(3)=XXP(3)+FLOAT(IC-5)
C CHECK FOR ATOM DUPLICATION
			     DO 590 I=1,LIST
      IF(NUMB(I)/64-JJ)590,570,590
  570			     DO 580 L=1,3
      IF(ABS(XYP(L)-XPR(L,I))-0.00001)580,590,590
  580 CONTINUE
      GO TO  680
  590 CONTINUE
C NEW ATOM, COMPUTE DISTANCE
      S=0.0
			     DO 600 K=1,3
  600 YP(K)=XYZ(K,II)-XYP(K)
			     DO 620 N=1,3
      XZP=0.0
			     DO 610 K=1,3
  610 XZP=XZP+YP(K)*G(K,N)
  620 S=S+XZP*YP(N)
      IF((S-SX)*(SN-S))680,680,630
  630 S=SQRT(S)
C CHECK IF LIST ARRAYS FULL
      IF(LIST-300)650,640,640
  640 WRITE(IOUT,16)
      GO TO 680
C ADD ATOM TO LIST
  650 LIST=LIST+1
      NUMB(LIST)=64*JJ+J
			     DO 660 K=1,3
  660 XPR(K,LIST)=XYP(K)
      NLST=NLST+1
      IDPNT(NLST)=LIST
      DLIST(NLST)=S
  680 CONTINUE
  690 CONTINUE
C BONDS GENERATED INTO LIST - SORT AND PRINT
      DO 698 J=1,NLST
      SMIN=DLIST(J)
      IJ=J
      DO 696 K=J+1,NLST
      IF(DLIST(K).LT.SMIN)THEN
	SMIN=DLIST(K)
	IJ=K
      END IF
  696 CONTINUE
      IF(IJ.NE.J)THEN
C SWAP ELEMENTS
	DLIST(IJ)=DLIST(J)
	DLIST(J)=SMIN
	K=IDPNT(IJ)
	IDPNT(IJ)=IDPNT(J)
	IDPNT(J)=K
      END IF
      JJ=IDPNT(J)
      IJ=NUMB(JJ)/64
      IF(MBODAN.GE.0)WRITE(IOUT,14)(TAG(K,IJ),K=1,2)
     1	       ,(XPR(K,JJ),K=1,3),SMIN,JJ
  698 CONTINUE
  700 CONTINUE
      IF(MBODAN.LE.1)RETURN
C COMPUTE BOND ANGLES
  720 WRITE(IOUT,18)TITLE
      IJ=LIST-1
C LOOP THROUGH ATOMS
			     DO 850 I=1,NATOM
      WRITE(IOUT,20)(TAG(K,I),K=1,2),(XYZ(J,I),J=1,3)
C FIRST LOOP THROUGH LIST
			     DO 850 J=1,IJ
			     DO 730 K=1,3
  730 XP(K)=XPR(K,J)-XYZ(K,I)
      S=0.0
			     DO 750 N=1,3
      XXP(N)=0.0
			     DO 740 K=1,3
  740 XXP(N)=XXP(N)+G(K,N)*XP(K)
  750 S=S+XXP(N)*XP(N)
      IF((S-SX)*(SN-S))850,850,760
  760 S=1.0/SQRT(S)
      NAT=NUMB(J)/64
      IK=J+1
C SECOND LOOP THROUGH LIST
			     DO 840 K=IK,LIST
      S2=0.0
			     DO 770 N=1,3
  770 YP(N)=XPR(N,K)-XYZ(N,I)
			     DO 790 N=1,3
      XZP=0.0
			     DO 780 M=1,3
  780 XZP=XZP+YP(M)*G(M,N)
  790 S2=S2+XZP*YP(N)
      IF((S2-SX)*(SN-S2))840,840,800
  800 S2=S/SQRT(S2)
      NTA=NUMB(K)/64
      CTH=0.0
			     DO 810 M=1,3
  810 CTH=CTH+XXP(M)*YP(M)
      CTH=CTH*S2
      TEMP=ABS(CTH)
      IF(TEMP-1.0)830,830,820
  820 CTH=CTH/TEMP
  830 CTH=180.0*ACOS(CTH)/PI
      WRITE(IOUT,22)(TAG(M,NAT),M=1,2),J,(TAG(M,NTA),M=1,2),K,CTH
  840 CONTINUE
  850 CONTINUE
      RETURN
      END
      SUBROUTINE ELVIB1(MELVIB,MBODAN)
C
C        THIS ROUTINE REVISED  JULY 28, 1970
C
C             ELLIPSOIDS OF VIBRATION CALCULATION
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'A.INC'
      INCLUDE 'C.INC'
      INCLUDE 'E.INC'
      COMMON/D/B(3,3),D(3,3),F(6)  ,X(3,3),XP(3),ATEMP(65),NUMB(300),
     1XPR(3,300),LIST,SN,SX,COEF(3),EIG(3),ROOT(3),YP(3),Y(3),ZTEMP(381)
C
   10 FORMAT(1H1,18A4//10X,23HELLIPSOIDS OF VIBRATION)
   12 FORMAT(///10X,14HELLIPSOID FOR 2A3//15X,24HEQUIVALENT ISOTROPIC B
     1= F8.4//20X,4HAXIS,10X,13HRMS AMPLITUDE,10X,12HANGLE WITH A,5X,12H
     2ANGLE WITH B,5X,12HANGLE WITH C/)
   14 FORMAT(///10X,2A3,25H IS NON-POSITIVE DEFINITE)
   16 FORMAT(20X,I3,12X,F11.6,3X,3F17.2)
   18 FORMAT(/12X,56HRMS DISPLACEMENTS AND THERMAL CORRECTIONS TO BONDS
     1FROM ,2A3,4H AT ,3F10.5	  /65X,9HDISTANCES,20X,12HDISPLACEMENT/
     2 15X,4HATOM,7X,1HX,8X,1HY,8X,1HZ,6X,3HRAW,4X,5HLOWER,5X,6HRIDING,
     3 2X,8HNON-CORR,2X,5HUPPER,3X,10HALONG BOND/)
   20 FORMAT(14X,2A3,9F9.5)
C
      WRITE(IOUT,10)TITLE
      TFC=0.5/PI**2
			     DO 300 I=1,3
      F(I)=G(I,I)
  300 Y(I)=SQRT(G(I,I))*RECELL(I)
      F(4)=2.0*G(1,2)
      F(5)=2.0*G(1,3)
      F(6)=2.0*G(2,3)
C	      START ATOM LOOP
			     DO 690 I=1,NATOM
      IF(ISOT(I).EQ.0)GO TO 480
C	      IGNORE ISOTROPIC ATOMS
      DO 320 K=1,2
	B(K,K)=BETA(K,I)
	IJ=K+1
	DO 310 L=IJ,3
	  IK=IJ+L
	  B(K,L)=BETA(IK,I)
	  B(L,K)=B(K,L)
310	CONTINUE
320   CONTINUE
      B(3,3)=BETA(3,I)
C	      B NOW CONTAINS MATRIX OF THERMAL COEF. OF ATOM I
C	      NOW COMPUTE D = G B
      DO 330 K=1,3
	DO 328 L=1,3
	  D(K,L)=0.0
	  DO 326 M=1,3
	    D(K,L)=D(K,L)+G(K,M)*B(M,L)
326	  CONTINUE
328	CONTINUE
330   CONTINUE
C	      SOLVE FOR EIGENVALUES OF D
      COEF(1)=D(1,1)+D(2,2)+D(3,3)
      COEF(2)=D(1,1)*D(2,2)+D(1,1)*D(3,3)+D(2,2)*D(3,3)
     1	     -D(1,2)*D(2,1)-D(1,3)*D(3,1)-D(2,3)*D(3,2)
      COEF(3)=D(1,1)*(D(2,2)*D(3,3)-D(2,3)*D(3,2))-D(1,2)*(D(2,1)*D(3,3)
     1 -D(3,1)*D(2,3))+D(1,3)*(D(2,1)*D(3,2)-D(3,1)*D(2,2))
C	      EIGENVALUE EQUATION IS
C	      0 = -X**3 + COEF(1)*X**2-COEF(2)*X + COEF(3)
      BISO=4.0/3.0*COEF(1)
      TEMP=COEF(2)-COEF(1)**2/3.0
      TEMP1=(-2.0*COEF(1)**3+9.0*COEF(1)*COEF(2)-27.0*COEF(3))/27.0
      TEMP2=-TEMP**3/27.0
      IF(0.25*TEMP1**2-TEMP2-1.0E-6.LE.0.0)GO TO 350
C	      ATOM NON-POSITIVE DEFINITE
  340 WRITE(IOUT,14)(TAG(K,I),K=1,2)
      GO TO 690
  350 IF(TEMP2.GT.1.0E-20)THEN
	TEMP2=-0.5*TEMP1/SQRT(TEMP2)
      ELSE
	TEMP2=0.0
      END IF
      TEMP3=ABS(TEMP2)
      IF(TEMP3.GT.1.0)TEMP2=TEMP2/TEMP3
      PHI=ACOS(TEMP2)/3.0
      TEMP=2.0*SQRT(ABS(TEMP/3.0))
      COEF(2)=TEMP*COS(PHI+2.0*PI/3.0)+COEF(1)/3.0
      COEF(3)=TEMP*COS(PHI-2.0*PI/3.0)+COEF(1)/3.0
      COEF(1)=TEMP*COS(PHI)+COEF(1)/3.0
C	      ARRANGE ROOTS IN ORDER OF INCREASING SIZE
      DO 390 N=1,2
	IK=N+1
	DO 380 L=IK,3
	  IF(COEF(N).gt.COEF(L))then
	    TEMP=COEF(N)
	    COEF(N)=COEF(L)
	    COEF(L)=TEMP
	  end if
  380	continue
  390 CONTINUE
C	      CALCULATE RMS DISPLACEMENTS
      DO 400 K=1,3
	IF(COEF(K).lt.0.0)go to 340
	ROOT(K)=SQRT(COEF(K)*TFC)
  400 continue
      WRITE(IOUT,12)(TAG(K,I),K=1,2),BISO
C	      START FINDING EIGENVECTORS
			     DO 410 K=1,3
			     DO 410 L=1,3
  410 X(K,L)=D(K,L)
			     DO 470 N=1,3
			     DO 420 L=1,3
  420 X(L,L)=D(L,L)-COEF(N)
C	      X IS NOW G B - LAMBDA I
      TEMP=1.0E-6
			     DO 430 ITMP=1,12
      TEMP=TEMP*0.01
			     DO 430 K=1,3
      M=K+1-3*(K/3)
      EIG(1)=RECELL(1)*(X(K,2)*X(M,3)-X(K,3)*X(M,2))
      EIG(2)=RECELL(2)*(X(K,3)*X(M,1)-X(K,1)*X(M,3))
      EIG(3)=RECELL(3)*(X(K,1)*X(M,2)-X(K,2)*X(M,1))
      FNC=EIG(1)**2+EIG(2)**2+EIG(3)**2
C	      CHECK FOR NULL VECTOR
      IF(FNC-TEMP)430,430,440
  430 CONTINUE
  440 FNC=SQRT(FNC+2.0*(EIG(1)*EIG(2)*RECELL(6)+EIG(1)*EIG(3)*RECELL(5)+
     1 EIG(2)*EIG(3)*RECELL(4)))
			     DO 460 K=1,3
      TEMP=FNC*Y(K)
      IF(ABS(TEMP).GT.1.0E-8)TEMP=EIG(K)/TEMP
      TEMP1=ABS(TEMP)
      IF(TEMP1-1.0)460,460,450
  450 TEMP=TEMP/TEMP1
  460 EIG(K)=180.0*ACOS(TEMP)/PI
C	      OUTPUT RESULTS
  470 WRITE(IOUT,16)N,ROOT(N),EIG
  480 IF(MELVIB-1)490,690,490
C	      COMPUTE RMS DISP. ALONG BONDS AND THERMAL CORRECTION TO BO
  490 IF(MBODAN)510,500,510
C	      BOND DISTANCES NOT AVAILABLE, COMPUTE AND SUPPRESS OUTPUT
  500 CALL BODAN1(-1)
      MBODAN=1
  510 WRITE(IOUT,18)(TAG(J,I),J=1,2),(XYZ(J,I),J=1,3)
C	      TRANSD IS MEAN SQUARE DISPLACEMENT OF ATOM I
      TRANSD=0.75*TFC*BISO
			     DO 680 J=1,LIST
      SUM=0.0
			     DO 520 K=1,3
  520 XP(K)=XPR(K,J)-XYZ(K,I)
			     DO 540 N=1,3
      YP(N)=0.0
			     DO 530 K=1,3
  530 YP(N)=YP(N)+G(K,N)*XP(K)
  540 SUM=SUM+YP(N)*XP(N)
      IF((SUM-SX)*(SN-SUM))680,680,550
  550 DIST=SQRT(SUM)
      IF(ISOT(I))570,560,570
  560 DISP=SQRT(ABS(BETA(1,I))*0.25*TFC)
      CORR1=0.50*ABS(BETA(1,I))*TFC
      GO TO 600
  570 PARA=0.0
			     DO 590 N=1,3
      XP(N)=0.0
			     DO 580 K=1,3
  580 XP(N)=XP(N)+YP(K)*B(N,K)
  590 PARA=PARA+XP(N)*YP(N)
      PARA=PARA*TFC/SUM
C	      PARA IS THE MEAN SQUARE DISPLACEMENT ALONG BOND
      DISP=SQRT(PARA)
      CORR1=ABS(TRANSD-PARA)
  600 K=NUMB(J)/64
      IF(ISOT(K))620,610,620
  610 CORR2=0.50*ABS(BETA(1,K))*TFC
      GO TO 670
  620 TEMP=0.0
			     DO 630 N=1,6
  630 TEMP=TEMP+F(N)*BETA(N,K)
      TEMP=TFC*ABS(TEMP)
      TEMP1=0.0
      L=MOD(NUMB(J),64)
			     DO 640 M=1,3
      COEF(M)=0.0
			     DO 640 N=1,3
  640 COEF(M)=COEF(M)+RMAT(N,M,L)*YP(N)
			     DO 650 M=1,3
  650 YP(M)=COEF(M)
			     DO 660 N=1,3
  660 TEMP1=TEMP1+YP(N)**2*BETA(N,K)
      TEMP1=TEMP1+2.0*(BETA(4,K)*YP(1)*YP(2)+BETA(5,K)*YP(1)*YP(3)+BETA(
     1 6,K)*YP(2)*YP(3))
      TEMP1=ABS(TEMP1)*TFC/SUM
      CORR2=ABS(TEMP1-TEMP)
  670 RIDE=DIST+0.5*ABS(CORR1-CORR2)/DIST
      UNCOR=DIST+0.5*(CORR1+CORR2)/DIST
      BL=DIST+0.5*(SQRT(CORR1)-SQRT(CORR2))**2/DIST
      BU=DIST+0.5*(SQRT(CORR1)+SQRT(CORR2))**2/DIST
      WRITE(IOUT,20)(TAG(L,K),L=1,2),(XPR(L,J),L=1,3),DIST,BL,RIDE,UNCOR
     1,BU,DISP
  680 CONTINUE
  690 CONTINUE
      RETURN
      END
