C FORTICON8 (VAX VERSION) QCPE 517 C C THE FOLLOWING CODE IS NOT THE ORIGINAL SOURCE CODE. THIS CODE C INCLUDES MODIFICATIONS TO WRITE OUT PERTINENT INFORMATION TO A C TEMPORARY FILE (FOR013) TO BE USED IN MO PLOTTING ROUTINES. C C ALL MODIFICATIONS ARE COMMENTED TO PROVIDE EASE IN IDENTIFYING C CHANGES. JJN 8-7-90 C C********************************************************************** FORT0001 C * FORT0002 C PROGRAM FORTICON8 COMPLETE FORTRAN VERSION OF ICON8 * FORT0003 C * FORT0004 C THE FOLLOWING SUBPROGRAMS, WHICH EXIST AS ASSEMBLER ROUTINES * FORT0005 C IN ICON8, ARE TRANSLATED INTO FORTRAN0 MATRIX, ABFNS, * FORT0006 C LOVLAP, GRMSCH, TRNFRM, DSUM, ROTATE, DOT, VECSUM, REDUCE, * FORT0007 C FULCHM, AND REDCHM. FORTICON INLUDES THESE AS WELL AS ALL * FORT0008 C THE FORTRAN SUBPROGRAMS OF ICON8. * FORT0009 C * FORT0010 C********************************************************************** FORT0011 C FORT0012 IMPLICIT REAL*8(A-H,O-Z) FORT0013 C FORT0014 C PROGRAM ICON FOR PERFORMING EXTENDED HUCKEL CALCULATIONS FORT0015 C WITH OR WITHOUT CHARGE ITERATION. FORT0016 C ** QCPE VERSION ** FORT0017 C FORT0018 C ** SAMPLE DECK ** FORT0019 C ....0....1....0....2....0....3....0....4....0....5....0....6....0 FORT0020 C FORT0021 C ETHYLENE FORT0022 C 4 2 2 FORT0023 C 0.92665 1.205 0.0 FORT0024 C -0.92665 1.205 0.0 FORT0025 C 0.92665 -1.205 0.0 FORT0026 C -0.92665 -1.205 0.0 FORT0027 C 0.0 0.67 0.0 FORT0028 C 0.0 -0.67 0.0 FORT0029 C C C FORT0030 C FORT0031 C C REVISED TO ALLOW MORE FLEXIBILITY IN NUMBER OF ATOMS C POSSIBLE. MAXATM (MAXIMUM NUMBER OF ATOMS) SET AT 500. C MAXIMUM NUMBER OF HEAVY ATOMS SET AT 250. MAXIMUM NUMBER C OF USER-DEFINED ATOMS SET TO 230. C JJN 8-28-90 C PARAMETER (MAXATM=500) PARAMETER (BB=250) PARAMETER (MXUSER=230) PARAMETER (MXUSR2=231) C C THE ABOVE PARAMETERS ARE DEFINED AS SUCH TO ELIMINATE THE NEED C TO SEARCH THROUGH THE CODE FOR EVERY OCCURRENCE OF THE NUMBERS C (AS I HAD TO DO WITH THE ORIGINAL CODE). THE PARAMETER STATEMENTS C ARE DEFINED IN EACH SUBROUTINE REQUIRING THEM SO IF THESE NEED TO C BE CHANGED AGAIN, JUST FIND THE PARAMETER STATEMENTS. ALSO, THERE C ARE A FEW PLACES WHERE THE USE OF THE PARAMETERS IS NOT ALLOWED AND C THE ACTUAL NUMBER IS USED. MOST NOTABLY IN DATA STATEMENTS. A QUICK C SEARCH THROUGH FOR THE VALUES OF MXUSER AND MXUSR2 WILL BE C NECESSARY. JJN 9-8-90 C COMMON/TITLE/AB(10) FORT0032 COMMON/CNTRL/CON,PEEP,COULH,NH,NA,NATM,KA,NELEC,METH, FORT0033 . IPRINT,IPUNCH,L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT0034 LOGICAL*1 L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT0035 COMMON/OUT/PRT(20),PUN(20),IOVPOP(24),IENRGY(24) FORT0036 LOGICAL*1 PRT,PUN FORT0037 INTEGER*2 IOVPOP,IENRGY FORT0038 COMMON/ATOM/AC(BB),SYMBOL(BB),VELEC(BB),NS(BB),NP(BB), . ND(BB),EXPS(BB),EXPP(BB),EXPD(BB),EXPD2(BB),C1(BB),C2(BB), . COULS(BB),COULP(BB),COULD(BB),X(MAXATM),Y(MAXATM),Z(MAXATM) INTEGER*2 AC,SYMBOL,VELEC FORT0042 COMMON/ITPARM/DAMP1,DAMP2,DAMP3,LAMPRI,DELTAC,SENSE,MAXCYC, FORT0043 . PRTCYC,ICYCLE,NCON,PARTIT,PRINTX,ITABLE(20) FORT0044 REAL*8 LAMPRI FORT0045 INTEGER*4 PRTCYC FORT0046 LOGICAL*1 PARTIT,PRINTX,ITABLE FORT0047 COMMON/ABC/AS2(5),BS2(5),CS2(5),AP2(5),BP2(5),CP2(5),AD2(5), FORT0048 . BD2(5),CD2(5),AS3(5),BS3(5),CS3(5),AP3(5),BP3(5),CP3(5), FORT0049 . AD3(5),BD3(5),CD3(5) FORT0050 COMMON/STARS/STAR,STAR2 FORT0051 INTEGER*2 STAR,STAR2 FORT0052 C FORT0053 C CALL INPUT TO READ IN AND PRINT OUT INPUT DATA. FORT0054 C FORT0055 1000 CALL INPUT(NATOM,NDIM,NTYPE) FORT0056 IF(IPRINT.LT.-3) GO TO 2000 FORT0057 C FORT0058 C CALCULATE MATRIX DIMENSIONS. FORT0059 C FORT0060 NC=NDIM*(NDIM+1)/2 FORT0061 NHS=NC+NC-NDIM FORT0062 NSS=NHS FORT0063 C FORT0064 C IF NO WAVEFUNCTIONS NEEDED, SET NHS=0. THIS HAS THE EFFECT FORT0065 C OF EQUIVALENCING THE H AND S MATRICES. FORT0066 C FORT0067 ONEMAT=.FALSE. FORT0068 IF(.NOT.ITERAT.AND.IPRINT.LT.-1) ONEMAT=.TRUE. FORT0069 IF(IPUNCH.NE.0) GO TO 600 FORT0070 DO 400 I=6,20 FORT0071 IF(PUN(I)) GO TO 600 FORT0072 400 CONTINUE FORT0073 GO TO 500 FORT0074 600 ONEMAT=.FALSE. FORT0075 500 IF(ONEMAT) NHS=0 FORT0076 IF(METH.LT.3) NTYPE=1 FORT0077 NMD=NTYPE*NTYPE FORT0078 NCL=NDIM FORT0079 IF(NDIM.LT.10) NCL=10 FORT0080 NOCC=(NDIM+1)/2 FORT0081 NHDG=1 FORT0082 IF(METH.GT.2.AND.L5) NHDG=NDIM FORT0083 C FORT0084 C CALL MATRIX TO ALLOCATE SPACE FOR MATRICES. FORT0085 C ORDER OF MATRICES0 H S MAD C SP PD MAXS MAXP MAXD FORT0086 C COUL0 SORB IOCC HDG FORT0087 C FORT0088 200 CALL MATRIX(13,NHS,NSS,NMD,NC,NDIM,NDIM,NDIM,NDIM,2*NDIM, FORT0089 . NCL,NDIM,NOCC,NHDG, NDIM,NDIM,NC,NATOM,NTYPE,NHDG,NC,NC, FORT0090 . NC,NC,NC,NC,NDIM) FORT0091 2000 CONTINUE FORT0092 GO TO 1000 FORT0093 END FORT0094 BLOCK DATA FORT0095 C FORT0096 C INITIALIZATION OF INTERNAL ATOMIC DATA. THERE ARE PROVISIONS FORT0097 C FOR 20 USER DEFINED ATOMS, 15 INTERNALLY DEFINED ATOMS, AND FORT0098 C SPACE FOR 5 MORE TO BE USED EITHER WAY. FORT0099 C FORT0100 C C THIS HAS BEEN MODIFIED TO ALLOW 230 USER DEFINED ATOMS, 15 C INTERNALLY DEFINED ATOMS AND SPACE FOR 5 MORE TO BE USED EITHER C WAY (TOTAL=250). JJN 8-28-90 C PARAMETER (MAXATM=500) PARAMETER (BB=250) PARAMETER (MXUSER=230) PARAMETER (MXUSR2=231) IMPLICIT REAL*8(A-H,O-Z) COMMON/ATOM/KEY(BB),SYMBOL(BB),VELEC(BB),NS(BB),NP(BB), . ND(BB),EXPS(BB),EXPP(BB),EXPD(BB),EXPD2(BB),C1(BB), . C2(BB),COULS(BB),COULP(BB),COULD(BB),X(MAXATM),Y(MAXATM), . Z(MAXATM) INTEGER*2 KEY,SYMBOL,VELEC COMMON/STARS/STAR,STAR2 INTEGER*2 STAR,STAR2 COMMON/START/NUSER DATA SYMBOL/230*'**',' C',' N',' O',' F','SI',' P',' S', . 'CL','LI','BE',' B','NA','MG','AL',' H',5*' '/ DATA VELEC/230*0,4,5,6,7,4,5,6,7,1,2,3,1,2,3,1,5*0/ DATA NS/230*0,4*2,4*3,3*2,3*3,1,5*0/ DATA EXPS/230*0.0D0,1.625D0,1.950D0,2.275D0,2.425D0,1.383D0, . 1.6D0,1.817D0,2.033D0,0.65D0,0.975D0,1.3D0,0.733D0,0.95D0, FORT0114 . 1.167D0,1.3D0,5*0.0D0/ DATA COULS/230*0.0D0,-21.4D0,-26.0D0,-32.3D0,-40.0D0,-17.3D0, . -18.6D0,-20.0D0,-30.0D0,-5.4D0,-10.0D0,-15.2D0,-5.1D0,-9.0D0, . -12.3D0,-13.6D0,5*0.0D0/ DATA NP/230*0,4*2,4*3,3*2,3*3,6*0/ DATA EXPP/230*0.0D0,1.625D0,1.950D0,2.275D0,2.425D0,1.383D0, . 1.6D0,1.817D0,2.033D0,0.65D0,0.975D0,1.3D0,0.733D0,0.95D0, FORT0121 . 1.167D0,6*0.0D0/ DATA COULP/230*0.0D0,-11.4D0,-13.4D0,-14.8D0,-18.1D0,-9.2D0, . -14.0D0,-13.3D0,-15.0D0,-3.5D0,-6.0D0,-8.5D0,-3.0D0,-4.5D0, FORT0124 . -6.5D0,6*0.0D0/ DATA ND/234*0,4*3,12*0/ DATA EXPD/234*0.0D0,1.383D0,1.4D0,1.5D0,2.033D0,12*0.0D0/ DATA COULD/234*0.0D0,-6.0D0,-7.0D0,-8.0D0,-9.0D0,12*0.0D0/ DATA C1/250*0.0D0/ DATA C2/250*0.0D0/ DATA EXPD2/250*0.0D0/ DATA STAR,STAR2 / ' *','**'/ FORT0132 C C PROVISION FOR 230 USER DEFINED ATOMS. JJN 8-28-90 C DATA NUSER/231/ C END FORT0134 SUBROUTINE INPUT(NATOM,NDIM,NTYPE) FORT0135 C FORT0136 C SUBROUTINE FOR READING IN AND PRINTING OUT INPUT DATA. FORT0137 C FORT0138 PARAMETER (MAXATM=500) PARAMETER (BB=250) PARAMETER (MXUSER=230) PARAMETER (MXUSR2=231) IMPLICIT REAL*8(A-H,O-Z) FORT0139 COMMON/TITLE/AB(10) FORT0140 COMMON/CNTRL/CON,PEEP,COULH,NH,NA,NATM,KA,NELEC,METH, FORT0141 . IPRINT,IPUNCH,L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT0142 LOGICAL*1 L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT0143 COMMON/OUT/PRT(20),PUN(20),IOVPOP(24),IENRGY(24) FORT0144 LOGICAL*1 PRT,PUN FORT0145 INTEGER*2 IOVPOP,IENRGY FORT0146 COMMON/ATOM/AC(BB),SYMBOL(BB),VELEC(BB),NS(BB),NP(BB), . ND(BB),EXPS(BB),EXPP(BB),EXPD(BB),EXPD2(BB),C1(BB), . C2(BB),COULS(BB),COULP(BB),COULD(BB),X(MAXATM),Y(MAXATM), . Z(MAXATM),SYMBL(MAXATM),ATS(89) INTEGER*2 AC,SYMBOL,SYMBL,VELEC,ATS COMMON/ITPARM/DAMP1,DAMP2,DAMP3,LAMPRI,DELTAC,SENSE,MAXCYC, FORT0151 . PRTCYC,ICYCLE,NCON,PARTIT,PRINTX,ITABLE(20) FORT0152 REAL*8 LAMPRI FORT0153 INTEGER*4 PRTCYC FORT0154 LOGICAL*1 PARTIT,PRINTX,ITABLE FORT0155 C FORT0156 C SINCE THE ITERNAL ATOMIC PARAMETERS ( EXPS, EXPP, ETC. ) ARE FORT0157 C NOT USED WHEN DOING CHARGE ITERATION ( METH >1 ) THE SPACE FORT0158 C ALLOCATED TO THEM CAN BE USED FOR THE VSIE CHARGE ITERATION FORT0159 C PARAMETERS. FORT0160 C FORT0161 DIMENSION AS1(MXUSER),BS1(MXUSER),CS1(MXUSER),AP1(MXUSER), . BP1(MXUSER),CP1(MXUSER),AD1(MXUSER),BD1(MXUSER),CD1(MXUSER) EQUIVALENCE (AS1(1),EXPS(MXUSR2)),(BS1(1),EXPP(MXUSR2)), . (CS1(1),EXPD(MXUSR2)),(AP1(1),EXPD2(MXUSR2)), . (BP1(1),C1(MXUSR2)),(CP1(1),C2(MXUSR2)), . (AD1(1),COULS(MXUSR2)),(BD1(1),COULP(MXUSR2)), . (CD1(1),COULD(MXUSR2)) COMMON/ABC/AS2(5),BS2(5),CS2(5),AP2(5),BP2(5),CP2(5),AD2(5), FORT0168 . BD2(5),CD2(5),AS3(5),BS3(5),CS3(5),AP3(5),BP3(5),CP3(5), FORT0169 . AD3(5),BD3(5),CD3(5) FORT0170 INTEGER*2 CHANGE(MXUSER) EQUIVALENCE (CHANGE(1),AS2(1)) FORT0172 REAL*4 MADS(MXUSER),MADP(MXUSER),MADD(MXUSER) EQUIVALENCE (MADS(1),NS(MXUSR2)),(MADP(1),NP(MXUSR2)), . (MADD(1),ND(MXUSR2)) COMMON/STARS/STAR,STAR2 FORT0176 INTEGER*2 STAR,STAR2 FORT0177 COMMON/START/NUSER FORT0178 DIMENSION EXTRA(9) FORT0179 EQUIVALENCE (X(1),EXTRA(1)) FORT0180 INTEGER*2 CONTIN FORT0181 EQUIVALENCE (AB(10),CONTIN) FORT0182 INTEGER*2 HYDROG FORT0183 DATA HYDROG/' H'/ FORT0184 DATA ATS/' H','HE','LI','BE',' B',' C',' N',' O',' F', . 'NE','NA','MG','AL','SI',' P',' S','CL','AR',' K','CA', . 'SC','TI',' V','CR','MN','FE','CO','NI','CU','ZN','GA', . 'GE','AS','SE','BR','KR','RB','SR',' Y','ZR','NB','MO', . 'TC','RU','RH','PD','AG','CD','IN','SN','SB','TE',' I', . 'XE','CS','BA','LA','CE','PR','ND','PM','SM','EU','GD', . 'TB','DY','HO','ER','TM','YB','LU','HF','TA',' W','RE', . 'OS','IR','PT','AU','HG','TL','PB','BI','PO','AT','RN', . 'FR','RA','AC'/ C FORT0185 C READ AND WRITE TITLE. FORT0186 C IF CONTIN IS EQUAL TO STAR THEN ANOTHER TITLE CARD WILL FORT0187 C BE READ AND PRINTED. HOWEVER ONLY THE FIRST IS STORED FORT0188 C FOR PRINTING LATER ON. GET YOUR GOODIES ON THE FIRST. FORT0189 C FORT0190 1000 READ(5,1,END=115) AB FORT0191 1 FORMAT(8A8,A6,A2) FORT0192 C WRITE TITLE TO DISK FILE 13. JJN 9-3-90 C WRITE(13,2730) AB 2730 FORMAT(8A8,A6,A2) C WRITE(6,2) AB FORT0193 2 FORMAT('1',T10,8A8,A6,A2) FORT0194 11 IF(CONTIN.NE.STAR) GO TO 9 FORT0195 READ(5,1) EXTRA,CONTIN FORT0196 WRITE(6,12) EXTRA,CONTIN FORT0197 12 FORMAT(T10,8A8,A6,A2) FORT0198 GO TO 11 FORT0199 9 CONTINUE FORT0200 C FORT0201 C READ PARAMETER CARD. FORT0202 C FORT0203 READ(5,3) NH,NA,KA,METH,IPRINT,IPUNCH,L1,L2,L3,L4,L5,CON, FORT0204 . PEEP,COULH,(PRT(I),I=1,20),(PUN(J),J=1,20) FORT0205 3 FORMAT(6I3,5L1,F5.2,2F6.3,40L1) FORT0206 C FORT0207 C INSERT DEFAULT PARAMETERS. FORT0208 C FORT0209 IF(CON.LT.1.E-05) CON=1.75D0 FORT0210 IF(PEEP.LT.1.E-05) PEEP=1.3D0 FORT0211 IF(COULH.GT.-1.E-05) COULH=-13.6D0 FORT0212 ITERAT=METH.NE.0 FORT0213 NATOM=NH+NA FORT0214 NATM=NATOM FORT0215 C FORT0216 C SET IPRINT OPTION. FORT0217 C FORT0218 IF(IPRINT.GT.1) GO TO 250 FORT0219 PRT(6)=.TRUE. FORT0220 PRT(7)=.TRUE. FORT0221 PRT(11)=.TRUE. FORT0222 PRT(17)=.TRUE. FORT0223 PRT(19)=.TRUE. FORT0224 IF(IPRINT.GT.0) GO TO 250 FORT0225 PRT(12)=.TRUE. FORT0226 PRT(14)=.TRUE. FORT0227 PRT(20)=.TRUE. FORT0228 IF(IPRINT.GT.-1) GO TO 250 FORT0229 PRT(13)=.TRUE. FORT0230 PRT(15)=.TRUE. FORT0231 PRT(16)=.TRUE. FORT0232 PRT(18)=.TRUE. FORT0233 IF(IPRINT.GT.-2) GO TO 250 FORT0234 PRT(10)=.TRUE. FORT0235 C FORT0236 C READ COORDINATES AND HEAVY ATOM CARD. FORT0237 C FORT0238 250 READ(5,5) (X(I),Y(I),Z(I),I=1,NATOM) FORT0239 5 FORMAT(3F15.6) FORT0240 READ(5,8) (AC(I),I=1,NA) FORT0241 8 FORMAT(40A2) FORT0242 C FORT0243 C READ AND DECODE ATOM DEFINITION CARDS. FORT0244 C FORT0245 JOHN=0 NDIM=NH FORT0246 NTYPE=NH FORT0247 NELEC=NH-KA FORT0248 K=NUSER FORT0249 NUSER2=BB IF(METH.GE.2) NUSER2=MXUSER DO 100 I=1,NA FORT0252 IF(NUSER.GT.NUSER2) GO TO 103 FORT0253 DO 102 J=NUSER,NUSER2 FORT0254 JSAVE=J FORT0255 IF(AC(I) .EQ. SYMBOL(J)) GO TO 101 FORT0256 102 CONTINUE FORT0257 C FORT0258 C PROVISION FOR USER SPECIFIED DATA. FORT0259 C FORT0260 IF(AC(I).EQ.STAR) GO TO 103 FORT0261 IF(AC(I).EQ.STAR2) GO TO 105 FORT0262 WRITE(6,6) I,AC(I) FORT0263 6 FORMAT(//,T10,'HEAVY ATOM',I3,' NOT RECOGNIZED. SYMBOL0',A2) FORT0264 IF(METH.GE.2) WRITE(6,13) FORT0265 13 FORMAT(/,T10,'REMEMBER IF USING METH > 1 ALL ATOMIC', FORT0266 . ' PARAMETERS MUST BE DEFINED BY THE USER.') FORT0267 115 REWIND 7 FORT0268 STOP FORT0269 103 NUSER=NUSER-1 FORT0270 105 READ(5,7) SYMBOL(NUSER),VELEC(NUSER),NS(NUSER),EXPS(NUSER), FORT0271 . COULS(NUSER),NP(NUSER),EXPP(NUSER),COULP(NUSER),ND(NUSER), FORT0272 . EXPD(NUSER),COULD(NUSER),C1(NUSER),EXPD2(NUSER),C2(NUSER) FORT0273 7 FORMAT(A2,I3,3(I3,2F6.3),F6.4,F6.3,F6.4) C C THIS SECTION WRITE OUT THE USER-DEFINED PARAMETERS TO DISK C FILE 18. THESE WILL BE USED FOR CONSTRUCTION OF THE PSI1 C INPUT FILE AS THESE PARAMETERS ARE REQUIRED FOR ELEMENTS C GREATER THAN ATOMIC NUMBER 18. JJN 9-8-90 C JOHN=JOHN+1 WRITE(18,767) SYMBOL(NUSER),VELEC(NUSER),NS(NUSER), . EXPS(NUSER),NP(NUSER),EXPP(NUSER),ND(NUSER),EXPD(NUSER), . C1(NUSER),EXPD2(NUSER),C2(NUSER) 767 FORMAT(A2,I3,I3,F6.3,I3,F6.3,I3,F6.3,F6.4,F6.3,F6.4) C C JSAVE=NUSER FORT0275 C FORT0276 C NORMALIZE USER SPECIFIED CONTRACTED D ORBITAL. FORT0277 C FORT0278 IF(C2(NUSER).EQ.0.) GO TO 101 FORT0279 S=(4.D0*EXPD(NUSER)*EXPD2(NUSER)/(EXPD(NUSER)+EXPD2(NUSER) FORT0280 . )**2)**(ND(NUSER)+.5D0) FORT0281 S=1.D0/DSQRT(C1(NUSER)**2+C2(NUSER)**2+(S+S)*C1(NUSER) FORT0282 . *C2(NUSER)) FORT0283 C1(NUSER)=S*C1(NUSER) FORT0284 C2(NUSER)=S*C2(NUSER) FORT0285 101 NELEC=NELEC+VELEC(JSAVE) FORT0286 C FORT0287 C AC, LATER REFERENCED AS KEY, IS A POINTER TO THE PARAMETER TABLES.FORT0288 C FORT0289 111 AC(I)=JSAVE FORT0290 NDIM=NDIM+4 FORT0291 IF(NP(JSAVE).EQ.0) NDIM=NDIM-3 FORT0292 IF(ND(JSAVE).NE.0) NDIM=NDIM+5 FORT0293 NTYPE=NTYPE+2 FORT0294 IF(NP(JSAVE).EQ.0) NTYPE=NTYPE-1 FORT0295 IF(ND(JSAVE).NE.0) NTYPE=NTYPE+1 FORT0296 100 CONTINUE FORT0297 C FORT0298 C READ IN CHARGE ITERATION PARAMETERS. SET DEFAULT VALUES. FORT0299 C FORT0300 IF(.NOT.ITERAT) GO TO 60 FORT0301 IF(K.NE.MXUSR2) GO TO 60 READ(5,61) DAMP1,DAMP2,DAMP3,LAMPRI,DELTAC,SENSE,MAXCYC, FORT0303 . PRTCYC,NCON,PARTIT FORT0304 61 FORMAT(6F10.5,3I5,4X,L1) FORT0305 IF(.NOT.PARTIT.OR.METH.EQ.2) GO TO 65 FORT0306 WRITE(6,66) FORT0307 66 FORMAT(///,T10,'PARTIAL ITERATION ( PARTIT = TRUE ) MAY', FORT0308 . ' ONLY BE USED IF METH = 2.') FORT0309 STOP FORT0310 65 IF(DELTAC.EQ.0.0D0) DELTAC=0.0001D0 FORT0311 IF(SENSE.EQ.0.0D0) SENSE=2.0D0 FORT0312 IF(MAXCYC.EQ.0) MAXCYC=25 FORT0313 IF(PRTCYC.EQ.0) PRTCYC=MAXCYC FORT0314 IF(NCON.EQ.0) NCON=3 FORT0315 IF(DAMP1.EQ.0.0D0) DAMP1=0.1D0 FORT0316 IF(METH.GE.3) GO TO 62 FORT0317 IF(DAMP2.EQ.0.0D0) DAMP2=0.25D0 FORT0318 IF(LAMPRI.EQ.0.0D0) LAMPRI=0.25D0 FORT0319 GO TO 63 FORT0320 62 IF(DAMP2.EQ.0.0D0) DAMP2=0.75D0 FORT0321 IF(LAMPRI.EQ.0.0D0) LAMPRI=0.75D0 FORT0322 63 IF(METH.LT.2) GO TO 60 FORT0323 DO 32 I=1,20 FORT0324 32 ITABLE(I)=.FALSE. FORT0325 NUSER2=MXUSR2-NUSER C FORT0327 C READ IN SYMBOLS OF ATOMS ON WHICH CHARGE ITERATION FORT0328 C IS TO BE PERFORMED. FORT0329 C FORT0330 IF(.NOT.PARTIT) GO TO 30 FORT0331 READ(5,31) CHANGE FORT0332 31 FORMAT(20A2) FORT0333 DO 33 I=1,NUSER2 FORT0334 J=MXUSR2-I DO 33 K=1,NUSER2 FORT0336 IF(SYMBOL(J).EQ.CHANGE(K)) ITABLE(J)=.TRUE. FORT0337 33 CONTINUE FORT0338 GO TO 34 FORT0339 30 DO 35 I=1,NUSER2 FORT0340 J=MXUSR2-I 35 ITABLE(J)=.TRUE. FORT0342 C FORT0343 C READ IN VSIE AND MADELUNG PARAMETERS. FORT0344 C FORT0345 34 DO 36 I=1,NUSER2 FORT0346 J=MXUSR2-I IF(.NOT.ITABLE(J)) GO TO 36 FORT0348 READ(5,37) AS1(I),BS1(I),CS1(I),MADS(I) FORT0349 37 FORMAT(4F10.8) FORT0350 IF(NP(J).EQ.0) GO TO 36 FORT0351 IF(ND(J).NE.0) GO TO 38 FORT0352 READ(5,37) AP1(I),BP1(I),CP1(I),MADP(I) FORT0353 GO TO 36 FORT0354 38 IF(NCON.EQ.3) GO TO 39 FORT0355 READ(5,37) AP1(I),BP1(I),CP1(I),MADP(I),AD1(I),BD1(I), FORT0356 . CD1(I),MADD(I) FORT0357 GO TO 36 FORT0358 39 READ(5,40) AS2(I),BS2(I),CS2(I),AS3(I),BS3(I),CS3(I) FORT0359 40 FORMAT(3F10.8) FORT0360 READ(5,37) AP1(I),BP1(I),CP1(I),MADP(I) FORT0361 READ(5,40) AP2(I),BP2(I),CP2(I),AP3(I),BP3(I),CP3(I) FORT0362 READ(5,37) AD1(I),BD1(I),CD1(I),MADD(I) FORT0363 READ(5,40) AD2(I),BD2(I),CD2(I),AD3(I),BD3(I),CD3(I) FORT0364 36 CONTINUE FORT0365 C FORT0366 C READ IN IOVPOP(I) AND IENRGY(I). INDIVIDULAL OVERLAP POPULATION FORT0367 C ANALYSES ARE PERFORMED FROM ORBITAL IOVPOP(N) TO ORBITAL FORT0368 C IOVPOP(N+1). INDIVIDUAL ENERGY MATRIX ANALYSES ARE PERFORMED FORT0369 C FROM ORBITAL IENRGY(N) TO ORBITAL IENRGY(N+1). FORT0370 C FORT0371 60 IF(L3) READ(5,67) IOVPOP FORT0372 67 FORMAT(24I3) FORT0373 IF(L4) READ(5,67) IENRGY FORT0374 C FORT0375 C PRINT OUT TYPE OF CALCULATION. FORT0376 C FORT0377 IF(METH.EQ.0) WRITE(6,90) FORT0378 IF(METH.EQ.1) WRITE(6,91) FORT0379 IF(METH.GE.2) WRITE(6,92) FORT0380 IF(METH.GT.2) WRITE(6,93) FORT0381 IF(L5) WRITE(6,94) FORT0382 90 FORMAT(///,T10,'EXTENDED HUCKEL CALCULATION.') FORT0383 91 FORMAT(///,T10,'EXTENDED HUCKEL CALCULATION WITH CHARGE', FORT0384 . ' ITERATION.',/,T10,'LINEAR CHARGE DEPENDENCE OF SENSE*CHARG'FORT0385 . ,'E FOR H(I,I)''S.') FORT0386 92 FORMAT(///,T10,'EXTENDED HUCKEL CALCULATION WITH CHARGE', FORT0387 . ' ITERATION.') FORT0388 93 FORMAT(T10,'MADELUNG CORRECTION INCLUDED.') FORT0389 94 FORMAT(T10,'WEIGHTED HIJ FORMULA USED.') FORT0390 C FORT0391 C PRINT OUT ATOMIC COORDINATES AND PARAMETERS. FORT0392 C FORT0393 IF(PRT(1)) GO TO 80 FORT0394 WRITE(6,74) FORT0395 74 FORMAT(///,T5,'ATOM',T17,'X',T29,'Y',T41,'Z',T56,'S',T76,'P',FORT0396 . T96,'D',T113,'CONTRACTED D'/T47,'N',T50,'EXP',T59,'COUL', FORT0397 . T67,'N',T70,'EXP',T79,'COUL',T87,'N',T90,'EXPD1',T99,'COUL', FORT0398 . T109,'C1',T118,'C2',T125,'EXPD2') FORT0399 IF(NH.EQ.0) GO TO 72 FORT0400 J=1 FORT0401 DO 76 I=1,NH FORT0402 76 WRITE(6,53) HYDROG,I,X(I),Y(I),Z(I),J,PEEP,COULH FORT0403 53 FORMAT(T4,A2,I3,3F12.5,3(I3,F8.4,F9.4),2F9.5,F8.4) FORT0404 C C ** THIS CHUNK WRITES ANY HYDROGEN ATOM COORDINATES TO DISK FILE 13 C THAT ARE DEFINED SPECIFICALLY AS HYDROGENS (I.E., NOT DEFINED AS C HEAVY ATOMS). THE HYDROGEN AS A HEAVY ATOM CASE COMES LATER. C IT ALSO WRITES THE NUMBER OF ATOMS AND THE NUMBER OF VALENCE C ORBITALS (ALSO MOLECULAR ORBITALS) FOR USE IN SETTING UP THE C PLOTTING FILES. JJN 8-28-90 C IF (NH.EQ.0) GO TO 72 WRITE(13,9998) NATOM,NDIM,NELEC,KA,JOHN 9998 FORMAT(I3,2I4,4I3) DO 9991 I=1,NH 9991 WRITE(13,9992) X(I),Y(I),Z(I) 9992 FORMAT(' 1',3F12.6) C C ** C 72 CONTINUE FORT0405 DO 151 I=1,NA FORT0406 KEYI=AC(I) FORT0407 INH=I+NH FORT0408 IF(NP(KEYI).NE.0) GO TO 152 FORT0409 WRITE(6,53) SYMBOL(KEYI),INH,X(INH),Y(INH),Z(INH),NS(KEYI), FORT0410 . EXPS(KEYI),COULS(KEYI) FORT0411 GO TO 151 FORT0412 152 IF(ND(KEYI).NE.0) GO TO 153 FORT0413 WRITE(6,53) SYMBOL(KEYI),INH,X(INH),Y(INH),Z(INH),NS(KEYI), FORT0414 . EXPS(KEYI),COULS(KEYI),NP(KEYI),EXPP(KEYI),COULP(KEYI) FORT0415 GO TO 151 FORT0416 153 IF(C2(KEYI).NE.0.0D0) GO TO 154 FORT0417 WRITE(6,53) SYMBOL(KEYI),INH,X(INH),Y(INH),Z(INH),NS(KEYI), FORT0418 . EXPS(KEYI),COULS(KEYI),NP(KEYI),EXPP(KEYI),COULP(KEYI), FORT0419 . ND(KEYI),EXPD(KEYI),COULD(KEYI) FORT0420 GO TO 151 FORT0421 154 WRITE(6,53) SYMBOL(KEYI),INH,X(INH),Y(INH),Z(INH),NS(KEYI), FORT0422 . EXPS(KEYI),COULS(KEYI),NP(KEYI),EXPP(KEYI),COULP(KEYI), FORT0423 . ND(KEYI),EXPD(KEYI),COULD(KEYI),C1(KEYI),C2(KEYI),EXPD2(KEYI)FORT0424 151 CONTINUE FORT0425 C C THIS CHUNK WRITES THE COORDINATES AND ATOMIC NUMBER C (FOR ELEMENTS LESS THAN OR EQUAL TO 89) TO DISK FILE C 13 FOR NON-HYDROGEN ATOMS. JJN 9-8-90 C C THE FIRST ITEMS WRITTEN OUT ARE THE NUMBER OF ATOMS AND THE C NUMBER OF VALENCE ORBITALS (I.E., MOLECULAR ORBITALS) FOR USE C IN CONSTRUCTING THE PLOTTING ROUTINE INPUT FILES. C IF (NH.NE.0) GO TO 10051 10049 WRITE(13,10050) NATOM,NDIM,NELEC,KA,JOHN 10050 FORMAT(I3,2I4,4I3) 10051 DO 9996 I=1,NA KEYI=AC(I) INH=I+NH DO 9995 J=1,89 IF (SYMBOL(KEYI).EQ.ATS(J)) SYMBL(KEYI)=J 9995 CONTINUE WRITE(13,9994) SYMBL(KEYI),X(INH),Y(INH),Z(INH) 9994 FORMAT(I2,3F12.6) 9996 CONTINUE C C C WRITE(6,160) KA,IPRINT,IPUNCH,CON FORT0426 160 FORMAT(///,T10,'CHARGE =',I3,8X,'IPRINT =',I3,8X,'IPUNCH =', FORT0427 . I3,8X,'HUCKEL CONSTANT =',F7.3) FORT0428 C FORT0429 C PRINT OUT ITERATION PARAMETERS. FORT0430 C FORT0431 IF(.NOT.ITERAT) GO TO 80 FORT0432 WRITE(6,81) DAMP1,DAMP2,DAMP3,LAMPRI,MAXCYC,PRTCYC, FORT0433 . SENSE,DELTAC FORT0434 81 FORMAT(/,T10,'DAMP1 =',F6.3,6X,'DAMP2 =',F6.3,6X,'DAMP3 =', FORT0435 . F6.3,6X,'LAMPRI =',F6.3,//,T10,'MAXCYC =',I3,8X,'PRTCYC =', FORT0436 . I3,8X,'SENSE =',F6.3,6X,'DELTAC =',F10.7) FORT0437 C FORT0438 C PRINT OUT VSIE PARAMETERS. FORT0439 C FORT0440 IF(METH.LT.2) GO TO 80 FORT0441 WRITE(6,82) FORT0442 82 FORMAT(///,' VSIE PARAMETERS',//,T10,'ATOM',T26,'A',T39,'B', FORT0443 . T52,'C') FORT0444 NUSER2=MXUSR2-NUSER DO 83 I=1,NUSER2 FORT0446 J=MXUSR2-I IF(.NOT.ITABLE(J)) GO TO 83 FORT0448 WRITE(6,84) SYMBOL(J),AS1(I),BS1(I),CS1(I) FORT0449 84 FORMAT(/,T11,A2,4X,3F13.5) FORT0450 IF(NP(J).EQ.0) GO TO 83 FORT0451 IF(ND(J).NE.0) GO TO 85 FORT0452 WRITE(6,86) AP1(I),BP1(I),CP1(I) FORT0453 86 FORMAT(T17,3F13.5) FORT0454 GO TO 83 FORT0455 85 IF(NCON.EQ.3) GO TO 87 FORT0456 WRITE(6,86) AP1(I),BP1(I),CP1(I),AD1(I),BD1(I),CD1(I) FORT0457 GO TO 83 FORT0458 87 WRITE(6,86) AS2(I),BS2(I),CS2(I),AS3(I),BS3(I),CS3(I),AP1(I),FORT0459 . BP1(I),CP1(I),AP2(I),BP2(I),CP2(I),AP3(I),BP3(I),CP3(I), FORT0460 . AD1(I),BD1(I),CD1(I),AD2(I),BD2(I),CD2(I),AD3(I),BD3(I), FORT0461 . CD3(I) FORT0462 83 CONTINUE FORT0463 80 WRITE(6,99) FORT0464 99 FORMAT(///) FORT0465 RETURN FORT0466 END FORT0467 SUBROUTINE MOVLAP(H,S,MAD,C,SP,PD,MAXS,MAXP,MAXD,COUL0,SORB, FORT0468 . IOCC,HDG, NDIM, ND1 ,NC,NATOM,NTYPE,NHDG) FORT0469 C FORT0470 C SUBROUTINE TO CALCULATE INTERATOMIC DISTANCES, OVERLAP FORT0471 C INTEGRALS, AND MADELUNG PARAMETERS. FORT0472 C FORT0473 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MXUSER=230) PARAMETER (MXUSR2=231) DIMENSION H(NDIM,NDIM),S(NDIM,NDIM),MAD(NTYPE,NTYPE),C(NC), . SP(NDIM),PD(NDIM),MAXS(NDIM),MAXP(NDIM),MAXD(NDIM), . COUL0(MXUSER),SORB(NATOM),IOCC(NDIM),HDG(NHDG) REAL*8 MAD FORT0478 LOGICAL*4 SP,PD FORT0479 INTEGER*4 COUL0,SORB C FORT0482 C COUL0 DIMENSIONED AT 20 FOR EASY READING DURING PROCESSING FORT0483 C OF DELETION INPUT. DELETIONS DONE IN SUBROUTINE DELETS. FORT0484 C FORT0485 PARAMETER (MAXATM=500) PARAMETER (BB=250) COMMON/TITLE/AB(10) FORT0486 COMMON/CNTRL/CON,PEEP,COULH,NH,NA,NATM,KA,NELEC,METH, FORT0487 . IPRINT,IPUNCH,LA,LB,L3,L4,L5,ONEMAT,ITERAT FORT0488 LOGICAL*1 LA,LB,L3,L4,L5,ONEMAT,ITERAT FORT0489 COMMON/OUT/PRT(20),PUN(20) FORT0490 LOGICAL*1 PRT,PUN FORT0491 COMMON/ATOM/KEY(BB),SYMBOL(BB),VELEC(BB),NS(BB),NP(BB), . ND(BB),EXPS(BB),EXPP(BB),EXPD(BB),EXPD2(BB),C1(BB), . C2(BB),COULS(BB),COULP(BB),COULD(BB),X(MAXATM),Y(MAXATM), . Z(MAXATM) INTEGER*2 SYMBOL,KEY,VELEC FORT0495 COMMON/STARS/STAR,STAR2 FORT0496 INTEGER*2 STAR,STAR2 FORT0497 COMMON /LOCLAP/ SK1,SK2,R,L1,L2,M,N1,N2,MAX FORT0498 REAL*4 MADS(MXUSER),MADP(MXUSER),MADD(MXUSER) EQUIVALENCE (MADS(1),NS(MXUSR2)),(MADP(1),NP(MXUSR2)), . (MADD(1),ND(MXUSR2)) DIMENSION PTR(9),DTR(25) FORT0502 DIMENSION A(MXUSER),B(MXUSER),A1(MXUSER),B1(MXUSER) LOGICAL*1 JGO FORT0504 EQUIVALENCE (PTR(3),CA),(PTR(8),CB) FORT0505 DATA SQRT3/1.7320508075688770/ FORT0506 DATA AUI/1.889644746D0/ FORT0507 DATA PTR(9)/0.D0/,DTR(12)/0.D0/,DTR(22)/0.D0/ FORT0508 NH1=NH+1 FORT0509 C FORT0510 C HYDROGEN-HYDROGEN OVERLAPS. FORT0511 C FORT0512 IF(NH.LE.1) GO TO 106 FORT0513 DO 107 I=2,NH FORT0514 IM1=I-1 FORT0515 DO 107 J=1,IM1 FORT0516 R=DSQRT((X(I)-X(J))**2+(Y(I)-Y(J))**2+(Z(I)-Z(J))**2) FORT0517 C FORT0518 C STORE OVERLAPS IN UPPER RIGHT TRIANGLE OF S(I,J). PUT FORT0519 C DISTANCES IN LOWER RIGHT TRIANGLE. FORT0520 C FORT0521 S(I,J)=R FORT0522 R=R*PEEP*AUI FORT0523 IF(R.GT.50) GO TO 105 FORT0524 SIGMA=(1.D0+R*(1.D0+R/3.D0))/DEXP(R) FORT0525 GO TO 107 FORT0526 105 SIGMA=0.D0 FORT0527 107 S(J,I)=SIGMA FORT0528 C FORT0529 C HEAVY ATOM-HEAVY ATOM OVERLAPS. LOCAL COORDINATE SYSTEM FORT0530 C CENTERED ON ATOM J. FILL IN UPPER RIGHT TRIANGLE OF S(I,J). FORT0531 C FORT0532 106 IORB=NH1 FORT0533 DO 130 I=1,NA FORT0534 INH=I+NH FORT0535 IM1=I-1 FORT0536 KEYI=KEY(I) FORT0537 MAXD(I)=ND(KEYI) FORT0538 MAXP(I)=NP(KEYI) FORT0539 MAXS(I)=NS(KEYI) FORT0540 SP(I)=EXPS(KEYI) .EQ. EXPP(KEYI) FORT0541 PD(I)=EXPP(KEYI) .EQ. EXPD(KEYI) FORT0542 IF(PD(I)) MAXP(I)=MAX0(MAXP(I),MAXD(I)) FORT0543 IF(SP(I)) MAXS(I)=MAX0(NS(KEYI),MAXP(I)) FORT0544 SORB(I)=IORB FORT0545 C FORT0546 C SORB(I) CONTAINS A POINTER TO THE S ORBITAL ON ATOM I. FORT0547 C FORT0548 IORBS=IORB FORT0549 IORB=IORB+4 FORT0550 IF(NP(KEYI).EQ.0) IORB=IORB-3 FORT0551 IF(ND(KEYI).NE.0) IORB=IORB+5 FORT0552 IF(NP(KEYI).EQ.0) GO TO 298 FORT0553 JD=IORB-1 FORT0554 JD1=JD-1 FORT0555 DO 280 JC = IORBS,JD1 FORT0556 ID=JC+1 FORT0557 DO 280 IC=ID,JD FORT0558 280 S(JC,IC)=0.D0 FORT0559 298 CONTINUE FORT0560 IF(I.EQ.1) GO TO 300 FORT0561 DO 131 J=1,IM1 FORT0562 KEYJ=KEY(J) FORT0563 JNH=J+NH FORT0564 JORBS=SORB(J) FORT0565 DELX=X(INH)-X(JNH) FORT0566 DELY=Y(INH)-Y(JNH) FORT0567 DELZ=Z(INH)-Z(JNH) FORT0568 RT2=DELX**2+DELY**2 FORT0569 R=DSQRT(RT2+DELZ**2) FORT0570 S(INH,JNH)=R FORT0571 C FORT0572 C STORE DISTANCES IN LOWER LEFT TRIANGLE OF S(I,J). FORT0573 C FORT0574 IF(R.GT.0.0D0) GO TO 102 FORT0575 ID=IORB-1 FORT0576 JD=SORB(J+1)-1 FORT0577 DO 103 IC=IORBS,ID FORT0578 DO 103 JC=JORBS,JD FORT0579 103 S(JC,IC)=0.0D0 FORT0580 GO TO 131 FORT0581 102 IF(RT2.GT.1.E-10) GO TO 135 FORT0582 CB=1.D0 FORT0583 SB=0.D0 FORT0584 SA=0.D0 FORT0585 GOTO 136 FORT0586 135 T=DSQRT(RT2) FORT0587 CB=DELX/T FORT0588 SB=DELY/T FORT0589 SA=T/R FORT0590 136 CA=DELZ/R FORT0591 C FORT0592 C THE TRANSFORMATION MATRICES ARE CALCULATED EXPLICITLY. FORT0593 C PTR IS THE MATRIX FOR PROJECTING THE X,Y,Z ORBITALS FORT0594 C ONTO THE LOCAL SYSTEM. THE ELEMENTS ARE ORDERED SO THAT FIRST FORT0595 C X THEN Y THEN Z IS PROJECTED ONTO THE Z' AXIS (SIGMA). FORT0596 C THEN THE 3 ARE PROJECTED ONTO THE X' AXIS AND THEN THE Y' (PI). FORT0597 C THE D ORBITALS ARE HANDLED SIMILARLY. THE ORDER OF PROJECTION FORT0598 C IS X2-Y2,Z2,XY,XZ,YZ FIRST ONTO Z2'(SIGMA)AND THEN ONTO XZ' AND FORT0599 C YZ'(PI). FINALLY THE 5 ORBITALS ARE PROJECTED ONTO X'2-Y'2 AND FORT0600 C THEN XY' (DELTA). FORT0601 C FORT0602 C THOSE PTR AND DTR WHICH ARE ZERO ARE INITIALIZED IN A DATA STATE- FORT0603 C MENT. CA AND CB HAVE BEEN EQUIVALENCED TO PTR(3) AND PTR(8) FORT0604 C RESPECTIVELY TO SAVE TIME. FORT0605 C FORT0606 PTR(1)= SA*CB FORT0607 PTR(2)= SA*SB FORT0608 C ... PTR(3)= CA FORT0609 PTR(4)= CA*CB FORT0610 PTR(5)= CA*SB FORT0611 PTR(6)= -SA FORT0612 PTR(7)= -SB FORT0613 C ... PTR(8)= CB FORT0614 C ... PTR(9)= 0.D0 FORT0615 IF(ND(KEYI)+ND(KEYJ).EQ.0) GO TO 180 FORT0616 CA2=CA**2 FORT0617 SA2=SA*SA FORT0618 CB2=CB*CB FORT0619 SB2=SB*SB FORT0620 CBSB= CB*SB FORT0621 CASA= CA*SA FORT0622 CB2SB2= CB2-SB2 FORT0623 DTR(1)= SQRT3*.5D0*SA2*CB2SB2 FORT0624 DTR(2)= 1.D0-1.5D0*SA2 FORT0625 DTR(3)= SQRT3*CBSB*SA2 FORT0626 DTR(4)= SQRT3*CASA*CB FORT0627 DTR(5)= SQRT3*CASA*SB FORT0628 DTR(6)= CASA*CB2SB2 FORT0629 DTR(7)= -SQRT3*CASA FORT0630 DTR(8)= 2.D0*CASA*CBSB FORT0631 DTR(9)= CB*(CA2-SA2) FORT0632 DTR(10)= SB*(CA2-SA2) FORT0633 DTR(11)= -2.D0*SA*CBSB FORT0634 C ... DTR(12)= 0.D0 FORT0635 DTR(13)= SA* CB2SB2 FORT0636 DTR(14)= -PTR(5) FORT0637 DTR(15)= PTR(4) FORT0638 IF(ND(KEYI)*ND(KEYJ).EQ.0) GO TO 180 FORT0639 DTR(16)=.5D0*(1.D0+CA2)*CB2SB2 FORT0640 DTR(17)= .5D0*SQRT3*SA2 FORT0641 DTR(18)= CBSB*(1.D0+CA2) FORT0642 DTR(19)= -CASA*CB FORT0643 DTR(20)= -CASA*SB FORT0644 DTR(21)= -2.D0*CA*CBSB FORT0645 C ... DTR(22)= 0.D0 FORT0646 DTR(23)= CA*CB2SB2 FORT0647 DTR(24)= PTR(2) FORT0648 DTR(25)= -PTR(1) FORT0649 180 R=R*AUI FORT0650 C FORT0651 C (S(I)!S(J)). FORT0652 C FORT0653 N2=NS(KEYJ) FORT0654 N1=NS(KEYI) FORT0655 L2=0 FORT0656 L1=0 FORT0657 M=0 FORT0658 MAX=MAXS(I)+MAXS(J) FORT0659 SK1=EXPS(KEYI) FORT0660 SK2=EXPS(KEYJ) FORT0661 CALL ABFNS(A,B) FORT0662 CALL LOVLAP(SIGMA,A,B) FORT0663 S(JORBS,IORBS)=SIGMA FORT0664 C FORT0665 C IF THE S EXPONENT OF ATOM I EQUALS THE P EXPONENT WE NEED FORT0666 C NOT CALCULATE THE A AND B FUNCTIONS AGAIN. FORT0667 C FORT0668 C (P(I)!S(J)). FORT0669 C FORT0670 JGO=.FALSE. FORT0671 IF(KEYI.EQ.KEYJ) GO TO 126 FORT0672 IF((.NOT.SP(I)).OR.(NP(KEYI).EQ.0)) GO TO 126 FORT0673 220 N1=NP(KEYI) FORT0674 L1=1 FORT0675 CALL LOVLAP(SIGMA,A,B) FORT0676 SIGMA=-SIGMA FORT0677 DO 200 IC=1,3 FORT0678 200 S(JORBS,IORBS+IC)=PTR(IC)*SIGMA FORT0679 IF(PD(I)) GO TO 221 FORT0680 IF(JGO) GO TO 217 FORT0681 GO TO 137 FORT0682 C FORT0683 C (D(I)!S(J)) CONDITIONALLY AT FIRST CHANCE. FORT0684 C FORT0685 221 N1=ND(KEYI) FORT0686 L1=2 FORT0687 168 CALL LOVLAP(SIGMA,A,B) FORT0688 IF(C2(KEYI).EQ.0.D0) GO TO 167 FORT0689 SK1=EXPD2(KEYI) FORT0690 CALL ABFNS(A1,B1) FORT0691 CALL LOVLAP(PART2,A1,B1) FORT0692 SIGMA=C1(KEYI)*SIGMA+C2(KEYI)*PART2 FORT0693 SK1=EXPD(KEYI) FORT0694 167 ID=IORBS+3 FORT0695 DO 201 IC=1,5 FORT0696 201 S(JORBS,ID+IC)=DTR(IC)*SIGMA FORT0697 C FORT0698 C CALCULATE (D(I)!P(J)) IF CAN USE SAME A'S AND B'S. FORT0699 C FORT0700 IF(SP(J)) GO TO 222 FORT0701 IF(JGO) GO TO 228 FORT0702 GO TO 137 FORT0703 222 N2=NP(KEYJ) FORT0704 L2=1 FORT0705 M=0 FORT0706 CALL LOVLAP(SIGMA,A,B) FORT0707 M=1 FORT0708 CALL LOVLAP(PI,A,B) FORT0709 IF(C2(KEYI).EQ.0.D0) GO TO 1169 FORT0710 SK1=EXPD2(KEYI) FORT0711 CALL LOVLAP(PART2,A1,B1) FORT0712 PI=C1(KEYI)*PI+C2(KEYI)*PART2 FORT0713 M=0 FORT0714 CALL LOVLAP(PART2,A1,B1) FORT0715 SK1=EXPD(KEYI) FORT0716 SIGMA=C1(KEYI)*SIGMA+C2(KEYI)*PART2 FORT0717 1169 PI=-PI FORT0718 ID=IORBS+3 FORT0719 DO 195 JC=1,3 FORT0720 DO 195 IC=1,5 FORT0721 195 S(JORBS+JC,ID+IC)=PTR(JC)*DTR(IC)*SIGMA+(PTR(JC+3)*DTR(IC+5) FORT0722 . +PTR(JC+6)*DTR(IC+10))*PI FORT0723 IF(JGO) GO TO 131 FORT0724 C FORT0725 C NOW TEST FOR DUPLICATE EXPONENTS ON ATOM J. FORT0726 C HOWEVER DO CALCULATIONS ANYHOW. FORT0727 C FORT0728 137 N1=NS(KEYI) FORT0729 L1=0 FORT0730 C FORT0731 C (S(I)!P(J)). FORT0732 C FORT0733 126 IF(SP(J)) GO TO 138 FORT0734 IF(NP(KEYJ).EQ.0) GO TO 210 FORT0735 MAX=MAXS(I)+MAXP(J) FORT0736 SK2=EXPP(KEYJ) FORT0737 CALL ABFNS(A,B) FORT0738 138 N2=NP(KEYJ) FORT0739 L2=1 FORT0740 M=0 FORT0741 CALL LOVLAP(SIGMA,A,B) FORT0742 DO 202 IC=1,3 FORT0743 202 S(JORBS+IC,IORBS)=PTR(IC)*SIGMA FORT0744 IF(SP(I)) GO TO 156 FORT0745 JGO=.TRUE. FORT0746 IF(ND(KEYJ).NE.0) GO TO 149 FORT0747 C FORT0748 C BRANCH TO TEST FOR EXPP(J) .EQ. EXPD(J). CALCULATE (S!D) ANYHOW. FORT0749 C RETURN WILL BE MADE TO THE NEXT STATEMENT. FORT0750 C FORT0751 C (P(I)!P(J)) EXPP(I) EQ,NE EXPS(I). FORT0752 C FORT0753 GO TO 646 FORT0754 146 N2=NP(KEYJ) FORT0755 L2=1 FORT0756 SK2=EXPP(KEYJ) FORT0757 646 IF(NP(KEYI).EQ.0) GO TO 210 FORT0758 SK1=EXPP(KEYI) FORT0759 C FORT0760 C THESE STATEMENTS USED ONLY IF HAVE ALREADY CALCULATED (S(I)!D(J)) FORT0761 C WHICH MEANS THAT SP(I) IS FALSE. FORT0762 C FORT0763 MAX=MAXP(I)+MAXP(J) FORT0764 CALL ABFNS(A,B) FORT0765 156 N1=NP(KEYI) FORT0766 L1=1 FORT0767 148 M=0 FORT0768 CALL LOVLAP(SIGMA,A,B) FORT0769 SIGMA=-SIGMA FORT0770 M=1 FORT0771 CALL LOVLAP(PI,A,B) FORT0772 DO 204 JC=1,3 FORT0773 DO 204 IC=JC,3 FORT0774 S(JORBS+JC,IORBS+IC)=PTR(JC)*PTR(IC)*SIGMA + (PTR(JC+3)* FORT0775 . PTR(IC+3)+PTR(JC+6)*PTR(IC+6))*PI FORT0776 204 S(JORBS+IC,IORBS+JC)=S(JORBS+JC,IORBS+IC) FORT0777 147 IF(ND(KEYJ).EQ.0) GO TO 210 FORT0778 C FORT0779 C BRANCH AROUND (S(I)!D(J)) IF ALREADY DONE. FORT0780 C FORT0781 IF(JGO) GO TO 160 FORT0782 C FORT0783 C (S(I)!D(J)). FORT0784 C FORT0785 N1=NS(KEYI) FORT0786 L1=0 FORT0787 149 N2=ND(KEYJ) FORT0788 L2=2 FORT0789 IF(PD(J)) GO TO 142 FORT0790 SK2=EXPD(KEYJ) FORT0791 MAX=MAXS(I)+MAXD(J) FORT0792 CALL ABFNS(A,B) FORT0793 142 M=0 FORT0794 CALL LOVLAP(SIGMA,A,B) FORT0795 IF(C2(KEYJ).EQ.0.D0) GO TO 151 FORT0796 SK2=EXPD2(KEYJ) FORT0797 CALL ABFNS(A1,B1) FORT0798 CALL LOVLAP(PART2,A1,B1) FORT0799 SIGMA=C1(KEYJ)*SIGMA+C2(KEYJ)*PART2 FORT0800 SK2=EXPD(KEYJ) FORT0801 151 JD=JORBS+3 FORT0802 DO 205 IC=1,5 FORT0803 205 S(JD+IC,IORBS)=DTR(IC)*SIGMA FORT0804 150 IF(JGO) GO TO 146 FORT0805 C FORT0806 C SP(I) IS TRUE IF HERE SO BRANCH AS WE ALSO HAVE D ON ATOM J. FORT0807 C FORT0808 GO TO 170 FORT0809 160 JGO=.FALSE. FORT0810 C FORT0811 C (P(I)!D(J)). FORT0812 C FORT0813 N2=ND(KEYJ) FORT0814 L2=2 FORT0815 IF(PD(J)) GO TO 178 FORT0816 SK2=EXPD(KEYJ) FORT0817 MAX=MAXP(I)+MAXD(J) FORT0818 CALL ABFNS(A,B) FORT0819 178 IF(C2(KEYJ).EQ.0.D0) GO TO 170 FORT0820 SK2=EXPD2(KEYJ) FORT0821 CALL ABFNS(A1,B1) FORT0822 SK2=EXPD(KEYJ) FORT0823 170 N1=NP(KEYI) FORT0824 L1=1 FORT0825 M=0 FORT0826 CALL LOVLAP(SIGMA,A,B) FORT0827 M=1 FORT0828 CALL LOVLAP(PI,A,B) FORT0829 IF(C2(KEYJ).EQ.0.D0) GO TO 171 FORT0830 SK2=EXPD2(KEYJ) FORT0831 CALL LOVLAP(PART2,A1,B1) FORT0832 PI=C1(KEYJ)*PI+C2(KEYJ)*PART2 FORT0833 M=0 FORT0834 CALL LOVLAP(PART2,A1,B1) FORT0835 SIGMA=C1(KEYJ)*SIGMA+C2(KEYJ)*PART2 FORT0836 SK2=EXPD(KEYJ) FORT0837 171 SIGMA=-SIGMA FORT0838 DO 206 IC=1,3 FORT0839 DO 206 JC=1,5 FORT0840 206 S(JD+JC,IORBS+IC)=DTR(JC)*PTR(IC)*SIGMA+(DTR(JC+5)*PTR(IC+3) FORT0841 . +DTR(JC+10)*PTR(IC+6))*PI FORT0842 C FORT0843 C (D(I)!D(J)). FORT0844 C FORT0845 IF(ND(KEYI).EQ.0) GO TO 210 FORT0846 MAX=MAXD(I)+MAXD(J) FORT0847 IF(PD(I)) GO TO 208 FORT0848 SK1=EXPD(KEYI) FORT0849 CALL ABFNS(A,B) FORT0850 IF(C2(KEYJ).EQ.0.D0) GO TO 208 FORT0851 SK2=EXPD2(KEYJ) FORT0852 CALL ABFNS(A1,B1) FORT0853 SK2=EXPD(KEYJ) FORT0854 208 N1=ND(KEYI) FORT0855 L1=2 FORT0856 M=0 FORT0857 CALL LOVLAP(SIGMA,A,B) FORT0858 M=1 FORT0859 CALL LOVLAP(PI,A,B) FORT0860 M=2 FORT0861 CALL LOVLAP(DELTA,A,B) FORT0862 CC=C2(KEYI) FORT0863 IF(C2(KEYJ).EQ.0.D0) GO TO 173 FORT0864 CC=C1(KEYJ)*CC FORT0865 SK2=EXPD2(KEYJ) FORT0866 CALL LOVLAP(PART2,A1,B1) FORT0867 DELTA=C1(KEYJ)*DELTA+C2(KEYJ)*PART2 FORT0868 M=1 FORT0869 CALL LOVLAP(PART3,A1,B1) FORT0870 PI=C1(KEYJ)*PI+C2(KEYJ)*PART3 FORT0871 M=0 FORT0872 CALL LOVLAP(PART4,A1,B1) FORT0873 SIGMA=C1(KEYJ)*SIGMA+C2(KEYJ)*PART4 FORT0874 SK2=EXPD(KEYJ) FORT0875 M=2 FORT0876 173 IF(C2(KEYI).EQ.0.D0) GO TO 172 FORT0877 IF(KEYI.EQ.KEYJ) GO TO 176 FORT0878 SK1=EXPD2(KEYI) FORT0879 CALL ABFNS(A1,B1) FORT0880 CALL LOVLAP(PART2,A1,B1) FORT0881 M=1 FORT0882 CALL LOVLAP(PART3,A1,B1) FORT0883 M=0 CALL LOVLAP(PART4,A1,B1) FORT0884 176 SIGMA=C1(KEYI)*SIGMA+CC*PART4 FORT0885 PI =C1(KEYI)*PI+CC*PART3 FORT0886 DELTA=C1(KEYI)*DELTA+CC*PART2 FORT0887 IF(C2(KEYJ).EQ.0.D0) GO TO 172 FORT0888 SK1=EXPD2(KEYI) FORT0889 SK2=EXPD2(KEYJ) FORT0890 CALL ABFNS(A1,B1) FORT0891 M=0 FORT0892 CALL LOVLAP(PART2,A1,B1) FORT0893 CC=C2(KEYI)*C2(KEYJ) FORT0894 SIGMA=SIGMA+CC*PART2 FORT0895 M=1 FORT0896 CALL LOVLAP(PART2,A1,B1) FORT0897 PI=PI+CC*PART2 FORT0898 M=2 FORT0899 CALL LOVLAP(PART2,A1,B1) FORT0900 DELTA=DELTA+CC*PART2 FORT0901 172 PI=-PI FORT0902 JD=JORBS+3 FORT0903 DO 211 IC=1,5 FORT0904 ID=IORBS+3 FORT0905 DO 211 JC=1,5 FORT0906 S(JD+JC,ID+IC) = DTR(IC)*DTR(JC)*SIGMA+(DTR(IC+5)*DTR(JC+5) FORT0907 . +DTR(IC+10)*DTR(JC+10))*PI+(DTR(IC+15)*DTR(JC+15)+DTR(IC+20) FORT0908 . *DTR(JC+20))*DELTA FORT0909 211 S(JD+IC,ID+JC)=S(JD+JC,ID+IC) FORT0910 C FORT0911 C FILLING IN OTHER HALF OF OVERLAPS FOR (J!I) AS NEEDED. FORT0912 C FORT0913 210 IF(KEYI.EQ.KEYJ) GO TO 213 FORT0914 N2=NS(KEYJ) FORT0915 L2=0 FORT0916 SK2=EXPS(KEYJ) FORT0917 M=0 FORT0918 JGO=.TRUE. FORT0919 IF(NP(KEYI).EQ.0) GO TO 131 FORT0920 IF(SP(I)) GO TO 215 FORT0921 MAX=MAXP(I)+MAXS(J) FORT0922 SK1=EXPP(KEYI) FORT0923 CALL ABFNS(A,B) FORT0924 GO TO 220 FORT0925 215 IF(PD(I)) GO TO 227 FORT0926 217 IF(ND(KEYI).EQ.0) GO TO 131 FORT0927 MAX=MAXD(I)+MAXS(J) FORT0928 SK1=EXPD(KEYI) FORT0929 CALL ABFNS(A,B) FORT0930 GO TO 221 FORT0931 227 IF(SP(J)) GO TO 131 FORT0932 N1=ND(KEYI) FORT0933 L1=2 FORT0934 SK1=EXPD(KEYI) FORT0935 228 IF(NP(KEYJ).EQ.0) GO TO 131 FORT0936 SK2=EXPP(KEYJ) FORT0937 MAX=MAXD(I)+MAXP(J) FORT0938 CALL ABFNS(A,B) FORT0939 IF(C2(KEYI).EQ.0.D0) GO TO 222 FORT0940 SK1=EXPD2(KEYI) FORT0941 CALL ABFNS(A1,B1) FORT0942 SK1=EXPD(KEYI) FORT0943 GO TO 222 FORT0944 213 IF(NP(KEYI).EQ.0) GO TO 131 FORT0945 DO 237 IC=1,3 FORT0946 237 S(JORBS,IORBS+IC)=-S(JORBS+IC,IORBS) FORT0947 IF(ND(KEYI).EQ.0) GO TO 131 FORT0948 DO 238 IC=4,8 FORT0949 S(JORBS,IORBS+IC)=S(JORBS+IC,IORBS) FORT0950 DO 238 JC=1,3 FORT0951 238 S(JORBS+JC,IORBS+IC)=-S(JORBS+IC,IORBS+JC) FORT0952 131 CONTINUE FORT0953 300 IF(NH.EQ.0) GO TO 130 FORT0954 N2=1 FORT0955 L2=0 FORT0956 M=0 FORT0957 SK2=PEEP FORT0958 DO 301 J=1,NH FORT0959 DELX=X(J)-X(INH) FORT0960 DELY=Y(J)-Y(INH) FORT0961 DELZ=Z(J)-Z(INH) FORT0962 RT2=DELX**2+DELY**2 FORT0963 R=DSQRT(RT2+DELZ**2) FORT0964 C FORT0965 C STORE DISTANCES IN LOWER LEFT TRIANGLE OF S(I,J). FORT0966 C FORT0967 S(INH,J)=R FORT0968 IF(RT2.GT.1.D-10) GO TO 303 FORT0969 CB=1.D0 FORT0970 SB=0.D0 FORT0971 SA=0.D0 FORT0972 GO TO 302 FORT0973 303 T=DSQRT(RT2) FORT0974 CB=DELX/T FORT0975 SB=DELY/T FORT0976 SA=T/R FORT0977 302 CA=DELZ/R FORT0978 R=R*AUI FORT0979 C FORT0980 C H(J)!S(I)). FORT0981 C FORT0982 N1=NS(KEYI) FORT0983 L1=0 FORT0984 MAX=1+MAXS(I) FORT0985 SK1=EXPS(KEYI) FORT0986 CALL ABFNS(A,B) FORT0987 CALL LOVLAP(SIGMA,A,B) FORT0988 S(J,IORBS)=SIGMA FORT0989 IF(NP(KEYI).EQ.0) GO TO 301 FORT0990 IF(SP(I)) GO TO 304 FORT0991 SK1=EXPP(KEYI) FORT0992 MAX=1+MAXP(I) FORT0993 CALL ABFNS(A,B) FORT0994 304 N1=NP(KEYI) FORT0995 L1=1 FORT0996 CALL LOVLAP(SIGMA,A,B) FORT0997 S(J,IORBS+3)=CA*SIGMA FORT0998 SIGMA=SIGMA*SA FORT0999 S(J,IORBS+2)=SB*SIGMA FORT1000 S(J,IORBS+1)=CB*SIGMA FORT1001 IF(ND(KEYI).EQ.0) GO TO 301 FORT1002 IF(PD(I)) GO TO 305 FORT1003 SK1=EXPD(KEYI) FORT1004 MAX=1+ND(KEYI) FORT1005 CALL ABFNS(A,B) FORT1006 305 N1=ND(KEYI) FORT1007 L1=2 FORT1008 CALL LOVLAP(SIGMA,A,B) FORT1009 IF(C2(KEYI).EQ.0.D0) GO TO 181 FORT1010 SK1=EXPD2(KEYI) FORT1011 CALL ABFNS(A1,B1) FORT1012 CALL LOVLAP(PART2,A1,B1) FORT1013 SK1=EXPD(KEYI) FORT1014 SIGMA=C1(KEYI)*SIGMA+C2(KEYI)*PART2 FORT1015 181 CONTINUE FORT1016 S(J,IORBS+5)=(1.D0-1.5D0*SA*SA)*SIGMA FORT1017 SIGMA=SIGMA*SQRT3*SA FORT1018 S(J,IORBS+4)=.5D0*SA*(CB*CB-SB*SB)*SIGMA FORT1019 S(J,IORBS+6)=CB*SB*SA*SIGMA FORT1020 SIGMA=SIGMA*CA FORT1021 S(J,IORBS+7)=CB*SIGMA FORT1022 S(J,IORBS+8)=SB*SIGMA FORT1023 301 CONTINUE FORT1024 130 CONTINUE FORT1025 C FORT1026 C CALL DELETS TO SET CERTAIN OVERLAP INTEGRALS = 0. FORT1027 C FORT1028 IF(.NOT.LB) GO TO 835 FORT1029 CALL DELETS(S,COUL0,SORB,NDIM) FORT1030 WRITE(6,2010) FORT1031 2010 FORMAT(///) FORT1032 C FORT1033 C CALCULATE INTERATOMIC MADELUNG PARAMETERS. FORT1034 C FORT1035 835 IF(METH.LT.3) GO TO 450 FORT1036 IC=1 FORT1037 DO 401 I=1,NA FORT1038 KEYI=KEY(I) FORT1039 RS=0.0D0 FORT1040 ID=IC FORT1041 MAD(ID,ID)=DBLE(MADS(MXUSR2-KEYI)) IF(NP(KEYI).EQ.0) GO TO 402 FORT1043 ID=IC+1 FORT1044 MAD(ID,ID)=DBLE(MADP(MXUSR2-KEYI)) IF(ND(KEYI).EQ.0) GO TO 403 FORT1046 ID=IC+2 FORT1047 MAD(ID,ID)=DBLE(MADD(MXUSR2-KEYI)) 403 M=IC+1 FORT1049 DO 404 K=M,ID FORT1050 K1=K-1 FORT1051 DO 404 L=IC,K1 FORT1052 CA=MAD(K,K) FORT1053 CB=MAD(L,L) FORT1054 SA=VALMAD(CA,CB,RS) FORT1055 MAD(K,L)=SA FORT1056 404 MAD(L,K)=SA FORT1057 402 IF(I.EQ.1) GO TO 401 FORT1058 IM1=I-1 FORT1059 JC=1 FORT1060 DO 406 J=1,IM1 FORT1061 KEYJ=KEY(J) FORT1062 RS=S(I,J)*AUI/27.21D0 FORT1063 JD=JC FORT1064 IF(NP(KEYJ).NE.0) JD=JC+1 FORT1065 IF(ND(KEYJ).NE.0) JD=JC+2 FORT1066 DO 407 K=IC,ID FORT1067 CA=MAD(K,K) FORT1068 DO 407 L=JC,JD FORT1069 CB=MAD(L,L) FORT1070 SA=VALMAD(CA,CB,RS) FORT1071 MAD(K,L)=SA FORT1072 407 MAD(L,K)=SA FORT1073 406 JC=JD+1 FORT1074 401 IC=ID+1 FORT1075 C FORT1076 C SET UP DISTANCE MATRIX FOR PRINTING. FORT1077 C STUFF ELEMENTS OF S INTO C TO GET THEM OUT OF THE WAY. FORT1078 C FORT1079 450 ISUB=1 FORT1080 C FORT1081 C ZERO DISTANCE ALONG DIAGONAL. FORT1082 C FORT1083 S(1,1)=0.D0 FORT1084 DO 1010 I=2,NATOM FORT1085 S(I,I)=0.D0 FORT1086 IM1=I-1 FORT1087 DO 1005 J=1,IM1 FORT1088 C(ISUB)=S(J,I) FORT1089 ISUB=ISUB+1 FORT1090 1005 S(J,I)=S(I,J) FORT1091 1010 CONTINUE FORT1092 IF(PRT(3)) GO TO 2004 FORT1093 WRITE(6,2000) FORT1094 2000 FORMAT('DISTANCE MATRIX') FORT1095 CALL PEGLEG(S,NATOM,NDIM) FORT1096 2004 IF(PUN(3)) WRITE(7,2050) ((S(I,J),I=1,NATOM),J=1,NATOM) FORT1097 2050 FORMAT(8F9.6) FORT1098 C FORT1099 C SET UP OVERLAP MATRIX FOR PRINTING. FORT1100 C REPLACE ELEMENTS IN OVERLAP MATRIX FROM C. FORT1101 C FORT1102 1015 S(1,1)=1.D0 FORT1103 ISUB=1 FORT1104 DO 1025 I=2,NDIM FORT1105 S(I,I)=1.D0 FORT1106 IM1=I-1 FORT1107 DO 1020 J=1,IM1 FORT1108 IF(I.GT.NATOM) GO TO 1020 FORT1109 S(J,I)=C(ISUB) FORT1110 ISUB=ISUB+1 FORT1111 1020 S(I,J)=S(J,I) FORT1112 1025 CONTINUE FORT1113 IF(PRT(4)) GO TO 2005 FORT1114 WRITE(6,2001) FORT1115 2001 FORMAT('OVERLAP MATRIX') FORT1116 CALL PEGLEG(S,NDIM,NDIM) FORT1117 2005 IF(PUN(4)) WRITE(7,2050) S FORT1118 C FORT1119 C PRINT OUT MADELUNG PARAMETERS. FORT1120 C FORT1121 IF(METH.LT.3) GO TO 460 FORT1122 IF(PRT(5)) GO TO 2006 FORT1123 WRITE(6,2002) FORT1124 2002 FORMAT('MADELUNG PARAMETERS') FORT1125 CALL PEGLEG(MAD,NTYPE,NTYPE) FORT1126 2006 IF(PUN(5)) WRITE(7,2050) MAD FORT1127 460 RETURN FORT1128 END FORT1129 SUBROUTINE DELETS(S,COUL0,SORB,NDIM) FORT1130 C FORT1131 C SUBROUTINE FOR SETTING CERTAIN OVERLAP INTEGRALS EQUAL TO ZERO. FORT1132 C FORT1133 IMPLICIT REAL*8(A-H,O-Z) FORT1134 PARAMETER (MXUSER=230) PARAMETER (MXUSR2=231) DIMENSION S(NDIM,NDIM),COUL0(MXUSER),SORB(NDIM) INTEGER*4 COUL0 FORT1136 INTEGER*2 SORB FORT1137 COMMON/CNTRL/CON,PEEP,COULH,NH,NA,NATM,KA,NELEC,METH,IPRINT, FORT1138 1 IPUNCH,L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT1139 LOGICAL*1 L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT1140 COMMON/OUT/PRT(20),PUN(20) FORT1141 LOGICAL*1 PRT,PUN FORT1142 LOGICAL*1 IERR FORT1143 DATA ORBTL/' ORBITAL'/,ATMPR/' ATOM '/ FORT1144 SORB(NA+1)=NDIM+1 FORT1145 IERR=.FALSE. FORT1146 C FORT1147 C READ IN NUMBERS INDICATING WHICH OVERLAP INTEGRALS ARE TO BE SET FORT1148 C TO ZERO. A POSITIVE NUMBER REFERS TO AN ATOM, A NEGATIVE ONE TO FORT1149 C AN ORBITAL. FORT1150 C FORT1151 10 READ(5,1000,END=300) COUL0 FORT1152 1000 FORMAT(20I3) FORT1153 DO 200 IDEL=1,19,2 FORT1154 I=COUL0(IDEL) FORT1155 J=COUL0(IDEL+1) FORT1156 C FORT1157 C TERMINATE ON ENCOUNTERING A ZERO. FORT1158 C FORT1159 IF(I.EQ.0.OR.J.EQ.0) GO TO 400 FORT1160 IF(IERR) GO TO 200 FORT1161 IABSV=IABS(I) FORT1162 JABSV=IABS(J) FORT1163 IF(I.GT.NH) GO TO 20 FORT1164 C FORT1165 C I REFERS TO ORBITAL (NEGATIVE) OR H ATOM (POSITIVE,LE NH). FORT1166 C FORT1167 ILOW=IABSV FORT1168 IHIGH=IABSV FORT1169 C FORT1170 C ERROR IF I OUT OF RANGE OF ORBITAL NUMBERS. FORT1171 C FORT1172 IF(IABSV.GT.NDIM) GO TO 160 FORT1173 GO TO 30 FORT1174 C FORT1175 C ERROR IF I OUT OF RANGE OF ATOMS. FORT1176 C FORT1177 20 IF(I.GT.NATM) GO TO 160 FORT1178 ILOW=SORB(IABSV-NH) FORT1179 IHIGH=SORB(IABSV-NH+1)-1 FORT1180 30 IF(J.GT.NH) GO TO 40 FORT1181 JLOW=JABSV FORT1182 JHIGH=JABSV FORT1183 C FORT1184 C CHECK TO SEE IF J IS IN RANGE. FORT1185 C FORT1186 IF(JABSV.GT.NDIM) GO TO 160 FORT1187 GO TO 50 FORT1188 40 IF(J.GT.NATM) GO TO 160 FORT1189 JLOW=SORB(JABSV-NH) FORT1190 JHIGH=SORB(JABSV-NH+1)-1 FORT1191 50 X1=ATMPR FORT1192 IF(I.LT.0) X1=ORBTL FORT1193 X2=ATMPR FORT1194 IF(J.LT.0) X2=ORBTL FORT1195 IF(.NOT.PRT(2)) WRITE(6,1002) X1,IABSV,X2,JABSV FORT1196 1002 FORMAT('0ALL S(I,J) SET TO ZERO BETWEEN',A8,I4,' AND',A8,I4,'.') FORT1197 C FORT1198 C J MUST BE LESS THAN OR EQUAL TO I SINCE WE ONLY HAVE A HALF FORT1199 C MATRIX AT THIS POINT. FORT1200 C FORT1201 IF(JHIGH.LE.IHIGH) GO TO 60 FORT1202 I=JHIGH FORT1203 JHIGH=IHIGH FORT1204 IHIGH=I FORT1205 I=JLOW FORT1206 JLOW=ILOW FORT1207 ILOW=I FORT1208 60 DO 100 I=ILOW,IHIGH FORT1209 DO 100 J=JLOW,JHIGH FORT1210 100 S(J,I)=0.D0 FORT1211 GO TO 200 FORT1212 160 IERR=.TRUE. FORT1213 WRITE(6,1003) COUL0 FORT1214 1003 FORMAT('0NUMBER OUT OF RANGE IN FOLLOWING DELETION CARD'/'0',20I5/ *'0NO FURTHER DELETIONS, BUT SCANNING FOR ZERO TO TERMINATE CARD RE *ADING.') FORT1217 200 CONTINUE FORT1218 GO TO 10 FORT1219 300 WRITE(6,1004) FORT1220 1004 FORMAT('0END OF FILE IN DELETION CARDS, NO FURTHER DELETIONS.') FORT1221 400 RETURN FORT1222 END FORT1223 DOUBLE PRECISION FUNCTION VALMAD(A,B,R) FORT1224 C FORT1225 C FUNCTION ROUTINE FOR CALCULATING MADELUNG PARAMETERS. FORT1226 C FORT1227 IMPLICIT REAL*8(A-H,O-Z) FORT1228 IF(A.LT.0.01D0.OR.B.LT.0.01D0) GO TO 1 FORT1229 AB=(A+B)/(2.0D0*A*B) FORT1230 VALMAD=1.0D0/DSQRT(R*R+AB*AB) FORT1231 GO TO 2 FORT1232 1 VALMAD=0.0D0 FORT1233 IF(R.GT.0.001D0) VALMAD=1.0D0/R FORT1234 2 RETURN FORT1235 END FORT1236 SUBROUTINE PEGLEG(A,N,NL) FORT1237 C FORT1238 C SUBROUTINE TO PRINT OUT MATRICES IN READABLE FORMAT. FORT1239 C FORT1240 IMPLICIT REAL*8(A-H,O-Z) FORT1241 DIMENSION A(NL,NL) FORT1242 NROW=N FORT1243 NCOL=N FORT1244 GO TO 10 FORT1245 ENTRY OUTMAT(A,NL,NR,NC) FORT1246 NROW=NR FORT1247 NCOL=NC FORT1248 10 KITE=0 FORT1249 20 LOW=KITE+1 FORT1250 KITE=KITE+14 FORT1251 IF(KITE.GT.NCOL) KITE=NCOL FORT1252 WRITE(6,1000) (I,I=LOW,KITE) FORT1253 1000 FORMAT(/5X,14I8,//) FORT1254 DO 30 I=1,NROW FORT1255 30 WRITE(6,1001) I,(A(I,J),J=LOW,KITE) FORT1256 1001 FORMAT(I5,2X,14F8.4) FORT1257 IF(KITE.LT.NCOL) GO TO 20 FORT1258 RETURN FORT1259 END FORT1260 SUBROUTINE PEGLEG2(AA,NN,NLL) C C ROUTINE FOR WRITING OUT THE EIGENVECTORS TO DISK FILE 13 FOR C PLOTTING ROUTINES. JJN 8-8-90 C IMPLICIT REAL*8(A-H,O-Z) DIMENSION AA(NLL,NLL) NROWW=NN NCOLL=NN WRITE(13,9998) ((AA(J,I),J=1,NROWW), I=1,NCOLL) 9998 FORMAT(8F10.6) RETURN END C C THIS CONCLUDES THE EIGENVECTOR WRITE ROUTINE TO DISK FILE 13. C SUBROUTINE HUCKEL(H,S,MAD,C,SP,PD,MAXS,MAXP,MAXD,COUL0,SORB,IOCC, FORT1261 1 HDG, NDIM, ND1 ,NC,NATOM,NTYPE,NHDG) FORT1262 C FORT1263 C SUBROUTINE TO 1) DETERMINE ORBITAL OCCUPATION NUMBERS, AND 2) SETUP FORT1264 C HUCKEL MATRIX. FORT1265 C FORT1266 IMPLICIT DOUBLE PRECISION (A-H,O-Z) FORT1267 DIMENSION H(NDIM,NDIM),S(NDIM,NDIM),MAD(NTYPE,NTYPE),C(NC), FORT1268 1 SP(NDIM),PD(NDIM),MAXS(NDIM),MAXP(NDIM),MAXD(NDIM),COUL0(NATOM), FORT1269 2 SORB(NDIM),IOCC(NDIM),HDG(NHDG) FORT1270 REAL*8 MAD FORT1271 REAL*4 IOCC FORT1272 PARAMETER (MAXATM=500) PARAMETER (BB=250) PARAMETER (MXUSER=230) PARAMETER (MXUSR2=231) COMMON/TITLE/AB(10) FORT1273 COMMON/CNTRL/CON,PEEP,COULH,NH,NA,NATM,KA,NELEC,METH,IPRINT, FORT1274 1 IPUNCH,L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT1275 LOGICAL*1 L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT1276 COMMON/OUT/PRT(20),PUN(20) FORT1277 LOGICAL*1 PRT,PUN FORT1278 COMMON/ATOM/KEY(BB),SYMBOL(BB),VELEC(BB),NS(BB),NP(BB),ND(BB), 1 EXPS(BB),EXPP(BB),EXPD(BB),EXPD2(BB),C1(BB),C2(BB),COULS(BB), 2 COULP(BB),COULD(BB),X(MAXATM),Y(MAXATM),Z(MAXATM) INTEGER*2 SYMBOL,KEY,VELEC FORT1282 COMMON/ITPARM/DAMP1,DAMP2,DAMP3,LAMPRI,DELTAC,SENSE,MAXCYC,PRTCYC, 1 ICYCLE,NCON,PARTIT,PRINTX,ITABLE(20) FORT1284 REAL*8 LAMPRI FORT1285 INTEGER*4 PRTCYC FORT1286 LOGICAL*1 PARTIT,PRINTX,ITABLE FORT1287 COMMON/STARS/STAR,STAR2 FORT1288 INTEGER*2 STAR,STAR2 FORT1289 INTEGER*4 ONE,TWO,STAR1,TRUE FORT1290 DATA ONE,TWO,STAR1/'1','2','*'/ FORT1291 IF(ONEMAT.AND.IPRINT.LT.-2) RETURN FORT1292 ICYCLE=1 FORT1293 C FORT1294 C SETUP DEFAULT ORBITAL OCCUPATIONS. FORT1295 C FORT1296 IF(L1) GO TO 3 FORT1297 IH=NELEC/2 FORT1298 JH=NDIM+1-IH FORT1299 DO 1 I=1,JH FORT1300 1 IOCC(I)=0.0 FORT1301 DO 2 I=JH,NDIM FORT1302 2 IOCC(I)=2.0 FORT1303 IF(IH+IH.NE.NELEC) IOCC(JH-1)=1.0 FORT1304 GO TO 500 FORT1305 C FORT1306 C PROVISION FOR READING IN USER SPECIFIED OCCUPATIONS. FORT1307 C ALSO PROVISION FOR NON-INTEGER OCCUPATIONS. FORT1308 C FORT1309 3 READ(5,2000) (MAXD(I),I=1,NDIM) FORT1310 2000 FORMAT(80A1) FORT1311 TEL=0.0D0 FORT1312 DO 4 I=1,NDIM FORT1313 J=NDIM+1-I FORT1314 IOCC(J)=0.0 FORT1315 IF(MAXD(I).EQ.ONE) IOCC(J)=1.0 FORT1316 IF(MAXD(I).EQ.TWO) IOCC(J)=2.0 FORT1317 IF(MAXD(I).EQ.STAR1) READ(5,900) IOCC(J) FORT1318 900 FORMAT(F15.8) FORT1319 4 TEL=TEL+DBLE(IOCC(J)) FORT1320 TEL=DABS(TEL-DFLOAT(NELEC)) FORT1321 IF(TEL.LT.0.001D0) GO TO 500 FORT1322 WRITE(6,2001) FORT1323 2001 FORMAT('0*** WARNING **** ORBITAL POPULATIONS INCONSISTENT WITHFORT1324 1 ASSUMED CHARGE ON MOLECULE',////,T10,'I',T25,'IOCC(I)',/) FORT1325 DO 99 I=1,NDIM FORT1326 99 WRITE(6,2002) I,IOCC(I) FORT1327 2002 FORMAT(T8,I3,T22,F12.8) FORT1328 STOP FORT1329 C FORT1330 C CALL GRMSCH TO ORTHOGONALISE BASIS SET. FORT1331 C FORT1332 500 CALL GRMSCH(S,C,NDIM) FORT1333 CON=.5D0*CON FORT1334 IF(.NOT.ITERAT) GO TO 15 FORT1335 DO 5 I=1,NDIM FORT1336 5 SORB(I)=0.D0 FORT1337 DO 10 I=1,NATOM FORT1338 X(I)=0.D0 FORT1339 Z(I)=0.D0 FORT1340 10 COUL0(I)=0.D0 FORT1341 15 IF(.NOT.ONEMAT) GO TO 25 FORT1342 C FORT1343 C IN ONE MATRIX CASE, STUFF DIAGONAL ELEMENTS OF S INTO SP. FORT1344 C FORT1345 DO 20 I=1,NDIM FORT1346 20 SP(I)=S(I,I) FORT1347 C FORT1348 C SETUP DIAGONAL ELEMENTS OF HUCKEL MATRIX IN H(I,J). FORT1349 C FORT1350 25 IF(NH.EQ.0) GO TO 35 FORT1351 ET=COULH FORT1352 DO 30 I=1,NH FORT1353 IF(ITERAT) ET=COULH-COUL0(I) FORT1354 30 C(I)=ET FORT1355 35 IH=NH+1 FORT1356 ET=0.D0 FORT1357 DO 50 I=1,NA FORT1358 KEYI=KEY(I) FORT1359 IF(ITERAT) ET=COUL0(I+NH) FORT1360 C(IH)=COULS(KEYI)-ET FORT1361 IH=IH+1 FORT1362 IF(NP(KEYI).EQ.0) GO TO 50 FORT1363 HH=COULP(KEYI)-ET FORT1364 JH=IH+2 FORT1365 ASSIGN 40 TO IL FORT1366 GO TO 42 FORT1367 40 IF(ND(KEYI).EQ.0) GO TO 50 FORT1368 HH=COULD(KEYI)-ET FORT1369 JH=IH+4 FORT1370 ASSIGN 50 TO IL FORT1371 42 DO 45 J=IH,JH FORT1372 45 C(J)=HH FORT1373 IH=JH+1 FORT1374 GO TO IL,(40,50) FORT1375 50 CONTINUE FORT1376 DO 55 I=1,NDIM FORT1377 55 H(I,I)=C(I) FORT1378 IF(NHDG.EQ.1) GO TO 59 FORT1379 DO 56 I=1,NDIM FORT1380 56 HDG(I)=C(I) FORT1381 C FORT1382 C SETUP OFF-DIAGONAL ELEMENTS OF HUCKEL MATRIX. FORT1383 C FORT1384 59 CNST=CON FORT1385 DO 58 I=2,NDIM FORT1386 IL=I-1 FORT1387 DO 58 J=1,IL FORT1388 HH=C(I)+C(J) FORT1389 IF(.NOT.L5) GO TO 58 FORT1390 ET=(C(I)-C(J))/HH FORT1391 ET=ET*ET FORT1392 CNST=CON+ET/2.0D0+ET*ET*(0.5D0-CON) FORT1393 58 H(I,J)=CNST*HH*S(I,J) FORT1394 IF(ONEMAT) GO TO 100 FORT1395 DO 60 I=2,NDIM FORT1396 IL=I-1 FORT1397 DO 60 J=1,IL FORT1398 60 H(J,I)=H(I,J) FORT1399 C FORT1400 C PRINT OUT HUCKEL MATRIX. PRINT OUT TITLE IF METH IS NOT FORT1401 C EQUAL TO ZERO. FORT1402 C FORT1403 806 IF(ICYCLE.GT.1) GO TO 800 FORT1404 IF(PRT(6)) GO TO 805 FORT1405 IF(METH.EQ.0) GO TO 801 FORT1406 GO TO 802 FORT1407 800 IF(ICYCLE.GE.10000) WRITE(6,701) AB FORT1408 701 FORMAT('RESULTS OF CALCULATION ',8A8,A6,A2,//) FORT1409 IF(PRT(7).OR..NOT.PRINTX) GO TO 805 FORT1410 801 WRITE(6,803) FORT1411 803 FORMAT('HUCKEL MATRIX') FORT1412 CALL PEGLEG(H,NDIM,NDIM) FORT1413 GO TO 805 FORT1414 802 WRITE(6,804) FORT1415 804 FORMAT('INPUT HUCKEL MATRIX') FORT1416 CALL PEGLEG(H,NDIM,NDIM) FORT1417 805 IF(ICYCLE.EQ.1.AND.PUN(6)) WRITE(7,825) H FORT1418 IF(ICYCLE.GT.1.AND.PUN(7).AND.PRINTX) WRITE(7,825) H FORT1419 825 FORMAT(8F9.5) FORT1420 C FORT1421 C IF CALCULATING ENERGY MATRIX, STORE H(I,I) IN X(I),Y(I),Z(I). FORT1422 C FORT1423 IF(ICYCLE.LE.MAXCYC.AND.METH.NE.0) GO TO 100 FORT1424 IF(NH.EQ.0) GO TO 369 FORT1425 DO 370 I=1,NH FORT1426 370 X(I)=H(I,I) FORT1427 369 IH=NH+1 FORT1428 JH=NH+1 FORT1429 DO 371 I=1,NA FORT1430 KEYI=KEY(I) FORT1431 X(IH)=H(JH,JH) FORT1432 JH=JH+1 FORT1433 IF(NP(KEYI).EQ.0) GO TO 371 FORT1434 Y(IH)=H(JH,JH) FORT1435 JH=JH+3 FORT1436 IF(ND(KEYI).EQ.0) GO TO 371 FORT1437 Z(IH)=H(JH,JH) FORT1438 JH=JH+5 FORT1439 371 IH=IH+1 FORT1440 C FORT1441 C CALL TRNFRM TO TRANSFORM HUCKEL MATRIX TO ORTHOGONAL BASIS SET. FORT1442 C THEN CALL GIVENS TO PERFORM DIAGONALIZATION. FORT1443 C FORT1444 100 IH=1 FORT1445 IF(ONEMAT) IH=2 FORT1446 CALL TRNFRM(S,H,C,COUL0,NDIM,SP,IH) FORT1447 IF(ONEMAT) GO TO 110 FORT1448 IH=NDIM FORT1449 GO TO 120 FORT1450 110 IH=-NDIM FORT1451 120 CALL GIVENS(NDIM,IH,NDIM,C,SP,COUL0,H) FORT1452 130 IF(ICYCLE.GE.10000) ITERAT=.FALSE. FORT1453 C FORT1454 C PRINT OUT TITLE, ENERGY LEVELS, AND OCCUPATION NUMBERS. FORT1455 C FORT1456 IF(ITERAT) GO TO 700 FORT1457 IF(METH.EQ.0) WRITE(6,701) AB FORT1458 IF(PRT(8)) GO TO 710 FORT1459 WRITE(6,702) (I,COUL0(I),IOCC(I),I=1,NDIM) FORT1460 702 FORMAT(////,57X,'ENERGY LEVELS (EV)',/,(/52X,'E(',I3,') =',F12.5, FORT1461 18X,F6.4)) FORT1462 710 IF(PUN(8)) WRITE(7,825) (COUL0(I),I=1,NDIM) FORT1463 700 IF(ONEMAT) GO TO 200 FORT1464 C FORT1465 C DIDDLE WITH C,H. FORT1466 C FORT1467 DO 160 J=1,NDIM FORT1468 DO 140 K=1,NDIM FORT1469 140 C(K)=H(K,IH) FORT1470 DO 155 I=1,NDIM FORT1471 ET=0.D0 FORT1472 DO 150 K=I,NDIM FORT1473 150 ET=ET+S(I,K)*C(K) FORT1474 155 H(I,IH)=ET FORT1475 160 IH=IH-1 FORT1476 K=1 FORT1477 DO 180 I=2,NDIM FORT1478 IL=I-1 FORT1479 DO 180 J=1,IL FORT1480 C(K)=S(I,J) FORT1481 180 K=K+1 FORT1482 200 IF(METH.GT.1.AND.ITERAT) GO TO 210 FORT1483 C FORT1484 C CALL OUTPUT FOR FINAL PRINT OUT OF RESULTS. FORT1485 C FORT1486 205 CALL OUTPUT(H,S,MAD,C,COUL0,SORB,IOCC,HDG,NDIM,NTYPE,NC,NHDG) FORT1487 IF(.NOT.ITERAT) GO TO 999 FORT1488 GO TO 220 FORT1489 C FORT1490 C IF DOING CHARGE ITERATIVE CALCULATION ( METH >1 ), CALL ITRATE FORT1491 C TO SETUP HUCKEL MATRIX. FORT1492 C FORT1493 210 CALL ITRATE(H,S,MAD,C,COUL0,SORB,IOCC,HDG,NDIM,NTYPE,NC,NHDG) FORT1494 IF(.NOT.ITERAT) GO TO 205 FORT1495 220 IF(ICYCLE.GT.MAXCYC) ICYCLE=10000 FORT1496 IF(METH.GT.1) GO TO 806 FORT1497 GO TO 15 FORT1498 999 RETURN FORT1499 END FORT1500 SUBROUTINE GIVENS (NX,NROOTX,NJX,A,B,ROOT,VECT) FORT1501 C FORT1502 C SUBROUTINE TO CALCULATE THE EIGENVALUES AND EIGENVECTORS FORT1503 C OF A REAL SYMMETRIC MATRIX. FORT1504 C FORT1505 C FORT1506 C THE PARAMETERS FOR THE ROUTINE ARE0 FORT1507 C FORT1508 C NX ORDER OF MATRIX. FORT1509 C FORT1510 C NROOTX NUMBER OF ROOTS FOR WHICH EIGENVECTORS ARE WANTED. FORT1511 C IF NO VECTORS ARE WANTED, MAKE NROOTX NEGATIVE. FORT1512 C FORT1513 C NJX ROW DIMENSION OF VECT ARRAY. SEE 'VECT' BELOW. FORT1514 C NJX MUST BE NOT LESS THAN NX. FORT1515 C FORT1516 C A MATRIX STORED BY COLUMNS IN PACKED UPPER TRIANGULAR FORT1517 C FORM, I.E. OCCUPYING NX*(NX+1)/2 CONSECUTIVE FORT1518 C LOCATIONS. FORT1519 C FORT1520 C B SCRATCH ARRAY USED BY GIVENS. MUST BE AT LEAST FORT1521 C NX*6 CELLS. FORT1522 C FORT1523 C ROOT ARRAY TO HOLD THE EIGENVALUES. MUST BE AT LEAST FORT1524 C NX CELLS LONG. THE ROOTS ARE ORDERED LARGEST FIRST FORT1525 C IN THIS ARRAY. FORT1526 C FORT1527 C VECT EIGENVECTOR ARRAY. EACH COLUMN WILL HOLD AN FORT1528 C EIGENVECTOR FOR THE CORRESPONDING ROOT. MUST BE FORT1529 C DIMENSIONED WITH 'NJX' ROWS AND AT LEAST 'NJX' FORT1530 C COLUMNS, UNLESS NO VECTORS ARE REQUESTED (NEGATIVE FORT1531 C NROOTX). IN THIS LATTER CASE, THE ARGUMENT VECT FORT1532 C IS JUST A DUMMY, AND THE STORAGE IS NOT USED. FORT1533 C THE EIGENVECTORS ARE NORMALIZED TO UNIT LENGTH. FORT1534 C FORT1535 C THE ARRAYS A AND B ARE DESTROYED BY THE COMPUTATION. THE FORT1536 C RESULTS APPEAR IN ROOT AND VECT. FORT1537 C FORT1538 C FOR PROPER FUNCTIONING OF THIS ROUTINE, THE RESULT OF A FLOATING FORT1539 C POINT UNDERFLOW SHOULD BE A ZERO. FORT1540 C FORT1541 C THE ORIGINAL REFERENCE TO THE GIVENS TECHNIQUE IS IN OAK RIDGE FORT1542 C REPORT NUMBER ORNL 1574 (PHYSICS), BY WALLACE GIVENS. FORT1543 C FORT1544 C THE METHOD AS PRESENTED IN THIS PROGRAM CONSISTS OF FOUR STEPS0 FORT1545 C FORT1546 C FIRST, THE INPUT MATRIX IS REDUCED TO TRIDIAGONAL FORM BY THE FORT1547 C HOUSEHOLDER TECHNIQUE (J. H. WILKINSON, COMP. J. 3, 23 (1960)). FORT1548 C THE EIGENVALUES OF THE TRIDIAGONAL MATRIX ARE THEN FOUND USING FORT1549 C THE QR TRANSFORM METHOD. SEE J. H. WILKINSON, THE ALGEBRAIC FORT1550 C EIGENVALUE PROBLEM(1965) FOR A DESCRIPTION OF THIS ALGORITHM. FORT1551 C THE EIGENVECTORS OF THE TRIDIAGONAL FORM ARE THEN EVALUATED FORT1552 C (J. H. WILKINSON, COMP. J. 1, 90 (1958)), BY THE METHOD OF FORT1553 C INVERSE ITERATION, FOR NONDEGENERATE MATRICES. FORT1554 C FOR MATRICES WITH DEGENERATE OR NEAR-DEGENERATE EIGENVALUES, FORT1555 C THE EIGENVECTORS ARE EVALUATED INSTEAD BY FURTHER QR TRANSFORMS. FORT1556 C THIS METHOD GIVES ORTHOGONAL VECTORS EVEN FOR DEGENERATE ROOTS. FORT1557 C FINALLY THE TRIDIAGONAL VECTORS ARE ROTATED TO VECTORS OF THE FORT1558 C ORIGINAL ARRAY (FIRST REFERENCE). FORT1559 C FORT1560 C THE INVERSE ITERATION PORTION OF THIS PROGRAM WAS ADAPTED FORT1561 C FROM THE QUANTUM CHEMISTRY PROGRAM EXCHANGE NUMBER 62.1, BY FORT1562 C FRANKLIN PROSSER. THE E