C INTERNATIONAL AVS CENTER C (This disclaimer must remain at the top of all files) C C WARRANTY DISCLAIMER C C This module and the files associated with it are distributed free of charge. C It is placed in the public domain and permission is granted for anyone to use, C duplicate, modify, and redistribute it unless otherwise noted. Some modules C may be copyrighted. You agree to abide by the conditions also included in C the AVS Licensing Agreement, version 1.0, located in the main module C directory located at the International AVS Center ftp site and to include C the AVS Licensing Agreement when you distribute any files downloaded from C that site. C C The International AVS Center, MCNC, the AVS Consortium and the individual C submitting the module and files associated with said module provide absolutely C NO WARRANTY OF ANY KIND with respect to this software. The entire risk as to C the quality and performance of this software is with the user. IN NO EVENT C WILL The International AVS Center, MCNC, the AVS Consortium and the individual C submitting the module and files associated with said module BE LIABLE TO C ANYONE FOR ANY DAMAGES ARISING FROM THE USE OF THIS SOFTWARE, INCLUDING, C WITHOUT LIMITATION, DAMAGES RESULTING FROM LOST DATA OR LOST PROFITS, OR ANY C SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES. C C This AVS module and associated files are public domain software unless C otherwise noted. Permission is hereby granted to do whatever you like with C it, subject to the conditions that may exist in copyrighted materials. Should C you wish to make a contribution toward the improvement, modification, or C general performance of this module, please send us your comments: why you C liked or disliked it, how you use it, and most important, how it helps your C work. We will receive your comments at avs@ncsc.org. C C Please send AVS module bug reports to avs@ncsc.org. C subroutine espwrapper(natm,norbs,nelecs,numgridptsx,numgridptsy, & numgridptsz,nat,xyz,gridptx,gridpty,gridptz,vec,ezs,ezp,ezd, & maxo,maxatm,maxgrid,elecsp,pt_q_espflag,nlast1,nfirst1) c c Modified January 14, 1993 to use the GRAPH input to generate c the electrostatic potential over a set of grid points. c Lee Bartolotti c c Electrostatic potential code taken from MOPAC 6.0, authored c by Dr. James J. P. Stewart. c c Notice of Public Domain nature of MOPAC c c 'This computer program is a work of the United States c Government and as such is not subject to protection by c copyright (17 U.S.C. # 105.) Any person who fraudulently c places a copyright notice or does any other act contrary c to the provisions of 17 U.S. Code 506(c) shall be subject c to the penalties provided therein. This notice shall not c be altered or removed from this software and is to be on c all reproductions.' c c implicit real*8 (a-h,o-z) include 'SIZES' COMMON /POTESP/ NESP COMMON /ABC/ CO(3,NUMATM),ian(NUMATM),NATOM COMMON /WORK1/ POTPT(3,MESP1), ES(MESP), ESP(MESP) COMMON /VECTOR/ C(MORB2*2+MAXORB*2) COMMON /MOLKST/ NCLOSE,NOPEN,FRACT,NLAST(NUMATM),NFIRST(NUMATM) COMMON /EXPONT/ ZS(107),ZP(107),ZD(107) integer pt_q_espflag dimension nat(natm), xyz(3,maxatm), 1 gridptx(numgridptsx), gridpty(numgridptsy), 2 gridptz(numgridptsz), elecsp(maxgrid), 3 vec(maxo,maxo), ezs(natm), ezp(natm), ezd(natm), 4 nlast1(natm),nfirst1(natm) norb2 = norbs*norbs natom = natm nclose = nelecs/2 nopen = (nelecs + 1)/2 - nclose fract = nopen nopen = nclose + nopen nespx = numgridptsx nespy = numgridptsy nespz = numgridptsz do i=1,natom ian(i) = nat(i) co(1,i) = xyz(1,i) co(2,i) = xyz(2,i) co(3,i) = xyz(3,i) zs(ian(i)) = ezs(i) zp(ian(i)) = ezp(i) zd(ian(i)) = ezd(i) nlast(i) = nlast1(i) nfirst(i) = nfirst1(i) end do k=0 do ix=1,nespx do iy=1,nespy do iz=1,nespz k=k+1 potpt(1,k) = gridptx(ix) potpt(2,k) = gridpty(iy) potpt(3,k) = gridptz(iz) end do end do end do nesp = k k=0 do i=1,norbs do j=1,norbs k=k+1 c(k) = vec(i,j) end do end do if(pt_q_espflag.ne.0) then call pt_q_esp(norbs) else call elesp endif do i=1,nesp elecsp(i) = esp(i) end do return end subroutine pt_q_esp(nc) implicit real*8 (a-h,o-z) include 'SIZES' COMMON /POTESP/ NESP COMMON /WORK1/ POTPT(3,MESP1), ES(MESP), ESP(MESP) COMMON /VECTOR/ C(MORB2*2+MAXORB*2) COMMON /ABC/ CO(3,NUMATM),ian(NUMATM),NATOM COMMON /MOLKST/ NCLOSE,NOPEN,FRACT,NLAST(NUMATM),NFIRST(NUMATM) COMMON /CORE / TORE(107) dimension vecs(MORB2*2+MAXORB*2), q(numatm), p(MORB2*2+MAXORB*2) dimension f(MORB2*2+MAXORB*2) data zero/0.0d00/, con/1.0e-10/ call densit(c,nc,nc,nclose,nopen,fract,p,1) k=0 do i=1,natom sum = zero do j=nfirst(i),nlast(i) k = k + j sum = sum + p(k) end do nan = ian(i) q(i) = tore(nan) - sum end do do i=1,nesp rij = sqrt((co(1,1) - potpt(1,i))**2 + (co(2,1) - potpt(2,i))**2 1 +(co(3,1) - potpt(3,i))**2 + con) esp(i) = q(1)/rij end do do j=2,natom do i=1,nesp rij = sqrt((co(1,j) - potpt(1,i))**2 + (co(2,j) - potpt(2,i))**2 1 +(co(3,j) - potpt(3,i))**2 + con) esp(i) = esp(i) + q(j)/rij end do end do return end BLOCK DATA ESPBLO IMPLICIT DOUBLE PRECISION (A-H, O-Z) INCLUDE 'SIZES' COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR), 1 EX(MAXPR),ESPI(MAXORB,MAXORB), 2 FV(0:8,821),FAC(0:7), 3 DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR), 4 OVL(MAXORB,MAXORB),FC(MAXPR*6) COMMON /CORE / TORE(107) DATA TF/33.D0,37.D0,41.D0/ DATA FAC/1.D0,1.D0,2.D0,6.D0,24.D0,120.D0,720.D0,5040.D0/ DATA TORE/1.D0,0.D0, 1 1.D0,2.D0,3.D0,4.D0,5.D0,6.D0,7.D0,0.D0, 2 1.D0,2.D0,3.D0,4.D0,5.D0,6.D0,7.D0,0.D0, 3 1.D0,2.D0,3.D0,4.D0,5.D0,6.D0,7.D0,8.D0,9.D0,10.D0,11.D0,2.D0, 4 3.D0,4.D0,5.D0,6.D0,7.D0,0.D0, 5 1.D0,2.D0,3.D0,4.D0,5.D0,6.D0,7.D0,8.D0,9.D0,10.D0,11.D0,2.D0, 6 3.D0,4.D0,5.D0,6.D0,7.D0,0.D0, 7 1.D0,2.D0,3.D0,4.D0,5.D0,6.D0,7.D0,8.D0,9.D0,10.D0, 8 11.D0,12.D0,13.D0,14.D0,15.D0,16.D0, 9 3.D0,4.D0,5.D0,6.D0,7.D0,8.D0,9.D0,10.D0,11.D0,2.D0, 1 3.D0,4.D0,5.D0,6.D0,7.D0,0.D0, 2 15*0.D0,1.D0,2.D0,1.D0,-2.D0,-1.D0,0.D0/ END SUBROUTINE ELESP IMPLICIT DOUBLE PRECISION (A-H,O-Z) C*********************************************************************** C ELESP LOADS THE STO-6G BASIS SET ONTO THE ATOMS, PERFOMS THE C DEORTHOGONALIZATION OF THE COEFFICIENTS AND EVALUATES THE C ELECTRONIC CONTRIBUTION TO THE ESP. IT WAS WRITTEN BY B.H.BESLER C AND K.M.MERZ IN FEB. 1989 AT UCSF. C C*********************************************************************** DOUBLE PRECISION NORM,OVL INCLUDE 'SIZES' COMMON/ESPF/ B(NUMATM), CESPM(MAXORB,MAXORB) COMMON /POTESP/ NESP COMMON /ABC/ CO(3,NUMATM),IAN(NUMATM),NATOM COMMON /WORK1/ POTPT(3,MESP1), ES(MESP), ESP(MESP) COMMON /STO6G/ ALLC(6,5,2),ALLZ(6,5,2) COMMON /VECTOR/ C(MORB2*2+MAXORB*2) COMMON /MOLKST/ NCLOSE,NOPEN,FRACT,NLAST(NUMATM),NFIRST(NUMATM) COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR), 1 EX(MAXPR),ESPI(MAXORB,MAXORB), 2 FV(0:8,821),FAC(0:7), 3 DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR), 4 OVL(MAXORB,MAXORB),FC(MAXPR*6) 6 /CORE / TORE(107) 7 /EXPONT/ ZS(107),ZP(107),ZD(107) * * END OF MINDO/3 COMMON BLOCKS * COMMON /INDX/ INDC(MAXORB) DIMENSION CESPM2(MAXORB,MAXORB) DIMENSION CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB) PI=4.D0*ATAN(1.D0) C C PUT STO-6G BASIS SET ON ATOM CENTERS C DO 10 I=-1,10 DEX(I)=DEX2(I) 10 CONTINUE DO 20 I=0,7 FAC(I)=1.D0/FAC(I) 20 CONTINUE DO 30 M=0,8 K=1 FV(M,1)=1.D0/(2.D0*M+1.D0) DO 30 T=0.05D0,41.D0,0.05D0 K=K+1 CALL FSUB(M,T,FVAL) FV(M,K)=FVAL 30 CONTINUE C C LOAD BASIS FUNCTIONS INTO ARRAYS C ICD = 6 CALL SETUPG NC=0 NPR=0 DO 80 I=1,NATOM IF (IAN(I) .LE. 2) THEN DO 40 J=1,ICD CC(NPR+J)=ALLC(J,1,1) EX(NPR+J)=ALLZ(J,1,1)*ZS(1)**2 CEN(NPR+J,1)=CO(1,I) CEN(NPR+J,2)=CO(2,I) CEN(NPR+J,3)=CO(3,I) IAM(NPR+J,1)=0 IAM(NPR+J,2)=0 FC(NPR+J)=I 40 CONTINUE NC=NC+1 NPR=NPR+ICD ELSE C DETERMINE PRINCIPAL QUANTUM NUMBER(NQN) C OF ORBITALS TO BE USED C NQN=2 IF(IAN(I) .GT. 10 .AND. IAN(I) .LE. 18) NQN=3 IF(IAN(I) .GT. 18 .AND. IAN(I) .LE. 36) NQN=4 IF(IAN(I) .GT. 36 .AND. IAN(I) .LE. 54) NQN=5 C DO 50 J=1,ICD CC(NPR+J)=ALLC(J,NQN,1) EX(NPR+J)=ALLZ(J,NQN,1)*ZS(IAN(I))**2 CEN(NPR+J,1)=CO(1,I) CEN(NPR+J,2)=CO(2,I) CEN(NPR+J,3)=CO(3,I) IAM(NPR+J,1)=0 IAM(NPR+J,2)=0 50 CONTINUE NC=NC+1 NPR=NPR+ICD DO 70 K=1,3 DO 60 J=1,ICD CC(NPR+J)=ALLC(J,NQN,2) EX(NPR+J)=ALLZ(J,NQN,2)*ZP(IAN(I))**2 CEN(NPR+J,1)=CO(1,I) CEN(NPR+J,2)=CO(2,I) CEN(NPR+J,3)=CO(3,I) IAM(NPR+J,1)=1 IAM(NPR+J,2)=K 60 CONTINUE NC=NC+1 NPR=NPR+ICD 70 CONTINUE ENDIF 80 CONTINUE C C CALCULATE NORMALIZATION CONSTANTS AND INCLUDE C THEM IN THE CONTRACTION COEFFICIENTS C DO 90 I=1,NPR NORM=(2.D0*EX(I)/PI)**0.75D0*(4.D0*EX(I))**(IAM(I,1)/2.D0)/ 1 SQRT(DEX(2*IAM(I,1)-1)) CC(I)=CC(I)*NORM 90 CONTINUE IPR=0 C C PERFORM SORT OF PRIMITIVES BY ANGULAR MOMENTUM C IS=0 IP=0 IPC=0 ISC=0 J=0 DO 100 I=1,NPR IF (IAM(I,1) .EQ. 0) THEN IS=IS+1 IND(IS)=I ENDIF 100 CONTINUE IP=IS DO 110 I=1,NPR IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 1) THEN IP=IP+1 IND(IP)=I ENDIF 110 CONTINUE DO 120 I=1,NPR IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 2) THEN IP=IP+1 IND(IP)=I ENDIF 120 CONTINUE DO 130 I=1,NPR IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 3) THEN IP=IP+1 IND(IP)=I ENDIF 130 CONTINUE DO 140 I=1,NC IN=I*ICD-ICD+1 IF (IAM(IN,1) .EQ. 0) THEN ISC=ISC+1 INDC(ISC)=I ENDIF 140 CONTINUE IPC=ISC DO 150 I=1,NC IN=I*ICD-ICD+1 IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 1) THEN IPC=IPC+1 INDC(IPC)=I ENDIF 150 CONTINUE DO 160 I=1,NC IN=I*ICD-ICD+1 IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 2) THEN IPC=IPC+1 INDC(IPC)=I ENDIF 160 CONTINUE DO 170 I=1,NC IN=I*ICD-ICD+1 IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 3) THEN IPC=IPC+1 INDC(IPC)=I ENDIF 170 CONTINUE DO 180 I=1,NPR TEMP(I)=CC(IND(I)) 180 CONTINUE DO 190 I=1,NPR CC(I)=TEMP(I) 190 CONTINUE DO 200 I=1,NPR TEMP(I)=EX(IND(I)) 200 CONTINUE DO 210 I=1,NPR EX(I)=TEMP(I) 210 CONTINUE DO 220 I=1,NPR TEMP(I)=CEN(IND(I),1) 220 CONTINUE DO 230 I=1,NPR CEN(I,1)=TEMP(I) 230 CONTINUE DO 240 I=1,NPR TEMP(I)=CEN(IND(I),2) 240 CONTINUE DO 250 I=1,NPR CEN(I,2)=TEMP(I) 250 CONTINUE DO 260 I=1,NPR TEMP(I)=CEN(IND(I),3) 260 CONTINUE DO 270 I=1,NPR CEN(I,3)=TEMP(I) 270 CONTINUE DO 280 I=1,NPR ITEMP(I)=IAM(IND(I),1) 280 CONTINUE DO 290 I=1,NPR IAM(I,1)=ITEMP(I) 290 CONTINUE DO 300 I=1,NPR ITEMP(I)=IAM(IND(I),2) 300 CONTINUE DO 310 I=1,NPR IAM(I,2)=ITEMP(I) 310 CONTINUE C CALCULATE OVERLAP MATRIX OF STO-6G FUNCTIONS C DO 320 J=1,NC CALL OVLP(J,1,IS,IP,NPR,NC,ICD) 320 CONTINUE C DO 330 J=1,NC DO 330 K=1,NC CESPM2(INDC(J),INDC(K))=OVL(J,K) 330 CONTINUE DO 340 J=1,NC DO 340 K=1,NC OVL(J,K)=CESPM2(J,K) 340 CONTINUE L=0 DO 350 I=1,NC DO 350 J=1,I L=L+1 CESP(L)=OVL(I,J) 350 CONTINUE C C DEORTHOGONALIZE THE COEFFICIENTS AND REFORM THE DENSITY MATRIX C CALL RSP(CESP,NC,1,TEMP,CESPML) DO 360 I=1,NC DO 360 J=1,I SUM=0.D0 DO 360 K=1,NC SUM=SUM+CESPML(I+(K-1)*NC)/SQRT(TEMP(K))*CESPML(J+(K-1)*N 1C) CESP(I+(J-1)*NC)=SUM CESP(J+(I-1)*NC)=SUM 360 CONTINUE CALL MULT(C,CESP,CESPML,NC) CALL DENSIT(CESPML,NC,NC,NCLOSE,NOPEN,FRACT,CESP,2) C C NOW CALCULATE THE ELECTRONIC CONTRIBUTION TO THE ELECTROSTATIC POT C L=0 DO 370 I=1,NC DO 370 J=1,I L=L+1 CESPM(I,J)=CESP(L) CESPM(J,I)=CESP(L) 370 CONTINUE IPX=(NPR-IS)/3 IPE=IS+IPX DO 380 I=1,NESP ES(I)=0.D0 380 CONTINUE CALL NAICAS(ISC,IS,IP,NPR,NC,IPE,IPX,ICD) CALL NAICAP(ISC,IS,IP,NPR,NC,IPE,IPX,ICD) C CALCULATE TOTAL ESP AND FORM ARRAYS FOR ESPFIT DO 400 I=1,NESP ESP(I)=0.D0 DO 390 J=1,NATOM RA=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2+(CO( 13,J)-POTPT(3,I))**2) ESP(I)=ESP(I)+TORE(IAN(J))/(RA) 390 CONTINUE ESP(I)=ESP(I)-ES(I) DO 400 J=1,NATOM RIJ=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2 1+(CO(3,J)-POTPT(3,I))**2) B(J)=B(J)+ESP(I)*1.D0/RIJ 400 CONTINUE C C IF REQUESTED WRITE OUT ELECTRIC POTENTIAL DATA TO C UNIT 21 C OPEN(21,STATUS='unknown') rewind 21 WRITE(21,'(I5)') NESP C *** THE FOLLOWING WRITE STATEMENT MODIFIED AS NOTED IN C ON THE CHEMISTRY MAILING LIST. C DO 410 I=1,NESP 410 WRITE(21,420) I, ESP(I),POTPT(1,I),POTPT(2,I),POTPT(3,I) 420 FORMAT(i4,1X,4E16.7) close(21) RETURN END DOUBLE PRECISION FUNCTION DEX2(M) IMPLICIT DOUBLE PRECISION (A-H,O-Z) IF(M .LT. 2) THEN DEX2=1 ELSE DEX2=1 DO 10 I=1,M,2 10 DEX2=DEX2*I ENDIF RETURN END C*********************************************************************** SUBROUTINE FSUB(N,X,FVAL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C*********************************************************************** C C CALCULATE THE FM(T). KINDLY SUPPLIED BY RUS PITZER AND CLEANED UP C BY K.M.MERZ IN FEB. 1989 AT UCSF. C C ON INPUT: N = INDEX C X = EXPONENT C ON OUTPUT: FVAL = VALUE OF THE FUNCTION C C FOR MORE DETAILS SEE: OBARA AND SAIKA J. CHEM. PHYS. 1986,84,3963 C*********************************************************************** DIMENSION FF(21),TERM(200) DATA A0, A1S2, PIE4, A1 1 /0.0D0,0.5D0,0.7853981633974483096156608D0,1.0D0/ DATA XSW /24.0D0/ E=A1S2*EXP(-X) FAC0=N FAC0=FAC0+A1S2 IF(X.GT.XSW) GO TO 50 C C USE POWER SERIES C 10 FAC=FAC0 TERM0=E/FAC SUM=TERM0 KU=(X-FAC0) IF(KU.LT.1) GO TO 30 C C SUM INCREASING TERMS FORWARDS C DO 20 K=1,KU FAC=FAC+A1 TERM0=TERM0*X/FAC SUM=SUM+TERM0 20 CONTINUE 30 I=1 FAC=FAC+A1 TERM(1)=TERM0*X/FAC SUMA=SUM+TERM(1) IF(SUM.EQ.SUMA) GO TO 90 40 I=I+1 FAC=FAC+A1 TERM(I)=TERM(I-1)*X/FAC SUM1=SUMA SUMA=SUMA+TERM(I) IF(SUM1-SUMA) 40,90,40 C C USE ASYMPTOTIC SERIES C 50 SUM=SQRT(PIE4/X) IF(N.EQ.0) GO TO 70 FAC=-A1S2 DO 60 K=1,N FAC=FAC+A1 SUM=SUM*FAC/X 60 CONTINUE 70 I=1 TERM(1)=-E/X SUMA=SUM+TERM(1) IF(SUM.EQ.SUMA) GO TO 90 FAC=FAC0 KU=(X+FAC0-A1) DO 80 I=2,KU FAC=FAC-A1 TERM(I)=TERM(I-1)*FAC/X SUM1=SUMA SUMA=SUMA+TERM(I) IF(SUM1.EQ.SUMA) GO TO 90 80 CONTINUE C C XSW SET TOO LOW. USE POWER SERIES. C GO TO 10 C C SUM DECREASING TERMS BACKWARDS C 90 SUM1=A0 DO 100 K=1,I SUM1=SUM1+TERM(I+1-K) 100 CONTINUE FF(N+1)=SUM+SUM1 C C USE RECURRENCE RELATION C IF(N.EQ.0) GOTO 120 DO 110 K=1,N FAC0=FAC0-A1 FF(N+1-K)=(E+X*FF(N+2-K))/FAC0 110 CONTINUE 120 FVAL=FF(N+1) RETURN END SUBROUTINE SETUPG IMPLICIT DOUBLE PRECISION (A-H,O-Z) INCLUDE 'SIZES' COMMON /STO6G/ ALLC(6,5,2),ALLZ(6,5,2) C SET-UP THE STEWART'S STO-6G EXPANSIONS C 1S ALLZ(1,1,1) =2.310303149D01 ALLZ(2,1,1) =4.235915534D00 ALLZ(3,1,1) =1.185056519D00 ALLZ(4,1,1) =4.070988982D-01 ALLZ(5,1,1) =1.580884151D-01 ALLZ(6,1,1) =6.510953954D-02 C ALLC(1,1,1) =9.163596280D-03 ALLC(2,1,1) =4.936149294D-02 ALLC(3,1,1) =1.685383049D-01 ALLC(4,1,1) =3.705627997D-01 ALLC(5,1,1) =4.164915298D-01 ALLC(6,1,1) =1.303340841D-01 C 2S ALLZ(1,2,1) =2.768496241D01 ALLZ(2,2,1) =5.077140627D00 ALLZ(3,2,1) =1.426786050D00 ALLZ(4,2,1) =2.040335729D-01 ALLZ(5,2,1) =9.260298399D-02 ALLZ(6,2,1) =4.416183978D-02 C ALLC(1,2,1) =-4.151277819D-03 ALLC(2,2,1) =-2.067024148D-02 ALLC(3,2,1) =-5.150303337D-02 ALLC(4,2,1) =3.346271174D-01 ALLC(5,2,1) =5.621061301D-01 ALLC(6,2,1) =1.712994697D-01 C 2P ALLZ(1,2,2) =5.868285913D00 ALLZ(2,2,2) =1.530329631D00 ALLZ(3,2,2) =5.475665231D-01 ALLZ(4,2,2) =2.288932733D-01 ALLZ(5,2,2) =1.046655969D-01 ALLZ(6,2,2) =4.948220127D-02 C ALLC(1,2,2) =7.924233646D-03 ALLC(2,2,2) =5.144104825D-02 ALLC(3,2,2) =1.898400060D-01 ALLC(4,2,2) =4.049863191D-01 ALLC(5,2,2) =4.012362861D-01 ALLC(6,2,2) =1.051855189D-01 C 3S ALLZ(1,3,1) =3.273031938D00 ALLZ(2,3,1) =9.200611311D-01 ALLZ(3,3,1) =3.593349765D-01 ALLZ(4,3,1) =8.636686991D-02 ALLZ(5,3,1) =4.797373812D-02 ALLZ(6,3,1) =2.724741144D-02 ALLC(1,3,1) =-6.775596947D-03 ALLC(2,3,1) =-5.639325779D-02 ALLC(3,3,1) =-1.587856086D-01 ALLC(4,3,1) =5.534527651D-01 ALLC(5,3,1) =5.015351020D-01 ALLC(6,3,1) =7.223633674D-02 C 3P ALLZ(1,3,2) =5.077973607D00 ALLZ(2,3,2) =1.340786940D00 ALLZ(3,3,2) =2.248434849D-01 ALLZ(4,3,2) =1.131741848D-01 ALLZ(5,3,2) =6.076408893D-02 ALLZ(6,3,2) =3.315424265D-02 ALLC(1,3,2) =-3.329929840D-03 ALLC(2,3,2) =-1.419488340D-02 ALLC(3,3,2) =1.639395770D-01 ALLC(4,3,2) =4.485358256D-01 ALLC(5,3,2) =3.908813050D-01 ALLC(6,3,2) =7.411456232D-02 C 4S ALLZ(1,4,1) = 1.365346 D+00 ALLZ(2,4,1) = 4.393213 D-01 ALLZ(3,4,1) = 1.877069 D-01 ALLZ(4,4,1) = 9.360270 D-02 ALLZ(5,4,1) = 5.052263 D-02 ALLZ(6,4,1) = 2.809354 D-02 ALLC(1,4,1) = 3.775056 D-03 ALLC(2,4,1) =-5.585965 D-02 ALLC(3,4,1) =-3.192946 D-01 ALLC(4,4,1) =-2.764780 D-02 ALLC(5,4,1) = 9.049199 D-01 ALLC(6,4,1) = 3.406258 D-01 C 4P ALLC(1,4,2) =-7.052075 D-03 ALLC(2,4,2) =-5.259505 D-02 ALLC(3,4,2) =-3.773450 D-02 ALLC(4,4,2) = 3.874773 D-01 ALLC(5,4,2) = 5.791672 D-01 ALLC(6,4,2) = 1.221817 D-01 ALLZ(1,4,2) = 1.365346 D+00 ALLZ(2,4,2) = 4.393213 D-01 ALLZ(3,4,2) = 1.877069 D-01 ALLZ(4,4,2) = 9.360270 D-02 ALLZ(5,4,2) = 5.052263 D-02 ALLZ(6,4,2) = 2.809354 D-02 C 5S ALLZ(1,5,1) = 7.701420258D-01 ALLZ(2,5,1) = 2.756268915D-01 ALLZ(3,5,1) = 1.301847480D-01 ALLZ(4,5,1) = 6.953441940D-02 ALLZ(5,5,1) = 4.002545502D-02 ALLZ(6,5,1) = 2.348388309D-02 ALLC(1,5,1) = 1.267447151D-02 ALLC(2,5,1) = 3.266734789D-03 ALLC(3,5,1) =-4.307553999D-01 ALLC(4,5,1) =-3.231998963D-01 ALLC(5,5,1) = 1.104322879D+00 ALLC(6,5,1) = 4.368498703D-01 C 5P ALLZ(1,5,2) = 7.701420258D-01 ALLZ(2,5,2) = 2.756268915D-01 ALLZ(3,5,2) = 1.301847480D-01 ALLZ(4,5,2) = 6.953441940D-02 ALLZ(5,5,2) = 4.002545502D-02 ALLZ(6,5,2) = 2.348388309D-02 ALLC(1,5,2) =-1.105673292D-03 ALLC(2,5,2) =-6.243132446D-02 ALLC(3,5,2) =-1.628476766D-01 ALLC(4,5,2) = 3.210328714D-01 ALLC(5,5,2) = 6.964579592D-01 ALLC(6,5,2) = 1.493146125D-01 RETURN END SUBROUTINE OVLP(IC,IESP,IS,IP,NPR,NC,ICD) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C*********************************************************************** C C OVLP CALCULATES THE OVERLAP INTEGRALS FOR A STO-6G BASIS SET. C THE RESULTING INTEGRALS ARE USED IN THE DEORTHOGONALIZATION C PROCESS. C THE CODE WAS WRITTEN BY B.H.BESLER AND K.M.MERZ IN FEB. 1989 C AT UCSF. C C ON INPUT: IC = LOOP INDEX C IESP = LOOP INDEX C IS = NUMBER OF S ORBITALS C IP = NUMBER OF P ORBITALS C NPR = NUMBER OF PRIMITIVES C NC = NUMBER OF CONTRACTED FUNCTIONS C C ON OUTPUT: OVL IS FILLED WITH THE OVERLAP INTEGRAL VALUE C C FOR FURTHER INFO SEE: OBARA & SAIKA J.CHEM.PHYS. 1986,84,3963 C*********************************************************************** DOUBLE PRECISION NAI,NAI1,NAI2 INCLUDE 'SIZES' COMMON /POTESP/ NESP COMMON /ABC/ CO(3,NUMATM),IAN(NUMATM),NATOM COMMON /WORK1/ POTPT(3,MESP1), ES(MESP), ESP(MESP) COMMON /EXPONT/ ZS(107),ZP(107),ZD(107) COMMON /STO6G/ ALLC(6,5,2),ALLZ(6,5,2) COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR), 1EX(MAXPR),ESPI(MAXORB,MAXORB),FV(0:8,821), 2FAC(0:7),DEX(-1:96),TF(0:2), 3TEMP(MAXPR),ITEMP(MAXPR),OVL(MAXORB,MAXORB),XDMY(MAXPR*6) COMMON/X/ DX(MAXPR),DY(MAXPR),DZ(MAXPR),F1(MAXPR,6),F2(MAXPR,6), 1TD(MAXPR),CE(MAXPR,6),U(MAXPR,6),EXS(MAXPR,6),EXPN(MAXPR,6), 2NAI(MAXPR,6),EWCX(MAXPR,6),EWCY(MAXPR,6),EWCZ(MAXPR,6),F0(MAXPR,6) 3,NAI1(MAXPR,6),NAI2(MAXPR,6) C C CALCULATE DISTANCE ARRAYS C PI=4.D0*ATAN(1.D0) IPR=IC*ICD-ICD+1 ISTART=IPR DO 10 I=ISTART,NPR DX(I)=CEN(IPR,1)-CEN(I,1) DY(I)=CEN(IPR,2)-CEN(I,2) DZ(I)=CEN(IPR,3)-CEN(I,3) TD(I)=DX(I)**2+DY(I)**2+DZ(I)**2 10 CONTINUE C C CALCULATE EXPONENT SUM C DO 20 I=ISTART,NPR DO 20 J=1,ICD EXS(I,J)=1.D0/(EX(IPR+J-1)+EX(I)) CE(I,J)=EX(IPR+J-1)*EX(I)*EXS(I,J) 20 CONTINUE C C CALCULATE EXPONENT WEIGHTED CENTERS C DO 30 I=ISTART,NPR DO 30 J=1,ICD EWCX(I,J)=(EX(I)*CEN(I,1)+EX(IPR+J-1) 1*CEN(IPR+J-1,1))*EXS(I,J) EWCY(I,J)=(EX(I)*CEN(I,2)+EX(IPR+J-1) 1*CEN(IPR+J-1,2))*EXS(I,J) EWCZ(I,J)=(EX(I)*CEN(I,3)+EX(IPR+J-1) 1*CEN(IPR+J-1,3))*EXS(I,J) 30 CONTINUE DO 40 I=1,NPR DO 40 J=1,ICD EXPN(I,J)=EXP(-CE(I,J)*TD(I)) NAI(I,J)=(PI*EXS(I,J))**1.5D0*EXPN(I,J) EXPN(I,J)=NAI(I,J) 40 CONTINUE C C CALCULATE (S||P) ESP INTEGRALS C IF((IAM(IPR,1) .EQ. 0) .AND. (IS .NE. IP)) THEN NP=IS+1 DO 80 I=NP,NPR DO 80 J=1,ICD GO TO (50,60,70),IAM(I,2) 50 NAI(I,J)=(EWCX(I,J)-CEN(I,1))*EXPN(I,J) go TO 80 60 NAI(I,J)=(EWCY(I,J)-CEN(I,2))*EXPN(I,J) GO TO 80 70 NAI(I,J)=(EWCZ(I,J)-CEN(I,3))*EXPN(I,J) 80 CONTINUE ENDIF C C CALCULATE (P||S) ESP INTEGRALS C IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN NP=IS+1 DO 120 I= ISTART,NPR DO 120 J=1,ICD GO TO (90,100,110),IAM(IPR+J-1,2) 90 NAI(I,J)=(EWCX(I,J)-CEN(IPR+J-1,1))*EXPN(I,J) GO TO 120 100 NAI(I,J)=(EWCY(I,J)-CEN(IPR+J-1,2))*EXPN(I,J) GO TO 120 110 NAI(I,J)=(EWCZ(I,J)-CEN(IPR+J-1,3))*EXPN(I,J) 120 CONTINUE ENDIF C C CALCULATE (P||P) ESP INTEGRALS C IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN DO 160 I=ISTART,NPR DO 160 J=1,ICD GO TO (130,140,150),IAM(I,2) 130 NAI(I,J)=(EWCX(I,J)-CEN(I,1))*NAI(I,J) IF(IAM(IPR+J-1,2) .EQ. IAM(I,2)) 1NAI(I,J)=NAI(I,J)+EXS(I,J)*0.5D0 2 *EXPN(I,J) GO TO 160 140 NAI(I,J)=(EWCY(I,J)-CEN(I,2))*NAI(I,J) IF(IAM(IPR+J-1,2) .EQ. IAM(I,2)) 1NAI(I,J)=NAI(I,J)+EXS(I,J)*0.5D0 2 *EXPN(I,J) GO TO 160 150 NAI(I,J)=(EWCZ(I,J)-CEN(I,3))*NAI(I,J) IF(IAM(IPR+J-1,2) .EQ. IAM(I,2)) 1NAI(I,J)=NAI(I,J)+EXS(I,J)*0.5D0 2 *EXPN(I,J) 160 CONTINUE ENDIF IPS=IC*ICD-ICD+1 DO 180 I=IC,NC JPS=I*ICD-ICD+1 OVL(IC,I)=0.D0 DO 170 J=JPS,JPS+ICD-1 DO 170 K=IPS,IPS+ICD-1 OVL(IC,I)=OVL(IC,I)+CC(J)*CC(K)*NAI(J,K-IPS+1) 170 CONTINUE OVL(I,IC)=OVL(IC,I) 180 CONTINUE RETURN END SUBROUTINE MULT(C,S,VECS,N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION C(N,*), S(N,*), VECS(N,*) *********************************************************************** * * MULT IS USED IN THE MULLIKEN ANALYSIS ONLY. IT PERFORMS THE * OPERATION:- * VECS=BACK-TRANSFORMED EIGENVECTORS * VECS = C*S C =UN-BACK-TRANSFORMED VECTORS * S =1/SQRT(OVERLAP MATRIX) * *********************************************************************** DO 20 I=1,N DO 20 J=1,N SUM=0.D0 DO 10 K=1,N 10 SUM=SUM+C(K,I)*S(J,K) 20 VECS(J,I)=SUM RETURN END SUBROUTINE DENSIT( C,MDIM, NORBS,NDUBL, NSINGL, FRACT, P,MODE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION P(*), C(MDIM,*) C*********************************************************************** C C DENSIT COMPUTES THE DENSITY MATRIX GIVEN THE EIGENVECTOR MATRIX, AND C INFORMATION ABOUT THE M.O. OCCUPANCY. C C INPUT: C = SQUARE EIGENVECTOR MATRIX, C IS OF SIZE MDIM BY MDIM C AND THE EIGENVECTORS ARE STORED IN THE TOP LEFT-HAND C CORNER. C NORBS = NUMBER OF ORBITALS C NDUBL = NUMBER OF DOUBLY-OCCUPIED M.O.S ( =0 IF UHF) C NSINGL= NUMBER OF SINGLY OR FRACTIONALLY OCCUPIED M.O.S. C MODE = 2 IF POSITRON EQUIVALENT IS NOT TO BE USED C C ON EXIT: P = DENSITY MATRIX C C*********************************************************************** C C SET UP LIMITS FOR SUMS C NL1 = BEGINING OF ONE ELECTRON SUM C NU1 = END OF SAME C NL2 = BEGINING OF TWO ELECTRON SUM C NU2 = END OF SAME C NORBS2=NORBS/2 NSINGL=MAX(NDUBL,NSINGL) IF(NDUBL.NE.0.AND.NSINGL .GT. NORBS2 .AND. MODE.NE.2) THEN C C TAKE POSITRON EQUIVALENT C SIGN=-1.D0 FRAC=2.D0-FRACT CONST=2.D0 NL2=NSINGL+1 NU2=NORBS NL1=NDUBL+1 NU1=NSINGL ELSE C C TAKE ELECTRON EQUIVALENT C SIGN=1.D0 FRAC=FRACT CONST=0.D0 NL2=1 NU2=NDUBL NL1=NDUBL+1 NU1=NSINGL ENDIF L=0 DO 40 I=1,NORBS DO 30 J=1,I L=L+1 SUM2=0.D0 SUM1=0.D0 DO 10 K=NL2,NU2 10 SUM2=SUM2+C(I,K)*C(J,K) SUM2=SUM2*2.D0 DO 20 K=NL1,NU1 20 SUM1=SUM1+C(I,K)*C(J,K) 30 P(L)=(SUM2+SUM1*FRAC)*SIGN 40 P(L)=CONST+P(L) RETURN END SUBROUTINE NAICAS(ISC,IS,IP,NPR,NC,IPE,IPX,ICD) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C*********************************************************************** C C THIS SUBROUTINE EVALUATES (S|S) , (S|P) TYPE NUCLEAR ATTRACTION C INTEGRALS FOR A STO-NG BASIS SET C WRITTEN BY B.H. BESLER AT FORD SCIENTIFIC RESEARCH LABS IN C DECEMBER 1989. C C ON INPUT: IC = LOOP INDEX OF THE GAUSSIAN C IESP = LOOP INDEX OF THE ESP POINT C IPE = INDEX OF LAST Px PRIMITIVE C IPX = NUMBER OF Px PRIMITIVES C IS = NUMBER OS S ORBITALS C ISC = NUMBER OF CONTRACTED S ORBITALS C IP = NUMBER OF P ORBITALS C NPR = NUMBER OF PRIMITIVES C NC = NUMBER OF CONTRACTED FUNCTIONS C C C FOR MORE INFO SEE: OBARA&SAIKA J.CHEM.PHYS. 1986,84,3963. C*********************************************************************** INCLUDE 'SIZES' DOUBLE PRECISION NAI,NAI1,NAI2 COMMON/ESPF/ B(NUMATM), CESPM(MAXORB,MAXORB) COMMON /INDX/ INDC(MAXORB) COMMON /POTESP/ NESP COMMON /ABC/ CO(3,NUMATM),IAN(NUMATM),NATOM COMMON /WORK1/ POTPT(3,MESP1), ES(MESP), ESP(MESP) COMMON /EXPONT/ ZS(107),ZP(107),ZD(107) COMMON /STO6G/ ALLC(6,5,2),ALLZ(6,5,2) COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR), 1EX(MAXPR),ESPI(MAXORB,MAXORB),FV(0:8,821), 2FAC(0:7),DEX(-1:96),TF(0:2), 3TEMP(MAXPR),ITEMP(MAXPR),OVL(MAXORB,MAXORB),EXSR(MAXPR,6) COMMON/X/ DX(MAXPR),DY(MAXPR),DZ(MAXPR),F1(MAXPR,6),F2(MAXPR,6), 1TD(MAXPR),CE(MAXPR,6),U(MAXPR,6),EXS(MAXPR,6),EXPN(MAXPR,6), 2NAI(MAXPR,6),EWCX(MAXPR,6),EWCY(MAXPR,6),EWCZ(MAXPR,6),F0(MAXPR,6) 3,NAI1(MAXPR,6),NAI2(MAXPR,6) C C CALCULATE DISTANCE ARRAYS C PI=4.D0*ATAN(1.D0) IPX2=2*IPX JSTART=1 NP=IS+1 DO 200 IC=JSTART,ISC IPR=IC*ICD-ICD+1 ISTART=IPR DO 20 I=ISTART,IPE DX(I)=CEN(IPR,1)-CEN(I,1) DY(I)=CEN(IPR,2)-CEN(I,2) DZ(I)=CEN(IPR,3)-CEN(I,3) TD(I)=DX(I)**2+DY(I)**2+DZ(I)**2 20 CONTINUE C C CALCULATE EXPONENT SUM C DO 30 I=ISTART,IPE DO 30 J=1,ICD EXSR(I,J)=EX(IPR+J-1)+EX(I) EXS(I,J)=1.D0/EXSR(I,J) CE(I,J)=EX(IPR+J-1)*EX(I)*EXS(I,J) EXPN(I,J)=EXP(-CE(I,J)*TD(I)) 30 CONTINUE C C CALCULATE EXPONENT WEIGHTED CENTERS C DO 40 I=ISTART,IPE DO 40 J=1,ICD EWCX(I,J)=(EX(I)*CEN(I,1)+EX(IPR+J-1) 1*CEN(IPR+J-1,1))*EXS(I,J) EWCY(I,J)=(EX(I)*CEN(I,2)+EX(IPR+J-1) 1*CEN(IPR+J-1,2))*EXS(I,J) EWCZ(I,J)=(EX(I)*CEN(I,3)+EX(IPR+J-1) 1*CEN(IPR+J-1,3))*EXS(I,J) 40 CONTINUE C C BEGIN LOOP OVER ESP POINTS C DO 180 IESP=1,NESP POTP1=POTPT(1,IESP) POTP2=POTPT(2,IESP) POTP3=POTPT(3,IESP) C C BEGIN LOOP OVER COMPONENTS OF CONTRACTED FUNCTION IC C DO 150 J=1,ICD C C CALCULATE DISTANCE BETWEEN EXPONENT WEIGHTED AND PROBE POINT C DO 50 I=ISTART,IPE U(I,J)=((EWCX(I,J)-POTP1)**2+(EWCY(I,J)-POTP2)**2+ 1 (EWCZ(I,J)-POTP3)**2)*EXSR(I,J) NAI(I,J)=SQRT(PI/U(I,J)) 50 CONTINUE C C CALCULATE ESP INTEGRALS C DO 70 I=ISTART,IPE IF(U(I,J) .LE. TF(0)) THEN IREF=DNINT(U(I,J)*20.D0) REF=0.05D0*IREF RES=U(I,J)-REF TERM=1.D0 F0(I,J)=0.D0 DO 60 K=0,6 F=FV(K,IREF+1) TS=F*TERM*FAC(K) TERM=-TERM*RES F0(I,J)=F0(I,J)+TS 60 CONTINUE ELSE F0(I,J)=NAI(I,J)*0.5D0 ENDIF 70 CONTINUE DO 90 I=NP,IPE IF(U(I,J) .LE. TF(1)) THEN IREF=DNINT(U(I,J)*20.D0) REF=0.05D0*IREF RES=U(I,J)-REF TERM1=1.D0 F1(I,J)=0.D0 DO 80 K=0,6 FI=FV(K+1,IREF+1) TS1=FI*TERM1*FAC(K) TERM1=-TERM1*RES F1(I,J)=F1(I,J)+TS1 80 CONTINUE ELSE F1(I,J)=NAI(I,J)*0.25D0/U(I,J) ENDIF 90 CONTINUE DO 100 I=ISTART,IS 100 U(I,J)=2.D0*PI*EXS(I,J)*EXPN(I,J)*F0(I,J) NP=IS+1 DO 110 I=NP,IPE NAI(I,J)=2.D0*PI*EXS(I,J)*EXPN(I,J)*F0(I,J) NAI1(I,J)=2.D0*PI*EXS(I,J)*EXPN(I,J)*F1(I,J) 110 CONTINUE C C CALCULATE (S||P) ESP INTEGRALS C IF((IAM(IPR,1) .EQ. 0) .AND. (IS .NE. IP)) THEN DO 120 I=NP,IPE 120 U(I,J)=(EWCX(I,J)-CEN(I,1))*NAI(I,J) 1-(EWCX(I,J)-POTP1)*NAI1(I,J) DO 130 I=IPE+1,IPE+1+IPX 130 U(I,J)=(EWCY(I-IPX,J)-CEN(I-IPX,2))*NAI(I-IPX,J) 1-(EWCY(I-IPX,J)-POTP2)*NAI1(I-IPX,J) DO 140 I=IPE+1+IPX,NPR 140 U(I,J)=(EWCZ(I-IPX2,J)-CEN(I-IPX2,3))*NAI(I-IPX2,J) 1-(EWCZ(I-IPX2,J)-POTP3)*NAI1(I-IPX2,J) ENDIF 150 CONTINUE IPS=IC*ICD-ICD+1 DO 170 I=IC,NC JPS=I*ICD-ICD+1 ESPI(I,IC)=0.D0 DO 160 J=JPS,JPS+ICD-1 DO 160 K=IPS,IPS+ICD-1 ESPI(I,IC)=ESPI(I,IC)+CC(J)*CC(K)*U(J,K-IPS+1) 160 CONTINUE ES(IESP)=ES(IESP)+2.D0*CESPM(INDC(I),INDC(IC))*ESPI(I,IC) 170 CONTINUE ES(IESP)=ES(IESP)-CESPM(INDC(IC),INDC(IC))*ESPI(IC,IC) 180 CONTINUE IESPS=0 C 200 CONTINUE RETURN END SUBROUTINE NAICAP(ISC,IS,IP,NPR,NC,IPE,IPX,ICD) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C*********************************************************************** C THIS ROUTINE EVALUATES (P|P) NUCLEAR ATTRACTION INTEGRALS OVER C C A STO-NG BASIS SET. C WRITTEN BY B.H. BESLER AT FORD SCIENTIFIC RESEARCH LABS IN C SEPT. 1989 C C ON INPUT: IC = LOOP INDEX OF THE GAUSSIAN C ICD = CONTRACTION DEPTH OF BASIS SET C IESP = LOOP INDEX OF THE ESP POINT C IS = NUMBER OS S PRIMITIVES C IPE = INDEX OF LAST PX PRIMITIVE C IPX = NUMBER OF PX PRIMITIVES C IS = NUMBER OS S PRIMITIVES C ISC = NUMBER OF CONTRACTED C NPR = NUMBER OF PRIMITIVES C NC = NUMBER OF CONTRACTED FUNCTIONS C C C FOR MORE INFO SEE: OBARA&SAIKA J.CHEM.PHYS. 1986,84,3963. C*********************************************************************** INCLUDE 'SIZES' DOUBLE PRECISION NAI,NAI1,NAI2 COMMON/ESPF/ B(NUMATM), CESPM(MAXORB,MAXORB) COMMON /INDX/ INDC(MAXORB) COMMON /POTESP/ NESP COMMON /ABC/ CO(3,NUMATM),IAN(NUMATM),NATOM COMMON /WORK1/ POTPT(3,MESP1), ES(MESP), ESP(MESP) COMMON /EXPONT/ ZS(107),ZP(107),ZD(107) COMMON /STO6G/ ALLC(6,5,2),ALLZ(6,5,2) COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR), 1EX(MAXPR),ESPI(MAXORB,MAXORB),FV(0:8,821), 2FAC(0:7),DEX(-1:96),TF(0:2), 3TEMP(MAXPR),ITEMP(MAXPR),OVL(MAXORB,MAXORB),EXSR(MAXPR,6) COMMON/X/ DX(MAXPR),DY(MAXPR),DZ(MAXPR),F1(MAXPR,6),F2(MAXPR,6), 1TD(MAXPR),CE(MAXPR,6),U(MAXPR,6),EXS(MAXPR,6),EXPN(MAXPR,6), 2NAI(MAXPR,6),EWCX(MAXPR,6),EWCY(MAXPR,6),EWCZ(MAXPR,6),F0(MAXPR,6) 3,NAI1(MAXPR,6),NAI2(MAXPR,6) COMMON/FP/ PF0(MAXHES),PF1(MAXHES),PF2(MAXHES),ID(MAXPAR), 1PEXS(MAXHES),PCE(MAXHES),PEXPN(MAXHES),PTD(MAXHES), 2PEWCX(MAXHES),PEWCY(MAXHES),PEWCZ(MAXHES),IRD(MAXHES) C SET NUMBER OF EQUALLY SPACED DUMPS IDN=10 C IDC=0 IPX2=2*IPX PI=4.D0*ATAN(1.D0) NP=IS+1 C SETUP INDEX ARRAY DO 10 I=NP,IPE IRD(I)=I-IS IRD(I+IPX)=I-IS IRD(I+IPX2)=I-IS 10 CONTINUE C C CALCULATE QUANTITIES INVARIANT WITH ESP POINT FOR C (P|P) ESP INTEGRALS C IL=L L=0 DO 30 I=NP,IPE DO 20 J=I,IPE L=L+1 PTD(L)=(CEN(I,1)-CEN(J,1))**2+(CEN(I,2)-CEN(J,2))**2+ 1(CEN(I,3)-CEN(J,3))**2 PEXS(L)=1.d0/(EX(I)+EX(J)) PCE(L)=EX(I)*EX(J)*PEXS(L) PEXPN(L)=EXP(-PCE(L)*PTD(L)) PEWCX(L)=(EX(I)*CEN(I,1)+EX(J)*CEN(J,1))*PEXS(L) PEWCY(L)=(EX(I)*CEN(I,2)+EX(J)*CEN(J,2))*PEXS(L) PEWCZ(L)=(EX(I)*CEN(I,3)+EX(J)*CEN(J,3))*PEXS(L) 20 CONTINUE C C SET UP OTHER INDEX ARRAY FOR PACKED SYMMETRIC ARRAY C STORAGE C ID(I-IS)=L-IPX 30 CONTINUE IESPS=0 C C LOOP OVER ESP PROBE POINTS C DO 250 IESP=IESPS+1,NESP POTP1=POTPT(1,IESP) POTP2=POTPT(2,IESP) POTP3=POTPT(3,IESP) C CALCULATE QUANTITY U C L=0 DO 60 I=NP,IPE DO 60 J=I,IPE L=L+1 PTD(L)=((PEWCX(L)-POTP1)**2+(PEWCY(L)-POTP2)**2+ 1 (PEWCZ(L)-POTP3)**2)/PEXS(L) PCE(L)=SQRT(PI/PTD(L)) 60 CONTINUE C C CALCULATE F0, F1, AND F2(U) USING TAYLOR SERIES C OR ASYMPTOTIC EXPANSION C IL=L L=0 DO 100 I=1,IL IF(PTD(I) .LE. TF(0)) THEN IREF=DNINT(PTD(I)*20.D0) REF=0.05D0*IREF RES=PTD(I)-REF TERM=1.D0 PF0(I)=0.D0 DO 70 K=0,6 F=FV(K,IREF+1) TS=F*TERM*FAC(K) TERM=-TERM*RES PF0(I)=PF0(I)+TS 70 CONTINUE ELSE PF0(I)=PCE(I)*0.5D0 ENDIF IF(PTD(I) .LE. TF(1)) THEN IREF=DNINT(PTD(I)*20.D0) REF=0.05D0*IREF RES=PTD(I)-REF TERM1=1.D0 PF1(I)=0.D0 DO 80 K=0,6 FI=FV(K+1,IREF+1) TS1=FI*TERM1*FAC(K) TERM1=-TERM1*RES PF1(I)=PF1(I)+TS1 80 CONTINUE ELSE PF1(I)=PCE(I)*0.25D0/PTD(I) ENDIF IF(PTD(I) .LE. TF(2)) THEN IREF=DNINT(PTD(I)*20.D0) REF=0.05D0*IREF RES=PTD(I)-REF TERM2=1.D0 PF2(I)=0.D0 DO 90 K=0,6 FII=FV(K+2,IREF+1) TS2=FII*TERM2*FAC(K) TERM2=-TERM2*RES PF2(I)=PF2(I)+TS2 90 CONTINUE ELSE PF2(I)=PCE(I)*0.375D0/(PTD(I)*PTD(I)) ENDIF 100 CONTINUE C C CALCULATE (S||S) TYPE INTEGRALS C DO 110 I=1,IL PF0(I)=2.D0*PI*PEXS(I)*PEXPN(I)*PF0(I) PTD(I)=PF0(I) PF1(I)=2.D0*PI*PEXS(I)*PEXPN(I)*PF1(I) PF2(I)=2.D0*PI*PEXS(I)*PEXPN(I)*PF2(I) 110 CONTINUE C DO 230 IC=ISC+1,NC IPR=IC*ICD-ICD+1 ISTART=IPR DO 200 J=1,ICD C C CALCULATE (P||S) ESP INTEGRALS C IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN DO 150 I=ISTART,NPR IN=IPR+J-1 IR=IRD(I)+ID(IRD(IN)) IR2=ID(IRD(I))+IRD(IN) IF(IR2 .LE. IR ) IR=IR2 GO TO (120,130,140),IAM(IN,2) 120 NAI2(I,J)=(PEWCX(IR)-CEN(IN,1))*PF1(IR)-PF2(IR)* 1 (PEWCX(IR)-POTP1) NAI(I,J)=(PEWCX(IR)-CEN(IN,1))*PF0(IR)-PF1(IR)* 1 (PEWCX(IR)-POTP1) GO TO 150 130 NAI2(I,J)=(PEWCY(IR)-CEN(IN,2))*PF1(IR)-PF2(IR)* 1 (PEWCY(IR)-POTP2) NAI(I,J)=(PEWCY(IR)-CEN(IN,2))*PF0(IR)-PF1(IR)* 1 (PEWCY(IR)-POTP2) GO TO 150 140 NAI2(I,J)=(PEWCZ(IR)-CEN(IN,3))*PF1(IR)-PF2(IR)* 1 (PEWCZ(IR)-POTP3) NAI(I,J)=(PEWCZ(IR)-CEN(IN,3))*PF0(IR)-PF1(IR)* 1 (PEWCZ(IR)-POTP3) 150 CONTINUE ENDIF C C CALCULATE (P||P) ESP INTEGRALS C IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN DO 190 I=ISTART,NPR IN=IPR+J-1 IR=IRD(I)+ID(IRD(IN)) IR2=ID(IRD(I))+IRD(IN) IF(IR2 .LE. IR ) IR=IR2 GO TO (160,170,180),IAM(I,2) 160 NAI(I,J)=(PEWCX(IR)-CEN(I,1))*NAI(I,J)-(PEWCX(IR)-P 1OTP1)* NAI2(I,J) IF(IAM(IN,2) .EQ. IAM(I,2)) NAI(I,J)=NAI(I,J)+PEXS( 1IR)* 0.5D0*(PTD(IR)-PF1(IR)) GO TO 190 170 NAI(I,J)=(PEWCY(IR)-CEN(I,2))*NAI(I,J)-(PEWCY(IR)-P 1OTP2)* NAI2(I,J) IF(IAM(IN,2) .EQ. IAM(I,2)) NAI(I,J)=NAI(I,J)+PEXS( 1IR)* 0.5D0*(PTD(IR)-PF1(IR)) GO TO 190 180 NAI(I,J)=(PEWCZ(IR)-CEN(I,3))*NAI(I,J)-(PEWCZ(IR)-P 1OTP3)* NAI2(I,J) IF(IAM(IN,2) .EQ. IAM(I,2)) NAI(I,J)=NAI(I,J)+PEXS( 1IR)* 0.5D0*(PTD(IR)-PF1(IR)) 190 CONTINUE ENDIF 200 CONTINUE C C FORM INTEGRALS OVER CONTRACTED FUNCTIONS C IPS=IC*ICD-ICD+1 DO 220 I=IC,NC JPS=I*ICD-ICD+1 ESPI(I,IC)=0.D0 DO 210 J=JPS,JPS+ICD-1 DO 210 K=IPS,IPS+ICD-1 ESPI(I,IC)=ESPI(I,IC)+CC(J)*CC(K)*NAI(J,K-IPS+1) 210 CONTINUE ES(IESP)=ES(IESP)+2.D0*CESPM(INDC(I),INDC(IC))*ESPI(I,IC) 220 CONTINUE ES(IESP)=ES(IESP)-CESPM(INDC(IC),INDC(IC))*ESPI(IC,IC) 230 CONTINUE IF(MOD(IESP,NESP/IDN) .EQ. 0) THEN JSTART=ISC*2 IDC=IDC+1 ENDIF 250 CONTINUE RETURN END SUBROUTINE RSP(A,N,MATZ,W,Z) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INCLUDE 'SIZES' DIMENSION A(*), W(N), Z(N,N) ******************************************************************* * * EISPACK DIAGONALIZATION ROUTINES: TO FIND THE EIGENVALUES AND * EIGENVECTORS (IF DESIRED) OF A REAL SYMMETRIC PACKED MATRIX. * ON INPUT- N IS THE ORDER OF THE MATRIX A, * A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC * PACKED MATRIX STORED ROW-WISE, * MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF ONLY * EIGENVALUES ARE DESIRED, OTHERWISE IT IS SET TO * ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND * EIGENVECTORS. * ON OUTPUT- W CONTAINS THE EIGENVALUES IN ASCENDING ORDER, * Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO, * ******************************************************************* * THIS SUBROUTINE WAS CHOSEN AS BEING THE MOST RELIABLE. (JJPS) C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C DIMENSION FV1(MAXHEV*4+MAXLIT*3),FV2(MAXHEV*4+MAXLIT*3) SAVE FIRST, EPS, ETA, NV LOGICAL FIRST DATA FIRST /.TRUE./ IF (FIRST) THEN FIRST=.FALSE. CALL EPSETA(EPS,ETA) ENDIF NV=(N*(N+1))/2 NM=N CALL TRED3(N,NV,A,W,FV1,FV2,EPS,EPS) IF (MATZ .NE. 0) GO TO 10 C ********** FIND EIGENVALUES ONLY ********** CALL TQLRAT(N,W,FV2,IERR,EPS) GO TO 40 C ********** FIND BOTH EIGENVALUES AND EIGENVECTORS ********** 10 DO 30 I = 1, N C DO 20 J = 1, N Z(J,I)=0.0D0 20 CONTINUE C Z(I,I)=1.0D0 30 CONTINUE C CALL TQL2(NM,N,W,FV1,Z,IERR,EPS) IF (IERR .NE. 0) GO TO 40 CALL TRBAK3(NM,N,NV,A,N,Z,EPS) C ********** LAST CARD OF RSP ********** 40 RETURN END SUBROUTINE EPSETA(EPS,ETA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C COMPUTE AND RETURN ETA, THE SMALLEST REPRESENTABLE NUMBER, C AND EPS IS THE SMALLEST NUMBER FOR WHICH 1+EPS.NE.1. C C ETA = 1.D0 10 IF((ETA/2.D0).EQ.0.D0) GOTO 20 IF(ETA.LT.1.D-38) GOTO 20 ETA = ETA / 2.D0 GOTO 10 20 EPS = 1.D0 30 IF((1.D0+(EPS/2.D0)).EQ.1.D0) GOTO 40 IF(EPS.LT.1.D-17) GOTO 40 EPS = EPS / 2.D0 GOTO 30 40 RETURN END SUBROUTINE TQL2(NM,N,D,E,Z,IERR,EPS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C ===== PROCESSED BY AUGMENT, VERSION 4N ===== C APPROVED FOR VAX 11/780 ON MAY 6,1980. J.D.NEECE C ----- LOCAL VARIABLES ----- C ----- GLOBAL VARIABLES ----- DIMENSION D(N), E(N), Z(NM,N) C ----- SUPPORTING PACKAGE FUNCTIONS ----- C ===== TRANSLATED PROGRAM ===== C C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2, C NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND C WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD. C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS C FULL MATRIX TO TRIDIAGONAL FORM. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX, C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY, C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN C THE IDENTITY MATRIX. C C ON OUTPUT- C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT C UNORDERED FOR INDICES 1,2,...,IERR-1, C C E HAS BEEN DESTROYED, C C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED C EIGENVALUES, C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C IERR = 0 IF (N .EQ. 1) GO TO 160 C DO 10 I = 2, N 10 E(I-1) = E(I) C F=0.0D0 B=0.0D0 E(N)=0.0D0 C DO 110 L = 1, N J = 0 H=EPS*(ABS (D(L))+ABS (E(L))) IF (B .LT. H) B=H C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT ********** DO 20 M = L, N IF (ABS (E(M)).LE.B) GO TO 30 C ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP ********** 20 CONTINUE C 30 IF (M .EQ. L) GO TO 100 40 IF (J .EQ. 30) GO TO 150 J = J + 1 C ********** FORM SHIFT ********** L1 = L + 1 G = D(L) P=(D(L1)-G)/(2.0D0*E(L)) R=SQRT (P*P+1.0D0) D(L)=E(L)/(P+SIGN (R,P)) H = G - D(L) C DO 50 I = L1, N 50 D(I) = D(I) - H C F = F + H C ********** QL TRANSFORMATION ********** P = D(M) C=1.0D0 S=0.0D0 MML = M - L C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** DO 90 II = 1, MML I = M - II G = C * E(I) H = C * P IF (ABS (P).LT.ABS (E(I))) GO TO 60 C = E(I) / P R=SQRT (C*C+1.0D0) E(I+1) = S * P * R S = C / R C=1.0D0/R GO TO 70 60 C = P / E(I) R=SQRT (C*C+1.0D0) E(I+1) = S * E(I) * R S=1.0D0/R C = C * S 70 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C ********** FORM VECTOR ********** DO 80 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 80 CONTINUE C 90 CONTINUE C E(L) = S * P D(L) = C * P IF (ABS (E(L)).GT.B) GO TO 40 100 D(L) = D(L) + F 110 CONTINUE C ********** ORDER EIGENVALUES AND EIGENVECTORS ********** DO 140 II = 2, N I = II - 1 K = I P = D(I) C DO 120 J = II, N IF (D(J) .GE. P) GO TO 120 K = J P = D(J) 120 CONTINUE C IF (K .EQ. I) GO TO 140 D(K) = D(I) D(I) = P C DO 130 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 130 CONTINUE C 140 CONTINUE C GO TO 160 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** 150 IERR = L 160 RETURN C ********** LAST CARD OF TQL2 ********** END SUBROUTINE TQLRAT(N,D,E2,IERR,EPS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C ===== PROCESSED BY AUGMENT, VERSION 4N ===== C APPROVED FOR VAX 11/780 ON MAY 6,1980. J.D.NEECE C ----- LOCAL VARIABLES ----- C ----- GLOBAL VARIABLES ----- DIMENSION D(N), E2(N) C ----- SUPPORTING PACKAGE FUNCTIONS ----- C ===== TRANSLATED PROGRAM ===== C C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT, C ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH. C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC C TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD. C C ON INPUT- C C N IS THE ORDER OF THE MATRIX, C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX, C C E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE C INPUT MATRIX IN ITS LAST N-1 POSITIONS. E2(1) IS ARBITRARY. C C ON OUTPUT- C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES, C C E2 HAS BEEN DESTROYED, C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C IERR = 0 IF (N .EQ. 1) GO TO 140 C DO 10 I = 2, N 10 E2(I-1) = E2(I) C F=0.0D0 B=0.0D0 E2(N)=0.0D0 C DO 120 L = 1, N J = 0 H=EPS*(ABS (D(L))+SQRT (E2(L))) IF (B .GT. H) GO TO 20 B = H C = B * B C ********** LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ********** 20 DO 30 M = L, N IF (E2(M) .LE. C) GO TO 40 C ********** E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP ********** 30 CONTINUE C 40 IF (M .EQ. L) GO TO 80 50 IF (J .EQ. 30) GO TO 130 J = J + 1 C ********** FORM SHIFT ********** L1 = L + 1 S=SQRT (E2(L)) G = D(L) P=(D(L1)-G)/(2.0D0*S) R=SQRT (P*P+1.0D0) D(L)=S/(P+SIGN (R,P)) H = G - D(L) C DO 60 I = L1, N 60 D(I) = D(I) - H C F = F + H C ********** RATIONAL QL TRANSFORMATION ********** G = D(M) IF (G.EQ.0.0D0) G=B H = G S=0.0D0 MML = M - L C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** DO 70 II = 1, MML I = M - II P = G * H R = P + E2(I) E2(I+1) = S * R S = E2(I) / R D(I+1) = H + S * (H + D(I)) G = D(I) - E2(I) / G IF (G.EQ.0.0D0) G=B H = G * P / R 70 CONTINUE C E2(L) = S * G D(L) = H C ********** GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST ********** IF (H.EQ.0.0D0) GO TO 80 IF (ABS (E2(L)).LE.ABS (C/H)) GO TO 80 E2(L) = H * E2(L) IF (E2(L).NE.0.0D0) GO TO 50 80 P = D(L) + F C ********** ORDER EIGENVALUES ********** IF (L .EQ. 1) GO TO 100 C ********** FOR I=L STEP -1 UNTIL 2 DO -- ********** DO 90 II = 2, L I = L + 2 - II IF (P .GE. D(I-1)) GO TO 110 D(I) = D(I-1) 90 CONTINUE C 100 I = 1 110 D(I) = P 120 CONTINUE C GO TO 140 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** 130 IERR = L 140 RETURN C ********** LAST CARD OF TQLRAT ********** END SUBROUTINE TRBAK3(NM,N,NV,A,M,Z,EPS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C ===== PROCESSED BY AUGMENT, VERSION 4N ===== C APPROVED FOR VAX 11/780 ON MAY 6,1980. J.D.NEECE C ----- LOCAL VARIABLES ----- C ----- GLOBAL VARIABLES ----- DIMENSION A(NV), Z(NM,M) C ===== TRANSLATED PROGRAM ===== C C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING C SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY TRED3. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT, C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS C USED IN THE REDUCTION BY TRED3 IN ITS FIRST C N*(N+1)/2 POSITIONS, C C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED, C C Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED C IN ITS FIRST M COLUMNS. C C ON OUTPUT- C C Z CONTAINS THE TRANSFORMED EIGENVECTORS C IN ITS FIRST M COLUMNS. C C NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C IF (M .EQ. 0) GO TO 50 IF (N .EQ. 1) GO TO 50 C DO 40 I = 2, N L = I - 1 IZ = (I * L) / 2 IK = IZ + I H = A(IK) IF (H.EQ.0.0D0) GO TO 40 C DO 30 J = 1, M S=0.0D0 IK = IZ C DO 10 K = 1, L IK = IK + 1 S = S + A(IK) * Z(K,J) 10 CONTINUE C ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ********** S = (S / H) / H IK = IZ C DO 20 K = 1, L IK = IK + 1 Z(K,J) = Z(K,J) - S * A(IK) 20 CONTINUE C 30 CONTINUE C 40 CONTINUE C 50 RETURN C ********** LAST CARD OF TRBAK3 ********** END SUBROUTINE TRED3(N,NV,A,D,E,E2,EPS,ETA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C ===== PROCESSED BY AUGMENT, VERSION 4N ===== C APPROVED FOR VAX 11/780 ON MAY 6,1980. J.D.NEECE C ----- LOCAL VARIABLES ----- C ----- GLOBAL VARIABLES ----- DIMENSION A(NV), D(N), E(N), E2(N) C ----- SUPPORTING PACKAGE FUNCTIONS ----- C ===== TRANSLATED PROGRAM ===== C C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS C A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX C USING ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT- C C N IS THE ORDER OF THE MATRIX, C C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT, C C A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC C INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL C ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS. C C ON OUTPUT- C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL C TRANSFORMATIONS USED IN THE REDUCTION, C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX, C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO, C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C ********** FOR I=N STEP -1 UNTIL 1 DO -- ********** TOL=ETA/EPS DO 100 II = 1, N I = N + 1 - II L = I - 1 IZ = ( I * L ) / 2 H=0.0D0 SCALE=0.0D0 DO 10 K = 1, L IZ = IZ + 1 D(K) = A(IZ) SCALE=SCALE+ABS( D(K) ) 10 CONTINUE C IF ( SCALE.NE.0.D0 ) GO TO 20 E(I)=0.0D0 E2(I)=0.0D0 GO TO 90 C 20 DO 30 K = 1, L D(K) = D(K) / SCALE H = H + D(K) * D(K) 30 CONTINUE C E2(I) = SCALE * SCALE * H F = D(L) G=-SIGN (SQRT (H),F) E(I) = SCALE * G H = H - F * G D(L) = F - G A(IZ) = SCALE * D(L) IF (L .EQ. 1) GO TO 90 F=0.0D0 C DO 70 J = 1, L G=0.0D0 JK = (J * (J-1)) / 2 C ********** FORM ELEMENT OF A*U ********** K = 0 40 K = K + 1 JK = JK + 1 G = G + A(JK) * D(K) IF ( K .LT. J ) GO TO 40 IF ( K .EQ. L ) GO TO 60 50 JK = JK + K K = K + 1 G = G + A(JK) * D(K) IF ( K .LT. L ) GO TO 50 C ********** FORM ELEMENT OF P ********** 60 CONTINUE E(J) = G / H F = F + E(J) * D(J) 70 CONTINUE C HH = F / (H + H) JK = 0 C ********** FORM REDUCED A ********** DO 80 J = 1, L F = D(J) G = E(J) - HH * F E(J) = G C DO 80 K = 1, J JK = JK + 1 A(JK) = A(JK) - F * E(K) - G * D(K) 80 CONTINUE C 90 D(I) = A(IZ+1) A(IZ+1)=SCALE*SQRT (H) 100 CONTINUE C RETURN C ********** LAST CARD OF TRED3 ********** END