c	Not vectorized
c       CHELP-NET ATOMIC CHARGES FROM AB INITIO WAVE FUNCTIONS
c
c	Modified for Grid Operations by Curt Breneman, Yale University
c	Department of Chemistry, 3/88.  (Currently of Rensselaer
c	Polytechnic Institute, Troy NY 12180.)
c
c
c        CHELPG
c
c
c       (NET ATOMIC) CHARGES FIT TO ELECTROSTATIC POTENTIALS
c
c
c
c	 Original CHELP code by:
c        M.M. FRANCL
c
c        L.E. CHIRLIAN
c
c
c        OCTOBER 1985
c
c        PRINCETON CHEMISTRY DEPARTMENT VAX 11/780
c
c        VMS 3.7
c
c
c        FEBRUARY 1988
c	 MODIFIED TO USE GAUSSIAN86 CHECKPOINT FILES
c	 Modified to use G88/90 checkpoint files 1/89
c	 YALE UNIVERSITY DEPARTMENT OF CHEMISTRY
c	 WIBERG GROUP VMS 4.5
c	 CURT BRENEMAN
c
c g90 version kjb 12/18/90
c Modification due to Kenneth J. Butenhof,  CWRU Cleveland Ohio:
C FOPEN REMOVE ,JUNK PARAMETER
c /ipo/ ipo(10) changed to /ipo/ ipo(5)
c "e=e+(value*coef_)"  where coef_ is undefined, thua e was unchanged
c  the line above it is the correct and indended statement that
c changes the value of e.  THIS BAD LINE WAS REMOVED (three lies above # 2441)
c 
c add definition of chkfile to first occurence of /out/
c in /io/ make sure all have ",ipunch" in common
c
c 12/7/90 add g88 cray /b/ common  --- from  grep of utilam.f
c  this modified common block uses SHLADF(2000) and is present
c every time /b/ uses the C3(2000),CF(2000).... form.
c When c3(6000) is used in /b/ no change was made.
c    Also changed the /bb/ reference to /b/  -- note
c that although the /b/ array is loaded by reference to the
c address of EXX(1) the subroutine uep, which contains the
c only reference to /bb/ accesses EXX  only through the common
c block reference.  AT best /bb/ is forced to lie on /b/ due
c to the previous definition of EXX in /b/, the old way, at worst,
c a new EXX ... is created (with no value).
c 
c 12/11/90 change double precision to real*8 -- for CRAY implementation
c	   CHANGE ALL DSQRT AND DABS TO SQRT AND ABS
c	   CHANGE DMAX1 TO AMAX1 (NOTE MAX1 IS TYPE INTEGER)
c
c 12/17/90 remove commented out junk in READIO
C 	   add z-matric read and print in READIO
c
C 12/18/90 add fclose
c***********************************************************************
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
      integer*4 handle1
c***	Unix Fix
c      ISTAT1=LIB$INIT_TIMER(HANDLE1)
c
c***
c
c
c        READ IN DATA FROM CHECKPOINT FILE
c
c
c
c# 26 "chelpg.for"
      handle1 = 0
c
c
c
c
c        SELECT POINTS FOR FITTING, BEGIN WITH SHELL OF RADIUS 2A AND
c
c        INCREASING BY .5A SELECT POINTS FROM THE ROUGHLY RADIAL
c        DISTRIBUTION WHICH ARE NOT ENCLOSED BY THE VAN DER WAALS
c        ENVELOPE OF THE MOLECULE UNTIL A PREDETERMINED MAXIMUM NUMBER
c        OF POINTS HAS BEEN REACHED
c
c# 33 "chelpg.for"
      call readin
c
c
c        CALCULATE THE ELECTROSTATIC POTENTIAL USING FIRST ORDER
c        HARTREE-FOCK PERTURBATION THEORY
c
c# 41 "chelpg.for"
      call ball
c
c
c        USING METHOD OF LAGRANGE MULTIPLIERS, FIT BY LEAST SQUARES THE
c        CHARGES TO THE ELECTROSTATIC POTENTIAL, CONSTRAINING THE FIT TO 
c        REPRODUCE THE TOTAL MOLECULAR CHARGE
c
c# 46 "chelpg.for"
      call ep
c
c        PRINT OUT TABLE OF RESULTS
c
c
c# 52 "chelpg.for"
      call fit
c
c
c***	TRACE-7
c      ISTAT1=LIB$SHOW_TIMER(HANDLE1)
c
c***
c# 56 "chelpg.for"
      call output
c
c
c
c
c# 61 "chelpg.for"
      end
c
c
c  ROUTINE TO SELECT POINTS FOR FITTING TO THE ELECTROSTATIC POTENTIAL.
c
c
c  POINTS WHICH LIE WITHIN THE VAN DER WAALS ENVELOPE OF THE MOLECULE
c        ARE EXCLUDED.
c
c
c
c      POINTS ARE INITIALLY SELECTED IN A CUBE AROUND THE MOLECULE WHICH
c     IS SCALED TO THE SIZE OF THE MOLECULE+RMAX. THIS IS PRESENTLY AN INPUT
c	 PARAMETER.  POINTS ARE THEN EXCLUDED IF THEY FALL WITHIN THE INPUT
c	 VDW RADIUS OF ANY OF THE ATOMS, OR, IF THEY FALL OUTSIDE
c	 A DESIGNATED DISTANCE (RMAX) FROM ALL OF THE ATOMS.   THE REMAINING
c	 POINTS ARE PACKED IN A SET OF THREE (X,Y,Z) VECTORS, AND SENT TO THE
c	 LAGRANGE LEAST-SQUARES FITTING ROUTINE.  THE ORIGINAL CHELP INPUT
c	 DECK IS AUGMENTED BY ADDING TWO FREE-FORMAT VARIABLES AT THE END.
c	 THE TWO NEW INPUT VARIABLES ARE 'RMAX' AND 'DELR', WHERE RMAX
c	 IS THE MAXIMUM DISTANCE A POINT CAN BE FROM ANY ATOM AND STILL
c	 BE CONSIDERED IN THE FIT, AND DELR IS THE DISTANCE BETWEEN POINTS
c	 IN THE GRID.  BOTH RMAX AND DELR ARE IN ANGSTROMS.
c
c
c	 CURT BRENEMAN AND TERESA LEPAGE
c	 YALE UNIVERSITY DEPARTMENT OF CHEMISTRY 3/88
c
c	 ORIGINAL CODE BY:
c
c        L.E. CHIRLIAN
c
c        M.M. FRANCL
c
c        APRIL 1985
c
c
c
      subroutine ball()
c
c
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
      parameter (npoints = 50000)
c+++
      common /io/ in, iout, ipunch
c+++
      common /mol/ natoms, icharg, multip, nae, nbe, nel, nbasis, ian(
     &401), atmchg(400), c(3, 400)
      common /ipo/ ipo(5)
      common /sphere/ radii(400), ntotp
c
c
      common /points/ p(3, npoints), maxpnts
c
c***	READ IN THE THE RMAX AND DELR VALUES IN ANGSTROMS.
c
      data ang2au / 1.889726878d0 /
c# 109 "chelpg.for"
      read(unit=in, fmt=*) rmax, delr
c***
c
c        CONVERT RADII TO AU
c
c
c
c# 110 "chelpg.for"
      write(unit=iout, fmt=*) ' RMAX = ', rmax, ' (ANGS), DELR = ', delr
     &, ' (ANGS).'
c# 115 "chelpg.for"
      delr = delr * ang2au
c
c	 WHILE CONVERTING THE VDW RADII TO AU, FIND THE EXTREMA OF THE
c	 MOLECULAR GEOMETRY.
c
c# 116 "chelpg.for"
      rmax = rmax * ang2au
c# 121 "chelpg.for"
      xmax = -50.0d0
      xmin = 50.0d0
      ymax = -50.0d0
      ymin = 50.0d0
      zmax = -50.0d0
c
c# 126 "chelpg.for"
      zmin = 50.0d0
c# 128 "chelpg.for"
      write(unit=iout, fmt=*) ' THERE ARE ', natoms,
     &' ATOMS TO CONSIDER.'
c# 129 "chelpg.for"
      do 10 i = 1, natoms
c
c# 130 "chelpg.for"
      radii(i) = radii(i) * ang2au
c# 132 "chelpg.for"
      if (c(1,i) .gt. xmax) xmax = c(1,i)
      if (c(1,i) .lt. xmin) xmin = c(1,i)
      if (c(2,i) .gt. ymax) ymax = c(2,i)
      if (c(2,i) .lt. ymin) ymin = c(2,i)
      if (c(3,i) .gt. zmax) zmax = c(3,i)
      if (c(3,i) .lt. zmin) zmin = c(3,i)
c
c
c# 138 "chelpg.for"
   10 continue
c# 140 "chelpg.for"
      write(unit=iout, fmt=*) ' XMAX = ', xmax, ' (AU), XMIN = ', xmin,
     &' (AU).'
c# 141 "chelpg.for"
      write(unit=iout, fmt=*) ' YMAX = ', ymax, ' (AU), YMIN = ', ymin,
     &' (AU).'
c
c        DETERMINE THE MINIMUM CUBE DIMENSIONS REQUIRED TO CONTAIN
c	 THE MOLECULE, INCLUDING THE MAXIMUM SELECTION RADIUS (RMAX)
c	 ON BOTH SIDES.
c
c# 142 "chelpg.for"
      write(unit=iout, fmt=*) ' ZMAX = ', zmax, ' (AU), ZMIN = ', zmin,
     &' (AU).'
c# 148 "chelpg.for"
      xrange = (xmax - xmin) + (2.0d0 * rmax)
      yrange = (ymax - ymin) + (2.0d0 * rmax)
c
c
c# 150 "chelpg.for"
      zrange = (zmax - zmin) + (2.0d0 * rmax)
c# 152 "chelpg.for"
      nxpts = int(xrange / delr)
      nypts = int(yrange / delr)
c
c# 154 "chelpg.for"
      nzpts = int(zrange / delr)
c# 156 "chelpg.for"
      write(unit=iout, fmt=*) ' NUMBER OF X POINTS REQUIRED = ', nxpts
      write(unit=iout, fmt=*) ' NUMBER OF Y POINTS REQUIRED = ', nypts
      write(unit=iout, fmt=*) ' NUMBER OF Z POINTS REQUIRED = ', nzpts
      maxposs = (nxpts * nypts) * nzpts
c
c
c	 RESET POINT COUNTER FOR NUMBER OF SELECTED POINTS
c
c# 160 "chelpg.for"
      write(unit=iout, fmt=*) ' TOTAL NUMBER OF POINTS CONSIDERED = ',
     &maxposs
c
c       LOOP OVER POSSIBLE POINTS
c
c# 165 "chelpg.for"
      ipoint = 0
c
c# 169 "chelpg.for"
      do 200 ii = 1, nxpts + 1
c
c# 171 "chelpg.for"
      p1 = (xmin - rmax) + (dble(ii - 1) * delr)
c
c# 173 "chelpg.for"
      do 200 jj = 1, nypts + 1
c
c# 175 "chelpg.for"
      p2 = (ymin - rmax) + (dble(jj - 1) * delr)
c
c# 177 "chelpg.for"
      do 200 kk = 1, nzpts + 1
c
c
c
c        IS THIS POINT WITHIN A VAN DER WAALS SPHERE OR OUTSIDE THE
c	 RMAX DISTANCE FROM ALL ATOMS?
c
c
c# 179 "chelpg.for"
      p3 = (zmin - rmax) + (dble(kk - 1) * delr)
c# 185 "chelpg.for"
      radmin = 50.0d0
      do 100 i = 1, natoms
      vrad = radii(i)
      dist = (((p1 - c(1,i)) ** 2) + ((p2 - c(2,i)) ** 2)) + ((p3 - c(3,
     &i)) ** 2)
c# 189 "chelpg.for"
      dist = sqrt(dist)
      if (dist .lt. vrad) goto 210
      if (dist .lt. radmin) radmin = dist
  100 continue
c
c
c        STORE POINTS (IN ATOMIC UNITS)
c
c
c
c# 193 "chelpg.for"
      if (radmin .gt. rmax) goto 210
c# 197 "chelpg.for"
      ipoint = ipoint + 1
      p(1,ipoint) = p1
      p(2,ipoint) = p2
      p(3,ipoint) = p3
      if (ipo(2) .eq. 1) write(unit=iout, fmt=*) 'POINT ', ipoint,
     &' X,Y,Z ', p1, p2, p3
  210 continue
c
c
c# 204 "chelpg.for"
  200 continue
c# 206 "chelpg.for"
      maxpnts = ipoint
      write(unit=iout, fmt=*)
     &' NUMBER OF POINTS SELECTED FOR FITTING : ', maxpnts
c# 208 "chelpg.for"
      return
c
c
c# 209 "chelpg.for"
      end
c
c
c        ROUTINE TO CALCULATE THE ELECTROSTATIC POTENTIAL FROM FIRST ORD
cER
c        PERTURBATION THEORY
c
c
c
c        M.M. FRANCL    APRIL 1985
c
c        MODIFIED VERSION OF A MEPHISTO ROUTINE
c
c        RESTRICTED TO CLOSED SHELL MOLECULES
c
c
c
      subroutine ep()
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
      parameter (npoints = 50000)
c
c
      integer*4 shella, shelln, shellt, aos, shellc, aon, handle
      character chkfil*40
      common /io/ in, iout, ipunch
c+++
      common /ipo/ ipo(5)
c
c===	Gaussian88 Modification for enlarged common /b/.
      common /mol/ natoms, icharg, multip, nae, nbe, nel, nbasis, ian(
     &401), atmchg(400), c(3, 400)
c==== Old G86 Version of common /b/
c      COMMON/B/EXX(1200),C1(1200),C2(1200),C3(1200),
c     $         X(400),Y(400),Z(400),JAN(400),SHELLA(400),SHELLN(400),
c     $         SHELLT(400),SHELLC(400),AOS(400),AON(400),NSHELL,MAXTYP
c+++
c      COMMON /B/ EXX(240),C1(240),C2(240),C3(240),X(80),Y(80),Z(80),
c
c     $           JAN(80),SHELLA(80),SHELLN(80),SHELLT(80),SHELLC(80)
c
c     $          ,AOS(80),AON(80),NSHELL,MAXTYP
c
c      COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS,IAN(101),
c
c     $             ATMCHG(100),C(3,100)
c
      common /b/ exx(6000), c1(6000), c2(6000), c3(6000), x(2000), y(
     &2000), z(2000), jan(2000), shella(2000), shelln(2000), shellt(2000
     &), shellc(2000), aos(2000), aon(2000), nshell, maxtyp
      common /points/ p(3, npoints), maxpnts
      common /elp/ elecp(npoints)
      common /charge/ coef_alpha(100000), coef_beta(100000), iuhf
c
c
      common /out/ ntitle(20, 3), chkfil, q(400), i6to5, rms, percent,
     &nlin, nend(3)
c
c
      dimension hpert(100000), index(1280)
c        DIVERT TO ROUTINE UEP IF WAVEFUNCTION IS UNRESTRICTED
c
c        HARTREE-FOCK WAVEFUNCTION
c
c
c
      data iptchg / 1.0 /
      data zero / 0.0 /
      data two / 2.0 /
      data vnucmax / 30.0 /
c# 257 "chelpg.for"
      if (iuhf .eq. 1) then
      call uep
      return
c
c
c# 260 "chelpg.for"
      end if
c
c
c        SET UP THE INDEXING TABLE FOR HPERT
c
c
c
c# 262 "chelpg.for"
      handle = 0
c# 266 "chelpg.for"
      do 100 i = 1, nbasis
      index(i) = ((i - 1) * i) / 2
c
c
c        BEGIN LOOP TO CALCULATE ELECTROSTATIC POTENTIAL
c
c
c
c# 268 "chelpg.for"
  100 continue
c# 272 "chelpg.for"
      nocc = nel / 2
c
c
c        START OF LOOP
c
c
c
c# 273 "chelpg.for"
      mvir = nocc + 1
c# 277 "chelpg.for"
      do 200 npnt = 1, maxpnts
      x1 = p(1,npnt)
      x2 = p(2,npnt)
c
c
c     CALCULATE THE ONE-ELECTRON INTEGRALS
c
c
c
c# 280 "chelpg.for"
      x3 = p(3,npnt)
c# 284 "chelpg.for"
      if (ipo(5) .eq. 1) then
      write(unit=iout, fmt=3010)
c***
c      ISTAT = LIB$INIT_TIMER(HANDLE)
c
c***
c# 286 "chelpg.for"
 3010 format(1x,18hTIME FOR INTEGRALS)
c
c
      end if
c
c
c***
c      IF (IPO(5).EQ.1) ISTAT = LIB$SHOW_TIMER(HANDLE)
c
c***
c
c
c# 292 "chelpg.for"
      call intgrl(hpert, x1, x2, x3, iptchg, i6to5)
c
c
c# 298 "chelpg.for"
      if (ipo(4) .eq. 1) call linout(hpert, nbasis, 0, 0)
c# 300 "chelpg.for"
      if (ipo(5) .eq. 1) then
      write(unit=iout, fmt=3000)
c***
c      ISTAT = LIB$INIT_TIMER(HANDLE)
c
c***
c# 302 "chelpg.for"
 3000 format(1x,18hTIME FOR TRANSFORM)
c
c
c     FORM THE HPERT MATRIX ELEMENTS
c
c
c
c# 306 "chelpg.for"
      end if
c# 310 "chelpg.for"
      e = zero
c
c
c        SUM OVER OCCUPIED MOS
c
c
c
c# 311 "chelpg.for"
      icoefi = - nbasis
c# 315 "chelpg.for"
      do 220 ii = 1, nocc
c
c
c        CALCULATE ELECTROSTATIC POTENTIAL
c
c
c
c# 316 "chelpg.for"
      icoefi = icoefi + nbasis
c# 320 "chelpg.for"
      do 221 ip = 1, nbasis
      cpi = coef_alpha(icoefi + ip)
c
c
c# 322 "chelpg.for"
      ipdex = index(ip)
c# 324 "chelpg.for"
      do 222 iq = 1, ip
      e = e + ((cpi * coef_alpha(icoefi + iq)) * hpert(ipdex + iq))
  222 continue
      do 223 iq = ip + 1, nbasis
      e = e + ((cpi * coef_alpha(icoefi + iq)) * hpert(ip + index(iq)))
c
c
c# 329 "chelpg.for"
  223 continue
c# 331 "chelpg.for"
  221 continue
c
c
c***
c      IF (IPO(5).EQ.1) ISTAT = LIB$SHOW_TIMER(HANDLE)
c
c***
c
c
c        CALCULATE NUCLEAR PART OF ELECTROSTATIC POTENTIAL
c
c
c
c# 332 "chelpg.for"
  220 continue
c# 340 "chelpg.for"
      vnuc = zero
      do 300 iatom = 1, natoms
      del1 = c(1,iatom) - x1
      del2 = c(2,iatom) - x2
      del3 = c(3,iatom) - x3
      ra = sqrt(((del1 * del1) + (del2 * del2)) + (del3 * del3))
      if (ra .eq. zero) then
      vnuc = vnucmax
      goto 310
      end if
      vnuc = vnuc + (ian(iatom) / ra)
  300 continue
c
c
c# 352 "chelpg.for"
  310 continue
c# 354 "chelpg.for"
      elecp(npnt) = (e * two) + (vnuc * iptchg)
      if (ipo(5) .eq. 1) write(unit=iout, fmt=*) 'E(', npnt, ') = ', e
  200 continue
      return
      end
c
c
c        ROUTINE TO USE METHOD OF LAGRANGE MULTIPLIERS TO OBTAIN BEST
c
c        LEAST SQUARE FIT WITH CONSTRAINTS
c
c
c
c        M.M. FRANCL
c
c        APRIL 1985
c
c
c
      subroutine fit()
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
      parameter (npoints = 50000)
      integer*4 which1
c
c
      character chkfil*40
      common /io/ in, iout, ipunch
      common /ipo/ ipo(5)
      common /elp/ e(npoints)
      common /points/ p(3, npoints), maxpnts
c+++
      common /out/ ntitle(20, 3), chkfil, x(400), i6to5, rms, percent,
     &nlin, nend(3)
c+++
c      COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS,IAN(101),
c
c     $             ATMCHG(100),C(3,100)
c
c
c
      common /mol/ natoms, icharg, multip, nae, nbe, nel, nbasis, ian(
     &401), atmchg(400), c(3, 400)
      dimension a(400, 400), y(400), is(2, 400), iad1(400), iad2(400)
c
c
c        DEBYE = CONVERSION FROM DEBYES TO AU
c
c
c
      dimension d(400), which1(3)
c
c
c        SET UP MATRIX OF LINEAR COEFFICIENTS, A
c
c
c
c        BEGIN LOOP OVER ROWS
c
c
c
c
c
c        BEGIN LOOP OVER COLUMNS
c
c
c
      data one / 1.0 /
      data zero / 0.0 /
      data debye / 0.393427328 /
      data maxdim / 400 /
      data au2cal / 627.51 /
      data half / 0.5 /
      data hundred / 100.0 /
      data nconstr / 1 /
c# 397 "chelpg.for"
      do 100 k = 1, natoms
c
c
      do 200 mu = 1, natoms
c# 403 "chelpg.for"
      sum = zero
      do 400 i = 1, maxpnts
      rik = (((p(1,i) - c(1,k)) ** 2) + ((p(2,i) - c(2,k)) ** 2)) + ((p(
     &3,i) - c(3,k)) ** 2)
c# 406 "chelpg.for"
      rik = sqrt(rik)
      rimu = (((p(1,i) - c(1,mu)) ** 2) + ((p(2,i) - c(2,mu)) ** 2)) + (
     &(p(3,i) - c(3,mu)) ** 2)
      rimu = sqrt(rimu)
      sum = sum + (one / (rik * rimu))
c
c
c# 411 "chelpg.for"
  400 continue
c# 413 "chelpg.for"
      a(k,mu) = sum
c
c
c        FILL OUT COLUMNS CORRESPONDING TO LAGRANGE MULTIPLIERS
c
c
c
c# 414 "chelpg.for"
  200 continue
c
c
c
c
c# 418 "chelpg.for"
      a(k,natoms + 1) = half
c
c
c        FILL OUT THE ROWS CORRESPONDING TO CONSTRAINTS
c
c
c
c# 421 "chelpg.for"
  100 continue
c# 425 "chelpg.for"
      do 500 mu = 1, natoms
c
c
c# 426 "chelpg.for"
      a(natoms + 1,mu) = one
c
c
c        FILL OUT THE BLOCK WHICH CONNECTS LAGRANGE MULTIPLIERS TO
c
c        CONSTRAINTS
c
c
c
c# 428 "chelpg.for"
  500 continue
c# 433 "chelpg.for"
      do 600 k = natoms + 1, natoms + nconstr
      do 600 mu = natoms + 1, natoms + nconstr
      a(k,mu) = zero
c
c
c****DEBUG*****
c
c
c
c# 436 "chelpg.for"
  600 continue
c# 440 "chelpg.for"
      if (ipo(3) .eq. 1) then
      write(unit=iout, fmt=*) 'A MATRIX'
      do 699 k = 1, natoms + nconstr
      write(unit=iout, fmt=1699) (a(k,mu),mu = 1, natoms + nconstr)
 1699 format(1x,10f10.4)
  699 continue
c***************
c
c
c
c        CONSTRUCT COLUMN VECTOR, Y
c
c
c
c# 446 "chelpg.for"
      end if
c# 451 "chelpg.for"
      do 700 k = 1, natoms
      sum = zero
      do 710 i = 1, maxpnts
      rik = (((p(1,i) - c(1,k)) ** 2) + ((p(2,i) - c(2,k)) ** 2)) + ((p(
     &3,i) - c(3,k)) ** 2)
      rik = sqrt(rik)
      sum = sum + (e(i) / rik)
  710 continue
      y(k) = sum
      if (ipo(3) .eq. 1) write(unit=iout, fmt=*) k, y(k)
c
c
c        CONSTRUCT THE PORTION OF Y CORRESPONDING TO LAGRANGE MULTIPLIER
cS
c
c
c
c
c# 461 "chelpg.for"
  700 continue
c
c
c
c
c# 466 "chelpg.for"
      y(natoms + 1) = dble(icharg)
c
c
c        SOLVE MATRIX EQUATION AX = Y;
c
c        WHERE X = (Q1,Q2, ... QN,L1,L2, ... ,LN)
c
c
c
c        X = A(INV)Y
c
c
c
c        INVERT A
c
c
c
c# 469 "chelpg.for"
      if (ipo(3) .eq. 1) write(unit=iout, fmt=*) 'COL VECTR Y', (y(kk),
     &kk = 1, natoms + nconstr)
c
c
c****DEBUG*****
c
c
c
c# 479 "chelpg.for"
      call inv(a, natoms + nconstr, is, iad1, iad2, d, maxdim)
c# 483 "chelpg.for"
      if (ipo(3) .eq. 1) then
      write(unit=iout, fmt=*) 'A INVERSE'
      do 799 k = 1, natoms + nconstr
      write(unit=iout, fmt=1699) (a(k,mu),mu = 1, natoms + nconstr)
  799 continue
c**************
c
c
c
c        PERFORM MATRIX MULTIPLICATION A(INV)Y
c
c
c
c# 488 "chelpg.for"
      end if
c
c
c# 493 "chelpg.for"
      call multay(a, y, x, natoms + nconstr, maxdim)
c# 495 "chelpg.for"
      if (ipo(3) .eq. 1) then
      write(unit=iout, fmt=*) 'CHARGES:  '
      do 899 i = 1, natoms
      write(unit=iout, fmt=*) ian(i), x(i)
  899 continue
c
c
c        COMPUTE RMS DEVIATION AND MEAN ABSOLUTE % DEVIATION
c
c
c
c# 500 "chelpg.for"
      end if
c# 504 "chelpg.for"
      rms = zero
      percent = zero
      do 800 i = 1, maxpnts
      eq = zero
      do 810 j = 1, natoms
      dist = (((p(1,i) - c(1,j)) ** 2) + ((p(2,i) - c(2,j)) ** 2)) + ((p
     &(3,i) - c(3,j)) ** 2)
      dist = sqrt(dist)
      eq = eq + (x(j) / dist)
  810 continue
      rms = rms + ((e(i) - eq) ** 2)
      percent = percent + abs(((e(i) - eq) / e(i)) * hundred)
      if (ipo(3) .eq. 1) write(unit=iout, fmt=*) 'ACTUAL,CALC ', e(i),
     &eq
c# 517 "chelpg.for"
  800 continue
      if (ipo(3) .eq. 1) write(unit=iout, fmt=*) 'SUM OF SQUARES ', rms
      rms = (sqrt(rms) * au2cal) / maxpnts
      percent = percent / maxpnts
      if (ipo(3) .eq. 1) write(unit=iout, fmt=*) 'RMS, %', rms, percent
      return
      end
c
c
      subroutine fmgen(f, t, m)
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
c
c
      common /io/ in, iout, ipunch
      dimension f(m)
c
c
      dimension ga(35)
c
c
      equivalence (oldsum, approx)
c
c
 2001 format(42h1FAILURE IN FMGEN FOR SMALL T:  IX.GT.50, /6h IX = ,i3,
     &7h,  T = ,e20.14)
c
c
 2002 format(37h1FAILURE IN FMGEN FOR INTERMEDIATE T,/6h  T = ,e20.14)
      data zero / 0.0e0 /
      data half / 0.5e0 /
      data one / 1.0e0 /
      data two / 2.0e0 /
      data ten / 10.0e0 /
      data pi / 3.14159265358979e0 /
      data f42 / 42.0e0 /
      data f80 / 80.0e0 /
c# 542 "chelpg.for"
      texp = zero
      if (t - f80) 2, 3, 3
    2 texp = exp(- t)
    3 continue
c***********************************************************************
c
c        0 .LT. T .LT. 10
c
c***********************************************************************
c
c# 546 "chelpg.for"
      if (t - ten) 10, 70, 70
c# 550 "chelpg.for"
   10 term = (half * ga(m)) * texp
      tx = one
      ix = m + 1
      sum = tx / ga(ix)
      oldsum = sum
   20 ix = ix + 1
      tx = tx * t
      if (ix - 35) 40, 40, 30
   30 write(unit=iout, fmt=2001) ix, t
      stop 'FMGEN'
   40 sum = sum + (tx / ga(ix))
      if (tol - abs((oldsum / sum) - one)) 50, 60, 60
   50 oldsum = sum
      goto 20
   60 f(m) = sum * term
c
c
c# 565 "chelpg.for"
      goto 160
c***********************************************************************
c
c        10 .LE. T .LT. 42
c
c***********************************************************************
c
c# 567 "chelpg.for"
   70 if (t - f42) 80, 150, 150
c# 571 "chelpg.for"
   80 a = float(m - 1)
      b = a + half
      a = a - half
      tx = one / t
      mm1 = m - 1
      approx = (rpitwo * sqrt(tx)) * (tx ** mm1)
      if (mm1) 90, 110, 90
   90 do 100 ix = 1, mm1
      b = b - one
  100 approx = approx * b
  110 fimult = (half * texp) * tx
      sum = zero
      if (fimult) 120, 140, 120
  120 fiprop = fimult / approx
      term = one
      sum = one
      notrms = int(t) + mm1
      do 130 ix = 2, notrms
      term = (term * a) * tx
      sum = sum + term
      if (abs((term * fiprop) / sum) - tol) 140, 140, 130
  130 a = a - one
      write(unit=iout, fmt=2002) t
      stop 'FMGEN'
  140 f(m) = approx - (fimult * sum)
c***********************************************************************
c
c        T .GE. 42
c
c***********************************************************************
c
c# 596 "chelpg.for"
      goto 160
c# 600 "chelpg.for"
  150 tx = float(m) - half
c***********************************************************************
c
c        RECUR DOWNWARDS TO F(1)
c
c***********************************************************************
c
c# 601 "chelpg.for"
      f(m) = (half * ga(m)) / (t ** tx)
c# 605 "chelpg.for"
  160 tx = t + t
      sum = float((m + m) - 3)
      mm1 = m - 1
      if (mm1) 170, 190, 170
  170 do 180 ix = 1, mm1
      f(m - ix) = ((tx * f((m - ix) + 1)) + texp) / sum
  180 sum = sum - two
c
c
c# 612 "chelpg.for"
  190 return
c
c
c# 614 "chelpg.for"
      entry fmset()
c# 616 "chelpg.for"
      ga(1) = sqrt(pi)
      tol = half
      do 200 i = 2, 35
      ga(i) = ga(i - 1) * tol
  200 tol = tol + one
      tol = 5.0e-09
      rpitwo = half * ga(1)
      return
      end
c
c
c     ROUTINE TO CALCULATE THE ELECTRON-CHARGE MATRIX ELEMENTS FOR THE
c
c     POLARIZATION POTENTIAL. CODE REVISED FROM THE ONE ELECTRON PACKAGE
c
c     AS IT EXISTED AUGUST, 1983.
c
c
c
c
c
c        REVISED BY M.M. FRANCL JANUARY 1984 FOR PRINCETON CHEMISTRY
c
c        DEPARTMENT VAX 11/780
c
c
c
c        REVISED TO BE COMPATIBLE WITH COMMON /B/ FROM GAUSSIAN 82
c
c       MAY 1984 M.M. FRANCL
c
c
c
c        REVISED TO USE ** BASIS SETS AND THOSE HAVING P ONLY SHELLS
c
c        JANUARY 1986 M.M. FRANCL
c
c
c
c	 REVISED FOR GAUSSIAN 86 CHECKPOINT FILES FOR YALE UNIVERSITY
c	FEBRUARY 1988 CURT BRENEMAN
c
c
      subroutine intgrl(h, x1, x2, x3, icharg, i6to5)
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
c
c
c+++
      integer*4 shella, shelln, shellt, shellc, aos, aon, shladf
c
c===  Gaussian88 Modification.  New Common /b/ size.
      common /mol/ natoms, jcharg, multip, nae, nbe, nel, nbasis, ian(
     &401), atmchg(400), c(3, 400)
c
c===	Old G86 common /b/
c      COMMON/B/EXX(1200),C1(1200),C2(1200),C3(400),CF(400),SHLADF(800),
c     $         X(400),Y(400),Z(400),JAN(400),SHELLA(400),SHELLN(400),
c     $         SHELLT(400),SHELLC(400),AOS(400),AON(400),NSHELL,MAXTYP
c
c+++
c      COMMON /B/ EXX(240),C1(240),C2(240),C3(80),CF(80),SHLADF(160),
c
c     $           X(80),Y(80),Z(80),
c
c     $           JAN(80),SHELLA(80),SHELLN(80),SHELLT(80),SHELLC(80)
c
c     $          ,AOS(80),AON(80),NSHELL,MAXTYP
c
c      COMMON /MOL/ NATOMS,JCHARG,MULTIP,NAE,NBE,NEL,NBASIS,IAN(101),
c
c     $             ATMCHG(100),C(3,100)
c
c g88/90
c      common /b/ exx(6000), c1(6000), c2(6000), c3(2000), cf(2000),
c     &shladf(4000), x(2000), y(2000), z(2000), jan(2000), shella(2000),
c     &shelln(2000), shellt(2000), shellc(2000), aos(2000), aon(2000),
c     &nshell, maxtyp
c g88 cray
      common /b/ exx(6000), c1(6000), c2(6000), c3(2000), cf(2000),
     &shladf(2000), x(2000), y(2000), z(2000), jan(2000), shella(2000),
     &shelln(2000), shellt(2000), shellc(2000), aos(2000), aon(2000),
     &nshell, maxtyp
      common /ipo/ ipo(5)
c
c
      common /io/ in, iout, ipunch
      dimension h(1)
      dimension renorm(10)
      dimension of(9), ox(9), tx(13), abx(5), aby(5), abz(5), absq(5), a
     &(5), b(5), f(5), apb(5), cpx(5), cpy(5), cpz(5), fm(5)
c
c
      dimension epn(100)
c
c
      common /h100/ ep00, ep10, ep20, ep30, ep40, ep50, ep60, ep70, ep80
     &, ep90, ep01, ep11, ep21, ep31, ep41, ep51, ep61, ep71, ep81, ep91
     &, ep02, ep12, ep22, ep32, ep42, ep52, ep62, ep72, ep82, ep92, ep03
     &, ep13, ep23, ep33, ep43, ep53, ep63, ep73, ep83, ep93, ep04, ep14
     &, ep24, ep34, ep44, ep54, ep64, ep74, ep84, ep94, ep05, ep15, ep25
     &, ep35, ep45, ep55, ep65, ep75, ep85, ep95, ep06, ep16, ep26, ep36
     &, ep46, ep56, ep66, ep76, ep86, ep96, ep07, ep17, ep27, ep37, ep47
     &, ep57, ep67, ep77, ep87, ep97, ep08, ep18, ep28, ep38, ep48, ep58
     &, ep68, ep78, ep88, ep98, ep09, ep19, ep29, ep39, ep49, ep59, ep69
     &, ep79, ep89, ep99
      dimension eep(100)
c
c SPDSTV
c     LOCAL VARIABLES.
c SPDSTV
c
c
      dimension max(6)
      dimension ag(6), csa(6), cpa(6), cda(6), bg(6), csb(6), cpb(6),
     &cdb(6), dpp(9)
      equivalence (of(1), of0), (of(2), of1), (of(3), of2), (of(4), of3)
     &, (of(5), of4), (of(6), of5), (of(7), of6), (of(8), of7), (of(9),
     &of8)
      equivalence (ox(1), ox0), (ox(2), ox1), (ox(3), ox2), (ox(4), ox3)
     &, (ox(5), ox4), (ox(6), ox5), (ox(7), ox6), (ox(8), ox7), (ox(9),
     &ox8)
      equivalence (a(2), a1), (a(3), a2), (a(4), a3), (a(5), a4)
      equivalence (b(2), b1), (b(3), b2), (b(4), b3), (b(5), b4)
      equivalence (t0, t01), (t1, t02), (t2, t03), (t3, t04), (t4, t05)
     &, (t5, t06), (t6, t07), (t7, t08), (t8, t09)
      equivalence (tx(10), t10), (tx(11), t11), (tx(12), t12), (tx(13),
     &t13)
      equivalence (tx(1), t0), (tx(2), t1), (tx(3), t2), (tx(4), t3), (
     &tx(5), t4), (tx(6), t5), (tx(7), t6), (tx(8), t7), (tx(9), t8)
      equivalence (t01, c001), (t02, c050), (t09, c054), (t13, c067), (
     &t08, c068), (t03, c074)
      equivalence (abx(2), abx1), (abx(3), abx2), (abx(4), abx3), (abx(5
     &), abx4)
      equivalence (abx1, ab004), (abx2, ab006), (abx3, ab023), (abx4,
     &ab029)
      equivalence (aby(2), aby1), (aby(3), aby2), (aby(4), aby3), (aby(5
     &), aby4)
      equivalence (aby1, ab007), (aby2, ab010), (aby3, ab032), (aby4,
     &ab035)
      equivalence (abz(2), abz1), (abz(3), abz2), (abz(4), abz3), (abz(5
     &), abz4)
      equivalence (abz1, ab002), (abz2, ab003), (abz3, ab011), (abz4,
     &ab017)
      equivalence (absq(2), absq1), (absq(3), absq2), (absq(4), absq3),
     &(absq(5), absq4)
      equivalence (apb(2), apb1), (apb(3), apb2), (apb(4), apb3), (apb(5
     &), apb4)
      equivalence (cpx(2), cpx1), (cpx(3), cpx2), (cpx(4), cpx3), (cpx(5
     &), cpx4)
      equivalence (cpy(2), cpy1), (cpy(3), cpy2), (cpy(4), cpy3), (cpy(5
     &), cpy4)
      equivalence (cpz(2), cpz1), (cpz(3), cpz2), (cpz(4), cpz3), (cpz(5
     &), cpz4)
      equivalence (f(2), f1), (f(3), f2), (f(4), f3), (f(5), f4)
      equivalence (fm(1), fm0), (fm(2), fm1), (fm(3), fm2), (fm(4), fm3)
     &, (fm(5), fm4)
      equivalence (fm0, d001)
c
c
      equivalence (eep(1), ep00)
c
c
c
c
c        CALL ROUTINE TO MODIFY COMMON /B/ IF P ONLY SHELLS ARE PRESENT
c
c
c
 2010 format(/1x,31hELECTRON-CHARGE MATRIX ELEMENTS/)
c
c
c***********************************************************************
c SPDSTV
c        INITIALIZE THIS SEGMENT.
c SPDSTV
c***********************************************************************
c SPDSTV
c
c SPDSTV
c     ******************************************************************
c SPDSTV
c     COMPUTE SIZE OF S T AND V ARRAYS
c
c     ******************************************************************
c SPDSTV
      data max / 1, 4, 9, 1, 4, 10 /
      data tx / 1.0e0, 0.5e0, 0.25e0, 0.125e0, 0.375e0, 0.625e-01,
     &0.1875e0, 0.75e0, 1.5e0, 2.25e0, 1.125e0, 0.0e0, 3.0e0 /
      data zero / 0.0 /
      data half / 0.5 /
      data one / 1.0 /
      data onept5 / 1.5 /
      data two / 2.0 /
      data three / 3.0 /
      data root3 / 1.732050808 /
      data pi / 3.14159265358979 /
      data antoau / 1.889726878d0 /
c# 753 "chelpg.for"
      call star(nbasis, shellt, shellc, aos, nshell, nostar)
c# 762 "chelpg.for"
      ntt = (nbasis * (nbasis + 1)) / 2
cC    IF(IGO(4) .NE. 0) I5OR6 = 0
c
c     ******************************************************************
c SPDSTV
c     INITIALIZE RENORM  USED TO NORMALIZE D FUNCTIONS
c
c     ******************************************************************
c SPDSTV
c# 763 "chelpg.for"
      i5or6 = 3
c# 768 "chelpg.for"
      do 100 i = 1, 10
  100 renorm(i) = one
      renorm(5) = root3
      renorm(8) = root3
c     ******************************************************************
c SPDSTV
c     CLEAR H ARRAY
c
c     ******************************************************************
c SPDSTV
c# 772 "chelpg.for"
      renorm(9) = root3
c# 776 "chelpg.for"
      do 50 i = 1, ntt
c     ******************************************************************
c SPDSTV
c     *  INITIALIZE THE VARIABLES USED BY ROUTINE FMGEN.               *
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 777 "chelpg.for"
   50 h(i) = zero
c# 781 "chelpg.for"
      call fmset
      do 95 i = 1, 5
   95 fm(i) = zero
      abx(1) = one
      aby(1) = one
      abz(1) = one
      a(1) = one
      b(1) = one
      f(1) = one
      cpx(1) = one
      cpy(1) = one
      cpz(1) = one
      apb(1) = one
c***********************************************************************
c SPDSTV
c        LOOP OVER SHELLS ISHELL AND JSHELL.
c SPDSTV
c***********************************************************************
c SPDSTV
c# 794 "chelpg.for"
      absq(1) = one
c# 798 "chelpg.for"
      do 1000 ishell = 1, nshell
      do 1000 jshell = 1, ishell
c     ******************************************************************
c SPDSTV
c     ZERO LOCATIONS
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 800 "chelpg.for"
      symfac = one
c# 804 "chelpg.for"
   80 continue
      do 9447 ji = 1, 100
      epn(ji) = zero
 9447 continue
      if (shellt(ishell) - shellt(jshell)) 120, 120, 110
  110 inew = jshell
      jnew = ishell
      la = shellt(jshell)
      lb = shellt(ishell)
      goto 200
  120 inew = ishell
      jnew = jshell
      la = shellt(ishell)
      lb = shellt(jshell)
  200 continue
      lap1 = la + 1
      lbp1 = lb + 1
      lamax = max(lap1 + i5or6)
      lbmax = max(lbp1 + i5or6)
      itype = (3 * lb) + la
      m = (la + lb) + 1
      nga = shelln(inew)
      ngb = shelln(jnew)
      ax = x(inew)
      bx = x(jnew)
      ay = y(inew)
      by = y(jnew)
      az = z(inew)
      bz = z(jnew)
      isha = shella(inew)
      ishb = shella(jnew)
      ishad = shladf(inew)
      ishbd = shladf(jnew)
      iaos = aos(inew)
c     ******************************************************************
c SPDSTV
c     OBTAIN INFORMATION ABOUT SHELLS INEW AND JNEW
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 838 "chelpg.for"
      jaos = aos(jnew)
c# 842 "chelpg.for"
      do 101 i = 1, nga
      n = (isha + i) - 1
      nd = (ishad + i) - 1
      if (maxtyp .le. 1) nd = 1
      ag(i) = exx(n)
      csa(i) = c1(n)
      cpa(i) = c2(n)
  101 cda(i) = c3(nd)
c# 851 "chelpg.for"
      do 102 i = 1, ngb
      n = (ishb + i) - 1
      nd = (ishbd + i) - 1
      bg(i) = exx(n)
      csb(i) = c1(n)
      cpb(i) = c2(n)
  102 cdb(i) = c3(nd)
c# 859 "chelpg.for"
      abx(2) = bx - ax
      aby(2) = by - ay
      abz(2) = bz - az
      rabsq = ((abx(2) * abx(2)) + (aby(2) * aby(2))) + (abz(2) * abz(2)
     &)
c# 863 "chelpg.for"
      absq(2) = rabsq
      do 103 i = 3, 5
      abx(i) = abx(i - 1) * abx(2)
      aby(i) = aby(i - 1) * aby(2)
      abz(i) = abz(i - 1) * abz(2)
  103 absq(i) = absq(i - 1) * absq(2)
      ab001 = one
      ab005 = abx1 * abz1
      ab008 = aby1 * abz1
      ab009 = abx1 * aby1
      ab012 = abx1 * abz2
      ab013 = abx2 * abz1
      ab014 = aby1 * abz2
      ab015 = (abx1 * aby1) * abz1
      ab016 = aby2 * abz1
      ab018 = abx1 * abz3
      ab019 = abx2 * abz2
      ab020 = aby1 * abz3
      ab021 = (abx1 * aby1) * abz2
      ab022 = aby2 * abz2
      ab024 = abx2 * aby1
      ab025 = abx1 * aby2
      ab026 = abx3 * abz1
      ab027 = (abx2 * aby1) * abz1
      ab028 = (abx1 * aby2) * abz1
      ab030 = abx3 * aby1
      ab031 = abx2 * aby2
      ab033 = aby3 * abz1
c***********************************************************************
c SPDSTV
c        LOOP OVER GAUSSIANS  (CONTRACTION LOOP).
c SPDSTV
c***********************************************************************
c SPDSTV
c# 891 "chelpg.for"
      ab034 = abx1 * aby3
c# 895 "chelpg.for"
      do 105 igauss = 1, nga
      aa = ag(igauss)
      do 105 jgauss = 1, ngb
      bb = bg(jgauss)
      aapbb = aa + bb
      apbb = one / aapbb
      f2 = ((two * aa) * bb) * apbb
      px = ((aa * ax) + (bb * bx)) * apbb
      py = ((aa * ay) + (bb * by)) * apbb
      pz = ((aa * az) + (bb * bz)) * apbb
      a(2) = one / aa
      b(2) = one / bb
      f(2) = f2
      apb(2) = apbb
      yx = pi * apbb
      exx1 = (half * f2) * rabsq
      if (exx1 - 80.0e0) 4172, 4173, 4173
 4173 exx1 = zero
      goto 4714
 4172 exx1 = exp(- exx1)
 4714 continue
      ov = (yx ** onept5) * exx1
      ovek = ((three * aa) * bb) * apbb
      ek = (((f2 * aa) * bb) * apbb) * ov
      ep = (two * yx) * exx1
      do 119 i = 3, 5
      a(i) = a(i - 1) * a(2)
      b(i) = b(i - 1) * b(2)
      apb(i) = apb(i - 1) * apb(2)
  119 f(i) = f(i - 1) * f(2)
      dpp(1) = csa(igauss) * csb(jgauss)
      dpp(2) = cpa(igauss) * csb(jgauss)
      dpp(3) = cda(igauss) * csb(jgauss)
      dpp(4) = csa(igauss) * cpb(jgauss)
      dpp(5) = cpa(igauss) * cpb(jgauss)
      dpp(6) = cda(igauss) * cpb(jgauss)
      dpp(7) = csa(igauss) * cdb(jgauss)
      dpp(8) = cpa(igauss) * cdb(jgauss)
      dpp(9) = cda(igauss) * cdb(jgauss)
      do 2132 i = 1, 9
      of(i) = dpp(i) * ov
 2132 ox(i) = dpp(i) * ek
      do 2139 i = 1, 100
 2139 eep(i) = zero
      c002 = (t02 * a1) * f1
      c006 = (t02 * b1) * f1
      c007 = ((t03 * a1) * b1) * f2
      c008 = ((t03 * a1) * b1) * f1
      c027 = t01 * a1
      c031 = ((t01 * a1) * b1) * f1
      c032 = (t02 * a1) * b1
      c051 = ((t02 * a1) * b1) * f2
      c012 = t02 * b1
      c013 = (t03 * b2) * f2
      c014 = (t03 * b2) * f1
      c036 = (t01 * b2) * f1
      c037 = t02 * b2
      c056 = (t01 * b1) * f1
      c030 = t01 * b1
      c018 = ((t04 * a1) * b2) * f2
      if (itype - 7) 3060, 3040, 3041
 3041 continue
      c003 = t02 * a1
      c004 = (t03 * a2) * f2
      c005 = (t03 * a2) * f1
      c009 = ((t04 * a2) * b1) * f3
      c010 = ((t05 * a2) * b1) * f2
      c011 = ((t04 * a2) * b1) * f2
      c017 = (t03 * a1) * b1
      c019 = ((t04 * a1) * b2) * f1
      c020 = ((t04 * a2) * b1) * f1
      c021 = ((t06 * a2) * b2) * f4
      c022 = ((t05 * a2) * b2) * f3
      c023 = ((t07 * a2) * b2) * f2
      c024 = ((t07 * a2) * b2) * f3
      c025 = ((t06 * a2) * b2) * f3
      c026 = ((t06 * a2) * b2) * f2
      c028 = (t01 * a2) * f1
      c029 = t02 * a2
      c033 = ((t08 * a2) * b1) * f2
      c034 = ((t09 * a2) * b1) * f1
      c035 = ((t02 * a2) * b1) * f1
      c040 = ((t02 * a1) * b2) * f1
      c041 = (t03 * a1) * b2
      c042 = (t03 * a2) * b1
      c043 = ((t02 * a2) * b2) * f3
      c044 = ((t10 * a2) * b2) * f2
      c045 = ((t08 * a2) * b2) * f1
      c046 = ((t11 * a2) * b2) * f2
      c047 = ((t05 * a2) * b2) * f2
      c048 = ((t03 * a2) * b2) * f1
      c049 = (t01 * a1) * f1
      c057 = ((t12 * a1) * b1) * f1
      c058 = t03 * a1
      c059 = t03 * b1
      c060 = ((t03 * a1) * b2) * f3
      c061 = (t04 * b2) * f2
      c062 = (t04 * b2) * f1
      c063 = ((t03 * a2) * b1) * f3
      c064 = ((t01 * a1) * b1) * f2
      c065 = (t09 * b1) * f1
      c066 = (t09 * a1) * f1
      c069 = (t04 * a2) * f2
      c070 = (t04 * a2) * f1
      c071 = ((t03 * a2) * b1) * f2
      c072 = (t08 * a1) * f1
      c073 = ((t03 * a1) * b2) * f2
      c075 = (t08 * b1) * f1
      c076 = ((t04 * a1) * b1) * f2
 3040 continue
      c015 = ((t04 * a1) * b2) * f3
      c016 = ((t05 * a1) * b2) * f2
      c038 = ((t08 * a1) * b2) * f2
      c039 = ((t09 * a1) * b2) * f1
      c040 = ((t02 * a1) * b2) * f1
      c052 = ((t02 * a1) * b1) * f1
      c053 = (t03 * b1) * f1
      c055 = (t03 * a1) * f1
 3060 continue
      cx = x1
      cy = x2
      cz = x3
      cpx(2) = px - cx
      cpy(2) = py - cy
      cpz(2) = pz - cz
      cp2 = ((cpx(2) * cpx(2)) + (cpy(2) * cpy(2))) + (cpz(2) * cpz(2))
      call fmgen(fm, aapbb * cp2, m)
      do 108 i = 3, 5
      cpx(i) = cpx(i - 1) * cpx(2)
      cpy(i) = cpy(i - 1) * cpy(2)
  108 cpz(i) = cpz(i - 1) * cpz(2)
      epan = ep * float(- icharg)
      do 2136 i = 1, 9
 2136 of(i) = dpp(i) * epan
      d002 = cpz1 * fm1
      d003 = cpz2 * fm2
      d004 = apb1 * fm1
      d005 = cpx1 * fm1
      d006 = (cpx1 * cpz1) * fm2
      d007 = cpx2 * fm2
      d008 = cpy1 * fm1
      d009 = (cpy1 * cpz1) * fm2
      d010 = (cpx1 * cpy1) * fm2
      d011 = cpy2 * fm2
      d012 = cpz3 * fm3
      d013 = (apb1 * cpz1) * fm2
      d014 = (cpx1 * cpz2) * fm3
      d015 = (apb1 * cpx1) * fm2
      d016 = (cpx2 * cpz1) * fm3
      d017 = (cpy1 * cpz2) * fm3
      d018 = (apb1 * cpy1) * fm2
      d019 = ((cpx1 * cpy1) * cpz1) * fm3
      d020 = (cpy2 * cpz1) * fm3
      d034 = cpx3 * fm3
      d035 = (cpx2 * cpy1) * fm3
      d036 = (cpx1 * cpy2) * fm3
c     ******************************************************************
c SPDSTV
c     *                                SS                              *
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 1051 "chelpg.for"
      d043 = cpy3 * fm3
c# 1055 "chelpg.for"
      ep00 = of0 * ((c001 * ab001) * d001)
c     ******************************************************************
c SPDSTV
c     *                                SP                              *
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 1056 "chelpg.for"
      if (itype) 3230, 3262, 3230
c# 1060 "chelpg.for"
 3230 continue
      ep01 = of3 * ((- ((c006 * ab002) * d001)) - ((c001 * ab001) * d002
     &))
c# 1062 "chelpg.for"
      ep03 = of3 * ((- ((c006 * ab004) * d001)) - ((c001 * ab001) * d005
     &))
c# 1063 "chelpg.for"
      ep06 = of3 * ((- ((c006 * ab007) * d001)) - ((c001 * ab001) * d008
     &))
c# 1064 "chelpg.for"
      if (itype - 7) 3240, 3242, 3241
c     ******************************************************************
c SPDSTV
c     *                                DD                              *
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 1065 "chelpg.for"
 3240 if (itype - 4) 3262, 3261, 3260
c# 1069 "chelpg.for"
 3241 continue
      d021 = cpz4 * fm4
      d022 = (apb1 * cpz2) * fm3
      d023 = apb2 * fm2
      d024 = (cpx1 * cpz3) * fm4
      d025 = ((apb1 * cpx1) * cpz1) * fm3
      d026 = (cpx2 * cpz2) * fm4
      d027 = (apb1 * cpx2) * fm3
      d028 = (cpy1 * cpz3) * fm4
      d029 = ((apb1 * cpy1) * cpz1) * fm3
      d030 = ((cpx1 * cpy1) * cpz2) * fm4
      d031 = ((apb1 * cpx1) * cpy1) * fm3
      d032 = (cpy2 * cpz2) * fm4
      d033 = (apb1 * cpy2) * fm3
      d037 = (cpx3 * cpz1) * fm4
      d038 = ((cpx2 * cpy1) * cpz1) * fm4
      d039 = ((cpx1 * cpy2) * cpz1) * fm4
      d040 = cpx4 * fm4
      d041 = (cpx3 * cpy1) * fm4
      d042 = (cpx2 * cpy2) * fm4
      d044 = (cpy3 * cpz1) * fm4
      d045 = (cpx1 * cpy3) * fm4
      d046 = cpy4 * fm4
      ep20 = of2 * (((((((c003 * ab001) * d001) + ((c004 * ab003) * d001
     &)) - ((c005 * ab001) * d001)) - ((c049 * ab002) * d002)) + ((c001
     & * ab001) * d003)) - ((c050 * ab001) * d004))
c# 1094 "chelpg.for"
      ep40 = of2 * (((((c004 * ab005) * d001) - ((c002 * ab004) * d002))
     & - ((c002 * ab002) * d005)) + ((c001 * ab001) * d006))
      ep50 = of2 * (((((((c003 * ab001) * d001) + ((c004 * ab006) * d001
     &)) - ((c005 * ab001) * d001)) - ((c049 * ab004) * d005)) + ((c001
     & * ab001) * d007)) - ((c050 * ab001) * d004))
c# 1098 "chelpg.for"
      ep70 = of2 * (((((c004 * ab008) * d001) - ((c002 * ab007) * d002))
     & - ((c002 * ab002) * d008)) + ((c001 * ab001) * d009))
      ep80 = of2 * (((((c004 * ab009) * d001) - ((c002 * ab004) * d008))
     & - ((c002 * ab007) * d005)) + ((c001 * ab001) * d010))
      ep90 = of2 * (((((((c003 * ab001) * d001) + ((c004 * ab010) * d001
     &)) - ((c005 * ab001) * d001)) - ((c049 * ab007) * d008)) + ((c001
     & * ab001) * d011)) - ((c050 * ab001) * d004))
c# 1104 "chelpg.for"
      ep21 = of5 * ((((((((((((((- ((c008 * ab002) * d001)) - ((c003 *
     &ab001) * d002)) - ((c009 * ab011) * d001)) + ((c010 * ab002) *
     &d001)) + ((c051 * ab003) * d002)) - ((c052 * ab001) * d002)) - ((
     &c006 * ab002) * d003)) + ((c053 * ab002) * d004)) - ((c004 * ab003
     &) * d002)) + ((c005 * ab001) * d002)) + ((c049 * ab002) * d003))
     & - ((c002 * ab002) * d004)) - ((c001 * ab001) * d012)) + ((c054 *
     &ab001) * d013))
c# 1108 "chelpg.for"
      ep41 = of5 * ((((((((((((- ((c009 * ab012) * d001)) + ((c011 *
     &ab004) * d001)) + ((c007 * ab005) * d002)) + ((c007 * ab003) *
     &d005)) - ((c008 * ab001) * d005)) - ((c006 * ab002) * d006)) - ((
     &c004 * ab005) * d002)) + ((c002 * ab004) * d003)) - ((c055 * ab004
     &) * d004)) + ((c002 * ab002) * d006)) - ((c001 * ab001) * d014))
     & + ((c050 * ab001) * d015))
c# 1112 "chelpg.for"
      ep51 = of5 * ((((((((((((- ((c008 * ab002) * d001)) - ((c003 *
     &ab001) * d002)) - ((c009 * ab013) * d001)) + ((c011 * ab002) *
     &d001)) + ((c051 * ab005) * d005)) - ((c006 * ab002) * d007)) + ((
     &c053 * ab002) * d004)) - ((c004 * ab006) * d002)) + ((c005 * ab001
     &) * d002)) + ((c049 * ab004) * d006)) - ((c001 * ab001) * d016))
     & + ((c050 * ab001) * d013))
c# 1116 "chelpg.for"
      ep71 = of5 * ((((((((((((- ((c009 * ab014) * d001)) + ((c011 *
     &ab007) * d001)) + ((c007 * ab008) * d002)) + ((c007 * ab003) *
     &d008)) - ((c008 * ab001) * d008)) - ((c006 * ab002) * d009)) - ((
     &c004 * ab008) * d002)) + ((c002 * ab007) * d003)) - ((c055 * ab007
     &) * d004)) + ((c002 * ab002) * d009)) - ((c001 * ab001) * d017))
     & + ((c050 * ab001) * d018))
c# 1120 "chelpg.for"
      ep81 = of5 * ((((((((- ((c009 * ab015) * d001)) + ((c007 * ab005)
     & * d008)) + ((c007 * ab008) * d005)) - ((c006 * ab002) * d010)) -
     &((c004 * ab009) * d002)) + ((c002 * ab004) * d009)) + ((c002 *
     &ab007) * d006)) - ((c001 * ab001) * d019))
c# 1123 "chelpg.for"
      ep91 = of5 * ((((((((((((- ((c008 * ab002) * d001)) - ((c003 *
     &ab001) * d002)) - ((c009 * ab016) * d001)) + ((c011 * ab002) *
     &d001)) + ((c051 * ab008) * d008)) - ((c006 * ab002) * d011)) + ((
     &c053 * ab002) * d004)) - ((c004 * ab010) * d002)) + ((c005 * ab001
     &) * d002)) + ((c049 * ab007) * d009)) - ((c001 * ab001) * d020))
     & + ((c050 * ab001) * d013))
c# 1127 "chelpg.for"
      ep22 = of8 * (((((((((((((((((((((((((((((((((((((c017 * ab001) *
     &d001) + ((c018 * ab003) * d001)) - ((c019 * ab001) * d001)) + ((
     &c057 * ab002) * d002)) + ((c003 * ab001) * d003)) - ((c058 * ab001
     &) * d004)) + ((c011 * ab003) * d001)) - ((c020 * ab001) * d001))
     & + ((c012 * ab001) * d003)) - ((c059 * ab001) * d004)) + ((c021 *
     &ab017) * d001)) - ((c022 * ab003) * d001)) - ((c060 * ab011) *
     &d002)) + ((c023 * ab001) * d001)) + ((c038 * ab002) * d002)) + ((
     &c013 * ab003) * d003)) - ((c061 * ab003) * d004)) - ((c014 * ab001
     &) * d003)) + ((c062 * ab001) * d004)) + ((c063 * ab011) * d002))
     & - ((c033 * ab002) * d002)) - ((c064 * ab003) * d003)) + ((c051 *
     &ab003) * d004)) + ((c031 * ab001) * d003)) - ((c052 * ab001) *
     &d004)) + ((c056 * ab002) * d012)) - ((c065 * ab002) * d013)) + ((
     &c004 * ab003) * d003)) - ((c005 * ab001) * d003)) - ((c049 * ab002
     &) * d012)) + ((c066 * ab002) * d013)) + ((c001 * ab001) * d021))
     & - ((c067 * ab001) * d022)) + ((c068 * ab001) * d023)) - ((c069 *
     &ab003) * d004)) + ((c070 * ab001) * d004))
c# 1136 "chelpg.for"
      ep42 = of8 * (((((((((((((((((((((((((((((c011 * ab005) * d001) -
     &((c008 * ab004) * d002)) - ((c008 * ab002) * d005)) + ((c012 *
     &ab001) * d006)) + ((c021 * ab018) * d001)) - ((c024 * ab005) *
     &d001)) - ((c015 * ab012) * d002)) - ((c015 * ab011) * d005)) + ((
     &c016 * ab002) * d005)) + ((c013 * ab003) * d006)) + ((c018 * ab004
     &) * d002)) - ((c014 * ab001) * d006)) + ((c063 * ab012) * d002))
     & - ((c071 * ab004) * d002)) - ((c051 * ab005) * d003)) + ((c007 *
     &ab005) * d004)) - ((c051 * ab003) * d006)) + ((c052 * ab001) *
     &d006)) + ((c056 * ab002) * d014)) - ((c006 * ab002) * d015)) + ((
     &c004 * ab005) * d003)) - ((c002 * ab004) * d012)) + ((c072 * ab004
     &) * d013)) - ((c002 * ab002) * d014)) + ((c001 * ab001) * d024))
     & - ((c054 * ab001) * d025)) - ((c069 * ab005) * d004)) + ((c055 *
     &ab002) * d015))
c# 1143 "chelpg.for"
      ep52 = of8 * (((((((((((((((((((((((((((((((((((((c017 * ab001) *
     &d001) + ((c018 * ab003) * d001)) - ((c019 * ab001) * d001)) + ((
     &c052 * ab002) * d002)) + ((c003 * ab001) * d003)) - ((c058 * ab001
     &) * d004)) + ((c011 * ab006) * d001)) - ((c020 * ab001) * d001))
     & - ((c052 * ab004) * d005)) + ((c012 * ab001) * d007)) - ((c059 *
     &ab001) * d004)) + ((c021 * ab019) * d001)) - ((c025 * ab003) *
     &d001)) - ((c060 * ab012) * d005)) + ((c013 * ab003) * d007)) - ((
     &c061 * ab003) * d004)) - ((c025 * ab006) * d001)) + ((c026 * ab001
     &) * d001)) + ((c073 * ab004) * d005)) - ((c014 * ab001) * d007))
     & + ((c062 * ab001) * d004)) + ((c063 * ab013) * d002)) - ((c071 *
     &ab002) * d002)) - ((c064 * ab005) * d006)) + ((c056 * ab002) *
     &d016)) - ((c006 * ab002) * d013)) + ((c004 * ab006) * d003)) - ((
     &c005 * ab001) * d003)) - ((c049 * ab004) * d014)) + ((c001 * ab001
     &) * d026)) - ((c050 * ab001) * d022)) - ((c069 * ab006) * d004))
     & + ((c070 * ab001) * d004)) + ((c002 * ab004) * d015)) - ((c050 *
     &ab001) * d027)) + ((c074 * ab001) * d023))
c# 1152 "chelpg.for"
      ep72 = of8 * (((((((((((((((((((((((((((((c011 * ab008) * d001) -
     &((c008 * ab007) * d002)) - ((c008 * ab002) * d008)) + ((c012 *
     &ab001) * d009)) + ((c021 * ab020) * d001)) - ((c024 * ab008) *
     &d001)) - ((c015 * ab014) * d002)) - ((c015 * ab011) * d008)) + ((
     &c016 * ab002) * d008)) + ((c013 * ab003) * d009)) + ((c018 * ab007
     &) * d002)) - ((c014 * ab001) * d009)) + ((c063 * ab014) * d002))
     & - ((c071 * ab007) * d002)) - ((c051 * ab008) * d003)) + ((c007 *
     &ab008) * d004)) - ((c051 * ab003) * d009)) + ((c052 * ab001) *
     &d009)) + ((c056 * ab002) * d017)) - ((c006 * ab002) * d018)) + ((
     &c004 * ab008) * d003)) - ((c002 * ab007) * d012)) + ((c072 * ab007
     &) * d013)) - ((c002 * ab002) * d017)) + ((c001 * ab001) * d028))
     & - ((c054 * ab001) * d029)) - ((c069 * ab008) * d004)) + ((c055 *
     &ab002) * d018))
c# 1159 "chelpg.for"
      ep82 = of8 * (((((((((((((((((((((((((c011 * ab009) * d001) - ((
     &c008 * ab004) * d008)) - ((c008 * ab007) * d005)) + ((c012 * ab001
     &) * d010)) + ((c021 * ab021) * d001)) - ((c015 * ab012) * d008))
     & - ((c015 * ab014) * d005)) + ((c013 * ab003) * d010)) - ((c025 *
     &ab009) * d001)) + ((c018 * ab004) * d008)) + ((c018 * ab007) *
     &d005)) - ((c014 * ab001) * d010)) + ((c063 * ab015) * d002)) - ((
     &c051 * ab005) * d009)) - ((c051 * ab008) * d006)) + ((c056 * ab002
     &) * d019)) + ((c004 * ab009) * d003)) - ((c002 * ab004) * d017))
     & - ((c002 * ab007) * d014)) + ((c001 * ab001) * d030)) - ((c069 *
     &ab009) * d004)) + ((c055 * ab004) * d018)) + ((c055 * ab007) *
     &d015)) - ((c050 * ab001) * d031))
c# 1165 "chelpg.for"
      ep92 = of8 * (((((((((((((((((((((((((((((((((((((c017 * ab001) *
     &d001) + ((c018 * ab003) * d001)) - ((c019 * ab001) * d001)) + ((
     &c052 * ab002) * d002)) + ((c003 * ab001) * d003)) - ((c058 * ab001
     &) * d004)) + ((c011 * ab010) * d001)) - ((c020 * ab001) * d001))
     & - ((c052 * ab007) * d008)) + ((c012 * ab001) * d011)) - ((c059 *
     &ab001) * d004)) + ((c021 * ab022) * d001)) - ((c025 * ab003) *
     &d001)) - ((c060 * ab014) * d008)) + ((c013 * ab003) * d011)) - ((
     &c061 * ab003) * d004)) - ((c025 * ab010) * d001)) + ((c026 * ab001
     &) * d001)) + ((c073 * ab007) * d008)) - ((c014 * ab001) * d011))
     & + ((c062 * ab001) * d004)) + ((c063 * ab016) * d002)) - ((c071 *
     &ab002) * d002)) - ((c064 * ab008) * d009)) + ((c056 * ab002) *
     &d020)) - ((c006 * ab002) * d013)) + ((c004 * ab010) * d003)) - ((
     &c005 * ab001) * d003)) - ((c049 * ab007) * d017)) + ((c001 * ab001
     &) * d032)) - ((c050 * ab001) * d022)) - ((c069 * ab010) * d004))
     & + ((c070 * ab001) * d004)) + ((c002 * ab007) * d018)) - ((c050 *
     &ab001) * d033)) + ((c074 * ab001) * d023))
c# 1174 "chelpg.for"
      ep23 = of5 * ((((((((((((- ((c008 * ab004) * d001)) - ((c003 *
     &ab001) * d005)) - ((c009 * ab012) * d001)) + ((c011 * ab004) *
     &d001)) + ((c051 * ab005) * d002)) - ((c006 * ab004) * d003)) + ((
     &c053 * ab004) * d004)) - ((c004 * ab003) * d005)) + ((c005 * ab001
     &) * d005)) + ((c049 * ab002) * d006)) - ((c001 * ab001) * d014))
     & + ((c050 * ab001) * d015))
c# 1178 "chelpg.for"
      ep43 = of5 * ((((((((((((- ((c009 * ab013) * d001)) + ((c007 *
     &ab006) * d002)) + ((c011 * ab002) * d001)) - ((c008 * ab001) *
     &d002)) + ((c007 * ab005) * d005)) - ((c006 * ab004) * d006)) - ((
     &c004 * ab005) * d005)) + ((c002 * ab004) * d006)) + ((c002 * ab002
     &) * d007)) - ((c001 * ab001) * d016)) - ((c055 * ab002) * d004))
     & + ((c050 * ab001) * d013))
c# 1182 "chelpg.for"
      ep53 = of5 * ((((((((((((((- ((c008 * ab004) * d001)) - ((c003 *
     &ab001) * d005)) - ((c009 * ab023) * d001)) + ((c010 * ab004) *
     &d001)) + ((c051 * ab006) * d005)) - ((c052 * ab001) * d005)) - ((
     &c006 * ab004) * d007)) + ((c053 * ab004) * d004)) - ((c004 * ab006
     &) * d005)) + ((c005 * ab001) * d005)) + ((c049 * ab004) * d007))
     & - ((c002 * ab004) * d004)) - ((c001 * ab001) * d034)) + ((c054 *
     &ab001) * d015))
c# 1186 "chelpg.for"
      ep73 = of5 * ((((((((- ((c009 * ab015) * d001)) + ((c007 * ab009)
     & * d002)) + ((c007 * ab005) * d008)) - ((c006 * ab004) * d009)) -
     &((c004 * ab008) * d005)) + ((c002 * ab007) * d006)) + ((c002 *
     &ab002) * d010)) - ((c001 * ab001) * d019))
c# 1189 "chelpg.for"
      ep83 = of5 * ((((((((((((- ((c009 * ab024) * d001)) + ((c007 *
     &ab006) * d008)) + ((c011 * ab007) * d001)) - ((c008 * ab001) *
     &d008)) + ((c007 * ab009) * d005)) - ((c006 * ab004) * d010)) - ((
     &c004 * ab009) * d005)) + ((c002 * ab004) * d010)) + ((c002 * ab007
     &) * d007)) - ((c001 * ab001) * d035)) - ((c055 * ab007) * d004))
     & + ((c050 * ab001) * d018))
c# 1193 "chelpg.for"
      ep93 = of5 * ((((((((((((- ((c008 * ab004) * d001)) - ((c003 *
     &ab001) * d005)) - ((c009 * ab025) * d001)) + ((c011 * ab004) *
     &d001)) + ((c051 * ab009) * d008)) - ((c006 * ab004) * d011)) + ((
     &c053 * ab004) * d004)) - ((c004 * ab010) * d005)) + ((c005 * ab001
     &) * d005)) + ((c049 * ab007) * d010)) - ((c001 * ab001) * d036))
     & + ((c050 * ab001) * d015))
c# 1197 "chelpg.for"
      ep24 = of8 * (((((((((((((((((((((((((((((c018 * ab005) * d001) +
     &((c008 * ab004) * d002)) + ((c008 * ab002) * d005)) + ((c003 *
     &ab001) * d006)) + ((c021 * ab018) * d001)) - ((c024 * ab005) *
     &d001)) - ((c060 * ab012) * d002)) + ((c073 * ab004) * d002)) + ((
     &c013 * ab005) * d003)) - ((c061 * ab005) * d004)) + ((c009 * ab012
     &) * d002)) - ((c011 * ab004) * d002)) - ((c051 * ab005) * d003))
     & + ((c007 * ab005) * d004)) + ((c006 * ab004) * d012)) - ((c075 *
     &ab004) * d013)) + ((c009 * ab011) * d005)) - ((c010 * ab002) *
     &d005)) - ((c051 * ab003) * d006)) + ((c052 * ab001) * d006)) + ((
     &c006 * ab002) * d014)) - ((c053 * ab002) * d015)) + ((c004 * ab003
     &) * d006)) - ((c005 * ab001) * d006)) - ((c049 * ab002) * d014))
     & + ((c002 * ab002) * d015)) + ((c001 * ab001) * d024)) - ((c054 *
     &ab001) * d025))
c# 1204 "chelpg.for"
      ep44 = of8 * (((((((((((((((((((((((((((((((((((c021 * ab019) *
     &d001) - ((c025 * ab006) * d001)) - ((c015 * ab013) * d002)) - ((
     &c025 * ab003) * d001)) + ((c026 * ab001) * d001)) + ((c018 * ab002
     &) * d002)) - ((c015 * ab012) * d005)) + ((c018 * ab004) * d005))
     & + ((c013 * ab005) * d006)) + ((c009 * ab013) * d002)) - ((c007 *
     &ab006) * d003)) + ((c076 * ab006) * d004)) - ((c011 * ab002) *
     &d002)) + ((c008 * ab001) * d003)) - ((c008 * ab001) * d004)) - ((
     &c051 * ab005) * d006)) + ((c006 * ab004) * d014)) - ((c053 * ab004
     &) * d015)) + ((c009 * ab012) * d005)) - ((c011 * ab004) * d005))
     & - ((c007 * ab003) * d007)) + ((c008 * ab001) * d007)) + ((c006 *
     &ab002) * d016)) + ((c076 * ab003) * d004)) - ((c053 * ab002) *
     &d013)) + ((c004 * ab005) * d006)) - ((c002 * ab004) * d014)) + ((
     &c055 * ab004) * d015)) - ((c002 * ab002) * d016)) + ((c001 * ab001
     &) * d026)) - ((c050 * ab001) * d027)) + ((c055 * ab002) * d013))
     & - ((c050 * ab001) * d022)) + ((c074 * ab001) * d023))
c# 1213 "chelpg.for"
      ep54 = of8 * (((((((((((((((((((((((((((((c018 * ab005) * d001) +
     &((c008 * ab004) * d002)) + ((c008 * ab002) * d005)) + ((c003 *
     &ab001) * d006)) + ((c021 * ab026) * d001)) - ((c024 * ab005) *
     &d001)) - ((c060 * ab013) * d005)) + ((c073 * ab002) * d005)) + ((
     &c013 * ab005) * d007)) - ((c061 * ab005) * d004)) + ((c009 * ab023
     &) * d002)) - ((c010 * ab004) * d002)) - ((c051 * ab006) * d006))
     & + ((c052 * ab001) * d006)) + ((c006 * ab004) * d016)) - ((c053 *
     &ab004) * d013)) + ((c009 * ab013) * d005)) - ((c011 * ab002) *
     &d005)) - ((c051 * ab005) * d007)) + ((c007 * ab005) * d004)) + ((
     &c006 * ab002) * d034)) - ((c075 * ab002) * d015)) + ((c004 * ab006
     &) * d006)) - ((c005 * ab001) * d006)) - ((c049 * ab004) * d016))
     & + ((c002 * ab004) * d013)) + ((c001 * ab001) * d037)) - ((c054 *
     &ab001) * d025))
c# 1220 "chelpg.for"
      ep74 = of8 * (((((((((((((((((((((((((c021 * ab021) * d001) - ((
     &c025 * ab009) * d001)) - ((c015 * ab015) * d002)) - ((c015 * ab012
     &) * d008)) + ((c018 * ab004) * d008)) + ((c013 * ab005) * d009))
     & + ((c009 * ab015) * d002)) - ((c007 * ab009) * d003)) + ((c076 *
     &ab009) * d004)) - ((c007 * ab005) * d009)) + ((c006 * ab004) *
     &d017)) - ((c053 * ab004) * d018)) + ((c009 * ab014) * d005)) - ((
     &c011 * ab007) * d005)) - ((c007 * ab008) * d006)) - ((c007 * ab003
     &) * d010)) + ((c008 * ab001) * d010)) + ((c006 * ab002) * d019))
     & + ((c004 * ab008) * d006)) - ((c002 * ab007) * d014)) + ((c055 *
     &ab007) * d015)) - ((c002 * ab002) * d019)) + ((c001 * ab001) *
     &d030)) - ((c050 * ab001) * d031))
c# 1226 "chelpg.for"
      ep84 = of8 * (((((((((((((((((((((((((c021 * ab027) * d001) - ((
     &c015 * ab013) * d008)) - ((c025 * ab008) * d001)) + ((c018 * ab002
     &) * d008)) - ((c015 * ab015) * d005)) + ((c013 * ab005) * d010))
     & + ((c009 * ab024) * d002)) - ((c007 * ab006) * d009)) - ((c011 *
     &ab007) * d002)) + ((c008 * ab001) * d009)) - ((c007 * ab009) *
     &d006)) + ((c006 * ab004) * d019)) + ((c009 * ab015) * d005)) - ((
     &c007 * ab005) * d010)) - ((c007 * ab008) * d007)) + ((c006 * ab002
     &) * d035)) + ((c076 * ab008) * d004)) - ((c053 * ab002) * d018))
     & + ((c004 * ab009) * d006)) - ((c002 * ab004) * d019)) - ((c002 *
     &ab007) * d016)) + ((c001 * ab001) * d038)) + ((c055 * ab007) *
     &d013)) - ((c050 * ab001) * d029))
c# 1232 "chelpg.for"
      ep94 = of8 * (((((((((((((((((((((((((c018 * ab005) * d001) + ((
     &c008 * ab004) * d002)) + ((c008 * ab002) * d005)) + ((c003 * ab001
     &) * d006)) + ((c021 * ab028) * d001)) - ((c025 * ab005) * d001))
     & - ((c060 * ab015) * d008)) + ((c013 * ab005) * d011)) - ((c061 *
     &ab005) * d004)) + ((c009 * ab025) * d002)) - ((c011 * ab004) *
     &d002)) - ((c051 * ab009) * d009)) + ((c006 * ab004) * d020)) - ((
     &c053 * ab004) * d013)) + ((c009 * ab016) * d005)) - ((c011 * ab002
     &) * d005)) - ((c051 * ab008) * d010)) + ((c006 * ab002) * d036))
     & - ((c053 * ab002) * d015)) + ((c004 * ab010) * d006)) - ((c005 *
     &ab001) * d006)) - ((c049 * ab007) * d019)) + ((c001 * ab001) *
     &d039)) - ((c050 * ab001) * d025))
c# 1238 "chelpg.for"
      ep25 = of8 * (((((((((((((((((((((((((((((((((((((c017 * ab001) *
     &d001) + ((c018 * ab006) * d001)) - ((c019 * ab001) * d001)) + ((
     &c052 * ab004) * d005)) + ((c003 * ab001) * d007)) - ((c058 * ab001
     &) * d004)) + ((c011 * ab003) * d001)) - ((c020 * ab001) * d001))
     & - ((c052 * ab002) * d002)) + ((c012 * ab001) * d003)) - ((c059 *
     &ab001) * d004)) + ((c021 * ab019) * d001)) - ((c025 * ab006) *
     &d001)) - ((c060 * ab013) * d002)) + ((c013 * ab006) * d003)) - ((
     &c061 * ab006) * d004)) - ((c025 * ab003) * d001)) + ((c026 * ab001
     &) * d001)) + ((c073 * ab002) * d002)) - ((c014 * ab001) * d003))
     & + ((c062 * ab001) * d004)) + ((c063 * ab012) * d005)) - ((c071 *
     &ab004) * d005)) - ((c064 * ab005) * d006)) + ((c056 * ab004) *
     &d014)) - ((c006 * ab004) * d015)) + ((c004 * ab003) * d007)) - ((
     &c005 * ab001) * d007)) - ((c049 * ab002) * d016)) + ((c001 * ab001
     &) * d026)) - ((c050 * ab001) * d027)) - ((c069 * ab003) * d004))
     & + ((c070 * ab001) * d004)) + ((c002 * ab002) * d013)) - ((c050 *
     &ab001) * d022)) + ((c074 * ab001) * d023))
c# 1247 "chelpg.for"
      ep45 = of8 * (((((((((((((((((((((((((((((c011 * ab005) * d001) -
     &((c008 * ab004) * d002)) - ((c008 * ab002) * d005)) + ((c012 *
     &ab001) * d006)) + ((c021 * ab026) * d001)) - ((c015 * ab023) *
     &d002)) - ((c024 * ab005) * d001)) + ((c016 * ab004) * d002)) - ((
     &c015 * ab013) * d005)) + ((c013 * ab006) * d006)) + ((c018 * ab002
     &) * d005)) - ((c014 * ab001) * d006)) + ((c063 * ab013) * d005))
     & - ((c051 * ab006) * d006)) - ((c071 * ab002) * d005)) + ((c052 *
     &ab001) * d006)) - ((c051 * ab005) * d007)) + ((c056 * ab004) *
     &d016)) + ((c007 * ab005) * d004)) - ((c006 * ab004) * d013)) + ((
     &c004 * ab005) * d007)) - ((c002 * ab004) * d016)) - ((c002 * ab002
     &) * d034)) + ((c001 * ab001) * d037)) + ((c072 * ab002) * d015))
     & - ((c054 * ab001) * d025)) - ((c069 * ab005) * d004)) + ((c055 *
     &ab004) * d013))
c# 1254 "chelpg.for"
      ep55 = of8 * (((((((((((((((((((((((((((((((((((((c017 * ab001) *
     &d001) + ((c018 * ab006) * d001)) - ((c019 * ab001) * d001)) + ((
     &c057 * ab004) * d005)) + ((c003 * ab001) * d007)) - ((c058 * ab001
     &) * d004)) + ((c011 * ab006) * d001)) - ((c020 * ab001) * d001))
     & + ((c012 * ab001) * d007)) - ((c059 * ab001) * d004)) + ((c021 *
     &ab029) * d001)) - ((c022 * ab006) * d001)) - ((c060 * ab023) *
     &d005)) + ((c023 * ab001) * d001)) + ((c038 * ab004) * d005)) + ((
     &c013 * ab006) * d007)) - ((c061 * ab006) * d004)) - ((c014 * ab001
     &) * d007)) + ((c062 * ab001) * d004)) + ((c063 * ab023) * d005))
     & - ((c033 * ab004) * d005)) - ((c064 * ab006) * d007)) + ((c051 *
     &ab006) * d004)) + ((c031 * ab001) * d007)) - ((c052 * ab001) *
     &d004)) + ((c056 * ab004) * d034)) - ((c065 * ab004) * d015)) + ((
     &c004 * ab006) * d007)) - ((c005 * ab001) * d007)) - ((c049 * ab004
     &) * d034)) + ((c066 * ab004) * d015)) + ((c001 * ab001) * d040))
     & - ((c067 * ab001) * d027)) + ((c068 * ab001) * d023)) - ((c069 *
     &ab006) * d004)) + ((c070 * ab001) * d004))
c# 1263 "chelpg.for"
      ep75 = of8 * (((((((((((((((((((((((((c011 * ab008) * d001) - ((
     &c008 * ab007) * d002)) - ((c008 * ab002) * d008)) + ((c012 * ab001
     &) * d009)) + ((c021 * ab027) * d001)) - ((c015 * ab024) * d002))
     & - ((c015 * ab013) * d008)) + ((c013 * ab006) * d009)) - ((c025 *
     &ab008) * d001)) + ((c018 * ab007) * d002)) + ((c018 * ab002) *
     &d008)) - ((c014 * ab001) * d009)) + ((c063 * ab015) * d005)) - ((
     &c051 * ab009) * d006)) - ((c051 * ab005) * d010)) + ((c056 * ab004
     &) * d019)) + ((c004 * ab008) * d007)) - ((c002 * ab007) * d016))
     & - ((c002 * ab002) * d035)) + ((c001 * ab001) * d038)) - ((c069 *
     &ab008) * d004)) + ((c055 * ab007) * d013)) + ((c055 * ab002) *
     &d018)) - ((c050 * ab001) * d029))
c# 1269 "chelpg.for"
      ep85 = of8 * (((((((((((((((((((((((((((((c011 * ab009) * d001) -
     &((c008 * ab004) * d008)) - ((c008 * ab007) * d005)) + ((c012 *
     &ab001) * d010)) + ((c021 * ab030) * d001)) - ((c015 * ab023) *
     &d008)) - ((c024 * ab009) * d001)) + ((c016 * ab004) * d008)) - ((
     &c015 * ab024) * d005)) + ((c013 * ab006) * d010)) + ((c018 * ab007
     &) * d005)) - ((c014 * ab001) * d010)) + ((c063 * ab024) * d005))
     & - ((c051 * ab006) * d010)) - ((c071 * ab007) * d005)) + ((c052 *
     &ab001) * d010)) - ((c051 * ab009) * d007)) + ((c056 * ab004) *
     &d035)) + ((c007 * ab009) * d004)) - ((c006 * ab004) * d018)) + ((
     &c004 * ab009) * d007)) - ((c002 * ab004) * d035)) - ((c002 * ab007
     &) * d034)) + ((c001 * ab001) * d041)) + ((c072 * ab007) * d015))
     & - ((c054 * ab001) * d031)) - ((c069 * ab009) * d004)) + ((c055 *
     &ab004) * d018))
c# 1276 "chelpg.for"
      ep95 = of8 * (((((((((((((((((((((((((((((((((((((c017 * ab001) *
     &d001) + ((c018 * ab006) * d001)) - ((c019 * ab001) * d001)) + ((
     &c052 * ab004) * d005)) + ((c003 * ab001) * d007)) - ((c058 * ab001
     &) * d004)) + ((c011 * ab010) * d001)) - ((c020 * ab001) * d001))
     & - ((c052 * ab007) * d008)) + ((c012 * ab001) * d011)) - ((c059 *
     &ab001) * d004)) + ((c021 * ab031) * d001)) - ((c025 * ab006) *
     &d001)) - ((c060 * ab024) * d008)) + ((c013 * ab006) * d011)) - ((
     &c061 * ab006) * d004)) - ((c025 * ab010) * d001)) + ((c026 * ab001
     &) * d001)) + ((c073 * ab007) * d008)) - ((c014 * ab001) * d011))
     & + ((c062 * ab001) * d004)) + ((c063 * ab025) * d005)) - ((c071 *
     &ab004) * d005)) - ((c064 * ab009) * d010)) + ((c056 * ab004) *
     &d036)) - ((c006 * ab004) * d015)) + ((c004 * ab010) * d007)) - ((
     &c005 * ab001) * d007)) - ((c049 * ab007) * d035)) + ((c001 * ab001
     &) * d042)) - ((c050 * ab001) * d027)) - ((c069 * ab010) * d004))
     & + ((c070 * ab001) * d004)) + ((c002 * ab007) * d018)) - ((c050 *
     &ab001) * d033)) + ((c074 * ab001) * d023))
c# 1285 "chelpg.for"
      ep26 = of5 * ((((((((((((- ((c008 * ab007) * d001)) - ((c003 *
     &ab001) * d008)) - ((c009 * ab014) * d001)) + ((c011 * ab007) *
     &d001)) + ((c051 * ab008) * d002)) - ((c006 * ab007) * d003)) + ((
     &c053 * ab007) * d004)) - ((c004 * ab003) * d008)) + ((c005 * ab001
     &) * d008)) + ((c049 * ab002) * d009)) - ((c001 * ab001) * d017))
     & + ((c050 * ab001) * d018))
c# 1289 "chelpg.for"
      ep46 = of5 * ((((((((- ((c009 * ab015) * d001)) + ((c007 * ab009)
     & * d002)) + ((c007 * ab008) * d005)) - ((c006 * ab007) * d006)) -
     &((c004 * ab005) * d008)) + ((c002 * ab004) * d009)) + ((c002 *
     &ab002) * d010)) - ((c001 * ab001) * d019))
c# 1292 "chelpg.for"
      ep56 = of5 * ((((((((((((- ((c008 * ab007) * d001)) - ((c003 *
     &ab001) * d008)) - ((c009 * ab024) * d001)) + ((c011 * ab007) *
     &d001)) + ((c051 * ab009) * d005)) - ((c006 * ab007) * d007)) + ((
     &c053 * ab007) * d004)) - ((c004 * ab006) * d008)) + ((c005 * ab001
     &) * d008)) + ((c049 * ab004) * d010)) - ((c001 * ab001) * d035))
     & + ((c050 * ab001) * d018))
c# 1296 "chelpg.for"
      ep76 = of5 * ((((((((((((- ((c009 * ab016) * d001)) + ((c007 *
     &ab010) * d002)) + ((c011 * ab002) * d001)) - ((c008 * ab001) *
     &d002)) + ((c007 * ab008) * d008)) - ((c006 * ab007) * d009)) - ((
     &c004 * ab008) * d008)) + ((c002 * ab007) * d009)) + ((c002 * ab002
     &) * d011)) - ((c001 * ab001) * d020)) - ((c055 * ab002) * d004))
     & + ((c050 * ab001) * d013))
c# 1300 "chelpg.for"
      ep86 = of5 * ((((((((((((- ((c009 * ab025) * d001)) + ((c011 *
     &ab004) * d001)) + ((c007 * ab009) * d008)) + ((c007 * ab010) *
     &d005)) - ((c008 * ab001) * d005)) - ((c006 * ab007) * d010)) - ((
     &c004 * ab009) * d008)) + ((c002 * ab004) * d011)) - ((c055 * ab004
     &) * d004)) + ((c002 * ab007) * d010)) - ((c001 * ab001) * d036))
     & + ((c050 * ab001) * d015))
c# 1304 "chelpg.for"
      ep96 = of5 * ((((((((((((((- ((c008 * ab007) * d001)) - ((c003 *
     &ab001) * d008)) - ((c009 * ab032) * d001)) + ((c010 * ab007) *
     &d001)) + ((c051 * ab010) * d008)) - ((c052 * ab001) * d008)) - ((
     &c006 * ab007) * d011)) + ((c053 * ab007) * d004)) - ((c004 * ab010
     &) * d008)) + ((c005 * ab001) * d008)) + ((c049 * ab007) * d011))
     & - ((c002 * ab007) * d004)) - ((c001 * ab001) * d043)) + ((c054 *
     &ab001) * d018))
c# 1308 "chelpg.for"
      ep27 = of8 * (((((((((((((((((((((((((((((c018 * ab008) * d001) +
     &((c008 * ab007) * d002)) + ((c008 * ab002) * d008)) + ((c003 *
     &ab001) * d009)) + ((c021 * ab020) * d001)) - ((c024 * ab008) *
     &d001)) - ((c060 * ab014) * d002)) + ((c073 * ab007) * d002)) + ((
     &c013 * ab008) * d003)) - ((c061 * ab008) * d004)) + ((c009 * ab014
     &) * d002)) - ((c011 * ab007) * d002)) - ((c051 * ab008) * d003))
     & + ((c007 * ab008) * d004)) + ((c006 * ab007) * d012)) - ((c075 *
     &ab007) * d013)) + ((c009 * ab011) * d008)) - ((c010 * ab002) *
     &d008)) - ((c051 * ab003) * d009)) + ((c052 * ab001) * d009)) + ((
     &c006 * ab002) * d017)) - ((c053 * ab002) * d018)) + ((c004 * ab003
     &) * d009)) - ((c005 * ab001) * d009)) - ((c049 * ab002) * d017))
     & + ((c002 * ab002) * d018)) + ((c001 * ab001) * d028)) - ((c054 *
     &ab001) * d029))
c# 1315 "chelpg.for"
      ep47 = of8 * (((((((((((((((((((((((((c021 * ab021) * d001) - ((
     &c025 * ab009) * d001)) - ((c015 * ab015) * d002)) - ((c015 * ab014
     &) * d005)) + ((c018 * ab007) * d005)) + ((c013 * ab008) * d006))
     & + ((c009 * ab015) * d002)) - ((c007 * ab009) * d003)) + ((c076 *
     &ab009) * d004)) - ((c007 * ab008) * d006)) + ((c006 * ab007) *
     &d014)) - ((c053 * ab007) * d015)) + ((c009 * ab012) * d008)) - ((
     &c011 * ab004) * d008)) - ((c007 * ab005) * d009)) - ((c007 * ab003
     &) * d010)) + ((c008 * ab001) * d010)) + ((c006 * ab002) * d019))
     & + ((c004 * ab005) * d009)) - ((c002 * ab004) * d017)) + ((c055 *
     &ab004) * d018)) - ((c002 * ab002) * d019)) + ((c001 * ab001) *
     &d030)) - ((c050 * ab001) * d031))
c# 1321 "chelpg.for"
      ep57 = of8 * (((((((((((((((((((((((((c018 * ab008) * d001) + ((
     &c008 * ab007) * d002)) + ((c008 * ab002) * d008)) + ((c003 * ab001
     &) * d009)) + ((c021 * ab027) * d001)) - ((c025 * ab008) * d001))
     & - ((c060 * ab015) * d005)) + ((c013 * ab008) * d007)) - ((c061 *
     &ab008) * d004)) + ((c009 * ab024) * d002)) - ((c011 * ab007) *
     &d002)) - ((c051 * ab009) * d006)) + ((c006 * ab007) * d016)) - ((
     &c053 * ab007) * d013)) + ((c009 * ab013) * d008)) - ((c011 * ab002
     &) * d008)) - ((c051 * ab005) * d010)) + ((c006 * ab002) * d035))
     & - ((c053 * ab002) * d018)) + ((c004 * ab006) * d009)) - ((c005 *
     &ab001) * d009)) - ((c049 * ab004) * d019)) + ((c001 * ab001) *
     &d038)) - ((c050 * ab001) * d029))
c# 1327 "chelpg.for"
      ep77 = of8 * (((((((((((((((((((((((((((((((((((c021 * ab022) *
     &d001) - ((c025 * ab010) * d001)) - ((c015 * ab016) * d002)) - ((
     &c025 * ab003) * d001)) + ((c026 * ab001) * d001)) + ((c018 * ab002
     &) * d002)) - ((c015 * ab014) * d008)) + ((c018 * ab007) * d008))
     & + ((c013 * ab008) * d009)) + ((c009 * ab016) * d002)) - ((c007 *
     &ab010) * d003)) + ((c076 * ab010) * d004)) - ((c011 * ab002) *
     &d002)) + ((c008 * ab001) * d003)) - ((c008 * ab001) * d004)) - ((
     &c051 * ab008) * d009)) + ((c006 * ab007) * d017)) - ((c053 * ab007
     &) * d018)) + ((c009 * ab014) * d008)) - ((c011 * ab007) * d008))
     & - ((c007 * ab003) * d011)) + ((c008 * ab001) * d011)) + ((c006 *
     &ab002) * d020)) + ((c076 * ab003) * d004)) - ((c053 * ab002) *
     &d013)) + ((c004 * ab008) * d009)) - ((c002 * ab007) * d017)) + ((
     &c055 * ab007) * d018)) - ((c002 * ab002) * d020)) + ((c001 * ab001
     &) * d032)) - ((c050 * ab001) * d033)) + ((c055 * ab002) * d013))
     & - ((c050 * ab001) * d022)) + ((c074 * ab001) * d023))
c# 1336 "chelpg.for"
      ep87 = of8 * (((((((((((((((((((((((((c021 * ab028) * d001) - ((
     &c025 * ab005) * d001)) - ((c015 * ab015) * d008)) - ((c015 * ab016
     &) * d005)) + ((c018 * ab002) * d005)) + ((c013 * ab008) * d010))
     & + ((c009 * ab025) * d002)) - ((c011 * ab004) * d002)) - ((c007 *
     &ab009) * d009)) - ((c007 * ab010) * d006)) + ((c008 * ab001) *
     &d006)) + ((c006 * ab007) * d019)) + ((c009 * ab015) * d008)) - ((
     &c007 * ab005) * d011)) + ((c076 * ab005) * d004)) - ((c007 * ab008
     &) * d010)) + ((c006 * ab002) * d036)) - ((c053 * ab002) * d015))
     & + ((c004 * ab009) * d009)) - ((c002 * ab004) * d020)) + ((c055 *
     &ab004) * d013)) - ((c002 * ab007) * d019)) + ((c001 * ab001) *
     &d039)) - ((c050 * ab001) * d025))
c# 1342 "chelpg.for"
      ep97 = of8 * (((((((((((((((((((((((((((((c018 * ab008) * d001) +
     &((c008 * ab007) * d002)) + ((c008 * ab002) * d008)) + ((c003 *
     &ab001) * d009)) + ((c021 * ab033) * d001)) - ((c024 * ab008) *
     &d001)) - ((c060 * ab016) * d008)) + ((c073 * ab002) * d008)) + ((
     &c013 * ab008) * d011)) - ((c061 * ab008) * d004)) + ((c009 * ab032
     &) * d002)) - ((c010 * ab007) * d002)) - ((c051 * ab010) * d009))
     & + ((c052 * ab001) * d009)) + ((c006 * ab007) * d020)) - ((c053 *
     &ab007) * d013)) + ((c009 * ab016) * d008)) - ((c011 * ab002) *
     &d008)) - ((c051 * ab008) * d011)) + ((c007 * ab008) * d004)) + ((
     &c006 * ab002) * d043)) - ((c075 * ab002) * d018)) + ((c004 * ab010
     &) * d009)) - ((c005 * ab001) * d009)) - ((c049 * ab007) * d020))
     & + ((c002 * ab007) * d013)) + ((c001 * ab001) * d044)) - ((c054 *
     &ab001) * d029))
c# 1349 "chelpg.for"
      ep28 = of8 * (((((((((((((((((((((((((c018 * ab009) * d001) + ((
     &c008 * ab004) * d008)) + ((c008 * ab007) * d005)) + ((c003 * ab001
     &) * d010)) + ((c021 * ab021) * d001)) - ((c025 * ab009) * d001))
     & - ((c060 * ab015) * d002)) + ((c013 * ab009) * d003)) - ((c061 *
     &ab009) * d004)) + ((c009 * ab012) * d008)) - ((c011 * ab004) *
     &d008)) - ((c051 * ab005) * d009)) + ((c006 * ab004) * d017)) - ((
     &c053 * ab004) * d018)) + ((c009 * ab014) * d005)) - ((c011 * ab007
     &) * d005)) - ((c051 * ab008) * d006)) + ((c006 * ab007) * d014))
     & - ((c053 * ab007) * d015)) + ((c004 * ab003) * d010)) - ((c005 *
     &ab001) * d010)) - ((c049 * ab002) * d019)) + ((c001 * ab001) *
     &d030)) - ((c050 * ab001) * d031))
c# 1355 "chelpg.for"
      ep48 = of8 * (((((((((((((((((((((((((c021 * ab027) * d001) - ((
     &c015 * ab024) * d002)) - ((c025 * ab008) * d001)) + ((c018 * ab007
     &) * d002)) - ((c015 * ab015) * d005)) + ((c013 * ab009) * d006))
     & + ((c009 * ab013) * d008)) - ((c007 * ab006) * d009)) - ((c011 *
     &ab002) * d008)) + ((c008 * ab001) * d009)) - ((c007 * ab005) *
     &d010)) + ((c006 * ab004) * d019)) + ((c009 * ab015) * d005)) - ((
     &c007 * ab009) * d006)) - ((c007 * ab008) * d007)) + ((c006 * ab007
     &) * d016)) + ((c076 * ab008) * d004)) - ((c053 * ab007) * d013))
     & + ((c004 * ab005) * d010)) - ((c002 * ab004) * d019)) - ((c002 *
     &ab002) * d035)) + ((c001 * ab001) * d038)) + ((c055 * ab002) *
     &d018)) - ((c050 * ab001) * d029))
c# 1361 "chelpg.for"
      ep58 = of8 * (((((((((((((((((((((((((((((c018 * ab009) * d001) +
     &((c008 * ab004) * d008)) + ((c008 * ab007) * d005)) + ((c003 *
     &ab001) * d010)) + ((c021 * ab030) * d001)) - ((c024 * ab009) *
     &d001)) - ((c060 * ab024) * d005)) + ((c073 * ab007) * d005)) + ((
     &c013 * ab009) * d007)) - ((c061 * ab009) * d004)) + ((c009 * ab023
     &) * d008)) - ((c010 * ab004) * d008)) - ((c051 * ab006) * d010))
     & + ((c052 * ab001) * d010)) + ((c006 * ab004) * d035)) - ((c053 *
     &ab004) * d018)) + ((c009 * ab024) * d005)) - ((c011 * ab007) *
     &d005)) - ((c051 * ab009) * d007)) + ((c007 * ab009) * d004)) + ((
     &c006 * ab007) * d034)) - ((c075 * ab007) * d015)) + ((c004 * ab006
     &) * d010)) - ((c005 * ab001) * d010)) - ((c049 * ab004) * d035))
     & + ((c002 * ab004) * d018)) + ((c001 * ab001) * d041)) - ((c054 *
     &ab001) * d031))
c# 1368 "chelpg.for"
      ep78 = of8 * (((((((((((((((((((((((((c021 * ab028) * d001) - ((
     &c015 * ab025) * d002)) - ((c025 * ab005) * d001)) + ((c018 * ab004
     &) * d002)) - ((c015 * ab015) * d008)) + ((c013 * ab009) * d009))
     & + ((c009 * ab015) * d008)) - ((c007 * ab009) * d009)) - ((c007 *
     &ab005) * d011)) + ((c006 * ab004) * d020)) + ((c076 * ab005) *
     &d004)) - ((c053 * ab004) * d013)) + ((c009 * ab016) * d005)) - ((
     &c007 * ab010) * d006)) - ((c011 * ab002) * d005)) + ((c008 * ab001
     &) * d006)) - ((c007 * ab008) * d010)) + ((c006 * ab007) * d019))
     & + ((c004 * ab008) * d010)) - ((c002 * ab007) * d019)) - ((c002 *
     &ab002) * d036)) + ((c001 * ab001) * d039)) + ((c055 * ab002) *
     &d015)) - ((c050 * ab001) * d025))
c# 1374 "chelpg.for"
      ep88 = of8 * (((((((((((((((((((((((((((((((((((c021 * ab031) *
     &d001) - ((c025 * ab006) * d001)) - ((c015 * ab024) * d008)) - ((
     &c025 * ab010) * d001)) + ((c026 * ab001) * d001)) + ((c018 * ab007
     &) * d008)) - ((c015 * ab025) * d005)) + ((c018 * ab004) * d005))
     & + ((c013 * ab009) * d010)) + ((c009 * ab024) * d008)) - ((c007 *
     &ab006) * d011)) + ((c076 * ab006) * d004)) - ((c011 * ab007) *
     &d008)) + ((c008 * ab001) * d011)) - ((c008 * ab001) * d004)) - ((
     &c051 * ab009) * d010)) + ((c006 * ab004) * d036)) - ((c053 * ab004
     &) * d015)) + ((c009 * ab025) * d005)) - ((c011 * ab004) * d005))
     & - ((c007 * ab010) * d007)) + ((c008 * ab001) * d007)) + ((c006 *
     &ab007) * d035)) + ((c076 * ab010) * d004)) - ((c053 * ab007) *
     &d018)) + ((c004 * ab009) * d010)) - ((c002 * ab004) * d036)) + ((
     &c055 * ab004) * d015)) - ((c002 * ab007) * d035)) + ((c001 * ab001
     &) * d042)) - ((c050 * ab001) * d027)) + ((c055 * ab007) * d018))
     & - ((c050 * ab001) * d033)) + ((c074 * ab001) * d023))
c# 1383 "chelpg.for"
      ep98 = of8 * (((((((((((((((((((((((((((((c018 * ab009) * d001) +
     &((c008 * ab004) * d008)) + ((c008 * ab007) * d005)) + ((c003 *
     &ab001) * d010)) + ((c021 * ab034) * d001)) - ((c024 * ab009) *
     &d001)) - ((c060 * ab025) * d008)) + ((c073 * ab004) * d008)) + ((
     &c013 * ab009) * d011)) - ((c061 * ab009) * d004)) + ((c009 * ab025
     &) * d008)) - ((c011 * ab004) * d008)) - ((c051 * ab009) * d011))
     & + ((c007 * ab009) * d004)) + ((c006 * ab004) * d043)) - ((c075 *
     &ab004) * d018)) + ((c009 * ab032) * d005)) - ((c010 * ab007) *
     &d005)) - ((c051 * ab010) * d010)) + ((c052 * ab001) * d010)) + ((
     &c006 * ab007) * d036)) - ((c053 * ab007) * d015)) + ((c004 * ab010
     &) * d010)) - ((c005 * ab001) * d010)) - ((c049 * ab007) * d036))
     & + ((c002 * ab007) * d015)) + ((c001 * ab001) * d045)) - ((c054 *
     &ab001) * d031))
c# 1390 "chelpg.for"
      ep29 = of8 * (((((((((((((((((((((((((((((((((((((c017 * ab001) *
     &d001) + ((c018 * ab010) * d001)) - ((c019 * ab001) * d001)) + ((
     &c052 * ab007) * d008)) + ((c003 * ab001) * d011)) - ((c058 * ab001
     &) * d004)) + ((c011 * ab003) * d001)) - ((c020 * ab001) * d001))
     & - ((c052 * ab002) * d002)) + ((c012 * ab001) * d003)) - ((c059 *
     &ab001) * d004)) + ((c021 * ab022) * d001)) - ((c025 * ab010) *
     &d001)) - ((c060 * ab016) * d002)) + ((c013 * ab010) * d003)) - ((
     &c061 * ab010) * d004)) - ((c025 * ab003) * d001)) + ((c026 * ab001
     &) * d001)) + ((c073 * ab002) * d002)) - ((c014 * ab001) * d003))
     & + ((c062 * ab001) * d004)) + ((c063 * ab014) * d008)) - ((c071 *
     &ab007) * d008)) - ((c064 * ab008) * d009)) + ((c056 * ab007) *
     &d017)) - ((c006 * ab007) * d018)) + ((c004 * ab003) * d011)) - ((
     &c005 * ab001) * d011)) - ((c049 * ab002) * d020)) + ((c001 * ab001
     &) * d032)) - ((c050 * ab001) * d033)) - ((c069 * ab003) * d004))
     & + ((c070 * ab001) * d004)) + ((c002 * ab002) * d013)) - ((c050 *
     &ab001) * d022)) + ((c074 * ab001) * d023))
c# 1399 "chelpg.for"
      ep49 = of8 * (((((((((((((((((((((((((c011 * ab005) * d001) - ((
     &c008 * ab004) * d002)) - ((c008 * ab002) * d005)) + ((c012 * ab001
     &) * d006)) + ((c021 * ab028) * d001)) - ((c015 * ab025) * d002))
     & - ((c015 * ab016) * d005)) + ((c013 * ab010) * d006)) - ((c025 *
     &ab005) * d001)) + ((c018 * ab004) * d002)) + ((c018 * ab002) *
     &d005)) - ((c014 * ab001) * d006)) + ((c063 * ab015) * d008)) - ((
     &c051 * ab009) * d009)) - ((c051 * ab008) * d010)) + ((c056 * ab007
     &) * d019)) + ((c004 * ab005) * d011)) - ((c002 * ab004) * d020))
     & - ((c002 * ab002) * d036)) + ((c001 * ab001) * d039)) - ((c069 *
     &ab005) * d004)) + ((c055 * ab004) * d013)) + ((c055 * ab002) *
     &d015)) - ((c050 * ab001) * d025))
c# 1405 "chelpg.for"
      ep59 = of8 * (((((((((((((((((((((((((((((((((((((c017 * ab001) *
     &d001) + ((c018 * ab010) * d001)) - ((c019 * ab001) * d001)) + ((
     &c052 * ab007) * d008)) + ((c003 * ab001) * d011)) - ((c058 * ab001
     &) * d004)) + ((c011 * ab006) * d001)) - ((c020 * ab001) * d001))
     & - ((c052 * ab004) * d005)) + ((c012 * ab001) * d007)) - ((c059 *
     &ab001) * d004)) + ((c021 * ab031) * d001)) - ((c025 * ab010) *
     &d001)) - ((c060 * ab025) * d005)) + ((c013 * ab010) * d007)) - ((
     &c061 * ab010) * d004)) - ((c025 * ab006) * d001)) + ((c026 * ab001
     &) * d001)) + ((c073 * ab004) * d005)) - ((c014 * ab001) * d007))
     & + ((c062 * ab001) * d004)) + ((c063 * ab024) * d008)) - ((c071 *
     &ab007) * d008)) - ((c064 * ab009) * d010)) + ((c056 * ab007) *
     &d035)) - ((c006 * ab007) * d018)) + ((c004 * ab006) * d011)) - ((
     &c005 * ab001) * d011)) - ((c049 * ab004) * d036)) + ((c001 * ab001
     &) * d042)) - ((c050 * ab001) * d033)) - ((c069 * ab006) * d004))
     & + ((c070 * ab001) * d004)) + ((c002 * ab004) * d015)) - ((c050 *
     &ab001) * d027)) + ((c074 * ab001) * d023))
c# 1414 "chelpg.for"
      ep79 = of8 * (((((((((((((((((((((((((((((c011 * ab008) * d001) -
     &((c008 * ab007) * d002)) - ((c008 * ab002) * d008)) + ((c012 *
     &ab001) * d009)) + ((c021 * ab033) * d001)) - ((c015 * ab032) *
     &d002)) - ((c024 * ab008) * d001)) + ((c016 * ab007) * d002)) - ((
     &c015 * ab016) * d008)) + ((c013 * ab010) * d009)) + ((c018 * ab002
     &) * d008)) - ((c014 * ab001) * d009)) + ((c063 * ab016) * d008))
     & - ((c051 * ab010) * d009)) - ((c071 * ab002) * d008)) + ((c052 *
     &ab001) * d009)) - ((c051 * ab008) * d011)) + ((c056 * ab007) *
     &d020)) + ((c007 * ab008) * d004)) - ((c006 * ab007) * d013)) + ((
     &c004 * ab008) * d011)) - ((c002 * ab007) * d020)) - ((c002 * ab002
     &) * d043)) + ((c001 * ab001) * d044)) + ((c072 * ab002) * d018))
     & - ((c054 * ab001) * d029)) - ((c069 * ab008) * d004)) + ((c055 *
     &ab007) * d013))
c# 1421 "chelpg.for"
      ep89 = of8 * (((((((((((((((((((((((((((((c011 * ab009) * d001) -
     &((c008 * ab004) * d008)) - ((c008 * ab007) * d005)) + ((c012 *
     &ab001) * d010)) + ((c021 * ab034) * d001)) - ((c024 * ab009) *
     &d001)) - ((c015 * ab025) * d008)) - ((c015 * ab032) * d005)) + ((
     &c016 * ab007) * d005)) + ((c013 * ab010) * d010)) + ((c018 * ab004
     &) * d008)) - ((c014 * ab001) * d010)) + ((c063 * ab025) * d008))
     & - ((c071 * ab004) * d008)) - ((c051 * ab009) * d011)) + ((c007 *
     &ab009) * d004)) - ((c051 * ab010) * d010)) + ((c052 * ab001) *
     &d010)) + ((c056 * ab007) * d036)) - ((c006 * ab007) * d015)) + ((
     &c004 * ab009) * d011)) - ((c002 * ab004) * d043)) + ((c072 * ab004
     &) * d018)) - ((c002 * ab007) * d036)) + ((c001 * ab001) * d045))
     & - ((c054 * ab001) * d031)) - ((c069 * ab009) * d004)) + ((c055 *
     &ab007) * d015))
c     ******************************************************************
c SPDSTV
c     *                                PD                              *
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 1428 "chelpg.for"
      ep99 = of8 * (((((((((((((((((((((((((((((((((((((c017 * ab001) *
     &d001) + ((c018 * ab010) * d001)) - ((c019 * ab001) * d001)) + ((
     &c057 * ab007) * d008)) + ((c003 * ab001) * d011)) - ((c058 * ab001
     &) * d004)) + ((c011 * ab010) * d001)) - ((c020 * ab001) * d001))
     & + ((c012 * ab001) * d011)) - ((c059 * ab001) * d004)) + ((c021 *
     &ab035) * d001)) - ((c022 * ab010) * d001)) - ((c060 * ab032) *
     &d008)) + ((c023 * ab001) * d001)) + ((c038 * ab007) * d008)) + ((
     &c013 * ab010) * d011)) - ((c061 * ab010) * d004)) - ((c014 * ab001
     &) * d011)) + ((c062 * ab001) * d004)) + ((c063 * ab032) * d008))
     & - ((c033 * ab007) * d008)) - ((c064 * ab010) * d011)) + ((c051 *
     &ab010) * d004)) + ((c031 * ab001) * d011)) - ((c052 * ab001) *
     &d004)) + ((c056 * ab007) * d043)) - ((c065 * ab007) * d018)) + ((
     &c004 * ab010) * d011)) - ((c005 * ab001) * d011)) - ((c049 * ab007
     &) * d043)) + ((c066 * ab007) * d018)) + ((c001 * ab001) * d046))
     & - ((c067 * ab001) * d033)) + ((c068 * ab001) * d023)) - ((c069 *
     &ab010) * d004)) + ((c070 * ab001) * d004))
c# 1440 "chelpg.for"
 3242 continue
      ep12 = of7 * (((((((((((((((c008 * ab002) * d001) - ((c012 * ab001
     &) * d002)) + ((c015 * ab011) * d001)) - ((c016 * ab002) * d001))
     & - ((c013 * ab003) * d002)) + ((c014 * ab001) * d002)) + ((c051 *
     &ab003) * d002)) - ((c052 * ab001) * d002)) - ((c056 * ab002) *
     &d003)) + ((c006 * ab002) * d004)) + ((c002 * ab002) * d003)) - ((
     &c001 * ab001) * d012)) + ((c054 * ab001) * d013)) - ((c055 * ab002
     &) * d004))
c# 1445 "chelpg.for"
      ep32 = of7 * (((((((((((((c008 * ab004) * d001) - ((c012 * ab001)
     & * d005)) + ((c015 * ab012) * d001)) - ((c013 * ab003) * d005)) -
     &((c018 * ab004) * d001)) + ((c014 * ab001) * d005)) + ((c051 *
     &ab005) * d002)) - ((c056 * ab002) * d006)) + ((c002 * ab004) *
     &d003)) - ((c001 * ab001) * d014)) - ((c055 * ab004) * d004)) + ((
     &c050 * ab001) * d015))
c# 1449 "chelpg.for"
      ep62 = of7 * (((((((((((((c008 * ab007) * d001) - ((c012 * ab001)
     & * d008)) + ((c015 * ab014) * d001)) - ((c013 * ab003) * d008)) -
     &((c018 * ab007) * d001)) + ((c014 * ab001) * d008)) + ((c051 *
     &ab008) * d002)) - ((c056 * ab002) * d009)) + ((c002 * ab007) *
     &d003)) - ((c001 * ab001) * d017)) - ((c055 * ab007) * d004)) + ((
     &c050 * ab001) * d018))
c# 1453 "chelpg.for"
      ep14 = of7 * (((((((((((((c015 * ab012) * d001) - ((c018 * ab004)
     & * d001)) - ((c013 * ab005) * d002)) + ((c007 * ab005) * d002)) -
     &((c006 * ab004) * d003)) + ((c053 * ab004) * d004)) + ((c007 *
     &ab003) * d005)) - ((c008 * ab001) * d005)) - ((c006 * ab002) *
     &d006)) + ((c002 * ab002) * d006)) - ((c001 * ab001) * d014)) + ((
     &c050 * ab001) * d015))
c# 1457 "chelpg.for"
      ep34 = of7 * (((((((((((((c015 * ab013) * d001) - ((c018 * ab002)
     & * d001)) - ((c013 * ab005) * d005)) + ((c007 * ab006) * d002)) -
     &((c008 * ab001) * d002)) - ((c006 * ab004) * d006)) + ((c007 *
     &ab005) * d005)) - ((c006 * ab002) * d007)) + ((c053 * ab002) *
     &d004)) + ((c002 * ab004) * d006)) - ((c001 * ab001) * d016)) + ((
     &c050 * ab001) * d013))
c# 1461 "chelpg.for"
      ep64 = of7 * (((((((((c015 * ab015) * d001) - ((c013 * ab005) *
     &d008)) + ((c007 * ab009) * d002)) - ((c006 * ab004) * d009)) + ((
     &c007 * ab008) * d005)) - ((c006 * ab002) * d010)) + ((c002 * ab007
     &) * d006)) - ((c001 * ab001) * d019))
c# 1464 "chelpg.for"
      ep15 = of7 * (((((((((((((c008 * ab002) * d001) - ((c012 * ab001)
     & * d002)) + ((c015 * ab013) * d001)) - ((c013 * ab006) * d002)) -
     &((c018 * ab002) * d001)) + ((c014 * ab001) * d002)) + ((c051 *
     &ab005) * d005)) - ((c056 * ab004) * d006)) + ((c002 * ab002) *
     &d007)) - ((c001 * ab001) * d016)) - ((c055 * ab002) * d004)) + ((
     &c050 * ab001) * d013))
c# 1468 "chelpg.for"
      ep35 = of7 * (((((((((((((((c008 * ab004) * d001) - ((c012 * ab001
     &) * d005)) + ((c015 * ab023) * d001)) - ((c016 * ab004) * d001))
     & - ((c013 * ab006) * d005)) + ((c014 * ab001) * d005)) + ((c051 *
     &ab006) * d005)) - ((c052 * ab001) * d005)) - ((c056 * ab004) *
     &d007)) + ((c006 * ab004) * d004)) + ((c002 * ab004) * d007)) - ((
     &c001 * ab001) * d034)) + ((c054 * ab001) * d015)) - ((c055 * ab004
     &) * d004))
c# 1472 "chelpg.for"
      ep65 = of7 * (((((((((((((c008 * ab007) * d001) - ((c012 * ab001)
     & * d008)) + ((c015 * ab024) * d001)) - ((c013 * ab006) * d008)) -
     &((c018 * ab007) * d001)) + ((c014 * ab001) * d008)) + ((c051 *
     &ab009) * d005)) - ((c056 * ab004) * d010)) + ((c002 * ab007) *
     &d007)) - ((c001 * ab001) * d035)) - ((c055 * ab007) * d004)) + ((
     &c050 * ab001) * d018))
c# 1476 "chelpg.for"
      ep17 = of7 * (((((((((((((c015 * ab014) * d001) - ((c018 * ab007)
     & * d001)) - ((c013 * ab008) * d002)) + ((c007 * ab008) * d002)) -
     &((c006 * ab007) * d003)) + ((c053 * ab007) * d004)) + ((c007 *
     &ab003) * d008)) - ((c008 * ab001) * d008)) - ((c006 * ab002) *
     &d009)) + ((c002 * ab002) * d009)) - ((c001 * ab001) * d017)) + ((
     &c050 * ab001) * d018))
c# 1480 "chelpg.for"
      ep37 = of7 * (((((((((c015 * ab015) * d001) - ((c013 * ab008) *
     &d005)) + ((c007 * ab009) * d002)) - ((c006 * ab007) * d006)) + ((
     &c007 * ab005) * d008)) - ((c006 * ab002) * d010)) + ((c002 * ab004
     &) * d009)) - ((c001 * ab001) * d019))
c# 1483 "chelpg.for"
      ep67 = of7 * (((((((((((((c015 * ab016) * d001) - ((c018 * ab002)
     & * d001)) - ((c013 * ab008) * d008)) + ((c007 * ab010) * d002)) -
     &((c008 * ab001) * d002)) - ((c006 * ab007) * d009)) + ((c007 *
     &ab008) * d008)) - ((c006 * ab002) * d011)) + ((c053 * ab002) *
     &d004)) + ((c002 * ab007) * d009)) - ((c001 * ab001) * d020)) + ((
     &c050 * ab001) * d013))
c# 1487 "chelpg.for"
      ep18 = of7 * (((((((((c015 * ab015) * d001) - ((c013 * ab009) *
     &d002)) + ((c007 * ab005) * d008)) - ((c006 * ab004) * d009)) + ((
     &c007 * ab008) * d005)) - ((c006 * ab007) * d006)) + ((c002 * ab002
     &) * d010)) - ((c001 * ab001) * d019))
c# 1490 "chelpg.for"
      ep38 = of7 * (((((((((((((c015 * ab024) * d001) - ((c018 * ab007)
     & * d001)) - ((c013 * ab009) * d005)) + ((c007 * ab006) * d008)) -
     &((c008 * ab001) * d008)) - ((c006 * ab004) * d010)) + ((c007 *
     &ab009) * d005)) - ((c006 * ab007) * d007)) + ((c053 * ab007) *
     &d004)) + ((c002 * ab004) * d010)) - ((c001 * ab001) * d035)) + ((
     &c050 * ab001) * d018))
c# 1494 "chelpg.for"
      ep68 = of7 * (((((((((((((c015 * ab025) * d001) - ((c018 * ab004)
     & * d001)) - ((c013 * ab009) * d008)) + ((c007 * ab009) * d008)) -
     &((c006 * ab004) * d011)) + ((c053 * ab004) * d004)) + ((c007 *
     &ab010) * d005)) - ((c008 * ab001) * d005)) - ((c006 * ab007) *
     &d010)) + ((c002 * ab007) * d010)) - ((c001 * ab001) * d036)) + ((
     &c050 * ab001) * d015))
c# 1498 "chelpg.for"
      ep19 = of7 * (((((((((((((c008 * ab002) * d001) - ((c012 * ab001)
     & * d002)) + ((c015 * ab016) * d001)) - ((c013 * ab010) * d002)) -
     &((c018 * ab002) * d001)) + ((c014 * ab001) * d002)) + ((c051 *
     &ab008) * d008)) - ((c056 * ab007) * d009)) + ((c002 * ab002) *
     &d011)) - ((c001 * ab001) * d020)) - ((c055 * ab002) * d004)) + ((
     &c050 * ab001) * d013))
c# 1502 "chelpg.for"
      ep39 = of7 * (((((((((((((c008 * ab004) * d001) - ((c012 * ab001)
     & * d005)) + ((c015 * ab025) * d001)) - ((c013 * ab010) * d005)) -
     &((c018 * ab004) * d001)) + ((c014 * ab001) * d005)) + ((c051 *
     &ab009) * d008)) - ((c056 * ab007) * d010)) + ((c002 * ab004) *
     &d011)) - ((c001 * ab001) * d036)) - ((c055 * ab004) * d004)) + ((
     &c050 * ab001) * d015))
c     ******************************************************************
c SPDSTV
c     *                                SD                              *
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 1506 "chelpg.for"
      ep69 = of7 * (((((((((((((((c008 * ab007) * d001) - ((c012 * ab001
     &) * d008)) + ((c015 * ab032) * d001)) - ((c016 * ab007) * d001))
     & - ((c013 * ab010) * d008)) + ((c014 * ab001) * d008)) + ((c051 *
     &ab010) * d008)) - ((c052 * ab001) * d008)) - ((c056 * ab007) *
     &d011)) + ((c006 * ab007) * d004)) + ((c002 * ab007) * d011)) - ((
     &c001 * ab001) * d043)) + ((c054 * ab001) * d018)) - ((c055 * ab007
     &) * d004))
 3260 continue
      ep02 = of6 * (((((((c012 * ab001) * d001) + ((c013 * ab003) * d001
     &)) - ((c014 * ab001) * d001)) + ((c056 * ab002) * d002)) + ((c001
     & * ab001) * d003)) - ((c050 * ab001) * d004))
c# 1516 "chelpg.for"
      ep04 = of6 * (((((c013 * ab005) * d001) + ((c006 * ab004) * d002))
     & + ((c006 * ab002) * d005)) + ((c001 * ab001) * d006))
      ep05 = of6 * (((((((c012 * ab001) * d001) + ((c013 * ab006) * d001
     &)) - ((c014 * ab001) * d001)) + ((c056 * ab004) * d005)) + ((c001
     & * ab001) * d007)) - ((c050 * ab001) * d004))
c# 1520 "chelpg.for"
      ep07 = of6 * (((((c013 * ab008) * d001) + ((c006 * ab007) * d002))
     & + ((c006 * ab002) * d008)) + ((c001 * ab001) * d009))
      ep08 = of6 * (((((c013 * ab009) * d001) + ((c006 * ab004) * d008))
     & + ((c006 * ab007) * d005)) + ((c001 * ab001) * d010))
      ep09 = of6 * (((((((c012 * ab001) * d001) + ((c013 * ab010) * d001
     &)) - ((c014 * ab001) * d001)) + ((c056 * ab007) * d008)) + ((c001
     & * ab001) * d011)) - ((c050 * ab001) * d004))
c     ******************************************************************
c SPDSTV
c     *                                PP                              *
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 1526 "chelpg.for"
      if (itype - 6) 3261, 3262, 3261
c# 1530 "chelpg.for"
 3261 continue
      ep10 = of1 * (((c002 * ab002) * d001) - ((c001 * ab001) * d002))
      ep30 = of1 * (((c002 * ab004) * d001) - ((c001 * ab001) * d005))
      ep60 = of1 * (((c002 * ab007) * d001) - ((c001 * ab001) * d008))
      ep11 = of4 * ((((((- ((c007 * ab003) * d001)) + ((c008 * ab001) *
     &d001)) + ((c006 * ab002) * d002)) - ((c002 * ab002) * d002)) + ((
     &c001 * ab001) * d003)) - ((c050 * ab001) * d004))
c# 1536 "chelpg.for"
      ep31 = of4 * ((((- ((c007 * ab005) * d001)) + ((c006 * ab002) *
     &d005)) - ((c002 * ab004) * d002)) + ((c001 * ab001) * d006))
      ep61 = of4 * ((((- ((c007 * ab008) * d001)) + ((c006 * ab002) *
     &d008)) - ((c002 * ab007) * d002)) + ((c001 * ab001) * d009))
      ep13 = of4 * ((((- ((c007 * ab005) * d001)) + ((c006 * ab004) *
     &d002)) - ((c002 * ab002) * d005)) + ((c001 * ab001) * d006))
      ep33 = of4 * ((((((- ((c007 * ab006) * d001)) + ((c008 * ab001) *
     &d001)) + ((c006 * ab004) * d005)) - ((c002 * ab004) * d005)) + ((
     &c001 * ab001) * d007)) - ((c050 * ab001) * d004))
c# 1544 "chelpg.for"
      ep63 = of4 * ((((- ((c007 * ab009) * d001)) + ((c006 * ab004) *
     &d008)) - ((c002 * ab007) * d005)) + ((c001 * ab001) * d010))
      ep16 = of4 * ((((- ((c007 * ab008) * d001)) + ((c006 * ab007) *
     &d002)) - ((c002 * ab002) * d008)) + ((c001 * ab001) * d009))
      ep36 = of4 * ((((- ((c007 * ab009) * d001)) + ((c006 * ab007) *
     &d005)) - ((c002 * ab004) * d008)) + ((c001 * ab001) * d010))
      ep66 = of4 * ((((((- ((c007 * ab010) * d001)) + ((c008 * ab001) *
     &d001)) + ((c006 * ab007) * d008)) - ((c002 * ab007) * d008)) + ((
     &c001 * ab001) * d011)) - ((c050 * ab001) * d004))
c# 1552 "chelpg.for"
 3262 continue
      do 2137 i = 1, 100
 2137 epn(i) = epn(i) + eep(i)
c     ******************************************************************
c SPDSTV
c     END OF LOOP OVER GAUSSIANS
c SPDSTV
c     STORE IN ARRAYS
c SPDSTV
c     ******************************************************************
c SPDSTV
c# 1555 "chelpg.for"
  105 continue
c# 1560 "chelpg.for"
      intc = 0
      do 500 j = 1, 10
      r3b = renorm(j)
      do 500 i = 1, 10
      r3a = r3b * renorm(i)
      intc = intc + 1
  500 epn(intc) = (epn(intc) * r3a) * symfac
      call reduc1(epn, lamax, lbmax, i6to5)
      call matfil(h, epn, aos, shellt, inew, jnew, lamax, lbmax, la, lb)
c
c
c        REFORMAT COMMON /B/ AND THE H ARRAY IF THIS BASIS CONTAINS
c
c        P ONLY SHELLS
c
c
c
c# 1569 "chelpg.for"
 1000 continue
c# 1574 "chelpg.for"
      if (ipo(4) .eq. 0) goto 1285
      write(unit=iout, fmt=*) 'DEBUG OF UNSTAR'
      call linout(h, nbasis, 0, 0)
c
c
c# 1577 "chelpg.for"
 1285 continue
c
c
c# 1579 "chelpg.for"
      call unstar(nbasis, shellt, shellc, aos, nshell, h, nostar)
c# 1581 "chelpg.for"
      if (ipo(4) .eq. 0) goto 1500
      write(unit=iout, fmt=2010)
      call linout(h, nbasis, 0, 0)
 1500 continue
      return
      end
      subroutine inv(a, n, is, iad1, iad2, d, mdm)
c     ******************************************************************
c
c     INVERSION OF SQUARE MATRIX A BY MEANS OF THE GAUSS-JORDAN
c
c     ALGORITHM
c
c
c
c     APRIL 72/RS9B
c
c     ******************************************************************
c
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
c
c
      dimension a(mdm, mdm), is(2, mdm), iad1(mdm), iad2(mdm), d(mdm)
c
c
      common /io/ in, iout, ipunch
c
c
c     ******************************************************************
c
 2000 format(37h WARNING FROM INV: MATRIX IS SINGULAR)
      data zero / 0.0d0 /
      data one / 1.0d0 /
      data small / 1.0d-20 /
c# 1603 "chelpg.for"
      do 1 l = 1, n
      is(1,l) = 0
    1 is(2,l) = 0
      do 9 ima = 1, n
      b = zero
      do 2 l = 1, n
      do 2 m = 1, n
      if ((is(1,l) .eq. 1) .or. (is(2,m) .eq. 1)) goto 2
      e = abs(a(l,m))
      if (e .lt. b) goto 8
      i = l
      k = m
    8 b = Amax1(b,e)
    2 continue
      is(1,i) = 1
      is(2,k) = 1
      iad1(k) = i
      iad2(i) = k
c.....PIVOT
c
c# 1621 "chelpg.for"
      b = a(i,k)
c# 1623 "chelpg.for"
      if (abs(b) .lt. small) goto 20
      a(i,k) = one / b
      do 6 l = 1, n
c.....KELLERZEILE
c
c# 1626 "chelpg.for"
      if (l .eq. k) goto 6
c# 1628 "chelpg.for"
      a(i,l) = - (a(i,l) / b)
    6 continue
      do 5 l = 1, n
      do 5 m = 1, n
c.....RECHTECK-REGEL
c
c# 1632 "chelpg.for"
      if ((l .eq. i) .or. (m .eq. k)) goto 5
c# 1634 "chelpg.for"
      a(l,m) = a(l,m) + (a(l,k) * a(i,m))
    5 continue
      do 11 l = 1, n
c.....PIVOT-SPALTE
c
c# 1637 "chelpg.for"
      if (l .eq. i) goto 11
c# 1639 "chelpg.for"
      a(l,k) = a(l,k) / b
   11 continue
c.....PERMUTATION DER ZEILEN, UM DIE NATUERLICHE ORDNUNG WIEDER HERZUSTE
c
c# 1641 "chelpg.for"
    9 continue
c# 1643 "chelpg.for"
      do 15 l = 1, n
      do 13 j = 1, n
      k = iad1(j)
   13 d(j) = a(k,l)
      do 14 j = 1, n
   14 a(j,l) = d(j)
c.....PERMUTATION DER SPALTEN
c
c# 1649 "chelpg.for"
   15 continue
c# 1651 "chelpg.for"
      do 16 l = 1, n
      do 17 j = 1, n
      k = iad2(j)
   17 d(j) = a(l,k)
      do 18 j = 1, n
   18 a(l,j) = d(j)
   16 continue
c
c
c     ERROR EXIT: MATRIX IS SINGULAR
c
c# 1658 "chelpg.for"
      return
c# 1661 "chelpg.for"
   20 write(unit=iout, fmt=2000)
      stop 'INV IN POLAR'
      end
c
c
c     GENERAL LINEAR MATRIX OUTPUT ROUTINE
c
c
c
c     KEY=0  MATRIX SYMMETRIC
c
c     KEY=1  MATRIX SQUARE ASYMMETRIC
c
c
c
c     IZERO=0  ZERO MATRIX ELEMENTS LESS THAN CUTOFF
c
c     IZERO=1  DO NOT ZERO MATRIX ELEMENTS
c
c
c
c     CUTOFF=1.0E-06
c
c
c
      subroutine linout(x, n, key, izero)
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
c
c
      common /io/ in, iout, ipunch
c
c
      dimension s(9), x(1)
c
c
c
c
c
c
      data cutoff / 1.0e-06 /
      data zero / 0.0e0 /
c# 1684 "chelpg.for"
      ia(i) = (i * (i - 1)) / 2
c# 1687 "chelpg.for"
      ilower = 1
  100 iupper = min0(ilower + 8,n)
      irange = min0((iupper - ilower) + 1,9)
      write(unit=iout, fmt=9000) (j,j = ilower, iupper)
      write(unit=iout, fmt=9010)
      do 160 i = 1, n
      k = 1
      do 150 j = ilower, iupper
      if (key) 110, 120, 110
  110 ij = (n * (j - 1)) + i
      goto 140
  120 ij = ia(i) + j
      if (i - j) 130, 140, 140
  130 ij = ia(j) + i
  140 s(k) = x(ij)
      if ((izero .eq. 0) .and. (abs(s(k)) .le. cutoff)) s(k) = zero
  150 k = k + 1
  160 write(unit=iout, fmt=9020) i, (s(j),j = 1, irange)
      write(unit=iout, fmt=9010)
      ilower = ilower + 9
      if (n - iupper) 170, 170, 100
  170 return
 9000 format(12x,8(i3,11x),i3)
 9010 format(/)
 9020 format(1x,i3,2x,9e14.6)
      end
c
c
c     GAUSSIAN 77/UCI
c
c
c
      subroutine matfil(a, aa, aos, shellt, inew, jnew, lamax, lbmax, la
     &, lb)
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
c
c
      integer aos(1), shellt(1)
c
c
      dimension a(1), aa(1)
c
c
c# 1722 "chelpg.for"
      lind(i) = (i * (i - 1)) / 2
c# 1724 "chelpg.for"
      istart = aos(inew)
      jstart = aos(jnew)
      ial = 0
      iau = 5
      ibl = 0
      ibu = 5
      ima = 0
      imb = 0
      if (shellt(inew) .eq. 2) ima = 1
c
c MATFIL
c# 1733 "chelpg.for"
      if (shellt(jnew) .eq. 2) imb = 1
c# 1735 "chelpg.for"
  120 intc = 0
      do 170 j = 1, lbmax
      do 170 i = 1, lamax
      intc = intc + 1
      if (((la .gt. 1) .and. (i .gt. ial)) .and. (i .lt. iau)) goto 170
      if ((la .eq. 1) .and. (i .eq. ima)) goto 170
      if (((lb .gt. 1) .and. (j .gt. ibl)) .and. (j .lt. ibu)) goto 170
      if ((lb .eq. 1) .and. (j .eq. imb)) goto 170
      ind = (istart + i) - 1
      jnd = (jstart + j) - 1
      if (ind - jnd) 130, 140, 150
  130 ij = lind(jnd) + ind
      goto 160
  140 ij = lind(ind + 1)
      goto 160
  150 ij = lind(ind) + jnd
  160 a(ij) = aa(intc)
  170 continue
      return
      end
c
c
c        MATRIX MULTIPLICATION ROUTINE
c
c
c
      subroutine multay(a, y, x, n, maxdim)
c
c
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
c
c
      dimension a(maxdim, maxdim), y(maxdim), x(maxdim)
c
c
      data zero / 0.0 /
c# 1765 "chelpg.for"
      do 200 irow = 1, n
      sum = zero
      do 100 jcol = 1, n
      sum = sum + (a(irow,jcol) * y(jcol))
  100 continue
      x(irow) = sum
  200 continue
      return
c
c
c# 1773 "chelpg.for"
      end
c
c
c
c
c        L.E. CHIRLIAN
c
c        APRIL 1985
c
c
c
c        A SUBROUTINE TO OUTPUT THE CHARGES AND OTHER PERTINANT
c
c        INFORMATION FROM THE CHELP PROGRAM
c
c
c	 Slightly Modified for CHELPG operations by Curt Breneman
c	 Yale University Department of Chemistry, 2/88
c
c
      subroutine output()
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
      parameter (npoints = 50000)
      integer*4 shella, shelln, shellt, shellc, aos, aon, shladf
c
c
      character chkfil*40
      common /io/ in, iout, ipunch
c+++
      common /ipo/ ipo(5)
c+++
c      COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NE,NBASIS,IAN(101),
c
c     1             ATMCHG(100),C(3,100)
c
      common /mol/ natoms, icharg, multip, nae, nbe, ne, nbasis, ian(401
     &), atmchg(400), c(3, 400)
      common /out/ ntitle(20, 3), chkfil, q(400), nd, rms, pd, nlin,
     &nend(3)
      common /points/ p(3, npoints), np
c
c
c	write(6,*) 'Debug --', nlin,nend(1),nend(2),nend(3),nwords
c
c
c        CALCULATE THE DIPOLE MOMENT FROM THE FITTED CHARGES
c
c
c
c
c
      data deb / 0.393427328 /
c# 1811 "chelpg.for"
      dipx = 0.
      dipy = 0.
      dipz = 0.
      do 99 i = 1, natoms
      dipx = dipx + (q(i) * c(1,i))
      dipy = dipy + (q(i) * c(2,i))
      dipz = dipz + (q(i) * c(3,i))
   99 continue
      dipx = dipx / deb
      dipy = dipy / deb
c
c
c     CALCULATE TOTAL DIPOLE MOMENT
c
c
c
c# 1821 "chelpg.for"
      dipz = dipz / deb
c
c
c        CREATE OUTPUT  (continued from READIN)
c
c# 1825 "chelpg.for"
      diptot = sqrt(((dipx ** 2) + (dipy ** 2)) + (dipz ** 2))
c# 1829 "chelpg.for"
c
c  old header printed out here  -- moved to readin
c
c
c***
c***********************************************************************
c                WRITE GEOMETRY
c
c***********************************************************************
c
c
c# 1858 "chelpg.for"
      write(unit=iout, fmt=180)
  180 format(/2x,36x,18hMOLECULAR GEOMETRY)
      write(unit=iout, fmt=190)
  190 format(/,/,17x,13hATOMIC NUMBER,8x,1hX,12x,1hY,12x,1hZ)
      do 30 i = 1, natoms
      write(unit=iout, fmt=200) ian(i), c(1,i), c(2,i), c(3,i)
  200 format(/,23x,i2,8x,f10.7,3x,f10.7,3x,f10.7)
   30 continue
      write(unit=iout, fmt=210) icharg
  210 format(/2x,37hTHE TOTAL CHARGE IS CONSTRAINED TO:  ,i3)
      write(unit=iout, fmt=240)
  240 format(/,36x,11hNET CHARGES)
      write(unit=iout, fmt=250)
  250 format(/,28x,13hATOMIC NUMBER,5x,6hCHARGE)
      write(unit=iout, fmt=260) (ian(i), q(i),i = 1, natoms)
      write(unit=6, fmt=101) diptot
  101 format(/,2x,40hTHE DIPOLE MOMENT OF THESE CHARGES IS:  ,f8.5)
  260 format(/,34x,i2,10x,f8.4)
      write(unit=iout, fmt=270) np
  270 format(/,2x,34hFIT TO ELECTROSTATIC POTENTIAL AT ,i6,7h POINTS)
      write(unit=iout, fmt=280) rms
  280 format(/,2x,30hROOT MEAN SQUARE DEVIATION IS ,f6.4,10h KCAL/MOLE)
      return
      end
c
c
c       WRITTEN BY M.M. FRANCL FOR THE
c
c        PRINCETON CHEMISTRY DEPARTMENT VAX 11/780 UNDER VMS 3.4.
c
c        MODIFIED BY L.E. CHIRLIAN UNDER VMS 3.7.
c
c
c	 MODIFIED FOR GAUSSIAN 86 BY CURT BRENEMAN (VMS 4.5)
c	 Modified for G88/90 by Curt Breneman, 2/89
c	 YALE UNIVERSITY DEPARTMENT OF CHEMISTRY.
c
c        THIS VERSION IS COMPATIBLE WITH GAUSSIAN 82 FROM CARNEGIE-
c
c        MELLON UNIVERSITY AND IS DESIGNED FOR THE INPUT OF MO AND
c
c        BASIS INFORMATION FROM CHECKPOINT FILES.  THIS VERSION IS
c
c        TO BE USED FOR THE DETERMINATION OF ATOMIC CHARGES FROM
c
c        ELECTROSTATIC POTENTIALS DETERMINED BY FIRST ORDER HARTREE
c
c        FOCK PERURBATION THEORY.
c
c
c
c        OLD LIMITATIONS:   NO MORE THAN 256 BASIS FUNCTIONS
c
c                      NO MORE THAN 80 SHELLS
c
c
c        NEW LIMITATIONS:   NO MORE THAN 1280 BASIS FUNCTIONS
c	               NO MORE THAN 400 SHELLS
c
c
      subroutine readin()
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
      integer*4 shella, shelln, shellt, shellc, aos, aon, shladf, filnum
      character savchkfil*40, chkfil*40
      parameter (numpts = 20)
c
      dimension line(20)
c
      common /io/ in, iout, ipunch
      common /zmat/ IANZ(400),IZ(400,4),B1(400),ALPHA(400),BETA(400),
     &                LB1(400),LALPHA(400),LBETA(400),NZ,NVAR
c
c===  Gaussian88 Modification.  New Common /b/ size.
      common /mol/ natoms, icharg, multip, nae, nbe, ne, nbasis,
     &             ian(401), atmchg(400), c(3, 400)
c%%%g88/g89
c      common /b/ exx(6000), c1(6000), c2(6000), c3(2000), cf(2000),
c     &shladf(4000), x(2000), y(2000), z(2000), jan(2000), shella(2000),
c     &shelln(2000), shellt(2000), shellc(2000), aos(2000), aon(2000),
c     &nshell, maxtyp
c g88 cray
      common /b/ exx(6000), c1(6000), c2(6000), c3(2000), cf(2000),
     &shladf(2000), x(2000), y(2000), z(2000), jan(2000), shella(2000),
     &shelln(2000), shellt(2000), shellc(2000), aos(2000), aon(2000),
     &nshell, maxtyp
c
      common /ipo/ ipo(5)
      common /out/ ntitle(20, 3), chkfil, q(400), nd, rms, pd, nlin,
     &nend(3)
      common /charge/ coef_alpha(100000), coef_beta(100000), iuhf
      common /sphere/ vdwr(400), npts
c
	DATA MAXZ/400/
	DATA BOHRTANG/0.52917706/
      data maxatm / 400 /
      data maxbas / 1280 /
      data maxpts / 50000 /
      data npo / 5 /
      data iunit / 3 /
      data iread / 2 /
      data ifind / 11 /
      data iblnk / '    ' /
c
      in = 5
      iout = 6
C print out header info:
c
      write(iout,'(/,/,12x,51(1h*))')
        itime=0
        idate=0
        call forjdate(idate,itime)
c PRINT  identifier and DATE
        write(unit=iout, fmt=170) idate,itime
 170	format(/17x,28hC  H  E  L  P  G  for  G 9 0,6x,a8,4x,a8)
c
      write(unit=iout, fmt=100)
  100 format(/,17x,41hCHARGES FROM ELECTROSTATIC POTENTIAL GRID)
      write(unit=iout, fmt=110)
  110 format(/,36x,9hCHELPGrid,/)
      write(unit=iout, fmt=111)
  111 format(/,15x,
     &50hGrid Modification - C. Breneman,  Yale Chem, 1989,/)
      WRITE(IOUT,'(1h ,30x,a)')' vsn 12/90'
      write(iout,'(/,/,12x,51(1h*))')
c
c       GET CHECKPOINT FILE NAME
c
	read(in,'(a40)') chkfil
c  	savchkfil=chkfil
c 1000 format(a40)
      write(unit=iout, fmt=150) chkfil
  150 format(/2x,18hCHECKPOINT FILE:  ,a40///)
c
c     INPUT INFORMATION
c
 1007 format(1h ,a40)
 1008 format(1h ,20a4)
 1010 format(20a4)
      do 1921 i = 1, 3
        nlin = i
        read(unit=in, fmt=1010) line
c        write(unit=iout, fmt=1008)line
        if (line(1) .eq. iblnk) then
          nlin = nlin - 1
          goto 192
        end if
        do 1922 j = 1, 20
          l = 20 - j
          if (line(l) .ne. iblnk) then
            nend(i) = l
            do 1923 k = 1, nend(i)
              ntitle(k,i) = line(k)
 1923       continue
            goto 1921
          end if
 1922   continue
 1921 continue
      nlin = nlin + 1
      do 24 i = 1, nlin
        write(unit=iout, fmt=1200) (ntitle(j,i),j = 1, nend(i))
 1200   format(2x,19a4)
   24 continue
c
c         READ IN PRINT OPTIONS
c
  192 continue
c
c        READ IN # OF D FUNCTIONS
c        NOTE: IF THE BASIS SET USES 5 D FUNCTION, ND MUST BE
c        SET EQUAL TO 1 TO ACCOMODATE THE INTEGRAL PACKAGE.
c        IF THE BASIS SET USES  6 D FUNCTIONS, ND IS SET EQUAL TO 0.
c
c# 2031 "chelpg.for"
3000    read(unit=in, fmt=*) (ipo(i),i = 1, npo)
c kjb start
c      write(unit=iout,fmt=*)(ipo(i),i=1,npo)
c kjb stop
      read(unit=in, fmt=*) nd
c kjb start
c      write(unit=iout,fmt=*)nd
c kjb stop
      if ((nd .ne. 5) .and. (nd .ne. 6)) then
        write(iout,*) '# OF D FUNCTIONS MUST BE 5 OR 6'
	write(iout,*) '  assuming 6 d funtions ...'
	nd=6
      end if
      if (nd .eq. 5) then
      nd = 1
      goto 15
      end if
      nd = 0
c
c        INITIATE FILEIO
c
c***   Different Fopen statement for Trace-7
c# 2048 "chelpg.for"
   15 continue
c
c# 2054 "chelpg.for"
c	chkfil=savchkfil
	write(iout,'(9h opening ,a,3h.../)') chkfil
c	call fclose(iunit,1)
C      call fopen(iunit, 5, chkfil, ialloc, JUNK)
      call fopen(iunit, 5, chkfil, ialloc)
c***********************************************************************
c
c        READ IN COMMON /MOL/
c
c      NWORDS = 1804
c      IFILENO = 30997
c
c      CALL FILEIO (IREAD,IFILENO,NWORDS,NATOMS,IALLOC)
c
c# 2057 "chelpg.for"
      iwwrit = ipo(1)
c# 2065 "chelpg.for"
      irwmol = 997
      lenmol = (4 * maxatm) + intowp(8 + maxatm)
c
c***********************************************************************
c
c        READ IN SPHERE DATA (VAN DER WAALS RADII, # OF POINTS TO FIT
c
c kjb start
c      write(unit=iout,fmt=*)' npts=',npts
c      write(unit=iout,fmt=*)' natoms,nbasis,c(1,3),c(3,3)=',
c     &    natoms,nbasis,c(1,3),c(3,3)
c kjb stop
c# 2068 "chelpg.for"
      ifileno= filnum(irwmol,iunit)
C      write(unit=iout,fmt=*)'#2068 ifileno=',ifileno
      call fileio(2,-ifileno, lenmol, natoms, 0)
c kjb start
c      write(unit=iout,fmt=*)'after:natoms,nbasis,lenmol=',
c     &      natoms,nbasis,lenmol
c      write(unit=iout,fmt=*)'c(xyz,i)',
c     &     (c(1,i),c(2,i),c(3,i),i=1,natoms)
c	write(unit=iout,fmt=*) 'atmchg',(atmchg(i),i=1,natom)
c kjb stop
c# 2076 "chelpg.for"
      read(unit=in, fmt=*) (vdwr(i),i = 1, natoms)
c kjb start
c      write(unit=iout,fmt=*)(vdwr(i),i=1,natoms)
c kjb stop
      read(unit=in, fmt=*) npts
c kjb start
c      write(unit=iout,fmt=*)' npts=',npts
c kjb stop
      if (npts .gt. maxpts) then
      stop 'MAXIMUM NUMBER OF POINTS MUST BE LESS THAN 50000'
      end if
c
c        SET MAXIMUM VAN DER WAALS RADII TO 4
c
c# 2084 "chelpg.for"
      vmax = 4.
      do 20 i = 1, natoms
      if (vdwr(i) .gt. vmax) then
      write(unit=iout, fmt=2500) i
 2500 format(3x,31hTHE VAN DER WAALS RADII OF ATOM,i3,15hIS OUT OF RANGE
     &)
      stop
      end if
   20 continue
c# 2094 "chelpg.for"
      read(unit=in, fmt=*) vfact
c kjb start
c      write(unit=iout,fmt=*)'vfact=',vfact
c kjb stop
      do 21 i = 1, natoms
      vdwr(i) = vdwr(i) * vfact
   21 continue
c***********************************************************************
c
c        READ IN BASIS SET INFORMATION (COMMON /B/)
c
c# 2102 "chelpg.for"
ckjb start
c      write(unit=iout,fmt=*)' #2102 nwords=',nwords
ckjb stop
      ifileno = 30506
      call fileio(ifind, ifileno, nwords, exx, 0)
c***********************************************************************
c
ckjb start
c      write(unit=iout,fmt=*)' #2104 nwords=',nwords
ckjb stop
      call fileio(iread, - ifileno, nwords, exx, 0)
c# 2106 "chelpg.for"
c      write(unit=iout, fmt=8000) (c(1,i), c(2,i), c(3,i),i = 1,400)
      write(unit=iout, fmt=8000) (c(1,i), c(2,i), c(3,i),i = 1, natoms)
 8000 format(/1x,11hCOORDINATES/(1x,3f12.6))
C	write(unit=iout,fmt='(/1x,3hxyz/(1x,3f12.6))')
C     &        (x(i),y(i),z(i),i=1,natoms)
c
      write(unit=iout, fmt=8020) natoms, icharg, multip, nae, nbe, ne,
     &nbasis, nshell, maxtyp
 8020 format(/1x,11hNATOMS   = ,i3/1x,9hICHARG = ,i3/1x,9hMULTIP =
     &,i3/1x,9hNAE    = ,i3/1x,9hNBE    = ,i3/1x,9hNE     = ,i3/1x,
     &9hNBASIS = ,i3/1x,9hNSHELL = ,i3/1x,9hMAXTYP = ,i3)
c
      write(unit=iout, fmt=8030) (ian(i),i = 1, natoms)
 8030 format(/1x,3hIAN/(1x,20i3))
c
      write(unit=iout, fmt=8050) (jan(i), shellt(i), shella(i),i = 1,
     &nshell)
 8050 format(/1x,18hCENTER TYPE SHELLA/(1x,3i7))
c
      write(unit=iout, fmt=8055) shella(nshell + 1)
 8055 format(1x,14x,i7)
c
      write(unit=iout, fmt=8060) (exx(i), c1(i), c2(i),i = 1, nshell)
 8060 format(/1x,12x,5hEXPON,8x,9hEXPCOF(S),8x,9hEXPCOF(P),/(1x,3e17.9))
c
      write(unit=iout, fmt=8070) (c3(i), cf(i),i = 1, nshell)
 8070 format(/1x,12x,9hEXPCOF(D),8x,9hEXPCOF(F),/(1x,2e17.9))
c
      write(unit=iout, fmt=3575)
 3575 format(/,15x,6hATOM #,5x,35hV.D.W. RADII (MULTIPLIED BY FACTOR))
c
        do 3500 i = 1, natoms
           write(unit=iout, fmt=4000) i, vdwr(i)
 4000      format(15x,i5,20x,f5.2)
c          WRITE (IOUT,4500)NPTS
 4500     format(/2x,23hNUMBER OF POINTS TO FIT,i6)
 3500     continue
c
c***********************************************************************
c
c        READ IN ALPHA MO COEFFICIENTS
c
c# 2145 "chelpg.for"
      ifileno = 30524
      call fileio(ifind, - ifileno, nwords, coef_alpha, 0)
c***********************************************************************
c
c        READ IN THE BETA MO COEFFICIENTS
c
c# 2147 "chelpg.for"
      call fileio(iread, - ifileno, nwords, coef_alpha, 0)
c# 2152 "chelpg.for"
      ifileno = 30526
      call fileio(ifind, ifileno, nwords, coef_beta, 0)
      if (nwords .eq. 0) then
      iuhf = 0
      goto 300
      end if
      iuhf = 1
c***********************************************************************
c# 2159 "chelpg.for"
      call fileio(iread, - ifileno, nwords, coef_beta, 0)
c***********************************************************************
c
c# 2161 "chelpg.for"
  300 continue
c# 2163 "chelpg.for"
	IFILENO=30507
C	NWORDS=3*MAXNZ +INTOWP(8*MAXNZ+2)
	CALL FILEIO(IFIND, - IFILENO, NWORDS, IANZ,0)
	CALL FILEIO(IREAD, - IFILENO, NWORDS, IANZ,0)
	write(iout,'(1h )')
	write(iout,'(15H Z MATRIX: for ,A/)')CHKFIL
        CALL ZPRINT(400,NZ,IANZ,IZ,B1,ALPHA,BETA,BOHRTANG)
	call fclose(iunit,1)
      return
      end
c
c
c        MODIFIED FOR POLARIZATION POTENTIAL CALCULATIONS
c
c        M.M. FRANCL  FEBRUARY 1984
c
c
c
      subroutine reduc1(x, lamax, lbmax, i6to5)
c
c
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
c
c
      dimension x(100), s(100), ind5(9), ind6(10)
c
c
c     ******************************************************************
c
c     ROUTINE REORDERS FROM ARRANGEMENT: S,Z,ZZ,X,XZ,XX,Y,YZ,XY,YY
c
c     TO                                 S,X,Y,Z,XX,YY,ZZ,XY,XZ,YZ
c
c     OR FROM   S,Z,X,Y TO S,X,Y,Z
c
c     THIS ENSURES LABELING COMPATIBILITY BETWEEN THE SP AND D PACKAGES
c
c     AT THE SAME TIME THE INTEGRALS ARE MOVED TO THE FIRST 1,4,10,16,
c
c     40 OR 100 LOCATIONS, DEPENDING ON THE SHELL QUANTUM NUMBERS
c
c     ******************************************************************
c
      data pt5 / 0.5 /
      data r3ov2 / 0.8660254040 /
      data ind5 / 1, 4, 7, 2, 3, 6, 9, 5, 8 /
      data ind6 / 1, 4, 7, 2, 6, 10, 3, 9, 5, 8 /
c# 2191 "chelpg.for"
      nword = lamax * lbmax
      if (nword - 1) 5, 180, 5
    5 intc = 0
      if (i6to5 .eq. 1) goto 40
   10 do 20 i = 1, lbmax
      isb = 10 * (ind6(i) - 1)
      do 20 j = 1, lamax
      isa = isb + ind6(j)
      intc = intc + 1
   20 s(intc) = x(isa)
c     ******************************************************************
c
c     ROUTINE TO REDUCE SIX D FUNCTIONS TO FIVE
c
c     ALSO REORDERS FROM : S,Z,ZZ,X,XZ,XX,Y,YZ,YX,YY
c
c                     TO   S,X,Y,Z,ZZ,XX-YY,XY,XZ,YZ
c
c     OR FROM S,Z,X,Y TO S,X,Y,Z
c
c     FOR COMPATIBILITY WITH SP PACKAGE
c
c     ******************************************************************
c
c# 2201 "chelpg.for"
      goto 160
c# 2209 "chelpg.for"
   40 do 150 i = 1, lbmax
c     IFB=0 FOR S,X,Y,Z,XY,XZ,YZ, IFB=1 FOR ZZ-RR, IFB=2 FOR XX-YY
c
c# 2210 "chelpg.for"
      isb = 10 * (ind5(i) - 1)
c# 2212 "chelpg.for"
      ifb = 0
      if (i .eq. 5) ifb = 1
      if (i .eq. 6) ifb = 2
   80 do 150 j = 1, lamax
      isa = isb + ind5(j)
      ifa = 0
      if (j .eq. 5) ifa = 1
      if (j .eq. 6) ifa = 2
  120 ihop = ((3 * ifb) + ifa) + 1
c
c
c     ******************************************************************
c
c     *                           (F,O,ZA2)                            *
c
c     ******************************************************************
c
c# 2221 "chelpg.for"
      goto (130, 122, 123, 124, 125, 126, 127, 128, 129), ihop
c# 2226 "chelpg.for"
  122 xx = zz1(x,isa,3,7)
c
c
c     ******************************************************************
c
c     *                           (F,O,XA2-YA2)                        *
c
c     ******************************************************************
c
c# 2227 "chelpg.for"
      goto 140
c# 2232 "chelpg.for"
  123 xx = xy1(x,isa,4)
c
c
c     ******************************************************************
c
c     *                           (ZB2,O,F)                            *
c
c     ******************************************************************
c
c# 2233 "chelpg.for"
      goto 140
c# 2238 "chelpg.for"
  124 xx = zz1(x,isa,30,70)
c
c
c     ******************************************************************
c
c     *                           (ZB2,O,ZA2)                          *
c
c     ******************************************************************
c
c# 2239 "chelpg.for"
      goto 140
c# 2244 "chelpg.for"
  125 xx = zz1(x,isa,30,70) - (pt5 * (zz1(x,isa + 3,30,70) + zz1(x,isa
     & + 7,30,70)))
c
c
c     ******************************************************************
c
c     *                           (ZB2,O,XA2-YA2)                      *
c
c     ******************************************************************
c
c# 2245 "chelpg.for"
      goto 140
c# 2250 "chelpg.for"
  126 xx = r3ov2 * (zz1(x,isa,30,70) - zz1(x,isa + 4,30,70))
c
c
c     ******************************************************************
c
c     *                           (XB2-YB2,O,F)                        *
c
c     ******************************************************************
c
c# 2251 "chelpg.for"
      goto 140
c# 2256 "chelpg.for"
  127 xx = xy1(x,isa,40)
c
c
c     ******************************************************************
c
c     *                           (XB2-YB2,O,ZA2)                      *
c
c     ******************************************************************
c
c# 2257 "chelpg.for"
      goto 140
c# 2262 "chelpg.for"
  128 xx = r3ov2 * (zz1(x,isa,3,7) - zz1(x,isa + 40,3,7))
c
c
c     ******************************************************************
c
c     *                           (XB2-YB2,O,XA2-YA2)                  *
c
c     ******************************************************************
c
c# 2263 "chelpg.for"
      goto 140
c# 2268 "chelpg.for"
  129 xx = r3ov2 * (xy1(x,isa,4) - xy1(x,isa + 40,4))
c
c
c     ******************************************************************
c
c     *                           (F,O,F)                              *
c
c     ******************************************************************
c
c# 2269 "chelpg.for"
      goto 140
c# 2274 "chelpg.for"
  130 xx = x(isa)
  140 intc = intc + 1
  150 s(intc) = xx
  160 do 170 i = 1, nword
  170 x(i) = s(i)
  180 return
      end
c
c
c        ROUTINE TO MODIFY COMMON /B/ TO THE EXPECTED FORMAT FOR INTGRL
c
c        FOR BASIS SETS HAVING P ONLY SHELLS, SUCH AS THE 6-31G** BASIS
c
c
c
      subroutine star(nbasis, shellt, shellc, aos, nshell, nostar)
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
c
c
      integer*4 shellc, shellt, aos
c
c
c        LOOP OVER SHELLS
c
c
c
      dimension shellc(2000), shellt(2000), aos(2000)
c# 2293 "chelpg.for"
      do 100 i = 1, nshell
      if ((shellt(i) .eq. 1) .and. (shellc(i) .eq. 1)) then
      nbasis = nbasis + 1
      nostar = 1
      do 200 j = i, nshell
      aos(j) = aos(j) + 1
  200 continue
      end if
  100 continue
      return
      end
c
c
c        ROUTINE TO CALCULATE THE ELECTROSTATIC POTENTIAL FROM FIRST ORD
cER
c        PERTURBATION THEORY
c
c
c
c        M.M. FRANCL    JULY 1985
c
c        MODIFIED VERSION OF A MEPHISTO ROUTINE
c
c        RESTRICTED TO UNRESTRICTED HARTREE-FOCK WAVEFUNCTIONS
c
c
c
      subroutine uep()
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
      parameter (npoints = 50000)
c
c
      integer*4 shella, shelln, shellt, aos, shellc, aon, handle
      common /io/ in, iout, ipunch
c+++
      common /ipo/ ipo(5)
c
c===	Gaussian88 Modification for enlarged common /b/.
      common /mol/ natoms, icharg, multip, nae, nbe, nel, nbasis, ian(
     &401), atmchg(400), c(3, 400)
c
c===	Old G86 Common /b/
c      COMMON/B/EXX(1200),C1(1200),C2(1200),C3(1200),
c     $         X(400),Y(400),Z(400),JAN(400),SHELLA(400),SHELLN(400),
c     $         SHELLT(400),SHELLC(400),AOS(400),AON(400),NSHELL,MAXTYP
c+++
c      COMMON /B/ EXX(240),C1(240),C2(240),C3(240),X(80),Y(80),Z(80),
c
c     $           JAN(80),SHELLA(80),SHELLN(80),SHELLT(80),SHELLC(80)
c
c     $          ,AOS(80),AON(80),NSHELL,MAXTYP
c
c      COMMON /MOL/ NATOMS,ICHARG,MULTIP,NAE,NBE,NEL,NBASIS,IAN(101),
c
c     $             ATMCHG(100),C(3,100)
c
      common /B/ exx(6000), c1(6000), c2(6000), c3(6000), x(2000), y(
c      common /BB/ exx(6000), c1(6000), c2(6000), c3(6000), x(2000), y(
     &2000), z(2000), jan(2000), shella(2000), shelln(2000), shellt(2000
     &), shellc(2000), aos(2000), aon(2000), nshell, maxtyp
      common /points/ p(3, npoints), maxpnts
      common /elp/ elecp(npoints)
      common /charge/ coef_alpha(100000), coef_beta(100000), iuhf
c
c
      common /out/ ntitle(20, 3), chkfil, q(400), i6to5, rms, percent,
     &nlin, nend(3)
c
c
      dimension hpert(100000), index(1280)
c
c
c        INTIALIZE TIMING
c
c
c
c
c
c        SET UP THE INDEXING TABLE FOR HPERT
c
c
c
      data iptchg / 1.0 /
      data zero / 0.0 /
      data two / 2.0 /
      data vnucmax / 30.0 /
c# 2351 "chelpg.for"
      handle = 0
c# 2355 "chelpg.for"
      do 100 i = 1, nbasis
      index(i) = ((i - 1) * i) / 2
c
c
c        BEGIN LOOP TO CALCULATE ELECTROSTATIC POTENTIAL
c
c
c
c# 2357 "chelpg.for"
  100 continue
c# 2361 "chelpg.for"
      do 200 npnt = 1, maxpnts
      x1 = p(1,npnt)
      x2 = p(2,npnt)
c
c
c     CALCULATE THE ONE-ELECTRON INTEGRALS
c
c
c
c# 2364 "chelpg.for"
      x3 = p(3,npnt)
c# 2368 "chelpg.for"
      if (ipo(5) .eq. 1) then
      write(unit=iout, fmt=3010)
c***
c      ISTAT = LIB$INIT_TIMER(HANDLE)
c
c***
c# 2370 "chelpg.for"
 3010 format(1x,18hTIME FOR INTEGRALS)
c
c
      end if
c
c
c***
c      IF (IPO(5).EQ.1) ISTAT = LIB$SHOW_TIMER(HANDLE)
c
c***
c
c
c# 2376 "chelpg.for"
      call intgrl(hpert, x1, x2, x3, iptchg, i6to5)
c
c
c# 2382 "chelpg.for"
      if (ipo(4) .eq. 1) call linout(hpert, nbasis, 0, 0)
c# 2384 "chelpg.for"
      if (ipo(5) .eq. 1) then
      write(unit=iout, fmt=3000)
c***
c      ISTAT = LIB$INIT_TIMER(HANDLE)
c
c***
c# 2386 "chelpg.for"
 3000 format(1x,18hTIME FOR TRANSFORM)
c
c
c     FORM THE HPERT MATRIX ELEMENTS
c
c
c
c        ALPHA CODE
c
c
c
c# 2390 "chelpg.for"
      end if
c# 2396 "chelpg.for"
      e = zero
c
c
c        SUM OVER OCCUPIED ALPHA MOS
c
c
c
c# 2397 "chelpg.for"
      icoefi = - nbasis
c# 2401 "chelpg.for"
      do 220 ii = 1, nae
c
c
c        CALCULATE ELECTROSTATIC POTENTIAL
c
c
c
c# 2402 "chelpg.for"
      icoefi = icoefi + nbasis
c# 2406 "chelpg.for"
      do 221 ip = 1, nbasis
      cpi = coef_alpha(icoefi + ip)
c
c
c# 2408 "chelpg.for"
      ipdex = index(ip)
c# 2410 "chelpg.for"
      do 222 iq = 1, ip
      e = e + ((cpi * coef_alpha(icoefi + iq)) * hpert(ipdex + iq))
  222 continue
      do 223 iq = ip + 1, nbasis
      e = e + ((cpi * coef_alpha(icoefi + iq)) * hpert(ip + index(iq)))
c
c
c# 2415 "chelpg.for"
  223 continue
c# 2417 "chelpg.for"
  221 continue
c
c
c        BETA CODE
c
c
c
c# 2418 "chelpg.for"
  220 continue
c
c
c        SUM OVER OCCUPIED BETA MOS
c
c
c
c# 2422 "chelpg.for"
      icoefi = - nbasis
c# 2426 "chelpg.for"
      do 420 ii = 1, nbe
c
c
c        CALCULATE ELECTROSTATIC POTENTIAL
c
c
c
c# 2427 "chelpg.for"
      icoefi = icoefi + nbasis
c# 2431 "chelpg.for"
      do 421 ip = 1, nbasis
      cpi = coef_beta(icoefi + ip)
c
c
c# 2433 "chelpg.for"
      ipdex = index(ip)
c# 2435 "chelpg.for"
      do 422 iq = 1, ip
      e = e + ((cpi * coef_beta(icoefi + iq)) * hpert(ipdex + iq))
  422 continue
      do 423 iq = ip + 1, nbasis
      e = e + ((cpi * coef_beta(icoefi + iq)) * hpert(ip + index(iq)))
c the follwing stmt was found active, it has been deactivated:  kjb
c      e = e + ((cpi * coef_) * hpert(ip + index(iq)))
c
c
c# 2441 "chelpg.for"
  423 continue
c# 2443 "chelpg.for"
  421 continue
c
c
c***
c      IF (IPO(5) .EQ. 1) ISTAT = LIB$SHOW_TIMER(HANDLE)
c
c***
c
c
c        CALCULATE NUCLEAR PART OF ELECTROSTATIC POTENTIAL
c
c
c
c# 2444 "chelpg.for"
  420 continue
c# 2452 "chelpg.for"
      vnuc = zero
      do 300 iatom = 1, natoms
      del1 = c(1,iatom) - x1
      del2 = c(2,iatom) - x2
      del3 = c(3,iatom) - x3
      ra = sqrt(((del1 * del1) + (del2 * del2)) + (del3 * del3))
      if (ra .eq. zero) then
      vnuc = vnucmax
      goto 310
      end if
      vnuc = vnuc + (ian(iatom) / ra)
  300 continue
c
c
c# 2464 "chelpg.for"
  310 continue
c# 2466 "chelpg.for"
      elecp(npnt) = e + (vnuc * iptchg)
      if (ipo(5) .eq. 1) write(unit=iout, fmt=*) 'E(', npnt, ') = ', e
  200 continue
      return
      end
c
c
c        ROUTINE TO REFORMAT COMMON/B/ AND THE H ARRAY FOR BASIS
c
c        SETS HAVING P ONLY SHELLS
c
c
c
      subroutine unstar(nbasis, shellt, shellc, aos, nshell, h, nostar)
c      implicit double precision (o-z, a-h)
	implicit real*8 (a-h,o-z)
c
c
      integer*4 shellc, shellt, aos
c
c
      dimension shellc(2000), shellt(2000), aos(2000), h(1)
c
c
c        LOOP OVER SHELLS
c
c
c
c# 2481 "chelpg.for"
      if (nostar .eq. 0) return
c# 2485 "chelpg.for"
      do 100 i = 1, nshell
c
c
c        REMOVE EXTRA ROWS AND COLUMNS
c
c
c
c        LOOP OVER ROWS
c
c
c
c# 2486 "chelpg.for"
      if ((shellt(i) .eq. 1) .and. (shellc(i) .eq. 1)) then
c# 2492 "chelpg.for"
      ibasis = aos(i)
      item = ((ibasis - 1) * ibasis) / 2
c
c
c        LOOP OVER COLUMNS
c
c
c
c# 2494 "chelpg.for"
      do 200 j = ibasis + 1, nbasis
c# 2498 "chelpg.for"
      newitem = ((j - 1) * j) / 2
      do 250 k = 1, ibasis - 1
      item = item + 1
      newitem = newitem + 1
      h(item) = h(newitem)
c
c
c        SKIP THE VALUE IN THE IBASIS TH COLUMN
c
c
c
c# 2503 "chelpg.for"
  250 continue
c
c
      newitem = newitem + 1
c# 2509 "chelpg.for"
      do 260 k = ibasis + 1, j
      item = item + 1
      newitem = newitem + 1
      h(item) = h(newitem)
  260 continue
  200 continue
c
c
c        RESTRUCTURE AOS TO ACCOUNT FOR THE S SHELL REMOVED
c
c
c
c# 2515 "chelpg.for"
      nbasis = nbasis - 1
c# 2519 "chelpg.for"
      do 300 iaos = i, nshell
      aos(iaos) = aos(iaos) - 1
c
c
c# 2521 "chelpg.for"
  300 continue
c# 2523 "chelpg.for"
      end if
  100 continue
      return
      end
c
c
      function xy1(x, i, iy)
c
c
      dimension x(100)
c
c
      data halfr3 / 0.8660254040 /
c# 2533 "chelpg.for"
      xy1 = halfr3 * (x(i) - x(i + iy))
      return
      end
c
c
      function zz1(x, i, ix, iy)
c
c
      dimension x(100)
c
c
      data half / 0.5 /
c# 2542 "chelpg.for"
      zz1 = x(i) - (half * (x(i + ix) + x(i + iy)))
      return
      end
c
c
        subroutine forjdate(idate,itime)
c	use system dependent calls to computer clock:
C IDATE  MM/DD/YY ITIME  HH-MM-SS  BOTH ARE BOOLEAN (USE A8)
        idate=date()
        itime=clock()
        return
        end