      PROGRAM ATOMCHAR
C***********************************************************************
C* Release September 1996                                              *
C*                                                                     *
C* ATOMCHAR: a computer routine for calculating the atomic charges in  *
C* a molecule and giving the input for ROTENER                         *
C*                           by M.NARDELLI                             *
C* Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica,  *
C* Chimica Fisica della Universita` degli Studi di Parma, Centro di    *
C* Studio per la Strutturistica Diffrattometrica del C.N.R.,           *
C* Viale delle Scienze, I-43100 Parma, Italy.                          *
C*                                                                     *
C* Two options are considered according to:                            *
C*  (1) the method of D.W. Smith (1990). J.Chem.Ed.,67, 559-562,       *
C*  (2) the method of J. Gasteiger & M. Marsili (1980). Tetrahedron,   *
C*      36, 3219-3228                                                  *
C*                                                                     *
C*                             Input data                              *
C*                          ----------------                           *
C* Input of PARST                                                      *
C*                          ----------------                           *
C* The program works interactively and gives data suitable for ROTENER *
C* calculations.                                                       *
C*                                                                     *
C*  Subroutines:                                                       *
C*               INPOU                                                 *
C*               MINUMA                                                *
C*               CONVE                                                 *
C*               LABL                                                  *
C*               VOLMET                                                *
C*               BONDI                                                 *
C*               ATMCHAR                                               *
C*               AVERQ                                                 *
C*               PEOE                                                  *
C*               ASSCOF                                                *
C*               ITPEOE                                                *
C*               BLOCK DATA                                            *
C*               CONNEC                                                *
C***********************************************************************
      CHARACTER*1 TIT(80),FRT(80),SDOL(6,500),MINUS(2),QS,TT,
     1 UNIM(3,5),II(5),EPSTR(50,40),BRV,SDL1(2,500),LA1(2,20),
     2 AAAA(6),LBA(2)
      CHARACTER*2 LA(20),SDL2(500)
      DIMENSION SC(3),SA(3),B(6),SB(6),S(3),ANAA(20)
      EQUIVALENCE (LA(1),LA1(1,1)),(SDL2(1),SDL1(1,1))
      COMMON/GRP1/V(3,500)
     1      /GRP2/SDOL
     2      /GRP3/O(3,3)
     3      /GRP4/MC(500,8)
     4      /GRP5/C(3),A(3)
     5      /GRP6/VO(3,500)
     6      /GRP7/CS(3),SN(3)
     7      /INP/IN,IO
     8      /EQP1/IC,NE,NT
     9      /EQP2/BRV,EPSTR
      COMMON/GRP10/LA
     1      /GRG1/SDL2
     2      /NATM/KY,NNS
      RD=57.295779513
      IN=13
      IO=15
      CALL INPOU
      DO 41 I=1,20
41    LA(I)='  '
C-----Read title, codes and parameters
      READ(IN,'(80A1)')TIT
      WRITE(IO,'(80A1)')TIT
      READ(IN,'(A1)')BRV
      READ(IN,*)N,D3,D4,NH1,NTN,LSP,LSL,NPR,NST,KL,KY,KC,KO,KB,KA,KT,
     1 KD,KE,KOR,NI
      IF(N.LE.500)GO TO 2
      N=500
      WRITE(IO,100)
100   FORMAT(/' There are too many atoms.'/
     1 ' The calculations are carried out on the first 500 atoms.'/)
2     IF(KY.EQ.0)GO TO 1
      READ(IN,*)NNS
      READ(IN,*)(LA(I),ANAA(I),I=1,NNS)
      DO 37 I=1,NNS
      MINUS(1)=LA1(1,I)
      MINUS(2)=LA1(2,I)
      CALL MINUMA(MINUS)
      LA1(1,I)=MINUS(1)
37    LA1(2,I)=MINUS(2)
1     READ(IN,*)C,SC,A,SA
      CALL VOLMET(RD,CSAR,SNAR)
C-----Read atom symbols,coordinates,thermal parameters & e.s.d.'s
      READ(IN,'(80A1)')FRT
      IF(KOR.EQ.0) GO TO 8
      GO TO (9,11,13),KL
8     GO TO (15,15,17,19),KL
9     DO 10 J=1,N
10    READ(IN,FRT)(SDOL(L,J),L=1,6),(V(I,J),I=1,3),(S(I),I=1,3),
     1 (B(I),I=1,6),(SB(I),I=1,6)
      GO TO 22
11    DO 12 J=1,N
12    READ(IN,FRT)(SDOL(L,J),L=1,6),(V(I,J),I=1,3),(B(I),I=1,6),
     1 (S(I),I=1,3),(SB(I),I=1,6)
      GO TO 22
13    DO 14 J=1,N
14    READ(IN,FRT)(SDOL(L,J),L=1,6),(V(I,J),S(I),I=1,3),
     1 (B(I),SB(I),I=1,6)
      GO TO 22
15    DO 16 J=1,N
16    READ(IN,FRT)(SDOL(L,J),L=1,6),(V(I,J),I=1,3),(S(I),I=1,3)
      GO TO 22
17    DO 18 J=1,N
18    READ(IN,FRT)(SDOL(L,J),L=1,6),(V(I,J),S(I),I=1,3)
      GO TO 22
19    DO 21 J=1,N
      READ(IN,101)(SDOL(L,J),L=1,6),(V(I,J),I=1,3),
     1 ((UNIM(I,K),K=1,5),I=1,3)
101   FORMAT(7X,6A1,2X,3F8.5,12X,3(5A1,3X))
      DO 21 I=1,3
      DO 20 JJ=1,5
20    II(JJ)=UNIM(I,JJ)
      CALL CONVE(II,DH)
21    S(I)=DH
22    DO 23 J=1,N
      MINUS(1)=SDOL(1,J)
      MINUS(2)=SDOL(2,J)
      CALL MINUMA(MINUS)
      SDOL(1,J)=MINUS(1)
23    SDOL(2,J)=MINUS(2)
C-----Matrix of the atom symbols
      DO 39 J=1,N
      MINUS(1)=SDOL(1,J)
      MINUS(2)=SDOL(2,J)
      CALL MINUMA(MINUS)
      SDOL(1,J)=MINUS(1)
      SDOL(2,J)=MINUS(2)
      DO 38 L=1,6
38    AAAA(L)=SDOL(L,J)
      CALL LABL(AAAA,LA1,LBA,KKD)
      SDL1(1,J)=LBA(1)
      SDL1(2,J)=LBA(2)
39    CONTINUE
C-----Calculated hydrogens
      IF(NH1.NE.0)THEN
      DO 6 I=1,NH1
6     READ(IN,'(A1)')TT
      ENDIF
C-----Torsion angles formed by non-bonded atoms
      IF(NTN.NE.0)THEN
      DO 7 I=1,NTN
7     READ(IN,'(A1)')TT
      ENDIF
C-----Planes
      IF(LSP.NE.0)THEN
      DO 40 I=1,LSP
      READ(IN,'(I2)')NT
      DO 40 J=12,36,12
      IF(NT.GT.J)READ(IN,'(A1)')TT
40    CONTINUE
      ENDIF
C-----Lines
      IF(LSL.NE.0)THEN
      DO 25 I=1,LSL
      READ(IN,'(I2)')NT
      DO 25 J=12,36,12
25    IF(NT.GT.12)READ(IN,'(A1)')TT
      ENDIF
C-----Ring puckering
      IF(NPR.NE.0)THEN
      DO 26 I=1,NPR
      READ(IN,'(I2)')NT
      DO 26 J=12,36,12
26    IF(NT.GT.12)READ(IN,'(A1)')TT
      ENDIF
C-----Atomic environments
      IF(NST.NE.0)THEN
      DO 27 I=1,NST
      READ(IN,'(1X,I2)')NT
      DO 27 J=12,36,12
27    IF(NT.GT.12)READ(IN,'(A1)')TT
      ENDIF
C-----Equivalent positions
      READ(IN,*)IC,NE,NT
      DO 28 I=1,NE
28    READ(IN,'(40A1)')(EPSTR(I,J),J=1,40)
C-----Orthogonalization matrix
      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 coordinates
      DO 4 J=1,N
      VO(1,J)=V(1,J)*O(1,1)+V(2,J)*O(1,2)+V(3,J)*O(1,3)
      VO(2,J)=V(2,J)*O(2,2)+V(3,J)*O(2,3)
4     VO(3,J)=V(3,J)*O(3,3)
C-----Printing of atomic coordinates
      WRITE(*,114)
114   FORMAT(' Do you want the fractional coordinates printed ? Y/N'/
     1 ' ?>')
      READ(*,'(A1)')QS
      IF(QS.EQ.'N'.OR.QS.EQ.'n')GO TO 29
      WRITE(IO,102)N
102   FORMAT(/' Number of atoms:  ',I3/)
      WRITE(IO,103)
103   FORMAT(' Atomic coordinates'/' Atom',9X,'X/a',10X,'Y/b',
     1 10X,'Z/c')
      DO 5 J=1,N
5     WRITE(IO,104)(SDOL(L,J),L=1,6),(V(I,J),I=1,3)
104   FORMAT(1X,6A1,4X,3(F8.5,5X))
29    WRITE(*,105)
105   FORMAT(' Do you want the orthogonal coordinates printed ? Y/N'/
     1 ' ?>')
      READ(*,'(A1)')QS
      IF(QS.EQ.'N'.OR.QS.EQ.'n')GO TO 31
      WRITE(IO,106)((O(I,J),J=1,3),I=1,3)
106   FORMAT(/' Orthogonal coordinates (Angstrom)'/
     1 1X,'Orthogonalization matrix:'/
     2 1X,'a  b cosgamma        c cosbeta     ',1X,3F10.5/
     3 1X,'0  b singamma  -c sinbeta cosalpha*',1X,3F10.5/
     4 1X,'0       0       c sinbeta sinalpha*',1X,3F10.5//
     5 1X,'Atom',8X,'X',11X,'Y',10X,'Z')
      DO 30 J=1,N
30    WRITE(IO,107)(SDOL(L,J),L=1,6),(VO(I,J),I=1,3)
107   FORMAT(1X,6A1,2X,3(F8.4,4X))
C-----Interatomic distances and connectivity matrix
31    KDC=1
      WRITE(*,108)
108   FORMAT(' Do you want the interatomic distances printed ? Y/N'/
     1 ' ?>')
      READ(*,'(A1)')QS
      IF(QS.EQ.'N'.OR.QS.EQ.'n') KDC=0
      CALL BONDI(N,NL,KDC)
      IF(KDC.EQ.1) WRITE(IO,109)NL
109   FORMAT(' Number of bond distances:  ',I4/)
C-----Printing of the connectivity matrix
      WRITE(*,110)
110   FORMAT(' Do you want the connectivity matrix printed ? Y/N'/
     1 ' ?>')
      READ(*,'(A1)')QS
      IF(QS.EQ.'N'.OR.QS.EQ.'n')GO TO 32
      WRITE(IO,111)
111   FORMAT(/' Connectivity matrix'/)
      CALL CONNEC(N)
C     DO 7 I=1,N
C7    WRITE(IO,'(I4,4X,8I5)')I,(MC(I,J),J=1,8)
C-----Calculation of electronegativities and atomic charges
32    WRITE(*,112)
112   FORMAT(' If you want charges:'/
     1       '    according to D.W. Smith,                key in  1'/
     1       '    according to J. Gasteiger & M. Marsili, key in  2'/
     2       ' 1/2 ?>')
      READ(*,*)KD
      GO TO (33,34),KD
33    CALL ATMCHAR(RD,N,TIT)
      GO TO 35
34    CALL PEOE(RD,N,TIT)
35    STOP
      END
 
      SUBROUTINE INPOU
C-----This routine defines the input and the output files
      CHARACTER*14 FILIN,FILOUT
      COMMON/INP/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)')FILOUT
      OPEN(UNIT=IO,FILE=FILOUT,FORM='FORMATTED')
      REWIND(UNIT=IO)
      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)
      COMMON/MIMA/MINU,MAIU
      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 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 LABL(AAAA,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/INP/IN,IO
     1      /MIMA/MINU,MAIU
     2      /NATM/KY,NNS
      IF(KY.EQ.1)THEN
      DO 1 J=1,NNS
      IF(AAAA(1).EQ.LA1(1,J).AND.AAAA(2).EQ.LA1(2,J))GO TO 3
1     CONTINUE
      DO 2 K=1,NNS
      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 VOLMET(RD,CSAR,SNAR)
C-----This routine calculates the unit cell volume, cos(alpha*),
C-----sin(alpha*)
      CHARACTER*1 QS
      DIMENSION T(3),CR(3)
      COMMON/GRP5/C(3),A(3)
     1      /GRP7/CS(3),SN(3)
     2      /INP/IN,IO
      WRITE(*,100)
100   FORMAT(' Do you want crystal data printed ? Y/N'/
     1 ' ?>')
      READ(*,'(A1)')QS
      IF(QS.EQ.'N'.OR.QS.EQ.'n')GO TO 1
      WRITE(IO,101)
101   FORMAT(/' Crystal data')
      WRITE(IO,102)(C(I),A(I),I=1,3)
102   FORMAT(' a = ',F7.4,25X,'alpha= ',F6.2/
     1       ' b = ',F7.4,25X,'beta = ',F6.2/
     2       ' c = ',F7.4,25X,'gamma= ',F6.2)
1     DO 3 I=1,3
      A(I)=A(I)/RD
2     CS(I)=COS(A(I))
      IF(ABS(CS(I)).LT.1E-6)CS(I)=0.
3     SN(I)=SIN(A(I))
      R=1.-CS(1)*CS(1)-CS(2)*CS(2)+CS(3)*(2.*CS(1)*CS(2)-CS(3))
      VL=C(1)*C(2)*C(3)*SQRT(R)
      T(1)=CS(2)*CS(3)-CS(1)
      T(2)=CS(1)*CS(3)-CS(2)
      T(3)=CS(1)*CS(2)-CS(3)
      CR(1)=C(2)*C(3)*SN(1)/VL
      CR(2)=C(1)*C(3)*SN(2)/VL
      CR(3)=C(1)*C(2)*SN(3)/VL
      CSAR=T(1)/(SN(2)*SN(3))
      IF(ABS(CSAR).LT.1E-6)CSAR=0.0
      SNAR=SQRT(1.-CSAR*CSAR)
      RETURN
      END
 
      SUBROUTINE BONDI(N,NL,KDC)
C-----This routine calculates bond distances and prints them on request
      CHARACTER*1 SDOL(6,500),KS1(2,2),STR(2),MINU(26),MAIU(26)
      CHARACTER*2 NOMAT(103),SDL2(500),KS(2)
      DIMENSION RG(2),IRAG(103)
      EQUIVALENCE (KS(1),KS1(1,1))
      COMMON/GRP1/V(3,500)
     1      /GRP2/SDOL
     2      /GRP4/MC(500,8)
     2      /GRP5/C(3),A(3)
     3      /GRP6/VO(3,500)
     4      /SMBL/NOMAT
     5      /ATRD/IRAG
     6      /INP/IN,IO
     7      /GRG1/SDL2
     5      /MIMA/MINU,MAIU
      NL=0
      IF(KDC.EQ.0)GO TO 7
      WRITE(IO,100)
100   FORMAT(/' Bond distances (Angstrom)')
      WRITE(IO,101)
101   FORMAT(25X,'Distance')
7     DO 1 J=1,N-1
      KK=0
      DO 1 K=J+1,N
      D1=0.0
      DO 2 I=1,3
      DR=VO(I,J)-VO(I,K)
      IF(ABS(DR).LT.1E-9) DR=0.0
2     D1=D1+DR**2
      D1=SQRT(D1)
      IF(D1.LT.1E-3) GO TO 1
      KS(1)=SDL2(J)
      KS(2)=SDL2(K)
C**************************
      STR(1)=SDOL(1,J)
      STR(2)=SDOL(2,J)
      CALL MINUMA(STR)
      KS1(1,1)=STR(1)
      KS1(2,1)=STR(2)
      STR(1)=SDOL(1,K)
      STR(2)=SDOL(2,K)
      CALL MINUMA(STR)
      KS1(1,2)=STR(1)
      KS1(2,2)=STR(2)
      DO 11 II=1,26
      IF(KS1(2,1).EQ.MAIU(II))GO TO 8
11    CONTINUE
      KS1(2,1)=' '
8     DO 9 II=1,26
      IF(KS1(2,2).EQ.MAIU(II))GO TO 10
9     CONTINUE
      KS1(2,2)=' '
C10   WRITE(IO,'(''KS='',A2,2X,A2)')KS
C**************************
10    DO 6 JJ=1,2
      DO 4 I=1,103
      IF(KS(JJ).EQ.NOMAT(I)) GO TO 5
4     CONTINUE
      WRITE(IO,'(''There is no '',A2,'' atom'')')KS(JJ)
      STOP
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
      MC(J,KK)=K
      NL=NL+1
      IF(KDC.EQ.0)GO TO 1
      WRITE(IO,102)(SDOL(L,J),L=1,6),(SDOL(L,K),L=1,6),D1
102   FORMAT(1X,6A1,'-',6A1,12X,F5.3)
1     CONTINUE
      RETURN
      END
 
      SUBROUTINE ATMCHAR(RD,N,TIT)
C-----This routine calculates the atomic charges by electronegativity
C-----equilibration according to D.W. Smith, J.Chem.Ed.,67, 559-562
C-----(1990)
      CHARACTER*1 SDOL(6,500),SATM(6,8),TIT(80),BRV,
     1 EPSTR(50,40)
      CHARACTER*2 NOMAT(103),SDL2(500),KS(2)
      DIMENSION X0B(8),XB(8),QB(8)
      COMMON/GRP1/V(3,500)
     1      /GRP2/SDOL
     2      /GRP4/MC(500,8)
     3      /GRP5/C(3),A(3)
     4      /GRP8/QT(500,12)
     5      /SMBL/NOMAT
     6      /ELNG/ELN(103)
     7      /INP/IN,IO
     8      /EQP1/IC,NE,NT
     9      /EQP2/BRV,EPSTR
     1      /GRG1/SDL2
      DO 28 I=1,N
      DO 28 J=1,12
28    QT(I,J)=0.0
      WRITE(IO,100)
100   FORMAT(//'Atomic electronegativities and charges (electrons)'/
     1 'according to D.W. Smith, J.Chem.Ed.,67, 559-562, (1990)'/)
      DO 1 I=1,N
      NA=0
      KS(1)=SDL2(I)
      DO 3 K=1,103
      IF(KS(1).EQ.NOMAT(K)) GO TO 4
3     CONTINUE
      GO TO 1
4     X0A=ELN(K)
      DO 29 J=1,8
      DO 29 L=1,6
29    SATM(L,J)=' '
      DO 5 J=1,8
      IF(MC(I,J).EQ.0) GO TO 9
      NA=NA+1
      KS(1)=SDL2(MC(I,J))
      DO 7 K=1,103
      IF(KS(1).EQ.NOMAT(K)) GO TO 8
7     CONTINUE
      GO TO 9
8     DO 15 L=1,6
15    SATM(L,J)=SDOL(L,MC(I,J))
      X0B(J)=ELN(K)
5     CONTINUE
9     IF(N.EQ.2)GO TO 16
      IF(NA.EQ.0.OR.NA.EQ.1) GO TO 1
16    IF(NA.EQ.0)GO TO 1
      SB=0
      DO 10 J=1,NA
10    SB=SB+X0B(J)
      SB=SB/NA
      SBB=0
      DO 11 J=1,NA
      DO 11 K=1,NA
      IF(K.EQ.J) GO TO 11
      SBB=SBB+X0B(J)/X0B(K)
11    CONTINUE
      DNT=SB+(X0A/(NA+1))*(2*NA+SBB)
      XA=((NA+1)*X0A*SB)/DNT
      QA=(XA-X0A)/X0A
      DO 12 J=1,NA
      SAB=0
      DO 13 K=1,NA
      IF(K.EQ.J) GO TO 13
      SAB=SAB+X0A*X0B(K)
13    CONTINUE
      XB(J)=(2*X0A*X0B(J)+SAB)/DNT
12    QB(J)=(XB(J)-X0B(J))/X0B(J)
C**********************
C     WRITE(IO,101)(SDOL(L,I),L=1,6),X0A,XA,QA
C101  FORMAT(/1X,6A1,3X,'X0A=',F5.2,5X,'XA=',F5.2,5X,'QA=',F6.3)
C     DO 14 K=1,NA
C14   WRITE(IO,102)(SATM(L,K),L=1,6),X0B(K),XB(K),QB(K)
C102  FORMAT(1X,6A1,3X,'X0B=',F5.2,5X,'XB=',F5.2,5X,'QB=',F6.3)
C**********************
      DO 23 J=1,N
      DO 17 L=1,6
      IF(SDOL(L,I).NE.SDOL(L,J)) GO TO 23
17    CONTINUE
      IF(QA.EQ.0.0)QA=1E-4
      DO 31 K=1,NA
      IF(QT(J,K).EQ.0.0)GO TO 32
31    CONTINUE
32    DO 33 KK=1,12
      IF(QT(J,KK).EQ.0.0) GO TO 34
33    CONTINUE
34    QT(J,KK)=QA
      IF(N.NE.2)GO TO 35
      QT(1,1)=QA
      QT(2,1)=QB(1)
      GO TO 36
35    DO 22 K=1,NA
      DO 21 JJ=1,N
      DO 18 L=1,6
      IF(SATM(L,K).NE.SDOL(L,JJ)) GO TO 21
18    CONTINUE
      DO 19 KK=1,12
      IF(QT(JJ,KK).EQ.0.0) GO TO 20
19    CONTINUE
      GO TO 22
20    IF(QB(K).EQ.0.0) QB(K)=1E-4
      QT(JJ,KK)=QB(K)
      GO TO 22
21    CONTINUE
22    CONTINUE
      GO TO 1
23    CONTINUE
1     CONTINUE
36    CALL AVERQ(RD,TIT,N,NA)
      STOP
      END
 
      SUBROUTINE AVERQ(RD,TIT,N,NA)
C-----Averages charge values
      CHARACTER*1 SDOL(6,500),TIT(80),QS,BRV,EPSTR(50,40)
      CHARACTER*2 LA(20)
      DIMENSION AV(500)
      COMMON/GRP1/V(3,500)
     1      /GRP2/SDOL
     2      /GRP5/C(3),A(3)
     3      /GRP8/QT(500,12)
     4      /INP/IN,IO
     5      /EQP1/IC,NE,NT
     6      /EQP2/BRV,EPSTR
     7      /GRP10/LA
     8      /NATM/KY,NNS
30    WRITE(*,103)
103   FORMAT(' Do you want the table of the possible charges on the',
     1 ' atoms ? Y/N'/' ?>')
      READ(*,'(A1)')QS
      WRITE(IO,'(1X,'' '')')
      DO 26 I=1,N
      IF(QS.EQ.'N'.OR.QS.EQ.'n') GO TO 38
      WRITE(IO,104)(SDOL(L,I),L=1,6),(QT(I,K),K=1,12)
104   FORMAT(1X,6A1,1X,12F7.3)
38    NA=0
      SS=0
      DO 24 K=1,12
      IF(QT(I,K).EQ.0.0) GO TO 24
      NA=NA+1
      SS=SS+QT(I,K)
24    CONTINUE
      IF(NA.EQ.0)GO TO 26
      AV(I)=SS/NA
26    CONTINUE
      SAV=0.0
      DO 35 I=1,N
35    SAV=SAV+AV(I)
C     WRITE(IO,105)SAV
C105  FORMAT(/'Sum of charges =',F7.3)
      SAV=SAV/N
      SS=0
      DO 36 I=1,N
      AV(I)=AV(I)-SAV
36    SS=SS+AV(I)
      WRITE(*,106)
106   FORMAT(' Do you want the table input for ROTENER ? Y/N'/
     1 ' ?>')
      READ(*,'(A1)')QS
      IF(QS.EQ.'N'.OR.QS.EQ.'n') GO TO 16
      WRITE(IO,107)
107   FORMAT(/'Data for ROTENER'/)
      DO 37 I=1,3
37    A(I)=A(I)*RD
      WRITE(IO,108)(TIT(I),I=1,80),BRV,KY,(C(I),I=1,3),(A(I),I=1,3)
108   FORMAT(80A1/A1/I1/3F9.4,3F8.2)
      WRITE(IO,'(I2)')NNS
      WRITE(IO,111)(LA(I),I=1,NNS)
111   FORMAT(20('''',A2,'''',','))
      WRITE(IO,109)N
109   FORMAT('(1X,6A1,3F9.5,F7.3)'/I3)
      DO 27 I=1,N
27    WRITE(IO,110)(SDOL(L,I),L=1,6),(V(K,I),K=1,3),AV(I)
110   FORMAT(1X,6A1,3F9.5,F7.3)
C-----Insertion of rotation fragments
      WRITE(IO,113)
113   FORMAT(7X,'Insert 1st fragment data:'/'N1,-180,180,5 (or 10)',
     1 15X,'[Format: Free]'/'AAAAA1AAAAA2AAAAA3AAAAA4',12X,
     2 '[Format: 4(6A1)]'/'AAAAA2AAAAA3AAAAA4AAAAA5AAAAA6...',3X,
     3 '[Format: 12(6A1)/12(6A1)]'/
     4 7X,'Insert 2nd fragment (if requested) data:'/'N2,-180,180,10',
     5 22X,'[Format: Free]'/'BBBBB1BBBBB2BBBBB3BBBBB4',12X,
     6 '[Format: 4(6A1)]'/'BBBBB2BBBBB3BBBBB4BBBBB5BBBBB6...',3X,
     7 '[Format: 12(6A1)/12(6A1)]')
C-----Insertion of equivalent positions
      WRITE(IO,'(I2,'','',I2,'','',I1)')IC,NE,NT
      DO 1 I=1,NE
1     WRITE(IO,'(40A1)')(EPSTR(I,J),J=1,40)
16    RETURN
      END
 
      SUBROUTINE PEOE(RD,N,TIT)
C-----Calculates the atomic charges in a molecule according to
C-----J. Gasteiger & M. Marsili (1980). Tetrahedron, 36, 3219-3228.
      CHARACTER*1 SDOL(6,500),SATM(6,8),TIT(80)
      CHARACTER*2 NOMA(9),SDL2(500),KS(2),NM1
      DIMENSION QB(8),CF(3),CFF(3,9)
      COMMON/GRP2/SDOL
     1      /GRP4/MC(500,8)
     2      /GRP8/QT(500,12)
     3      /INS1/SATM
     4      /INS2/QA,QB
     5      /GRP9/COEF(3,14)
     6      /GR10/NOMA
     7      /INP/IN,IO
     8      /GRG1/SDL2
      WRITE(IO,100)
100   FORMAT(//'Atomic charges (electrons) according to'/
     1 'J. Gasteiger & M. Marsili (1980). Tetrahedron, 36, 3219-3228'/)
      DO 28 I=1,N
      DO 28 J=1,12
28    QT(I,J)=0.0
      DO 1 I=1,N
      NA=0
      KS(1)=SDL2(I)
      DO 3 K=1,9
      IF(KS(1).EQ.NOMA(K))GO TO 4
3     CONTINUE
      WRITE(IO,101)KS(1)
101   FORMAT('Atom ',1A2,' cannot be considered.'/
     1 'Repeat the calculation using the Smith''s method.')
      STOP
4     NM1=NOMA(K)
      NN1=0
      DO 13 L=1,8
      IF(MC(I,L).EQ.0) GO TO 14
      NN1=NN1+1
13    CONTINUE
C-----Coefficients for the central atom
14    CALL ASSCOF(NM1,NN1,CF)
      DO 12 KK=1,3
12    CFF(KK,1)=CF(KK)
      DO 29 J=1,8
      DO 29 L=1,6
29    SATM(L,J)=' '
      DO 5 J=1,8
      IF(MC(I,J).EQ.0)GO TO 9
      NA=NA+1
      KS(1)=SDL2(MC(I,J))
      DO 7 K=1,9
      IF(KS(1).EQ.NOMA(K)) GO TO 8
7     CONTINUE
      WRITE(IO,101) KS(1)
      STOP
8     DO 15 L=1,6
15    SATM(L,J)=SDOL(L,MC(I,J))
      DO 11 K=1,N
      IF(MC(I,J).NE.K)GO TO 11
      NN=0
      DO 10 L=1,8
      IF(MC(K,L).EQ.0)GO TO 18
      NN=NN+1
10    CONTINUE
18    NM1=KS(1)
      CALL ASSCOF(NM1,NN,CF)
16    DO 17 KK=1,3
17    CFF(KK,J+1)=CF(KK)
      GO TO 5
11    CONTINUE
5     CONTINUE
9     IF(N.EQ.2)GO TO 26
      IF(NA.EQ.0.OR.NA.EQ.1) GO TO 1
26    IF(NA.EQ.0)GO TO 1
      CALL ITPEOE(CFF,NA)
      DO 23 J=1,N
      DO 24 L=1,6
      IF(SDOL(L,I).NE.SDOL(L,J)) GO TO 23
24    CONTINUE
      IF(QA.EQ.0.0)QA=1E-4
      DO 31 K=1,NA
      IF(QT(J,K).EQ.0.0)GO TO 32
31    CONTINUE
32    DO 33 KK=1,12
      IF(QT(J,KK).EQ.0.0) GO TO 34
33    CONTINUE
34    QT(J,KK)=QA
      IF(N.NE.2)GO TO 35
      QT(1,1)=QA
      QT(2,1)=QB(1)
      GO TO 36
35    DO 22 K=1,NA
      DO 21 JJ=1,N
      DO 25 L=1,6
      IF(SATM(L,K).NE.SDOL(L,JJ)) GO TO 21
25    CONTINUE
      DO 19 KK=1,12
      IF(QT(JJ,KK).EQ.0.0) GO TO 20
19    CONTINUE
      GO TO 22
20    IF(QB(K).EQ.0.0) QB(K)=1E-4
      QT(JJ,KK)=QB(K)
      GO TO 22
21    CONTINUE
22    CONTINUE
      GO TO 1
23    CONTINUE
1     CONTINUE
36    CALL AVERQ(RD,TIT,N,NA)
      RETURN
      END
 
      SUBROUTINE ASSCOF(NMA,NN,CF)
C-----Assignes the coefficients of the expression for orbital
C-----electronegativity on charge
      CHARACTER*2 NMA,NOMA(9)
      DIMENSION CF(3)
      COMMON/GRP9/COEF(3,14)
     1      /GR10/NOMA
      IF(NMA.EQ.'H ')KK=1
      IF(NMA.NE.'C ')GO TO 1
      IF(NN.EQ.4)KK=2
      IF(NN.EQ.3)KK=3
      IF(NN.EQ.2)KK=4
1     IF(NMA.NE.'N ')GO TO 2
      IF(NN.EQ.3)KK=5
      IF(NN.EQ.2)KK=6
      IF(NN.EQ.1)KK=7
2     IF(NMA.NE.'O ')GO TO 3
      IF(NN.EQ.2)KK=8
      IF(NN.EQ.1)KK=9
3     IF(NMA.EQ.'F ')KK=10
      IF(NMA.EQ.'CL')KK=11
      IF(NMA.EQ.'BR')KK=12
      IF(NMA.EQ.'I ')KK=13
      IF(NMA.EQ.'S ')KK=14
      DO 4 I=1,3
4     CF(I)=COEF(I,KK)
      RETURN
      END
 
      SUBROUTINE ITPEOE(CFF,NA)
C-----This subroutine calculates the atomic charges by iterative
C-----partial equalization of orbital electronegativity
      CHARACTER*1 SDOL(6,500),SATM(6,8)
      CHARACTER*2 SDL2(500)
      DIMENSION CFF(3,9),SM(9),Q(9),CHI(3,9),Q0(9),QB(8)
      COMMON/GRP2/SDOL
     1      /INS2/QA,QB
     2      /INS1/SATM
     3      /GRG1/SDL2
      DO 1 J=1,NA+1
      SM(J)=0.0
      Q0(J)=0.0
      DO 1 K=1,3
      SM(J)=SM(J)+CFF(K,J)
1     CHI(K,J)=CFF(K,J)
      NN=0
2     CONTINUE
      NN=NN+1
      Q(1)=0.0
      DO 3 J=2,NA+1
      AA=CHI(1,J)-CHI(1,1)
C     IF(SM(1).GT.SM(J))AA=-AA
3     Q(1)=Q(1)+AA/SM(1)
      Q(1)=Q(1)*0.5**NN
      DO 4 J=2,NA+1
      AA=CHI(1,1)-CHI(1,J)
4     Q(J)=(AA/SM(1))*0.5**NN
      DO 5 J=1,NA+1
5     Q0(J)=Q0(J)+Q(J)
      DO 6 J=1,NA+1
6     CHI(1,J)=CFF(1,J)+CFF(2,J)*Q0(J)+CFF(3,J)*Q0(J)**2
      DO 8 J=1,NA+1
      IF(ABS(Q(J)).GT.0.0005)GO TO 2
8     CONTINUE
      QA=Q0(1)
      DO 10 J=2,NA+1
10    QB(J-1)=Q0(J)
      RETURN
      END
 
      BLOCK DATA
      CHARACTER*1 MINU(26),MAIU(26)
      CHARACTER*2 NOMAT(103),NOMA(9)
      COMMON/SMBL/NOMAT
     1      /ATRD/IRAG(103)
     2      /ELNG/ELN(103)
     3      /GRP9/COEF(3,14)
     4      /GR10/NOMA
     5      /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/
      DATA ELN/2.1,0.0,1.0,1.5,2.0,2.5,3.0,3.5,4.0,0.0,0.9,1.2,1.5,
     1 1.8,2.1,2.5,3.0,0.0,0.8,1.0,1.3,1.5,1.6,1.6,1.5,1.8,1.8,1.8,
     2 1.9,1.6,1.6,1.8,2.0,2.4,2.8,0.0,0.8,1.0,1.2,1.4,1.6,1.8,1.9,
     3 2.2,2.2,2.2,1.9,1.7,1.7,1.8,1.9,2.1,2.5,0.0,0.7,0.9,1.1,1.1,
     4 1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.2,1.2,1.2,1.2,1.2,1.2,1.3,1.5,
     5 1.7,1.9,2.2,2.2,2.2,2.4,1.9,1.8,1.8,1.9,2.0,2.2,0.0,0.7,0.9,
     6 1.1,1.3,1.5,1.7,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3,1.3/
      DATA ((COEF(I,J),I=1,3),J=1,14)/7.17,6.24,-0.56,7.98,9.18,1.88,
     1 8.79,9.32,1.51,10.39,9.45,0.73,11.54,10.82,1.36,12.87,11.15,
     2 0.85,15.58,11.7,-0.27,14.18,12.92,1.39,17.07,13.79,0.47,14.66,
     3 13.85,2.31,11.00,9.69,1.35,10.08,8.47,1.16,9.90,7.96,0.96,
     4 10.14,9.13,1.38/
      DATA NOMA/'H ','C ','N ','O ','F ','CL','BR','I ','S '/
      END
 
      SUBROUTINE CONNEC(N)
C-----From the numerical connectivity matrix gives the same matrix
C-----for the atomic labels
      CHARACTER*1 SDOL(6,500),ATS(6,8)
      DIMENSION NC(8)
      COMMON/GRP2/SDOL
     1      /GRP4/MC(500,8)
     2      /INP/IN,IO
      DO 8 I=1,N
      DO 5 J=1,8
      NC(J)=0
      DO 4 K=1,N
      IF(MC(I,J).NE.K) GO TO 2
      NC(J)=K
      DO 1 L=1,6
1     ATS(L,J)=SDOL(L,K)
      GO TO 5
2     IF(MC(I,J).NE.0)GO TO 4
      DO 3 JJ=J,8
      DO 3 L=1,6
3     ATS(L,JJ)=' '
      GO TO 6
4     CONTINUE
      WRITE(IO,100)
100   FORMAT(' There is something wrong with the atom labels')
      GO TO 7
5     CONTINUE
6     IF(NC(2).EQ.0)THEN
      WRITE(IO,101)I,(SDOL(L,I),L=1,6),NC(1),(ATS(L,1),L=1,6)
101   FORMAT(I3,'(',6A1,'):',2X,I3,'(',6A1,') ')
      ELSEIF(NC(3).EQ.0)THEN
      WRITE(IO,102)I,(SDOL(L,I),L=1,6),(NC(J),(ATS(L,J),L=1,6),J=1,2)
102   FORMAT(I3,'(',6A1,'):',2X,2(I3,'(',6A1,') '))
      ELSEIF(NC(4).EQ.0)THEN
      WRITE(IO,103)I,(SDOL(L,I),L=1,6),(NC(J),(ATS(L,J),L=1,6),J=1,3)
103   FORMAT(I3,'(',6A1,'):',2X,3(I3,'(',6A1,') '))
      ELSEIF(NC(5).EQ.0)THEN
      WRITE(IO,104)I,(SDOL(L,I),L=1,6),(NC(J),(ATS(L,J),L=1,6),J=1,4)
104   FORMAT(I3,'(',6A1,'):',2X,4(I3,'(',6A1,') '))
      ELSEIF(NC(6).EQ.0)THEN
      WRITE(IO,105)I,(SDOL(L,I),L=1,6),(NC(J),(ATS(L,J),L=1,6),J=1,5)
105   FORMAT(I3,'(',6A1,'):',2X,5(I3,'(',6A1,') '))
      ELSEIF(NC(7).EQ.0)THEN
      WRITE(IO,106)I,(SDOL(L,I),L=1,6),(NC(J),(ATS(L,J),L=1,6),J=1,6)
106   FORMAT(I3,'(',6A1,'):',2X,6(I3,'(',6A1,') '))
      ELSEIF(NC(8).EQ.0)THEN
      WRITE(IO,107)I,(SDOL(L,I),L=1,6),(NC(J),(ATS(L,J),L=1,6),J=1,7)
107   FORMAT(I3,'(',6A1,'):',2X,7(I3,'(',6A1,') '))
      ELSE
      WRITE(IO,108)I,(SDOL(L,I),L=1,6),(NC(J),(ATS(L,J),L=1,6),J=1,8)
108   FORMAT(I3,'(',6A1,'):',2X,8(I3,'(',6A1,') '))
      ENDIF
8     CONTINUE
7     RETURN
      END
