      Program ROTENER
C**********************************************************************
C* Release September 1996                                             *
C*                                                                    *
C*                      Program ROTENER                               *
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* Literature: J.Chem.Soc.Dalton, 1990, p.179; ibid. 1991, p.203      *
C*             Gazz. Chim. Ital., (1991), 121, 433-439                *
C*             QCPE BULLETIN. (1991), 11, 61                          *
C*             Bull. Chem. Soc. Jpn.(1994). 67, 189-195               *
C*             J.Chem.Inform. Computer Sci. (1994). 34, 972-975       *
C*             Acta Cryst.(1994). C50, 1323-1326                      *
C*                                                                    *
C* This program calculates the difference molecular potential energy  *
C* after rotations of one or two molecular fragments about given      *
C* directions. The first two atoms of each set to be rotated define   *
C* the directions about which the rotations occur.                    *
C* Two options are available: isolated molecule or molecule in the    *
C* crystal. In the latter case the crystal environment is considered  *
C* static.                                                            *
C* The van der Waals energy is calculated only if the atoms involved  *
C* in the contacts are: H, C, N, O, F, Si, P, S, Cl, Fe, Br, Sn, I,   *
C* Pb.                                                                *
C* The Coulombic energy is considered only if the atomic partial      *
C* charges are given with non zero values.                            *
C* The options isolated molecule and molecule in the crystal (static  *
C* environment are available.                                         *
C*                                                                    *
C*                         INPUT                                      *
C* 1. TITLE (A80)                                                     *
C* 2. Bravais symbol (A1)                                             *
C* 3. KY code for atomic species. KY=0 if not given, KY=1 if given    *
C* 4. a, b, c, alpha, beta, gamma (free FORMAT)                       *
C* 5. NNS = number of atomic species (free FORMAT)                    *
C* 6. Atomic species between apices (free FORMAT)                     *
C* 7. Format of atomic coordinates (fractional) and charges           *
C*    (electrons) [e.g.(1X,6A1,1X,3F9.5,F7.3)]                        *
C* 8. N (total number of atoms)                                       *
C* 9. Atom, x, y, z, q (6A1,3F8.5,F7.3) (max number of atoms: 400)    *
C*10. The number of the atoms of the first rotating fragment,theta    *
C*    minimum (rotation angle, degr), theta maximum, delta theta      *
C*    (step of theta)(free format).                                   *
C*11. The labels of four atoms defining a torsion angle about the     *
C*    first direction of rotation (4(6A1))                            *
C*12. The labels of the atoms of the first rotating fragment. The     *
C*    origin of the reference orthogonal system of axes is on the     *
C*    first atom and the rotation occurs counterclock-wise about the  *
C*    direction defined by the first two atoms (12(6A1)/12(6A1))      *
C*  N.B. If the number of atoms is EQUAL to 12 or 24..., put a blank  *
C*       record after their list                                      *          *
C*13. The number of the atoms of the second rotating fragment,theta   *
C*    minimum (rotation angle, degr), theta maximum, delta theta      *
C*    (step of theta) (free format).                                  *
C*14. The labels of four atoms defining a torsion angle about the     *
C*    second direction of rotation (4(6A1)).                          *
C*15. The labels of the atoms of the second rotating fragment. The    *
C*    rotation occurs counterclock-wise about the direction defined   *
C*    by the first two atoms (12(6A1)/12(6A1)).                       *
C*16. Data for environment in the crystal (not given if isolated      *
C*    molecule is considered: IC,NE,NT (Free Format)                  *
C*    IC=1 if the space group is acentric                             *
C*      =-1 if the space group is centric                             *
C*    NE=Number of the equivalent positions (the general, x y z, the  *
C*       centrosymmetric and those generated by the Bravais lattice   *
C*       must be omitted                                              *
C*    NT=Maximum translation required (suggested 1 when the set of    *
C*       atoms in the general position is the nearest to the origin,  *
C*       2 in the other cases)                                        *
C*17. Equivalent position: one equivalent position per record (20A1)  *
C*    Examples:  1/2-X,1/2+Y,1/2+Z  (not 0.5-X, etc...)               *
C*               Y-X,-X,1/3+Z       (not 0.3333+Z)                    *
C*                                                                    *
C* N.B.- If the second fragment is part of the firts one the names of *
C*    the atoms of the latter fragment must be present also among     *
C*    those for the first one.                                        *
C*                                                                    *
C*       The program interactively ask for the number of fragments to *
C*    be considered and for isolated molecule or in the crystal.      *
C*    For rotations of 360 degr. and one fragment, steps of 5 degr.   *
C*    are suggested, while with two fragments, steps of 10 degr. are  *
C*    advisable. Of course the computing time depends also upon the   *
C*    number of atoms.                                                *
C*                                                                    *
C*                         OUTPUT                                     *
C*                                                                    *
C* In user defined files:                                             *
C*                                                                    *
C* Title                                                              *
C*                                                                    *
C* On request: crystal data and atomic coordinates.                   *
C*                                                                    *
C*                                                                    *
C* In the case of one rotating fragment:                              *
C* =====================================                              *
C*                                                                    *
C* A table giving the rotation angles, the torsion angles, the        *
C* difference energy with respect the unrotated molecule, in kcal/mol *
C* and kJ/mol. In addition, interatomic contacts less than a given    *
C* value (requested interactively; suggested 2-2.5 A) and energies    *
C* corresponding to each rotation of the fragment are also quoted.    *
C*                                                                    *
C* On request: tables of data for two-dimensional graphic routines,   *
C* the names of the files are given by the user.                      *
C*                                                                    *
C*                                                                    *
C* In the case of two fragments:                                      *
C* ============================                                       *
C*                                                                    *
C* A table giving the difference energy (with respect the unrotated   *
C* molecule) for each couple of rotation angles.                      *
C*                                                                    *
C* A line-printer plot of the difference energy values as function of *
C* the rotation angles.                                               *
C*                                                                    *
C* On request: a table giving the T1, T2, DE data suitable as input   *
C* for a tridimensional plotting routine.                             *
C*                                                                    *
C*                       -----------------                            *
C* WARNING                                                            *
C*   GETTIM and GETDAT are intrinsec routines of the IBM Professional *
C*   FORTRAN Compiler PROFORT Version 1.0 [(C) Copyright IBM Corp     *
C*   1984, (C) Copyright Ryan-McFarland Corp 1984] library used to    *
C*   get time and date. If a different compiler is used, the          *
C*   following sentences of the Main must be changed or deleted.      *
C*  At the beginning:                                                 *
C*    C-----Read the initial time                                     *
C*          CALL GETTIM(IHR,IMN,ISC,IHS)                              *
C*          HS=IHS                                                    *
C*          TM0=(IHR*3600)+(IMN*60)+ISC+HS/100.0                      *
C*  At the end:                                                       *
C*    C-----Print computing time and date                             *
C*          CALL GETDAT(IYR,IMT,IDY)                                  *
C*          CALL GETTIM(IHR,IMN,ISC,IHS)                              *
C*          HS=IHS                                                    *
C*          TM=(IHR*3600)+(IMN*60)+ISC+HS/100.0                       *
C*          IF(TM.LT.TM0)TM=TM+86400                                  *
C*          DTM=TM-TM0                                                *
C*          WRITE(IO,'(/'' Total computing time: '',F9.2,'' s'')')DTM *
C*          WRITE(IO,140)                                             *
C*          WRITE(IO,140)IDY,MONTH(IMT),IYR,IHR,IMN,ISC,IHS           *
C*    140   FORMAT(I4,'/',A4,'/',I4,3X,I4,':',I2,':',I2,'.',I2/)      *
C*                                                                    *
C**********************************************************************
      CHARACTER*1 TIT(80),FRT(80),AT(500,6),ATR1(500,6),MINUS(2),LR,
     1 AT1(4,6),AT2(4,6),ATR2(500,6),LP,QS,QL(2),BRV,AAA1(500,6),
     2 AAA2(500,6),AFR1(500,6),AFR2(500,6),STPE(38,20),SDL1(2,500),
     3 LA1(2,20),AAAA(6),LBA(2),CTEMP(6)
      CHARACTER*14 FILGRF,FILVIS,FILSRF
      CHARACTER*4 MONTH(12)
      CHARACTER*2 SDL2(500),LA(20)
      INTEGER*2 IHR,IMN,ISC,IHS,IYR,IMT,IDY
      DIMENSION XRO1(500,3),XRO2(500,3),XB(500,3),QB(500),
     1 R1(4,3),TAU1(180),TAU2(180),ENRG(72,72),IO1(72),IO2(72),
     2 IENRG(72,72),XYZ(5184,3),R2(4,3),RR1(3),RR2(3),XO1(500,3),
     3 XO2(500,3),XO3(500,3),XOO1(500,3),XOO2(500,3),X1K(3),
     4 TH(500),DEJ1(500),NT1(4),NT2(4),XF1(500,3),XXO(500,3),
     5 XBO(500,3),NA4(500),QF(500)
      DOUBLE PRECISION BI1(3,3),BI2(3,3)
      EQUIVALENCE (LA(1),LA1(1,1)),(SDL2(1),SDL1(1,1))
      COMMON/GRP1/O(3,3)
     1      /GRP2/AT
     2      /GRP3/MC(400,8)
     3      /GRP5/XO(500,3)
     4      /GRP6/NA1(500),NA2(500),NA3(500)
     5      /GRP7/Q1(500),Q2(500),Q3(500)
     6      /GRP8/DMAX
     7      /GRP10/X(500,3),Q(500)
     8      /LBL1/C(3),A(3)
     9      /LBL3/IC,NE,NT
      COMMON/LBL4/STPE
     2      /INPO/IN,IO
     3      /GRG1/SDL2
     4      /GRG2/LA
      DATA QL/'N','n'/
      DATA MONTH/'Jan.','Feb.','Mar.','Apr.','May.','Jun.','Jul.',
     1 'Aug.','Sep.','Oct.','Nov.','Dec.'/
C-----Read the initial time
      CALL GETTIM(IHR,IMN,ISC,IHS)
      HS=IHS
      TM0=(IHR*3600)+(IMN*60)+ISC+HS/100.0
C-----Constants and I/O files
      RD=0.0174532925
      IN=13
      IO=15
      CALL INPOUT
C-----Read title, codes and parameters
      READ(IN,100)TIT
100   FORMAT(80A1)
      WRITE(IO,101)TIT
101   FORMAT(1X,80A1)
      READ(IN,'(A1)')BRV
      READ(IN,*)KY
      READ(IN,*)C,A
      CALL ORTHOG(C,A)
      DO 60 I=1,20
60    LA(I)='  '
      IF(KY.EQ.0)GO TO 59
      READ(IN,*)NNS
      READ(IN,*)(LA(I),I=1,NNS)
      DO 57 I=1,NNS
      MINUS(1)=LA1(1,I)
      MINUS(2)=LA1(2,I)
      CALL MINUMA(MINUS)
      LA1(1,I)=MINUS(1)
57    LA1(2,I)=MINUS(2)
59    READ(IN,100) FRT
      WRITE(*,102)
102   FORMAT(' Do you want crystal data and atomic coordinates',
     1 ' printed? (Y/N)'/' ?>')
      READ(*,'(A1)')LP
      IF(LP.EQ.QL(1).OR.LP.EQ.QL(2))GO TO 1
      WRITE(IO,103)(C(I),A(I),I=1,3)
103   FORMAT(' Crystal data'/' a = ',F7.4,25X,'alpha= ',F6.2/
     1 ' b = ',F7.4,25X,'beta = ',F6.2/' c = ',F7.4,25X,'gamma= ',F6.2)
      WRITE(IO,104)((O(I,J),J=1,3),I=1,3)
104   FORMAT(/' Orthogonalization matrix'/3(1X,3F10.5/))
      WRITE(IO,105)
105   FORMAT(' Atomic coordinates'//' Atom',6X,'X/a',6X,'Y/b',
     1 6X,'Z/c',12X,'X',8X,'Y',8X,'Z'/)
1     READ(IN,*)N
      DO 4 I=1,N
      READ(IN,FRT)(AT(I,L),L=1,6),(X(I,J),J=1,3),Q(I)
      DO 2 K=1,5,2
      MINUS(1)=AT(I,K)
      MINUS(2)=AT(I,K+1)
      CALL MINUMA(MINUS)
      AT(I,K)=MINUS(1)
2     AT(I,K+1)=MINUS(2)
C-----Orthogonal coordinates
      DO 3 J=1,3
      XO(I,J)=0.0
      DO 3 K=1,3
3     XO(I,J)=XO(I,J)+O(J,K)*X(I,K)
C-----Matrix of the atom symbols
      DO 58 L=1,6
58    AAAA(L)=AT(I,L)
      CALL LABL(AAAA,KY,NNS,LA1,LBA,KKD)
      SDL1(1,I)=LBA(1)
      SDL1(2,I)=LBA(2)
      IF(LP.EQ.QL(1).OR.LP.EQ.QL(2))GO TO 4
      WRITE(IO,106)(AT(I,L),L=1,6),(X(I,J),J=1,3),(XO(I,K),K=1,3)
106   FORMAT(1X,6A1,2X,3F9.5,5X,3F9.4)
4     CONTINUE
C-----Set up the connectivity matrix MC(400,8)
      CALL BONDI(N)
C-----Read in the rotation parameters
      WRITE(*,107)
107   FORMAT(' Key in the number of fragments to be rotated (1 or 2)'/
     1 ' ?>')
      READ(*,*)KC
      WRITE(*,108)
108   FORMAT(' Key in 0 if the molecule is considered isolated '/
     1       '        1 if the molecule is considered in the crystal'/
     2 ' ?>')
      READ(*,*)KDE
      CALL READAT(NG1,T1,T3,DT1,AT1,ATR1)
      IF(KC.EQ.2) CALL READAT(NG2,T2,T4,DT2,AT2,ATR2)
C-----Consider crystal environment
      IF(KDE.EQ.1)THEN
      READ(IN,*)IC,NE,NT
      DO 5 I=1,NE
5     READ(IN,'(20A1)')(STPE(I,L),L=1,20)
      NEE=NE
C-----Environment of the first rotating fragment
      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)
C-----Environment of the second rotating fragment
      IF(KC.EQ.2)THEN
      NR2=NG2-2
      DO 7 I=3,NG2
      DO 7 L=1,6
7     AFR2(I-2,L)=ATR2(I,L)
      NE=NEE
      CALL CONTACT(BRV,N,NR2,AFR2,NB,AAA2,XBO,QB)
      ENDIF
      ENDIF
      T1=T1*RD
      T3=T3*RD
      DT1=DT1*RD
      IF(KC.EQ.2)THEN
      T2=T2*RD
      T20=T2
      T4=T4*RD
      DT2=DT2*RD
      ENDIF
C-----First rotating fragment
      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
      IF(KC.EQ.2)THEN
C-----Second rotating fragment
      DO 15 I=1,NG2
      DO 14 J=1,N
      DO 12 L=1,6
      IF(ATR2(I,L).NE.AT(J,L)) GO TO 14
12    CONTINUE
      DO 13 K=1,3
13    XO2(I,K)=XO(J,K)
      NA2(I)=J
      Q2(I)=Q(J)
      GO TO 15
14    CONTINUE
      WRITE(*,109)(ATR2(I,L),L=1,6)
      STOP
15    CONTINUE
      ENDIF
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
      IF(KC.EQ.2)THEN
      DO 17 I=3,NG2
      IF(NA2(I).EQ.J)GO TO 19
17    CONTINUE
      ENDIF
      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-----Atoms involved in the torsion angles
      CALL ATORS(N,AT1,NG1,NT1,NA1,R1)
      IF(KC.EQ.2)CALL ATORS(N,AT2,NG2,NT2,NA2,R2)
C-----Add crystal environment to the fix fragment
      IF(KDE.EQ.1)THEN
      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)
      IF(KC.EQ.2)THEN
      NI=0
      DO 23 I=1,NB
      DO 23 J=1,NL
      DO 20 L=1,3
      IF(ABS(XBO(I,L)-XXO(J,L)).GT.0.0001)GO TO 23
20    CONTINUE
      NI=NI+1
      DO 22 K=I,NB-NI
      DO 21 L=1,3
      XBO(K,L)=XBO((K+1),L)
      QB(K)=QB(K+1)
21    AAA2(K,L)=AAA2((K+1),L)
      DO 22 L=4,6
22    AAA2(K,L)=AAA2((K+1),L)
23    CONTINUE
      NB=NB-NI
      CALL ADD(N,NB,AAA2,XO3,XBO,Q3,QB,NG3)
C-----Add environmental atoms to the atom symbol matrix
      IF(N.GT.500) WRITE(*,'(/'' WARNING''/'' The total number '',
     1 ''of atoms is greater than 500''/)')
      DO 64 I=1,N
      DO 63 L=1,6
63    CTEMP(L)=AT(I,L)
      CALL LABL(CTEMP,KY,NNS,LA1,LBA,KKD)
      SDL1(1,I)=LBA(1)
64    SDL1(2,I)=LBA(2)
      ENDIF
      ENDIF
C-----Calculation of the molecular energy for the unrotated molecule
      KD=0
      CALL ENERMOL(KC,KD,NG1,NG2,NG3,XO1,XO2,XO3,E00)
      KD=1
      CALL ENERMOL(KC,KD,NG1,NG2,NG3,XO1,XO2,XO3,E01)
      IF(KC.EQ.2)THEN
      KD=2
      CALL ENERMOL(KC,KD,NG1,NG2,NG3,XO1,XO2,XO3,E02)
      E0=E00+E01+E02
      ELSE
      E0=E00+E01
      ENDIF
      E0J=E0*4.184
      IF(KDE.EQ.0)WRITE(IO,110)E0,E0J
110   FORMAT(/' Non-bonded energy for the isolated molecule:'/
     1 ' E0  = ',F10.2,' kcal/mol'/' E0  = ',F10.2,' kJ/mol')
C-----Ask for a upper limit value of delta(E)
      WRITE(*,111)
111   FORMAT(' Key in the maximum value of delta(E), kJ/mol'/' ?>')
      READ(*,*)DEMAX
      IF(KC.EQ.1)THEN
      WRITE(*,112)
112   FORMAT(/' Key in the maximum value of the distance (A)',
     1 ' to have the contacts quoted'/' ?>')
      READ(*,*)DMAX
C-----Write the head of the output
      WRITE(IO,113)
113   FORMAT(/' Rotation of the fragment')
      WRITE(IO,114)((ATR1(I,L),L=1,6),I=1,NG1)
114   FORMAT(12(1X,6A1)/12(1X,6A1))
      WRITE(IO,115)((ATR1(I,L),L=1,6),I=1,2),((AT1(I,L),L=1,6),I=1,4)
115   FORMAT(' about the ',6A1,'-',6A1,' direction [Phi(1)]'/
     1 ' Torsion  ',3(6A1,'-'),6A1,2X,'[Tau(1)]'/)
      WRITE(IO,116)
116   FORMAT(' Rotations are counter-clockwise'/)
      WRITE(IO,117)((AT1(I,L),L=1,6),I=1,4)
117   FORMAT(' Non-bonded difference molecular energy'//' Rotation',7X,
     1 'Torsion',16X,'DE'/'  angle',8X,2(6A1,'-'),
     2 5X,'kcal/mol',5X,'kJ/mol'/15X,2('-',6A1)/)
      ELSE
      WRITE(IO,113)
      WRITE(IO,114)((ATR1(I,L),L=1,6),I=1,NG1)
      WRITE(IO,115)((ATR1(I,L),L=1,6),I=1,2),((AT1(I,L),L=1,6),I=1,4)
      WRITE(IO,113)
      WRITE(IO,114)((ATR2(I,L),L=1,6),I=1,NG2)
      WRITE(IO,118)((ATR2(I,L),L=1,6),I=1,2),((AT2(I,L),L=1,6),I=1,4)
118   FORMAT(' about the ',6A1,'-',6A1,' direction [Phi(2)]'/
     1 ' Torsion  ',3(6A1,'-'),6A1,2X,'[Tau(2)]'/)
      WRITE(IO,116)
      ENDIF
C-----Ask for the complete list of data
      IF(KC.EQ.2)THEN
      WRITE(*,119)
119   FORMAT(' Do you want the complete list of data? (Y/N)'/' ?>')
      READ(*,'(A1)')LP
      IF(LP.EQ.QL(1).OR.LP.EQ.QL(2)) GO TO 24
C-----Write the head of the table giving the complete list of the data
      WRITE(IO,120)
120   FORMAT(/' Rotation angles',5X,'Torsion angles',9X,
     1 'Difference energy'/2X,'Phi(1)',2X,'Phi(2)',5X,'Tau(1)',2X,
     2 'Tau(2)',10X,'kcal/mol  kJ/mol')
C-----Ask for the data to use with graphics
24    WRITE(*,121)
121   FORMAT(' Do you want the list of T1, T2, DE data '/
     1 ' to use as input with plotting routines? (Y/N)'/' ?>')
      READ(*,'(A1)')QS
      ENDIF
      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)
      IF(KC.EQ.2)THEN
