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

c
sumx = zero
do 11 i=1,maxpnts
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)
sumx = sumx + 1.0d0/dist
11 continue
c
sumx2 = zero
do 12 i=1,maxpnts
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)
sumx2 = sumx2 + (dist * dist)
12 continue
sx(j) = resid*resid*maxpnts / (maxpnts*sumx2 - sumx*sumx) sx(j) = dsqrt( sx(j) )
C	here
c	write(iout,*) 's q',j,sx(j)
10 continue
c
c
return
end
subroutine multay(a,y,x,n,maxdim)
c
c	matrix multiplication routine
c
implicit real*8 (a-h,o-z)
c
dimension a(maxdim,maxdim),y(maxdim),x(maxdim) c
data zero/0.0/
c
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
end
subroutine constraints(nconstr,which1)
c
c	subroutine to determine if one of the dipole moment constraints
c	is degenerate
c
c	M.M. Francl
C	April 1985
c
implicit real*8 (a-h,o-z)
integer*4 which1
c
common /mol/ natoms,icharge,multip,nae,nbe,nel,nbasis,ian(101), $	atmchg(100),c(3,100)
dimension which1(3)
data tol/1.0e-6/
c
nconstr = 4
c
c	are all coordinates zero?
c
do 200 i=1,3
which1(i) = 0
do 100 j=1,natoms
if (c(i,j) .gt. tol) goto 200
100 continue
c
c	they're all zero (or close enough)
c
which1(i) = 1
nconstr = nconstr - 1
200 continue
return
end
SUBROUTINE INV(A,N,IS,IAD1,IAD2,D,MDM)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C	******************************************************************
C	INVERSION OF SQUARE MATRIX A BY MEANS OF THE GAUSS-JORDAN
C	ALGORITHM
C
C	APRIL 72/RS9B
C	******************************************************************
DIMENSION A(MDM,MDM),IS(2,MDM),IAD1(MDM),IAD2(MDM),D(MDM) C
COMMON/IO/IN,IOUT,IPUNCH
C
DATA ZERO/0.0D0/, ONE/1.0D0/, SMALL/1.0D-20/ C
2000 FORMAT(' WARNING FROM INV: MATRIX IS SINGULAR') C	******************************************************************
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=dABS(A(L,M))
IF(E.LT.B) GOTO 8
I=L
K=M
8 b=dMAX1(b,e)
2 CONTINUE
IS(1,I)=1
IS(2,K)=1
IAD1(K)=I
IAD2(I)=K
B=A(I,K)
C.....PIVOT
IF(dABS(B).LT. SMALL) GOTO 20
A(I,K)=ONE/B
DO 6 L=1,N
IF(L.EQ.K) GOTO 6
C.....KELLERZEILE
A(I,L)=-A(I,L)/B
6 CONTINUE
DO 5 L=1,N
DO 5 M=1,N
IF(L.EQ.I.OR.M.EQ.K) GOTO 5
C.....RECHTECK-REGEL
A(L,M)=A(L,M)+A(L,K)*A(I,M)
5 CONTINUE
DO 11 L=1,N
IF(L.EQ.I) GOTO 11
C.....PIVOT-SPALTE
A(L,K)=A(L,K)/B
11 CONTINUE
9 CONTINUE
C.....PERMUTATION DER ZEILEN, UM DIE NATUERLICHE ORDNUNG WIEDER HERZUSTE 
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)
15 CONTINUE
C.....PERMUTATION DER SPALTEN
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
RETURN
C
C	ERROR EXIT: MATRIX IS SINGULAR
20 WRITE(IOUT,2000)
stop 'inv in polar'
END

c
c
subroutine out
c
c
c	L.E. Chirlian
c	April 1985
c
c	A subroutine to output the charges and other pertinant
c	information from the CHELP program
c
implicit real*8 (a-h,o-z)
integer*4 shella,shelln,shellt,shellc,aos,aon,shlAdf character*8 title1,title2,title3
character*40 arcfil
c
common /io/ in,iout
common /mol/ natoms,icharge,multip,nae,nbe,ne,nbasis,ian(101), 1	atmchg(100),c(3,100)
common /xout/ title1(10),title2(10),title3(10),arcfil, 1	dipx, dipy, dipz, q(104), nd, rms, pd
common /points/ p(3,50000), maxpnts
c
c	create output
c
write (iout,100)
100 format (/,/,22X,'Charges from Electrostatic Potentials') 
write (iout,110)
110 format (/,38x,'CHELP')
c
c	write title
1010 format(9a8,a7)
if (title2(1) .eq. ' ') then
goto 22
end if
if (title3(1) .eq. ' ') then
goto 23
end if
write (6,1113)title1
write (6,1113)title2
write (6,1113)title3
goto 24
22 continue
write (6,1113)title1
1113 format (9a8,a7)
goto 24
23 continue
write (6,1113)title1
write (6,1113)title2
goto 24
24 continue

