      program ESPopt 
      IMPLICIT REAL*8(A-H,O-Z)
C    
C     ESP fit and POP fit restart program  
C     1995. 5.23  by S.N.
C
C     ifsel=1   ESP FIT  restart
C     ifsel=2   POP FIT  restart   Not currently used
C
      read(4,*) ifsel
      if(ifsel.eq.1)then
        call ESFIT(1)
      else if(ifsel.eq.2)then
        call ESFIT(2)
      else
      endif
      stop
      end
C23456789012345678901234567890123456789012345678901234567890123456789012
      SUBROUTINE DUMDIM16
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (I345=345,I346=I345+1,I347=I345+2)
      COMMON/SEKIBN/IX(17000)
      COMMON/COMRPL/IR(I347*47+6+47)
      COMMON/LIMIT /LMT(15)
      COMMON/INT2D /INT(274)
      COMMON/COMIMS/IMS(I346)
      COMMON/COMIOP/IOP(15)
      COMMON/CHCOR /H(I345*I346/2)
      COMMON/COMINV/V(I345*I346)
      COMMON/COMDC /D(I345*I346)
      COMMON/COMP12/P(I345*I346)
      COMMON/COMP34/Q(I345*I346)
C     WRITE(6,*) 'dummy subroutine'
      RETURN
      END
C
C
      SUBROUTINE ESFIT(IREST)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (I345=345,I346=I345+1,I347=I345+2)
      PARAMETER (I200=200,I340=340,I100=70)
C
C     ---------------------------------------------------------
C     IREST = 0 NORMAL RUNNING CONDITION AND NORMAL START
C           = 1 RESTART JOB  (-9 <= ESTAT < 0)
C%SN        = 2 RESTART JOB for bond/atom polarizability fit  (ESTAT < -9)
C     ---------------------------------------------------------
C     IVARY = 0 THE VARIABLE OF THE CENTER IS VARYING
C           = 1 THE VARIABLE OF THE CENTER IS FROZEN
C           = 2 THE VARIABLE OF THE CENTER IS EQUAL TO ONE OF 
C               THE OTHER CENTRE    
C     ---------------------------------------------------------
C
      COMMON/CGSCF3/JGSCF3,ITG,CUTG,SFTG,ESFACT,IESTAT
      COMMON/COMP12/CXNP(I340),CYNP(I340),CZNP(I340),
     & QUPOT(I340,I340)
      COMMON/CHCOR/IATEQ(I200,4),JATEQ(I200,4),NEQ(5),IUNIQ,NVAR,IFIX
     * ,IATOM(I200),IVARY(I200,4),XKPT(I200),YKPT(I200),ZKPT(I200)
     * ,QKPT (I200),VAR (I200*4)
C%SN ,FCAL(1000000)
      COMMON/COMINV/V(1)
      COMMON/COMP34/FCAL(1)
      DIMENSION VV(I340*I340)
      DIMENSION W(I340*I340*I100)
C
C%SN    REST_1 : WRITE restart data to   REST_1
C%SN    REST_2 : READ  restart data from REST_2
C
      OPEN(1,file='REST_1',ACCESS='SEQUENTIAL',FORM='FORMATTED')
      IF(IREST.GE.1)THEN
        OPEN(2,file='REST_2',ACCESS='SEQUENTIAL',FORM='FORMATTED')
      ENDIF
      IF(IREST.GE.2)THEN
        OPEN(3,file='REST_3',ACCESS='SEQUENTIAL',FORM='FORMATTED')
      ENDIF
C
      CALL ESREAD(IREST)
      CLOSE(UNIT=2)
      CLOSE(UNIT=3)
      CALL SETUP(VAR)
      I0= 1+NVAR*(NVAR+1)/2
      I1=I0+NVAR*(NVAR+1)/2
      I2=I1+NVAR
      I3=I2+NVAR
      I4=I3+NVAR
      I5=I4+NVAR
      J1=I5+NVAR
      J2=J1+IESTAT
      IBIG=340*340
      IF(NVAR.GT.70)
     &   STOP 'MEMORY OVER  NVAR'
      IF(IESTAT.GT.IBIG)
     &   STOP 'MEMORY OVER  IESTAT'
      CALL SETBCK(VAR)
      CALL SUMSQ (VAR,FCAL)
      CALL OPTSSQ(VAR,FCAL,W,
     *           V,V(I0),V(I1),V(I2),V(I3),V(I4),V(I5),V(J1),VV)
      CALL DPQUAD
      CALL POTOUT(VAR,FCAL,IREST)
      RETURN
      END
C
C
      SUBROUTINE ESREAD(IREST)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (I345=345,I346=I345+1,I347=I345+2)
      PARAMETER (I200=200,I340=340)
      COMMON/CONSTS/NATM,NDET,NSHL,NOCA,NOCB,NOCC,MDET,NGTO,NATT,NOVER
      COMMON/CGSCF3/JGSCF3,ITG,CUTG,SFTG,ESFACT,IESTAT
      COMMON/COMP12/CXNP(I340),CYNP(I340),CZNP(I340),
     & QUPOT(I340,I340)
      COMMON/CVLIST/CX(3200),CY(3200),CZ(3200),CV(3200),NPOINT
      COMMON/CHCOR/IATEQ(I200,4),JATEQ(I200,4),NEQ(5),IUNIQ,NVAR,IFIX
     * ,IATOM(I200),IVARY(I200,4),XKPT(I200),YKPT(I200),ZKPT(I200)
     * ,QKPT(I200)
      COMMON/SEKIBN/DIAG(I345),FCHRG(3)
      COMMON/PDPOL/IBON(I200),JBON(I200),IEQBON(I200),IBPAIR,IPTYP,
     & FPCHRG,PKPT(I200*2),NVARP,IVARP(I200),IESTAH,IATID,IATM(I200),
     & IEQATM(I200),IVARA(I200),NVARA,NVART,NVARAN,VECISO,IHTYP,NVARH,
     & ALPTRN(I200),ALPANI(I200),ALPATM(I200)
      DIMENSION XYZKPT(I200,4)
      EQUIVALENCE (XYZKPT(1,1),XKPT(1))
C
C%SN
C%SN    JUNE '92 <
C%SN
C%SN    ESTAT  IREST=0   : start
C%SN                 1   : restart
C%SN                 2   : restart; bond/atom polarizability fit
C%SN
      IATID=0
      IBPAIR=0
      IPTYP=0
      FPCHRG=1.0D0
      VECISO=1.0D0
      IF(IREST.EQ.0) THEN
        IUNIQ=0
        DO 100 I=1,NATM
        IF(DABS(CV(I)).LT.1.D-3) GO TO 100
        IUNIQ=IUNIQ+1
        IATOM(IUNIQ)=CV  (I)+0.01D0
        XKPT (IUNIQ)=CX  (I)
        YKPT (IUNIQ)=CY  (I)
        ZKPT (IUNIQ)=CZ  (I)
        QKPT (IUNIQ)=DIAG(I)
        IVARY(IUNIQ,1)=1
        IVARY(IUNIQ,2)=1
        IVARY(IUNIQ,3)=1
        IVARY(IUNIQ,4)=0
  100   CONTINUE
        IF(IUNIQ .GT.I200 ) STOP '** PLEASE INCREASE I200  **'
C%SN
      ELSE IF(IREST.EQ.1) THEN
        READ (2,'(8X,I5,E24.12)') IUNIQ,FCHRG(2)
        READ (2,'(I2,F16.10,I2,F16.10,I2,F16.10,I2,F16.10,I2)') 
     *   (IATOM(I),XKPT(I),IVARY(I,1),YKPT(I),
     *   IVARY(I,2),ZKPT(I),IVARY(I,3),QKPT(I),
     *   IVARY(I,4),I=1,IUNIQ)
        READ (2,*) IESTAT
        READ (2,'(4F15.10)') 
     *  (CXNP(I),CYNP(I),CZNP(I),QUPOT(I,1),I=1,IESTAT)
C%SN
      ELSE IF(IREST.EQ.2) THEN
