C ********** PROGRAM XHYDEX : APRIL 1986 : A. GUY ORPEN
C **********
C **********                               SCHOOL OF CHEMISTRY
C **********                               THE UNIVERSITY
C **********                               BRISTOL BS8 1TS
C **********                               U.K.
C **********
      INTEGER*2 IE(20000)
      INTEGER ED(10)
      COMMON EA(10),EB(10),EC(10),X1,X2,X3,COSA,COSB,COSC,HRAD(3),AX(3),
     *HSD(3),A1,A2,A3,XAT(150),YAT(150),ZAT(150),IR,ISRCH,ED,ILIST(150),
     *NSFAC(150)
      DIMENSION NATOM(150),ITIT(40),J(77),MN(46),HX(3),A(36),XTOL(3),
     *D(99),CW(20),FVAR(15),XH(26),YH(26),ZH(26),IAC(26),G(3),RX(3),
     *CX(20),CY(20),STEPMX(3),FMIN(26),HR(10,3),HS(10,3),XLOC(10),
     *YLOC(10),ZLOC(10),FLOC(10),NAME(10),ISR(10),ISC(30),AZ(9),
     *E(2809),AY(9),W(3),CR(10),P(9),AI(9),AP(3),XPLOT(3)
      DATA MN/4HTITL,4HCELL,4HTERM,4HEDGE,4HFACB,4HPARM,4HFVAR,4HEND ,
     *4HEND+,4H    ,4HHLOC,4HMATR,4HPLOT,4HCONT,4HAXIS,4HABSC,4HNLIM,
     *4HXLIM,4HSCAL,4HNUSF,4HFACE,4HHKLF,4HMERG,4HOMIT,4HEEES,4HTANG,
     *4HPHIX,4HWGHT,4HAFIX,4HANIS,4HUNDO,4HDFIX,4HL.S.,4HAFLS,4HBOND,
     *4HLIST,4HGRID,4HFMAP,4HNMAP,4HILSF,4HPLAN,4HFRPS,4HLATT,4HSYMM,
     *4HSFAC,4HUNIT/,IEND/0/,IA/0/,KL/0/,ILOC/0/,IHL/0/,
     *IF/1/,ISCR/0/,IC/0/,NN/3/,NIN/5/,NOUT/6/,MAXCAL/250/,G/3*0.E0/,
     *NAME/3HHL1,3HHL2,3HHL3,3HHL4,3HHL5,3HHL6,3HHL7,3HHL8,3HHL9,
     *4HHL10/,IFAIL/0/,
     *IPLOT/0/,P/1.,0.,0.,0.,1.,0.,0.,0.,1./,AY/9*0.0/,NC/7/,DC/1000./,
     *D/20.,50.,100.,200.,300.,500.,700.,92*10000./
     *,NJ/46/
C
      RADIAN=ATAN2(0.E0,-1.E0)
      RADIAN=RADIAN/1.8E2
      DO 94 I=1,9
      EA(I)=44.8E3
      EB(I)=2.04E0
      EC(I)=125.0E0
      CR(I)=0.77
94    ED(I)=6
      CR(10)=0.43
      EA(10)=6600.
      EB(10)=4.08
      EC(10)=49.2
      ED(10)=0
      DO 97 I=1,150
97    NSFAC(I)=1
996   ISRCH=0
      DO 500 I=1,3
500   HSD(I)=0.05E0
3     READ(NIN,100,END=99)LABL,(J(I),I=1,76)
100   FORMAT(A4,76A1)
      DO 1 K=1,NJ
      IF(LABL.EQ.MN(K)) GO TO 2
1     CONTINUE
      GO TO 14
13    WRITE(NOUT,101)
101   FORMAT(5X,11H**  IGNORED)
      GO TO 3
2     WRITE(NOUT,102)LABL,(J(I),I=1,76)
      IF(K.GT.14)GO TO 13
102   FORMAT(5X,12H INPUT .... ,A4,76A1)
      IF(K.EQ.1)GO TO 4
      CALL DECODE(A,N,J)
      GO TO (4,5,6,7,8,9,10,11,12,13,135,231,241,251),K
C  ***   TITL CARD  ***
4      DO 41 I=1,40
41    ITIT(I)=J(I)
      GO TO 3
C  ***  CELL CARD  ***
5     IF(N.GE.6 .AND. N.LT.8) GO TO 51
98    WRITE(NOUT,103)LABL
103   FORMAT(5X,7H ** BAD,2X,A4,5H CARD)
      GO TO 999
51    I=1
      IF(N.EQ.7)I=2
      AX(1)=A(I)
      AX(2)=A(I+1)
      AX(3)=A(I+2)
      COSA=A(I+3)*RADIAN
      COSB=A(I+4)*RADIAN
      COSC=A(I+5)*RADIAN
      A1=SIN(COSA)
      A2=SIN(COSB)
      A3=SIN(COSC)
      COSA=COS(COSA)
      COSB=COS(COSB)
      COSC=COS(COSC)
      X1=(COSA*COSB-COSC)/(A1*A2)
      X2=(1.-X1*X1)
      X2=SQRT(X2)
      AZ(1)=AX(1)*A2*X2
      AZ(2)=0.
      AZ(3)=0.
      AZ(4)=-AX(1)*A2*X1
      AZ(5)=AX(2)*A1
      AZ(6)=0.
      AZ(7)=AX(1)*COSB
      AZ(8)=AX(2)*COSA
      AZ(9)=AX(3)
      X1=2.E0*AX(2)*AX(3)*COSA
      X2=2.E0*AX(1)*AX(3)*COSB
      X3=2.E0*AX(1)*AX(2)*COSC
      A1=AX(1)*AX(1)
      A2=AX(2)*AX(2)
      A3=AX(3)*AX(3)
      GO TO 3