c	write checkpoint file name
c
write (iout,150)arcfil
150 format(/2x,'checkpoint file: ',a40)
c
c	print date
c
c************************************************************************** c	write geometry
c************************************************************************** c
write (iout,180)
180 Format (/2x,36x,'MOLECULAR GEOMETRY') 
write (iout,190)
190 Format (/,/,17x,'Atomic Number',8x,'X',12x,'Y',12x,'Z') 
do 30 i=1,natoms
write (iout,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 (iout,210)icharge
210 format (/2x,'The total charge is constrained to: ',i1) 
write (iout,230)dipx,dipy,dipz
230 format (/2x,'Dipole moment constrained to: x='f9.5,2x,'y=',f9.5 
1	,2x,'z=',f9.5)
write (iout,240)
240 format (/,36x,'NET CHARGES')
write (iout,250)
250 format (/,28x,'Atomic number',5x,'Charge') 
write (iout,260)(ian(i),q(i),i=1,natoms) 260 format (/,34x,i2,8x,f8.4)
write (iout,270)maxpnts
270 format(/,2x,'Fit to electrostatic potential at ',i5,' points') 
write (iout,280)rms
280 format (/,2x,'Root mean square deviation is ',f10.2,' kcal/mole') 
return
end


SUBROUTINE INTGRL (H,X1,X2,X3,ICHARG,I6TO5) C
C	ROUTINE TO CALCULATE THE ELECTRON-CHARGE MATRIX ELEMENTS FOR THE
C	POLARIZATION POTENTIAL. CODE REVISED FROM THE ONE ELECTRON PACKAGE
C	AS IT EXISTED AUGUST, 1983.
C
C
C	REVISED BY M.M. FRANCL JANUARY 1984 FOR PRINCETON CHEMISTRY
C	DEPARTMENT VAX 11/780
C
C	REVISED TO BE COMPATIBLE WITH COMMON /B/ FROM GAUSSIAN 82
C	MAY 1984 M.M. FRANCL
C
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER*4 SHELLA,SHELLN,SHELLT,SHELLC,AOS,AON,SHLADF C
COMMON /IO/ IN,IOUT
COMMON /IPO/ IPO(10)
COMMON /GG/ GA(35),FM(5),rpitwo,tol
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),NSHELL,MAXTYP
COMMON /MOL/ NATOMS,JCHARG,MULTIP,NAE,NBE,NEL,NBASIS,IAN(101), $	ATMCHG(100),C(3,100)
C
DIMENSION H(45150)
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) 
DIMENSION EPN(100)
C
COMMON/H100/
$EP00,EP10,EP20,EP30,EP40,EP50,EP60,EP70,EP80,EP90,	SPDSTV
$EP01,EP11,EP21,EP31,EP41,EP51,EP61,EP71,EP81,EP91,	SPDSTV
$EP02,EP12,EP22,EP32,EP42,EP52,EP62,EP72,EP82,EP92,	SPDSTV
$EP03,EP13,EP23,EP33,EP43,EP53,EP63,EP73,EP83,EP93,	SPDSTV
$EP04,EP14,EP24,EP34,EP44,EP54,EP64,EP74,EP84,EP94,	SPDSTV
$EP05,EP15,EP25,EP35,EP45,EP55,EP65,EP75,EP85,EP95,	SPDSTV
$EP06,EP16,EP26,EP36,EP46,EP56,EP66,EP76,EP86,EP96,	SPDSTV
$EP07,EP17,EP27,EP37,EP47,EP57,EP67,EP77,EP87,EP97,	SPDSTV
$EP08,EP18,EP28,EP38,EP48,EP58,EP68,EP78,EP88,EP98,	SPDSTV
$EP09,EP19,EP29,EP39,EP49,EP59,EP69,EP79,EP89,EP99	SPDSTV
C
DIMENSION EEP(100)	SPDSTV
DIMENSION MAX(6)
C	SPDSTV
C	LOCAL VARIABLES.	SPDSTV
C
DIMENSION AG(6),CSA(6),CPA(6),CDA(6),	SPDSTV
$	BG(6),CSB(6),CPB(6),CDB(6),	SPDSTV
$	DPP(9)	SPDSTV
EQUIVALENCE(OF0,OF(1)),(OF1,OF(2)),(OF2,OF(3)),	SPDSTV
$	(OF3,OF(4)),(OF4,OF(5)),(OF5,OF(6)),	SPDSTV
$	(OF6,OF(7)),(OF7,OF(8)),(OF8,OF(9))	SPDSTV
EQUIVALENCE(OX0,OX(1)),(OX1,OX(2)),(OX2,OX(3)),	SPDSTV
$	(OX3,OX(4)),(OX4,OX(5)),(OX5,OX(6)),	SPDSTV
$	(OX6,OX(7)),(OX7,OX(8)),(OX8,OX(9))	SPDSTV
EQUIVALENCE(A1,A(2)),(A2,A(3)),(A3,A(4)),(A4,A(5))	SPDSTV
EQUIVALENCE(B1,B(2)),(B2,B(3)),(B3,B(4)),(B4,B(5))	SPDSTV
EQUIVALENCE(T01,T0),(T02,T1),(T03,T2),	SPDSTV
$	(T04,T3),(T05,T4),(T06,T5),	SPDSTV
$	(T07,T6),(T08,T7),(T09,T8)	SPDSTV
EQUIVALENCE(T10,TX(10)),(T11,TX(11)),(T12,TX(12)),(T13,TX(13))	SPDSTV
EQUIVALENCE(T0,TX(1)),(T1,TX(2)),(T2,TX(3)),	SPDSTV
$	(T3,TX(4)),(T4,TX(5)),(T5,TX(6)),	SPDSTV
$	(T6,TX(7)),(T7,TX(8)),(T8,TX(9))	SPDSTV
EQUIVALENCE(C001,T01),(C050,T02),(C054,T09),	SPDSTV
$	(C067,T13),(C068,T08),(C074,T03)	SPDSTV
EQUIVALENCE(ABX1,ABX(2)),(ABX2,ABX(3)),	SPDSTV
$	(ABX3,ABX(4)),(ABX4,ABX(5))	SPDSTV
EQUIVALENCE(AB004,ABX1),(AB006,ABX2),(AB023,ABX3),(AB029,ABX4)	SPDSTV
EQUIVALENCE(ABY1,ABY(2)),(ABY2,ABY(3)),	SPDSTV
$	(ABY3,ABY(4)),(ABY4,ABY(5))	SPDSTV
EQUIVALENCE(AB007,ABY1),(AB010,ABY2),(AB032,ABY3),(AB035,ABY4)	SPDSTV
EQUIVALENCE(ABZ1,ABZ(2)),(ABZ2,ABZ(3)),	SPDSTV
$	(ABZ3,ABZ(4)),(ABZ4,ABZ(5))	SPDSTV
EQUIVALENCE(AB002,ABZ1),(AB003,ABZ2),(AB011,ABZ3),(AB017,ABZ4)	SPDSTV
EQUIVALENCE(ABSQ1,ABSQ(2)),(ABSQ2,ABSQ(3)),	SPDSTV
$	(ABSQ3,ABSQ(4)),(ABSQ4,ABSQ(5))	SPDSTV
EQUIVALENCE(APB1,APB(2)),(APB2,APB(3)),	SPDSTV
$	(APB3,APB(4)),(APB4,APB(5))	SPDSTV
EQUIVALENCE(CPX1,CPX(2)),(CPX2,CPX(3)),	SPDSTV
$	(CPX3,CPX(4)),(CPX4,CPX(5))	SPDSTV
EQUIVALENCE(CPY1,CPY(2)),(CPY2,CPY(3)),	SPDSTV
$	(CPY3,CPY(4)),(CPY4,CPY(5))	SPDSTV
EQUIVALENCE(CPZ1,CPZ(2)),(CPZ2,CPZ(3)),	SPDSTV
$	(CPZ3,CPZ(4)),(CPZ4,CPZ(5))	SPDSTV
EQUIVALENCE(F1,F(2)),(F2,F(3)),(F3,F(4)),(F4,F(5))	SPDSTV
EQUIVALENCE(FM0,FM(1)),(FM1,FM(2)),(FM2,FM(3)),(FM3,FM(4)),	SPDSTV
$	(FM4,FM(5))	SPDSTV
EQUIVALENCE (D001,FM0)	SPDSTV
EQUIVALENCE(EP00,EEP(1))	SPDSTV
C
DATA MAX/1,4,9,1,4,10/	SPDSTV
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/,HALF/0.5/,ONE/1.0/,ONEPT5/1.5/,TWO/2.0/,THREE/3.0/, *ROOT3/1.732050808/,PI/3.14159265358979/ 
DATA ANTOAU /1.889726878D0/
C
2010 FORMAT(/1X,'ELECTRON-CHARGE MATRIX ELEMENTS'/) C
C*********************************************************************** SPDSTV C	INITIALIZE THIS SEGMENT.	SPDSTV
C*********************************************************************** SPDSTV C	SPDSTV
C	****************************************************************** SPDSTV
C	COMPUTE SIZE OF S T AND V ARRAYS
C	****************************************************************** SPDSTV
NTT=(NBASIS*(NBASIS+1))/2
I5OR6=3
CC IF(IGO(4) .NE. 0) I5OR6 = 0
C	****************************************************************** SPDSTV
C	INITIALIZE RENORM USED TO NORMALIZE D FUNCTIONS
C	****************************************************************** SPDSTV
DO 100 I=1,10	SPDSTV
100 RENORM(I)=ONE	SPDSTV
RENORM(5)=ROOT3	SPDSTV
RENORM(8)=ROOT3	SPDSTV
RENORM(9)=ROOT3	SPDSTV
C	****************************************************************** SPDSTV
C	CLEAR H ARRAY
C	****************************************************************** SPDSTV
DO 50 I=1,NTT	SPDSTV
50 H(I)=ZERO
C	****************************************************************** SPDSTV
C	* INITIALIZE THE VARIABLES USED BY ROUTINE FMGEN.	* SPDSTV
C	****************************************************************** SPDSTV
CALL FMSET
DO 95 I=1,5
95 FM(I)=ZERO
ABX(1)=ONE	SPDSTV
ABY(1)=ONE	SPDSTV
ABZ(1)=ONE	SPDSTV
A(1)=ONE	SPDSTV
B(1)=ONE	SPDSTV
F(1)=ONE	SPDSTV
CPX(1)=ONE	SPDSTV
CPY(1)=ONE	SPDSTV
CPZ(1)=ONE	SPDSTV
APB(1)=ONE	SPDSTV
ABSQ(1)=ONE	SPDSTV
C*********************************************************************** SPDSTV C	LOOP OVER SHELLS ISHELL AND JSHELL.	SPDSTV
C*********************************************************************** SPDSTV 
DO 1000 ISHELL=1,NSHELL	SPDSTV
DO 1000 JSHELL=1,ISHELL	SPDSTV
SYMFAC = ONE
C	****************************************************************** SPDSTV
C	ZERO LOCATIONS	SPDSTV
C	****************************************************************** SPDSTV
80 CONTINUE
DO 9447 JI=1,100	SPDSTV
EPN(JI)=ZERO	SPDSTV
9447 CONTINUE	SPDSTV
IF(SHELLT(ISHELL)-SHELLT(JSHELL))120,120,110	SPDSTV
110 INEW=JSHELL	SPDSTV
JNEW=ISHELL	SPDSTV
LA=SHELLT(JSHELL)	SPDSTV
LB=SHELLT(ISHELL)	SPDSTV
GO TO 200	SPDSTV
120 INEW=ISHELL	SPDSTV
JNEW=JSHELL	SPDSTV
LA=SHELLT(ISHELL)	SPDSTV
LB=SHELLT(JSHELL)	SPDSTV
200 CONTINUE	SPDSTV
LAP1=LA+1	SPDSTV
LBP1=LB+1	SPDSTV
LAMAX=MAX(LAP1+I5OR6)	SPDSTV
LBMAX=MAX(LBP1+I5OR6)	SPDSTV
ITYPE=3*LB+LA	SPDSTV
M=LA+LB+1	SPDSTV
NGA=SHELLN(INEW)	SPDSTV
NGB=SHELLN(JNEW)	SPDSTV
AX=X(INEW)	SPDSTV
BX=X(JNEW)	SPDSTV
AY=Y(INEW)	SPDSTV
BY=Y(JNEW)	SPDSTV
AZ=Z(INEW)	SPDSTV
BZ=Z(JNEW)	SPDSTV
ISHA=SHELLA(INEW)	SPDSTV
ISHB=SHELLA(JNEW)	SPDSTV
ISHAD = SHLADF(INEW)
ISHBD = SHLADF(JNEW)
IAOS=AOS(INEW)	SPDSTV
JAOS=AOS(JNEW)	SPDSTV
C	****************************************************************** SPDSTV
C	OBTAIN INFORMATION ABOUT SHELLS INEW AND JNEW	SPDSTV
C	****************************************************************** SPDSTV
DO 101 I=1,NGA	SPDSTV
N=ISHA+I-1	SPDSTV
ND = ISHAD + I -1
IF (MAXTYP .LE. 1) ND=1
AG(I)=EXX(N)	SPDSTV
CSA(I)=C1(N)	SPDSTV
CPA(I)=C2(N)	SPDSTV
101 CDA(I)=C3(ND)	SPDSTV
DO 102 I=1,NGB	SPDSTV
N=ISHB+I-1	SPDSTV
ND = ISHBD + I -1
BG(I)=EXX(N)	SPDSTV
CSB(I)=C1(N)	SPDSTV
CPB(I)=C2(N)	SPDSTV
102 CDB(I)=C3(ND)	SPDSTV
ABX(2)=BX-AX	SPDSTV
ABY(2)=BY-AY	SPDSTV
ABZ(2)=BZ-AZ	SPDSTV
RABSQ=ABX(2)*ABX(2)+ABY(2)*ABY(2)+ABZ(2)*ABZ(2)	SPDSTV
ABSQ(2)=RABSQ	SPDSTV
DO 103 I=3,5	SPDSTV
ABX(I)=ABX(I-1)*ABX(2)	SPDSTV
ABY(I)=ABY(I-1)*ABY(2)	SPDSTV
ABZ(I)=ABZ(I-1)*ABZ(2)	SPDSTV
103 ABSQ(I)=ABSQ(I-1)*ABSQ(2)	SPDSTV
AB001=ONE	SPDSTV
AB005=ABX1*ABZ1	SPDSTV
AB008=ABY1*ABZ1	SPDSTV
AB009=ABX1*ABY1	SPDSTV
AB012=ABX1*ABZ2	SPDSTV
AB013=ABX2*ABZ1	SPDSTV
AB014=ABY1*ABZ2	SPDSTV
AB015=ABX1*ABY1*ABZ1	SPDSTV
AB016=ABY2*ABZ1	SPDSTV
AB018=ABX1*ABZ3	SPDSTV
AB019=ABX2*ABZ2	SPDSTV
AB020=ABY1*ABZ3	SPDSTV
AB021=ABX1*ABY1*ABZ2	SPDSTV
AB022=ABY2*ABZ2	SPDSTV
AB024=ABX2*ABY1	SPDSTV
AB025=ABX1*ABY2	SPDSTV
AB026=ABX3*ABZ1	SPDSTV
AB027=ABX2*ABY1*ABZ1	SPDSTV
AB028=ABX1*ABY2*ABZ1	SPDSTV
AB030=ABX3*ABY1	SPDSTV
AB031=ABX2*ABY2	SPDSTV
AB033=ABY3*ABZ1	SPDSTV
AB034=ABX1*ABY3	SPDSTV
C*********************************************************************** SPDSTV C	LOOP OVER GAUSSIANS (CONTRACTION LOOP).	SPDSTV
C*********************************************************************** SPDSTV 
DO 105 IGAUSS=1,NGA	SPDSTV
AA=AG(IGAUSS)	SPDSTV
DO 105 JGAUSS=1,NGB	SPDSTV
BB=BG(JGAUSS)	SPDSTV
AAPBB=AA+BB	SPDSTV
APBB=ONE/AAPBB	SPDSTV
F2=TWO*AA*BB*APBB	SPDSTV
PX=(AA*AX+BB*BX)*APBB	SPDSTV
PY=(AA*AY+BB*BY)*APBB	SPDSTV
PZ=(AA*AZ+BB*BZ)*APBB	SPDSTV
A(2)=ONE/AA	SPDSTV
B(2)=ONE/BB	SPDSTV
F(2)=F2	SPDSTV
APB(2)=APBB	SPDSTV
YX=PI*APBB	SPDSTV
EXX1=HALF*F2*RABSQ	SPDSTV
IF(EXX1-80.0E0)4172,4173,4173
4173 EXX1=ZERO
GO TO 4714	SPDSTV
4172 EXX1=EXP(-EXX1)
4714 CONTINUE	SPDSTV