      PROGRAM HYDROGEN
C**********************************************************************
C* January 1999                                                       *
C*                     Program HYDROGEN                               *
C*                                                                    *
C* Modelling Hydroxyl and Water Hydrogens: How to localize hydroxyl   *
C* and water hydrogens when diffraction data are not good enough.     *
C*                                                                    *
C*                      by M. Nardelli                                *
C* Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica  *
C* e Chimica Fisica della Universita' di Parma, Centro di Studio per  *
C* la Strutturistica Diffrattometrica del CNR.                        *
C*           Viale delle Scienze 78, I-43100 Parma, Italy.            *
C*           E-mail: nardelli@ipruniv.cce.unipr.it                    *
C*                                                                    *
C* This program calculates the hydroxyl and water H coordinates on    *
C* the basis of hydrogen bonding geometry and energy considerations.  *
C*                                                                    *
C*                        INPUT                                       *
Cp  - Input of PARST97 without the hydrogens to be calculated.        *
C*  - File containing the atom partial charges with the following     *
C*    structure:                                                      *
C*       Title (80A1)                                                 *
C*       N,Y0,HD (Free format) N  = number of hydroxyls or waters to  *
C*                                  consider (max. 50)                *
C*                             Y0 = coefficient for Gaussian factor   *
C*                                  suggested value: 5.0              *
C*                             HD = O-H distance, e.g. 0.938 A for    *
C*                                  neutron, 0.850 A for X ray        *
C*       FRT              format of data to be read e.g. (6A1,F7.3)   *
C*       OOOOOO  x.xxx    (OOOOOO = label of the hydroxyl or water    *
C*       OOOOOO  x.xxx    oxygen with the atomic partial charge,      *
C*       OOOOOO  x.xxx    x.xxxx, of the attached hydrogen,           *
C*       ......  .....    N of these entries)                         *
C*       AAAAAA  y.yyy    (AAAAAA = atom label with its atomic        *
C*       AAAAAA  y.yyy     partial charge. These atom labels must be  *
C*       ......  .....     equal to those of the PARST97 input)       *
C*                                                                    *
C*                            OUTPUT                                  *
C*     Hydroxyl or water hydrogen coordinates and hydrogen bonds      *
C*     On request: geometric data and non bonded energy profiles      *
C*                                                                    *
C*  N.B. - The program ask for the problem: hydroxyl or water, and    *
C*         for the names of the input and output files.               *
C*                                                                    *
C**********************************************************************
      CHARACTER*1 BRV,FRT(80),STPE(38,20),LTP(7),MINUS(2),UNIM(3,5),
     1 II(5),LA1(2,20),AT(500,6),AA(6),O1L(6),ATO(50,6),AAAA(6),QS,
     2 SDL1(2,500),LBA(2),ATH(6),ATHO(50,6),ACT(6),ACCT(50,6),STR(20),
     3 ATR1(3,6),ATW(3,6),ATA(2,6),ATAA(50,2,6),ATHW(50,2,6)
      CHARACTER*2 LA(20),SDL2(500)
      DIMENSION SC(3),SA(3),S(3),U(6),SU(6),QH(50),IATO(50),XS(50,3),
     1 XHX(50,3),ISO(50),XHH(3),XOCT(3),XCCT(50,3),IA(2,200),BB(3,3),
     2 FB(3),XHN(2,3),XON(2,3),XOAA(50,2,3),XOA(2,3),
     3 XHW(50,2,3),XHOW(50,2,3),XA(2,3),XAA(50,2,3),LE1(24,3,3),
     4 F1(24,3)
      EQUIVALENCE (LA(1),LA1(1,1)),(SDL2(1),SDL1(1,1))
      COMMON/GRP1/C(3),A(3)
     1      /GRP3/XO(500,3)
     2      /GRP4/AT
     3      /GRP5/MC(400,8)
     4      /GRP10/LA
     5      /GRP11/ANAA(20)
     6      /GRP12/NNS
     7      /CDO/KHW,KKD1,KKD2
     8      /LTT/KBB,BB
      COMMON/COOR/X(500,3),Q(500)
     1      /NEW1/O(3,3),OI(3,3)
     2      /INS4/CS(3),SN(3)
     3      /LBL4/STPE
     4      /LBL3/IC,NE,NT
     5      /INPOU/IN,IO,IR
     6      /GRG1/SDL2
     7      /EQVP/LE1,F1,FB
      DATA LTP /'P','A','B','C','I','F','R'/
      RD=57.295779513
      IN=13
      IO=15
      IR=16
      CALL INPOUT
C-----Questions concerning data and printings
      WRITE(*,'('' If it is a hydroxyl problem key in  0''/
     1 '' If it is a water problem key in     1''/'' ?>'')')
      READ(*,*)KHW
      WRITE(*,'('' Do you want geometric data and non bonded '',
     1 ''energy profiles?''/'' Y/N:>'')')
      READ(*,'(A1)')QS
      KKD1=1
      IF(QS.EQ.'N'.OR.QS.EQ.'n')KKD1=0
C     WRITE(*,'('' Do you want the data of profiles suitable '',
C    1 ''for graphic software?''/'' Y/N:>'')')
C     READ(*,'(A1)')QS
C     KKD2=1
C     IF(QS.EQ.'N'.OR.QS.EQ.'n')KKD2=0
      KKD2=0
C-----Read title, codes and parameters
      READ(IN,'(80A1)') FRT
      WRITE(IO,'(80A1)') FRT
      READ(IN,'(80A1)') FRT
      BRV=FRT(1)
      DO 1 I=1,7
      IF(BRV.EQ.LTP(I))GO TO 2
1     CONTINUE
      WRITE(IO,100)BRV
100   FORMAT(/' The symbol ',A1,'for lattice type is uncorrect.'/
     1 ' Symbol P has been assumed'/)
      BRV=LTP(1)
2     READ(IN,*)N,D3,DM,NH1,NTN,LSP,LSL,NPR,NST,KL,KY,KC,KO,KB,KA,KT,
     1 KD,KE,KOR,NI
      IF(N.LE.500)GO TO 3
      N=500
      WRITE(IO,101)
101   FORMAT(/'There are too many atoms.'/
     1 ' The calculations are carried out on the first 500 atoms.'/)
3     NCS=2
      IF(KY.EQ.0) GO TO 5
      READ(IN,*)NNS,Z,AL
      READ(IN,*)(LA(I),ANAA(I),I=1,NNS)
      DO 4 I=1,NNS
      MINUS(1)=LA1(1,I)
      MINUS(2)=LA1(2,I)
      CALL MINUMA(NCS,MINUS)
      LA1(1,I)=MINUS(1)
4     LA1(2,I)=MINUS(2)
5     READ(IN,*)C,SC,A,SA
C-----Read atom symbols,coordinates,thermal parameters & s.u.'s
      IF(KL.EQ.4) GO TO 17
      READ(IN,'(80A1)') FRT
      IF(KOR.EQ.0) GO TO 6
      GO TO (7,9,11),KL
6     GO TO (13,13,15,17),KL
7     DO 8 J=1,N
8     READ(IN,FRT)(AT(J,L),L=1,6),(X(J,I),I=1,3),(S(I),I=1,3),
     1 (U(I),I=1,6),(SU(I),I=1,6)
      GO TO 20
9     DO 10 J=1,N
10    READ(IN,FRT)(AT(J,L),L=1,6),(X(J,I),I=1,3),(U(I),I=1,6),
     1 (S(I),I=1,3),(SU(I),I=1,6)
      GO TO 20
11    DO 12 J=1,N
12    READ(IN,FRT)(AT(J,L),L=1,6),(X(J,I),S(I),I=1,3),
     1 (U(I),SU(I),I=1,6)
      GO TO 20
13    DO 14 J=1,N
14    READ(IN,FRT)(AT(J,L),L=1,6),(X(J,I),I=1,3),(S(I),I=1,3)
      GO TO 20
15    DO 16 J=1,N
16    READ(IN,FRT)(AT(J,L),L=1,6),(X(J,I),S(I),I=1,3)
      GO TO 20
17    DO 19 J=1,N
      READ(IN,'(7X,6A1,2X,3F8.5,12X,3(5A1,3X))')(AT(J,L),L=1,6),
     1 (X(J,I),I=1,3),((UNIM(I,K),K=1,5),I=1,3)
      DO 19 I=1,3
      DO 18 JJ=1,5
18    II(JJ)=UNIM(I,JJ)
      CALL CONVE(II,DH)
19    S(I)=DH
20    CONTINUE
      NSC=2
      DO 22 I=1,N
      DO 21 K=1,5,2
      MINUS(1)=AT(I,K)
      MINUS(2)=AT(I,K+1)
      CALL MINUMA(NSC,MINUS)
      AT(I,K)= MINUS(1)
21    AT(I,K+1)= MINUS(2)
22    CONTINUE
      DO 24 I=1,N
      DO 23 L=1,6
23    AAAA(L)=AT(I,L)
      CALL LABL(AAAA,KY,NNS,LA1,LBA,KKD)
      SDL1(1,I)=LBA(1)
      SDL1(2,I)=LBA(2)
24    CONTINUE
C-----Omit unecessary lines
      CALL DMLN(NH1,NTN,LSP,LSL,NPR,NST)
C-----Read equivalent positions
      READ(IN,*)IC,NE,NT
      DO 25 I=1,NE
      READ(IN,'(20A1)')STR
      CALL MINUMA(20,STR)
      DO 25 L=1,20
25    STPE(I,L)=STR(L)
      NE=NE+1
C-----Orthogonalization of coordinates
      CALL ORTHO(N,RD)
      CALL INVMAT3(O,OI,ICD)
C-----Sets up the connectivity matrix
      CALL BONDI(N)
C-----Read OH oxygens, Y0, O-H distance and H charges
      READ(IR,'(80A1)')FRT
      READ(IR,*)NO,Y0,HD
      IF(HD.LT.0.5.OR.HD.GT.1.5)THEN
      WRITE(IO,'(''There is something wrong with the O-H bond '',
     1 ''distance'')')
      STOP
      ENDIF
      IF(NO.GT.50)THEN
      NO=50
      WRITE(IO,'(''The number of hydroxyl groups or water molecules''/
     1 ''is greater than 50''/''Calculations will be carried on the '',
     2 ''first 50s''/)')
      ENDIF
      READ(IR,'(80A1)')FRT
      DO 26 I=1,NO
26    READ(IR,FRT)(ATO(I,L),L=1,6),QH(I)
      QH(I)=QH(I)
C-----Read atom charges
      NN=0
27    NN=NN+1
      READ(IR,FRT,ERR=30,END=30)(AA(L),L=1,6),QQ
      MINUS(1)=AA(1)
      MINUS(2)=AA(2)
      CALL MINUMA(2,MINUS)
      AA(1)=MINUS(1)
      AA(2)=MINUS(2)
      DO 29 I=1,N
      DO 28 L=1,6
      IF(AA(L).NE.AT(I,L))GO TO 29
28    CONTINUE
      Q(I)=QQ
      GO TO 90
29    CONTINUE
      WRITE(IO,'(''There is no correspondence between the atom '',
     1 ''labels and the two input files'')')
      STOP
90    GO TO 27
30    NN=NN-1
      IF(NN.NE.N)THEN
      WRITE(IO,'(''The number of the atom charges is not equal to '',
     1 ''the number of the atoms''/''Check your data'')')
      STOP
      ENDIF
      DO 33 J=1,NO
      DO 32 I=1,N
      DO 31,L=1,6
      IF(ATO(J,L).NE.AT(I,L))GO TO 32
31    CONTINUE
      IATO(J)=I
      GO TO 33
32    CONTINUE
33    CONTINUE
C-----Determination of the starting H coordinates for hydroxyl
      IF(KHW.EQ.1)GO TO 52
      DO 37 I=1,NO
      JO1=IATO(I)
      DO 34 L=1,6
34    O1L(L)=AT(IATO(I),L)
      IF(KKD1.EQ.1.AND.KHW.EQ.0)WRITE(IO,'(/19(''*'')/''* HYDROXYL '',
     1 6A1,'' *''/19(''*'')/)')O1L
      KKR=0
      CALL STARTH(KY,BRV,RD,N,Y0,HD,JO1,O1L,ATH,XHH,ISS,ACT,XOCT,KKR)
      IF(KKR.EQ.1)GO TO 37
      DO 35 L=1,6
      ATHO(I,L)=ATH(L)
35    ACCT(I,L)=ACT(L)
      DO 36 K=1,3
      XHX(I,K)=XHH(K)
36    XCCT(I,K)=XOCT(K)
37    ISO(I)=ISS
C-----Adds the hydroxyl hydrogens to the atom list and connectivity
C-----matrix
      DO 41 I=1,NO
      N=N+1
      DO 38 L=1,6
38    AT(N,L)=ATHO(I,L)
      DO 39 J=1,3
      DO 39 K=1,3
      XS(I,K)=XHX(I,K)
      X(N,K)=XHX(I,K)
39    XO(N,J)=XO(N,J)+O(J,K)*XHX(I,K)
      DO 88 L=1,8
      IF(MC(IATO(I),L).EQ.0)THEN
      MC(IATO(I),L)=N
      GO TO 89
      ENDIF
88    CONTINUE
89    MC(N,1)=IATO(I)
      DO 40 L=2,8
40    MC(N,L)=0
41    Q(N)=QH(I)
C-----Matrix of the atom symbols
      DO 43 I=1,N
      DO 42 L=1,6
42    AAAA(L)=AT(I,L)
      CALL LABL(AAAA,KY,NNS,LA1,LBA,KKD)
      SDL1(1,I)=LBA(1)
      SDL1(2,I)=LBA(2)
43    CONTINUE
C-----Hydrogen bonding with the starting H coordinates
      IF(KKD1.EQ.1)WRITE(IO,'(''Hydrogen bonding parameters with '',
     1 ''the starting H coordinates'')')
      DO 46 I=1,NO
      JO1=IATO(I)
      ISS=ISO(I)
      NN=N-NO+I
      DO 44 L=1,6
44    ACT(L)=ACCT(I,L)
      DO 45 K=1,3
45    XOCT(K)=XCCT(I,K)
46    CALL HBPAR(RD,JO1,NN,ACT,XOCT)
C-----Calculates the non bonded energy profiles
      DO 49 IK=1,3
      DO 49 I=1,NO
      ISS=ISO(I)
      DO 47 L=1,6
      ATR1(1,L)=AT(ISS,L)
      ATR1(2,L)=ATO(I,L)
      ATR1(3,L)=ATHO(I,L)
47    ACT(L)=ACCT(I,L)
      DO 48 K=1,3
48    XOCT(K)=XCCT(I,K)
      JO1=IATO(I)
      NNH=N-NO+I
      CALL ENRPRF(KY,BRV,RD,N,Y0,JO1,ISS,NNH,ATR1,ACT,XOCT)
49    CONTINUE
      WRITE(IO,'(/''Coordinates from geometric environments'')')
      DO 50 I=1,NO
      J=N-NO+I
50    WRITE(IO,'(6A1,3X,3F9.5)')(AT(J,L),L=1,6),(XS(I,K),K=1,3)
      WRITE(IO,'(/''Coordinates after profile refinement'')')
      DO 51 I=1,NO
      J=N-NO+I
51    WRITE(IO,'(6A1,3X,3F9.5)')(AT(J,L),L=1,6),(X(J,K),K=1,3)
      GO TO 80
C-----Determination of the starting H coordinates for water
52    DO 58 I=1,NO
      JO1=IATO(I)
      DO 53 L=1,6
53    O1L(L)=AT(IATO(I),L)
      IF(KKD1.EQ.1.AND.KHW.EQ.1)WRITE(IO,'(/19(''*'')/''* WATER    '',
     1 6A1,'' *''/19(''*'')/)')O1L
      KKK=0
      CALL STARTW(KY,BRV,RD,N,Y0,HD,JO1,O1L,ATA,XA,XOA,ATW,XHN,XON,KKK)
      IF(KKK.EQ.1)THEN
      WRITE(IO,'(''It is not possible to find suitable directions '',
     1 ''for hydrogens'')')
      STOP
      ENDIF
      DO 54 L=1,6
      ATHW(I,1,L)=ATW(1,L)
54    ATHW(I,2,L)=ATW(3,L)
      DO 56 J=1,2
      DO 55 K=1,3
      XHW(I,J,K)=XHN(J,K)
      XHOW(I,J,K)=XON(J,K)
      XAA(I,J,K)=XA(J,K)
55    XOAA(I,J,K)=XOA(J,K)
56    CONTINUE
      DO 57 L=1,6
      ATAA(I,1,L)=ATA(1,L)
57    ATAA(I,2,L)=ATA(2,L)
58    CONTINUE
C-----Add water hydrogens to the atom list and connectivity matrix
      DO 62 I=1,NO
      DO 62 J=1,2
      N=N+1
      Q(N)=QH(I)
      DO 59 L=1,6
59    AT(N,L)=ATHW(I,J,L)
      DO 60 K=1,3
      X(N,K)=XHW(I,J,K)
60    XO(N,K)=XHOW(I,J,K)
      DO 86 L=1,8
      IF(MC(IATO(I),L).EQ.0)THEN
      MC(IATO(I),L)=N
      GO TO 87
      ENDIF
86    CONTINUE
87    MC(N,1)=IATO(I)
      DO 61 L=2,8
61    MC(N,L)=0
62    CONTINUE
C-----Matrix of the atom symbols
      DO 64 I=1,N
      DO 63 L=1,6
63    AAAA(L)=AT(I,L)
      CALL LABL(AAAA,KY,NNS,LA1,LBA,KKD)
      SDL1(1,I)=LBA(1)
64    SDL1(2,I)=LBA(2)
C-----Hydrogen bonding with the starting bonding coordinates
      IF(KKD1.EQ.1)WRITE(IO,'(/''Hydrogen bonding parameters with '',
     1 ''the starting H coordinates'')')
      NN=N-2*NO
      DO 70 I=1,NO
      JO1=IATO(I)
      NN=NN+1
      DO 70 J=1,2
      NN=NN-1+J
      DO 65 L=1,6
65    ACT(L)=ATAA(I,J,L)
      DO 69 K=1,3
