      SUBROUTINE POLAR
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
C**********************************************************************
C
C   POLAR SETS UP THE CALCULATION OF THE MOLECULAR ELECTRIC RESPONSE
C   PROPERTIES BY FFHPOL.
C
C**********************************************************************
      character*2 ELEMNT
      common /TITLES/ KOMENT,TITLE
      common /POLVOL/ POLVOL(107)
      common /KEYWRD/ KEYWRD
      common /GEOKST/ NATOMS,LABELS(NUMATM),
     1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
      common /GEOVAR/ NVAR,LOC(2,MAXPAR),IDUMY,XPARAM(MAXPAR)
      common /SCFTYP/ EMIN,LIMSCF
      common /TIME  / TIME0
      common /ELEMTS/ ELEMNT(107)
      common /CORE  / CORE(107)
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      common /GEOSYM/ NDEP, LOCPAR(MAXPAR), IDEPFN(MAXPAR),
     1                    LOCDEP(MAXPAR)
      common /GEOM  / GEO(3,NUMATM)
      common /LAST  / LAST
      common /COORD / COORD(3,NUMATM)
      common /EULER / TVEC(3,3),IDTVEC
      common /SCRACH/ RXYZ(MPACK),XDUMY(MAXPAR**2-MPACK)
      dimension GRAD(MAXPAR),ROTVEC(3,3),DIPVEC(3),
     1          TEMPV(3,3)
      character  KEYWRD*241, TYPE*7, KOMENT*81, TITLE*81
      logical LET, LIMSCF
      TYPE=' MNDO  '
      LIMSCF=.FALSE.
      LET=(INDEX(KEYWRD,'LET').NE.0)
      if (INDEX(KEYWRD,'MINDO') .NE. 0) TYPE='MINDO/3'
      if (INDEX(KEYWRD,'AM1') .NE. 0)    TYPE='  AM1  '
      write(MIOB,10)
   10 format('1',20('*'),' FINITE-FIELD POLARIZABILITIES ',
     1        20('*'),//,
     2      '    THE F-F METHOD IS PERFORMED USING BOTH AN ENERGY',/,
     3      '    AND DIPOLE MOMENT EXPANSION.  THESE RESULTS ARE',/,
     4      '    LISTED BELOW AS "E4" AND "DIP", RESPECTIVELY.')
      call GMETRY(GEO,COORD)
C
C  Orient the molecule with the moments of inertia.
C  This is done to ensure a unique, reproduceable set of directions.
C  If LET is specified, the input orientation will be used.
C
      if (.NOT.LET) then
         MASS = 1.0D00
         call AXIS (COORD,NUMAT,A,B,C,SUMW,MASS,ROTVEC)
         write(MIOB,20)
   20    format(/' ROTATION MATRIX FOR ORIENTATION OF MOLECULE:'/)
         do 40 I = 1,3
            write(MIOB,30) (ROTVEC(I,J),J=1,3)
   30       format(5X,3F12.6)
   40    continue
C
C  ROTATE ATOMS
C
         do 70 I = 1,NUMAT
            do 60 J = 1,3
               SUM = 0.0D00
               do 50 K = 1,3
                  SUM = SUM + COORD(K,I)*ROTVEC(K,J)
   50          continue
               GEO(J,I) = SUM
   60       continue
   70    continue
         do 90 I = 1,NUMAT
            do 80 J = 1,3
               COORD(J,I) = GEO(J,I)
   80       continue
   90    continue
         write(MIOB,'(//10X,''CARTESIAN COORDINATES '',/)')
         write(MIOB,'(4X,''NO.'',7X,''ATOM'',9X,''X'',
     %                9X,''Y'',9X,''Z'',/)')
         L=0
         do 100 I=1,NUMAT
            if (NAT(I).EQ.99.OR.NAT(I).EQ.107) go to 100
            L=L+1
            write(MIOB,'(I6,8X,A2,4X,3F10.4)')
     %                             L,ELEMNT(NAT(I)),(COORD(J,L),J=1,3)
  100    continue
C
C  if POLYMER, ROTATE TVEC
C  (BEWARE:  THE POLYMER SECTIONS MAY NOT WORK YET)
C
         if (IDTVEC.GT.0) then
            do 130 I = 1,IDTVEC
               do 120 J = 1,3
                  SUM = 0.0D00
                  do 110 K = 1,3
                     SUM = SUM + TVEC(K,I)*ROTVEC(K,J)
  110             continue
                  TEMPV(J,I) = SUM
  120          continue
  130       continue
            do 150 I = 1,3
               do 140 J = 1,IDTVEC
                  TVEC(I,J) = TEMPV(I,J)
  140          continue
  150       continue
            write(MIOB,160) ((TVEC(J,I),J=1,3),I=1,IDTVEC)
  160       format(/' NEW TRANSLATION VECTOR:'/,
     1           ' ',3(3F15.5))
         end if
      end if
C
      LAST=1
      NA(1)=99
C
C  SET UP THE VARIABLES IN XPARAM AND LOC, THESE ARE IN CARTESIAN
C  COORDINATES.
C
      NDEP=0
      NUMAT=0
      SUMX=0.D0
      SUMY=0.D0
      SUMZ=0.D0
      do 180 I=1,NATOMS
         if ((LABELS(I).NE.99).AND.(LABELS(I).NE.107)) then
            NUMAT=NUMAT+1
            LABELS(NUMAT)=LABELS(I)
            SUMX=SUMX+COORD(1,NUMAT)
            SUMY=SUMY+COORD(2,NUMAT)
            SUMZ=SUMZ+COORD(3,NUMAT)
            do 170 J=1,3
  170       GEO(J,NUMAT)=COORD(J,NUMAT)
         end if
  180 continue
      SUMX=SUMX/NUMAT
      SUMY=SUMY/NUMAT
      SUMZ=SUMZ/NUMAT
      SUMMAX=0.D0
      ATPOL=0.D0
      do 190 I=1,NUMAT
         if (LABELS(I).NE.107) then
            ATPOL=ATPOL+POLVOL(LABELS(I))
         end if
         GEO(1,I)=GEO(1,I)-SUMX
         if (SUMMAX.LT.ABS(GEO(1,I))) SUMMAX=ABS(GEO(1,I))
         GEO(2,I)=GEO(2,I)-SUMY
         if (SUMMAX.LT.ABS(GEO(2,I))) SUMMAX=ABS(GEO(2,I))
         GEO(3,I)=GEO(3,I)-SUMZ
         if (SUMMAX.LT.ABS(GEO(3,I))) SUMMAX=ABS(GEO(3,I))
  190 continue
C
      NVAR=0
      NATOMS = NUMAT
      call COMPFG(GEO, .TRUE., HEAT0, .TRUE., GRAD, .FALSE.)
      write(MIOB,200) HEAT0
  200 format(//' ENERGY OF "REORIENTED" SYSTEM WITHOUT FIELD:',
     1        F20.10)
C
      call FFHPOL (HEAT0,ATPOL,DIPVEC)
C
      return
      end
