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

EP05=OF6*(+C012*AB001*D001+C013*AB006*D001-C014*AB001*D001+C056*AB SPDSTV $004*D005+C001*AB001*D007-C050*AB001*D004)	SPDSTV
EP07=OF6*(+C013*AB008*D001+C006*AB007*D002+C006*AB002*D008+C001*AB SPDSTV $001*D009)	SPDSTV
EP08=OF6*(+C013*AB009*D001+C006*AB004*D008+C006*AB007*D005+C001*AB SPDSTV $001*D010)	SPDSTV
EP09=OF6*(+C012*AB001*D001+C013*AB010*D001-C014*AB001*D001+C056*AB SPDSTV $007*D008+C001*AB001*D011-C050*AB001*D004)	SPDSTV
IF(ITYPE-6)3261,3262,3261	SPDSTV
C	****************************************************************** SPDSTV
C	*	PP	* SPDSTV
C	****************************************************************** SPDSTV
3261 CONTINUE	SPDSTV
EP10=OF1*(+C002*AB002*D001-C001*AB001*D002)	SPDSTV
EP30=OF1*(+C002*AB004*D001-C001*AB001*D005)	SPDSTV
EP60=OF1*(+C002*AB007*D001-C001*AB001*D008)	SPDSTV
EP11=OF4*(-C007*AB003*D001+C008*AB001*D001+C006*AB002*D002-C002*AB SPDSTV $002*D002+C001*AB001*D003-C050*AB001*D004)	SPDSTV
EP31=OF4*(-C007*AB005*D001+C006*AB002*D005-C002*AB004*D002+C001*AB SPDSTV $001*D006)	SPDSTV
EP61=OF4*(-C007*AB008*D001+C006*AB002*D008-C002*AB007*D002+C001*AB SPDSTV $001*D009)	SPDSTV
EP13=OF4*(-C007*AB005*D001+C006*AB004*D002-C002*AB002*D005+C001*AB SPDSTV $001*D006)	SPDSTV
EP33=OF4*(-C007*AB006*D001+C008*AB001*D001+C006*AB004*D005-C002*AB SPDSTV $004*D005+C001*AB001*D007-C050*AB001*D004)	SPDSTV
EP63=OF4*(-C007*AB009*D001+C006*AB004*D008-C002*AB007*D005+C001*AB SPDSTV $001*D010)	SPDSTV
EP16=OF4*(-C007*AB008*D001+C006*AB007*D002-C002*AB002*D008+C001*AB SPDSTV $001*D009)	SPDSTV
EP36=OF4*(-C007*AB009*D001+C006*AB007*D005-C002*AB004*D008+C001*AB SPDSTV $001*D010)	SPDSTV
EP66=OF4*(-C007*AB010*D001+C008*AB001*D001+C006*AB007*D008-C002*AB SPDSTV $007*D008+C001*AB001*D011-C050*AB001*D004)	SPDSTV
3262 CONTINUE
DO 2137 I=1,100	SPDSTV
2137 EPN(I)=EPN(I)+EEP(I)	SPDSTV
105 CONTINUE	SPDSTV
C	****************************************************************** SPDSTV
C	END OF LOOP OVER GAUSSIANS	SPDSTV
C	STORE IN ARRAYS	SPDSTV
C	****************************************************************** SPDSTV
INTC=0	SPDSTV
DO 500 J=1,10	SPDSTV
R3B=RENORM(J)	SPDSTV
DO 500 I=1,10	SPDSTV
R3A=R3B*RENORM(I)	SPDSTV
INTC=INTC+1	SPDSTV
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) 1000 CONTINUE	SPDSTV
IF(IPO(4).EQ.0) GOTO 1500
WRITE(IOUT,2010)
CALL LINOUT(H,NBASIS,0,0)
1500 CONTINUE
RETURN
END	SPDSTV
SUBROUTINE MATFIL(A,AA,AOS,SHELLT,INEW,JNEW,LAMAX,LBMAX,LA,LB) C
C	GAUSSIAN 77/UCI
C
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER AOS(1000), SHELLT(1000)
C
DIMENSION A(45150),AA(100)
C
LIND(I)=(I*(I-1))/2
C
ISTART=AOS(INEW)	MATFIL
JSTART=AOS(JNEW)	MATFIL
IAL = 0
IAU = 5
IBL = 0
IBU = 5
IMA = 0
IMB = 0
IF(SHELLT(INEW) .EQ. 2) IMA = 1
IF(SHELLT(JNEW) .EQ. 2) IMB = 1
C	MATFIL
120 INTC=0	MATFIL
DO 170 J=1,LBMAX	MATFIL
DO 170 I=1,LAMAX	MATFIL
INTC=INTC+1	MATFIL
IF( LA.GT.1 .AND. I.GT.IAL .AND. I.LT.IAU )GO TO 170	MATFIL
IF( LA.EQ.1 .AND. I.EQ.IMA	)GO TO 170	MATFIL
IF( LB.GT.1 .AND. J.GT.IBL .AND. J.LT.IBU )GO TO 170	MATFIL
IF( LB.EQ.1 .AND. J.EQ.IMB	)GO TO 170	MATFIL
IND=ISTART+I-1	MATFIL
JND=JSTART+J-1	MATFIL
IF(IND-JND)130,140,150	MATFIL
130 IJ=LIND(JND)+IND	MATFIL
GO TO 160	MATFIL
140 IJ=LIND(IND+1)	MATFIL
GO TO 160	MATFIL
150 IJ=LIND(IND)+JND	MATFIL
160 A(IJ)=AA(INTC)	MATFIL
170 CONTINUE	MATFIL
RETURN	MATFIL
END	MATFIL
SUBROUTINE REDUC1(X,LAMAX,LBMAX,I6TO5)
C
C	MODIFIED FOR POLARIZATION POTENTIAL CALCULATIONS
C	M.M. FRANCL FEBRUARY 1984
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X(100),S(100),IND5(9),IND6(10) 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
C	******************************************************************
C	ROUTINE REORDERS FROM ARRANGEMENT: S,Z,ZZ,X,XZ,XX,Y,YZ,XY,YY
C	TO	S,X,Y,Z,XX,YY,ZZ,XY,XZ,YZ
C	OR FROM S,Z,X,Y TO S,X,Y,Z
C	THIS ENSURES LABELING COMPATIBILITY BETWEEN THE SP AND D PACKAGES
C	AT THE SAME TIME THE INTEGRALS ARE MOVED TO THE FIRST 1,4,10,16,
C	40 OR 100 LOCATIONS, DEPENDING ON THE SHELL QUANTUM NUMBERS
C	******************************************************************
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)
GO TO 160
C	******************************************************************
C	ROUTINE TO REDUCE SIX D FUNCTIONS TO FIVE
C	ALSO REORDERS FROM : S,Z,ZZ,X,XZ,XX,Y,YZ,YX,YY
C	TO S,X,Y,Z,ZZ,XX-YY,XY,XZ,YZ
C	OR FROM S,Z,X,Y TO S,X,Y,Z
C	FOR COMPATIBILITY WITH SP PACKAGE
C	******************************************************************
40 DO 150 I=1,LBMAX
ISB=10*(IND5(I)-1)
C	IFB=0 FOR S,X,Y,Z,XY,XZ,YZ, IFB=1 FOR ZZ-RR, IFB=2 FOR XX-YY
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
GOTO(130,122,123,124,125,126,127,128,129),IHOP C
C	******************************************************************
C	*	(F,O,ZA2)	*
C	******************************************************************
122 XX=ZZ1(X,ISA,3,7)
GO TO 140
C
C	******************************************************************
C	*	(F,O,XA2-YA2)	*
C	******************************************************************
123 XX=XY1(X,ISA,4)
GO TO 140
C
C	******************************************************************
C	*	(ZB2,O,F)	*
C	******************************************************************
124 XX=ZZ1(X,ISA,30,70)
GO TO 140
C
C	******************************************************************
C	*	(ZB2,O,ZA2)	*
C	******************************************************************
125 XX=ZZ1(X,ISA,30,70)-PT5*(ZZ1(X,ISA+3,30,70)+ZZ1(X,ISA+7,30,70)) 
GO TO 140
C
C	******************************************************************
C	*	(ZB2,O,XA2-YA2)	*
C	******************************************************************
126 XX=R3OV2*(ZZ1(X,ISA,30,70)-ZZ1(X,ISA+4,30,70)) 
GO TO 140
C
C	******************************************************************
C	*	(XB2-YB2,O,F)	*
C	******************************************************************
127 XX=XY1(X,ISA,40)
GO TO 140
C
C	******************************************************************
C	*	(XB2-YB2,O,ZA2)	*
C	******************************************************************
128 XX=R3OV2*(ZZ1(X,ISA,3,7)-ZZ1(X,ISA+40,3,7)) 
GO TO 140
C
C	******************************************************************
C	*	(XB2-YB2,O,XA2-YA2)	*
C	******************************************************************
129 XX=R3OV2*(XY1(X,ISA,4)-XY1(X,ISA+40,4)) 
GO TO 140
C
C	******************************************************************
C	*	(F,O,F)	*
C	******************************************************************
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
FUNCTION ZZ1(X,I,IX,IY)
C
DIMENSION X(100)
C
DATA HALF/0.5/
C
ZZ1=X(I)-HALF*(X(I+IX)+X(I+IY))
RETURN
END
FUNCTION XY1(X,I,IY)
C
DIMENSION X(100)
C
DATA HALFR3/0.8660254040/
C
XY1=HALFR3*(X(I)-X(I+IY))
RETURN
END
SUBROUTINE FMGEN(T,M)
C
IMPLICIT REAL*8 (A-H,O-Z)
COMMON/IO/IN,IOUT
COMMON/GG/ GA(35),FM(5),rpitwo,tol
C
C
EQUIVALENCE (APPROX,OLDSUM)
C
DATA ZERO/0.0E0/, HALF/0.5E0/, ONE/1.0E0/, TWO/2.0E0/, TEN/10.0E0/ $	,PI/3.14159265358979E0/, F42/42.0E0/, F80/80.0E0/
C
2001 FORMAT(42H1FAILURE IN FMGEN FOR SMALL T: IX.GT.50, / 
$ 6H IX = ,I3,7H, T = ,E20.14)
2002 FORMAT(37H1FAILURE IN FMGEN FOR INTERMEDIATE T,/ 
$ 6H T = ,E20.14)
C
TEXP=ZERO
IF(T-F80)2,3,3
2 TEXP=EXP(-T)
3 CONTINUE
IF(T-TEN)10,70,70
C*********************************************************************** C	0 .LT. T .LT. 10
C*********************************************************************** 
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(IOUT,2001)IX,T
STOP 'FMGEN'
40 SUM=SUM+TX/GA(IX)
IF(TOL-ABS(OLDSUM/SUM-ONE))50,60,60
50 OLDSUM=SUM
GO TO 20
60 FM(M)=SUM*TERM
GO TO 160
C
70 IF(T-F42)80,150,150
C*********************************************************************** C	10 .LE. T .LT. 42
C*********************************************************************** 
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(IOUT,2002)T
stop 'fmgen'
140 FM(M)=APPROX-FIMULT*SUM
GO TO 160
C*********************************************************************** C	T .GE. 42
C*********************************************************************** 
150 TX=FLOAT(M)-HALF
FM(M)=HALF*GA(M)/(T**TX)
C*********************************************************************** C	RECUR DOWNWARDS TO FM(1)
C*********************************************************************** 
160 TX=T+T
SUM=FLOAT(M+M-3)
MM1=M-1
IF(MM1)170,190,170
170 DO 180 IX=1,MM1
FM(M-IX)=(TX*FM(M-IX+1)+TEXP)/SUM
180 SUM=SUM-TWO
190 RETURN
C
ENTRY FMSET
C
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
SUBROUTINE LINOUT(X,N,KEY,IZERO)
C
C	GENERAL LINEAR MATRIX OUTPUT ROUTINE
C
C	KEY=0 MATRIX SYMMETRIC
C	KEY=1 MATRIX SQUARE ASYMMETRIC
C
C	IZERO=0 ZERO MATRIX ELEMENTS LESS THAN CUTOFF
C	IZERO=1 DO NOT ZERO MATRIX ELEMENTS
C
C	CUTOFF=1.0E-06
C
implicit real*8 (a-h,o-z)
COMMON/IO/IN,IOUT
C
DIMENSION S(9),X(45150)
C
DATA CUTOFF/1.0E-06/
DATA ZERO/0.0E0/
C
IA(I)=(I*(I-1))/2
C
C
ILOWER=1
100 IUPPER=MIN0(ILOWER+8,N)
IRANGE=MIN0(IUPPER-ILOWER+1,9)
WRITE (IOUT,9000) (J,J=ILOWER,IUPPER)
WRITE (IOUT,9010)
DO 160 I=1,N
K=1
DO 150 J=ILOWER,IUPPER
IF(KEY)110,120,110
110 IJ=N*(J-1)+I
GO TO 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 (IOUT,9020) I,(S(J),J=1,IRANGE) 
WRITE (IOUT,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

subroutine select
c
c	routine to select points from a (hopefully) randomly distributed set
c	for fitting to the electrostatic potential
c
c	Points which lie within the van der Waals envelope of the molecule
c	are excluded. FACTOR controls the size of the envelope.
c
c	M.M. Francl
c	April 1985
c
implicit real*8 (a-h,o-z)
integer*4 seed
c
common /io/ in,iout
common /mol/ natoms,icharg,multip,nae,nbe,nel,nbasis,ian(101), $	atmchg(100),c(3,100)
common /ipo/ ipo(10)
common /sphere/ vdw(104),ntotp,npoints,sph(3,3074) common /points/ p(3,50000), maxpnts
common /envelope/ factor,factor2,xlen,ylen,zlen c
dimension vdwout(104)
c
c
data seed/12345/, half/0.5/
data ang2au /1.889726878/
c
c************************************************************ c	convert radii to au, and size for inner and outer
c	cutoffs
c
do 10 i=1,natoms
vdwout(i) = vdw(i) *ang2au * factor2
vdw(i) = vdw(i) * ang2au * factor
10 continue
c
c************************************************************ c
c	initialize the seed value, set the random seed value
c
seed = int(secnds(0.0))
call srand(seed)
c************************************************************ C
c	begin the loop
c
do 300 ipoint=1,maxpnts
c
c	loop over shells
c
200 continue
c
c	generate a random point within the boundaries of the box
c
C
p1 = rand() * (rand()-half) * xlen * ang2au p2 = rand() * (rand()-half) * ylen * ang2au p3 = rand() * (rand()-half) * zlen * ang2au C
c
c	is this point within the outer surface of any of the atoms?
c
iflag = 0
do 100 i=1,natoms
vrad = vdwout(i)
dist = (p1 - c(1,i))**2 + (p2 - c(2,i))**2 + (p3 - c(3,i))**2 dist = dsqrt(dist)
if (dist .lt. vrad) iflag = 1
100 continue
if (iflag.eq.0) goto 200
cc
c
c	is this point within the inner surface, ie van der Waals sphere?
c
do 110 i=1,natoms
vrad = vdw(i)
dist = (p1 - c(1,i))**2 + (p2 - c(2,i))**2 + (p3 - c(3,i))**2 dist = dsqrt(dist)
if (dist .lt. vrad) goto 200
110 continue
c
c	store points (in atomic units)
c
p(1,ipoint) = p1
p(2,ipoint) = p2
p(3,ipoint) = p3
cc	if (ipo(2) .eq. 1)
cc	$ write(iout,*) 'point ',ipoint,' x,y,z ',p1,p2,p3
if (ipo(2) .eq. 1)
$ write(iout,*) p1,p2,p3
c
300 continue
return
end