69    XOCT(K)=XOAA(I,J,K)
70    CALL HBPAR(RD,JO1,NN,ACT,XOCT)
C-----Calculates the non bonded energy profiles
      DO 77 IK=1,3
      NN=N-2*NO
      DO 77 I=1,NO
      NN=NN+1
      ISS=NN+1
      JO1=IATO(I)
      DO 71 L=1,6
      ATR1(2,L)=ATO(I,L)
71    ACT(L)=ATAA(I,1,L)
      DO 72 K=1,3
72    XOCT(K)=XOAA(I,1,K)
      DO 73 L=1,6
      ATR1(1,L)=ATHW(I,2,L)
73    ATR1(3,L)=ATHW(I,1,L)
      CALL ENRPRF(KY,BRV,RD,N,Y0,JO1,ISS,NN,ATR1,ACT,XOCT)
      NN=NN+1
      ISS=NN-1
      DO 74 L=1,6
74    ACT(L)=ATAA(I,2,L)
      DO 75 K=1,3
75    XOCT(K)=XOAA(I,2,K)
      DO 76 L=1,6
      ATR1(1,L)=ATHW(I,1,L)
76    ATR1(3,L)=ATHW(I,2,L)
      CALL ENRPRF(KY,BRV,RD,N,Y0,JO1,ISS,NN,ATR1,ACT,XOCT)
77    CONTINUE
      WRITE(IO,'(/''Coordinates from geometric environments'')')
      DO 78 I=1,NO
78    WRITE(IO,'(6A1,3X,3F9.5)')((ATHW(I,J,L),L=1,6),(XHW(I,J,K),
     1 K=1,3),J=1,2)
      WRITE(IO,'(/''Coordinates after profile refinement'')')
      NN=N-2*NO+1
      DO 79 I=NN,N
79    WRITE(IO,'(6A1,3X,3F9.5)')(AT(I,L),L=1,6),(X(I,K),K=1,3)
80    NH2=0
      DO 81 I=1,N
      IF(SDL2(I).EQ.'H ')THEN
      NH2=NH2+1
      IA(1,NH2)=I
      IA(2,NH2)=MC(I,1)
      ENDIF
81    CONTINUE
      WRITE(IO,102)
102   FORMAT(/'Hydrogen Bonds'//'D-H...A',17X,'D-H',5X,'D...A',4X,
     1 'H...A',3X,'D-H...A',2X,'H-D...A',1X,'H...A...D'/)
      DO 82 L=1,3
82    FB(L)=0.0
      CALL HDRB(N,SDL2,NH2,IA)
      IF(KBB.EQ.0)GO TO 85
      DO 84 L=1,KBB
      DO 83 K=1,3
83    FB(K)=BB(L,K)
      CALL HDRB(N,SDL2,NH2,IA)
84    CONTINUE
85    CLOSE (IN)
      CLOSE (IO)
      CLOSE (IR)
      STOP
      END
 
      SUBROUTINE ADD(N,NL,AAAA,XO3,XXO,Q3,QF,NG3)
      CHARACTER*1 AT(500,6),AAAA(500,6)
      DIMENSION XO3(500,3),XXO(500,3),Q3(500),QF(500)
      COMMON/GRP4/AT
     1      /GRP6/NA1(500),NA3(500)
     2      /INPOU/IN,IO,IR
      DO 2 I=1,NL
      DO 1 K=1,3
1     XO3(I+NG3,K)=XXO(I,K)
      Q3(I+NG3)=QF(I)
      NA3(I+NG3)=I+N
      DO 2 L=1,6
2     AT(I+N,L)=AAAA(I,L)
      NG3=NG3+NL
      N=N+NL
      RETURN
      END
 
      SUBROUTINE BONDI(N)
C-----This routine calculates bond distances and set up the the
C-----connectivity matrix MC
      CHARACTER*1 AT(500,6)
      CHARACTER*2 NOMAT(103),KS(2),SDL2(500)
      DIMENSION RG(2),IRAG(103)
      COMMON/GRP4/AT
     2      /GRP5/MC(400,8)
     3      /GRP3/XO(500,3)
     4      /SMBL/NOMAT
     5      /ATRD/IRAG
     6      /GRG1/SDL2
     7      /INPOU/IN,IO,IR
      DO 1 J=1,N
      KK=0
      DO 1 K=1,N
      D1=0.0
      DO 2 I=1,3
      DR=XO(J,I)-XO(K,I)
      IF(ABS(DR).LT.1E-9)THEN
      DR=0.0
      ELSE
      DR=DR**2
      ENDIF
2     D1=D1+DR
      IF(D1.GT.1E-2)D1=SQRT(D1)
      IF(D1.LT.1E-3) GO TO 1
      KS(1)=SDL2(J)
      KS(2)=SDL2(K)
      IF (KS(1).EQ.NOMAT(1).AND.KS(2).EQ.NOMAT(1))GO TO 1
      DO 6 JJ=1,2
      DO 4 I=1,103
      IF(KS(JJ).EQ.NOMAT(I)) GO TO 5
4     CONTINUE
      WRITE(*,'(/'' There is something wrong with the atom symbols''/
     1 '' Check the atom labels''/)')
      GO TO 1
5     RG(JJ)=IRAG(I)
6     CONTINUE
      DSTC=(RG(1)+RG(2))/100.
      IF(D1.GT.(DSTC+0.5))GO TO 1
      KK=KK+1
      IF (KK.GT.8) GO TO 1
      MC(J,KK)=K
1     CONTINUE
      RETURN
      END
 
      SUBROUTINE BNDAN(RD,X,THT)
C-----This routine calculates a bond angle
      DIMENSION X(3,3),D(3,3),DT(3)
      DO 1 I=1,3
      D(1,I)=X(1,I)-X(2,I)
      D(2,I)=X(3,I)-X(2,I)
1     D(3,I)=X(3,I)-X(1,I)
      DO 3 I=1,3
      DT(I)=0.
      DO 2 J=1,3
2     DT(I)=DT(I)+D(I,J)**2
3     IF(DT(I).GT.1E-2)DT(I)=SQRT(DT(I))
      DO 4 I=1,3
      DO 4 J=1,3
4     D(I,J)=D(I,J)/DT(I)
      CST=0.
      DO 5 K=1,3
5     CST=CST+(D(1,K)*D(2,K))
      IF(ABS(CST).GT.1.) THEN
      RS=1.0
      CST=SIGN(RS,CST)
      ENDIF
      SNT=SQRT(ABS(1.-CST*CST))
      THT=RD*ATAN2(SNT,CST)
      RETURN
      END
 
      SUBROUTINE CONTACT(BRV,N,NR,AFR,NL,AAAA,XXO,QF)
C-----Calculates the coordinates of the environment of the rotating
C-----atoms
      CHARACTER*1 BRV,STPE(38,20),AAAA(500,6),AFR(3,6),
     1 AT(500,6)
      DIMENSION XXO(500,3),QF(500),LE1(24,3,3),F1(24,3),FB(3)
      COMMON/GRP1/C(3),A(3)
     1      /LBL3/IC,NE,NT
     2      /LBL4/STPE
     3      /LBL8/XX(500,3)
     4      /NEW1/O(3,3),OI(3,3)
     5      /GRP4/AT
     6      /GRP3/XO(500,3)
     7      /COOR/X(500,3),Q(500)
     8      /INS4/CS(3),SN(3)
     9      /INPOU/IN,IO,IR
     1      /EQVP/LE1,F1,FB
      RD=57.295779513
      D3=2.0
      DM=4.0
C-----Intermolecular contacts less than DM=4.0 A
      NL=0
      CALL INTERC(N,DM,NL,BRV,NR,AFR,AAAA)
C-----File of coordinates for rotener
      CALL ELIMRA(NL,AAAA)
      DO 7 I=1,NL
      DO 6 J=1,N
      DO 4 L=1,6
      IF(AAAA(I,L).NE.AT(J,L))GO TO 6
4     CONTINUE
      DO 5 JJ=1,3
      XXO(I,JJ)=0.0
      DO 5 K=1,3
5     XXO(I,JJ)=XXO(I,JJ)+O(JJ,K)*XX(I,K)
      QF(I)=Q(J)
      GO TO 7
6     CONTINUE
      WRITE(IO,'(''There is some error in the atom labels'')')
      GO TO 8
7     CONTINUE
8     RETURN
      END
 
      SUBROUTINE CONVE(II,DX)
      DIMENSION IX(5)
      CHARACTER*1 II(5),IA(10)
      DATA IA/'0','1','2','3','4','5','6','7','8','9'/
      DO 2 I=1,5
      IX(I)=0
      DO 3 K=1,10
      IF(IA(K).EQ.II(I)) GO TO 1
3     CONTINUE
1     IX(I)=K-1
2     CONTINUE
      IF(II(5).NE.' ') DX=IX(1)*0.1+IX(2)*0.01+IX(3)*0.001+
     1 IX(4)*0.0001+IX(5)*0.00001
      IF(II(5).EQ.' ') DX=IX(1)*0.01+IX(2)*0.001+IX(3)*0.0001+
     1 IX(4)*0.00001
      IF(II(4).EQ.' ') DX=IX(1)*0.001+IX(2)*0.0001+IX(3)*0.00001
      IF(II(3).EQ.' ') DX=IX(1)*0.0001+IX(2)*0.00001
      IF(II(2).EQ.' ') DX=IX(1)*0.00001
      IF(II(1).EQ.' ') DX=0.0
      RETURN
      END
 
      SUBROUTINE COORDH(RD,ATM1,X,DT,W,T,XHH)
C-----Calculate the hydroxyl-H coordinate
      CHARACTER*1 ATM1(4,6)
      DIMENSION R(3,3),RI(3,3),X(4,3),XO(4,3),XO1(4,3),XO2(4,3),
     1 X3(4,3),XO3(4,3),XO4(3,3),DX(3),DT(4),W(4),T(4),XHH(3)
      COMMON/GRP1/C(3),A(3)
     1      /NEW1/O(3,3),OI(3,3)
     2      /INPOU/IN,IO,IR
     3      /CDO/KHW,KKD1,KKD2
C-----Calculation of cartesian local coordinates of the fragment
C-----from the internal coordinates
      CALL OH(RD,DT,W,T,XO1)
C-----Calculation of the orthogonalization matrix
C-----Orthogonalization of atomic coordinates
      DO 2 I=1,3
      DO 2 J=1,3
      XO(I,J)=0.0
      DO 2 K=1,3
2     XO(I,J)=XO(I,J)+O(J,K)*X(I,K)
C-----Transfer of the origin
      DO 5 I=1,4
      DO 5 J=1,3
5     XO2(I,J)=XO1(I,J)+XO(1,J)
C-----Rotate the system of axes
      CALL ROTMAT(XO,R)
      DO 8 I=1,3
      DO 8 K=1,3
      XO4(I,K)=0.0
      DO 8 J=1,3
8     XO4(I,K)=XO4(I,K)+R(K,J)*XO(I,J)
      DO 3 I=1,3
3     DX(I)=XO4(1,I)-XO2(1,I)
      DO 4 I=1,4
      DO 4 J=1,3
4     XO2(I,J)=XO2(I,J)+DX(J)
      CALL INVMAT3(R,RI,ICD)
      DO 7 I=1,4
      DO 7 K=1,3
      XO3(I,K)=0.0
      DO 7 J=1,3
7     XO3(I,K)=XO3(I,K)+RI(K,J)*XO2(I,J)
C-----Calculation of fractional coordinates of H
      CALL INVMAT3(O,OI,ICD)
      DO 6 I=1,4
      DO 6 K=1,3
      X3(I,K)=0.0
      DO 6 J=1,3
6     X3(I,K)=X3(I,K)+OI(K,J)*XO3(I,J)
C-----Print fractional coordinates
      DO 9 K=1,3
9     XHH(K)=X3(4,K)
      IF(KKD1.EQ.1)WRITE(15,104)(ATM1(4,L),L=1,6),(X3(4,K),K=1,3)
104   FORMAT(/'Atom',9X,'X/a',7X,'Y/b',7X,'Z/c'/6A1,3X,3F10.5/)
C     CALL DSTANG4(ATM1,XO3,RD)
      RETURN
      END
 
      SUBROUTINE DECEQP(STR,LE,F)
C-----This routine decodes the equivalent positions
      CHARACTER*1 STR(20),COO(13)
      INTEGER LE(3,3)
      DIMENSION F(3)
      COMMON/INPOU/IN,IO,IR
      DATA COO/'X','Y','Z',',','/','+','-','1','2','3','4','5','6'/
      DO 1 I=1,3
      F(I)=0.0
      DO 1 J=1,3
1     LE(I,J)=0
      N=0
      DO 3 I=1,20
      IF(STR(I).EQ.' '.AND.STR(I+1).EQ.' ')GO TO 4
      IF(STR(I).NE.COO(4))GO TO 3
      N=N+1
      IF(N.EQ.1)NA=I
      IF(N.EQ.2)NB=I
3     CONTINUE
4     DO 17 K=1,3
      GO TO (5,6,7),K
5     N1=1
      N2=NA-1
      GO TO 8
6     N1=NA+1
      N2=NB-1
      GO TO 8
7     N1=NB+1
      N2=20
8     KS=1
      KD=0
      DO 16 I=N1,N2
      IF(STR(I).EQ.' '.AND.STR(I+1).EQ.' ')RETURN
      IF(STR(I).EQ.' ')GO TO 16
      DO 9 J=1,13
      IF(STR(I).EQ.COO(J))GO TO 10
9     CONTINUE
      GO TO 16
10    IF(J.EQ.4)GO TO 16
      IF(J.EQ.7)KS=-1
      IF(J.EQ.5)KD=1
      IF(J.GT.3)GO TO 12
      DO 11 L=1,3
11    IF(K.EQ.L)LE(L,J)=1*KS
      KS=1
      GO TO 16
12    IF(J.LT.8)GO TO 16
      IF(KD.EQ.1)GO TO 14
      DO 13 L=1,3
13    IF(K.EQ.L)F(L)=(J-7)*KS
      KS=1
      GO TO 16
14    DO 15 L=1,3
      IF(K.EQ.L)F(L)=F(L)/(J-7)
15    CONTINUE
      KD=0
16    CONTINUE
17    CONTINUE
      RETURN
      END
 
      SUBROUTINE DISCON(AD2,DM,NL,AAAA)
C-----Set data and calculates contacts involving atoms translated in
C-----equivalent positions
      CHARACTER*1 STPE(38,20),AD2(6),AAAA(500,6)
      INTEGER LE1(24,3,3),FT(3),LE(3,3)
      DIMENSION F1(24,3),FB(3),F(3)
      COMMON/COOA/X1(3),X2(3)
     1      /LBL3/IC,NE,NT
     2      /GRP1/C(3),A(3)
     3      /INS4/CS(3),SN(3)
     4      /EQVP/LE1,F1,FB
     6      /LBL4/STPE
     7      /LBL8/XX(500,3)
     8      /INPOU/IN,IO,IR
      KM=2*NT+1
      DO 3 I=1,NE
      DO 3 I1=1,KM
      FT(1)=I1-NT-1
      DO 3 I2=1,KM
      FT(2)=I2-NT-1
      DO 3 I3=1,KM
      FT(3)=I3-NT-1
      DO 1 J=1,3
      F(J)=F1(I,J)+FB(J)+FT(J)
      DO 1 K=1,3
1     LE(J,K)=LE1(I,J,K)
      CALL DSTCHB(DM,AD2,LE,F,NL,AAAA)
      IF(IC.NE.-1) GO TO 3
      DO 2 J=1,3
      F(J)=-F1(I,J)+FB(J)+FT(J)
      DO 2 K=1,3
2     LE(J,K)=-LE1(I,J,K)
      CALL DSTCHB(DM,AD2,LE,F,NL,AAAA)
3     CONTINUE
      RETURN
      END
 
      SUBROUTINE DIST(X1,X2,D)
C-----This routine calculates the distance between two atoms
      DIMENSION X1(3),X2(3)
      D=0.0
      DO 2 K=1,3
      AA=ABS(X1(K)-X2(K))
      IF(AA.GT.1E-9)GO TO 1
      AA=0.0
      GO TO 2
1     AA=AA*AA
2     D=D+AA
      IF(D.GT.1E-2)D=SQRT(D)
      RETURN
      END
 
      SUBROUTINE DLTLN(LL,KK)
C-----Skip plane, line, puckering, stereo lines
      CHARACTER*1 TIT(4)
      COMMON/INPOU/IN,IO,IR
      DO 1 L=1,LL
      IF(KK.EQ.0)READ(IN,'(I2)')NP
      IF(KK.EQ.1)READ(IN,'(1X,I2)')NP
      AA=NP
      BB=AA/12.0
      IF(BB.GT.1.00)THEN
      IF(BB.LE.2.00)NP=1
      IF(BB.GT.2.00.AND.BB.LE.3.00)NP=2
      IF(BB.GT.3.00.AND.BB.LE.4.00)NP=3
      IF(BB.GT.4.00.AND.BB.LE.5.00)NP=4
      IF(BB.GT.5.00)NP=5
      ELSE
      GO TO 1
      ENDIF
      DO 2 J=1,NP
2     READ(IN,'(4A1)')TIT
1     CONTINUE
      RETURN
      END
 
      SUBROUTINE DMLN(NH1,NTN,LSP,LSL,NPR,NST)
C-----Skip lines giving unnecessary data
      CHARACTER*1 TIT(4)
      COMMON/INPOU/IN,IO,IR
      NT=NH1+NTN
      IF(NT.NE.0)THEN
      DO 1 J=1,NT
1     READ(IN,'(4A1)')TIT
      ENDIF
      NT=LSP+LSL+NPR
      IF(NT.NE.0)THEN
      KK=0
      CALL DLTLN(NT,KK)
      ENDIF
      IF(NST.NE.0)THEN
      KK=1
      CALL DLTLN(NST,KK)
      ENDIF
      RETURN
      END
 
      SUBROUTINE DSTCHB(DM,AD2,LE,F,NL,AAAA)
