Date: Fri, 12 May 1995 15:40:12 -0400
From: mfrancl (Francl Michelle M)
To: srusso
Subject: chelp.f
X-UIDL: 800307735.006

C
C	CHELP
C
C	(Net Atomic) Charges fit to ELectrostatic Potentials
C
C	M.M. Francl
C	L.E. Chirlian
C
c	platform:
c	runs on SGI under Iris 4.0.5
c
c	ab initio interface:
c	currently uses Spartan archive files as input for
c	wavefunction
c
c	code history:
c
C	April 1985
C	Princeton Chemistry Department VAX 11/780
C	VMS 3.7
c
C	modified 6/91 mmf Bryn Mawr College
c	to do point selection via random number algorithm
cc
c	modified 2/94 mmf Bryn Mawr College
c	to eliminate dipole constraint
c
c
c	modified 6/94 mmf Bryn Mawr College
c	to fit to selected subset of charges
c
implicit real*8 (a-h,o-z)
c
c	read in data from checkpoint file
c
call readin
c
c	choose random points within a user-selected box until a pre-set
c	maximum has been reached
c
call select
c
c	calculate the electrostatic potential using first order Hartree-Fock
c	perturbation theory
c
call ep
c
c	using method of Lagrange multipliers, fit by least squares the charges
c	to the electrostatic potential, constraining the fit to reproduce the
c	total molecular charge and the molecular dipole moment
c
call fit
c
c	print out table of results
c
call out
c
end
c
SUBROUTINE READIN
C
C	written by M.M. Francl for the
C	Princeton Chemistry Department VAX 11/780 under VMS 3.4.
c	modified by L.E. Chirlian under VMS 3.7.
C	This version is compatible with GAUSSIAN 82 from Carnegie-
C	Mellon University and is designed for the input of MO and
C	basis information from checkpoint files. This version is
C	to be used for the determination of atomic charges from
C	electrostatic potentials determined by first order Hartree
c	Fock perurbarion theory.
c
c	modified 5/29/91 to read in size of VDW envelope (FACTOR)
c	and step size between shells (DELR)
c	mmf
cc	modified 6/3/91 to deal with random # point selection
c	mmf
c
c	modified 1993 to use Spartan archive file C.Carey
C
C	limitations: no more than 256 basis functions
C	no more than 80 shells
C
implicit real*8 (a-h,o-z)
character*8 title1, title2, title3
integer*4 shella,shelln,shellt,shellc,aos,aon,shlAdf character*40 arcfil
character*80 temp
character*8 word,beta