C  ***  TERM CARD - SEARCH FOR TERMINAL HYDRIDE  ****
6     ISRCH=1
104   FORMAT(1H1,5X,28H SEARCH FOR TERMINAL HYDRIDE,2X,40A1//)
      HRAD(1)=1.6
      IF(N.LT.1) GO TO 98
      ILIST(1)=INT(A(1))
      IF(N.GT.1)HRAD(1)=A(2)
      IF(N.GT.2)HSD(1)=A(3)
      IF(N.GE.4)MAXCAL=INT(A(4))
      GO TO 3
C  ***  EDGE CARD SEARCH FOR EDGE-BRIDGING HYDRIDE  ***
7     IF(N.LT.2)GO TO 98
      HRAD(2)=1.85
      HRAD(1)=1.85
106   FORMAT(1H1,5X,33H SEARCH FOR EDGE BRIDGING HYDRIDE,2X,40A1//)
      ISRCH=2
      ILIST(1)=INT(A(1))
      ILIST(2)=INT(A(2))
      IF(N.GT.2)HRAD(1)=A(3)
      HRAD(2)=HRAD(1)
      IF(N.GE.4)HRAD(2)=A(4)
      IF(N.GE.5)HSD(1)=A(5)
      HSD(2)=HSD(1)
      IF(N.GE.6)HSD(2)=A(6)
      IF(N.GE.7)MAXCAL=INT(A(7))
      GO TO 3
C  ***  FACB CARD - SEARCH FOR FACE -BRIDGING HYDRIDE  ***
8     IF(N.LT.3) GO TO 98
107   FORMAT(1H1,5X,33H SEARCH FOR FACE BRIDGIN HYDRIDE,2X,40A1//)
      DO 1000 I=1,3
1000  HRAD(I)=1.95
      ISRCH=3
      ILIST(1)=INT(A(1))
      ILIST(2)=INT(A(2))
      ILIST(3)=INT(A(3))
      IF(N.GE.4)HRAD(1)=A(4)
      HRAD(2)=HRAD(1)
      IF(N.GE.5)HRAD(2)=A(5)
      HRAD(3)=HRAD(2)
      IF(N.GE.6)HRAD(3)=A(6)
      IF(N.GE.7)HSD(1)=A(7)
      HSD(2)=HSD(1)
      IF(N.GE.8)HSD(2)=A(8)
      HSD(3)=HSD(1)
      IF(N.GE.9)HSD(3)=A(9)
      IF(N.GE.10)MAXCAL=INT(A(10))
      GO TO 3
C  ***  PARM CARD - INPUT P.E.PARAMETERS ***
9     IF(N.LT.5) GO TO 13
      IZ=INT(A(1))
      EA(IZ)=A(2)
      EB(IZ)=A(3)
      EC(IZ)=A(4)
      ED(IZ)=INT(A(5))
      IF(N.GT.5)CR(IZ)=A(6)
      GO TO 3
C  ***  FVAR CARD  ***
10      IF(N.LT.2) GO TO 3
      IF=N
      IF(IF.GT.15)IF=15
      DO91 I=1,IF
93      IF(ABS(A(I)).LT.5.E0) GO TO 92
      A(I)=A(I)-10.E0
      GO TO 93
92    FVAR(I)=A(I)
91      CONTINUE
      GO TO 3
C  ***   HLOC CARD (LOCATED FROM PREVIOUS JOB(S)  ***
135   IF(N.LT.5) GO TO 98
      IF(ILOC.GE.10)GO TO 98
      IZ=INT(A(4))
      IF(IZ.LE.0)GO TO 98
      IF(N.LT.4+IZ)GO TO 98
      ILOC=ILOC+1
      KL=ILOC
       DO 1353 I=1,3
      HS(ILOC,I)=0.05
1353  HR(ILOC,I)=1.85
      IA=IA+1
      IF(IA.GE.150)GO TO 98
      XAT(IA)=A(1)
      YAT(IA)=A(2)
      ZAT(IA)=A(3)
      NATOM(IA)=NAME(ILOC)
      NSFAC(IA)=10
      ISR(ILOC)=IZ
      DO 1351 I=1,IZ
      ISCR=ISCR+1
1351  ISC(ISCR)=INT(A(4+I))
      IF(N.GE.IZ+IZ+5)HS(ILOC,1)=A(IZ+IZ+5)
      HS(ILOC,2)=HS(ILOC,1)
      IF(N.GE.IZ+IZ+6)HS(ILOC,2)=A(IZ+IZ+6)
      HS(ILOC,3)=HS(ILOC,2)
      IF(N.GE.IZ+IZ+7)HS(ILOC,3)=A(IZ+IZ+7)
      DO 1354 I=1,IZ
      IF(N.GE.IZ+I+4)HR(ILOC,I)=A(IZ+I+4)
      IF(I.EQ.1)GO TO 1354
      IF(N.LT.IZ+I+4)HR(ILOC,I)=HR(ILOC,I-1)
1354  CONTINUE
      GO TO 3
C ***   MATR CARD ***
231   IF(N.LT.9)GO TO 98
      DO 232 I=1,9
232   P(I)=A(I)
      AY(1)=P(1)*AZ(1)+P(2)*AZ(4)+P(3)*AZ(7)
      AY(2)=P(1)*AZ(2)+P(2)*AZ(5)+P(3)*AZ(8)
      AY(3)=P(1)*AZ(3)+P(2)*AZ(6)+P(3)*AZ(9)
      AY(4)=P(4)*AZ(1)+P(5)*AZ(4)+P(6)*AZ(7)
      AY(5)=P(4)*AZ(2)+P(5)*AZ(5)+P(6)*AZ(8)
      AY(6)=P(4)*AZ(3)+P(5)*AZ(6)+P(6)*AZ(9)
      AY(7)=P(7)*AZ(1)+P(8)*AZ(4)+P(9)*AZ(7)
      AY(8)=P(7)*AZ(2)+P(8)*AZ(5)+P(9)*AZ(8)
      AY(9)=P(7)*AZ(3)+P(8)*AZ(6)+P(9)*AZ(9)
      CALL INV(AY,AI,DET)
      IF(ABS(DET).LT.0.00001)GO TO 98
      GO TO 3
C  ***  PLOT  CARD  ***
241   IF(N.LT.5)GO TO 98
      IPLOT=1
      DO 242 I=1,3
242   XPLOT(I)=A(I)
      ISRCH=INT(A(4))
      IF(ISRCH.LE.0)GO TO 98
      IF(N.LT.4+ISRCH)GO TO 98
      DO 245 I=1,ISRCH
245   ILIST(I)=INT(A(4+I))
      IR=ISRCH
      DO 246 I=1,IA
      DO 247 IJ=1,ISRCH
      IF(ILIST(IJ).EQ.I)GO TO 246
247   CONTINUE
      IR=IR+1
      ILIST(IR)=I
246   CONTINUE
      HRAD(1)=1.85
      HSD(1)=0.05
      IF(N.GE.5+ISRCH+ISRCH)HSD(1)=A(6+ISRCH+ISRCH)
      HSD(2)=HSD(1)
      IF(N.GE.6+ISRCH+ISRCH)HSD(2)=A(ISRCH+ISRCH+7)
      HSD(3)=HSD(2)
      IF(N.GE.7+ISRCH+ISRCH)HSD(3)=A(ISRCH+ISRCH+8)
      DO 1352 I=1,ISRCH
      IF(N.GE.ISRCH+I+5)HRAD(I)=A(I+ISRCH+5)
      IF(I.EQ.1)GO TO 1352
      IF(N.LT.ISRCH+I+5)HRAD(I)=HRAD(I-1)
1352  CONTINUE
      SCAL=1.0
      SP=3.0
      IF(N.GT.4+ISRCH)SP=A(5+ISRCH)
      IF(SP.GT.4.0)SCAL=4./SP
      GO TO 3
C  ***   CONT CARD  ***
251   IF(N.LT.3)GO TO 98
      NC=INT(A(1))
      IF(N.LT.2+NC)GO TO 98
      DC=A(2)
      DO 243 I=1,NC
      D(I)=A(I+2)
243   CONTINUE
      GO TO 3
C  ***   END  CARD  ***
11    IEND=0
      GO TO 99
C  ***  END+ CARD  ***
12    IEND=1
      GO TO 99
C  ***  ATOM CARDS INPUT  ***
14    WRITE(NOUT,102)LABL,(J(I),I=1,76)
      CALL DECODE(A,N,J)
      IF(N.LT.4) GO TO 13
      IA=IA+1
      IF(IA.GT.150) GO TO 98
      NATOM(IA)=LABL
      DO 141 I=2,4
      IF(ABS(A(I)).LT.5.E0) GO TO 141
      Z=A(I)/10.E0
      IZ=INT(Z)
      IF(IABS(IZ).GT.IF) GO TO  98
      ZA=0.E0
      IF(A(I).LT.0.E0) ZA=1.E0
      ZP=A(I)-10.E0*FLOAT(IZ)
      IZ=IABS(IZ)
      A(I)=ZP*(FVAR(IZ) -ZA)
141   CONTINUE
      NSFAC(IA)=INT(A(1))
      IF(NSFAC(IA).GE.10)GO TO 13
      XAT(IA)=A(2)
      YAT(IA)=A(3)
      ZAT(IA)=A(4)
      GO TO 3
C  ***  SEARCH SECTION OF PROGRAM  ***
C  ***  REFINEMENT OF HLOC POSITIONS
99    IF(IPLOT.GT.0)GO TO 9998
      IF(ILOC.LE.0)GO TO 136
      ISCR=0
      DO 138 I=1,3
      XTOL(I)=0.001/AX(I)
138   STEPMX(I)=XTOL(I)*5.
      FEST=0.0
      WRITE(NOUT,137)
      WRITE(NOUT,133)
991   IR=0
      IC=IC+1
      IF(IC.GT.ILOC*5)GO TO 144
      RSRCH=4.0
      IHL=IHL+1
      IF(IHL.LE.ILOC)GO TO 1393
      IHL=1
      ISCR=0
1393  ISRCH=ISR(IHL)
      DO 1392 I=1,ISRCH
      ISCR=ISCR+1
      IR=IR+1
      ILIST(IR)=ISC(ISCR)
      HSD(I)=HS(IHL,I)
      HRAD(I)=HR(IHL,I)
1392  CONTINUE
      IZ=IA-ILOC
      IZL=IZ+IHL
      DO 139 I=1,IZ
      CALL RAD(XAT(I),YAT(I),ZAT(I),XAT(IZL),YAT(IZL),ZAT(IZL),R)
      DO 1395 IJ=1,ISRCH
      IF(I.NE.ILIST(IJ))GO TO 1395
      IF(R.GT.HRAD(IJ)+0.45)GO TO 1391
      GO TO 139
1395  CONTINUE
      IF(R.GT.RSRCH)GO TO 139
      IR=IR+1
      ILIST(IR)=I
139   CONTINUE
      GO TO 1396
1391  WRITE(NOUT,134)IHL
134   FORMAT(//5X,8H HYDRIDE,I4,20H  NOT CLOSE TO METAL)
      STOP
1396   DO 143 I=1,ILOC
      IF(I.EQ.IHL)GO TO 143
      IR=IR+1
      ILIST(IR)=I+IZ
143   CONTINUE
      HX(1)=XAT(IZL)
      HX(2)=YAT(IZL)
      HX(3)=ZAT(IZL)
      CALL MIN(NN,HX,F,G,XTOL,STEPMX,FEST,MAXCAL,IFAIL)
      XAT(IZL)=HX(1)
      YAT(IZL)=HX(2)
      ZAT(IZL)=HX(3)
      WRITE(NOUT,132)IHL,HX(1),HX(2),HX(3),F
      IF(IC.LE.ILOC*4)GO TO 991
      NF=2
      CALL FUNCT(NN,NF,HX,F,G)
      FLOC(IHL)=F
      GO TO 991
144   DO 145  I=1,ILOC
      IZL=I+IZ
      XLOC(I)=XAT(IZL)
      YLOC(I)=YAT(IZL)
      ZLOC(I)=ZAT(IZL)
      HX(1)=XLOC(I)
      HX(2)=YLOC(I)
      HX(3)=ZLOC(I)
145   CONTINUE
      GO TO 9998
C  ****   INITIAL SEARCH SECTION ***
136   K=1
      DO 95 I=1,IA
      IF(NSFAC(I).GT.10)GO TO 95
      K=MAX0(K,NSFAC(I))
95    CONTINUE
      IF(ISRCH.EQ.0) GO TO 199
      RSRCH=4.0
      IF(ISRCH.EQ.1)WRITE(NOUT,104)ITIT
      IF(ISRCH.EQ.2)WRITE(NOUT,106)ITIT
      IF(ISRCH.EQ.3)WRITE(NOUT,107)ITIT
      DO 165 I=1,3
165   RX(I)=0.E0
      DO 151 I=1,ISRCH
      IR=ILIST(I)
      RX(1)=XAT(IR)+RX(1)
      RX(2)=YAT(IR)+RX(2)
      RX(3)=RX(3)+ZAT(IR)
151   CONTINUE
      IF(ISRCH.EQ.1) GO TO158
      DO 161  I=1,3
161   RX(I)=RX(I)/FLOAT(ISRCH)
      DO 166  I=1,ISRCH
      IZ=ILIST(I)
      CALL RAD(RX(1),RX(2),RX(3),XAT(IZ),YAT(IZ),ZAT(IZ),R)
      IF(R.LT.(HRAD(I)+5.*HSD(I)))GO TO 166
      WRITE(NOUT,120)NATOM(IZ)
120   FORMAT(5X,A4,2X,24H IS IMPOSSIBLE TO BRIDGE)
      GO TO 199
166   CONTINUE
158   IR=ISRCH
      DO 152  I=1,IA
      DO 171 IZ=1,ISRCH
      IF(I.EQ.ILIST(IZ)) GO TO 152
171   CONTINUE
      CALL RAD(RX(1),RX(2),RX(3),XAT(I),YAT(I),ZAT(I),R)
      IF(R.GT.RSRCH) GO TO 152
      IR=IR+1
      ILIST(IR)=I
      IF(IR.LE.15) GO TO 152
      RSRCH=RSRCH-0.25E0
      GO TO 158
152   CONTINUE
      IF(IR.GT.ISRCH+3) GO TO 156
      RSRCH=RSRCH+0.2E0
      GO TO 158
156   DO 181 I=1,3
      XTOL(I)=0.02/AX(I)
181   STEPMX(I)=XTOL(I)*5.E0
      FEST=0.E0
C  ***  GENERATE INITIAL HYDRIDE POSITIONS
      DO 211 I=1,26
211   IAC(I)=0
      Z=1.E0
      K=0
      NH=0
      IH=0
163   DO 153 I=1,13
      DO 154 IZ=1,3
      HX(IZ)=RX(IZ)
      IF(IZ.EQ.I)HX(IZ)=HX(IZ)+Z*RSRCH/AX(I)
154   CONTINUE
      NH=NH+1
      IF(I.GT.9)GO TO 202
      IF(I.GT.3)GO TO 201
      GO TO 203
201   IZ=I-3
      IF(IZ.GT.3)GO TO 204
      DO 205 IJ=1,3
      IF(IJ.NE.IZ)HX(IJ)=HX(IJ)+Z*RSRCH/(AX(IJ)*SQRT(2.E0))
205   CONTINUE
      GO TO 203
204   IZ=IZ-3
      DO 206 IJ=1,3
      IF(IJ.EQ.IZ)GO TO 206
      Z=-Z
      HX(IJ)=HX(IJ)+Z*RSRCH/(AX(IJ)*SQRT(2.E0))
206   CONTINUE
      GO TO 203
202   IZ=I-9
      DX=Z*RSRCH/SQRT(3.E0)
      Z1=DX
      IF(IZ.EQ.2)Z1=-Z1
      IF(IZ.EQ.3)DX=-DX
      IF(IZ.LT.4)GO TO 207
      Z1=-Z1
      DX=-DX
207   HX(1)=HX(1)+Z*RSRCH/(SQRT(3.E0)*AX(1))
      HX(2)=HX(2)+Z1/AX(2)
      HX(3)=HX(3)+DX/AX(3)
203   CALL MIN(NN,HX,F,G,XTOL,STEPMX,FEST,MAXCAL,IFAIL)
      IH=IH+1
      IF(IFAIL.NE.0)IAC(IH)=2
111   FORMAT(//5X,29H MINIMISED POTENTIAL ENERGY =,F10.5,
     */5X,18H HYDRIDE POSITION ,3(F9.5,2X))
114   FORMAT(//5X,39H H - OTHER ATOM DISTANCES FOR HYDRIDE #,I4)
      XH(IH)=HX(1)
      YH(IH)=HX(2)
      ZH(IH)=HX(3)
      FMIN(IH)=F
153   CONTINUE
      IF(Z.LT.0.E0) GO TO 998
      Z=-Z
      GO TO 163
C  ***  EXAMINE LIST OF SOLUTIONS ***
998    WRITE(NOUT,126)RSRCH,(HRAD(I),HSD(I),I=1,ISRCH)
      DO 196 I=1,3
       XTOL(I)=0.1*XTOL(I)
196   STEPMX(I)=STEPMX(I)*0.1
      K=0
126   FORMAT(//5X,42H ATTEMPTED HOLE LOCATION :SEARCH RADIUS = ,F6.3
     *,22H : M-H BOND LENGTHS = ,3(F6.3,2H (,F4.3,2H )),//5X,8H ATTEMPT,
     *5X,2H X,10X,2H Y,10X,2H Z,7X,8H MINIMUM,13X,9H COMMENTS,/)
      Z=1.0
      DO 11010 I=1,ISRCH
11010 Z=AMIN1(Z,HRAD(I))
      DO 155 I=1,IH
      IF(FMIN(I).GT.999.)FMIN(I)=999.
      IF(IAC(I).NE.0)GO TO 213
      DO 157 IK=1,IR
      IF(IAC(I).NE.0)GO TO 213
      IZ=ILIST(IK)
      CALL RAD(XH(I),YH(I),ZH(I),XAT(IZ),YAT(IZ),ZAT(IZ),R)
      IF(IK.GT.ISRCH) GO TO 212
      IF(ABS(R-HRAD(IK)).LE.5.0*HSD(IK))GO TO 157
      IAC(I)=1
      NBAD=NATOM(IZ)
      GO TO 157
212   IF(R.GE.(Z+0.15))GO TO 157
      IAC(I)=1
      NBAD=NATOM(IZ)
157   CONTINUE
      IF(IAC(I).EQ.0)K=K+1
213   IF(IAC(I).EQ.1)WRITE(NOUT,127)I,XH(I),YH(I),ZH(I),FMIN(I),NBAD,R
      IF(IAC(I).EQ.0)WRITE(NOUT,128)I,XH(I),YH(I),ZH(I),FMIN(I)
      IF(IAC(I).EQ.2)WRITE(NOUT,125)I,XH(I),YH(I),ZH(I),FMIN(I)
155   CONTINUE
127   FORMAT(5X,I5,4(4X,F8.4),10X,' BAD',2X,A4,2X,' DISTANCE =',F8.4)
128   FORMAT(5X,I5,4(4X,F8.4),10X,' **  GOOD SOLUTION')
125   FORMAT(5X,I5,4(4X,F8.4),10X,' BAD CONVERGENCE - UNACCEPTABLE')
115   FORMAT(5X,8H H......,A4,5X,F8.4)
117   FORMAT(5X,8H H------,A4,5X,F8.4)
      WRITE(NOUT,122)K
      IF(K.LE.0)GO TO 199
      NH=1
      RSRCH=RSRCH+0.5
997   DO 190 I=1,26
      IF(IAC(I).NE.0)GO TO 190
      Z=FMIN(I)
      IZ=I
      GO TO 214
190   CONTINUE
214   DO 191 I=1,IH
      IF(IAC(I).NE.0) GO TO 191
      IF(Z.LE.FMIN(I))GO TO 191
      Z=FMIN(I)
      IZ=I
191   CONTINUE
      HX(1)=XH(IZ)
      HX(2)=YH(IZ)
      HX(3)=ZH(IZ)
      DO 210 I=1,26
      CALL RAD(HX(1),HX(2),HX(3),XH(I),YH(I),ZH(I),R)
      IF(R.GE.0.25)GO TO 210
      IAC(I)=1
      K=K-1
210   CONTINUE
      DO 193 I=1,IA
      DO 194 IZ=1,IR
      IF(I.EQ.ILIST(IZ))GO TO 193
194   CONTINUE
      CALL RAD(HX(1),HX(2),HX(3),XAT(I),YAT(I),ZAT(I),R)
      IF(R.GT.RSRCH)GO TO 193
      IR=IR+1
      ILIST(IR)=I
193   CONTINUE
      CALL MIN(NN,HX,F,G,XTOL,STEPMX,FEST,MAXCAL,IFAIL)
      WRITE(NOUT,111)F,HX
      NF=2
      CALL FUNCT(NN,NF,HX,F,G)
124   FORMAT(1H1)
      KL=KL+1
      IF(KL.LE.10)GO TO 1242
      WRITE(NOUT,1243)
1243  FORMAT(47H SORRY, TEN HOLES ALREADY LOCATED, STOPPING NOW)
1242  XLOC(KL)=HX(1)
      YLOC(KL)=HX(2)
      ZLOC(KL)=HX(3)
      FLOC(KL)=F
      DO 1244 I=1,3
      HS(KL,I)=HSD(I)
1244  HR(KL,I)=HRAD(I)
      IF(K.LE.0)GO TO 199
      NH=NH+1
      GO TO 997
199   IF(KL.LE.0)GO TO 244
      WRITE(NOUT,1241)NH
      GO TO 236
1241  FORMAT(//5X,I4,35H UNIQUE HYDRIDE POSITION(S) LOCATED,/1H1)
137   FORMAT(//5X,34H JOINT HYDRIDE LOCATION REFINEMENT)
131   FORMAT(//5X,30H REFINED HYDRIDE POSITION DATA)
133   FORMAT(//5X,3H H#,10X,4H X/A,7X,4H Y/B,7X,4H Z/C,7X,5H P.E.)
132   FORMAT(5X,I3,6X,4(2X,F9.4))
244   IF(IPLOT.EQ.0)STOP
248   XP=XPLOT(1)*AY(1)+XPLOT(2)*AY(2)+XPLOT(3)*AY(3)
      YP=XPLOT(1)*AY(4)+XPLOT(2)*AY(5)+XPLOT(3)*AY(6)
      ZP=XPLOT(1)*AY(7)+XPLOT(2)*AY(8)+XPLOT(3)*AY(9)
      AP(3)=ZP
      FMAX=1000.1
      IK=0
      IJ=0
      DO 238 I=1,53
      IJ=IJ+1
      AP(1)=XP+FLOAT(IJ)*(SP/26.)-SP
      KJ=0
      DO 238 K=1,53
      NF=3
      KJ=KJ+1
      IK=IK+1
      AP(2)=YP+FLOAT(KJ)*(SP/26.)-SP
      HX(1)=AP(1)*AI(1)+AP(2)*AI(2)+AP(3)*AI(3)
      HX(2)=AP(1)*AI(4)+AP(2)*AI(5)+AP(3)*AI(6)
      HX(3)=AP(1)*AI(7)+AP(2)*AI(8)+AP(3)*AI(9)
      CALL FUNCT(NN,NF,HX,F,G)
      E(IK)=F
      IF(ABS(F).LT.FMAX)FMAX=ABS(F)
238   CONTINUE
      Z=9.1/FMAX
      WRITE(NOUT,230)Z
230   FORMAT(/5X,19H P.E. MULTIPLIED BY ,F9.4)
      L=1
      I=1
      DO 910 IJ=1,53
      DO 910 K=1,53
      E(I)=E(I)*Z
      IF(E(I).GT.999.1)E(I)=999.1
      IE(L)=E(I)
      IF(MOD(IJ,2)+MOD(K,2).EQ.0)L=L+1
910   I=I+1
      WRITE(6,912)(IE(I),I=1,676)
912    FORMAT(1H1,(//1X,26(1X,I3)))
      WRITE(NOUT,233)SCAL
C  INITIALISE PLOTTER
C     CALL BRKPLT(8,300.)
      SCAL=SP*2.*SCAL
      DO 911 I=1,4
      CALL XPLOD(SCAL,2,0.,0.)
      CALL XPLOD(SCAL,0,0.,SCAL)
      CALL XPLOD(SCAL,0,SCAL,SCAL)
      CALL XPLOD(SCAL,0,SCAL,0.)
      CALL XPLOD(SCAL,0,0.,0.)
911   CONTINUE
233   FORMAT(//5X,12H PLOT SCALE ,F6.4,18H   INCHES/ANGSTROM)
C
C INTERPOLATE
C
      CX(10)=0.
      DO 22 I=1,10
22    CY(I)=0.
      DO 914 I=3,5
      CY(I)=100.
914   CY(I+4)=-100.
      CY(1)=-100.
      DO 915 I=1,9
915   CX(I)=CY(10-I)
      DO 30 M=1,NC
      NA=-1
      N=54
16    N=N+1
      IF(N.GT.2701)GOTO 21
      IF(MOD(N+1,53).LT.3)GOTO 16
      CW(1)=E(N)
      CW(3)=E(N+1)
      CW(5)=E(N+54)
      CW(7)=E(N+53)
      CW(9)=CW(1)
      CW(2)=0.0625*(9.*(CW(1)+CW(3))-E(N-1)-E(N+2))
      CW(4)=0.0625*(9.*(CW(3)+CW(5))-E(N-52)-E(N+107))
      CW(6)=0.0625*(9.*(CW(5)+CW(7))-E(N+55)-E(N+52))
      CW(8)=0.0625*(9.*(CW(7)+CW(9))-E(N+106)-E(N-53))
      DO 17 I=1,9
17    CW(I)=CW(I)-D(M)
      CW(10)=(4.*(CW(2)+CW(4)+CW(6)+CW(8))-CW(1)-CW(3)-CW(5)-CW(7))/12.
      CP=FLOAT(N/53)*200.-99.5
      Q=FLOAT(MOD(N,53))*200.-299.5
      IJ=10
      L=0
18    L=L+1
      K=L+1
      DO 20 NI=1,3
      IF(CW(IJ)*CW(K).GT.0.)GOTO 19
      Z=CW(K)/(CW(K)-CW(IJ))
      T=1.-Z
      NA=MIN0(19997,NA+2)
      IE(NA)=INT(Z*CX(IJ)+T*CX(K)+CP)
      IE(NA+1)=INT(Z*CY(IJ)+T*CY(K)+Q)
      IF(MOD(NA,4).NE.3)GOTO 19
      IF(IE(NA).NE.IE(NA-2))GOTO 19
      IF(IE(NA+1).EQ.IE(NA-1))NA=NA-4
19    I=IJ
      IJ=K
      K=L
20    L=I
      IF(L-8)18,16,16
C
C PLOT CONTOUR
C
21    IF(NA.LT.0)GO TO 30
      N=NA
23    IF(IE(N).LT.1)GOTO 24
      IF(IE(N+1).LT.1)GOTO 24
      IF(IE(N).GT.9999)GOTO 24
      IF(IE(N+1).GT.9999)GOTO 24
      N=N-2
      IF(N.GT.1)GOTO 23
24    NX=IE(N)
      NY=IE(N+1)
      IJ=2
      K=0
      GOTO 28
25    N=-1
      IE(NA+2)=NX
      IE(NA+3)=NY
26    N=N+2
      IF(NX.NE.IE(N))GOTO 26
      IF(NY.NE.IE(N+1))GOTO 26
      IF(N.GT.NA)GOTO 21
      L=8*((N+3)/4)-N-4
      NX=IE(L)
      NY=IE(L+1)
      IE(L)=IE(NA-2)
      IE(L+1)=IE(NA-1)
      IE(N)=IE(NA)
      IE(N+1)=IE(NA+1)
      NA=NA-4
      IF(D(M).GT.DC)K=2-K
      IJ=K
28    CP=0.0001*FLOAT(NX)
      Q=0.0001*FLOAT(NY)
      CP=(CP+0.02)*SCAL
      Q=(Q+0.02)*SCAL
      CALL XPLOD(SCAL,IJ,CP,Q)
      GOTO 25
30    CONTINUE
      DO 261 I=1,IA
      CP=XAT(I)*AY(1)+YAT(I)*AY(2)+ZAT(I)*AY(3)
      Q=XAT(I)*AY(4)+YAT(I)*AY(5)+ZAT(I)*AY(6)
      IF(SP.LT.ABS(XP-CP) .OR. SP.LT.ABS(YP-Q))GO TO 261
      CP=CP+SP-XP+0.025
      Q=Q+SP-YP+0.025
      CALL XPLOD(SCAL,2,CP,Q)
      CP=CP-0.05
      CALL XPLOD(SCAL,0,CP,Q)
      Q=Q-0.05
      CALL XPLOD(SCAL,0,CP,Q)
      CP=CP+0.05
      CALL XPLOD(SCAL,0,CP,Q)
      Q=Q+0.05
      CALL XPLOD(SCAL,0,CP,Q)
      CP=CP-0.025
      Q=Q-0.025
      IZ=NSFAC(I)
      DO 262 IK=1,IA
      IF(I.EQ.IK) GO TO 262
      CALL RAD(XAT(I),YAT(I),ZAT(I),XAT(IK),YAT(IK),ZAT(IK),R)
      K=NSFAC(IK)
      IF(R.GT.CR(IZ)+CR(K)+0.25)GO TO 262
      IJ=2
      CALL XPLOD(SCAL,IJ,CP,Q)
      AP(1)=XAT(IK)*AY(1)+YAT(IK)*AY(2)+ZAT(IK)*AY(3)
      AP(2)=XAT(IK)*AY(4)+YAT(IK)*AY(5)+ZAT(IK)*AY(6)
      AP(1)=AP(1)+SP-XP
      AP(2)=AP(2)+SP-YP
      IJ=0
      CALL XPLOD(SCAL,IJ,AP(1),AP(2))
262   CONTINUE
261   CONTINUE
      IPLOT=0
236   IF(IEND.EQ.0)GO TO 9998
      IEND=0
      GO TO 996
122   FORMAT(//5X,I5,34H  SUITABLE HYDRIDE POSITIONS FOUND)
9998  IF(IPLOT.GT.0 .AND. KL.LE.0)GO TO 248
      IF(KL.LE.0)STOP
      WRITE(NOUT,131)
      WRITE(NOUT,133)
      DO 221 I=1,KL
      WRITE(NOUT,132)I,XLOC(I),YLOC(I),ZLOC(I),FLOC(I)
221   CONTINUE
      IZ=IA+1
      DO 220   I=1,KL
      WRITE(NOUT,114)I
      IF(ILOC.GT.0)IZ=IA-ILOC+I
      DO 222 IJ=1,IA
      IF(IJ.EQ.IZ)GO TO 222
      CALL RAD(XLOC(I),YLOC(I),ZLOC(I),XAT(IJ),YAT(IJ),ZAT(IJ),R)
      IF(R.LT.4.5)WRITE(NOUT,115)NATOM(IJ),R
222   CONTINUE
220   CONTINUE
      IF(ILOC.GT.0)STOP
      IF(KL.LE.1)STOP
      WRITE(NOUT,227)
      DO 224 I=1,KL
      DO 225 IJ=I,KL
      IF(I.EQ.IJ)GO TO 225
      CALL RAD(XLOC(I),YLOC(I),ZLOC(I),XLOC(IJ),YLOC(IJ),ZLOC(IJ),R)
      WRITE(NOUT,226)I,IJ,R
225   CONTINUE
224   CONTINUE
226   FORMAT(5X,I4,8H .......,I4,5X,F7.4)
227   FORMAT(//5X,25H INTER-HYDRIDE DISTANCES,//7X,10H HOLE NOS.//)
      KL=0
      IF(IPLOT.GT.0)GO TO 248
999   STOP
      END
C
C SUBROUTINE USING CALCOMP-TYPE PLOTTING COMMAND "PLOT"
      SUBROUTINE XPLOD(SCAL,N,X,Y)
      P=25.4*(AMAX1(0.,AMIN1(X,SCAL)))
      Q=25.4*(AMAX1(0.,AMIN1(Y,SCAL)))
      CALL PLOT(P,Q,N)
      RETURN
      END
C
      SUBROUTINE RAD(X1,Y1,Z1,X2,Y2,Z2,R)
      INTEGER ED(10)
      COMMON EA(10),EB(10),EC(10),X ,Y ,Z ,COSA,COSB,COSC,HRAD(3),AX(3),
     *HSD(3),A1,A2,A3,XAT(150),YAT(150),ZAT(150),IR,ISRCH,ED,ILIST(150),
     *NSFAC(150)
      DX=X1-X2
      DY=Y1-Y2
      DZ=Z1-Z2
      R=DX*(A1*DX+DY*Z+DZ*Y)+DY*(A2*DY+DZ*X)+DZ*A3*DZ
      R=SQRT(R)
      RETURN
      END
C
C
      SUBROUTINE FUNCT(N,NF,XC,FC,GC)
      INTEGER ED(10)
      DIMENSION XC(3),GC(3)
      COMMON EA(10),EB(10),EC(10),X ,Y ,Z ,COSA,COSB,COSC,HRAD(3),AX(3),
     *HSD(3),A1,A2,A3,XAT(150),YAT(150),ZAT(150),IR,ISRCH,ED,ILIST(150),
     *NSFAC(150)
      FC=0.E0
      DO 1 I=1,3
1     GC(I)=0.E0
      DO 2 I=1,IR
      IZ=ILIST(I)
      IS=NSFAC(IZ)
      CALL RAD(XC(1),XC(2),XC(3),XAT(IZ),YAT(IZ),ZAT(IZ),R)
      IF(NF.EQ.3)GO TO 4
      IF(R.GT.0.5)GO TO 6
      NF=-NF
      RETURN
6     DX=XC(1)-XAT(IZ)
      DY=XC(2)-YAT(IZ)
      DZ=XC(3)-ZAT(IZ)
      PX=((DY*Z+DZ*Y)/2.E0+DX*A1)/R
      PY=(DY*A2+(DX*Z+DZ*X)/2.E0)/R
      PZ=(DZ*A3+(DX*Y+DY*X)/2.E0)/R
4     IF(I.LE.ISRCH)GO TO 3
      IF(R.GT.6.) GO TO 2
      IF(NF.NE.3 .OR. R.GT.0.9)GO TO 5
      R=0.9
5     Z1=EA(IS)*EXP(-EB(IS)*R)
      Z1=Z1/R**ED(IS)
      Z2=-EC(IS)/R**6
      FC=Z1+Z2+FC
      IF(NF.EQ.3)GO TO 2
      Z3=-Z1*(EB(IS)+FLOAT(ED(IS))/R)-6.*Z2/R
      GC(1)=Z3*PX+GC(1)
      GC(2)=Z3*PY+GC(2)
      GC(3)=Z3*PZ+GC(3)
      GO TO 2
3     IF(NF.EQ.2)GO TO 2
      H=(R-HRAD(I))
      FC=FC+H*H/(HSD(I)*HSD(I))
      IF(NF.EQ.3)GO TO 2
      GC(1)=GC(1)+PX*H*2.E0/(HSD(I)*HSD(I))
      GC(2)=GC(2)+PY*H*2.E0/(HSD(I)*HSD(I))
      GC(3)=GC(3)+PZ*H*2.E0/(HSD(I)*HSD(I))
2     CONTINUE
      NF=1
      RETURN
      END
C
      SUBROUTINE DECODE(A,N,JC)
      DIMENSION JC(77),A(36),M(13),L(3)
      DATA M/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H.,1H-,1H /
      L(1)=0
      L(2)=0
      L(3)=0
      S=1.E0
      N=1
      JC(77)=M(13)
      DO 1 I=1,36
    1 A(I)=0.E0
      DO 11 I=1,77
      DO 2 J=1,12
      IF(JC(I).EQ.M(J))GOTO 7
    2 CONTINUE
      IF(L(2))3,11,3
    3 IF(L(3)-1)4,5,4
    4 A(N)=AINT(A(N)*10.E0**L(1)+0.1E-1)
    5 A(N)=A(N)*S
      N=N+1
      S=1.E0
      DO 6 K=1,3
    6 L(K)=0
      GOTO 11
    7 L(2)=1
      IF(J-11)10,9,8
    8 S=-1.E0
      GOTO 11
    9 A(N)=A(N)*10.E0**L(1)
      L(1)=0
      L(3)=1
      GOTO 11
   10 L(1)=L(1)+1
      A(N)=A(N)+FLOAT(J-1)/10.E0**L(1)
   11 CONTINUE
      N=N-1
      RETURN
      END
C
      SUBROUTINE MIN(N,X,F,G,XTOL,STEPMX,FEST,MAXCAL,IFAIL)
      DIMENSION OLDX(3),STEPMX(3),XTOL(3),X(3),G(3),P(3)
      NOUT=6
      NJ=N+1
      NCALL=0
      NF=1
      SCAL=1.E0
      CALL FUNCT(N,NF,X,F,G)
      IF(NF.LT.0)GO TO 23
      NCALL=NCALL+1
      ITER=0
      GG=0.E0
5     OLDF=F
      DO 20 I=1,N
20    OLDX(I)=X(I)
      DO 3 J=1,NJ
      ITER=ITER+1
      DO 17 I=1,N
17    GG=GG+G(I)*G(I)
      IF(GG.GT.0.1E-6) GO TO 21
      GO TO 31
21    IF(J.GT.1) GO TO 15
      DO 1 I=1,N
1      P(I)=-G(I)
      GO TO 14
15    BETA=GG/OLDGG
      DO 16 I=1,N
16    P(I)=-G(I)+BETA*P(I)
14    YB=F
      PP=0.E0
      VB=0.E0
      DO 2 I=1,N
      PP=P(I)*P(I)+PP
2     VB=VB+P(I)*G(I)
      IF(VB.GE.0.E0) GO TO 4
      AK=2.E0*(FEST-F)/VB
      H=1.E0/SQRT(PP)
      IF(AK.LT.0.E0 .OR. AK.GT.H)AK=H
      H=AK
      DO 32 I=1,N
      STEP=ABS(H*P(I)*SCAL)
      IF(STEP.GT.STEPMX(I))SCAL=SCAL*STEPMX(I)/STEP
32    CONTINUE
      H=H*SCAL
9     YA=YB
      VA=VB
      DO 6 I=1,N
6     X(I)=X(I)+H*P(I)
      CALL FUNCT(N,NF,X,F,G)
      IF(NF.LT.0)GO TO 23
      NCALL=NCALL+1
      IF(NCALL.GT.MAXCAL) GO TO 23
      VB=0.E0
      YB=F
      DO 7 I=1,N
7     VB=VB+G(I)*P(I)
      IF(VB.GE.0.E0 .OR. YB.GE.YA) GO TO 8
      GO TO 9
8     T=0.E0
13    IF(SQRT(ABS(VB))*SQRT(ABS(VA)).GE.1.E16)GO TO 24
      Z=3.E0*(YA-YB)/H+VA+VB
      W=Z*Z-VA*VB
      IF(W.LE.0.E0)GO TO 4
      W=SQRT(W)
      AK=H*(VB+W-Z)
      Z=VB-VA+2.E0*W
      AK=AK/Z
25    DO 10 I=1,N
10    X(I)=X(I)+(T-AK)*P(I)
      CALL FUNCT(N,NF,X,F,G)
      IF(NF.LT.0)GO TO 23
      NCALL=NCALL+1
      IF(NCALL.GT.MAXCAL)GO TO 23
      IF(F.LE.YA .AND. F.LE.YB) GO TO 4
      DO 131 I=1,N
      Z=ABS((T-AK)*P(I))
      IF(Z.GT.0.1E0*XTOL(I)) GO TO 132
131   CONTINUE
      GO TO 4
24    AK=H/2.
      GO TO 25
132   VC=0.E0
      DO 11 I=1,N
11    VC=VC+G(I)*P(I)
      IF(VC.GE.0.E0 .OR. F.GE.YA ) GO TO 12
      YA=F
      VA=VC
      H=AK
      GO TO 13
12    YB=F
      VB=VC
      H=H-AK
      T=0.E0
      GO TO 13
4     OLDGG=GG
      GG=0.E0
3     CONTINUE
      DO 19 I=1,N
      Z=OLDX(I)-X(I)
      IF(ABS(Z).GE.XTOL(I))GO TO 18
19    CONTINUE
      GO TO 31
18    IF(F.LT.OLDF) GO TO 5
      DO 22 I=1,N
22    X(I)=OLDX(I)
      IFAIL=1
      RETURN
31    IFAIL=0
      RETURN
23    IFAIL=1
      RETURN
      END
      SUBROUTINE INV(A,AI,DET)
C
C   SUBROUTINE FOR THE INVERSION OF A 3*3 MATRIX
C      INPUT:  A(1) A(2) A(3)
C              A(4) A(5) A(6)
C              A(7) A(8) A(9)
C      OUTPUT:  AI(1) AI(2) AI(3)
C               AI(4) AI(5) AI(6)
C               AI(7) AI(8) AI(9)
C
      DIMENSION A(9), AI(9)
      BC1 = A(5) * A(9) - A(8) * A(6)
      BC2 = A(6) * A(7) - A(4) * A(9)
      BC3 = A(4) * A(8) - A(5) * A(7)
      DET = A(1) * BC1 + A(2) * BC2 + A(3) * BC3
      IF (ABS(DET) .LE. 1.E-6) GO TO 20
      DETI = 1.0 / DET
      AI(1) = BC1 * DETI
      AI(4) = BC2 * DETI
      AI(7) = BC3 * DETI
      AI(2) = (A(8) * A(3) - A(2) * A(9)) * DETI
      AI(5) = (A(1) * A(9) - A(3) * A(7)) * DETI
      AI(8) = (A(7) * A(2) - A(1) * A(8)) * DETI
      AI(3) = (A(2) * A(6) - A(5) * A(3)) * DETI
      AI(6) = (A(4) * A(3) - A(1) * A(6)) * DETI
      AI(9) = (A(1) * A(5) - A(2) * A(4)) * DETI
      GO TO 25
  20  DET = 0.0
  25  RETURN
      END