C%SN
C%SN IPTYP  : Type of polarizability
C%SN             0 : transvers + anisotropy 
C%SN             1 : transvers  
C%SN             2 : anisotropy 
C%SN             3 : atom + anisotropy
C%SN             4 : atom
C%SN FPCHRG : electric field point charge
C%SN VECISO : vector length for atom/transvers induced dipole
C%SN IHTYP  : Inclusion of hyperpolarizabilities
C%SN             0 : No
C%SN             1 : Yes
C%SN
C%SN IBPAIR : # of bond pair for bond polarizability
C%SN IATID  : # of atom for atom polarizability 
C%SN
        READ (3,*) IPTYP,FPCHRG,VECISO,IHTYP
        WRITE(6,*) '  '
        WRITE(6,'(A25,I5,2F8.3,I5)') 'IPTYP,FPCHRG,VECISO,IHTYP',
     &   IPTYP,FPCHRG,VECISO,IHTYP
        READ (3,*) IBPAIR,IATID
        IF(IHTYP.EQ.0)THEN
          WRITE(6,*) ' Polarizability Fitting'
          IF(IBPAIR.GT.0)THEN
            READ (3,*) (IBON(I),JBON(I),IEQBON(I),I=1,IBPAIR)
          ENDIF
          IF(IATID.GT.0)THEN
            READ (3,*) (IATM(I),IEQATM(I),I=1,IATID)
          ENDIF
        ELSE
          WRITE(6,*) ' Hyperpolarizability Fitting'
          IF(IBPAIR.GT.0)THEN
            READ (3,*) (IBON(I),JBON(I),IEQBON(I),ALPTRN(I),
     &        ALPANI(I),I=1,IBPAIR)
          ENDIF
          IF(IATID.GT.0)THEN
            READ (3,*) (IATM(I),IEQATM(I),ALPATM(I),I=1,IATID)
          ENDIF
        ENDIF
        WRITE(6,*) '  '
        WRITE(6,*) ' Bond/Atom Polarizability Fit '
        IF(IPTYP.EQ.1)THEN
          WRITE(6,*) '  TYPE  transvers'
        ELSE IF(IPTYP.EQ.2) THEN
          WRITE(6,*) '  TYPE  anisotropy'
        ELSE IF(IPTYP.EQ.3) THEN
          WRITE(6,*) '  TYPE  atom and anisotropy'
        ELSE IF(IPTYP.EQ.4) THEN
          WRITE(6,*) '  TYPE  atom'
        ELSE 
          WRITE(6,*) '  TYPE  transvers and anisotropy'
        ENDIF
        IF(IBPAIR.NE.0)THEN
          WRITE(6,*) '  '
          WRITE(6,*) '  No. Bond  i -- j  (similarity)'
          DO I=1,IBPAIR,2
            II=I+1
            WRITE(6,'(2(I5,2X,I5,A3,I5,5X,A1,I2,A1,8X))') 
     &       I,IBON(I),' --',JBON(I),'(',IEQBON(I),')',
     &       II,IBON(II),' --',JBON(II),'(',IEQBON(II),')' 
          ENDDO
        ENDIF
        IF(IATID.NE.0)THEN
          WRITE(6,*) '  '
          WRITE(6,*) '  No. Atom  (similarity) ' 
          DO I=1,IATID
            WRITE(6,'(I5,2X,I5,5X,A1,I2,A1)') 
     &       I,IATM(I),'(',IEQATM(I),')'
          ENDDO
        ENDIF
        WRITE(6,*) '  '
C%SN
C%SN    read first data 
C%SN
        READ (2,'(8X,I5,E24.12)') IUNIQ,FCHRG(2)
        READ (2,'(I2,F16.10,I2,F16.10,I2,F16.10,I2,F16.10,I2)') 
     *   (IATOM(I),XKPT(I),IVARY(I,1),YKPT(I),
     *   IVARY(I,2),ZKPT(I),IVARY(I,3),QKPT(I),
     *   IVARY(I,4),I=1,IUNIQ)
        READ (2,'(8X,I5)') IESTAH
        IESTAT=IESTAH*IESTAH
        DO I=1,IESTAH
          READ (2,'(4F15.10)') 
     *     CXNP(I),CYNP(I),CZNP(I),QUPOT(I,1)
        ENDDO
C%SN
C%SN    read data 2-IESTAH
C%SN
        DO I=2,IESTAH
          DO J=1,IUNIQ+2
            READ (2,'(A4)') DUMMY
          ENDDO 
          DO J=1,IESTAH
            READ (2,'(45X,F15.10)') QUPOT(J,I)
          ENDDO
        ENDDO
        CALL PREPOL
      ENDIF
C%SN
C%SN  JUNE '92 >
C%SN
  999 DO 300 K=1,4
      NEQ(K)=0
      DO 301 I=1,IUNIQ
  301 IF(IVARY(I,K).NE.1.AND.ABS(IVARY(I,K)).LT.8) IVARY(I,K)=0
      DO 310 I=2,IUNIQ
      IF(IVARY(I,K).NE.0                       ) GO TO 310
C     IF(DABS(XYZKPT(I,K)            ).LT.1.D-5) GO TO 310
      DO 311 J=1,I-1
      IF(IVARY(J,K).NE.0                       ) GO TO 311
      IF(DABS(XYZKPT(J,K)            ).LT.1.D-5) GO TO 311
      IF(DABS(XYZKPT(I,K)-XYZKPT(J,K)).GT.1.D-5) GO TO 312
      NEQ(K)=NEQ(K)+1
      IATEQ(NEQ(K),K)=I
      JATEQ(NEQ(K),K)=J
      IVARY(I,K)=2
      IVARY(J,K)=0
                                                 GO TO 311
  312 IF(DABS(XYZKPT(I,K)+XYZKPT(J,K)).GT.1.D-5) GO TO 311
      NEQ(K)=NEQ(K)+1
      IATEQ(NEQ(K),K)=I
      JATEQ(NEQ(K),K)=J
      IVARY(I,K)=-2
      IVARY(J,K)= 0
  311 CONTINUE
  310 CONTINUE
  300 CONTINUE
C
C%SN
C%SN April '92 <
C%SN Find Fixed Charge
C%SN
      IALL=0
      IFIX=0
      DO I=1,IUNIQ
       IF(IVARY(I,4).EQ.0) THEN
        IFIX=I
        IALL=IALL+1
       ENDIF
      ENDDO
      IF(IALL.EQ.0)THEN
       WRITE(6,*) '**ALL THE CHARGES ARE FIXED**'
      ELSE IF(IALL.EQ.1) THEN
       IVARY(IFIX,4)=0
      ELSE
       IVARY(IFIX,4)=3
       DO I=1,NEQ(4)
        IF(JATEQ(I,4).EQ.IFIX) THEN
         IF(IVARY(IATEQ(I,4),4).EQ.2) IVARY(IATEQ(I,4),4)=3
         IF(IVARY(IATEQ(I,4),4).EQ.-2) IVARY(IATEQ(I,4),4)=-3
        ENDIF
       ENDDO
      ENDIF
C%SN April '92 >
  220 IF(FCHRG(2).LT.0.D0) ICHRG=FCHRG(2)-0.01D0
      IF(FCHRG(2).GE.0.D0) ICHRG=FCHRG(2)+0.01D0
                  FCHRG(2)=ICHRG
C%SN April '92 <
C%SN Find Fixed Charge
C%SN
      CHARGE=0.D0
      NTST=0
      DO I=1,IUNIQ
       IF(IVARY(I,4).NE.3.AND.IVARY(I,4).NE.-3) THEN
        CHARGE=CHARGE+QKPT(I)
       ELSE
        IF(IVARY(I,4).EQ.3)NTST=NTST+1
        IF(IVARY(I,4).EQ.-3)NTST=NTST-1
       ENDIF
      ENDDO
      IF(NTST.NE.0)THEN
       QFIX=(FCHRG(2)-CHARGE)/NTST
       DO I=1,IUNIQ
        IF(IVARY(I,4).EQ.3) QKPT(I)=QFIX
        IF(IVARY(I,4).EQ.-3) QKPT(I)=-1.0D0*QFIX
       ENDDO
      ENDIF
C%SN April '92 >
      WRITE(6,1)
      WRITE(6,'(9X,2I7,2X,F15.7,I2,F15.7,I2,F15.7,I2,F15.7,I2))')
     *(I,IATOM(I),XKPT(I),IVARY(I,1),YKPT(I),IVARY(I,2),
     *            ZKPT(I),IVARY(I,3),QKPT(I),IVARY(I,4),I=1,IUNIQ)
      WRITE(6,'(13X,''TOTAL'',51X,F22.7)') FCHRG(2)
      WRITE(6,'(12X,42(''* ''))')
    1 FORMAT(/12X,42(2H* )/13X,'ATOM   NUCLEAR',18X,
     *     'COORDINATES (AU)',21X,'INITIAL'/20X,'CHARGE',
     *  10X,'X',14X,'Y',14X,'Z',15X,'CHARGE'/12X,42(2H* ))
      RETURN
      END
C
C
      SUBROUTINE SETUP(VAR)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (I345=345,I346=I345+1,I347=I345+2)
      PARAMETER (I200=200,I340=340)
      COMMON/CONSTS/NATM,NDET,NSHL,NOCA,NOCB,NOCC,MDET,NGTO,NATT,NOVER
      COMMON/CVLIST/CX(3200),CY(3200),CZ(3200),CV(3200),NPOINT
      COMMON/CHCOR/IATEQ(I200,4),JATEQ(I200,4),NEQ(5),IUNIQ,NVAR,IFIX
     * ,IATOM(I200),IVARY(I200,4),XKPT(I200),YKPT(I200),ZKPT(I200)
     * ,QKPT (I200)
      COMMON/SEKIBN/DIAG(I345),FCHRG(3)
      COMMON/PDPOL/IBON(I200),JBON(I200),IEQBON(I200),IBPAIR,IPTYP,
     & FPCHRG,PKPT(I200*2),NVARP,IVARP(I200),IESTAH,IATID,IATM(I200),
     & IEQATM(I200),IVARA(I200),NVARA,NVART,NVARAN,VECISO,IHTYP,NVARH,
     & ALPTRN(I200),ALPANI(I200),ALPATM(I200)
      DIMENSION XYZKPT(I200,4)
      DIMENSION JEQBON(I200),JEQATM(I200)
      EQUIVALENCE (XYZKPT(1,1),XKPT(1))
      DIMENSION VAR(1)
      NVARP =0
      NVARQ =0
      NVAR  =0
      DO 1000 I=1,IUNIQ
      IF(IVARY(I,4).NE.0.AND.IVARY(I,4).LT.8) GO TO 1000
      NVARQ=NVARQ+1
      NVAR =NVAR +1
      VAR(NVAR) = QKPT(I)
 1000 CONTINUE
      DO 1100 K=1,3
      DO 1100 I=1,IUNIQ
      IF(IVARY(I,K).NE.0.AND.IVARY(I,K).LT.8) GO TO 1100
      NVAR=NVAR+1
      VAR(NVAR) = XYZKPT(I,K)
 1100 CONTINUE
      NVARC=NVAR-NVARQ