C
COMMON /IO/ IN,IOUT
common /mol/ natoms,icharge,multip,nae,nbe,ne,nbasis,ian(101), $	atmchg(100),c(3,100)
common /b/ exx(1000),c1(1000),c2(1000),c3(1000),cf(1000), $	shladf(1000),x(1000),y(1000),z(1000),
$	jan(1000),shella(1000),shelln(1000),shellt(1000),
$	shellc(1000),
$	aos(1000),aon(1000),nshells,maxtyp
common /ipo/ ipo(10)
common /xout/ title1(10),title2(10),title3(10),arcfil, 1	dipx, dipy, dipz, q(104), nd,rms,percent
common /charge/ coef_alpha(9000),coef_beta(9000),iUHF common /sphere/ vdwr(104), npts, nvect,vect(3,3074) common /points/ p(3,50000), maxpnts
common /envelope/ factor,factor2,xlen,ylen,zlen common /pot/ eigen(300),eigenII(300)
common /tobefit/ isfit(100),nfit
C
DATA MAXbas /300/
DATA NPO/ 10/
data iunit/ 3/, iread/ 2/, ifind/11/
in = 5
iout = 6
C
C	checkpoint file name
C
READ(IN,1000) ARCFIL
1000 FORMAT(a40)
C
C	INPUT INFORMATION
C
read (in,1010) (title1(i),i=1,10)
1010 format(9a8,a7)
read (in,1010)(title2(i),i=1,10)
if (title2(1) .eq. ' ') then
goto 24
end if
read (in,1010) (title3(i),i=1,10)
if (title3(1) .eq. ' ') then
goto 24
end if
24 continue
c
C
C	read in print options
C
3000 READ(IN,*) (IPO(I),I=1,NPO)
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
c	0.
c
read(in,*) nd
if (nd .ne. 5 .and. nd .ne. 6) then
stop '# of d functions must be 5 or 6'
end if
if (nd .eq. 5) then
nd = 1
goto 15
end if
nd = 0
15 continue
c
c	read in dipole moment - TURNED OFF NOW
c	the three components are read in seperately
c
cdip	read (in,*)dipx,dipy,dipz
c
C*********************************************************************** C
C open archive
C
open(iunit, file=arcfil,status='OLD')
open = FALSE
C
IWWRIT = IPO(1)
C
C*********************************************************************** C
C open SPARTAN Archive
C
C read banner and title - discard
C
read(iunit,1030) temp
read(iunit,1030) temp
C
C*********************************************************************** C
C read #atoms, #basis fns, #shells, #gaussians, molecular charge, and C spin multiplicity.
C (NORB is the #of orbitals stored on archive. for open shells: total C alpha = total beta)
C
read(iunit,9002) natoms,nbasis,nshells,ngauss,icharge, $	multip,norb,irs
C
C*********************************************************************** C
C more than 300 basis fns
C
if (nbasis.gt.300) stop 'exceeds masbas' C
C read filler line
C
read(iunit,1020) word
C
C read atomic numbers and cartesian coords - these are in bohrs C
do 10 i=1,natoms
read(iunit,*) ian(i),(c(j,i),j=1,3)
10 continue
C
C read filler line
C
read(iunit,1020) word
C
C*********************************************************************** C
C read basis set info
C
icnt=1
do 29 i=1,nshells
read(iunit,*) shellt(i),shelln(i),shella(i),jan(i), $	shellc(i)
x(i)=c(1,jan(i))
y(i)=c(2,jan(i))
z(i)=c(3,jan(i))
if (shellt(i).lt.2) then
shladf(i)=0
else
shladf(i)=shella(i)
endif
if (i.eq.1) then
aos(i)=icnt
else
if ( ((ian(jan(i)) .le. 2) .and. (shellt(i) .eq. 0)) $ .or. ((ian(jan(i)) .gt. 2) .and. (ian(jan(i)) .le. 18) $ .and. (shellt(i) .le. 1)) .or. ((ian(jan(i)) .gt. 18) $ .and. (ian(jan(i)) .le. 54) .and. (shellt(i) .le. 2)) ) then 
aos(i)=icnt
else
aos(i)=aos(i-1)
endif
endif
if (shellt(i).eq.0) icnt=icnt + 1
if (shellt(i).eq.1) icnt=icnt + 4
if (shellt(i).eq.2) icnt=icnt + 6
if (shellt(i).eq.3) icnt=icnt + 10
29	continue
C
C*********************************************************************** C exx- gaussian exponent
C c1 - a coefficient
C c2 - p coefficient
C c3 - d coefficient
C cf - f coefficient
C
do 30 i=1,ngauss
read(iunit,*) exx(i)
read(iunit,*) c1(i),c2(i),c3(i),cf(i)
30 continue
C
C*********************************************************************** C read filler line
read(iunit,1020) word
C
C read hf and mp2 energies - discard
read(iunit,*)ehf,emp2
C
C read filler line
read(iunit,1020) word
C
C*********************************************************************** C
C read in molecular energies (to discard) and the alpha coefficients C
read(iunit,*) (eigen(i),i=1,nbasis)
C
do 40 imo=1,norb
read(iunit,*) (coef_alpha(nbasis*(imo-1)+j),j=1,nbasis) 40 continue
iUHF=0
c
C*********************************************************************** C
C read in molecular energies (to discard ) and the beta MO coefficients C
read(iunit,1020,END=80) word
C
if (word.eq.'BETA ') then
read(iunit,*) (eigenII(i),i=1,nbasis)
iUHF=1
C
do 41 imo=1,norb
read(iunit,*) (coef_beta(nbasis*(imo-1)+j),j=1,nbasis) 41 continue
open=TRUE
endif
80 continue
C
close(iunit,status='KEEP')
C
C*********************************************************************** C
C calculations
C
ne=0
do 99 i=1,natoms
ne=ne + ian(i)
99 continue
ne=ne - icharge
C
itmp=ne - multip + 1
nae=itmp/2 + mulitp - 1
nbe=itmp/2
C
maxtyp=0
do 999 i=1,nshells
if (shellt(i).gt.maxtyp) maxtyp=shellt(i) 999 continue
C
C*********************************************************************** C
C read in sphere data (van der Waals radii, # of points to fit) C
read (in,*) (vdwr(i),i=1,natoms)
read (in,*) maxpnts
if (maxpnts .gt. 49999) then
stop 'Maximum number of points must be less than 49,999' end if
c
c	set maximum van der waals radii to 4
c
vmax=4.
do 20 i=1,natoms
if (vdwr(i) .gt. vmax) then
write (iout, 2500) i
2500	format (3x, 'The van der waals radii of atom', i3,
1	'is out of range')
stop
end if
20 continue
c
C************************************************************************ C
C read in the size of VDW envelope and step size - STEP SIZES MUST BE in Angstroms C
read(in,*) factor,factor2,xlen,ylen,zlen C
C*********************************************************************** c
c	read in the numbers of charges to be fit and the atoms numbers
c	at which the charges should be centered
c
read(in,*) nfit
read(in,*) (isfit(i),i=1,nfit)
c
c********************************************************************* 
IF(IWWRIT .NE. 1) GOTO 170
WRITE(IOUT,8000)(C(1,I),C(2,I),C(3,I),I=1,natoms) 8000 FORMAT(/1X,'Coordinates'/(1X,3F12.6)) 
WRITE(IOUT,8020) natoms,icharge,multip,nae,nbe,ne,nbasis $ ,nshells,maxtyp
8020 FORMAT(/1X,'natoms = ',I3
$/1X,'icharge = ',I3
$/1X,'multip = ',I3
$/1X,'nae = ',I3
$/1X,'nbe = ',I3
$/1X,'ne	= ',I3
$/1X,'nbasis = ',I3
$/1x,'nshells = ',i3
$/1x,'maxtyp = ',i3)
WRITE(IOUT,8030) (ian(I),I=1,natoms)
8030 FORMAT(/1X,'ian'/(1X,20I3))
WRITE(IOUT,8050) (jan(I),shellt(I),shella(I),I=1,NSHELL) 8050 FORMAT(/1X,'center type shella'/(1X,3I7)) 
WRITE(IOUT,8055) shella(NSHELL+1)
8055 FORMAT(1X,14X,I7)
WRITE(IOUT,8060) (EXX(I),C1(I),C2(I),I=1,nshells) 8060 FORMAT(/1X,12X,'EXPON',8X,'EXPCOF(S)',8X,'EXPCOF(P)', 
$/(1X,3E17.9))
write(iout,8070) (c3(i),cf(i),i=1,nshells) 8070 FORMAT(/1X,12X,'EXPCOF(D)',8X,'EXPCOF(F)', 
$/(1X,2E17.9))
write (iout,3010) nd, dipx, dipy, dipz
3010 format(/1x,'# of d functions',i1,/,6x,'dipole-x',f6.3,6x,' 
1	dipole-y',f6.3,6x,'dipole-z'f6.3)
write (iout,3575)
3575 format(/,15x,'atom #',15x,'v.d.w. radii') 
do 3500 i=1,natoms
write (iout,4000)i,vdwr(i)
4000 format (15x,i5,20x,f5.2)
3500 continue
write (iout,4500) maxpts
4500 format(/2x,'number of points to fit',i6) 
170 CONTINUE

C
C****************************************************************************** C
return
C
1020 format(A8)
1030 format(A80)
9002 format(8I4,1X)
END
C
C****************************************************************************** c
subroutine ep
c
c	routine to calculate the electrostatic potential from first order
c	perturbation theory
c
c	M.M. Francl April 1985
c	modified version of a MEPhisto routine
c	restricted to closed shell molecules
c	uep is version for unrestricted Hartree-Fock wavefunctions
c
implicit real*8 (a-h,o-z)
integer*4 shella,shelln,shellt,aos,shellc,aon,handle C
common /io/ in,iout
common /ipo/ ipo(10)
common /b/ exx(1000),c1(1000),c2(1000),c3(2500), $	x(1000),y(1000),z(1000),
$	jan(1000),shella(1000),shelln(1000),shellt(1000),
$	shellc(1000),aos(1000),aon(1000),nshells,maxtyp
common /mol/ natoms,icharge,multip,nae,nbe,nel,nbasis,ian(101), $	atmchg(100),c(3,100)
common /points/ p(3,50000),maxpnts
common /elp/ elecp(50000)
common /charge/ coef_alpha(9000),coef_beta(9000),iUHF common /xout/ title1(10),title2(10),title3(10),arcfil,dipole(3), $	q(104),i6to5,rms,percent
C
dimension hpert(45150),index(300)
C
data iptchg/1.0/
DATA ZERO/0.0/, TWO/2.0/, vnucmax/30.0/
C
c	divert to routine uep if wavefunction is unrestricted
c	Hartree-Fock wavefunction
c
if (iUHF .eq. 1) then
call uep
return
end if
c
c	initialize timing
c
handle = 0
C
c	set up the indexing table for hpert
c
do 100 i=1,nbasis
index(i) = (i-1)*i/2
100 continue
C
c	begin loop to calculate electrostatic potential
C
nocc = nel / 2
c
do 200 npnt=1,maxpnts
x1 = p(1,npnt)
x2 = p(2,npnt)
x3 = p(3,npnt)
C
C	CALCULATE THE ONE-ELECTRON INTEGRALS
C
if (ipo(7).eq.1) then
write(iout,3010)
3010 format(1x,'Time for integrals')
end if
c
CALL INTGRL (HPERT,X1,X2,X3,IPTCHG,i6to5) c
if (ipo(4).eq.1) call linout (hpert,nbasis,0,0) c
if (ipo(7).eq.1) then
write(iout,3000)
3000 format(1x,'Time for transform')
end if
C
C	FORM THE HPERT MATRIX ELEMENTS
C
e = zero
ICOEFI = -NBASIS
C
C	sum over occupied MOs
C
DO 220 II=1,NOCC
ICOEFI = ICOEFI + NBASIS
c
c	Calculate electrostatic potential
c
DO 221 IP=1,NBASIS
cpi = coef_alpha(icoefi+ip)
ipdex = index(ip)
c
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)) 223 continue
c
221 continue
220 CONTINUE
C
C	Calculate nuclear part of electrostatic potential
C
vnuc = zero
do 300 iatom=1,natoms
del1 = c(1,iatom) - x1
del2 = c(2,iatom) - x2
del3 = c(3,iatom) - x3
ra = dsqrt(del1*del1 + del2*del2 + del3*del3) if (ra.eq.zero) then
vnuc=vnucmax
goto 310
end if
vnuc = vnuc + ian(iatom) / ra
300 continue
310 continue
c
elecp(npnt) = (e * two + vnuc * iptchg)
if (ipo(7) .eq. 1) write(iout,*) 'e(',npnt,') = ',e 200 CONTINUE
END
subroutine uep
c
c	routine to calculate the electrostatic potential from first order
c	perturbation theory
c
c	M.M. Francl July 1985
c	modified version of a MEPhisto routine
c	restricted to unrestricted Hartree-Fock wavefunctions
c
implicit real*8 (a-h,o-z)
integer*4 shella,shelln,shellt,aos,shellc,aon,handle C
common /io/ in,iout
common /ipo/ ipo(10)
common /b/ exx(1000),c1(1000),c2(1000),c3(2500), $	x(1000),y(1000),z(1000),
$	jan(1000),shella(1000),shelln(1000),shellt(1000),
$	shellc(1000),aos(1000),aon(1000),nshells,maxtyp
common /mol/ natoms,icharge,multip,nae,nbe,nel,nbasis,ian(101), $	atmchg(100),c(3,100)
common /points/ p(3,50000),maxpnts
common /elp/ elecp(50000)
common /charge/ coef_alpha(9000),coef_beta(9000),iUHF common /xout/ title1(10),title2(10),title3(10),arcfil,dipole(3), $	q(104),i6to5,rms,percent
C
dimension hpert(45150),index(300)
C
data iptchg/1.0/
DATA ZERO/0.0/, TWO/2.0/, vnucmax/30.0/
C
c	intialize timing
c
handle = 0
C
c	set up the indexing table for hpert
c
do 100 i=1,nbasis
index(i) = (i-1)*i/2
100 continue
C
c	begin loop to calculate electrostatic potential
C
do 200 npnt=1,maxpnts
x1 = p(1,npnt)
x2 = p(2,npnt)
x3 = p(3,npnt)
C
C	CALCULATE THE ONE-ELECTRON INTEGRALS
C
if (ipo(7).eq.1) then
write(iout,3010)
3010 format(1x,'Time for integrals')
end if
c
CALL INTGRL (HPERT,X1,X2,X3,IPTCHG,i6to5) c
if (ipo(4).eq.1) call linout (hpert,nbasis,0,0) c
if (ipo(7).eq.1) then
write(iout,3000)
3000 format(1x,'Time for transform')
end if
C
C	FORM THE HPERT MATRIX ELEMENTS
C
c	alpha code
c
e = zero
ICOEFI = -NBASIS
C
C	sum over occupied alpha MOs
C
DO 220 II=1,nae
ICOEFI = ICOEFI + NBASIS
c
c	Calculate electrostatic potential
c
DO 221 IP=1,NBASIS
cpi = coef_alpha(icoefi+ip)
ipdex = index(ip)
c
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)) 223 continue
c
221 continue
220 CONTINUE
C
c	beta code
c
ICOEFI = -NBASIS
C
C	sum over occupied beta MOs
C
DO 420 II=1,nbe
ICOEFI = ICOEFI + NBASIS
c
c	Calculate electrostatic potential
c
DO 421 IP=1,NBASIS
cpi = coef_beta(icoefi+ip)
ipdex = index(ip)
c
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)) e = e + cpi * coef_* hpert(ip+index(iq)) 423 continue
c
421 continue
420 CONTINUE
C
C	Calculate nuclear part of electrostatic potential
C
vnuc = zero
do 300 iatom=1,natoms
del1 = c(1,iatom) - x1
del2 = c(2,iatom) - x2
del3 = c(3,iatom) - x3
ra = dsqrt(del1*del1 + del2*del2 + del3*del3) if (ra.eq.zero) then
vnuc=vnucmax
goto 310
end if
vnuc = vnuc + ian(iatom) / ra
300 continue
310 continue
c
elecp(npnt) = (e * two + vnuc * iptchg)
if (ipo(7) .eq. 1) write(iout,*) 'e(',npnt,') = ',e 200 CONTINUE
END
subroutine fit
c
c	routine to use method of Lagrange multipliers to obtain best
c	least square fit with constraints
c
c	M.M. Francl
c	April 1985
c
implicit real*8 (a-h,o-z)
integer*4 which1
character*40 arcfil
c
common /io/ in,iout
common /ipo/ ipo(10)
common /elp/ e(50000)
common /points/ p(3,50000),maxpnts
common /xout/ title1(10),title2(10),title3(10),arcfil,dipole(3), $	x(104),i6to5,rms,percent
common /mol/ nats,icharge,multip,nae,nbe,nel,nbasis,ian(101), $	atmchg(100),c(3,100)
common /tobefit/ isfit(100),nfit
c
dimension a(104,104),y(104),is(2,104),iad1(104),iad2(104) dimension d(104),which1(3),sx(100)
c
c	debye = conversion from debyes to au
c
data one/1.0/, zero/0.0/, debye/0.393427328/, maxdim/104/ data au2cal/627.51/, half/0.5/, hundred/100.0/ c
c	-DIPOLE CONSTRAINT TURNED OFF-
c
c	determine number of dipole constraints to impose. If molecule
c	is planar, since the charges are constrained to be placed at
c	atomic centers, by definition the component of the dipole in
c	the direction perpendicular to the plane will be zero
c
c	call constraints (nconstr,which1)
nconstr = 1
c
c	set up number of atoms to be fit
c
natoms = nfit
c
c	set up matrix of linear coefficients, A
c
c	begin loop over rows
c
do 100 k=1,natoms
kat = isfit(k)
c
c	begin loop over columns
c
do 200 mu=1,natoms
muat = isfit(mu)
c
sum = zero
do 400 i=1,maxpnts
rik = (p(1,i)-c(1,kat))**2 + (p(2,i)-c(2,kat))**2 + $	(p(3,i)-c(3,kat))**2
rik = dsqrt(rik)
rimu = (p(1,i)-c(1,muat))**2 + (p(2,i)-c(2,muat))**2 + $	(p(3,i)-c(3,muat))**2
rimu = dsqrt(rimu)
sum = sum + one / (rik * rimu)
400 continue
c
a(k,mu) = sum
200 continue
c
c	fill out columns corresponding to Lagrange multipliers
c
a(k,natoms+1) = half
c
c	disregard coordinate information if corresponding dipole
c	constraint is not being used
c
cdip	jcon = natoms + 1
cdip	do 300 icon=1,3
cdip	if (which1(icon) .eq. 1) goto 300
cdip	jcon = jcon + 1
cdip	a(k,jcon) = c(icon,k) * half
cdip 300 continue
c
100 continue
c
c	fill out the rows corresponding to constraints
c
do 500 mu=1,natoms
a(natoms+1,mu) = one
c
c	don't include data for rows corresponding to dipole moment
c	components which are not being considered explicitly as
c	constraints
c
cdip	jcon = natoms + 1
cdip	do 550 icon=1,3
cdip	if (which1(icon) .eq. 1) goto 550
cdip	jcon = jcon + 1
cdip	a(jcon,mu) = c(icon,mu)
cdip 550 continue
c
500 continue
c
c	fill out the block which connects Lagrange multipliers to
c	constraints
c
cdip	do 600 k=natoms+1,natoms+nconstr
cdip	do 600 mu=natoms+1,natoms+nconstr
cdip	a(k,mu) = zero
cdip 600 continue
c
a(natoms+1,natoms+1) = zero
c
c****debug*****
c
if (ipo(3) .eq. 1) then
write(iout,*) 'A matrix'
do 699 k=1,natoms+nconstr
write(iout,1699) (a(k,mu),mu=1,natoms+nconstr) 1699 format(1x,10f10.4)
699 continue
end if
c***************
c
c	construct column vector, Y
c
do 700 k=1,natoms
kat = isfit(k)
sum = zero
do 710 i=1,maxpnts
rik = (p(1,i)-c(1,kat))**2 + (p(2,i)-c(2,kat))**2 + $	(p(3,i)-c(3,kat))**2
rik = dsqrt(rik)
sum = sum + e(i) / rik
710 continue
y(k) = sum
if (ipo(3) .eq. 1) write(iout,*) k,y(k)
700 continue
c
c	construct the portion of y corresponding to Lagrange multipliers
c
c	convert dipole moment to atomic units and reverse it, the
c	convention for dipole moments in Gaussian 82 is opposite that of
c	the actual charge distribution
c
y(natoms+1) = dfloat(icharge)
c
c	disregard dipole moment components which are not being used for
c	explicit constraints
c
cdip	jcon = natoms + 1
c	do 750 icon=1,3
c	if (which1(icon) .eq. 1) goto 750
c	jcon = jcon + 1
c	y(jcon) = -dipole(icon) * debye
c 750 continue
c
if (ipo(3) .eq. 1)
$ write(iout,*) 'col vectr y', (y(kk),kk=1,natoms+nconstr) c
c	solve matrix equation Ax = y;
c	where x = (q1,q2, ... qn,l1,l2, ... ,ln)
c
c	x = A(inv)y
c
c	invert A
c
call inv(a,natoms+nconstr,is,iad1,iad2,d,maxdim) c
c****debug*****
c
if (ipo(3) .eq. 1) then
write(iout,*) 'A inverse'
do 799 k=1,natoms+nconstr
write(iout,1699) (a(k,mu),mu=1,natoms+nconstr) 799 continue
end if
c**************
c
c	perform matrix multiplication A(inv)y
c
call multay(a,y,x,natoms+nconstr,maxdim) c
if (ipo(3) .eq. 1) then
write(iout,*) 'Charges: '
do 899 i=1,natoms
write(iout,*) isfit(i),ian(isfit(i)),x(i) 899 continue
end if
c	compute rms deviation and mean absolute % deviation
c	of electrostatic potential
c
rms = zero
percent = zero
do 800 i=1,maxpnts
eq = zero
do 810 j=1,natoms
jat = isfit(j)
dist = (p(1,i)-c(1,jat))**2 + (p(2,i)-c(2,jat))**2 + $	(p(3,i)-c(3,jat))**2
dist = dsqrt(dist)
eq = eq + x(j) / dist
810 continue
rms = rms + (e(i) - eq)**2
percent = percent + dabs((e(i) - eq) / e(i) * hundred) cc	if (ipo(3) .eq. 1) write(iout,*) 'actual,calc ',e(i),eq
800 continue
if (ipo(3) .eq. 1) write(iout,*) 'sum of squares ',rms resid = dsqrt( rms / (maxpnts - 2))
if (ipo(3) .eq. 1) write(iout,*) 's y/x ',resid rms = dsqrt(rms / maxpnts) * au2cal
percent = percent / maxpnts
if (ipo(3) .eq. 1) write(iout,*) 'rms, %',rms,percent c
c	calculate the standard deviations on the charges
c
do 10 j=1,natoms
jat = isfit(j)