C-----Modify the coordinates of the atoms of the second fragment,
C-----if this fragment is part of the first one
      DO 30 I=1,NG1
      DO 29 J=3,NG2
      DO 26 L=1,6
      IF(ATR2(J,L).NE.ATR1(I,L)) GO TO 29
26    CONTINUE
      DO 28 K=1,3
      X1K(K)=0.0
      DO 27 KK=1,3
27    X1K(K)=X1K(K)+BI1(K,KK)*XOO1(I,KK)
28    XO2(J,K)=X1K(K)+RR1(K)
      GO TO 30
29    CONTINUE
30    CONTINUE
      ENDIF
C-----The torsion angle for the first fragment is calculated
      CALL FRTORS(NG1,NT1,NA1,AT1,NN1,TAU1,XRO1,R1)
      KD=1
      CALL ENERMOL(KC,KD,NG1,NG2,NG3,XRO1,XRO2,XO3,E1)
      IF(KC.EQ.1) THEN
      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
      WRITE(IO,122)TH(NN1),TAU1(NN1),DE,DEJ1(NN1)
122   FORMAT(F7.1,7X,F8.1,3X,2F12.2)
      WRITE(16,123)TH(NN1)
123   FORMAT(' Rotation angle:',F7.1,' degr.'/)
      ELSE
      NE2=0
      NN2=0
      T2=T20
      CALL NEWCOO(NG2,XO2,BI2,XOO2,RR2)
      ENDIF