C%SN
      WRITE(6,1) NVARQ,NVARC
    1 FORMAT(/10X,'NO. OF VARYING CHARGES   =',I5
     *       /10X,'NO OF VARYING COORDINATES=',I5)
C%SN  <
C%SN  for bond polarizability
C%SN  
      K=0
      IF(IBPAIR.GT.0)THEN
        DO I=1,IBPAIR
          JEQBON(I)=IEQBON(I)
        ENDDO
        DO I=1,IBPAIR
          IF(JEQBON(I).NE.-100) THEN
            DO J=I+1,IBPAIR
              IF(JEQBON(I).EQ.JEQBON(J)) JEQBON(J)=-100 
            ENDDO
            K=K+1
            IVARP(K)=JEQBON(I)
          ENDIF
        ENDDO
        NVART=0
        NVARAN=0
        IF(IPTYP.EQ.0.OR.IPTYP.EQ.1)THEN
          NVART=K
        WRITE(6,5) NVART
    5   FORMAT(/10X,'NO OF VARING TRANSVERSE POLARIZABIRITY=',I5)
        ENDIF
        IF(IPTYP.EQ.0.OR.IPTYP.EQ.2.OR.IPTYP.EQ.3)THEN
          NVARAN=K
        WRITE(6,2) NVARAN
    2   FORMAT(/10X,'NO OF VARING ANISOTROPY POLARIZABIRITY=',I5)
        ENDIF
        NVARP=NVART+NVARAN
        DO I=1,NVARP
          NVAR=NVAR+1
          VAR(NVAR) = 0.001D0  
        ENDDO
      ENDIF  
C%SN  <
C%SN  for atom polarizability
C%SN  
      K=0
      NVARA=0
      IF(IATID.GT.0)THEN
        DO I=1,IATID
          JEQATM(I)=IEQATM(I)
        ENDDO
        DO I=1,IATID
          IF(JEQATM(I).NE.-100) THEN
            DO J=I+1,IATID
              IF(JEQATM(I).EQ.JEQATM(J)) JEQATM(J)=-100 
            ENDDO
            K=K+1
            IVARA(K)=JEQATM(I)
          ENDIF
        ENDDO
        NVARA=K
        DO I=1,NVARA
          NVAR=NVAR+1
          VAR(NVAR) = 0.002D0  
        ENDDO
        WRITE(6,4) NVARA
    4   FORMAT(/10X,'NO OF VARING ATOM POLARIZABIRITY=',I5)
      ENDIF  
      WRITE(6,3) NVAR
    3 FORMAT(/10X,'TOTAL NO OF VARIABLES    =',I5)
C%SN
C%SN  >
      RETURN
C
      ENTRY SETBCK(VAR)
C
      ixxx=ixxx+1
      NRST  =0
      DO 2000 I=1,IUNIQ
      IF(IVARY(I,4).NE.0.AND.IVARY(I,4).LT.8) GO TO 2000
      NRSTQ=NRSTQ+1
      NRST=NRST+1
      QKPT(I)=VAR(NRST)
 2000 CONTINUE
      DO 2100 K=1,3
      DO 2100 I=1,IUNIQ
      IF(IVARY(I,K).NE.0.AND.IVARY(I,K).LT.8) GO TO 2100
      NRST=NRST+1
      XYZKPT(I,K)=VAR(NRST)
 2100 CONTINUE
C
C%SN     IF(NRST.NE.NVARQ+NVARC) STOP 'NRST.NE.NVARQ+NVARC'
      DO 2500 K=1,4
      IF(NEQ(K).EQ.0) GO TO 2500
      DO 2510 I=1,NEQ(K)
      J=IATEQ(I,K)
                          XYZKPT(J,K)= XYZKPT(JATEQ(I,K),K)
 2510 IF(IVARY(J,K).LT.0) XYZKPT(J,K)=-XYZKPT(J         ,K)
 2500 CONTINUE
      DO I=1,IUNIQ
       IF(IVARY(I,4).EQ.-9.OR.IVARY(I,4).EQ.-8)THEN
        QKPT(I)=QKPT(I+1)*(-1.0D0)
       ENDIF
      ENDDO
C%SN MAY '92 <
      CHARGE=0.D0
      NTST=0
      DO I=1,IUNIQ
       IF(IVARY(I,4).NE.3.AND.IVARY(I,4).NE.-3) THEN
        CHARGE=CHARGE+QKPT(I)
       ELSE
        IF(IVARY(I,4).EQ.3)NTST=NTST+1
        IF(IVARY(I,4).EQ.-3)NTST=NTST-1
       ENDIF
      ENDDO
      IF(NTST.NE.0)THEN
       QFIX=(FCHRG(2)-CHARGE)/NTST
       DO I=1,IUNIQ
        IF(IVARY(I,4).EQ.3) QKPT(I)=QFIX
        IF(IVARY(I,4).EQ.-3) QKPT(I)=-1.0D0*QFIX
       ENDDO
      ENDIF
C%SN MAY '92 >
      KA=0
      IF(IBPAIR.GT.0)THEN
        IF(IPTYP.EQ.1.OR.IPTYP.EQ.2.OR.IPTYP.EQ.3)THEN
          K=NVARP
          KA=IBPAIR
        ELSE
          K=NVARP/2
          KA=IBPAIR*2
        ENDIF
        DO I=1,K
          NRST=NRST+1
          DO J=1,IBPAIR
            IF(IEQBON(J).EQ.IVARP(I))THEN
              PKPT(J)=VAR(NRST)
            ENDIF
          ENDDO
        ENDDO
        IF(IPTYP.EQ.0)THEN
          DO I=1,K
            NRST=NRST+1
            DO J=1,IBPAIR
              IF(IEQBON(J).EQ.IVARP(I))THEN
                PKPT(IBPAIR+J)=VAR(NRST)
              ENDIF
            ENDDO
          ENDDO
        ENDIF
      ENDIF
C%SN
C%SN July '92 >
      IF(IATID.GT.0)THEN
        DO I=1,NVARA
          NRST=NRST+1
          DO J=1,IATID
            IF(IEQATM(J).EQ.IVARA(I))THEN
              PKPT(KA+J)=VAR(NRST)
            ENDIF
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
C
C
      SUBROUTINE POTOUT(VAR,FCAL,IREST)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (I200=200,I340=340)
      COMMON/CGSCF3/JGSCF3,ITG,CUTG,SFTG,ESFACT,IESTAT
      COMMON/COMP12/CXNP(I340),CYNP(I340),CZNP(I340),
     & QUPOT(I340,I340)
      COMMON/CHCOR/IATEQ(I200,4),JATEQ(I200,4),NEQ(5),IUNIQ,NVAR,IFIX
     * ,IATOM(I200),IVARY(I200,4),XKPT(I200),YKPT(I200),ZKPT(I200)
     * ,QKPT (I200)
      COMMON/CONSTR/DIPCAL(3),QUCAL(3),CMAS(3),DIPSUM
      COMMON/PDPOL/IBON(I200),JBON(I200),IEQBON(I200),IBPAIR,IPTYP,
     & FPCHRG,PKPT(I200*2),NVARP,IVARP(I200),IESTAH,IATID,IATM(I200),
     & IEQATM(I200),IVARA(I200),NVARA,NVART,NVARAN,VECISO,IHTYP,NVARH,
     & ALPTRN(I200),ALPANI(I200),ALPATM(I200)
      DIMENSION VAR(1),FCAL(1)
      DIMENSION ALPH(I200)
      CHARACTER*4 DLABL(3)/ 'X  =', 'Y  =', 'Z  ='/
      CHARACTER*5 QLABL(3)/' AA =',' BB =',' CC ='/
C
      CHARGE=0.D0
      DO 20 I=1,IUNIQ
   20 CHARGE=CHARGE+QKPT(I)
      WRITE(6,1)
      WRITE(7,'(9X,2I7,2X,F15.7,I2,F15.7,I2,F15.7,I2,F15.7,I2))')
     *(I,IATOM(I),XKPT(I),IVARY(I,1),YKPT(I),IVARY(I,2),
     *            ZKPT(I),IVARY(I,3),QKPT(I),IVARY(I,4),I=1,IUNIQ)
      WRITE(6,'(9X,2I7,2X,F15.7,I2,F15.7,I2,F15.7,I2,F15.7,I2))')
     *(I,IATOM(I),XKPT(I),IVARY(I,1),YKPT(I),IVARY(I,2),
     *            ZKPT(I),IVARY(I,3),QKPT(I),IVARY(I,4),I=1,IUNIQ)
      WRITE(6,'(13X,''TOTAL'',51X,F22.7)') CHARGE
      WRITE(6,'(12X,42(''* ''))')
    1 FORMAT(/12X,42(2H* )/13X,'ATOM   NUCLEAR',18X,
     *     'COORDINATES (AU)',22X,'FINAL'/20X,'CHARGE',
     *   9X,'X',16X,'Y',16X,'Z',13X,'CHARGE'/12X,42(2H* ))
      WRITE(6,2) ( DLABL(I),DIPCAL(I),QLABL(I),QUCAL(I),DLABL(I),
     *             CMAS(I)*0.529177D0, I=1,3 ) , DIPSUM
    2 FORMAT(/20X,'DIPOLE',12X,'QUADRUPOLE',13X,'CMASS'
     *       /20X,'(DEBYE)',10X,'(BUCKINGHAM)',12X,'(ANG)'
     * /3(16X,A5,F10.5,5X,A5,F10.5,5X,A5,F10.5/),16X,'TOT  ',F10.5/)