C-----This routine calculates interatomic distances involving
C-----atoms translated in equivalent positions
      CHARACTER*1 AD2(6),AAAA(500,6)
      INTEGER LE(3,3)
      DIMENSION F(3),X(3),D(3)
      COMMON/COOA/X1(3),X2(3)
     1      /LBL3/IC,NE,NT
     2      /GRP1/C(3),A(3)
     3      /INS4/CS(3),SN(3)
     4      /LBL8/XX(500,3)
     5      /INPOU/IN,IO,IR
      DO 1 I=1,3
      X(I)=F(I)
      DO 1 J=1,3
1     X(I)=X(I)+LE(I,J)*X2(J)
      DO 2 I=1,3
2     IF(ABS(X(I)-X2(I)).GT.1E-4)GO TO 3
      GO TO 8
3     DA=0
      DO 4 I=1,3
      D(I)=X1(I)-X(I)
      D(I)=D(I)*C(I)
4     DA=DA+D(I)**2
      IF(DA.LT.1E-2)GO TO 8
      DA=DA+2*(D(1)*D(2)*CS(3)+D(1)*D(3)*CS(2)+D(2)*D(3)*CS(1))
      IF(DA.GT.1E-2)DA=SQRT(DA)
      IF(DA.GT.DM) GO TO 8
C-----Check for duplicates
      DO 11 I=1,NL
      DO 9 L=1,6
      IF(AD2(L).NE.AAAA(NL,L))GO TO 11
9     CONTINUE
      DO 10 K=1,3
      IF(X(K).NE.XX(NL,K))GO TO 11
10    CONTINUE
      GO TO 8
11    CONTINUE
      NL=NL+1
      IF(NL.LE.500)GO TO 5
      WRITE(IO,101)
101   FORMAT('There are too many atoms whose environment is ',
     1 'requested'/' Please, reduce their number'/)
      NL=500
      GO TO 8
5     DO 6 K=1,3
6     XX(NL,K)=X(K)
      DO 7 L=1,6
7     AAAA(NL,L)=AD2(L)
8     RETURN
      END
 
      SUBROUTINE DSTN(KY,NNS,LA1,LE,F,NS4,ATT,ATM4,X4,DS4,SSTR)
C-----This routine calculates an interatomic distance: DA
C-----atoms translated in equivalent positions
      CHARACTER*1 ATT(6),ATM4(20,6),STR(38),SSTR(20,38),AT1(2),LBA(2),
     1 LA1(2,20)
      CHARACTER*2 AT2
      INTEGER LE(3,3)
      DIMENSION F(3),X(3),D(3),X4(50,3),DS4(50)
      EQUIVALENCE (AT2,AT1(1))
      COMMON/COOA/X1(3),X2(3)
     2      /GRP1/C(3),A(3)
     3      /INS4/CS(3),SN(3)
     5      /INPOU/IN,IO,IR
      DO 1 I=1,3
      X(I)=F(I)
      DO 1 J=1,3
1     X(I)=X(I)+LE(I,J)*X2(J)
      CALL LABL(ATT,KY,NNS,LA1,LBA,KKD)
      AT1(1)=LBA(1)
      AT1(2)=LBA(2)
      DO 2 I=1,3
2     IF(ABS(X(I)-X2(I)).GT.1E-4)GO TO 3
      RETURN
3     DA=0
      DO 4 I=1,3
      D(I)=X1(I)-X(I)
      D(I)=D(I)*C(I)
4     DA=DA+D(I)**2
      IF(DA.LT.1E-2)RETURN
      DA=DA+2*(D(1)*D(2)*CS(3)+D(1)*D(3)*CS(2)+D(2)*D(3)*CS(1))
      DA=SQRT(DA)
      IF(DA.GT.3.8)RETURN
      IF((AT2.EQ.'O '.OR.AT2.EQ.'N '.OR.AT2.EQ.'F ').AND.DA.LT.3.5)
     1 GO TO 10
      IF((AT2.EQ.'CL'.OR.AT2.EQ.'S '.OR.AT2.EQ.'BR'.OR.AT2.EQ.'I ').
     1 AND.DA.LT.3.8)GO TO 10
      RETURN
10    NS4=NS4+1
      DS4(NS4)=DA
      CALL EQPDEC(LE,F,STR)
      DO 9 LL=1,38
9     SSTR(NS4,LL)=STR(LL)
      DO 5 K=1,3
5     X4(NS4,K)=X(K)
      DO 6 L=1,6
6     ATM4(NS4,L)=ATT(L)
      DO 7 L=1,6
      IF(ATM4(NS4,L).EQ.' ')THEN
      ATM4(NS4,L)='"'
      GO TO 8
      ENDIF
7     CONTINUE
8     RETURN
      END
 
      SUBROUTINE ELIMRA(NL,AAAA)
C-----Eliminate repeated atoms
      CHARACTER*1 AAAA(500,6),AAA1(500,6)
      DIMENSION X1(500,3)
      COMMON/LBL8/XX(500,3)
      DO 6 I=1,NL-1
      DO 5 J=I+1,NL
      DO 1 L=1,6
      IF(AAAA(I,L).NE.AAAA(J,L))GO TO 5
1     CONTINUE
      DO 2 K=1,3
      IF(ABS(XX(I,K)-XX(J,K)).GT.0.0001)GO TO 5
2     CONTINUE
      DO 3 L=1,6
3     AAAA(J,L)=' '
      DO 4 K=1,3
4     XX(J,K)=0.0
5     CONTINUE
6     CONTINUE
      NN=0
      DO 11 I=1,NL
      DO 7 L=1,6
      IF(AAAA(I,1).NE.' ') GO TO 8
7     CONTINUE
      GO TO 11
8     NN=NN+1
      DO 9 L=1,6
9     AAA1(NN,L)=AAAA(I,L)
      DO 10 K=1,3
10    X1(NN,K)=XX(I,K)
11    CONTINUE
      DO 13 I=1,NN
      DO 12 L=1,6
12    AAAA(I,L)=AAA1(I,L)
      DO 13 K=1,3
13    XX(I,K)=X1(I,K)
      NL=NN
      RETURN
      END
 
      SUBROUTINE ENERMOL(KD,Y0,NG1,NG3,XO1,XO3,ET)
C-----This routine calculates the interatomic contacts and the
C-----van der Waals and Coulombic energy, EC, with a gaussian correction
C-----for hydrogen bonding [EC=EC*(1+Y0*EXP(-CG*(d-1.8)**2))], where
C-----Y0=5.0, CG=ln2/wt=2.773, wt=1/2(half-height amplitude)
      CHARACTER*1 AT(500,6)
      CHARACTER*2 KS(2),SDL2(500)
      DIMENSION XO1(500,3),XO3(500,3),X1(3),X2(3)
      COMMON/GRP4/AT
     1      /GRP5/MC(400,8)
     2      /GRP6/NA1(500),NA3(500)
     3      /GRP7/Q1(500),Q3(500)
     4      /GRG1/SDL2
     5      /INPOU/IN,IO,IR
      ET=0
      IF(KD.EQ.0)THEN
      DO 5 I=1,NG3-1
      DO 5 J=I+1,NG3
      DO 2 K=1,8
      IF(NA3(J).EQ.MC(NA3(I),K))GO TO 5
      DO 1 L=1,8
      IF(NA3(J).EQ.MC(MC(NA3(I),K),L))GO TO 5
1     CONTINUE
2     CONTINUE
      DO 3 K=1,3
      X1(K)=XO3(I,K)
3     X2(K)=XO3(J,K)
      CALL DIST(X1,X2,D)
      IF(D.LT.0.1)GO TO 5
      IF(D.GT.7.0)GO TO 5
      KS(1)=SDL2(NA3(I))
      KS(2)=SDL2(NA3(J))
      CALL ENERPOT(KS,D,E,IND)
      IF(IND.EQ.1)E=0
      EQ=(Q3(I)*Q3(J)/D)*331.979
      ED=0.0
      ET=ET+E
5     CONTINUE
      ELSEIF(KD.EQ.1)THEN
      DO 12 I=1,NG3
      DO 12 J=3,NG1
      DO 10 K=1,3
      X1(K)=XO3(I,K)
10    X2(K)=XO1(J,K)
      CALL DIST(X1,X2,D)
      IF(D.LT.0.1)GO TO 12
      IF(D.GT.7.0)GO TO 12
      DIE=1.0+Y0*EXP(-2.773*(D-1.8)**2)
      KS(1)=SDL2(NA3(I))
      KS(2)=SDL2(NA1(J))
      KKAA=1
      IF(KS(1).EQ.'O '.OR.
     1   KS(1).EQ.'N '.OR.KS(1).EQ.'F '.OR.KS(1).EQ.'S '.OR.
     2   KS(1).EQ.'SI'.OR.KS(1).EQ.'Si'.OR.KS(1).EQ.'P '.OR.
     3   KS(1).EQ.'CL'.OR.KS(1).EQ.'Cl'.OR.KS(1).EQ.'BR'.OR.
     4   KS(1).EQ.'Br'.OR.KS(1).EQ.'I ')KKAA=0
      CALL ENERPOT(KS,D,E,IND)
      IF(IND.EQ.1)E=0
      EQ=(Q3(I)*Q1(J)/D)*331.979
      EQ=EQ*DIE
      ED=0.0
      IF(D.LT.2.0.AND.KKAA.EQ.1)ED=1.0*(EXP(-10.0*(D-2.0))-1)
      ET=ET+E+EQ+ED
12    CONTINUE
      ENDIF
      RETURN
      END
 
      SUBROUTINE ENERPOT(KS,D,E,IND)
C-----This routine calculates the van der Waals potential energy for a
C-----non bonding contact between two atoms (only H, C, N, O, F, Si, P,
C-----S, Cl, Fe, Br, Sn, I, Pb).
C-----Coefficients from K.Mirsky (Computing in Crystallography, 1978)
C-----for  H, C, N, O, F, S, Cl; from P. Johnston et al. (J.Cryst.
C-----Spectr.Res.,18, 403, 1988) for P, Fe I; from Ahmed et al. (Acta
C-----Cryst. B27, 867, 1971) for Si, Sn, Pb; from Eliel et al. (Con-
C-----formational Analysis, Chap.7, Interscience Publ.,1965) for Br.
      CHARACTER*2 ATM(14),KS(2)
      DIMENSION COA(14),COB(14),COC(14)
      DATA ATM/'H ','C ','N ','O ','F ','SI','P ','S ','CL','FE','BR',
     1 'SN','I ','PB'/
      DATA COA/4900,71600,42000,77700,42000,42000,525963,235000,4580,
     1 70301,359377,42000,38994,42000/
      DATA COB/4.29,3.68,3.78,4.18,4.15,2.95,4.0,3.49,2.262,3.18,3.484,
     1 2.85,2.5,3.06/
      DATA COC/29,421,259,259.4,148,1160,1531.54,2346,2980,277.8,
     1 3436.3,1401,12451,912/
      IF(D.EQ.0.)RETURN
      IND=0
      DO 1 I=1,14
      IF(KS(1).EQ.ATM(I)) GO TO 2
1     CONTINUE
      IND=1
      RETURN
2     K1=I
      DO 3 I=1,14
      IF(KS(2).EQ.ATM(I)) GO TO 4
3     CONTINUE
      IND=1
      RETURN
4     K2=I
      A=COA(K1)
      B=COB(K1)
      C=COC(K1)
      F=(COC(K1)+COC(K2))/COC(1)
      IF(K1.NE.K2) GO TO 5
      GO TO 6
5     A=SQRT(COA(K1)*COA(K2))
      B=(COB(K1)+COB(K2))/2
      C=SQRT(COC(K1)*COC(K2))
6     E=A*EXP(-B*D)-C*D**(-6)+F*D**(-12)
      RETURN
      END
 
      SUBROUTINE ENRPRF(KY,BRV,RD,N,Y0,JO1,ISS,NNH,ATR1,ACT,XOCT)
C-----Calculates the non bonded energy profile for rotation of a
C-----hydroxyl group about the C-O bond, then find the best orientation
C-----of it according to the position of the minimum in the profile
      CHARACTER*1 BRV,STPE(38,20),SDL1(2,500),LA1(2,20),LB1(6,2),
     1 LBL(6,40),AT(500,6),ACT(6),ATR1(3,6)
      CHARACTER*2 LA(20),SDL2(500)
      DIMENSION XOCT(3),FB(3),F1(24,3),LE1(24,3,3),XP(3),YP(3),CP(3),
     1 XON(3,3),XRO(3,3),XR(3,3)
      EQUIVALENCE (LA(1),LA1(1,1)),(LBL(1,1),LB1(1,1)),(SDL2(1),
     1 SDL1(1,1))
      COMMON/GRP1/C(3),A(3)
     1      /GRP3/XO(500,3)
     2      /GRP4/AT
     3      /GRP10/LA
     5      /GRP12/NNS
     6      /GRP5/MC(400,8)
     7      /CDO/KHW,KKD1,KKD2
      COMMON/COOR/X(500,3),Q(500)
     1      /NEW1/O(3,3),OI(3,3)
     2      /INS4/CS(3),SN(3)
     3      /LBL4/STPE
     4      /LBL3/IC,NE,NT
     5      /INPOU/IN,IO,IR
     6      /COOA/X1(3),X2(3)
     7      /EQVP/LE1,F1,FB
     9      /PRF1/TH(100),DEJ1(100)
     1      /PRF2/NN1
      IF(KKD1.EQ.1.AND.KHW.EQ.0)WRITE(IO,'(/19(''*'')/''* HYDROXYL '',
     1 6A1,'' *''/19(''*'')/)')(AT(JO1,L),L=1,6)
      IF(KKD1.EQ.1.AND.KHW.EQ.1)WRITE(IO,'(/19(''*'')/''* WATER    '',
     1 6A1,'' *''/19(''*'')/)')(AT(JO1,L),L=1,6)
      KDK=0
67    NTT=N
      CALL ROTH(KY,BRV,NTT,KDK,Y0,ATR1)
C-----Finds the minimum of the profile
      CALL MINPROF(NN1,TH,DEJ1,XP,YP)
      IF(KDK.EQ.0)THEN
      IF(KKD1.EQ.1)WRITE(IO,'(''Non bonded energy profile for '',
     1 ''rotation about the '',6A1,''-'',6A1,''bond''/)')
     2 (AT(ISS,L),L=1,6),(AT(JO1,L),L=1,6)
      IF(KKD1.EQ.1)THEN
      DO 212 I=1,NN1/3
212   WRITE(IO,'(3(F6.1,F11.3,6X))')TH(I),DEJ1(I),TH(I+24),DEJ1(I+24),
     1 TH(I+48),DEJ1(I+48)
      WRITE(IO,'(46X,F6.1,F11.3)')TH(73),DEJ1(73)
      ENDIF
      ENDIF
      DO 69 K=1,3
69    XP(K)=XP(K)+180
      CALL MINPAR(XP,YP,CP,XM,YM,ICD)
      IF(KKD1.EQ.1)WRITE(IO,'(''('',I1,'')'',5X,''Minimum of profile''/
     1 )')KDK
      IF(ICD.EQ.1)THEN
      IF(KKD1.EQ.1)WRITE(15,'(''Xm='',F8.2/''Ym='',F8.2/)')XP(2),YP(2)
      GO TO 59
      ELSE
      IF(KKD1.EQ.1)WRITE(IO,'(''Xm='',F8.2,3X,''Ym='',F8.2)')XM,YM
      ENDIF
59    CONTINUE
      NF=3
      FI=0.75*(XM/RD)
      DO 60 K=1,3
      XON(1,K)=XO(ISS,K)
      XON(2,K)=XO(JO1,K)
60    XON(3,K)=XO(NNH,K)
      CALL ROTFRG(NF,XON,FI,XRO)
      DO 61 I=1,NF
      DO 61 J=1,3
      XR(I,J)=0.0
      DO 61 K=1,3
61    XR(I,J)=XR(I,J)+OI(J,K)*XRO(I,K)
      FI=FI*RD
      IF(KKD1.EQ.1)WRITE(IO,105)(AT(ISS,L),L=1,6),(AT(JO1,L),L=1,6),FI