C-----Start the loop of the second fragment rotations
      IF(KC.EQ.2)THEN
31    NE2=NE2+1
      CALL ROTCOO(NG2,XOO2,T2,BI2,XRO2,RR2)
      CALL FRTORS(NG2,NT2,NA2,AT2,NN2,TAU2,XRO2,R2)
C-----Calculate the non bonded energy for fragment 2
      KD=2
      CALL ENERMOL(KC,KD,NG1,NG2,NG3,XRO1,XRO2,XO3,E2)
      DE=(E1+E2)-(E01+E02)
      DEJ=DE*4.184
      IF(DEJ.GT.DEMAX)DEJ=DEMAX
      DE=DEJ/4.184
      TH1=T1/RD
      TH2=T2/RD
      IF(LP.EQ.QL(1).OR.LP.EQ.QL(2)) GO TO 32
C-----Printing of the complete table of energy calculations
      WRITE(IO,124)TH1,TH2,TAU1(NN1),TAU2(NN2),DE,DEJ
124   FORMAT(2(F7.1,1X),3X,2F8.1,6X,2F10.1)
32    ENRG(NE1,NE2)=DEJ
      NX=NX+1
      XYZ(NX,1)=TH1
      XYZ(NX,2)=TH2
      XYZ(NX,3)=DEJ