C%SN
      J=0
      IF(NVART.GT.0)THEN
        WRITE(6,*) ' '
        WRITE(6,*) '      Parameter of Transverse polarizability'
        WRITE(6,*) '          No. Bond  i -- j  (parm. trans.)'
        DO I=1,IBPAIR
          J=J+1
          WRITE(6,'(8X,I5,2X,I5,A3,I5,5X,A1,F10.4,A1,8X)') 
     &     I,IBON(I),' --',JBON(I),'(',PKPT(J),')' 
        ENDDO
      ENDIF
      IF(NVARAN.GT.0)THEN
        WRITE(6,*) ' '
        WRITE(6,*) '      Parameter of Anisotropy polarizability'
        WRITE(6,*) '          No. Bond  i -- j  (parm. aniso.)'
        DO I=1,IBPAIR
          J=J+1
          WRITE(6,'(8X,I5,2X,I5,A3,I5,5X,A1,F10.4,A1,8X)') 
     &     I,IBON(I),' --',JBON(I),'(',PKPT(J),')' 
        ENDDO
      ENDIF
      IF(NVARA.GT.0)THEN
        WRITE(6,*) ' '
        WRITE(6,*) '      Parameter of Atom polarizability'
        WRITE(6,*) '          No.     Atom  (parm. atom) ' 
        DO I=1,IATID
          J=J+1
          WRITE(6,'(8X,I5,2X,I5,5X,A1,F10.4,A1)') 
     &     I,IATM(I),'(',PKPT(J),')'
        ENDDO
      ENDIF
      WRITE(6,*) ' '
      SSQ = 0.D0
      VSQ = 0.D0
      SMIN=99999.D0
      SMAX=0.0D0
      SMAXN=0.D-5
      KMAX=0
      KMAXN=0
      IF(IBPAIR.EQ.0.AND.IATID.EQ.0)THEN
        DO I=1,IESTAT
          SSQ = SSQ + FCAL(I)*FCAL(I)
          VSQ = VSQ + QUPOT(I,1)*QUPOT(I,1)
          IF(DABS(FCAL(I)).LT.SMIN) SMIN=DABS(FCAL(I))
          IF(DABS(FCAL(I)).GT.SMAX) THEN
            SMAX=DABS(FCAL(I))
            KMAX=I
          ENDIF
          IF(DABS(FCAL(I)).LT.SMAX.AND.DABS(FCAL(I)).GT.SMAXN) THEN
            SMAXN=DABS(FCAL(I))
            KMAXN=I
          ENDIF
        ENDDO
      ELSE
        KK=0
        WRITE(6,*) '   #      Sigma(kcal/mol)    S(%)   Max.#    Maximum 
     & (kcal/mol)'
        DO I=1,IESTAH
          SSQE=0.0D0
          VSQE=0.0D0
          SMAXE=0.0D0
          KMAXE=0
          DO J=1,IESTAH
            KK=KK+1
            FF =  FCAL(KK)*FCAL(KK)
            QQ =  QUPOT(J,I)*QUPOT(J,I)
            SSQE= SSQE+ FF
            VSQE= VSQE+ QQ
            SSQ = SSQ + FF
            VSQ = VSQ + QQ
            IF(DABS(FCAL(KK)).LT.SMIN) SMIN=DABS(FCAL(KK))
            IF(DABS(FCAL(KK)).GT.SMAX) THEN
              SMAX=DABS(FCAL(KK))
              KMAX=KK
            ENDIF
            IF(DABS(FCAL(KK)).LT.SMAX.AND.DABS(FCAL(KK)).GT.SMAXN) THEN
              SMAXN=DABS(FCAL(KK))
              KMAXN=KK
            ENDIF
            IF(DABS(FCAL(KK)).GT.SMAXE) THEN
              SMAXE=DABS(FCAL(KK))
              KMAXE=J
            ENDIF
          ENDDO
          SS = DSQRT(SSQE/DFLOAT(IESTAH))
          RR   = SS/DSQRT(VSQE/DFLOAT(IESTAH))/627.5251387752D0
          WRITE(6,'(I5,2F15.4,I5,F15.4)') I,SS,RR*100,KMAXE,SMAXE
        ENDDO
      ENDIF
      SIGMA = DSQRT(SSQ/DFLOAT(IESTAT))
      RMS   = SIGMA/DSQRT(VSQ/DFLOAT(IESTAT))/627.5251387752D0
C%SN                    AU TO KCAL    627.5251387752
      WRITE(6,1080) SIGMA,RMS*100,SMIN,KMAX,SMAX,KMAXN,SMAXN
 1080 FORMAT( 10X,'STANDARD DEVIATION(SIGMA) =',F15.4,'  kcal/mol'
     *       /10X,'PERCENT  DEVIATION   (S)  =',F15.4,'   %      '  
     *       /10X,'Minimum  DEVIATION        =',F15.4,'  kcal/mol'  
     *       /10X,'Maximum  DEVIATION ',I6,' =',F15.4,'  kcal/mol'  
     *       /10X,'Next Max DEVIATION ',I6,' =',F15.4,'  kcal/mol')  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IF(IREST.NE.2)THEN
        WRITE(1,*) IUNIQ,CHARGE
        WRITE(1,'(I2,F16.10,I2,F16.10,I2,F16.10,I2,F16.10,I2)') 
     *   (IATOM(I),XKPT(I),IVARY(I,1),YKPT(I),IVARY(I,2),
     *        ZKPT(I),IVARY(I,3),QKPT(I),IVARY(I,4),I=1,IUNIQ)
        WRITE(1,*) IESTAT
        WRITE(1,'(4F15.10)') 
     *   (CXNP(I),CYNP(I),CZNP(I),QUPOT(I,1),I=1,IESTAT)
C%SN        REWIND 1
      ELSE
        KK=0
        DO I=1,IESTAH
          WRITE(1,*) IUNIQ,CHARGE
          WRITE(1,'(I2,F16.10,I2,F16.10,I2,F16.10,I2,F16.10,I2)') 
     *     (IATOM(K),XKPT(K),IVARY(K,1),YKPT(K),IVARY(K,2),
     *     ZKPT(K),IVARY(K,3),QKPT(K),IVARY(K,4),K=1,IUNIQ)
          WRITE(1,*) IESTAH
          DO J=1,IESTAH
            KK=KK+1
            QQ=QUPOT(J,I)*627.5251387752D0
            WRITE(1,'(5F15.9)') 
     *       CXNP(J),CYNP(J),CZNP(J),QQ,FCAL(KK)
          ENDDO
        ENDDO
      ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      RETURN
      END
C
C
      SUBROUTINE SUMSQ(VAR,F)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (I200=200,I340=340)
      COMMON/CHCOR/IATEQ(I200,4),JATEQ(I200,4),NEQ(5),IUNIQ,NVAR,IFIX
     * ,IATOM(I200),IVARY(I200,4),XKPT(I200),YKPT(I200),ZKPT(I200)
     * ,QKPT (I200)
      COMMON/CGSCF3/JGSCF3,ITG,CUTG,SFTG,ESFACT,IESTAT
      COMMON/COMP12/CXNP(I340),CYNP(I340),CZNP(I340),
     & QUPOT(I340,I340)
      COMMON/PDPOL/IBON(I200),JBON(I200),IEQBON(I200),IBPAIR,IPTYP,
     & FPCHRG,PKPT(I200*2),NVARP,IVARP(I200),IESTAH,IATID,IATM(I200),
     & IEQATM(I200),IVARA(I200),NVARA,NVART,NVARAN,VECISO,IHTYP,NVARH,
     & ALPTRN(I200),ALPANI(I200),ALPATM(I200)
      COMMON/PRPOL/ XMID(I200),YMID(I200),ZMID(I200),
     &          TRNX1(I340,I200),TRNY1(I340,I200),TRNZ1(I340,I200),
     &          TRNX2(I340,I200),TRNY2(I340,I200),TRNZ2(I340,I200),
     &          BONX(I200),BONY(I200),BONZ(I200),BONLEN(I200),
     &          FLDLEN(I340,I200),ANIFAC(I340,I200),PALEN(I340,I200),
     &          SIGNA(I340,I200),FTMP(I340),
     &          ATMX1(I340,I200),ATMY1(I340,I200),ATMZ1(I340,I200),
     &          ATMX2(I340,I200),ATMY2(I340,I200),ATMZ2(I340,I200) 
      DIMENSION VAR(1),F(1)
      DIMENSION DLTAQ(I200)
      IF(IBPAIR.EQ.0.AND.IATID.EQ.0)THEN
        DO 100 I=1,IESTAT
        SUMP =0.D0
        DO 110 J=1,IUNIQ
        XC     =CXNP(I)-XKPT(J)
        YC     =CYNP(I)-YKPT(J)
        ZC     =CZNP(I)-ZKPT(J)
        RA     =DSQRT(XC*XC+YC*YC+ZC*ZC)
  110   SUMP   =SUMP-QKPT(J)/RA
  100   F(I)   =(QUPOT(I,1) -SUMP)*627.5251387752D0