105   FORMAT(/'After rotation about ',6A1,'-',6A1,' by',F7.1,' deg'/
     1 12X,'x',8X,'y',8X,'z')
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,'''',1X,3F9.5)')(AT(NNH,L),L=1,6),
     1 (XR(3,K),K=1,3)
64    DO 66 K=1,3
      X(NNH,K)=XR(3,K)
66    XO(NNH,K)=XRO(3,K)
C-----Hydrogen bonding parameters
      IF(KKD1.EQ.1)WRITE(IO,'(/''Hydrogen bonding parameters with '',
     1 ''the refined H coordinates'')')
      CALL HBPAR(RD,JO1,NNH,ACT,XOCT)
      IF(ABS(FI).LT.0.2)GO TO 68
      KDK=KDK+1
      IF(KDK.EQ.9)GO TO 68
      GO TO 67
68    RETURN
      END
 
      SUBROUTINE EQPBRV(BRV,KB,B)
C-----Decodes the equivalent position strings and takes account
C-----of the Bravais translation lattice
      CHARACTER*1 STR(20),BRV,STPE(38,20)
      INTEGER LE1(24,3,3),LE(3,3)
      DIMENSION B(3,3),F1(24,3),F(3),FB(3)
      COMMON/LBL3/IC,NE,NT
     1      /LBL4/STPE
     2      /EQVP/LE1,F1,FB
      DO 1 I=1,24
      DO 1 J=1,3
      F1(I,J)=0.0
      DO 1 K=1,3
1     LE1(I,J,K)=0
      IF(NE.EQ.1)GO TO 4
      DO 3 I=2,NE
      DO 2 L=1,20
2     STR(L)=STPE(I-1,L)
      CALL DECEQP(STR,LE,F)
      DO 3 J=1,3
      F1(I,J)=F(J)
      LE1(1,J,J)=1
      DO 3 K=1,3
3     LE1(I,J,K)=LE(J,K)
      GO TO 6
4     DO 5 J=1,3
5     LE1(1,J,J)=1
6     CONTINUE
      DO 7 I=1,3
      DO 7 J=1,3
7     B(I,J)=0.0
      KB=0
      IF(BRV.EQ.'P') RETURN
      IF(BRV.NE.'A') GO TO 8
      B(1,2)=1./2.
      B(1,3)=1./2.
      KB=1
      RETURN
8     IF(BRV.NE.'B') GO TO 9
      B(1,1)=1./2.
      B(1,3)=1./2.
      KB=1
      RETURN
9     IF(BRV.NE.'C') GO TO 10
      B(1,1)=1./2.
      B(1,2)=1./2.
      KB=1
      RETURN
10    IF(BRV.NE.'I') GO TO 11
      B(1,1)=1./2.
      B(1,2)=1./2.
      B(1,3)=1./2.
      KB=1
      RETURN
11    IF(BRV.NE.'R') GO TO 12
      B(1,1)=1./3.
      B(2,2)=1./3.
      B(2,3)=1./3.
      B(1,2)=2*B(1,1)
      B(1,3)=2*B(1,1)
      B(2,1)=2*B(1,1)
      KB=2
      RETURN
12    B(1,2)=1./2.
      B(1,3)=1./2.
      B(2,1)=1./2.
      B(2,3)=1./2.
      B(3,1)=1./2.
      B(3,2)=1./2.
      KB=3
      RETURN
      END
 
      SUBROUTINE EQPDEC(LE,F,STR)
C-----This routine encodes the matrix and the vector describing
C-----an equivalent position into an alphanumeric string
      CHARACTER*1 STR(38),COO(6)
      INTEGER NT(3),LE(3,3),LE1(3)
      DIMENSION F(3),F1(3)
      DATA COO/'X','Y','Z','1','2','3'/
      DO 12 I=1,38
12    STR(I)=' '
      STR(13)=','
      STR(26)=','
      DO 1 I=1,3
      NT(I)=F(I)
1     F1(I)=F(I)-NT(I)
      DO 8 I=1,3
      IF(I.EQ.1) N=0
      IF(I.EQ.2) N=13
      IF(I.EQ.3) N=26
      DO 2 J=1,3
      LE1(J)=LE(I,J)
      IF(LE1(J).LT.0)STR(2*J-1+N)='-'
      IF(LE1(J).GT.0)STR(2*J-1+N)='+'
      IF(IABS(LE1(J)).EQ.1)STR(2*J+N)=COO(J)
2     CONTINUE
      IF(ABS(F(I)).LT.1E-4) GO TO 8
      IF(ABS(F1(I)).LT.1E-4) GO TO 6
      IF(F(I).LT.0.)STR(7+N)='-'
      IF(F(I).LT.0.)STR(11+N)='-'
      IF(F(I).GT.0.)STR(7+N)='+'
      IF(F(I).GT.0.)STR(11+N)='+'
      STR(9+N)='/'
      IF(ABS(F1(I)).GT.0.6) GO TO 3
      STR(8+N)='1'
      IF(ABS(F1(I)).LT.0.2)STR(10+N)='6'
      IF(ABS(F1(I)).LT.0.3.AND.ABS(F1(I)).GT.0.2)STR(10+N)='4'
      IF(ABS(F1(I)).LT.0.4.AND.ABS(F1(I)).GT.0.3)STR(10+N)='3'
      IF(ABS(F1(I)).LT.0.6.AND.ABS(F1(I)).GT.0.4)STR(10+N)='2'
      GO TO 6
3     IF(ABS(F1(I)).GT.0.7) GO TO 4
      STR(8+N)='2'
      STR(10+N)='3'
      GO TO 6
4     IF(ABS(F1(I)).GT.0.8) GO TO 5
      STR(8+N)='3'
      STR(10+N)='4'
      GO TO 6
5     STR(8+N)='5'
      STR(10+N)='6'
6     IF(NT(I).EQ.0) GO TO 8
      DO 7 J=1,3
      IF(IABS(NT(I)).EQ.J.AND.NT(I).GT.0)STR(11+N)='+'
      IF(IABS(NT(I)).EQ.J.AND.NT(I).LT.0)STR(11+N)='-'
7     IF(IABS(NT(I)).EQ.J)STR(12+N)=COO(J+3)
8     CONTINUE
      IF(STR(12).EQ.' ') STR(11)=' '
      IF(STR(25).EQ.' ') STR(24)=' '
      IF(STR(38).EQ.' ') STR(37)=' '
      DO 11 K=1,4
      DO 10 I=1,37
      IF(STR(I).NE.' ') GO TO 10
      DO 9 J=I,37
9     STR(J)=STR(J+1)
      STR(38)=' '
10    CONTINUE
11    CONTINUE
      RETURN
      END
 
      SUBROUTINE HBPAR(RD,JO1,N,ACT,XOCT)
C-----Calculates and prints the H-bonding geometric parameters
      CHARACTER*1 AT(500,6),ACT(6)
      DIMENSION XD(3,3),XB(3),XC(3),XOCT(3),RIS(3)
      COMMON/GRP3/XO(500,3)
     1      /GRP4/AT
     2      /INPOU/IN,IO,IR
     3      /CDO/KHW,KKD1,KKD2
      DO 1 K=1,3
      XD(1,K)=XO(JO1,K)
      XD(2,K)=XO(N,K)
1     XD(3,K)=XOCT(K)
      DO 4 K=1,3
      XB(K)=XD(1,K)
4     XC(K)=XD(2,K)
      CALL DIST(XB,XC,DK)
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,''-'',6A1,6X,F9.3)')
     1 (AT(JO1,L),L=1,6),(AT(N,L),L=1,6),DK
      DO 2 K=1,3
      XB(K)=XD(1,K)
2     XC(K)=XD(3,K)
      CALL DIST(XB,XC,DK)
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,''...'',6A1,6X,F7.3)')
     1 (AT(JO1,L),L=1,6),(ACT(L),L=1,6),DK
      DO 3 K=1,3
      XB(K)=XD(2,K)
3     XC(K)=XD(3,K)
      CALL DIST(XB,XC,DK)
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,''...'',6A1,6X,F7.3)')
     1 (AT(N,L),L=1,6),(ACT(L),L=1,6),DK
      CALL BNDAN(RD,XD,T)
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,''-'',6A1,''-'',6A1,F8.1)')
     1 (AT(JO1,L),L=1,6),(AT(N,L),L=1,6),(ACT(L),L=1,6),T
      DO 5 K=1,3
      RIS(K)=XD(1,K)
      XD(1,K)=XD(2,K)
5     XD(2,K)=RIS(K)
      CALL BNDAN(RD,XD,T)
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,''-'',6A1,''-'',6A1,F8.1)')
     1 (AT(N,L),L=1,6),(AT(JO1,L),L=1,6),(ACT(L),L=1,6),T
      DO 6 K=1,3
      RIS(K)=XD(2,K)
      XD(2,K)=XD(3,K)
6     XD(3,K)=RIS(K)
      CALL BNDAN(RD,XD,T)
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,''-'',6A1,''-'',6A1,F8.1/)')
     1 (AT(N,L),L=1,6),(ACT(L),L=1,6),(AT(JO1,L),L=1,6),T
      RETURN
      END
 
      SUBROUTINE HBPAR1(RD,XD,DT,ANG)
C-----Calculates the three distances and angles involved in a
C-----hydrogen bond
      Dimension XD(3,3),DT(3),ANG(3),XB(3),XC(3),RIS(3)
      COMMON/INPOU/IN,IO,IR
      RD=57.295779513
      DO 1 K=1,3
      XB(K)=XD(1,K)
1     XC(K)=XD(2,K)
      CALL DIST(XB,XC,DT1)
      DT(1)=DT1
      DO 2 K=1,3
      XB(K)=XD(1,K)
2     XC(K)=XD(3,K)
      CALL DIST(XB,XC,DT1)
      DT(2)=DT1
      DO 3 K=1,3
      XB(K)=XD(2,K)
3     XC(K)=XD(3,K)
      CALL DIST(XB,XC,DT1)
      DT(3)=DT1
      CALL BNDAN(RD,XD,T)
      ANG(1)=T
      DO 4 K=1,3
      RIS(K)=XD(1,K)
      XD(1,K)=XD(2,K)
4     XD(2,K)=RIS(K)
      CALL BNDAN(RD,XD,T)
      ANG(2)=T
      DO 5 K=1,3
      RIS(K)=XD(2,K)
      XD(2,K)=XD(3,K)
5     XD(3,K)=RIS(K)
      CALL BNDAN(RD,XD,T)
      ANG(3)=T
      RETURN
      END
 
      SUBROUTINE HCOO(D13,XO,XH,XOH)
C-----This routine calculates the coordinates of atom X lying on the
C-----line joining atom A to atom B at a requested distance from atom A
C-----and nearest to atom B
      DIMENSION XO(3,3),DDC(3),X1(3),X2(3),XH(3),XOH(3),X(3,3)
      COMMON/NEW1/O(3,3),OI(3,3)
      D12=0.0
      DO 1 K=1,3
      DDC(K)=XO(2,K)-XO(1,K)
1     D12=D12+DDC(K)**2
      D12=SQRT(D12)
      DO 2 K=1,3
      DDC(K)=(DDC(K)/D12)*D13
      X1(K)=XO(1,K)+DDC(K)
2     X2(K)=XO(1,K)-DDC(K)
      DH1=0.0
      DH2=0.0
      DO 3 K=1,3
      DH1=DH1+(XO(2,K)-X1(K))**2
3     DH2=DH2+(XO(2,K)-X2(K))**2
      DH1=SQRT(DH1)
      DH2=SQRT(DH2)
      DO 4 K=1,3
      XO(3,K)=X1(K)
4     IF(DH1.GT.DH2)XO(3,K)=X2(K)
      DH=0.0
      DO 5 K=1,3
5     DH=DH+(XO(1,K)-XO(3,K))**2
      DH=SQRT(DH)
      DO 6 J=1,3
      X(3,J)=0.0
      DO 6 K=1,3
6     X(3,J)=X(3,J)+OI(J,K)*XO(3,K)
      DO  7 K=1,3
      XH(K)=X(3,K)
7     XOH(K)=XO(3,K)
      RETURN
      END
 
      SUBROUTINE HDRB(N,SDL2,NH2,IA)
C-----Finds the hydrogen bonds
      CHARACTER*1 AT(500,6),STR(38)
      CHARACTER*2 SDL2(500)
      DIMENSION XD(3,3),DT(3),ANG(3),XP(3),XOT(3),IA(2,200),F(3)
      INTEGER FT(3),LE(3,3)
      COMMON/GRP3/XO(500,3)
     1      /GRP4/AT
     2      /COOR/X(500,3),Q(500)
     3      /EQVP/LE1(24,3,3),F1(24,3),FB(3)
     5      /LBL3/IC,NE,NT
     6      /NEW1/O(3,3),OI(3,3)
     7      /INPOU/IN,IO,IR
      KM=2*NT+1
      DO 12 J=1,NH2
      DO 11 I=1,N
      IF(SDL2(I).EQ.'O '.OR.SDL2(I).EQ.'N '.OR.SDL2(I).EQ.'F '.OR.
     1 SDL2(I).EQ.'CL'.OR.SDL2(I).EQ.'BR'.OR.SDL2(I).EQ.'I '.OR.
     2 SDL2(I).EQ.'S ')GO TO 1
      GO TO 11
1     CONTINUE
      DO 10 L=1,NE
      DO 10 I1=1,KM
      FT(1)=I1-NT-1
      DO 10 I2=1,KM
      FT(2)=I2-NT-1
      DO 10 I3=1,KM
      IKD=0
      FT(3)=I3-NT-1
      DO 2 JJ=1,3
      F(JJ)=F1(L,JJ)+FB(JJ)+FT(JJ)
      DO 2 K=1,3
2     LE(JJ,K)=LE1(L,JJ,K)
3     DO 4 LL=1,3
      XP(LL)=F(LL)
      DO 4 JJ=1,3
4     XP(LL)=XP(LL)+LE(LL,JJ)*X(I,JJ)
      DO 5 JJ=1,3
      XOT(JJ)=0.0
      DO 5 KK=1,3
5     XOT(JJ)=XOT(JJ)+O(JJ,KK)*XP(KK)
      D=0.0
      DO 6 KK=1,3
6     D=D+(XO(IA(1,J),KK)-XOT(KK))**2
      D=SQRT(D)
      IF(D.GT.4.0.OR.D.LT.1.5)GO TO 8
      DO 7 K=1,3
      XD(1,K)=XO(IA(2,J),K)
      XD(2,K)=XO(IA(1,J),K)
7     XD(3,K)=XOT(K)
      CALL HBPAR1(RD,XD,DT,ANG)
      IF(ANG(1).LT.90)GO TO 8
      IF(DT(1).GT.1.2.OR.DT(3).GT.4.0)GO TO 8
      IF(SDL2(IA(2,J)).EQ.'C '.AND.DT(3).GT.2.5)GO TO 8
      IF((SDL2(I).EQ.'S '.OR.SDL2(I).EQ.'CL').AND.DT(3).GT.3.05)
     1 GO TO 8
      IF(SDL2(I).EQ.'Br'.AND.DT(3).GT.3.15) GO TO 8
      IF(SDL2(I).EQ.'I '.AND.DT(3).GT.3.35) GO TO 8
      IF((SDL2(I).EQ.'O '.OR.SDL2(I).EQ.'N '.OR.SDL2(I).EQ.'F ')
     1 .AND.DT(3).GT.2.7)GO TO 8
      CALL EQPDEC(LE,F,STR)
      WRITE(IO,'(6A1,''-'',6A1,''...'',6A1,3(F6.3,3X),3(F7.2,1X))')
     1 (AT(IA(2,J),LL),LL=1,6),(AT(IA(1,J),LL),LL=1,6),
     2 (AT(I,LL),LL=1,6),DT,ANG
      IF(STR(1).EQ.'+')STR(1)=' '
      DO 13 LL=1,38
      IF(STR(LL).EQ.'X')STR(LL)='x'
      IF(STR(LL).EQ.'Y')STR(LL)='y'
      IF(STR(LL).EQ.'Z')STR(LL)='z'
      IF(STR(LL).EQ.','.AND.STR(LL+1).EQ.'+')STR(LL+1)=' '
13    CONTINUE
      WRITE(IO,'(6A1,''at '',38A1/)')(AT(I,LL),LL=1,6),STR
8     IF(IKD.EQ.1)GO TO 10
      IF(IC.NE.-1)GO TO 10
      DO 9 JJ=1,3
      F(JJ)=-F1(L,JJ)+FB(JJ)+FT(JJ)
      DO 9 K=1,3
9     LE(JJ,K)=-LE1(L,JJ,K)
      IKD=1
      GO TO 3
10    CONTINUE
11    CONTINUE
12    CONTINUE
      RETURN
      END
 
      SUBROUTINE INTERC(N,DM,NL,BRV,NR,AFR,AAAA)
C-----This routine calculates all the interatomic contacts,
C-----involving atoms of different equivalent positions in the same
C-----or in adjacent cells, less than a given maximum value
      CHARACTER*1 BRV,STPE(38,20),AT(500,6),AD2(6),AFR(3,6),
     1 AAAA(500,6)
      INTEGER LE1(24,3,3)
      DIMENSION B(3,3),F1(24,3),FB(3)
      COMMON/GRP1/C(3),A(3)
     1      /LBL3/IC,NE,NT
     2      /LBL4/STPE
     3      /LBL8/XX(500,3)
     4      /INS4/CS(3),SN(3)
     5      /GRP4/AT
     6      /COOR/X(500,3),Q(500)
     7      /EQVP/LE1,F1,FB
     9      /COOA/X1(3),X2(3)
     1      /INPOU/IN,IO,IR
      CALL EQPBRV(BRV,KB,B)
      DO 22 II=1,NR
      DO 21 I=1,N
      DO 14 L=1,6
      IF(AT(I,L).NE.AFR(II,L))GO TO 21
14    CONTINUE
      DO 20 J=1,N
      DO 15 L=1,6
15    AD2(L)=AT(J,L)
      DO 16 K=1,3
      X1(K)=X(I,K)
16    X2(K)=X(J,K)
      DO 17 L=1,3
17    FB(L)=0.0
      CALL DISCON(AD2,DM,NL,AAAA)
      IF(KB.EQ.0)GO TO 20
      DO 19 L=1,KB
      DO 18 K=1,3
18    FB(K)=B(L,K)
      CALL DISCON(AD2,DM,NL,AAAA)
19    CONTINUE
20    CONTINUE
21    CONTINUE
22    CONTINUE
      RETURN
      END
 
      SUBROUTINE INPOUT
C-----Assignes the input and output files
      CHARACTER*10 STR1
      CHARACTER*8,STR2
      CHARACTER*14 FILIN,FILCHAR
      COMMON/INPOU/IN,IO,IR
      STR1='Key in the'
      STR2='put file'
      WRITE(*,'(1X,A10,'' in'',A8/'' ?>'')')STR1,STR2
      READ(*,'(A14)')FILIN
      OPEN(UNIT=IN,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=IN)
      WRITE(*,'(1X,A10,'' out'',A8/'' ?>'')')STR1,STR2
      READ(*,'(A14)')FILIN
      OPEN(UNIT=IO,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=IO)
      WRITE(*,'(1X,A10,'' name of the file of the atom charges''/
     1 '' ?>'')')STR1
      READ(*,'(A14)')FILCHAR
      OPEN(UNIT=IR,FILE=FILCHAR,FORM='FORMATTED')
      REWIND(UNIT=IR)
      RETURN
      END
 
      SUBROUTINE INVMAT3(A1,B1,ICD)
C-----Inverts a 3x3 matrix by partition method
      DOUBLE PRECISION A(3,3),B(3,3),D,AM1,AM2,AM3,AM4,X
      DIMENSION A1(3,3),B1(3,3)
      ICD=0
      DO 3 I=1,3
      DO 3 J=1,3
3     A(I,J)=A1(I,J)
      D=A(1,1)*A(2,2)-A(1,2)*A(2,1)
      IF(DABS(D).LT.1D-10) GO TO 1
      AM1=A(1,2)*A(2,3)-A(1,3)*A(2,2)
      AM2=A(1,3)*A(2,1)-A(1,1)*A(2,3)
      AM3=A(2,1)*A(3,2)-A(2,2)*A(3,1)
      AM4=A(1,2)*A(3,1)-A(1,1)*A(3,2)
      X=A(3,3)+(AM1*A(3,1)+AM2*A(3,2))/D
      IF(DABS(X).LT.1D-10) GO TO 1
      B(3,3)=1D0/X
      X=X*D
      B(1,3)=AM1/X
      B(1,1)=(B(1,3)*AM3+A(2,2))/D
      B(1,2)=(B(1,3)*AM4-A(1,2))/D
      B(2,3)=AM2/X
      B(2,1)=(B(2,3)*AM3-A(2,1))/D
      B(2,2)=(B(2,3)*AM4+A(1,1))/D
      B(3,1)=AM3/X
      B(3,2)=AM4/X
      DO 4 I=1,3
      DO 4 J=1,3
4     B1(I,J)=B(I,J)
      GO TO 2
1     WRITE(15,100)
100   FORMAT(' The matrix cannot be inverted')
      ICD=1
C     STOP
2     RETURN
      END
 
      SUBROUTINE LABLAT(CHR1,CHR2)
C-----This routine transforms the second character of an atom label
C-----in a low case character when this is part of the atomic symbol
      CHARACTER*1 MINU(26),MAIU(26),CHR1,CHR2,LA1(2,20)
      CHARACTER*2 LA(20)
      EQUIVALENCE (LA1(1,1),LA(1))
      COMMON/MIMA/MINU,MAIU
     1      /GRP10/LA
     2      /GRP12/NNS
      DO 1 I=1,NNS
      IF(CHR1.EQ.LA1(1,I).AND.CHR2.EQ.LA1(2,I))GO TO 2
1     CONTINUE
      RETURN
2     DO 3 J=1,26
      IF(CHR2.EQ.MAIU(J)) CHR2=MINU(J)
3     CONTINUE
      RETURN
      END
 
      SUBROUTINE LABL(AAAA,KY,NS,LA1,LBA,KKD)
C-----This routine finds the atom species from its label
      CHARACTER*1 AAAA(6),LA1(2,20),LBA(2),MINU(26),MAIU(26),MINUS(2)
      COMMON/INPOU/IN,IO,IR
     1      /MIMA/MINU,MAIU
      KKD=0
      DO 7 I=1,NS
      MINUS(1)=LA1(1,I)
      MINUS(2)=LA1(2,I)
      CALL MINUMA(2,MINUS)
      LA1(1,I)=MINUS(1)
7     LA1(2,I)=MINUS(2)
      IF(KY.EQ.1)THEN
      DO 1 J=1,NS
      IF(AAAA(1).EQ.LA1(1,J).AND.AAAA(2).EQ.LA1(2,J))GO TO 3
1     CONTINUE
      DO 2 K=1,NS
      IF(AAAA(1).EQ.LA1(1,K))GO TO 4
2     CONTINUE
      WRITE(IO,100)AAAA
100   FORMAT('Atom ',6A1,' is not among the given species'/
     1 'The first two characters are assumed')
      LBA(1)=AAAA(1)
      LBA(2)=AAAA(2)
      KKD=1
      RETURN
3     LBA(1)=LA1(1,J)
      LBA(2)=LA1(2,J)
      RETURN
4     LBA(1)=LA1(1,K)
      LBA(2)=' '
      RETURN
      ENDIF
      LBA(1)=AAAA(1)
      DO 5 I=1,26
      IF(AAAA(2).EQ.MAIU(I))GO TO 6
5     CONTINUE
      LBA(2)=' '
      RETURN
6     LBA(2)=AAAA(2)
      RETURN
      END
 
      SUBROUTINE MINPAR(X,Y,C,XM,YM,ICD)
C-----This routine calculates the coordinates of the minimum of a
C-----three-points parabole
      DIMENSION X(3),Y(3),C(3),A(3,3),B(3,3)
C-----Sets up matrix A
      DO 2 I=1,3
      A(I,1)=X(I)**2
      A(I,2)=X(I)
2     A(I,3)=1.D0
C-----Inverts matrix A
      CALL INVMAT3(A,B,ICD)
      IF(ICD.EQ.1)THEN
      XM=X(2)-180
      YM=Y(2)
      GO TO 4
      ENDIF
C-----Determines coefficients of the parabolic equation
      DO 3 I=1,3
      C(I)=0.0
      DO 3 J=1,3
3     C(I)=C(I)+B(I,J)*Y(J)
      XM=-C(2)/(2*C(1))
      YM=C(1)*XM**2+C(2)*XM+C(3)
      XM=XM-180
4     RETURN
      END
 
      SUBROUTINE MINPROF(N,XX,YY,X,Y)
C-----Read the profile and determine the minimum
      DIMENSION XX(100),YY(100),X(3),Y(3)
      COMMON/INPOU/IN,IO,IR
      YM=YY(1)
      ID=1
      DO 1 I=1,N
      IF(YY(I).LT.YM)THEN
      YM=YY(I)
      ID=I
      ENDIF
1     CONTINUE
      IF(ID.EQ.1)THEN
      X(1)=-185
      X(2)=-180
      X(3)=-175
      Y(1)=YY(2)
      Y(2)=YY(1)
      Y(3)=YY(2)
      RETURN
      ENDIF
      DO 2 I=1,3
      X(I)=XX(ID-2+I)
2     Y(I)=YY(ID-2+I)
      RETURN
      END
 
      SUBROUTINE MINUMA(N,STR)
C-----This routine transforms the lower case characters into upper
C-----case characters, in a string of N characters
      CHARACTER*1 MINU(26),MAIU(26),STR(N)
      COMMON/MIMA/MINU,MAIU
      DO 2 I=1,N
      DO 1 J=1,26
      IF(STR(I).NE.MINU(J)) GO TO 1
      STR(I)=MAIU(J)
      GO TO 2
1     CONTINUE
2     CONTINUE
      RETURN
      END
 
      SUBROUTINE NEWCOO(NG,XNO,BI,XOO,RR)
C-----Modifies coordinates to have them referred to a Cartesian system
C-----of axis with the origin at the first atom, with the z axis along
C-----the first to the second atom, the y axis perpendicular to the z
C-----axis and lying on the plane defined by the first three atoms, and
C-----the x axis perpendicular to y and z.
      DIMENSION XNO(500,3),X1(500,3),XOO(500,3),RR(3)
      DOUBLE PRECISION V,XR(3),B(3,3),BI(3,3)
      DO 1 I=1,NG
      DO 1 J=1,3
1     X1(I,J)=XNO(I,J)-XNO(1,J)
C-----Building up matrix B
      DO 2 I=1,3
      XR(I)=X1(2,I)
2     RR(I)=XNO(1,I)
      CALL NORM(XR)
      DO 3 I=1,3
3     B(3,I)=XR(I)
      V=0
      DO 4 K=1,3
4     V=V+X1(3,K)*B(3,K)
      DO 5 K=1,3
5     XR(K)=X1(3,K)-V*B(3,K)
      CALL NORM(XR)
      DO 6 I=1,3
6     B(2,I)=XR(I)
      XR(1)=B(2,2)*B(3,3)-B(3,2)*B(2,3)
      XR(2)=B(3,1)*B(2,3)-B(2,1)*B(3,3)
      XR(3)=B(2,1)*B(3,2)-B(3,1)*B(2,2)
      CALL NORM(XR)
      DO 7 I=1,3
7     B(1,I)=XR(I)
      DO 8 I=1,NG
      DO 8 J=1,3
      XOO(I,J)=0.0
      DO 8 K=1,3
8     XOO(I,J)=XOO(I,J)+B(J,K)*X1(I,K)
      BI(1,1)=B(2,2)*B(3,3)-B(2,3)*B(3,2)
      BI(2,1)=B(3,1)*B(2,3)-B(2,1)*B(3,3)
      BI(3,1)=B(2,1)*B(3,2)-B(3,1)*B(2,2)
      BI(1,2)=B(3,2)*B(1,3)-B(1,2)*B(3,3)
      BI(2,2)=B(1,1)*B(3,3)-B(3,1)*B(1,3)
      BI(3,2)=B(3,1)*B(1,2)-B(1,1)*B(3,2)
      BI(1,3)=B(1,2)*B(2,3)-B(2,2)*B(1,3)
      BI(2,3)=B(2,1)*B(1,3)-B(1,1)*B(2,3)
      BI(3,3)=B(1,1)*B(2,2)-B(1,2)*B(2,1)
      RETURN
      END
 
      SUBROUTINE NORM(X)
C-----This routine normalizes a 3-dimensional vector
      DOUBLE PRECISION X(3)
      S=0.0
      DO 1 I=1,3
1     S=S+X(I)**2
      IF(S.GT.1E-2)S=SQRT(S)
      DO 2 I=1,3
2     X(I)=X(I)/S
      RETURN
      END
 
      SUBROUTINE NORM1(X)
C-----This routine normalizes a 3-dimensional vector
      DIMENSION X(3)
      S=0.0
      DO 1 I=1,3
1     S=S+X(I)**2
      IF(S.GT.1E-2)S=SQRT(S)
      DO 2 I=1,3
2     X(I)=X(I)/S
      RETURN
      END
 
      SUBROUTINE OH(RD,D,W,T,XO)
C-----This routine calculates the local orthogonal coordinates
C-----in a system of four atoms whose internal coordinates are
C-----given
      DIMENSION Q(0:2,0:100),A(0:8,0:100),XO(4,3),D(0:3),W(0:3),T(0:3)
      DO 1 K=1,7
1     A(K,0)=0.
      A(0,0)=1.
      A(4,0)=1.
      A(8,0)=1.
      Q(0,0)=0.
      Q(1,0)=0.
      Q(2,0)=0.
      DO 2 I=0,3
      N=I
      IF(I.EQ.0) GO TO 6
      IF(I.GT.2) GO TO 4
      IF(I.LT.2) GO TO 3
      T(I)=0.0
      M=1
      GO TO 5
3     W(I)=0.
      T(I)=180.0-T(I)
      M=0
      GO TO 5
4     M=N-1
5     N=M
      B0=-COS(W(I)/RD)
      B1=-SIN(W(I)/RD)
      B5=-SIN(T(I)/RD)
      B8=COS(T(I)/RD)
      B3=-B1*B8
      B4=B0*B8
      B6=B1*B5
      B7=-B0*B5
      B2=0.
      A(0,I)=A(0,M)*B0+A(1,M)*B3+A(2,M)*B6
      A(1,I)=A(0,M)*B1+A(1,M)*B4+A(2,M)*B7
      A(2,I)=A(0,M)*B2+A(1,M)*B5+A(2,M)*B8
      A(3,I)=A(3,M)*B0+A(4,M)*B3+A(5,M)*B6
      A(4,I)=A(3,M)*B1+A(4,M)*B4+A(5,M)*B7
      A(5,I)=A(3,M)*B2+A(4,M)*B5+A(5,M)*B8
      A(6,I)=A(6,M)*B0+A(7,M)*B3+A(8,M)*B6
      A(7,I)=A(6,M)*B1+A(7,M)*B4+A(8,M)*B7
      A(8,I)=A(6,M)*B2+A(7,M)*B5+A(8,M)*B8
      Q(0,I)=Q(0,M)+A(0,I)*D(I)
      Q(1,I)=Q(1,M)+A(3,I)*D(I)
      Q(2,I)=Q(2,M)+A(6,I)*D(I)
6     DO 7 K=1,3
7     XO(I+1,K)=Q(K-1,I)
2     CONTINUE
      RETURN
      END
 
      SUBROUTINE ORTHO(N,RD)
C-----Orthogonalize atomic coordinates
      COMMON/INPOU/IN,IO,IR
     1      /GRP1/C(3),A(3)
     2      /COOR/X(500,3),Q(500)
     3      /GRP3/XO(500,3)
     4      /NEW1/O(3,3),OI(3,3)
     5      /INS4/CS(3),SN(3)
C-----Calculation of the orthogonalization matrix
      DO 1 I=1,3
      AA=A(I)/RD
      SN(I)=SIN(AA)
1     CS(I)=COS(AA)
      CSAR=(CS(2)*CS(3)-CS(1))/(SN(2)*SN(3))
      SNAR=SQRT(1.0-CSAR**2)
      O(1,1)=C(1)
      O(1,2)=C(2)*CS(3)
      O(1,3)=C(3)*CS(2)
      O(2,1)=0.0
      O(2,2)=C(2)*SN(3)
      O(2,3)=-C(3)*SN(2)*CSAR
      O(3,1)=0.0
      O(3,2)=0.0
      O(3,3)=C(3)*SN(2)*SNAR
C-----Orthogonalization of atomic coordinates
      DO 2 I=1,N
      DO 2 J=1,3
      XO(I,J)=0.0
      DO 2 K=1,3
2     XO(I,J)=XO(I,J)+O(J,K)*X(I,K)
      RETURN
      END
 
      SUBROUTINE ROTCOO(NG,XOO,FI,BI,XRO,RR)
C-----Rotates the coordinates of a fragment: XOO --> XRO
      DIMENSION XOO(500,3),XRO(500,3),X1(500,3),RT(3,3),RR(3)
      DOUBLE PRECISION BI(3,3)
C-----Building up rotation matrix
      RT(1,1)=COS(FI)
      RT(1,2)=-SIN(FI)
      RT(1,3)=0.0
      RT(2,1)=SIN(FI)
      RT(2,2)=COS(FI)
      RT(2,3)=0.0
      RT(3,1)=0.0
      RT(3,2)=0.0
      RT(3,3)=1.0
C-----Applying rotation matrix
      DO 1 I=1,NG
      DO 1 J=1,3
      X1(I,J)=0.0
      DO 1 K=1,3
1     X1(I,J)=X1(I,J)+RT(J,K)*XOO(I,K)
      DO 2 I=1,NG
      DO 2 J=1,3
      XRO(I,J)=0.0
      DO 2 K=1,3
2     XRO(I,J)=XRO(I,J)+BI(J,K)*X1(I,K)
      DO 3 I=1,NG
      DO 3 J=1,3
      XRO(I,J)=XRO(I,J)+RR(J)
3     CONTINUE
      RETURN
      END
 
      SUBROUTINE ROTFRG(N,X,FI,XO)
C-----This routine modifies the coordinates so to have them referred
C-----to a system of orthogonal axes with the origin at the first atom,
C-----the z axis directed from the first to the second atom, the y
C-----axis normal to z and lying in the plane defined by the first
C-----three atoms and the x axis normal to y and z.
      DIMENSION X(3,3),XO(3,3),XR(3),X1(3,3),X2(3,3),B(3,3),
     1 RT(3,3),BI(3,3)
      DO 1 I=1,N
      DO 1 J=1,3
1     X1(I,J)=X(I,J)-X(1,J)
C-----Building up matrix B
      XR(1)=X1(2,1)
      XR(2)=X1(2,2)
      XR(3)=X1(2,3)
      CALL NORM1(XR)
      DO 2 I=1,3
2     B(3,I)=XR(I)
      V=0
      DO 3 K=1,3
3     V=V+X1(3,K)*B(3,K)
      DO 4 K=1,3
4     XR(K)=X1(3,K)-V*B(3,K)
      CALL NORM1(XR)
      DO 5 I=1,3
5     B(2,I)=XR(I)
      XR(1)=B(2,2)*B(3,3)-B(3,2)*B(2,3)
      XR(2)=B(3,1)*B(2,3)-B(2,1)*B(3,3)
      XR(3)=B(2,1)*B(3,2)-B(3,1)*B(2,2)
      CALL NORM1(XR)
      DO 6 I=1,3
6     B(1,I)=XR(I)
      DO 7 I=1,N
      DO 7 J=1,3
      X2(I,J)=0.0
      DO 7 K=1,3
7     X2(I,J)=X2(I,J)+B(J,K)*X1(I,K)
C-----Building up rotation matrix
      RT(1,1)=COS(FI)
      RT(1,2)=-SIN(FI)
      RT(1,3)=0.0
      RT(2,1)=SIN(FI)
      RT(2,2)=COS(FI)
      RT(2,3)=0.0
      RT(3,1)=0.0
      RT(3,2)=0.0
      RT(3,3)=1.0
C-----Applying rotation matrix
      DO 8 I=1,N
      DO 8 J=1,3
      XO(I,J)=0.0
      DO 8 K=1,3
8     XO(I,J)=XO(I,J)+RT(J,K)*X2(I,K)
      CALL INVMAT3(B,BI,ICD)
      DO 9 I=1,N
      DO 9 J=1,3
      X1(I,J)=0.0
      DO 9 K=1,3
9     X1(I,J)=X1(I,J)+BI(J,K)*XO(I,K)
      DO 10 I=1,N
      DO 10 J=1,3
10    XO(I,J)=X1(I,J)+X(1,J)
      RETURN
      END
 
      SUBROUTINE ROTH(KY,BRV,N,KDK,Y0,ATR1)
C-----Calculates the energy profile by rotating H about the C-O bond
      CHARACTER*1 AT(500,6),ATR1(3,6),BRV,
     1 AAA1(500,6),AFR1(3,6),STPE(38,20),SDL1(2,500),LA1(2,20),
     2 LBA(2),CTEMP(6)
      CHARACTER*14 FILGRF
      CHARACTER*2 SDL2(500),LA(20)
      DIMENSION XRO1(500,3),RR1(3),XO1(500,3),XO3(500,3),XOO1(500,3),
     1 XF1(500,3),XXO(500,3),QF(500),LE1(24,3,3),F1(24,3),FB(3)
      DOUBLE PRECISION BI1(3,3)
      EQUIVALENCE (LA(1),LA1(1,1)),(SDL2(1),SDL1(1,1))
      COMMON/NEW1/O(3,3),OI(3,3)
     1      /GRP4/AT
     2      /GRP5/MC(400,8)
     3      /GRP3/XO(500,3)
     4      /GRP6/NA1(500),NA3(500)
     5      /GRP7/Q1(500),Q3(500)
     7      /COOR/X(500,3),Q(500)
     8      /GRP1/C(3),A(3)
     9      /LBL3/IC,NE,NT
      COMMON/LBL4/STPE
     1      /INPOU/IN,IO,IR
     2      /GRG1/SDL2
     3      /GRP10/LA
     4      /EQVP/LE1,F1,FB
     6      /PRF1/TH(100),DEJ1(100)
     7      /PRF2/NN1
     8      /GRP12/NNS
     9      /INS4/CS(3),SN(3)
     1      /CDO/KHW,KKD1,KKD2
C-----Constants and I/O files
      RD=0.0174532925
C-----Rotation parameters
      KC=1
      KDE=1
      NG1=3
      T1=-180.0
      T3=180.0
      DT1=5.0
C-----Crystal environment of the rotating hydrogen
      NR1=NG1-2
      DO 6 I=3,NG1
      DO 6 L=1,6
6     AFR1(I-2,L)=ATR1(I,L)
      CALL CONTACT(BRV,N,NR1,AFR1,NL,AAA1,XXO,QF)
      T1=T1*RD
      T3=T3*RD
      DT1=DT1*RD
C-----Rotating hydrogen
      DO 11 I=1,NG1
      DO 10 J=1,N
      DO 8 L=1,6
      IF(ATR1(I,L).NE.AT(J,L)) GO TO 10
8     CONTINUE
      DO 9 K=1,3
      XO1(I,K)=XO(J,K)
9     XF1(I,K)=X(J,K)
      NA1(I)=J
      Q1(I)=Q(J)
      GO TO 11
10    CONTINUE
      WRITE(*,109)(ATR1(I,L),L=1,6)
109   FORMAT(' Atom ',6A1,' does not exist in the molecule.'/
     1 ' Check your data. Calculation must stop.')
      STOP
11    CONTINUE
C-----Unrotating fragment
      KL=0
      DO 19 J=1,N
      DO 16 I=3,NG1
      IF(NA1(I).EQ.J)GO TO 19
16    CONTINUE
      KL=KL+1
      NA3(KL)=J
      Q3(KL)=Q(J)
      DO 18 K=1,3
18    XO3(KL,K)=XO(J,K)
19    CONTINUE
      NG3=KL
C-----Add crystal environment to the fix fragment
      CALL ADD(N,NL,AAA1,XO3,XXO,Q3,QF,NG3)
C-----Add environmental atoms to the atom symbol matrix
      DO 62 I=1,N
      DO 61 L=1,6
61    CTEMP(L)=AT(I,L)
      CALL LABL(CTEMP,KY,NNS,LA1,LBA,KKD)
      SDL1(1,I)=LBA(1)
62    SDL1(2,I)=LBA(2)
C-----Calculation of the molecular energy for the unrotated molecule
      KD=0
      CALL ENERMOL(KD,Y0,NG1,NG3,XO1,XO3,E00)
      KD=1
      CALL ENERMOL(KD,Y0,NG1,NG3,XO1,XO3,E01)
      E0=E00+E01
      E0J=E0*4.184
      DEMAX=2000
      CALL NEWCOO(NG1,XO1,BI1,XOO1,RR1)
      NE1=0
      NN1=0
      NX=0
C-----Start the loop of the first fragment rotations
25    NE1=NE1+1
      CALL ROTCOO(NG1,XOO1,T1,BI1,XRO1,RR1)
      NN1=NN1+1
      KD=1
      CALL ENERMOL(KD,Y0,NG1,NG3,XRO1,XO3,E1)
      TH(NN1)=T1/RD
      EJ=E1*4.184
      E01J=E01*4.184
      DE=E1-E01
      DEJ1(NN1)=EJ-E01J
      IF(DEJ1(NN1).GT.DEMAX)DEJ1(NN1)=DEMAX
      DE=DEJ1(NN1)/4.184
      T1=T1+DT1
      IF(T1.GT.T3)GO TO 34
      GO TO 25
C-----End of the rotation loop
34    IF(KDK.EQ.0.AND.KKD2.EQ.1)THEN
      WRITE(*,137)
137   FORMAT(' Key in the name of the file for GRAFIT'/' ?>')
      READ(*,'(A14)')FILGRF
      OPEN(UNIT=17,FILE=FILGRF,FORM='FORMATTED')
      REWIND(UNIT=17)
      DO 55 I=1,NN1
55    WRITE(17,138)TH(I),DEJ1(I)
138   FORMAT(F11.1,1X,F15.3)
      CLOSE(17)
      ENDIF
      RETURN
      END
 
      SUBROUTINE ROTMAT(X,B1)
C-----This routine built up a matrix which transform the orthogonal
C-----coordinates of three atoms 1,2,3 in a new orthogonal system
C-----having atom 1 at the origin, x along the 1-2 direction,
C-----y orthogonal to x and in the 1,2,3 plane, z normal to the x
C-----and y axes
      DOUBLE PRECISION B(3,3),RIS
      DIMENSION X(4,3),X1(3,3),R2(3),R3(3),B1(3,3)
      DO 1 I=1,3
      DO 1 J=1,3
1     X1(I,J)=X(I,J)-X(1,J)
      X1(2,1)=-X1(2,1)
      X1(2,2)=-X1(2,2)
      X1(2,3)=-X1(2,3)
      AA=0.0
      DO 2 K=1,3
2     AA=AA+X1(2,K)**2
      IF(AA.GT.1E-2)AA=SQRT(AA)
      DO 3 K=1,3
3     B(3,K)=X1(2,K)/AA
      V=0.0
      DO 4 K=1,3
4     V=V+X1(3,K)*B(3,K)
      DO 5 K=1,3
5     R2(K)=V*B(3,K)
      DO 6 K=1,3
6     R3(K)=X1(3,K)-R2(K)
      AA=0.0
      DO 7 K=1,3
7     AA=AA+R3(K)**2
      IF(AA.GT.1E-2)AA=SQRT(AA)
      DO 8 K=1,3
8     B(2,K)=R3(K)/AA
      B(1,1)=B(2,2)*B(3,3)-B(2,3)*B(3,2)
      B(1,2)=B(2,3)*B(3,1)-B(2,1)*B(3,3)
      B(1,3)=B(2,1)*B(3,2)-B(2,2)*B(3,1)
      DO 9 K=1,3
      RIS=-B(1,K)
      B(1,K)=B(3,K)
9     B(3,K)=RIS
      DO 10 I=1,3
      DO 10 J=1,3
10    B1(I,J)=B(I,J)
      RETURN
      END
 
      SUBROUTINE STARTH(KY,BRV,RD,N,Y0,HD,JO1,O1L,ATH,XHH,ISS,ACT,XOCT,
     1 KKR)
C-----Considers the environment of a hydroxyl oxygen and finds its
C-----coordinates
      CHARACTER*1 BRV,STPE(38,20),SDL1(2,500),LA1(2,20),LB1(6,2),
     1 LBL(6,40),AT(500,6),ATT(6),ATM4(20,6),ACT(6),ATM1(4,6),
     2 O1L(6),SSTR(20,38),ATH(6)
      CHARACTER*1 AT5(50,6),SST5(20,38)
      CHARACTER*2 LA(20),SDL2(500)
      DIMENSION IS1(50),DS1(50),IS2(50),DS2(50),IS3(50),DS3(50),
     1 XD(3,3),R(4,3),DS4(50),XHH(3),X4(50,3),XCT(3),XOCT(3),
     2 B(3,3),FB(3),F1(24,3),LE1(24,3,3),D(4),W(4),T(4),XAT(4,3)
      DIMENSION DS5(50),X5(50,3),XO5(50,3),AN5(50)
      EQUIVALENCE (LA(1),LA1(1,1)),(LBL(1,1),LB1(1,1)),(SDL2(1),
     1 SDL1(1,1))
      COMMON/GRP1/C(3),A(3)
     1      /GRP3/XO(500,3)
     2      /GRP4/AT
     3      /GRP10/LA
     5      /GRP12/NNS
     6      /CDO/KHW,KKD1,KKD2
     7      /GRG1/SDL2
     8      /LTT/KB,B
      COMMON/COOR/X(500,3),Q(500)
     1      /NEW1/O(3,3),OI(3,3)
     2      /INS4/CS(3),SN(3)
     3      /LBL4/STPE
     4      /LBL3/IC,NE,NT
     5      /INPOU/IN,IO,IR
     6      /COOA/X1(3),X2(3)
     7      /EQVP/LE1,F1,FB
     9      /PRF1/TH(100),DEJ1(100)
     1      /PRF2/NN1
C-----Calculates distances of O1 from the atoms of the set
C-----and choose those in the 1.2/1.6 and 2.4/3.0 A ranges
C-----NS1=number, IS1(10)=positions in atom tab., DS1(10)=distances
C-----in the first range; NS2, IS2, DS2 the same in the second range
      NS1=0
      NS2=0
      DO 2 I=1,N
      IF(I.EQ.JO1)GO TO 2
      DT=0.0
      DO 1 K=1,3
1     DT=DT+(XO(I,K)-XO(JO1,K))**2
      IF(DT.LE.1.0.OR.DT.GE.12.25)GO TO 2
      IF(DT.LT.2.56)THEN
      NS1=NS1+1
      IS1(NS1)=I
      DS1(NS1)=SQRT(DT)
      ELSEIF(DT.GT.5.76.AND.(SDL2(I).EQ.'O '.OR.SDL2(I).EQ.'N '.OR.
     1 SDL2(I).EQ.'F '.OR.SDL2(I).EQ.'CL'.OR.SDL2(I).EQ.'S '))THEN
      NS2=NS2+1
      IS2(NS2)=I
      DS2(NS2)=SQRT(DT)
      ENDIF
2     CONTINUE
      IF(NS1.EQ.0)THEN
      WRITE(IO,'(''There is no atom bound to '',6A1)')(AT(JO1,L),L=1,6)
      KKR=1
      RETURN
      ENDIF
      IF(NS2.EQ.0)THEN
      IF(KKD1.EQ.1)WRITE(IO,'(''There are no '',
     1 ''acceptable contacts of the O...O or O...N or O...F''/
     2 ''type in the input set'')')
      ENDIF
      IF(KKD1.EQ.1)THEN
      DO 3 I=1,NS1
3     WRITE(IO,'(6A1,''-'',6A1,5X,F7.3,'' A'')')(AT(JO1,L),L=1,6),
     1 (AT(IS1(I),L),L=1,6),DS1(I)
      ENDIF
      IF(NS2.NE.0)THEN
      DO 4 I=1,NS2
4     IF(KKD1.EQ.1)WRITE(IO,'(6A1,''...'',6A1,3X,F7.3)')
     1 (AT(JO1,L),L=1,6),(AT(IS2(I),L),L=1,6),DS2(I)
      ENDIF
C-----Choose the most probable O1-C1 distance
      DD=DS1(1)
      IM1=1
      DO 5 I=1,NS1
      IF(ABS(DS1(I)-1.4).LT.ABS(DD-1.4))THEN
      DD=DS1(I)
      IM1=I
      ENDIF
5     CONTINUE
C-----Consider the intermolecular contacts
      CALL EQPBRV(BRV,KB,B)
      DO 7 K=1,3
7     X1(K)=X(JO1,K)
      NS4=0
      DO 13 J=1,N
      DO 8 L=1,6
8     ATT(L)=AT(J,L)
      DO 9 K=1,3
9     X2(K)=X(J,K)
      DO 10 L=1,3
10    FB(L)=0.0
      CALL TRSL(KY,NNS,LA1,ATT,ATM4,NS4,X4,DS4,SSTR)
      IF(KB.EQ.0)GO TO 13
      DO 12 L=1,KB
      DO 11 K=1,3
11    FB(K)=B(L,K)
      CALL TRSL(KY,NNS,LA1,ATT,ATM4,NS4,X4,DS4,SSTR)
12    CONTINUE
13    CONTINUE
      NS5=0
      IF(NS2.EQ.0)GO TO 17
C-----Transfer intramolecular data into common files (NS5)
      DO 16 I=1,NS2
      IF(DS2(I).GT.3.8)GO TO 16
      NS5=NS5+1
      DO 14 L=1,6
14    AT5(NS5,L)=AT(IS2(I),L)
      DS5(NS5)=DS2(I)
      DO 15 K=1,3
      X5(NS5,K)=X(IS2(I),K)
15    XO5(NS5,K)=XO(IS2(I),K)
      SST5(NS5,1)=' '
      SST5(NS5,2)='x'
      SST5(NS5,3)=','
      SST5(NS5,4)='y'
      SST5(NS5,5)=','
      SST5(NS5,6)='z'
16    CONTINUE
17    IF(NS4.EQ.0)GO TO 22
C-----Transfer intermolecular data into common files (NS5)
      DO 21 I=1,NS4
      IF(DS4(I).GT.3.8)GO TO 21
      NS5=NS5+1
      DO 18 L=1,6
18    AT5(NS5,L)=ATM4(I,L)
      DS5(NS5)=DS4(I)
      DO 19 K=1,3
      X5(NS5,K)=X4(I,K)
      XO5(NS5,K)=0.0
      DO 19 J=1,3
19    XO5(NS5,K)=XO5(NS5,K)+O(K,J)*X4(I,J)
      DO 20 L=1,38
      IF(SSTR(I,1).EQ.'+')SSTR(I,1)=' '
      IF(SSTR(I,L).EQ.'+'.AND.SSTR(I,L-1).EQ.',')SSTR(I,L)=' '
      IF(SSTR(I,L).EQ.'X')SSTR(I,L)='x'
      IF(SSTR(I,L).EQ.'Y')SSTR(I,L)='y'
      IF(SSTR(I,L).EQ.'Z')SSTR(I,L)='z'
20    SST5(NS5,L)=SSTR(I,L)
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,3X,20A1,F9.4)')(AT5(NS5,L),L=1,6),
     1 (SST5(NS5,L),L=1,20),DS5(NS5)
21    CONTINUE
22    CONTINUE
      IF(NS5.EQ.0)THEN
      WRITE(IO,'('' There are no contacts suitable for determining '',
     1 ''the hydrogen coordinates'')')
      KKK=1
      RETURN
      ENDIF
      NAN=0
      ANDM=100
      DO 23 K=1,3
      XD(1,K)=XO(IS1(IM1),K)
23    XD(2,K)=XO(JO1,K)
      DO 25 I=1,NS5
      DO 24 K=1,3
24    XD(3,K)=XO5(I,K)
      CALL BNDAN(RD,XD,T1)
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,''-'',6A1,''-'',6A1,2X,F7.2)')
     1 (AT(IS1(IM1),L),L=1,6),(AT(JO1,L),L=1,6),(AT5(I,L),L=1,6),T1
      DTA=ABS(T1-109.5)
      CDT=COS(DTA*0.0174532925)
      ANDX=DS5(I)/CDT
      IF(KKD1.EQ.1)WRITE(IO,'(F7.3,F9.2,F7.2,F10.4)')DS5(I),T1,DTA,
     1 ANDX
C-----Choses the best angle C1-O1...O2 (nearest to 109.5) with the
C-----shortest Distance
      IF(ABS(ANDX).LT.ABS(ANDM))THEN
      ANDM=ANDX
      IANDM=I
      ANM=T1
      ENDIF
      NAN=NAN+1
      AN5(NAN)=T1
25    CONTINUE
      IF(KKD1.EQ.1)THEN
      WRITE(IO,'(''Best angle''/2F7.3,F10.4,F8.2)')
     1 DS1(IM1),DS5(IANDM),ANDM,ANM
      WRITE(IO,'(''NAN = '',I2)')NAN
      ENDIF
      DO 6 K=1,3
      XCT(K)=X5(IANDM,K)
6     XOCT(K)=XO5(IANDM,K)
      DO 37 L=1,6
37    ACT(L)=AT5(IANDM,L)
C-----Calculate distances of C1 from an atom that could be C2
C-----and choose one in the 1.1/1.6 A range
C-----NS3=number, IS3(10)=position in atom tab., DS3(10)=distances
      NS3=0
      DO 27 I=1,N
      IF(I.EQ.JO1)GO TO 27
      DT=0.0
      DO 26 K=1,3
26    DT=DT+(XO(I,K)-XO(IS1(IM1),K))**2
      IF(DT.LE.1.2.OR.DT.GE.2.6)GO TO 27
      NS3=NS3+1
      IS3(NS3)=I
      DS3(NS3)=SQRT(DT)
27    CONTINUE
      IF(NS3.EQ.0)THEN
      WRITE(IO,'(''There is no atom bound to '',6A1)')(AT(IS1(IM1),L),
     1 L=1,6)
      KKR=1
      RETURN
      ENDIF
      IF(KKD1.EQ.1)THEN
      WRITE(IO,'(/''Bonds formed by '',6A1)')(AT(IS1(IM1),L),L=1,6)
      DO 28 I=1,NS3
28    WRITE(IO,'(6A1,''-'',6A1,3X,F7.3)')(AT(IS1(IM1),L),L=1,6),
     1 (AT(IS3(I),L),L=1,6),DS3(I)
      ENDIF
C-----Chooses the most acceptable C1-C2 distance
      DD=DS1(1)
      IM3=1
      DO 29 I=1,NS3
      IF(ABS(DS3(I)-1.5).LT.ABS(DD-1.5))IM3=I
29    CONTINUE
C-----Calculates the C2-C1-O1 angle: T1
      DO 30 K=1,3
      XD(1,K)=XO(JO1,K)
      XD(2,K)=XO(IS1(IM1),K)
30    XD(3,K)=XO(IS3(IM3),K)
      CALL BNDAN(RD,XD,T1)
C-----Calculates the C1-O1...O2 angle: T2
      DO 31 K=1,3
      XD(1,K)=XO(IS1(IM1),K)
      XD(2,K)=XO(JO1,K)
31    XD(3,K)=XOCT(K)
      CALL BNDAN(RD,XD,T2)
C-----Calculates the C2-C1-O1...O2 torsion angle: TRS
      DO 32 K=1,4
      R(1,K)=XO(IS3(IM3),K)
      R(2,K)=XO(IS1(IM1),K)
      R(3,K)=XO(JO1,K)
32    R(4,K)=XOCT(K)
      CALL TORS(RD,R,TAU)
C-----Calculates the starting hydroxyl-H coordinates
      DO 33 L=1,6
      ATM1(1,L)=AT(IS3(IM3),L)
      ATM1(2,L)=AT(IS1(IM1),L)
      ATM1(3,L)=O1L(L)
33    ATM1(4,L)=' '
      ATM1(4,1)='H'
      ATM1(4,3)=O1L(1)
      ATM1(4,4)=O1L(2)
      ATM1(4,5)=O1L(3)
      ATM1(4,6)=O1L(4)
      DO 34 K=1,3
      XAT(1,K)=X(IS3(IM3),K)
      XAT(2,K)=X(IS1(IM1),K)
      XAT(3,K)=X(JO1,K)
34    XAT(4,K)=0.0
      D(1)=0.0
      D(2)=DS3(IM3)
      D(3)=DS1(IM1)
      D(4)=HD
      W(1)=0.0
      W(2)=0.0
      W(3)=T1
      W(4)=109.5
      DO 35 J=1,3
35    T(J)=0.0
      T(4)=TAU
      CALL COORDH(RD,ATM1,XAT,D,W,T,XHH)
      ISS=IS1(IM1)
      DO 36 L=1,6
36    ATH(L)=ATM1(4,L)
      RETURN
      END
 
      SUBROUTINE STARTW(KY,BRV,RD,N,Y0,HD,JO1,O1L,ATA,XA,XOA,ATW,
     1 XHN,XON,KKK)
C-----Considers the environment of a water oxygen and finds its
C-----coordinates
      CHARACTER*1 BRV,STPE(38,20),SDL1(2,500),LA1(2,20),LB1(6,2),
     1 LBL(6,40),AT(500,6),ATT(6),ATM4(20,6),O1L(6),SSTR(20,38),
     2 AT2(20,6),SST1(20,38),ATA(2,6),ATTH(2,6),ATW(3,6),LR,
     3 AT22(20,6),SST2(20,38),A1(20,6),S1(20,38)
      CHARACTER*2 LA(20),SDL2(500)
      DIMENSION IS1(50),DS1(50),DS2(50),XD(3,3),DS4(50),X4(50,3),
     1 B(3,3),FB(3),F1(24,3),LE1(24,3,3),AN2(40),XE(40,3),XOE(40,3),
     2 IAN(40,2),XA(2,3),XOA(2,3),XOO(3,3),XH(3),XOH(3),XHN(2,3),
     3 XON(2,3),XOW(3,3),XAQ(3),IANDM(2),DS22(50),XEE(40,3),DB(50),
     4 XOEE(40,3),XB(50,3),XOB(50,3),NA(50)
      EQUIVALENCE (LA(1),LA1(1,1)),(LBL(1,1),LB1(1,1)),(SDL2(1),
     1 SDL1(1,1))
      COMMON/GRP1/C(3),A(3)
     1      /GRP3/XO(500,3)
     2      /GRP4/AT
     3      /GRP10/LA
     4      /GRP12/NNS
     5      /GRG1/SDL2
     6      /LTT/KB,B
     7      /CDO/KHW,KKD1,KKD2
      COMMON/COOR/X(500,3),Q(500)
     1      /NEW1/O(3,3),OI(3,3)
     2      /INS4/CS(3),SN(3)
     3      /LBL4/STPE
     4      /LBL3/IC,NE,NT
     5      /INPOU/IN,IO,IR
     6      /COOA/X1(3),X2(3)
     7      /EQVP/LE1,F1,FB
     8      /PRF1/TH(100),DEJ1(100)
     9      /PRF2/NN1
C-----Calculates distances of O1 from the atoms of the set
C-----and choose those in the 1.0/3.8 A range
C-----NS1=number, IS1(10)=positions in atom tab., DS1(10)=distances
      KII=0
      NS1=0
      NS4=0
      DO 3 I=1,N
      IF(I.EQ.JO1)GO TO 3
      DT=0.0
      DO 1  K=1,3
1     DT=DT+(XO(I,K)-XO(JO1,K))**2
      DT=SQRT(DT)
      IF(DT.GE.3.8)GO TO 3
      IF((SDL2(I).EQ.'O '.OR.SDL2(I).EQ.'N '.OR.SDL2(I).EQ.'F ').AND.
     1 DT.LT.3.5)GO TO 2
      IF((SDL2(I).EQ.'CL'.OR.SDL2(I).EQ.'S '.OR.SDL2(I).EQ.'BR'.OR.
     1 SDL2(I).EQ.'I ').AND.DT.LT.3.8)GO TO 2
      GO TO 3
2     NS1=NS1+1
      IS1(NS1)=I
      DS1(NS1)=DT
3     CONTINUE
C-----Consider the intermolecular contacts
      CALL EQPBRV(BRV,KB,B)
      DO 4 K=1,3
4     X1(K)=X(JO1,K)
      DO 10 J=1,N
      DO 5 L=1,6
5     ATT(L)=AT(J,L)
      DO 6 K=1,3
6     X2(K)=X(J,K)
      DO 7 L=1,3
7     FB(L)=0.0
      CALL TRSL(KY,NNS,LA1,ATT,ATM4,NS4,X4,DS4,SSTR)
      IF(KB.EQ.0)GO TO 10
      DO 9 L=1,KB
      DO 8 K=1,3
8     FB(K)=B(L,K)
      CALL TRSL(KY,NNS,LA1,ATT,ATM4,NS4,X4,DS4,SSTR)
9     CONTINUE
10    CONTINUE
      NS2=0
C-----Blank SST1(20,38)
      DO 45 I=1,20
      DO 45 L=1,38
45    SST1(I,L)=' '
      IF(NS1.EQ.0)GO TO 14
C-----Transfer intramolecular data into common files (NS2)
      DO 13 I=1,NS1
      IF(DS1(I).GT.3.8)GO TO 13
      NS2=NS2+1
      NSS=NS2
      DO 11 L=1,6
      AT2(NS2,L)=AT(IS1(I),L)
11    AT22(NSS,L)=AT2(NS2,L)
      DS2(NS2)=DS1(I)
      DS22(NSS)=DS2(NS2)
      DO 12 K=1,3
      XE(NS2,K)=X(IS1(I),K)
      XEE(NSS,K)=XE(NS2,K)
      XOE(NS2,K)=XO(IS1(I),K)
12    XOEE(NSS,K)=XOE(NS2,K)
      SST1(NS2,1)=' '
      SST1(NS2,2)='x'
      SST1(NS2,3)=','
      SST1(NS2,4)='y'
      SST1(NS2,5)=','
      SST1(NS2,6)='z'
      DO 48 L=1,6
48    SST2(NSS,L)=SST1(NS2,L)
13    CONTINUE
14    IF(NS4.EQ.0)GO TO 19
C-----Transfer intermolecular data into common files (NS2)
      DO 18 I=1,NS4
      IF(DS4(I).GT.3.8)GO TO 18
      NS2=NS2+1
      NSS=NS2
      DO 15 L=1,6
      AT2(NS2,L)=ATM4(I,L)
15    AT22(NSS,L)=AT2(NS2,L)
      DS2(NS2)=DS4(I)
      DS22(NSS)=DS2(NS2)
      DO 16 K=1,3
      XE(NS2,K)=X4(I,K)
      XEE(NSS,K)=XE(NS2,K)
      XOE(NS2,K)=0.0
      DO 16 J=1,3
      XOE(NS2,K)=XOE(NS2,K)+O(K,J)*X4(I,J)
16    XOEE(NSS,K)=XOE(NS2,K)
      DO 17 L=1,38
      IF(SSTR(I,1).EQ.'+')SSTR(I,1)=' '
      IF(SSTR(I,L).EQ.'+'.AND.SSTR(I,L-1).EQ.',')SSTR(I,L)=' '
      IF(SSTR(I,L).EQ.'X')SSTR(I,L)='x'
      IF(SSTR(I,L).EQ.'Y')SSTR(I,L)='y'
      IF(SSTR(I,L).EQ.'Z')SSTR(I,L)='z'
      SST1(NS2,L)=SSTR(I,L)
17    SST2(NSS,L)=SST1(NS2,L)
18    CONTINUE
19    IF(KKD1.EQ.1)THEN
      WRITE(IO,'(''Atomic contacts, (NS2) = '',I3)')NS2
      DO 201 I=1,NS2
 201  WRITE(IO,'(6A1,''...'',6A1,1X,20A1,F9.3)')(AT(JO1,L),L=1,6),
     1 (AT2(I,L),L=1,6),(SST1(I,L),L=1,20),DS2(I)
      ENDIF
      CONTINUE
      IF(NS2.LT.2)THEN
      WRITE(IO,'('' There are no contacts suitable for determining '',
     1 ''the hydrogens coordinates'')')
      KKK=1
      RETURN
      ENDIF
C-----Look for forbidden directions
      DO 20 K=1,3
20    XAQ(K)=XO(JO1,K)
      IF(NS2.LT.3)GO TO 21
      CALL FRBDR(BRV,N,JO1,HD,NS2,AT2,DS2,XE,XOE,SST1)
      IF(KKD1.EQ.1)THEN
      WRITE(IO,'(''Number of contacts after elimination of '',
     1 ''forbidden directions:'',I3)')NS2
      DO 202 I=1,NS2
202   WRITE(IO,'(6A1,20A1,F7.3)')(AT2(I,L),L=1,6),(SST1(I,L),L=1,20),
     1 DS2(I)
      ENDIF
21    IF(NS2.LT.2)THEN
C-----The number of contacts is less than 2
      NS2=NSS
      DO 47 I=1,NS2
      DS2(I)=DS22(I)
      DO 46 K=1,3
      XE(I,K)=XEE(I,K)
46    XOE(I,K)=XOEE(I,K)
      DO 47 L=1,6
      AT2(I,L)=AT22(I,L)
47    SST1(I,L)=SST2(I,L)
      ENDIF
      IF(KKD1.EQ.1)THEN
      WRITE(IO,'(''Total number of contacts (NS2) ='',I2)')NS2
      DO 500 I=1,NS2
500   WRITE(IO,'(6A1,1X,20A1,F7.3)')(AT2(I,L),l=1,6),(SST1(I,L),L=1,20),
     1 DS2(I)
      ENDIF
C-----The number of contacts is just 2
      IF(NS2.EQ.2)GO TO 33
      IF(NS2.LT.3)GO TO 28
C-----Check the AA...OW...BB angles
      IF(KKD1.EQ.1)WRITE(IO,'(''Angles:'')')
      NAN=0
      ANDM=100
      DO 22 K=1,3
22    XD(2,K)=XO(JO1,K)
      DO 24 I=1,NS2-1
      DO 24 J=I+1,NS2
      DO 23 K=1,3
      XD(1,K)=XOE(I,K)
23    XD(3,K)=XOE(J,K)
      CALL BNDAN(RD,XD,T1)
      IF(KKD1.EQ.1)WRITE(IO,'(6A1,''-'',6A1,''-'',6A1,2X,F7.2)')
     1 (AT2(I,L),L=1,6),(AT(JO1,L),L=1,6),(AT2(J,L),L=1,6),T1
      DSA=DS2(I)
      DSB=DS2(J)
      IF(AT2(I,1).EQ.'C'.AND.AT2(I,2).EQ.'L')DSA=DS2(I)-0.400
      IF(AT2(J,1).EQ.'C'.AND.AT2(J,2).EQ.'L')DSB=DS2(J)-0.400
      IF(AT2(I,1).EQ.'B'.AND.AT2(I,2).EQ.'R')DSA=DS2(I)-0.600
      IF(AT2(J,1).EQ.'B'.AND.AT2(J,2).EQ.'R')DSB=DS2(J)-0.600
      IF(AT2(I,1).EQ.'I'.AND.(AT2(I,2).NE.'N'.OR.AT2(I,2).NE.'R'))
     1 DSA=DS2(I)-0.750
      IF(AT2(J,1).EQ.'I'.AND.(AT2(J,2).NE.'N'.OR.AT2(J,2).NE.'R'))
     1 DSA=DS2(I)-0.750
      DTA=ABS(T1-109.5)
      CDT=COS(DTA*0.0174532925)
      ANDX=(DSA+DSB)/CDT
      IF(KKD1.EQ.1)WRITE(IO,'(2F7.3,F6.2,2F10.4)')DSA,DSB,
     1 DTA,CDT,ANDX
C-----Choses the best angle (nearest to 109.5) with the shortest
C-----Distances
      IF(ABS(ANDX).LT.ABS(ANDM))THEN
      ANDM=ANDX
      IANDM(1)=I
      IANDM(2)=J
      ANM=T1
      ENDIF
C     IF(T1.LT.85.0.OR.T1.GT.125)GO TO 24
      NAN=NAN+1
      AN2(NAN)=T1
      IAN(NAN,1)=I
      IAN(NAN,2)=J
24    CONTINUE
      KI=0
      DO 49 I=1,NAN
      IF(AN2(I).LT.70.OR.AN2(I).GT.150)THEN
      KI=KI+1
      ENDIF
49    CONTINUE
      IF(KI.EQ.NAN)GO TO 28
C-----All the angles are in the 70-150 deg range
      IF(KKD1.EQ.1)WRITE(IO,'(''Best angle:''/2F7.3,F10.4,F8.2)')
     1 DS2(IANDM(1)),DS2(IANDM(2)),ANDM,ANM
      IF(KKD1.EQ.1)WRITE(IO,'(''Number of angles (NAN) = '',I2)')NAN
      IF(NAN.EQ.0)GO TO 28
      DO 27 I=1,2
      DO 26 K=1,3
      XA(I,K)=XE(IANDM(I),K)
26    XOA(I,K)=XOE(IANDM(I),K)
      DO 27 L=1,6
27    ATA(I,L)=AT2(IANDM(I),L)
C     IF(KKD1.EQ.1)WRITE(IO,'(6A1,3X,3F8.4)')((ATA(I,L),L=1,6),
C    1 (XA(I,K),K=1,3),I=1,2)
      GO TO 61
C-----Contacts are put in increasing order
28    DO 32 I=1,NS2
      DO 32 J=1,NS2-1
      IF(DS2(J).GT.DS2(J+1))THEN
      RS=DS2(J)
      DS2(J)=DS2(J+1)
      DS2(J+1)=RS
      DO 29 K=1,3
      RS=XE(J,K)
      XE(J,K)=XE(J+1,K)
      XE(J+1,K)=RS
      RS=XOE(J,K)
      XOE(J,K)=XOE(J+1,K)
29    XOE(J+1,K)=RS
      DO 30 L=1,6
      LR=AT2(J,L)
      AT2(J,L)=AT2(J+1,L)
30    AT2(J+1,L)=LR
      DO 31 L=1,20
      LR=SST1(J,L)
      SST1(J,L)=SST1(J+1,L)
31    SST1(J+1,L)=LR
      ENDIF
32    CONTINUE
C-----Omit repeated values
      NJ=1
      DO 50 I=1,NS2-1
      DO 50 J=I,NS2
      SX1=XE(I,1)+XE(I,2)+XE(I,3)
      SX2=XE(J,1)+XE(J,2)+XE(J,3)
      IF(ABS(SX1-SX2).GT.0.5)GO TO 50
      NJ=NJ+1
      NA(NJ)=J
50    Continue
      N1=0
      DO 55 I=1,NS2
      DO 51 J=1,NJ
      IF(I.EQ.NA(J))GO TO 55
51    CONTINUE
      N1=N1+1
      DB(N1)=DS2(I)
      DO 52 L=1,6
52    A1(N1,L)=AT2(I,L)
      DO 53 L=1,20
53    S1(N1,L)=SSTR(I,L)
      DO 54 K=1,3
      XB(N1,K)=XE(I,K)
54    XOB(N1,K)=XOE(I,K)
55    CONTINUE
      DO 59 I=1,N1
      DS2(I)=DB(I)
      DO 56 L=1,6
56    AT2(I,L)=A1(I,L)
      DO 57 L=1,20
57    SSTR(I,L)=S1(I,L)
      DO 58 K=1,3
      XE(I,K)=XB(I,K)
58    XOE(I,K)=XOB(I,K)
59    CONTINUE
      NS2=N1
      IF(NS2.LT.2.AND.KII.EQ.0)THEN
      KII=KII+1
      GO TO 21
      ENDIF
      IF(KKD1.EQ.1)THEN
      WRITE(IO,'(''Distances from '',6A1)')(AT(JO1,L),L=1,6)
      DO 203 I=1,NS2
 203  WRITE(IO,'(6A1,20A1,F7.3)')(AT2(I,L),L=1,6),(SST1(I,L),L=1,20),
     1 DS2(I)
      ENDIF
      IF(NS2.LT.2)THEN
      WRITE(IO,'(''It is not possible to locate the hydrogens of '',
     1 ''water '',6A1/''Change your input data omitting it'')')
     2 (AT(JO1,L),L=1,6)
      STOP
      ENDIF
C-----The two shortest contacts are assumed
33    DO 35 I=1,2
      DO 34 K=1,3
      XA(I,K)=XE(I,K)
34    XOA(I,K)=XOE(I,K)
      DO 35 L=1,6
35    ATA(I,L)=AT2(I,L)
61    IF(KKD1.EQ.1)THEN
      WRITE(IO,'(2(6A1,3X,3F8.4/))')(ATA(1,L),L=1,6),(XA(1,K),K=1,3),
     1 (ATA(2,L),L=1,6),(XA(2,K),K=1,3)
      ENDIF
C-----Calculates the coordinates of the two water hydrogens on the
C-----lines of the best contacts
36    ATTH(1,1)='H'
      ATTH(2,1)='H'
      ATTH(1,2)='1'
      ATTH(2,2)='2'
      DO 40 I=1,2
      DT=0.0
      DO 37 K=1,3
      XOO(1,K)=XO(JO1,K)
37    XOO(2,K)=XOA(I,K)
      CALL HCOO(HD,XOO,XH,XOH)
      DO 38 K=1,3
      XHN(I,K)=XH(K)
      XON(I,K)=XOH(K)
38    DT=DT+(XO(JO1,K)-XON(I,K))**2
      DT=SQRT(DT)
      DO 39 L=3,6
39    ATTH(I,L)=AT(JO1,L-2)
40    CONTINUE
C-----Print distances and angle of the water molecule
      D1=0.0
      D2=0.0
      DO 41 K=1,3
      XD(1,K)=XON(1,K)
      XD(2,K)=XO(JO1,K)
      XD(3,K)=XON(2,K)
      D1=D1+(XON(1,K)-XO(JO1,K))**2
41    D2=D2+(XON(2,K)-XO(JO1,K))**2
      D1=SQRT(D1)
      D2=SQRT(D2)
      CALL BNDAN(RD,XD,T1)
      DO 42 K=1,3
      XD(1,K)=XOA(1,K)
      XD(2,K)=XO(JO1,K)
42    XD(3,K)=XOA(2,K)
      CALL BNDAN(RD,XD,T1)
      DO 43 L=1,6
      ATW(1,L)=ATTH(1,L)
      ATW(2,L)=AT(JO1,L)
43    ATW(3,L)=ATTH(2,L)
      DO 44 K=1,3
      XOW(1,K)=XON(1,K)
      XOW(2,K)=XO(JO1,K)
44    XOW(3,K)=XON(2,K)
      CALL WATANG(RD,HD,ATW,XOW,XHN,XON)
      RETURN
      END
 
      SUBROUTINE TORS(RD,R,TAU)
C-----This routine calculates the torsion angles formed by the atoms
C-----1-2-3-4.The angle is positive when the 1-2 bond, view down the
C-----2-3 bond, will eclipse the 3-4 bond when rotates less than 180
C-----degrees in a clockwise direction ("right-hand rule",Klyne,W.&
C-----Prelog,V. (1960), Experientia, 16, 521).
      REAL M(3,3)
      DIMENSION R(4,3),DLT(3,3),DST(3)
      DO 1 J=1,3
      DLT(1,J)=R(2,J)-R(1,J)
      DLT(2,J)=R(2,J)-R(3,J)
1     DLT(3,J)=R(4,J)-R(3,J)
      DO 2 I=1,3
      DST(I)=SQRT(DLT(I,1)**2+DLT(I,2)**2+DLT(I,3)**2)
      IF(DST(I).EQ.0.) GO TO 4
2     CONTINUE
      DO 3 I=1,3
      DO 3 J=1,3
3     M(I,J)=DLT(I,J)/DST(I)
      CF1=M(1,1)*M(2,1)+M(1,2)*M(2,2)+M(1,3)*M(2,3)
      CF2=M(2,1)*M(3,1)+M(2,2)*M(3,2)+M(2,3)*M(3,3)
      IF(ABS(CF1).GT..98.OR.ABS(CF2).GT..98 )GO TO 4
      SF1=SQRT(1.-CF1**2)
      SF2=SQRT(1.-CF2**2)
      P=SF1*SF2
      IF(ABS(P).LT.1.E-5)GO TO 4
      STAU=(M(3,1)*(M(2,2)*M(1,3)-M(1,2)*M(2,3))
     1 +M(3,2)*(M(1,1)*M(2,3)-M(1,3)*M(2,1))
     2 +M(3,3)*(M(2,1)*M(1,2)-M(1,1)*M(2,2)))/P
      CTAU=((M(2,2)*M(1,3)-M(1,2)*M(2,3))*( M(3,2)*M(2,3)-M(2,2)*M(3,3)
     1 )+(M(2,3)*M(1,1)-M(2,1)*M(1,3))*(M(3,3)*M(2,1)-M(3,1)*M(2,3))
     2 +(M(2,1)*M(1,2)-M(1,1)*M(2,2))*(M(3,1)*M(2,2)-M(2,1)*M(3,2)))/P
      TAU1=ATAN2(STAU,CTAU)
      TAU=RD*TAU1
      IF(TAU.GT.180.)TAU=TAU-360.
4     RETURN
      END
 
      SUBROUTINE TRSL(KY,NNS,LA1,ATT,ATM4,NS4,X4,DS4,SSTR)
C-----Calculate translation coordinates
      CHARACTER*1 ATT(6),ATM4(20,6),SSTR(20,38),LA1(2,20)
      INTEGER LE1(24,3,3),FT(3),LE(3,3)
      DIMENSION F1(24,3),FB(3),F(3),X4(50,3),DS4(50)
      COMMON/COOA/X1(3),X2(3)
     1      /LBL3/IC,NE,NT
     2      /GRP1/C(3),A(3)
     3      /INS4/CS(3),SN(3)
     4      /EQVP/LE1,F1,FB
     6      /INPOU/IN,IO,IR
      KM=2*NT+1
      DO 3 I=1,NE
      DO 3 I1=1,KM
      FT(1)=I1-NT-1
      DO 3 I2=1,KM
      FT(2)=I2-NT-1
      DO 3 I3=1,KM
      FT(3)=I3-NT-1
      DO 1 J=1,3
      F(J)=F1(I,J)+FB(J)+FT(J)
      DO 1 K=1,3
1     LE(J,K)=LE1(I,J,K)
      CALL DSTN(KY,NNS,LA1,LE,F,NS4,ATT,ATM4,X4,DS4,SSTR)
      IF(IC.NE.-1) GO TO 3
      DO 2 J=1,3
      F(J)=-F1(I,J)+FB(J)+FT(J)
      DO 2 K=1,3
2     LE(J,K)=-LE1(I,J,K)
      CALL DSTN(KY,NNS,LA1,LE,F,NS4,ATT,ATM4,X4,DS4,SSTR)
3     CONTINUE
      RETURN
      END
 
      SUBROUTINE WATANG(RD,HD,ATW,XOW,XH,XOH)
C-----Calculates the cooordinates of the two water hydrogens in such
C-----a way as to have the two O-H bonds forming an angle of 107.7 deg
C-----without changing the plane of the molecule
      CHARACTER*1 AT(4,6),ATW(3,6)
      DIMENSION XH(2,3),XOH(2,3),XD(3,3),DDC(3),XO(4,3),XOW(3,3),X(3)
      COMMON/NEW1/O(3,3),OI(3,3)
     1      /INPOU/IN,IO,IR
      DO 1 I=1,3
      DO 1 L=1,6
1     AT(I,L)=ATW(I,L)
      DO 2 I=1,3
      DO 2 K=1,3
2     XO(I,K)=XOW(I,K)
      DO 3 K=1,3
3     XO(4,K)=(XO(1,K)+XO(3,K))/2
      AT(4,1)='H'
      AT(4,2)='-'
      AT(4,3)='H'
      AT(4,4)='M'
      AT(4,5)='I'
      AT(4,6)='D'
      DO 4 I=1,3
      DO 4 K=1,3
4     XD(I,K)=XO(I,K)
      DD1=0
      DD2=0
      DO 5 K=1,3
      DD1=DD1+(XO(1,K)-XO(2,K))**2
5     DD2=DD2+(XO(3,K)-XO(2,K))**2
      DD1=SQRT(DD1)
      DD2=SQRT(DD2)
      CALL BNDAN(RD,XD,T)
      D14=0.0
      D24=0.0
      DO 6 K=1,3
      DDC(K)=XO(1,K)-XO(4,K)
      D14=D14+DDC(K)**2
6     D24=D24+(XO(2,K)-XO(4,K))**2
      D14=SQRT(D14)
      D24=SQRT(D24)
      DDD=D24*TAN(53.85/RD)
      DLT=D14-DDD
      DO 7 K=1,3
      DDC(K)=DDC(K)/D14
7     XOH(1,K)=XO(1,K)-DDC(K)*DLT
      D14=0.0
      DO 8 K=1,3
      DDC(K)=XO(3,K)-XO(4,K)
8     D14=D14+DDC(K)**2
      D14=SQRT(D14)
      DLT=D14-DDD
      DO 9 K=1,3
      DDC(K)=DDC(K)/D14
9     XOH(2,K)=XO(3,K)-DDC(K)*DLT
      D15=0.0
      DO 10 K=1,3
      DDC(K)=XO(2,K)-XOH(1,K)
10    D15=D15+DDC(K)**2
      D15=SQRT(D15)
      DL1=D15-HD
      DO 11 K=1,3
      DDC(K)=DDC(K)/D15
11    XOH(1,K)=XOH(1,K)+DDC(K)*DL1
      D15=0.0
      DO 12 K=1,3
      DDC(K)=XO(2,K)-XOH(2,K)
12    D15=D15+DDC(K)**2
      D15=SQRT(D15)
      DL1=D15-HD
      DO 13 K=1,3
      DDC(K)=DDC(K)/D15
13    XOH(2,K)=XOH(2,K)+DDC(K)*DL1
      DO 14 I=1,2
      DO 14 J=1,3
      XH(I,J)=0.0
      DO 14 K=1,3
14    XH(I,J)=XH(I,J)+OI(J,K)*XOH(I,K)
      DO 15 K=1,3
      XD(1,K)=XOH(1,K)
      XD(2,K)=XO(2,K)
15    XD(3,K)=XOH(2,K)
      DO 17 I=1,3,2
      DO 16 L=1,6
      IF(AT(I,L).EQ.' ')THEN
      AT(I,L)='"'
      IF(I.EQ.1)GO TO 17
      IF(I.EQ.3)GO TO 18
      ENDIF
16    CONTINUE
17    CONTINUE
18    CALL BNDAN(RD,XD,T)
      DD1=0
      DD2=0
      DO 19 K=1,3
      X(K)=0.0
      DD1=DD1+(XOH(1,K)-XO(2,K))**2
19    DD2=DD2+(XOH(2,K)-XO(2,K))**2
      DD1=SQRT(DD1)
      DD2=SQRT(DD2)
      CALL BNDAN(RD,XD,T)
      RETURN
      END
 
      BLOCK DATA
      CHARACTER*1 MINU(26),MAIU(26)
      CHARACTER*2 NOMAT(103)
      DIMENSION IRAG(103)
      COMMON/SMBL/NOMAT
     1      /ATRD/IRAG
     2      /MIMA/MINU,MAIU
      DATA MINU/'a','b','c','d','e','f','g','h','i','j','k','l','m',
     1 'n','o','p','q','r','s','t','u','v','w','x','y','z'/
      DATA MAIU/'A','B','C','D','E','F','G','H','I','J','K','L','M',
     1 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
      DATA NOMAT/'H ','HE','LI','BE','B ','C ','N ','O ','F ','NE',
     1 'NA','MG','AL','SI','P ','S ','CL','AR','K ','CA','SC',
     2 'TI','V ','CR','MN','FE','CO','NI','CU','ZN','GA','GE',
     3 'AS','SE','BR','KR','RB','SR','Y ','ZR','NB','MO','TC',
     4 'RU','RH','PD','AG','CD','IN','SN','SB','TE','I ','XE',
     5 'CS','BA','LA','CE','PR','ND','PM','SM','EU','GD','TB',
     6 'DY','HO','ER','TM','YB','LU','HF','TA','W ','RE','OS',
     7 'IR','PT','AU','HG','TL','PB','BI','PO','AT','RN','FR',
     8 'RA','AC','TH','PA','U ','NP','PU','AM','CM','BK','CF',
     9 'ES','FM','MD','NO','LW'/
      DATA IRAG/37,93,134,90,82,77,75,73,72,131,154,130,118,
     1 111,106,102,99,174,196,174,144,136,117,128,122,130,126,
     2 140,138,131,126,122,119,116,114,189,211,192,162,148,
     3 139,133,135,140,145,148,153,148,144,141,138,135,133,
     4 209,225,198,169,165,165,164,165,166,185,161,159,159,
     5 158,157,156,170,156,157,143,137,128,126,132,150,150,
     6 149,148,147,146,168,150,214,230,204,188,165,150,142,
     7 165,166,185,161,160,160,160,160,160,160,160/
      END
 
      SUBROUTINE FRBDR(BRV,N,JO1,HD,NS2,AT2,DS2,XE,XOE,SSTR)
C-----Omit atoms too near (<3.0 A) to a water hydrogen
      CHARACTER*1 AT2(20,6),SSTR(20,38),AT2D(50,6),AT(500,6),
     1 A1(20,6),S1(20,20),BRV
      DIMENSION DS2(50),XE(40,3),XOE(40,3),XOO(3,3),XH(3),XOH(3),
     1 FB(3),X1(3),X2(3),X3(3),ID(20),D1(50),XA(40,3),XOA(40,3),
     2 B(3,3)
      COMMON/INPOU/IN,IO,IR
     1      /LTT/KB,B
     2      /COOR/X(500,3),Q(500)
     3      /GRP3/XO(500,3)
     4      /GRP4/AT
      CALL EQPBRV(BRV,KB,B)
      NO=0
      DO 1 K=1,3
      FB(K)=0.0
1     XOO(1,K)=XO(JO1,K)
      DO 12 I=1,NS2
      DO 2 K=1,3
      X3(K)=XE(I,K)
2     XOO(2,K)=XOE(I,K)
      CALL HCOO(HD,XOO,XH,XOH)
      DO 11 J=1,N
      IF(J.EQ.JO1)GO TO 11
      DO 3 L=1,6
      IF(AT(J,L).NE.AT2(I,L))GO TO 4
3     CONTINUE
      GO TO 11
4     DT=0.0
C-----Compare intramolecular  H...A contacts
      DO 5 K=1,3
      X1(K)=XH(K)
      X2(K)=X(J,K)
5     DT=DT+(XOH(K)-XO(J,K))**2
      DT=SQRT(DT)
      IF(DT.LT.2.4)THEN
      NO=NO+1
      IF(NO.GT.49)RETURN
      ID(NO)=I
      DO 6 L=1,6
6     AT2D(NO,L)=AT2(I,L)
      GO TO 12
      ENDIF
C-----Compare intermolecular  H...A contacts
      ICT=0
      CALL TRS1(X1,X2,X3,ICT)
      IF(ICT.EQ.1)THEN
      NO=NO+1
      IF(NO.GT.49)RETURN
      ID(NO)=I
      DO 7 L=1,6
7     AT2D(NO,L)=AT2(I,L)
      GO TO 12
      ENDIF
      IF(KB.EQ.0)GO TO 11
      DO 10 L=1,KB
      DO 8 K=1,3
8     FB(K)=B(L,K)
      CALL TRS1(X1,X2,X3,ICT)
      IF(ICT.EQ.1)THEN
      NO=NO+1
      IF(NO.GT.49)RETURN
      ID(NO)=I
      DO 9 LL=1,6
9     AT2D(NO,LL)=AT2(I,LL)
      GO TO 12
      ENDIF
10    CONTINUE
11    CONTINUE
12    CONTINUE
      N1=0
      DO 17 I=1,NS2
      DO 13 J=1,NO
      IF(I.EQ.ID(J))GO TO 17
13    CONTINUE
      N1=N1+1
      D1(N1)=DS2(I)
      DO 14 L=1,6
14    A1(N1,L)=AT2(I,L)
      DO 15 L=1,20
15    S1(N1,L)=SSTR(I,L)
      DO 16 K=1,3
      XA(N1,K)=XE(I,K)
16    XOA(N1,K)=XOE(I,K)
17    CONTINUE
      DO 21 I=1,N1
      DS2(I)=D1(I)
      DO 18 L=1,6
18    AT2(I,L)=A1(I,L)
      DO 19 L=1,20
19    SSTR(I,L)=S1(I,L)
      DO 20 K=1,3
      XE(I,K)=XA(I,K)
20    XOE(I,K)=XOA(I,K)
21    CONTINUE
      NS2=N1
      RETURN
      END
 
 
      SUBROUTINE TRS1(X1,X2,X3,ICT)
C-----Set data and calculates contacts involving atoms translated in
C-----equivalent positions
      CHARACTER*1 STPE(38,20)
      INTEGER LE1(24,3,3),FT(3),LE(3,3)
      DIMENSION F1(24,3),FB(3),F(3),X1(3),X2(3),X3(3)
      COMMON/LBL3/IC,NE,NT
     1      /GRP1/C(3),A(3)
     2      /INS4/CS(3),SN(3)
     3      /EQVP/LE1,F1,FB
     4      /LBL4/STPE
     5      /INPOU/IN,IO,IR
      KM=2*NT+1
      DO 3 I=1,NE
      DO 3 I1=1,KM
      FT(1)=I1-NT-1
      DO 3 I2=1,KM
      FT(2)=I2-NT-1
      DO 3 I3=1,KM
      FT(3)=I3-NT-1
      DO 1 J=1,3
      F(J)=F1(I,J)+FB(J)+FT(J)
      DO 1 K=1,3
1     LE(J,K)=LE1(I,J,K)
      CALL DST1(X1,X2,X3,LE,F,ICT)
      IF(IC.NE.-1) GO TO 3
      DO 2 J=1,3
      F(J)=-F1(I,J)+FB(J)+FT(J)
      DO 2 K=1,3
2     LE(J,K)=-LE1(I,J,K)
      CALL DST1(X1,X2,X3,LE,F,ICT)
3     CONTINUE
      RETURN
      END
 
      SUBROUTINE DST1(X1,X2,X3,LE,F,ICT)
C-----This routine calculates interatomic distances involving
C-----atoms translated in equivalent positions
      INTEGER LE(3,3)
      DIMENSION F(3),X(3),D(3),X1(3),X2(3),X3(3)
      COMMON/LBL3/IC,NE,NT
     1      /GRP1/C(3),A(3)
     2      /INS4/CS(3),SN(3)
     3      /INPOU/IN,IO,IR
      DO 1 I=1,3
      X(I)=F(I)
      DO 1 J=1,3
1     X(I)=X(I)+LE(I,J)*X2(J)
      DO 2 I=1,3
2     IF(ABS(X(I)-X2(I)).GT.1E-4)GO TO 3
      RETURN
3     DO 5 K=1,3
      IF((X(K)-X3(K)).GT.1E-4)GO TO 6
5     CONTINUE
      RETURN
6     DA=0
      DO 4 I=1,3
      D(I)=X1(I)-X(I)
      D(I)=D(I)*C(I)
4     DA=DA+D(I)**2
      IF(DA.LT.1E-2)RETURN
      DA=DA+2*(D(1)*D(2)*CS(3)+D(1)*D(3)*CS(2)+D(2)*D(3)*CS(1))
      IF(DA.GT.1E-2)THEN
      DA=SQRT(DA)
      ELSE
      RETURN
      ENDIF
      IF(DA.LT.2.4)THEN
      ICT=1
      ENDIF
      RETURN
      END