C-----Rounding off the angular values
      IF(IO1(NE1).GE.0)IO1(NE1)=TH1+0.5
      IF(IO1(NE1).LT.0)IO1(NE1)=TH1-0.5
      IF(IO2(NE2).GE.0)IO2(NE2)=TH2+0.5
      IF(IO2(NE2).LT.0)IO2(NE2)=TH2-0.5
      T2=T2+DT2
      IF(T2.GT.T4) GO TO 33
      GO TO 31
C-----End of the second fragment rotation loop
33    ENDIF
      T1=T1+DT1
      IF(T1.GT.T3)GO TO 34
      GO TO 25
C-----End of the first fragment rotation loop
34    IF(KC.EQ.2)THEN
      WRITE(IO,'(1X)')
C-----Scaling of the energy data
      EMX=ENRG(1,1)
      EMN=ENRG(1,1)
      DO 35 I=1,NE1
      DO 35 J=1,NE2
      IF(ENRG(I,J).GT.EMX)EMX=ENRG(I,J)
      IF(ENRG(I,J).LT.EMN)EMN=ENRG(I,J)
35    CONTINUE
      EINT=EMX-EMN
      FT=1
      IF(ABS(EINT).LE.10)FT=10
      IF(ABS(EINT).GT.100)FT=0.1
      IF(ABS(EINT).GT.1000)FT=0.01
      IF(ABS(EINT).GT.10000)FT=0.001
      IF(ABS(EINT).GT.100000)FT=0.0001
      WRITE(IO,125)FT