C                  AU TO KCAL    627.5251387752
      ELSE 
C%SN
C%SN   calculation of electrostatic potential for induced dipoles
C%SN
C%SN   bond/atom polarizability
C%SN
        KK=0
        DO I=1,IESTAH
          DO J=1,IESTAH
            KK=KK+1  
            F(KK)=QUPOT(J,I)+FTMP(J)
          ENDDO
        ENDDO
C%SN
C%SN     transverse, induced diople length is VECISO 
C%SN
        JFIRS=0
        IF(IBPAIR.NE.0)THEN
          KK=0
          IF(IPTYP.EQ.0.OR.IPTYP.EQ.1)THEN
            JFIRS=IBPAIR
            DO K=1,IESTAH
              DO J=1,IBPAIR
C#####
                FLD1=FPCHRG/(FLDLEN(K,J)*FLDLEN(K,J))
                IF(IHTYP.EQ.0)THEN
                  DLTAQ(J)=-FLD1*PKPT(J)/VECISO
                ELSE 
                  FLD2=FLD1*FLD1*FLD1
                  DLTAQ(J)=-FLD1*ALPTRN(J)/VECISO
                  DLTAQ(J)=DLTAQ(J)-FLD2*ALPTRN(J)*PKPT(J)/VECISO
                ENDIF
C#####
              ENDDO
              DO I=1,IESTAH
                KK=KK+1
                DO J=1,IBPAIR
                  FLEN1=DSQRT((CXNP(I)-TRNX1(K,J))**2+
     &                        (CYNP(I)-TRNY1(K,J))**2+
     &                        (CZNP(I)-TRNZ1(K,J))**2)
                  FLEN2=DSQRT((CXNP(I)-TRNX2(K,J))**2+
     &                        (CYNP(I)-TRNY2(K,J))**2+
     &                        (CZNP(I)-TRNZ2(K,J))**2)
                  F(KK)=F(KK)-DLTAQ(J)/FLEN1+DLTAQ(J)/FLEN2
                ENDDO
              ENDDO
            ENDDO
          ENDIF
C%SN
C%SN     anisotoropy, induced dipole length .eq. bond length
C%SN
          IF(IPTYP.EQ.0.OR.IPTYP.EQ.2.OR.IPTYP.EQ.3)THEN
            KK=0
            KFIRS=JFIRS
            JFIRS=JFIRS+IBPAIR
            DO K=1,IESTAH
              DO J=1,IBPAIR
                JJ=KFIRS+J
C#####
                FLD1=FPCHRG*ANIFAC(K,J)/(FLDLEN(K,J)*FLDLEN(K,J))
                IF(IHTYP.EQ.0)THEN
                  DLTAQ(J)=-FLD1*SIGNA(K,J)*PKPT(JJ)/BONLEN(J)
                ELSE 
                  DLTAQ(J)=-FLD1*SIGNA(K,J)*ALPANI(J)/BONLEN(J)
                  FLD2=FLD1*FLD1*FLD1
                  DLTAQ(J)=DLTAQ(J)-FLD2*SIGNA(K,J)*ALPANI(J)*PKPT(JJ)
     &              /BONLEN(J)
                ENDIF
C#####
C               DLTAQ(J)=-FPCHRG*SIGNA(K,J)*ANIFAC(K,J)
C    &                 *PKPT(JJ)/(FLDLEN(K,J)**2*BONLEN(J))
              ENDDO
              DO I=1,IESTAH
                KK=KK+1
                DO J=1,IBPAIR
                  FLEN1=PALEN(I,IBON(J))
                  FLEN2=PALEN(I,JBON(J))
                  F(KK)=F(KK)-DLTAQ(J)/FLEN1+DLTAQ(J)/FLEN2
                ENDDO
              ENDDO
            ENDDO
          ENDIF
        ENDIF
C%SN
C%SN     atom, induced diople length is VECISO 
C%SN
        IF(IATID.NE.0)THEN
          KK=0
          DO K=1,IESTAH
            DO J=1,IATID
              JJ=JFIRS+J
C#####
                FLD1=FPCHRG/(PALEN(K,J)*PALEN(K,J))
                IF(IHTYP.EQ.0)THEN
                  DLTAQ(J)=-FLD1*PKPT(JJ)/VECISO
                ELSE
                  DLTAQ(J)=-FLD1*ALPATM(J)/VECISO
                  FLD2=FLD1*FLD1*FLD1
                  DLTAQ(J)=DLTAQ(J)+FLD2*ALPATM(J)*PKPT(JJ)/VECISO
                ENDIF
C#####
C             DLTAQ(J)=-FPCHRG*PKPT(JJ)/(PALEN(K,J)*PALEN(K,J)
C    &          *VECISO)
            ENDDO
            DO I=1,IESTAH
              KK=KK+1
              DO J=1,IATID
                FLEN1=DSQRT((CXNP(I)-ATMX1(K,J))**2+
     &                      (CYNP(I)-ATMY1(K,J))**2+
     &                      (CZNP(I)-ATMZ1(K,J))**2)
                FLEN2=DSQRT((CXNP(I)-ATMX2(K,J))**2+
     &                      (CYNP(I)-ATMY2(K,J))**2+
     &                      (CZNP(I)-ATMZ2(K,J))**2)
                F(KK)=F(KK)-DLTAQ(J)/FLEN1+DLTAQ(J)/FLEN2
              ENDDO
            ENDDO
          ENDDO
C%SN
        ENDIF
        KK=0
        DO K=1,IESTAH
          DO I=1,IESTAH
            KK=KK+1
            F(KK)=F(KK)*627.5251387752D0
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
C
C
      SUBROUTINE DPQUAD
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (I200=200)
C
C      ROUTINE TO CALCULATE THE DIPOLE & QUADRUPOLE MOMENTS 
C
      COMMON/CHCOR/IATEQ(I200,4),JATEQ(I200,4),NEQ(5),IUNIQ,NVAR,IFIX
     * ,IATOM(I200),IVARY(I200,4),XKPT(I200),YKPT(I200),ZKPT(I200)
     * ,QKPT (I200)
      COMMON/CONSTR/DIPCAL(3),QUCAL(3),XC,YC,ZC,DIPSUM
      DIMENSION AIN(3,3),S(3,3),D(3)
      DIMENSION CO(3,50),WT(40)
C
      DATA (WT(I),I=1,40) /          1.0079D0, 4.0026D0,
     *                     6.9410D0, 9.0122D0,10.8100D0,12.0110D0,
     *                    14.0067D0,15.9994D0,18.9984D0,20.1700D0,
     *                    22.9897D0,24.3050D0,26.9815D0,28.0860D0,
     *                    30.9738D0,32.0600D0,35.4530D0,39.9480D0,
     *                    39.0900D0,40.0800D0, 20*50.D0 /
C
      DIPCAL(1)=0.D0
      DIPCAL(2)=0.D0
      DIPCAL(3)=0.D0
      DO 200 I =1,IUNIQ
      DIPCAL(1) = DIPCAL(1)+QKPT(I)*XKPT(I)
      DIPCAL(2) = DIPCAL(2)+QKPT(I)*YKPT(I)
  200 DIPCAL(3) = DIPCAL(3)+QKPT(I)*ZKPT(I)
      DIPSUM=0.D0
      DO 210 I=1,3
      DIPCAL(I) = DIPCAL(I)*2.54158059D0
  210 DIPSUM=DIPSUM+DIPCAL(I)**2
      DIPSUM=DSQRT(DIPSUM)
C
C CALCULATES THE CENTER OF MASS OF THE MOLECULE
C
      XC   = 0.D0
      YC   = 0.D0
      ZC   = 0.D0
      SUMX = 0.D0
      SUMY = 0.D0
      SUMZ = 0.D0
      SUM  = 0.D0
      DO 3 I =1,IUNIQ
      IF(IATOM(I).LE.0) GO TO 3
      SUMX = SUMX + XKPT(I) * WT(IATOM(I))
      SUMY = SUMY + YKPT(I) * WT(IATOM(I))
      SUMZ = SUMZ + ZKPT(I) * WT(IATOM(I))
      SUM  = SUM  +           WT(IATOM(I))
    3 CONTINUE
      XC   = SUMX/SUM
      YC   = SUMY/SUM
      ZC   = SUMZ/SUM
C
C MOVES THE ORIGIN TO THE CENTER OF MASS.
C
      DO 5 I=1,IUNIQ
      CO(1,I) = XKPT(I) - XC
      CO(2,I) = YKPT(I) - YC
    5 CO(3,I) = ZKPT(I) - ZC
