|
|
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
|