125   FORMAT(' The following energy values have been multiplied by',
     1 F8.4)
C-----Rounding off the energy values
      DO 36 I=1,NE1
      DO 36 J=1,NE2
      IF(ENRG(I,J).GE.0.0)IENRG(I,J)=ENRG(I,J)*FT+0.5
      IF(ENRG(I,J).LT.0.0)IENRG(I,J)=ENRG(I,J)*FT-0.5
36    CONTINUE
C-----Printing of the three-dimensional energy table
      IF(NE2.GT.20)GO TO 38
      WRITE(IO,126)(IO2(I),I=1,NE2)
126   FORMAT(/'  phi2³',20I4)
      WRITE(IO,127)
127   FORMAT(' _____³',80('_')/' phi1 ³')
      WRITE(IO,128)
128   FORMAT(6X,'³')
      DO 37 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=1,NE2)
129   FORMAT(I5,1X,'³',20I4)
37    WRITE(IO,128)
      GO TO 50
38    IF(NE2.GT.40)GO TO 41
      WRITE(IO,126)(IO2(I),I=1,20)
      WRITE(IO,127)
      WRITE(IO,128)
      DO 39 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=1,20)
39    WRITE(IO,128)
      WRITE(IO,126)(IO2(I),I=21,NE2)
      WRITE(IO,127)
      WRITE(IO,128)
      DO 40 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=21,NE2)
40    WRITE(IO,128)
      GO TO 50
41    IF(NE2.GT.60)GO TO 45
      WRITE(IO,126)(IO2(I),I=1,20)
      WRITE(IO,127)
      WRITE(IO,128)
      DO 42 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=1,20)
42    WRITE(IO,128)
      WRITE(IO,126)(IO2(I),I=21,40)
      WRITE(IO,127)
      WRITE(IO,128)
      DO 43 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=21,40)
43    WRITE(IO,128)
      WRITE(IO,126)(IO2(I),I=41,NE2)
      WRITE(IO,127)
      WRITE(IO,128)
      DO 44 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=41,NE2)
44    WRITE(IO,128)
45    WRITE(IO,126)(IO2(I),I=1,20)
      WRITE(IO,127)
      WRITE(IO,128)
      DO 46 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=1,20)
46    WRITE(IO,128)
      WRITE(IO,126)(IO2(I),I=21,40)
      WRITE(IO,127)
      WRITE(IO,128)
      DO 47 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=21,40)
47    WRITE(IO,128)
      WRITE(IO,126)(IO2(I),I=41,60)
      WRITE(IO,127)
      WRITE(IO,128)
      DO 48 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=41,60)
48    WRITE(IO,128)
      WRITE(IO,126)(IO2(I),I=61,NE2)
      WRITE(IO,127)
      WRITE(IO,128)
      DO 49 I=1,NE1
      WRITE(IO,129)IO1(I),(IENRG(I,K),K=61,NE2)
49    WRITE(IO,128)
50    IF(QS.EQ.QL(1).OR.QS.EQ.QL(2))GO TO 51
      WRITE(*,130)
130   FORMAT(' Key in the name of the file for graphic routines'/' ?>')
      READ(*,'(A14)')FILSRF
      OPEN(UNIT=17,FILE=FILSRF,FORM='FORMATTED')
      REWIND(UNIT=17)
      WRITE(17,131)((XYZ(I,K),K=1,3),I=1,NX)
131   FORMAT(2F8.1,1X,F15.2)
51    CONTINUE
      ELSE
      WRITE(IO,132)DMAX
132   FORMAT(/' Contacts less than',F5.2,' A ','after rotations of the')
      WRITE(IO,133)((ATR1(I,L),L=1,6),I=1,NG1)
133   FORMAT(12(1X,6A1)/12(1X,6A1))
      WRITE(IO,134)((ATR1(I,L),L=1,6),I=1,2)