C
C CALCULATES THE MOMENTS OF INERTIA AND THE PRINCIPAL AXES OF
C ROTATION OF THE MOLECULE.  IT THEN REORIENTS THE MOLECULE ALONG 
C THE PRINCIPAL AXES OF ROTATION (WITH THE ORIGIN AT THE CENTER OF MASS)
c IN PREPARATION FOR CALCULATION OF QUADRUPOLE MOMENT COMPONENTS.
C
      SXX=0.D0
      SYY=0.D0
      SZZ=0.D0
      SXY=0.D0
      SXZ=0.D0
      SYZ=0.D0
      DO 10 I=1,IUNIQ
      IF(IATOM(I).LE.0) GO TO 10
      XX=CO(1,I)*CO(1,I)
      YY=CO(2,I)*CO(2,I)
      ZZ=CO(3,I)*CO(3,I)
      SXX=SXX + WT(IATOM(I))*(YY + ZZ)
      SYY=SYY + WT(IATOM(I))*(XX + ZZ)
      SZZ=SZZ + WT(IATOM(I))*(XX + YY)
      SXY=SXY + WT(IATOM(I))*CO(1,I)*CO(2,I)
      SXZ=SXZ + WT(IATOM(I))*CO(1,I)*CO(3,I)
      SYZ=SYZ + WT(IATOM(I))*CO(2,I)*CO(3,I)
   10 CONTINUE
      AIN(1,1)= SXX
      AIN(1,2)=-SXY
      AIN(1,3)=-SXZ
      AIN(2,1)= AIN(1,2)
      AIN(2,2)= SYY
      AIN(2,3)=-SYZ
      AIN(3,1)= AIN(1,3)
      AIN(3,2)= AIN(2,3)
      AIN(3,3)= SZZ
C
C CALCULATE PRINCIPAL AXES OF INERTIA -----
C EV ARE PRINCIPLE MOMENTS OF INERTIA -----
C
      CALL JACOBI(AIN,S,3,3)
C
C ARRANGE SMALLEST TO LARGEST -----
C
      DO 30 I=1,2
      I1=I + 1
      DO 30 K=I1,3
      IF(AIN(I,I)-AIN(K,K))30,30,20
   20 RET=AIN(I,I)
      AIN(I,I)=AIN(K,K)
      AIN(K,K)=RET
      DO 25 L=1,3
      D(L)=S(L,I)
      S(L,I)=S(L,K)
   25 S(L,K)=D(L)
   30 CONTINUE
C
      DO 35 I=1,IUNIQ
      XT=S(1,1)*CO(1,I) + S(2,1)*CO(2,I) + S(3,1)*CO(3,I)
      YT=S(1,2)*CO(1,I) + S(2,2)*CO(2,I) + S(3,2)*CO(3,I)
      ZT=S(1,3)*CO(1,I) + S(2,3)*CO(2,I) + S(3,3)*CO(3,I)
      CO(1,I)=XT
      CO(2,I)=YT
   35 CO(3,I)=ZT
C
C MOMENT OF INERTIA WITH THE PRINCIPAL AXES
C      UNIT = 10**(-39) GM * CM **2 -----
C
C      DATA CONST/6.0228D0/
C
C     XPI=AIN(1,1)/CONST
C     YPI=AIN(2,2)/CONST
C     ZPI=AIN(3,3)/CONST
C
C
C CALCULATES THE  COMPONENTS OF THE QUADROPOLE
C MOMENT ALONG THE PRINCIPAL ROTATIONAL AXES.
C
      QXX=0.D0
      QYY=0.D0
      QZZ=0.D0
      DO 7 I=1,IUNIQ
      X2 = CO(1,I)*CO(1,I)
      Y2 = CO(2,I)*CO(2,I)
      Z2 = CO(3,I)*CO(3,I)
      RQ = X2 + Y2 + Z2
      QXX = QXX+QKPT(I)* ( 3.D0*CO(1,I)*CO(1,I) - RQ )
      QYY = QYY+QKPT(I)* ( 3.D0*CO(2,I)*CO(2,I) - RQ )
    7 QZZ = QZZ+QKPT(I)* ( 3.D0*CO(3,I)*CO(3,I) - RQ )
C
C CONVERT THE QUADRUPOLE TO BUCKINGHAM ----
C
      QUCAL(1) = QXX* 1.34504D0 *0.5D0
      QUCAL(2) = QYY* 1.34504D0 *0.5D0
      QUCAL(3) = QZZ* 1.34504D0 *0.5D0
      RETURN
      END
C
C
      SUBROUTINE JACOBI(HAM,X,NDIM,N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION HAM(NDIM,NDIM),X(NDIM,NDIM)
      DO 100 I=1,N
      DO 100 J=1,N
  100 X(J,I)=0.D0
      DO 110 I=1,N
  110 X(I,I)=1.D0
 1000 HMAX=0.D0
      DO 200 I=2,N
      DO 200 J=1,I-1
  200 HMAX=DMAX1(HMAX,DABS(HAM(J,I)))
      IF(HMAX.LT.1.D-10) RETURN
      THRESH=0.1D0*HMAX
      DO 2000 I=2,N
      DO 2000 J=1,I-1
      IF(DABS(HAM(J,I)).LT.THRESH) GO TO 2000
      TAN2=0.5D0*(HAM(I,I)-HAM(J,J))/HAM(J,I)
      IF(TAN2.LT.0.D0) THEN
      TANT=TAN2-DSQRT(1+TAN2*TAN2)
                       ELSE
      TANT=TAN2+DSQRT(1+TAN2*TAN2)
                       ENDIF
      COST=1.D0/DSQRT(1+TANT*TANT)
      SINT=TANT*COST
      DIF=(HAM(I,I)-HAM(J,J))*(COST*COST-SINT*SINT)
     *     -4*HAM(J,I)*SINT*COST
      DO 2100 K=1,N
      IF(K.EQ.J.OR.K.EQ.I) GO TO 2100
      HAMKJ=HAM(K,J)
      HAM(K,J)= HAMKJ*COST+HAM(K,I)*SINT
      HAM(K,I)=-HAMKJ*SINT+HAM(K,I)*COST
      HAM(J,K)= HAM(K,J)
      HAM(I,K)= HAM(K,I)
 2100 CONTINUE
      HAMII=HAM(I,I)
      HAMJJ=HAM(J,J)
      HAM(J,J)=   HAMJJ*COST*COST + 2*HAM(J,I)*SINT*COST
     *          + HAMII*SINT*SINT 
      HAM(I,I)=   HAMJJ*SINT*SINT - 2*HAM(J,I)*SINT*COST
     *          + HAMII*COST*COST
      HAM(J,I)=HAM(J,I)*(COST*COST-SINT*SINT)
     *                         + (HAMII-HAMJJ)*SINT*COST
      HAM(I,J)=HAM(J,I)
      DO 2200 K=1,N
      XKI=X(K,I)
      X(K,I)=XKI*COST-X(K,J)*SINT
 2200 X(K,J)=XKI*SINT+X(K,J)*COST
 2000 CONTINUE
      GO TO 1000
      END
      
C
C
      SUBROUTINE OPTSSQ(X,F,XJAC,XJTJ,WORK,W1,W2,W3,W4,W5,F1,F2)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (I200=200)
      COMMON/CGSCF3/JGSCF3,ITG,CUTG,SFTG,ESFACT,IESTAT
      COMMON/CHCOR/IATEQ(I200,4),JATEQ(I200,4),NEQ(5),IUNIQ,NVAR,IFIX
      REAL*8 PARM(4)/1.2D0,2.D0,120.D0,1.D-4/
      DATA AX/0.1D0/,MAXFN/2000/
      DIMENSION X(1),F(1),XJAC(IESTAT,1),XJTJ(1),WORK(1)
      DIMENSION W1(1),W2(1),W3(1),W4(1),W5(1),F1(1),F2(1)
C
C     ----- MACHINE DEPENDENT CONSTANTS -----
C
      DATA IEPS/9/,IDELTA/5/,NSIG/6/,KSIG/16/
      EPS    = 10.D0**(-IEPS)
      DELTA  = 10.D0**(-IDELTA)
      PREC   = 10.D0**(-KSIG-1)
      REL    = 10.D0**(-KSIG/2)
      RELCON = 10.D0**(-NSIG)
C
      IF(NVAR.LE.0.OR.PARM(2).LE.1.D0.OR.PARM(1).LE.0.D0) GO TO 9305
C
      IER = 0
      AL     = PARM(1)
      ERL2 = 1.D10
      IBAD = -99
      ISW = 1
      INFER = 0
      IER = 0
      CALL SETBCK(X)
      CALL SUMSQ (X,F1)
      IEVAL = 1
      SSQOLD = 0.D0
      DO 975 I=1,IESTAT
  975 SSQOLD = SSQOLD+F1(I)*F1(I)
      WRITE(6,3) IEVAL,SSQOLD,(X(I),I=1,MIN0(NVAR,8))
    3 FORMAT(//'  IEVAL',7X,'SUM.SQ',8X,'NORM OF GRAD.',6X,'VARIABLES'
     *       /I5,F17.7,22X,8F11.5)
                                                     GO TO 1100
C
C********** MAIN ITERATION LOOP **********
C
 1000 SSQOLD = SSQ
      IF(INFER.GT.0.OR.IJAC.GE.NVAR.OR.ICOUNT.GT.0) GO TO 1100
C
C     ----- RANK ONE UPDATE TO JACOBIAN -----
C
      IJAC = IJAC+1
      DSQ = 0.D0
      DO 35 J=1,NVAR
   35 DSQ = DSQ+W2(J)*W2(J)
      IF(DSQ.LE.0.D0) GO TO 1100
      DO 50 I=1,IESTAT
      G = F(I)-F2(I)
      DO 40 J=1,NVAR
   40 G = G+XJAC(I,J)*W2(J)
      G = G/DSQ
      DO 45 J=1,NVAR
   45 XJAC(I,J) = XJAC(I,J)-G*W2(J)
   50 CONTINUE
                                                     GO TO 1200
C
C     ----- JACOBIAN BY INCREMENTING X -----
C
 1100 IJAC = 0
      DO 75 J=1,NVAR
      XDABS = DABS(X(J))
      HH = REL*(DMAX1(XDABS,AX))
      XHOLD = X(J)
      X(J) = X(J)+HH
      CALL SETBCK(X)
      CALL SUMSQ (X,F1)
      IEVAL = IEVAL+1
      X(J) = XHOLD
      IF(ISW.EQ.1) GO TO 65
C
C     ----- CENTRAL DIFFERENCES -----
C
      X(J) = XHOLD-HH
      CALL SETBCK(X)
      CALL SUMSQ (X,F2)
      IEVAL = IEVAL+1
      X(J) = XHOLD
      RHH = 0.5D0/HH
      DO 60 I=1,IESTAT
   60 XJAC(I,J) = (F1(I)-F2(I))*RHH
      GO TO 75
C
C     ----- FORWARD DIFFERENCES -----
C
   65 RHH = 1.D0/HH
      DO 70 I=1,IESTAT
   70 XJAC(I,J) = (F1(I)-F(I))*RHH
   75 CONTINUE
C
C     ----- CALCULATE GRADIENT -----
C
 1200 ERL2X = ERL2
      ERL2 = 0.D0
      DO 90 J=1,NVAR
      SUM = 0.D0
      DO 91 I=1,IESTAT
   91 SUM = SUM+XJAC(I,J)*F(I)
      W1(J) = SUM
   90 ERL2 = ERL2+SUM*SUM
      ERL2 = DSQRT(ERL2)
C
C    -----  CONVERGENCE TEST FOR NORM OF GRADIENT -----
C
      IF(IJAC  .GT.  0    ) GO TO 95
      IF(ERL2*2.LE.DELTA  ) INFER = INFER+4
      IF(ERL2*2.LE.PARM(4)) ISW = 2
C
C     ----- CALCULATE THE LOWER SUPER TRIANGE OF
C           JACOBIAN (TRANSPOSED) * JACOBIAN -----
C
   95 L=0
      DO 100 I=1,NVAR
      DO 100 J=1,I
      L=L+1
      SUM = 0.D0
      DO 110 K=1,IESTAT
  110 SUM = SUM+XJAC(K,I)*XJAC(K,J)
  100 XJTJ(L) = SUM
C
C     ----- CONVERGENCE CHECKS -----
C
      IF(INFER.GT.  0  )                                  GO TO 9000
      IF(IEVAL.GE.MAXFN)                                  GO TO 9290
C
C     ----- COMPUTE SCALING VECTOR -----
C
      K=0
      DO 115 J=1,NVAR
      K=K+J
  115 W3(J) = XJTJ(K)
C
C     ----- ADD L-M FACTOR TO DIAGONAL -----
C
      ICOUNT = 0
 1300 K = 0
      DO 150 I=1,NVAR
      DO 145 J=1,I
      K = K+1
  145 WORK(K) = XJTJ(K)
      WORK(K) = WORK(K)+W3(I)*AL
  150 W2(I) = W1(I)
C
C     ----- CHOLESKY DECOMPOSITION -----
C
C           (1) DECOMPOSE A -----
C
      CALL LUDECP (WORK,WORK,NVAR,G,XHOLD,IER)
      IF(IER .EQ.0) GO TO 160
      IER = 0
      IF(IJAC.GT.0)                                        GO TO 1100
      IF(IBAD.LE.0)                                        GO TO 1900
      IF(IBAD.GE.2)                                        GO TO 9310
                                                           GO TO 1400
C
C           (2) ELIMINATION -----
C
  160 CALL LUELMP (WORK,W2,NVAR,W2)
      IF(IBAD.NE.-99) IBAD = 0
C
C     ----- CALCULATE SUM OF SQUARES -----
C
      DO 170 J=1,NVAR
  170 W4(J) = X(J)-W2(J)
      CALL SETBCK(W4)
      CALL SUMSQ (W4,F1)
      IEVAL = IEVAL+1
      SSQ = 0.D0
      DO 175 I=1,IESTAT
  175 SSQ = SSQ+F1(I)*F1(I)
      IF(SSQ.LE.SSQOLD*1.000001D0) GO TO 205
C
C     ----- INCREASE PARAMETER AND TRY AGAIN -----
C
 1400 ICOUNT = ICOUNT+1
      AL = AL*PARM(2)**2
      IF(IJAC  .EQ.0                 ) GO TO  195
      IF(ICOUNT.GE.4.OR.AL.GT.PARM(3)) GO TO  200
  195 IF(               AL.LE.PARM(3))                     GO TO 1300
      IF(IBAD  .EQ.1                 )                     GO TO 9310
                                                           GO TO 9300
  200 AL = AL/PARM(2)**8
                                                           GO TO 1100
C
C     ----- ADJUST MARQUARDT PARAMETER -----
C
  205 IF(ICOUNT.EQ.0) AL = AL/PARM(2)
      IF(ERL2X.LE.0.D0) GO TO 210
      G = ERL2/ERL2X
      IF(ERL2.LT.ERL2X) AL = AL*DMAX1(1.D0/PARM(2),G)
      IF(ERL2.GT.ERL2X) AL = AL*DMIN1(PARM(2),G)
  210 AL = DMAX1(AL,PREC)
C
C************ ONE ITERATION CYCLE COMPLETED ************
C
      DO 220 J=1,NVAR
  220 X(J) = W4(J)
      DO 225 I=1,IESTAT
      F2(I) = F (I)
  225 F (I) = F1(I)
C
      WRITE(6,1) IEVAL,SSQ,ERL2*2,(X(I),I=1,MIN0(NVAR,8))
    1 FORMAT(I5,F17.7,F19.6,3X,8F11.5)
C
C     ----- RELATIVE CONVERGENCE TEST FOR X -----
C
      IF(ICOUNT.GT.0.OR.IJAC.GT.0)                        GO TO 1000
      XDIFMX=0.D0
      DO 230 J=1,NVAR
      XDIF = DABS(W2(J))/DMAX1(DABS(X(J)),AX)
      XDIFMX=DMAX1(XDIFMX,XDIF)
      IF(XDIF.GT.RELCON) GO TO 235
  230 CONTINUE
      INFER = 1
C
C     ----- RELATIVE CONVERGENCE TEST FOR SSQ -----
C
  235 SQDIF = DABS(SSQ-SSQOLD)/DMAX1(SSQOLD,AX)
      IF(SQDIF.LE.EPS) INFER = INFER+2
                                                           GO TO 1000
C
C     ----- SINGULAR DECOMPOSITION -----
C
 1900 IF(IBAD) 1910,1920,1930
C
C     ----- CHECK TO SEE IF CURRENT ITERATE HAS CYCLED BACK TO
C           THE LAST SINGULAR POINT -----
C
 1920 DO 250 J=1,NVAR
      XHOLD = W5(J)
      IF(DABS(X(J)-XHOLD).GT.RELCON*DMAX1(AX,DABS(XHOLD)))GO TO 1910
  250 CONTINUE
                                                           GO TO 9295
C
C     ----- UPDATE THE BAD X VALUES -----
C
 1910 DO 260 J=1,NVAR
  260 W5(J) = X(J)
      IBAD = 1
C
C     ----- REPLACE ZEROES ON HESSIAN DIAGONAL -----
C
 1930 IZRO = 0
      DO 285 J=1,NVAR
      IF(W3(J).GT.0.D0) GO TO 285
      IZRO = IZRO+1
      W3(J) = 1.D0
  285 CONTINUE
      IF(IZRO.LT.NVAR)                                     GO TO 1300
                                                           IER = 38
                                                           GO TO 9000
C
C     ----- ERROR TERMINATION -----
C
 9290 IER = IER+1
 9295 IER = IER+1
 9300 IER = IER+1
 9305 IER = IER+1
 9310 IER = IER+129
 9000 CONTINUE
           IF(IER.EQ.129) THEN
                                 WRITE(6, 8)
      ELSE IF(IER.EQ.130) THEN
                                 WRITE(6,18)
                                 RETURN
      ELSE IF(IER.EQ.131) THEN
                                 WRITE(6,28)
      ELSE IF(IER.EQ.132) THEN
                                 WRITE(6,38)
      ELSE IF(IER.NE. 0 ) THEN
                                 WRITE(6,48) IER
                          ENDIF
      WRITE(6,2) NSIG,XDIFMX,IEPS,SQDIF,IDELTA,ERL2*2
C
    2 FORMAT(//28X,'*** RESULTS OF LEVENBERG-MARQUARDT ALGORITHM ***',
     *//22X,'CONVERGENCE LIMIT FOR THE VARIABLES =',' 1.D-0',I1,E15.3
     * /22X,'CONVERGENCE LIMIT FOR THE SSQ       =',' 1.D-0',I1,E15.3
     * /22X,'CONVERGENCE LIMIT FOR THE NORM OF G =',' 1.D-0',I1,E15.3)
    8 FORMAT(/10X,'FAILED TO RECOVER FROM SINGULAR JACOBIAN')
   18 FORMAT(/10X,'ONE OF M,N,PARM,IS INCORRECT')
   28 FORMAT(/10X,'MARQUARDT PARAMETER EXCEEDED')
   38 FORMAT(/10X,'THE SOLUTION CYCLED BACKTO SINGULARITY')
   48 FORMAT(/10X,'ERROR INDICATING PARAMETER      =',I5)
      RETURN
      END
C
C
      SUBROUTINE LUELMP (A,B,N,X)
      IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION A(1),B(1),X(1)
C
C     ----- SOLUTION OF LY = B -----
C
      IP=1
      IW = 0
      DO 15 I=1,N
      T=B(I)
      IM1 = I-1
      IF(IW .EQ. 0) GO TO 9
      IP=IP+IW-1
      DO 5 K=IW,IM1
      T = T-A(IP)*X(K)
      IP=IP+1
    5 CONTINUE
      GO TO 10
    9 IF(T .NE. 0.D0) IW = I
      IP = IP+IM1
   10 X(I)=T*A(IP)
      IP=IP+1
   15 CONTINUE
C
C     ----- SOLUTION OF UX = Y -----
C
      N1 = N+1
      DO 30 I = 1,N
      II = N1-I
      IP=IP-1
      IS=IP
      IQ=II+1
      T=X(II)
      IF(N.LT.IQ) GO TO 25
      KK = N
      DO 20 K=IQ,N
      T = T - A(IS) * X(KK)
      KK = KK-1
      IS = IS-KK
   20 CONTINUE
   25 X(II)=T*A(IS)
   30 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE LUDECP (A,UL,N,D1,D2,IER)
      IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION A(1),UL(1)
C
      kount=0
      D1=1.D0
      D2=0.D0
      RN = 1.D0/(FLOAT(N)*16.D0)
      IP = 1
      IER=0
      DO 45 I = 1,N
      IQ = IP
      IR = 1
      DO 40 J = 1,I
      X = A(IP)
      IF(J .EQ. 1) GO TO 10
      DO 5  K=IQ,IP1
      X = X - UL(K) * UL(IR)
      IR = IR+1
    5 CONTINUE
   10 IF(I.NE.J) GO TO 30
      D1 = D1*X
      IF(A(IP) + X*RN .LE. A(IP)) GO TO 50
   15 IF(DABS(D1).LE.1.D0) GO TO 20
      D1 = D1 * 0.0625D0
      D2 = D2 + 4.D0
      GO TO 15
   20 IF(DABS(D1) .GE. 0.0625D0) GO TO 25
      D1 = D1 * 16.D0
      D2 = D2 - 4.D0
      GO TO 20
   25 UL(IP) = 1.D0/DSQRT(X)
      GO TO 35
   30 UL(IP) = X * UL(IR)
   35 IP1 = IP
      IP = IP+1
      IR = IR+1
   40 CONTINUE
   45 CONTINUE
      GO TO 60
   50 IER = 129
   60 CONTINUE
      RETURN
      END
      SUBROUTINE PREPOL
C%SN
C%SN  preparation for bond polarizability fit
C%SN
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (I200=200,I340=340)
      COMMON/CHCOR/IATEQ(I200,4),JATEQ(I200,4),NEQ(5),IUNIQ,NVAR,IFIX
     * ,IATOM(I200),IVARY(I200,4),XKPT(I200),YKPT(I200),ZKPT(I200)
     * ,QKPT (I200)
      COMMON/CGSCF3/JGSCF3,ITG,CUTG,SFTG,ESFACT,IESTAT
      COMMON/COMP12/CXNP(I340),CYNP(I340),CZNP(I340),
     & QUPOT(I340,I340)
      COMMON/PDPOL/IBON(I200),JBON(I200),IEQBON(I200),IBPAIR,IPTYP,
     & FPCHRG,PKPT(I200*2),NVARP,IVARP(I200),IESTAH,IATID,IATM(I200),
     & IEQATM(I200),IVARA(I200),NVARA,NVART,NVARAN,VECISO,IHTYP,NVARH,
     & ALPTRN(I200),ALPANI(I200),ALPATM(I200)
      COMMON/PRPOL/ XMID(I200),YMID(I200),ZMID(I200),
     &          TRNX1(I340,I200),TRNY1(I340,I200),TRNZ1(I340,I200),
     &          TRNX2(I340,I200),TRNY2(I340,I200),TRNZ2(I340,I200),
     &          BONX(I200),BONY(I200),BONZ(I200),BONLEN(I200),
     &          FLDLEN(I340,I200),ANIFAC(I340,I200),PALEN(I340,I200),
     &          SIGNA(I340,I200),FTMP(I340),
     &          ATMX1(I340,I200),ATMY1(I340,I200),ATMZ1(I340,I200),
     &          ATMX2(I340,I200),ATMY2(I340,I200),ATMZ2(I340,I200) 
C           
C     midpoint of bond   XMID,YMID,ZMID
C  
      DO I=1,IBPAIR
        XMID(I)=0.5D0*(XKPT(IBON(I))+XKPT(JBON(I)))
        YMID(I)=0.5D0*(YKPT(IBON(I))+YKPT(JBON(I)))
        ZMID(I)=0.5D0*(ZKPT(IBON(I))+ZKPT(JBON(I)))
      ENDDO
C  
C     electrostatic potential for fixed parameters
C  
      DO I=1,IESTAH
        SUMP=0.D0
        DO J=1,IUNIQ
          XC     =CXNP(I)-XKPT(J)
          YC     =CYNP(I)-YKPT(J)
          ZC     =CZNP(I)-ZKPT(J)
          PALEN(I,J)=DSQRT(XC*XC+YC*YC+ZC*ZC)
          SUMP   =SUMP-QKPT(J)/PALEN(I,J)
          FLDL2=PALEN(I,J)*2.D0/VECISO
          ATMX1(I,J)=XKPT(J)-XC/FLDL2
          ATMX2(I,J)=XKPT(J)+XC/FLDL2
          ATMY1(I,J)=YKPT(J)-YC/FLDL2
          ATMY2(I,J)=YKPT(J)+YC/FLDL2
          ATMZ1(I,J)=ZKPT(J)-ZC/FLDL2
          ATMZ2(I,J)=ZKPT(J)+ZC/FLDL2
C         write(6,*) 'ATM XYZ'
C         write(6,'(3f10.5)') ATMX1(I,J),ATMY1(I,J),ATMZ1(I,J)
C         write(6,'(3f10.5)') ATMX2(I,J),ATMY2(I,J),ATMZ2(I,J)
        ENDDO
        FTMP(I)   = -SUMP
      ENDDO
C  
C     bond length 
C  
      DO I=1,IBPAIR
        BONX(I)=XKPT(IBON(I))-XKPT(JBON(I))
        BONY(I)=YKPT(IBON(I))-YKPT(JBON(I))
        BONZ(I)=ZKPT(IBON(I))-ZKPT(JBON(I))
        BONLEN(I)=DSQRT(BONX(I)**2+BONY(I)**2+BONZ(I)**2)
      ENDDO
C%SN
C%SN    electric field length & anisotoropy factor & 
C%SN    x y z coordinates of transverse induced dipoles
C%SN
      DO I=1,IESTAH
        DO J=1,IBPAIR
          FLDX=CXNP(I)-XMID(J)
          FLDY=CYNP(I)-YMID(J)
          FLDZ=CZNP(I)-ZMID(J)
          FLDLEN(I,J)=DSQRT(FLDX**2+FLDY**2+FLDZ**2)
          IF(PALEN(I,IBON(J)).GE.PALEN(I,JBON(J)))
     &                                               THEN
            ANIFAC(I,J)=(-BONX(J)*FLDX-BONY(J)*FLDY
     &               -BONZ(J)*FLDZ)/(FLDLEN(I,J)*BONLEN(J))
            SIGNA(I,J)=1.0D0
          ELSE
            ANIFAC(I,J)=(BONX(J)*FLDX+BONY(J)*FLDY
     &               +BONZ(J)*FLDZ)/(FLDLEN(I,J)*BONLEN(J))
            SIGNA(I,J)=-1.0D0
          ENDIF
          FLDL2=FLDLEN(I,J)*2.D0/VECISO
          TRNX1(I,J)=XMID(J)-FLDX/FLDL2
          TRNX2(I,J)=XMID(J)+FLDX/FLDL2
          TRNY1(I,J)=YMID(J)-FLDY/FLDL2
          TRNY2(I,J)=YMID(J)+FLDY/FLDL2
          TRNZ1(I,J)=ZMID(J)-FLDZ/FLDL2
          TRNZ2(I,J)=ZMID(J)+FLDZ/FLDL2
C         write(6,*) 'TRN XYZ'
C         write(6,'(3f10.5)') TRNX1(I,J),TRNY1(I,J),TRNZ1(I,J)
C         write(6,'(3f10.5)') TRNX2(I,J),TRNY2(I,J),TRNZ2(I,J)
        ENDDO
      ENDDO
      RETURN
      END