134   FORMAT(' fragment about the ',6A1,'-',6A1,' direction'/)
      REWIND(UNIT=16)
      DO 52 I=1,10000
      READ(16,'(80A1)',END=53)TIT
52    WRITE(IO,'(80A1)')TIT
53    WRITE(*,135)
135   FORMAT(' Do you want data suitable for graphic software? (Y/N)'/
     1 ' ?>')
      READ(*,'(A1)')LR
      IF(LR.EQ.QL(1).OR.LR.EQ.QL(2))GO TO 56
      WRITE(*,136)
136   FORMAT('   Key in:'/'   0 if D(E) is a function of the rotation',
     1 ' angle'/'   1 if D(E) is a function of the torsion angle'/
     2 ' ?>')
      READ(*,*)KLD
      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
      IF(KLD.EQ.1)GO TO 54
      WRITE(17,138)TH(I),DEJ1(I)
138   FORMAT(F11.1,1X,F15.3)
      GO TO 55
54    WRITE(17,138)TAU1(I),DEJ1(I)
55    CONTINUE
      WRITE(*,139)
139   FORMAT(' Key in the name of the file for VISION'/' ?>')
      READ(*,'(A14)')FILVIS
      OPEN(UNIT=18,FILE=FILVIS,FORM='FORMATTED')
      REWIND(UNIT=18)
      WRITE(18,'(8F10.4)')(DEJ1(I),I=1,NN1)
56    CONTINUE
      ENDIF
C-----Print computing time and date
      CALL GETDAT(IYR,IMT,IDY)
      CALL GETTIM(IHR,IMN,ISC,IHS)
      HS=IHS
      TM=(IHR*3600)+(IMN*60)+ISC+HS/100.0
      IF(TM.LT.TM0)TM=TM+86400
      DTM=TM-TM0
      WRITE(IO,'(/'' Total computing time: '',F9.2,'' s'')')DTM
      WRITE(IO,140)
      WRITE(IO,140)IDY,MONTH(IMT),IYR,IHR,IMN,ISC,IHS
140   FORMAT(I4,'/',A4,'/',I4,3X,I4,':',I2,':',I2,'.',I2/)
      CLOSE(IN)
      CLOSE(IO)
      CLOSE(16)
      CLOSE(17)
      CLOSE(18)
      STOP
      END
 
      SUBROUTINE INPOUT
C-----Assign input and output files
      CHARACTER*14 FILIN
      COMMON/INPO/IN,IO
      WRITE(*,100)
100   FORMAT(' Key in the input file'/' ?>')
      READ(*,'(A14)')FILIN
      OPEN(UNIT=IN,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=IN)
      WRITE(*,101)
101   FORMAT(' Key in the output file'/' ?>')
      READ(*,'(A14)')FILIN
      OPEN(UNIT=IO,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=IO)
      OPEN(UNIT=16,STATUS='SCRATCH',FORM='FORMATTED')
      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)
      COMMON/INPO/IN,IO
     1      /MIMA/MINU,MAIU
      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 ATORS(N,AT1,NG1,NT1,NA1,R1)
C-----Considers the atoms involved in the torsion angles
      CHARACTER*1 AT(500,6),AT1(4,6)
      DIMENSION R1(4,3),XO(500,3),NT1(4),NA1(500)
      COMMON/GRP2/AT
     1      /GRP5/XO
      DO 5 I=1,4
      DO 2 J=1,N
      DO 1 L=1,6
      IF(AT1(I,L).NE.AT(J,L))GO TO 2
1     CONTINUE
      GO TO 3
2     CONTINUE
      WRITE(*,100)(AT1(I,L),L=1,6)
100   FORMAT(' Atom ',6A1,' does not exists in the molecule.'/
     1 ' Check your data. Calculation must stop.')
      STOP
3     DO 4 K=1,3
4     R1(I,K)=XO(J,K)
      NT1(I)=J
5     CONTINUE
      DO 6 I=3,NG1
      IF(NT1(1).EQ.NA1(I)) GO TO 7
6     CONTINUE
      RETURN
7     IRIS=NT1(1)
      NT1(1)=NT1(4)
      NT1(4)=IRIS
      IRIS=NT1(2)
      NT1(2)=NT1(3)
      NT1(3)=IRIS
      DO 8 K=1,3
      RIS=R1(1,K)
      R1(1,K)=R1(4,K)
      R1(4,K)=RIS
      RIS=R1(2,K)
      R1(2,K)=R1(3,K)
8     R1(3,K)=RIS
      RETURN
      END
 
      SUBROUTINE FRTORS(NG1,NT1,NA1,AT1,NN1,TAU1,XRO1,R1)
C-----Considers the atoms defining the torsion angle involved in the
C-----rotating fragment
      CHARACTER*1 AT1(4,6)
      DIMENSION NT1(4),NA1(500),TAU1(180),R1(4,3),XRO1(500,3)
      DO 1 J=1,NG1
      IF(NT1(4).EQ.NA1(J))GO TO 2
1     CONTINUE
      WRITE(*,100) (AT1(NT1(4),L),L=1,6)
100   FORMAT(' Atom ',6A1,' does not exists in the molecule.'/
     1 ' Check your data. Calculation must stop.')
      STOP
2     DO 3 K=1,3
3     R1(4,K)=XRO1(J,K)
      CALL TORS(R1,TAU)
      NN1=NN1+1
      TAU1(NN1)=TAU
      RETURN
      END
 
      SUBROUTINE ENERMOL(KC,KD,NG1,NG2,NG3,XO1,XO2,XO3,ET)
C-----This routine calculates the interatomic contacts and the
C-----van der Waals energy
      CHARACTER*1 ATOM(2,2),AT(500,6)
      CHARACTER*2 KS(2),SDL2(500)
      DIMENSION XO1(500,3),XO2(500,3),XO3(500,3),X1(3),X2(3)
      COMMON/GRP2/AT
     1      /GRP3/MC(400,8)
     2      /GRP4/ATOM
     3      /GRP6/NA1(500),NA2(500),NA3(500)
     4      /GRP7/Q1(500),Q2(500),Q3(500)
     5      /GRP8/DMAX
     6      /GRG1/SDL2
      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.10.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
      ET=ET+E+EQ
5     CONTINUE
      ELSEIF(KD.EQ.1)THEN
      DO 12 I=1,NG3
      DO 12 J=3,NG1
      DO 7 K=1,8
      IF(NA1(J).EQ.MC(NA3(I),K))GO TO 12
      DO 6 L=1,8
      IF(NA1(J).EQ.MC(MC(NA3(I),K),L))GO TO 12
6     CONTINUE
7     CONTINUE
      IF(KC.EQ.2)GO TO 9
      DO 8 L=3,NG2
      IF(NA1(J).EQ.NA2(L))GO TO 12
8     CONTINUE
9     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.10.0)GO TO 12
      KS(1)=SDL2(NA3(I))
      KS(2)=SDL2(NA1(J))
      CALL ENERPOT(KS,D,E,IND)
      IF(IND.EQ.1)E=0
      EQ=(Q3(I)*Q1(J)/D)*331.979
      ET=ET+E+EQ
      IF(KC.EQ.2)GO TO 12
      IF(D.LE.DMAX)WRITE(16,100)(AT(NA3(I),L),L=1,6),
     1 (AT(NA1(J),L),L=1,6),D,E,ET
100   FORMAT(1X,6A1,'...',6A1,5X,'d =',F8.5,8X,'E =',F10.3,3X,'ET=',
     1 F10.3)
12    CONTINUE
      DO 19 I=3,NG1-1
      DO 19 J=I+1,NG1
      DO 14 K=1,8
      IF(NA1(J).EQ.MC(NA1(I),K))GO TO 19
      DO 13 L=1,8
      IF(NA1(J).EQ.MC(MC(NA1(I),K),L))GO TO 19
13    CONTINUE
14    CONTINUE
      IF(KC.EQ.2)GO TO 16
      DO 15 L=3,NG2
      IF(NA1(J).EQ.NA2(L))GO TO 19
15    CONTINUE
16    DO 17 K=1,3
      X1(K)=XO1(I,K)
17    X2(K)=XO1(J,K)
      CALL DIST(X1,X2,D)
      IF (D.LT.0.1) GO TO 19
      IF(D.GT.10.0)GO TO 19
      KS(1)=SDL2(NA1(I))
      KS(2)=SDL2(NA1(J))
      CALL ENERPOT(KS,D,E,IND)
      IF(IND.EQ.1)E=0
      EQ=(Q1(I)*Q1(J)/D)*331.979
      ET=ET+E+EQ
19    CONTINUE
      IF(KC.EQ.1) GO TO 36
      ELSEIF(KD.EQ.2)THEN
      DO 24 I=1,NG3
      DO 24 J=3,NG2
      DO 21 K=1,8
      IF(NA2(J).EQ.MC(NA3(I),K))GO TO 24
      DO 20 L=1,8
      IF(NA2(J).EQ.MC(MC(NA3(I),K),L))GO TO 24
20    CONTINUE
21    CONTINUE
      DO 22 K=1,3
      X1(K)=XO3(I,K)
22    X2(K)=XO2(J,K)
      CALL DIST(X1,X2,D)
      IF (D.LT.0.1) GO TO 24
      IF(D.GT.10.0)GO TO 24
      KS(1)=SDL2(NA3(I))
      KS(2)=SDL2(NA2(J))
      CALL ENERPOT(KS,D,E,IND)
      IF(IND.EQ.1)E=0
      EQ=(Q3(I)*Q2(J)/D)*331.979
      ET=ET+E+EQ
24    CONTINUE
      DO 30 I=3,NG1
      DO 30 J=3,NG2
      DO 26 K=1,8
      IF(NA2(J).EQ.MC(NA1(I),K))GO TO 30
      DO 25 L=1,8
      IF(NA2(J).EQ.MC(MC(NA1(I),K),L))GO TO 30
25    CONTINUE
26    CONTINUE
      DO 27 L=3,NG2
      IF(NA1(I).EQ.NA2(L))GO TO 30
27    CONTINUE
      DO 28 K=1,3
      X1(K)=XO1(I,K)
28    X2(K)=XO2(J,K)
      CALL DIST(X1,X2,D)
      IF (D.LT.0.1) GO TO 30
      IF(D.GT.10.0)GO TO 30
      KS(1)=SDL2(NA1(I))
      KS(2)=SDL2(NA2(J))
      CALL ENERPOT(KS,D,E,IND)
      IF(IND.EQ.1)E=0
      EQ=(Q1(I)*Q2(J)/D)*331.979
      ET=ET+E+EQ
30    CONTINUE
      DO 35 I=3,NG2-1
      DO 35 J=I+1,NG2
      DO 32 K=1,8
      IF(NA2(J).EQ.MC(NA2(I),K))GO TO 35
      DO 31 L=1,8
      IF(NA2(J).EQ.MC(MC(NA2(I),K),L))GO TO 35
31    CONTINUE
32    CONTINUE
      DO 33 K=1,3
      X1(K)=XO2(I,K)
33    X2(K)=XO2(J,K)
      CALL DIST(X1,X2,D)
      IF (D.LT.0.1) GO TO 35
      IF(D.GT.10.0)GO TO 35
      KS(1)=SDL2(NA2(I))
      KS(2)=SDL2(NA2(J))
      CALL ENERPOT(KS,D,E,IND)
      IF(IND.EQ.1)E=0
      EQ=(Q2(I)*Q2(J)/D)*331.979
      ET=ET+E+EQ
35    CONTINUE
36    ENDIF
      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
      D=SQRT(D)
      RETURN
      END
 
      SUBROUTINE READAT(NG,T1,T2,DT,ATT,ATR)
C-----This routine reads the data defining the rotation of the fragment
      CHARACTER *1 ATT(4,6),MINUS(2),ATR(500,6)
      COMMON/INPO/IN,IO
      READ(IN,*)NG,T1,T2,DT
      READ(IN,'(4(6A1))')((ATT(I,L),L=1,6),I=1,4)
      DO 1 I=1,4
      DO 1 K=1,3,2
      MINUS(1)=ATT(I,K)
      MINUS(2)=ATT(I,K+1)
      CALL MINUMA(MINUS)
      ATT(I,K)=MINUS(1)
1     ATT(I,K+1)=MINUS(2)
      READ(IN,'(12(6A1))')((ATR(I,L),L=1,6),I=1,NG)
      DO 2 J=1,NG
      DO 2 K=1,3,2
      MINUS(1)=ATR(J,K)
      MINUS(2)=ATR(J,K+1)
      CALL MINUMA(MINUS)
      ATR(J,K)=MINUS(1)
2     ATR(J,K+1)=MINUS(2)
      RETURN
      END
 
      SUBROUTINE ORTHOG(C,A)
C-----This routine calculates the orthogonalization matrix, O(3,3)
      DIMENSION C(3),A(3),SN(3),CS(3)
      COMMON/GRP1/O(3,3)
      RD=0.0174532925
      DO 1 I=1,3
      AR=A(I)*RD
      SN(I)=SIN(AR)
1     CS(I)=COS(AR)
      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
      RETURN
      END
 
      SUBROUTINE MINUMA(STR)
C-----This routine transforms the lower case characters into upper
C-----case characters, in a string of two characters
      CHARACTER*1 MINU(26),MAIU(26),STR(2)
      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'/
      DO 2 I=1,2
      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 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
      S=SQRT(S)
      DO 2 I=1,3
2     X(I)=X(I)/S
      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)
      CALL INVMAT3(B,BI)
      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 INVMAT3(B,BI)
C-----This routine inverts a 3x3 matrix
      DOUBLE PRECISION B(3,3),BI(3,3)
      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 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 TORS(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
      RR=57.295779
      TAU1=ATAN2(STAU,CTAU)
      TAU=RR*TAU1
      IF(TAU.GT.180.)TAU=TAU-360.
4     RETURN
      END
 
      SUBROUTINE BONDI(N)
C-----This routine calculates bond distances and set up the the
C-----connectivity matrix MC
      CHARACTER*1 ATOM(2,2),AT(500,6)
      CHARACTER*2 NOMAT(103),KS(2),SDL2(500)
      DIMENSION RG(2),IRAG(103)
      COMMON/GRP2/AT
     2      /GRP3/MC(400,8)
     3      /GRP5/XO(500,3)
     4      /SMBL/NOMAT
     5      /ATRD/IRAG
     6      /GRG1/SDL2
     7      /INPO/IN,IO
      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
      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
 
      BLOCK DATA
      CHARACTER*2 NOMAT(103)
      COMMON/SMBL/NOMAT
     1      /ATRD/IRAG(103)
      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 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),QS,AAAA(500,6),MINUS(2),LTP(7),
     1 AFR(500,6),AT(500,6)
      DIMENSION XXO(500,3),QF(500),NS(3)
      COMMON/LBL1/C(3),A(3)
     1      /LBL3/IC,NE,NT
     2      /LBL4/STPE
     3      /LBL8/XX(500,3)
     4      /GRP1/O(3,3)
     5      /GRP2/AT
     6      /GRP5/XO(500,3)
     7      /GRP10/X(500,3),Q(500)
     8      /INS4/CS(3),SN(3)
     9      /INPO/IN,IO
      DATA LTP /'P','A','B','C','I','F','R'/
      PG=3.14159265
      RD=57.295779513
      DO 1 I=1,7
      IF(BRV.EQ.LTP(I))GO TO 2
1     CONTINUE
      WRITE(IO,102) BRV
102   FORMAT(/' The symbol ',A1,'for lattice type is uncorrect.'/
     1 ' Symbol P has been assumed'/)
      BRV=LTP(1)
2     CONTINUE
      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 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 STR(20),BRV,STPE(38,20),AT(500,6),ATOM(2,2),
     1 AD2(6),AFR(500,6),AAAA(500,6)
      INTEGER LE1(24,3,3),LE(3,3)
      DIMENSION B(3,3),F1(24,3),F(3),FB(3)
      COMMON/LBL1/C(3),A(3)
     1      /LBL3/IC,NE,NT
     3      /LBL4/STPE
     4      /LBL8/XX(500,3)
     5      /INS4/CS(3),SN(3)
     6      /GRP2/AT
     7      /GRP10/X(500,3),Q(500)
     8      /EQVP/LE1,F1,FB
     9      /COOA/X1(3),X2(3)
      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
      NE=NE+1
      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') GO TO 13
      IF(BRV.NE.'A') GO TO 8
      B(1,2)=1./2.
      B(1,3)=1./2.
      KB=1
      GO TO 13
8     IF(BRV.NE.'B') GO TO 9
      B(1,1)=1./2.
      B(1,3)=1./2.
      KB=1
      GO TO 13
9     IF(BRV.NE.'C') GO TO 10
      B(1,1)=1./2.
      B(1,2)=1./2.
      KB=1
      GO TO 13
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
      GO TO 13
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
      GO TO 13
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
13    CONTINUE
      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 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      /LBL1/C(3),A(3)
     3      /INS4/CS(3),SN(3)
     4      /EQVP/LE1,F1,FB
     5      /LBL4/STPE
     6      /LBL8/XX(500,3)
     7      /INPO/IN,IO
      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 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      /LBL1/C(3),A(3)
     3      /INS4/CS(3),SN(3)
     4      /LBL8/XX(500,3)
     5      /INPO/IN,IO
      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))
      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 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/INPO/IN,IO
      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.' ')GO TO 18
      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
18    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 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/GRP2/AT
     1      /GRP6/NA1(500),NA2(500),NA3(500)
     2      /INPO/IN,IO
      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
