| http://server.ccl.net/cca/software/SOURCES/FORTRAN/EHT/forticon8.f.shtml |
|
CCL forticon8.f | |||
|
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 EIGENVALUE SUBROUTINE (EVQR) WAS WRITTEN FORT1563
C BY WALTER NIELSEN. FORT1564
C ROY GORDON, SEPT. 1969 FORT1565
C FORT1566
C AN EXCELLENT PRESENTATION OF THE GIVENS TECHNIQUE IS FOUND IN FORT1567
C J. M. ORTEGA'S ARTICLE IN 'MATHEMATICS FOR DIGITAL COMPUTERS,' FORT1568
C VOLUMD 2, ED. BY RALSTON AND WILF, WILEY (1967), PAGE 94. FORT1569
C FORT1570
IMPLICIT DOUBLE PRECISION(A-H,R-Z) FORT1571
COMMON /VECTOR/ FACT,IDIF FORT1572
DIMENSION B(NX,6),A(1),ROOT(NX),VECT(NJX,NROOTX) FORT1573
C FORT1574
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * FORT1575
C * * FORT1576
C * USERS PLEASE NOTE0 TWO PARAMETERS, ETA AND THETA, SHOULD BE * FORT1577
C * ADJUSTED BY THE USER FOR HIS PARTICULAR MACHINE. * FORT1578
C * * FORT1579
C * ETA IS AN INDICATION OF THE PRECISION OF THE FLOATING POINT * FORT1580
C * REPRESENTATION ON THE COMPUTER BEING USED (ROUGHLY 10**(-M), * FORT1581
C * WHERE M IS THE NUMBER OF DECIMALS OF PRECISION ). * FORT1582
C * * FORT1583
C * THETA IS AN INDICATION OF THE RANGE OF NUMBERS THAT CAN BE * FORT1584
C * EXPRESSED IN THE FLOATING POINT REPRESENTATION (ROUGHLY THE * FORT1585
C * LARGEST NUMBER). * FORT1586
C * * FORT1587
C * SOME RECOMMENDED VALUES FOLLOW. * FORT1588
C * * FORT1589
C * FOR IBM 7094, UNIVAC 1108, ETC. (27-BIT BINARY FRACTION, * FORT1590
C * 8-BIT BINARY EXPONENT), ETA=1.E-8, THETA=1.E37. * FORT1591
C * FOR CONTROL DATA 3600 (36-BIT BINARY FRACTION, 11-BIT BINARY * FORT1592
C * EXPONENT), ETA=1.E-11, THETA=1.E307. * FORT1593
C * FOR CONTROL DATA 6600 (48-BIT BINARY FRACTION, 11-BIT BINARY * FORT1594
C * EXPONENT), ETA=1.E-14, THETA=1.E307. * FORT1595
C * FOR IBM 360/50 AND 360/65 DOUBLE PRECISION (56-BIT HEXA- * FORT1596
C * DECIMAL FRACTION, 7-BIT HEXADECIMAL EXPONENT), ETA=1.E-16, * FORT1597
C * THETA=1.E75. * FORT1598
C * * FORT1599
C * OTHER PARAMETERS WHICH MUST BE ADJUSTED ARE0 * FORT1600
C * * FORT1601
C * DEL1 = ETA/1.D2, DELTA = ETA**2*1.D2, SMALL = ETA**2/1.D2, * FORT1602
C * DELBIG = THETA*DELTA/1.D3, THETA1 = 1.D3/THETA, EMAG = ETA, * FORT1603
C * TOLER = 1.D2*DSQRT(ETA) * FORT1604
C * * FORT1605
C * TOLER IS A FACTOR USED TO DETERMINE IF ANY ROOTS ARE CLOSE * FORT1606
C * ENOUGH TOGETHER TO BE CONSIDERED DEGENERATE FOR PURPOSES OF * FORT1607
C * CALCULATING EIGENVECTORS. FOR THE MATRIX NORMED TO UNITY, IF * FORT1608
C * THE DIFFERENCE BETWEEN TWO ROOTS IS LESS THAN TOLER, THEN THE * FORT1609
C * QR TRANSFORMATION IS USED TO FORM THE EIGENVECTORS. * FORT1610
C * * FORT1611
C * EMAG IS A TOLERANCE FOR NEGLIGIBLE ELEMENTS IN THE QR * FORT1612
C * ITERATION FOR EIGENVECTORS FOR DEGENERATE EIGENVALUES. * FORT1613
C * * FORT1614
C * IN THE FOLLOWING ROUTINE, ETA = 1.D-16 AND THETA = 1.D75. * FORT1615
C * * FORT1616
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * FORT1617
C FORT1618
DATA ETA/1.D-16/,THETA/1.D38/,DEL1/1.D-18/,DELTA/1.D-30/, FORT1619
*SMALL/1.D-34/,DELBIG/1.D05/,THETA1/1.D-38/,TOLER/1.D-6/, FORT1620
*EMAG/1.D-16/ FORT1621
N = NX FORT1622
FLOATN = DFLOAT(N) FORT1623
NROOT = IABS(NROOTX) FORT1624
IF (NROOT.EQ.0) GO TO 1001 FORT1625
IF (N-1) 1001,1003,105 FORT1626
1003 ROOT(1) = A(1) FORT1627
IF (NROOTX.GT.0) VECT(1,1) = 1.D0 FORT1628
GO TO 1001 FORT1629
105 CONTINUE FORT1630
C FORT1631
C NSIZE IS THE NUMBER OF ELEMENTS IN THE PACKED ARRAY. FORT1632
C FORT1633
NSIZE = (N*(N+1))/2 FORT1634
NM1 = N-1 FORT1635
NM2 = N-2 FORT1636
NP1 = N+1 FORT1637
C FORT1638
C COMPUTE TRACE. FORT1639
C FORT1640
TRACE = 0.D0 FORT1641
JUMP = 1 FORT1642
DO 1 J=2,NP1 FORT1643
TRACE = TRACE + A(JUMP) FORT1644
1 JUMP = JUMP + J FORT1645
TRACE = TRACE/FLOATN FORT1646
C FORT1647
C SUBTRACT TRACE FROM DIAGONAL ELEMENTS TO GIVE A MORE RELIABLE NORM FORT1648
C WHEN THERE ARE LARGE DIAGONAL ELEMENTS. FORT1649
C FORT1650
JUMP = 1 FORT1651
DO 2 J=2,NP1 FORT1652
A(JUMP) = A(JUMP) - TRACE FORT1653
2 JUMP = JUMP + J FORT1654
C FORT1655
C SCALE MATRIX TO EUCLIDEAN NORM OF 1. SCALE FACTOR IS ANORM. FORT1656
C FORT1657
FACTOR = 0.D0 FORT1658
DO 70 I=1,NSIZE FORT1659
70 FACTOR = DMAX1(FACTOR,DABS(A(I))) FORT1660
IF (FACTOR.NE.0.D0) GO TO 72 FORT1661
C FORT1662
C NULL MATRIX. FIX UP ROOTS AND VECTORS, THEN EXIT. FORT1663
C FORT1664
DO 78 I=1,NROOT FORT1665
IF (NROOTX.LT.0) GO TO 78 FORT1666
DO 77 J=1,N FORT1667
77 VECT(J,I) = 0.D0 FORT1668
VECT(I,I) = 1.D0 FORT1669
78 ROOT(I) = 0.D0 FORT1670
GO TO 1001 FORT1671
72 ANORM = 0.D0 FORT1672
86 SCALE = 1.D0/FACTOR FORT1673
DO 80 I=1,NSIZE FORT1674
80 ANORM = ANORM + (A(I)*SCALE)**2 FORT1675
ANORM = ANORM+ANORM FORT1676
C FORT1677
C SUBTRACT DIAGONAL CONTRIBUTIONS WHICH WERE COUNTED TWICE. FORT1678
C FORT1679
JUMP = 1 FORT1680
DO 81 J=2,NP1 FORT1681
ANORM = ANORM -(A(JUMP)*SCALE)**2 FORT1682
81 JUMP = JUMP + J FORT1683
83 ANORM = FACTOR*DSQRT(ANORM) FORT1684
SCALE = 1.D0/ANORM FORT1685
DO 91 I=1,NSIZE FORT1686
91 A(I) = A(I)*SCALE FORT1687
ALIMIT = 1.D0 FORT1688
C FORT1689
C TRIDIAGONALIZATION OF SYMMETRIC MATRIX. FORT1690
C FORT1691
ID = 0 FORT1692
IA = 1 FORT1693
IF (NM2.EQ.0) GO TO 201 FORT1694
DO 200 J=1,NM2 FORT1695
C FORT1696
C J COUNTS ROW OF A MATRIX TO BE DIAGONALIZED. IA INDICATES START OF FORT1697
C NON-CODIAGONAL ELEMENTS IN THE ROW. ID IS THE INDEX OF CODIAGONAL FORT1698
C ELEMENT ON THE ROW BEING CODIAGONALIZED. FORT1699
C FORT1700
IA = IA+J+2 FORT1701
ID = ID+J+1 FORT1702
JP2 = J + 2 FORT1703
J1 = J + 1 FORT1704
C FORT1705
C FIND LIMITS FOR BAND OF SIGNIFICANT MATRIX ELEMENTS. FORT1706
C FORT1707
LIMIT = J1 FORT1708
II = IA FORT1709
DO 99 I=JP2,N FORT1710
B(I,5) = A(II) FORT1711
IF (DABS(B(I,5)).GT.DEL1) LIMIT = I FORT1712
99 II = II + I FORT1713
DTEMP = A(ID) FORT1714
IF (LIMIT.GT.J1) GO TO 110 FORT1715
C FORT1716
C NO TRANSFORMATION NECESSARY IF ALL THE NON-CODIAGONAL FORT1717
C ELEMENTS ARE TINY. FORT1718
C FORT1719
120 B(J,1) = DTEMP FORT1720
A(ID) = 0.D0 FORT1721
GO TO 200 FORT1722
C FORT1723
C SUM SQUARES OF SIGNIFICANT NON-CODIAGONAL ELEMENTS OF ROW J. FORT1724
C FORT1725
110 IDIF = LIMIT -JP2 FORT1726
SUM = DOT(B(JP2,5),B(JP2,5)) FORT1727
C FORT1728
C NOW COMPLETE THE SUM OF OFF-DIAGONAL SQUARES. FORT1729
C FORT1730
SUM = DSQRT(SUM + DTEMP**2) FORT1731
C FORT1732
C NEW CODIAGONAL ELEMENT. FORT1733
C FORT1734
B(J,1) = -DSIGN(SUM,DTEMP) FORT1735
C FORT1736
C FIRST NON-ZERO ELEMENT OF THIS W-VECTOR. FORT1737
C FORT1738
B(J+1,2) = DSQRT((1.D0 + DABS(DTEMP)/SUM)*5.D-1) FORT1739
C FORT1740
C FORM REST OF THE W-VECTOR ELEMENTS. FORT1741
C FORT1742
TEMP = DSIGN(5.D-1/(B(J+1,2)*SUM),DTEMP) FORT1743
II = IA FORT1744
DO 130 I=JP2,LIMIT FORT1745
B(I,2) = A(II)*TEMP FORT1746
130 II = II + I FORT1747
C FORT1748
C FORM P-VECTOR AND SCALAR. P-VECTOR = A-MATRIX*W-VECTOR. FORT1749
C SCALAR = W-VECTOR*P-VECTOR. FORT1750
C FORT1751
DAK = 0.D0 FORT1752
C FORT1753
C IC IS THE LOCATION OF THE NEXT DIAGONAL ELEMENT. I RUNS OVER THE FORT1754
C NON-ZERO P-ELEMENTS. CASES FOR I LESS THAN LIMIT. FORT1755
C FORT1756
IC = ID + 1 FORT1757
LIMLES = LIMIT - 1 FORT1758
DO 188 I=J1,LIMLES FORT1759
C FORT1760
C FORM FIRST PART OF P ELEMENT THEN MOVE IC TO TOP OF NEXT FORT1761
C A-MATRIX 'ROW'. FORT1762
C FORT1763
IDIF = I - J1 FORT1764
DTEMP = DOT(B(J1,2),A(IC)) FORT1765
IC = IC + I FORT1766
C FORT1767
C COMPLETE P ELEMENT. CHANGE INCREMENTING MODE AT DIAGONAL ELEMENT. FORT1768
C FORT1769
178 IP1 = I + 1 FORT1770
JJ = IC + IDIF FORT1771
DTEMP = DTEMP + DSUM(B(N,1),A(JJ),IP1,LIMIT) FORT1772
C FORT1773
C BUILD UP THE K-SCALAR (AK). FORT1774
C FORT1775
DAK = DAK + DTEMP*B(I,2) FORT1776
188 B(I,1) = DTEMP FORT1777
C FORT1778
C CASE FOR I = LIMIT. FORT1779
C FORT1780
IDIF = LIMIT - J1 FORT1781
DTEMP = DOT(B(J1,2),A(IC)) FORT1782
DAK = DAK + DTEMP*B(LIMIT,2) FORT1783
B(LIMIT,1) = DTEMP FORT1784
IDIF = LIMIT - J1 FORT1785
C FORT1786
C TEST TO SEE IF ANY I VALUES REMAIN. DO REMAINING VALUES. FORT1787
C FORT1788
IF (LIMIT.EQ.N) GO TO 190 FORT1789
IC = IC + LIMIT FORT1790
LIMLO = LIMIT + 1 FORT1791
DO 189 I=LIMLO,N FORT1792
B(I,1) = DOT(B(J1,2),A(IC)) FORT1793
B(I,2) = 0.D0 FORT1794
189 IC = IC + I FORT1795
C FORT1796
C FORM THE Q-VECTOR. FORT1797
C FORT1798
190 FACT = -DAK FORT1799
CALL VECSUM(B(J1,1),B(J1,2)) FORT1800
C FORT1801
C TRANSFORM THE REST OF THE A-MATRIX. JJ INDICATES START-1 OF THE FORT1802
C REST OF THE A-MATRIX. MOVE W-VECTOR INTO THE OLD A-MATRIX LOCATIONS FORT1803
C TO SAVE SPACE. I RUNS OVER THE SIGNIFICANT ELEMENTS OF THE W-VECTOR. FORT1804
C FORT1805
JJ = ID FORT1806
DO 160 I=J1,N FORT1807
A(JJ) = B(I,2) FORT1808
IF (I.GT.LIMIT) GO TO 161 FORT1809
B2 = B(I,2) FORT1810
FACT = -B2 - B2 FORT1811
IDIF = I - J1 FORT1812
CALL VECSUM(A(JJ+1),B(J1,1)) FORT1813
161 B1 = B(I,1) FORT1814
FACT = -B1 - B1 FORT1815
IDIF = MIN0(I,LIMIT) - J1 FORT1816
CALL VECSUM(A(JJ+1),B(J1,2)) FORT1817
160 JJ = JJ + I FORT1818
C FORT1819
C STORE AWAY LIMIT FOR LATER USE IN BACK TRANSFORMATION. MOVE LAST FORT1820
C CODIAGONAL ELEMENT OUT INTO ITS PROPER PLACE. FORT1821
C FORT1822
200 B(J,6) = LIMIT FORT1823
201 CONTINUE FORT1824
B(NM1,1) = A(NSIZE-1) FORT1825
A(NSIZE-1) = 0.D0 FORT1826
C FORT1827
C USE QR TRANSFORM METHOD TO FIND EIGENVALUES OF THE TRIDIAGONAL FORT1828
C MATRIX. MOVE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX INTO FORT1829
C ROOT ARRAY. THIS IS A MORE CONVENIENT INDEXING POSITION. ALSO, FORT1830
C PUT SQUARE OF CODIAGONAL ELEMENTS IN THIRD N ELEMENTS. FORT1831
C FORT1832
JUMP = 0 FORT1833
DO 320 J=1,NM1 FORT1834
JUMP = JUMP + J FORT1835
ROOT(J) = A(JUMP) FORT1836
320 B(J,3) = B(J,1)**2 FORT1837
ROOT(N) = A(NSIZE) FORT1838
CALL EVQR(ROOT,B(1,3),N,30,SMALL) FORT1839
C FORT1840
C ROOT NOW CONTAINS THE SHIFTED AND SCALED EIGENVALUES. STORE FORT1841
C EIGENVALUES FOR POSSIBLE LATER USE AS SHIFTS IN EVALUATING FORT1842
C EIGENVECTORS FOR DEGENERATE MATRICES. FORT1843
C FORT1844
DO 325 J=1,N FORT1845
325 B(J,2) = ROOT(J) FORT1846
C FORT1847
C SORT THE EIGENVALUES INTO DESCENDING ALGEBRAIC ORDER. FORT1848
C FORT1849
DO 330 I=1,NM1 FORT1850
IP1 = I + 1 FORT1851
DO 330 J=IP1,N FORT1852
IF (ROOT(I).GE.ROOT(J)) GO TO 330 FORT1853
TEMP = ROOT(I) FORT1854
ROOT(I) = ROOT(J) FORT1855
ROOT(J) = TEMP FORT1856
330 CONTINUE FORT1857
C FORT1858
C QUIT NOW IF NO VECTORS WERE REQUESTED. OTHERWISE, TEST FOR FORT1859
C DEGENERACY OR NEAR DEGENERACY OF EIGENVALUES FOR WHICH EIGENVECTORS FORT1860
C WERE REQUESTED. IF ONLY ONE VECTOR REQUESTED, DEGENERACY DOESN'T FORT1861
C MATTER. FORT1862
C FORT1863
IF (NROOTX.LT.0) GO TO 1002 FORT1864
IF (NROOTX.EQ.1) GO TO 807 FORT1865
NTOP = NROOT - 1 FORT1866
DO 400 I=1,NTOP FORT1867
IF (DABS(ROOT(I+1)-ROOT(I)).LE.TOLER) GO TO 410 FORT1868
400 CONTINUE FORT1869
C FORT1870
C NEXT STATEMENT IS REACHED IF ALL EIGENVALUES FOR WHICH EIGENVECTORS FORT1871
C WERE REQUESTED ARE WELL SEPARATED. FORT1872
C FORT1873
GO TO 807 FORT1874
C FORT1875
C THE FOLLOWING IS REACHED IF THERE ARE ANY DEGENERATE CLUSTERS FORT1876
C OF EIGENVALUES. USE FURTHER QR TRANSFORMS TO EVALUATE FORT1877
C THE EIGENVECTORS OF THE TRIDIAGONAL MATRIX. THIS METHOD FORT1878
C GIVES ORTHOGONAL EIGENVECTORS EVEN WHEN THE EIGENVALUES ARE FORT1879
C DEGENERATE. HOWEVER, IT TAKES MORE ARITHMETIC THAN THE METHOD FORT1880
C OF INVERSE ITERATION (AT LEAST AT LARGE N). FORT1881
C FORT1882
C PUT DIAGONAL ELEMENTS OF TRIDIAGONAL MATRIX INTO ROOT. PUT OFF- FORT1883
C DIAGONAL ELEMENTS INTO B(I,3). FORT1884
C FORT1885
410 JUMP = 0 FORT1886
DO 440 J=1,NM1 FORT1887
JUMP = JUMP + J FORT1888
ROOT(J) = A(JUMP) FORT1889
440 B(J,3) = B(J,1) FORT1890
C FORT1891
C LAST DIAGONAL ELEMENT. FORT1892
C FORT1893
ROOT(N) = A(NSIZE) FORT1894
C FORT1895
C INITIALIZE VECTORS TO A UNIT MATRIX. FORT1896
C FORT1897
DO 450 I=1,N FORT1898
DO 445 J=1,N FORT1899
445 VECT(J,I) = 0.D0 FORT1900
450 VECT(I,I) = 1.D0 FORT1901
C FORT1902
C FORM EIGENVECTORS OF TRIDIAGONAL MATRIX FOR DEGENERATE FORT1903
C MATRICES AND TRANSPOSE THE VECTORS. FORT1904
C FORT1905
CALL QRTN(ROOT,B(1,3),VECT,B(1,2),N,25,EMAG,NJX) FORT1906
DO 456 I=1,NM1 FORT1907
IP1 = I + 1 FORT1908
DO 455 J=IP1,N FORT1909
FLIP = VECT(I,J) FORT1910
VECT(I,J) = VECT(J,I) FORT1911
455 VECT(J,I) = FLIP FORT1912
456 CONTINUE FORT1913
C FORT1914
C IF ROOTS WERE NOT LOCATED IN DESCENDING ORDER, INTERCHANGE ROOTS FORT1915
C AND VECTORS. FORT1916
C FORT1917
ITOP = NROOT-1 FORT1918
DO 480 I=1,ITOP FORT1919
IP1 = I+1 FORT1920
DO 480 J=IP1,N FORT1921
IF (ROOT(I).GE.ROOT(J)) GO TO 480 FORT1922
TEMP = ROOT(I) FORT1923
ROOT(I) = ROOT(J) FORT1924
ROOT(J) = TEMP FORT1925
DO 470 K=1,N FORT1926
TEMP = VECT(K,I) FORT1927
VECT(K,I) = VECT(K,J) FORT1928
470 VECT(K,J) = TEMP FORT1929
480 CONTINUE FORT1930
C FORT1931
C DEGENERATE VECTORS ARE NOW COMPLETE. FORT1932
C FORT1933
GO TO 940 FORT1934
C FORT1935
C EIGENVECTORS OF TRIDIAGONAL MATRIX FOR NONDEGENERATE MATRICES. FORT1936
C INITIALIZE VECTOR ARRAY. FORT1937
C FORT1938
807 CONTINUE FORT1939
DO 705 I=1,NROOT FORT1940
DO 15 J=1,N FORT1941
15 VECT(J,I) = 1.D0 FORT1942
705 CONTINUE FORT1943
DO 700 I=1,NROOT FORT1944
C FORT1945
C USE INVERSE ITERATION TO FIND VECTORS. FORT1946
C FORT1947
701 AROOT = ROOT(I) FORT1948
ELIM1 = A(1) - AROOT FORT1949
ELIM2 = B(1,1) FORT1950
JUMP = 1 FORT1951
DO 750 J=1,NM1 FORT1952
JUMP = JUMP + J + 1 FORT1953
C FORT1954
C GET THE CORRECT PIVOT EQUATION FOR THIS STEP. FORT1955
C FORT1956
IF (DABS(ELIM1).LE.DABS(B(J,1))) GO TO 760 FORT1957
C FORT1958
C FIRST (ELIM1) EQUATION IS THE PIVOT THIS TIME. CASE 1. FORT1959
C FORT1960
B(J,2) = ELIM1 FORT1961
B(J,3) = ELIM2 FORT1962
B(J,4) = 0.D0 FORT1963
TEMP = B(J,1)/ELIM1 FORT1964
ELIM1 = A(JUMP) - AROOT - TEMP*ELIM2 FORT1965
ELIM2 = B(J+1,1) FORT1966
GO TO 755 FORT1967
C FORT1968
C SECOND EQUATION IS THE PIVOT THIS TIME. CASE 2. FORT1969
C FORT1970
760 B(J,2) = B(J,1) FORT1971
B(J,3) = A(JUMP) - AROOT FORT1972
B(J,4) = B(J+1,1) FORT1973
TEMP = 1.D0 FORT1974
IF (DABS(B(J,1)).GT.THETA1) TEMP = ELIM1/B(J,1) FORT1975
ELIM1 = ELIM2 - TEMP*B(J,3) FORT1976
ELIM2 = -TEMP*B(J+1,1) FORT1977
C FORT1978
C SAVE FACTOR FOR THE SECOND ITERATION. FORT1979
C FORT1980
755 B(J,5) = TEMP FORT1981
750 CONTINUE FORT1982
B(N,2) = ELIM1 FORT1983
B(N,3) = 0.D0 FORT1984
B(N,4) = 0.D0 FORT1985
B(NM1,4) = 0.D0 FORT1986
ITER = 1 FORT1987
C FORT1988
C BACK SUBSTITUTE TO GET THIS VECTOR. FORT1989
C FORT1990
790 L = N + 1 FORT1991
DO 780 J=1,N FORT1992
L = L - 1 FORT1993
786 CONTINUE FORT1994
ELIM1 = VECT(L,I)-VECT(L+1,I)*B(L,3)-VECT(L+2,I)*B(L,4) FORT1995
C FORT1996
C IF OVERFLOW IS CONCEIVABLE, SCALE THE VECTOR DOWN. THIS APPROACH FORT1997
C IS USED TO AVOID MACHINE-DEPENDENT AND SYSTEM-DEPENDENT CALLS TO FORT1998
C OVERFLOW ROUTINES. FORT1999
C FORT2000
IF (DABS(ELIM1).GT.DELBIG) GO TO 782 FORT2001
TEMP = B(L,2) FORT2002
IF (DABS(B(L,2)).LT.DELTA) TEMP = DELTA FORT2003
VECT(L,I) = ELIM1/TEMP FORT2004
GO TO 780 FORT2005
782 DO 784 K=1,N FORT2006
784 VECT(K,I) = VECT(K,I)/DELBIG FORT2007
GO TO 786 FORT2008
780 CONTINUE FORT2009
GO TO (820,900), ITER FORT2010
C FORT2011
C SECOND ITERATION. FORT2012
C FORT2013
820 ITER = ITER + 1 FORT2014
890 ELIM1 = VECT(1,I) FORT2015
DO 830 J=1,NM1 FORT2016
IF (B(J,2).EQ.B(J,1)) GO TO 840 FORT2017
C FORT2018
C CASE ONE. FORT2019
C FORT2020
850 VECT(J,I) = ELIM1 FORT2021
ELIM1 = VECT(J+1,I) - ELIM1*B(J,5) FORT2022
GO TO 830 FORT2023
C FORT2024
C CASE TWO. FORT2025
C FORT2026
840 VECT(J,I) = VECT(J+1,I) FORT2027
ELIM1 = ELIM1 - VECT(J+1,I)*TEMP FORT2028
830 CONTINUE FORT2029
VECT(N,I) = ELIM1 FORT2030
GO TO 790 FORT2031
C FORT2032
C NORMALIZE THE VECTOR. FORT2033
C FORT2034
900 ELIM1 = 0.D0 FORT2035
DO 904 J=1,N FORT2036
904 ELIM1 = DMAX1(DABS(VECT(J,I)),ELIM1) FORT2037
TEMP = 0.D0 FORT2038
DO 910 J=1,N FORT2039
ELIM2 = VECT(J,I)/ELIM1 FORT2040
910 TEMP = TEMP + ELIM2**2 FORT2041
TEMP = 1.D0/(DSQRT(TEMP)*ELIM1) FORT2042
DO 920 J=1,N FORT2043
VECT(J,I) = VECT(J,I)*TEMP FORT2044
IF (DABS(VECT(J,I)).LT.DEL1) VECT(J,I) = 0.D0 FORT2045
920 CONTINUE FORT2046
700 CONTINUE FORT2047
C FORT2048
C ROTATE THE CODIAGONAL VECTORS INTO VECTORS OF ORIGINAL ARRAY. FORT2049
C LOOP OVER ALL THE TRANSFORMATION VECTORS. FORT2050
C FORT2051
940 IF (NM2.EQ.0) GO TO 1002 FORT2052
JUMP = NSIZE - NP1 FORT2053
IM = NM1 FORT2054
DO 950 I=1,NM2 FORT2055
LIMIT = IDINT(B(IM-1,6)) FORT2056
J1 = JUMP FORT2057
C FORT2058
C MOVE A TRANSFORMATION VECTOR OUT INTO BETTER INDEXING POSITION. FORT2059
C FORT2060
DO 955 J=IM,LIMIT FORT2061
B(J,2) = A(J1) FORT2062
955 J1 = J1 + J FORT2063
IDIF = LIMIT - IM FORT2064
C FORT2065
C MODIFY ALL REQUESTED VECTORS. FORT2066
C FORT2067
DO 960 K=1,NROOT FORT2068
C FORT2069
C FORM SCALAR PRODUCT OF TRANSFORMATION VECTOR WITH EIGENVECTOR. FORT2070
C FORT2071
TMP = DOT(B(IM,2),VECT(IM,K)) FORT2072
FACT = -TMP - TMP FORT2073
CALL VECSUM(VECT(IM,K),B(IM,2)) FORT2074
960 CONTINUE FORT2075
JUMP = JUMP - IM FORT2076
950 IM = IM - 1 FORT2077
1002 CONTINUE FORT2078
C FORT2079
C RESTORE ROOTS TO THEIR PROPER SIZE AND ADD BACK TRACE. FORT2080
C FORT2081
DO 95 I=1,N FORT2082
95 ROOT(I) = ROOT(I)*ANORM + TRACE FORT2083
1001 RETURN FORT2084
END FORT2085
SUBROUTINE EVQR(A,B,N,M,TOL) FORT2086
C FORT2087
C SUBROUTINE FOR PERFORMING A QR TRANSFORM ON A REAL SYMMETRIC FORT2088
C TRIDIAGONAL MATRIX. FORT2089
C FORT2090
C ON ENTERING, A CONTAINS THE N DIAGONAL ELEMENTS OF THE TRIDIAGONAL FORT2091
C MATRIX. B CONTAINS THE SQUARES OF THE N-1 OFF-DIAGONAL ELEMENTS. FORT2092
C ITERATION IS CONTINUED UNTIL THE SQUARES OF THE OFF-DIAGONAL FORT2093
C ELEMENTS ARE LESS THAN TOL. TYPICALLY LESS THAN TWO ITERATIONS PER FORT2094
C EIGENVALUE ARE REQUIRED. THUS THE UPPER LIMIT M TO THE NUMBER OF FORT2095
C ITERATIONS PER EIGENVALUE MAY BE SAFELY SET AT 20 OR SO. FORT2096
C FORT2097
IMPLICIT REAL*8(A-H,O-Z) FORT2098
DIMENSION A(1),B(1) FORT2099
NX = N FORT2100
NI = 1 FORT2101
SH = 0.D0 FORT2102
C FORT2103
C K COUNTS THE NUMBER OF ITERATIONS PER EIGENVALUE. FORT2104
C FORT2105
K = 0 FORT2106
IF (NX-2) 50,60,85 FORT2107
C FORT2108
C EACH NEW ITERATION BEGINS HERE. FORT2109
C FORT2110
100 K = K + 1 FORT2111
IF (K.LT.M) GO TO 101 FORT2112
WRITE (6,1000) K FORT2113
1000 FORMAT(35H NO CONVERGENCE OF QR ALGORITHM IN ,I4,11H ITERATIONS)FORT2114
CALL EXIT FORT2115
C FORT2116
C SOLVE THE TWO BY TWO IN THE LOWER RIGHT CORNER AND USE SMALLER ROOT FORT2117
C AS THE NEXT SHIFT. FORT2118
C FORT2119
101 AT = A(NX) + A(NX-1) FORT2120
ST = AT*5.D-1 FORT2121
DISC = AT**2 - 4.D0*(A(NX)*A(NX-1)-B(NX-1)) FORT2122
IF (DISC.LE.0.D0) GO TO 15 FORT2123
ST = ST - DSIGN(DSQRT(DISC),ST)*5.D-1 FORT2124
C FORT2125
C INCREASE THE TOTAL SHIFT BY THE TEMPORARY SHIFT. FORT2126
C FORT2127
15 SH = SH + ST FORT2128
C FORT2129
C THIS LOOP SUBTRACTS THE TEMPORARY SHIFT FROM THE DIAGONAL ELEMENTS. FORT2130
C FORT2131
DO 20 I=1,NX FORT2132
20 A(I) = A(I) - ST FORT2133
C FORT2134
C INITIALIZE. FORT2135
C FORT2136
G = A(NI) FORT2137
PS = G**2 FORT2138
RS = PS + B(NI) FORT2139
SX = B(NI)/RS FORT2140
CXS = 1.D0 FORT2141
CX = PS/RS FORT2142
U = SX*(G+A(NI+1)) FORT2143
A(NI) = G + U FORT2144
NTOP = NX - 2 FORT2145
C FORT2146
C THIS LOOP COMPLETES ONE ITERATION, THAT IS ONE QR TRANSFORM. FORT2147
C FORT2148
DO 10 I=NI,NTOP FORT2149
C FORT2150
C G IS THE GAMMA IN THE NOTATION OF WILKINSON. FORT2151
C FORT2152
G = A(I+1) - U FORT2153
IF (CX.GT.TOL) GO TO 12 FORT2154
PS = B(I)*CXS FORT2155
GO TO 16 FORT2156
12 PS = G**2/CX FORT2157
16 RS = PS + B(I+1) FORT2158
C FORT2159
C ROTATE AN OFF-DIAGONAL ELEMENT. FORT2160
C FORT2161
B(I) = SX*RS FORT2162
SX = B(I+1)/RS FORT2163
CXS = CX FORT2164
CX = PS/RS FORT2165
U = SX*(G+A(I+2)) FORT2166
C FORT2167
C ROTATE A DIAGONAL ELEMENT. FORT2168
C FORT2169
A(I+1) = G + U FORT2170
10 CONTINUE FORT2171
C FORT2172
C COMPUTE THE LAST DIAGONAL ELEMENT. FORT2173
C FORT2174
A(NX) = A(NX) - U FORT2175
C FORT2176
C COMPUTE THE LAST OFF-DIAGONAL ELEMENT. FORT2177
C FORT2178
IF (CX.GT.TOL) GO TO 112 FORT2179
PS = B(NTOP+1)*CXS FORT2180
GO TO 116 FORT2181
112 PS = ((A(NX))**2)/CX FORT2182
116 B(NTOP+1) = SX*PS FORT2183
C FORT2184
C END OF ONE ITERATION. FORT2185
C FORT2186
85 IT = NX FORT2187
C FORT2188
C CHECK UPWARD THROUGH THE OFF-DIAGONAL ELEMENTS TO FIND THOSE LESS FORT2189
C THAN TOL. IF NO OFF-DIAGONAL ELEMENTS LESS THAN TOL ARE FOUND, FORT2190
C PERFORM ANOTHER ITERATION. FORT2191
C FORT2192
30 IT = IT-1 FORT2193
IF (DABS(B(IT)).LE.TOL) GO TO 40 FORT2194
IF (IT-NI) 100,100,30 FORT2195
C FORT2196
C BRANCH ACCORDING TO WHETHER THE MATRIX ISOLATED BY THE SMALL FORT2197
C OFF-DIAGONAL ELEMENT IS OF DIMENSION ONE, TWO, OR MORE. FORT2198
C FORT2199
40 IF (NX-IT-2) 50,60,70 FORT2200
C FORT2201
C EXTRACT THE EIGENVALUE OF A ONE BY ONE MATRIX BY ADDING BACK FORT2202
C THE SHIFT. FORT2203
C FORT2204
50 A(NX) = A(NX) + SH FORT2205
C FORT2206
C DECREASE THE SIZE OF THE PORTION OF THE MATRIX AFFECTED BY FORT2207
C LATER ITERATIONS. FORT2208
C FORT2209
NX = NX-1 FORT2210
C FORT2211
C RESET THE ITERATION COUNTER. FORT2212
C FORT2213
K = 1 FORT2214
GO TO 80 FORT2215
C FORT2216
C EXTRACT THE EIGENVALUES FROM A TWO BY TWO MATRIX. FORT2217
C FORT2218
60 AL = B(NX-1) FORT2219
AM = 5.D-1*(A(NX-1)-A(NX)) FORT2220
AMS = AM**2 FORT2221
SAM = DSIGN(1.D0,AM) FORT2222
AN = DSQRT(AL+AM**2) FORT2223
CX = (AN+DABS(AM))/(2.D0*AN) FORT2224
SX = B(NX-1)/(4.D0*AN**2*CX) FORT2225
TA = A(NX-1) FORT2226
TB = A(NX) FORT2227
TC = B(NX-1) FORT2228
I = NX FORT2229
C FORT2230
C ROTATE THE DIAGONAL ELEMENTS AND THE OFF-DIAGONAL ELEMENTS. FORT2231
C FORT2232
A(NX-1) = TA*CX+TB*SX+TC*SAM/AN+SH FORT2233
A(NX) = TA*SX+TB*CX-TC*SAM/AN+SH FORT2234
B(NX-1) = 4.D0*AMS*CX*SX-DABS(AM)*TC/AN+TC*(CX-SX)**2 FORT2235
C FORT2236
C RESET THE ITERATION COUNTER. FORT2237
C FORT2238
K = 1 FORT2239
C FORT2240
C DECREASE THE SIZE OF THE PORTION OF THE MATRIX AFFECTED BY THE LATER FORT2241
C ITERATIONS. FORT2242
C FORT2243
NX = NX-2 FORT2244
GO TO 80 FORT2245
C FORT2246
C THE NEXT STATEMENT IS REACHED WHEN THE PORTION OF THE MATRIX FORT2247
C ISOLATED IS GREATER THAN TWO BY TWO. IT CHANGES THE LOWER LIMIT FORT2248
C OF THE ITERATION SO THAT ONLY THIS PORTION WILL BE AFFECTED BY FORT2249
C SUBSEQUENT ROTATIONS UNTIL ALL ITS EIGENVALUES ARE FOUND. FORT2250
C FORT2251
70 NI = IT + 1 FORT2252
C FORT2253
C TRANSFER TO BEGINNING OF ANOTHER ITERATION. FORT2254
C FORT2255
GO TO 85 FORT2256
C FORT2257
C NEXT STATEMENT IS REACHED AFTER EITHER ONE OR TWO EIGENVALUES FORT2258
C HAVE JUST BEEN FOUND. IT TRANSFERS IF ALL THE EIGENVALUES IN FORT2259
C THIS PORTION OF THE MATRIX HAVE BEEN FOUND. FORT2260
C FORT2261
80 IF (NX.LT.NI) GO TO 90 FORT2262
C FORT2263
C BRANCH ACCORDING TO WHETHER ONE, TWO, OR MORE EIGENVALUES REMAIN FORT2264
C TO BE FOUND. FORT2265
C FORT2266
95 IF (NX-NI-1) 50,60,85 FORT2267
C FORT2268
C THE NEXT STATEMENT IS REACHED WHEN ALL EIGENVALUES IN THIS PART FORT2269
C OF THE MATRIX HAVE BEEN FOUND. IT RETURNS IF THIS IS THE LAST FORT2270
C PART OF THE MATRIX. FORT2271
C FORT2272
90 IF (NI.EQ.1) RETURN FORT2273
C FORT2274
C ENLARGE THE PORTION OF THE MATRIX BEING TREATED TO INCLUDE THE FORT2275
C BEGINNING OF THE MATRIX. FORT2276
C FORT2277
NI = 1 FORT2278
GO TO 95 FORT2279
END FORT2280
SUBROUTINE QRTN (A,B,V,EIG,N,M,TOL,NJX) FORT2281
C FORT2282
C SUBROUTINE FOR PERFORMING A QR TRANSFORM ON A REAL SYMMETRIC FORT2283
C TRIDIAGONAL MATRIX. FORT2284
C FORT2285
C N IS THE DIMENSION OF THE MATRIX. A CONTAINS THE N DIAGONAL FORT2286
C ELEMENTS OF THE TRIDIAGONAL MATRIX. B CONTAINS THE N-1 OFF- FORT2287
C DIAGONAL ELEMENTS. M IS THE MAXIMUM NUMBER OF ITERATIONS ( SAY FORT2288
C 20 ). NJX IS THE PHYSICAL ROW DIMENSION OF THE EIGENVECTOR FORT2289
C MATRIX V. FORT2290
C FORT2291
C THE EIGENVALUES ARE ASSUMED KNOWN AND PLACED IN THE FIRST N FORT2292
C ELEMENTS OF EIG. FORT2293
C FORT2294
C N VECTORS V ( EACH OF LENGTH N ) ARE TRANSFORMED INTO THE FORT2295
C BASIS IN WHICH THE MATRIX IS DIAGONAL. FORT2296
C FORT2297
IMPLICIT REAL*8 (A-H,O-Z) FORT2298
DIMENSION A(1),B(1),V(1),EIG(1) FORT2299
NX = N FORT2300
NNM1 = NJX*(NX-1) FORT2301
NI = 1 FORT2302
C FORT2303
C SET INITIAL TOTAL SHIFT. FORT2304
C FORT2305
SH = 0 FORT2306
IF (NX-2) 50,60,1 FORT2307
C FORT2308
C K COUNTS THE NUMBER OF ITERATIONS PER EIGENVALUE. FORT2309
C FORT2310
1 K=0 FORT2311
C FORT2312
C SET INITIAL TEMPORARY SHIFT. FORT2313
C FORT2314
98 ST = EIG(NX)-SH FORT2315
C FORT2316
C CHECK FOR SMALL OFF-DIAGONAL ELEMENTS. FORT2317
C FORT2318
IT = NX FORT2319
99 IT = IT-1 FORT2320
IF (DABS(B(IT)).LE.TOL) GO TO 40 FORT2321
IF (IT.GT.NI) GO TO 99 FORT2322
C FORT2323
C NO SMALL OFF-DIAGONAL ELEMENTS FOUND. ITERATE. EACH NEW FORT2324
C ITERATION BEGINS HERE. FORT2325
C FORT2326
100 K = K + 1 FORT2327
IF (K.LT.M) GO TO 11 FORT2328
WRITE (6,1000) K FORT2329
1000 FORMAT(35H NO CONVERGENCE OF QR ALGORITHM IN ,I4,11H ITERATIONS)
CALL EXIT FORT2331
11 IF (K.EQ.1) GO TO 15 FORT2332
C FORT2333
C SOLVE THE TWO BY TWO IN THE LOWER RIGHT CORNER AND USE SMALLER ROOT FORT2334
C AS THE NEXT SHIFT. FORT2335
C FORT2336
12 AT = A(NX) + A(NX-1) FORT2337
ST = AT*5.D-1 FORT2338
DISC = AT**2-4.D0*(A(NX)*A(NX-1)-B(NX-1)**2) FORT2339
IF (DISC.LE.0.D0) GO TO 15 FORT2340
ST = ST-DSIGN(DSQRT(DISC),ST)*5.D-1 FORT2341
C FORT2342
C INCREASE THE TOTAL SHIFT BY THE TEMPORARY SHIFT. FORT2343
C FORT2344
15 SH = SH + ST FORT2345
C FORT2346
C THIS LOOP SUBTRACTS THE TEMPORARY SHIFT FROM THE DIAGONAL ELEMENTS. FORT2347
C FORT2348
DO 20 I=1,NX FORT2349
20 A(I) = A(I) - ST FORT2350
R = DSQRT(A(NI)**2+B(NI)**2) FORT2351
S = B(NI)/R FORT2352
CS = S FORT2353
C = A(NI)/R FORT2354
U = (S**2)*(A(NI)+A(NI+1)) FORT2355
A(NI) = A(NI) + U FORT2356
CALL ROTATE(V(NI),C,S,NJX,NNM1) FORT2357
NTOP = NX - 2 FORT2358
C FORT2359
C THIS LOOP COMPLETES ONE ITERATION, THAT IS ONE QR TRANSFORM. FORT2360
C FORT2361
DO 10 I=NI,NTOP FORT2362
C FORT2363
C G IS THE GAMMA AND Q IS THE P IN THE NOTATION OF WILKINSON. FORT2364
C FORT2365
G = A(I+1) - U FORT2366
Q = C*A(I+1) - CS*B(I) FORT2367
R = DSQRT(Q**2+B(I+1)**2) FORT2368
C FORT2369
C ROTATE AN OFF-DIAGONAL ELEMENT. FORT2370
C FORT2371
B(I) = S*R FORT2372
C FORT2373
C FIND THE NEW SINE AND COSINE FOR THE JACOBI ROTATION. THEN FORT2374
C COMPUTE A NEW U. FORT2375
C FORT2376
S = B(I+1)/R FORT2377
CS = C*S FORT2378
C = Q/R FORT2379
U = (S**2)*(G+A(I+2)) FORT2380
C FORT2381
C ROTATE A DIAGONAL ELEMENT. FORT2382
C FORT2383
A(I+1) = G + U FORT2384
C FORT2385
C ROTATE THE VECTORS. FORT2386
C FORT2387
CALL ROTATE (V(I+1),C,S,NJX,NNM1) FORT2388
10 CONTINUE FORT2389
C FORT2390
C COMPUTE THE LAST OFF DIAGONAL ELEMENT. FORT2391
C FORT2392
B(NTOP+1) = S*(C*A(NX)-CS*B(NTOP+1)) FORT2393
C FORT2394
C COMPUTE THE LAST DIAGONAL ELEMENT. FORT2395
C FORT2396
A(NX) = A(NX) - U FORT2397
C FORT2398
C END OF ONE ITERATION. FORT2399
C FORT2400
85 IT = NX FORT2401
C FORT2402
C CHECK UPWARD THROUGH THE OFF DIAGONAL ELEMENTS TO FIND THOSE LESS FORT2403
C THAN TOL. IF NO OFF-DIAGONAL ELEMENTS LESS THAN TOL ARE FOUND, FORT2404
C PERFORM ANOTHER ITERATION. FORT2405
C FORT2406
30 IT = IT - 1 FORT2407
IF (DABS(B(IT)).LE.TOL) GO TO 40 FORT2408
IF (IT-NI) 100,100,30 FORT2409
C FORT2410
C BRANCH ACCORDING TO WHETHER THE MATRIX ISOLATED BY THE SMALL FORT2411
C OFF-DIAGONAL ELEMENT IS OF DIMENSION ONE, TWO, OR MORE. FORT2412
C FORT2413
40 IF (NX-IT-2) 50,60,70 FORT2414
C FORT2415
C EXTRACT THE EIGENVALUE OF A ONE BY ONE MATRIX BY ADDING BACK FORT2416
C THE SHIFT. FORT2417
C FORT2418
50 A(NX) = A(NX) + SH FORT2419
C FORT2420
C DECREASE THE SIZE OF THE PORTION OF THE MATRIX AFFECTED BY FORT2421
C LATER ITERATIONS. FORT2422
C FORT2423
NX = NX - 1 FORT2424
C FORT2425
C RESET THE ITERATION COUNTER. FORT2426
C FORT2427
K = 0 FORT2428
GO TO 80 FORT2429
C FORT2430
C EXTRACT THE EIGENVALUES FROM A TWO BY TWO MATRIX AND PERFORM THE FORT2431
C CORRESPONDING ROTATIONS ON THE VECTORS. FORT2432
C FORT2433
60 AL = -B(NX-1) FORT2434
AM = 5.D-1*(A(NX-1)-A(NX)) FORT2435
AN = DSQRT(AL**2+AM**2) FORT2436
C = DSQRT((AN+DABS(AM))/(2.D0*AN)) FORT2437
S = DSIGN(5.D-1,AM)*AL/(AN*C) FORT2438
TA = A(NX-1) FORT2439
TB = A(NX) FORT2440
TC = B(NX-1) FORT2441
CX = C**2 FORT2442
SX = S**2 FORT2443
CS = C*S FORT2444
C FORT2445
C ROTATE THE DIAGONAL ELEMENTS, THE OFF-DIAGONAL ELEMENTS, AND FORT2446
C THE VECTORS. FORT2447
C FORT2448
A(NX-1) = TA*CX+TB*SX-2.D0*TC*CS+SH FORT2449
A(NX) = TA*SX+TB*CX+2.D0*TC*CS+SH FORT2450
B(NX-1) = 2.D0*AM*CS+TC*(CX-SX) FORT2451
I = NX-1 FORT2452
S = -S FORT2453
CALL ROTATE (V(I),C,S,NJX,NNM1) FORT2454
C FORT2455
C RESET THE ITERATION COUNTER. FORT2456
C FORT2457
K = 0 FORT2458
C FORT2459
C DECREASE THE SIZE OF THE PORTION OF THE MATRIX AFFECTED BY THE FORT2460
C LATER ITERATIONS. FORT2461
C FORT2462
NX = NX-2 FORT2463
GO TO 80 FORT2464
C FORT2465
C THE NEXT STATEMENT IS REACHED WHEN THE PORTION OF THE MATRIX FORT2466
C ISOLATED IS GREATER THAN TWO BY TWO. IT CHANGES THE LOWER LIMIT FORT2467
C OF THE ITERATION SO THAT ONLY THIS PORTION WILL BE AFFECTED BY FORT2468
C SUBSEQUENT ROTATIONS UNTIL ALL ITS EIGENVALUES ARE FOUND. FORT2469
C FORT2470
70 NI = IT + 1 FORT2471
C FORT2472
C TRANSFER TO BEGINNING OF ANOTHER ITERATION. FORT2473
C FORT2474
GO TO 100 FORT2475
C FORT2476
C NEXT STATEMENT IS REACHED AFTER EITHER ONE OR TWO EIGENVALUES FORT2477
C HAVE JUST BEEN FOUND. IT TRANSFERS IF ALL THE EIGENVALUES IN FORT2478
C THIS PORTION OF THE MATRIX HAVE BEEN FOUND. FORT2479
C FORT2480
80 IF (NX.LT.NI) GO TO 90 FORT2481
C FORT2482
C BRANCH ACCORDING TO WHETHER ONE, TWO, OR MORE EIGENVALUES REMAIN FORT2483
C TO BE FOUND. FORT2484
C FORT2485
95 IF (NX-NI-1) 50,60,98 FORT2486
C FORT2487
C THE NEXT STATEMENT IS REACHED WHEN ALL EIGENVALUES IN THIS PART FORT2488
C OF THE MATRIX HAVE BEEN FOUND. IT RETURNS IF THIS IS THE LAST FORT2489
C PART OF THE MATRIX. FORT2490
C FORT2491
90 IF (NI.EQ.1) RETURN FORT2492
C FORT2493
C ENLARGE THE PORTION OF THE MATRIX BEING TREATED TO INCLUDE THE FORT2494
C BEGINNING OF THE MATRIX. FORT2495
C FORT2496
NI = 1 FORT2497
GO TO 95 FORT2498
END FORT2499
SUBROUTINE ITRATE(H,U,MAD,C,E,W,IOCC,HDG,NDIM,NTYPE,NC,NHDG) FORT2500
C FORT2501
C SUBROUTINE TO SETUP HUCKEL MATRIX WHEN USING CHARGE ITERATION FORT2502
C OPTION ( METH = 2 OR 3 ). FORT2503
C FORT2504
IMPLICIT REAL*8(A-H,O-Z) FORT2505
DIMENSION H(NDIM,NDIM),U(NDIM,NDIM),MAD(NTYPE,NTYPE),C(NC), FORT2506
1 E(NDIM),W(NDIM),IOCC(NDIM),HDG(NHDG) FORT2507
REAL*8 MAD FORT2508
REAL*4 IOCC FORT2509
PARAMETER (MAXATM=500)
PARAMETER (BB=250)
PARAMETER (MXUSER=230)
PARAMETER (MXUSR2=231)
COMMON/TITLE/AB(10) FORT2510
COMMON/CNTRL/CON,PEEP,COULH,NH,NA,NATM,KA,NELEC,METH,IPRINT, FORT2511
1 IPUNCH,L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT2512
LOGICAL*1 L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT2513
COMMON/OUT/PRT(20),PUN(20) FORT2514
LOGICAL*1 PRT,PUN FORT2515
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 FORT2519
COMMON/ITPARM/DAMP1,DAMP2,DAMP3,LAMPRI,DELTAC,SENSE,MAXCYC,PRTCYC,
1 ICYCLE,NCON,PARTIT,PRINTX,ITABLE(20) FORT2521
REAL*8 LAMPRI FORT2522
INTEGER*4 PRTCYC FORT2523
LOGICAL*1 PARTIT,PRINTX,ITABLE FORT2524
C FORT2525
C SINCE THE INTERNAL ATOMIC PARAMETERS ( EXPS,EXPP, ETC. ) ARE FORT2526
C NOT USED WHEN DOING CHARGE ITERATION, THE SPACE ALLOCATED TO FORT2527
C THEM HAS BEEN USED FOR THE VSIE PARAMETERS. FORT2528
C FORT2529
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), FORT2535
1 BD2(5),CD2(5),AS3(5),BS3(5),CS3(5),AP3(5),BP3(5),CP3(5),AD3(5), FORT2536
2 BD3(5),CD3(5) FORT2537
REAL*4 ZS(MXUSER),ZP(MXUSER),ZD(MXUSER)
EQUIVALENCE (ZS(1),NS(MXUSR2)),(ZP(1),NP(MXUSR2)),
. (ZD(1),ND(MXUSR2))
DIMENSION IL(3),JL(3) FORT2540
EQUIVALENCE (PEEP,ADJUST),(COULH,SIGNN),(SENSE,DENSE), FORT2541
1 (VELEC(81),DCHG),(PRT(1),NMIN),(PUN(1),NLN) FORT2542
DIMENSION F(6) FORT2543
EQUIVALENCE (A1,F(1)),(A2,F(3)),(A3,F(5)) FORT2544
REAL*8 LAMBDA FORT2545
QON=DFLOAT(KA)/DFLOAT(NELEC) FORT2546
IF(ICYCLE.GT.1) GO TO 696 FORT2547
WRITE(6,750) FORT2548
750 FORMAT('1') FORT2549
DO 5 I=1,NDIM FORT2550
NMIN=I FORT2551
IF(IOCC(I).GT.0.0001) GO TO 697 FORT2552
5 CONTINUE FORT2553
C FORT2554
C IF CHARGE ITERATION WITH MADELUNG CORRECTION ( METH>2 ) IS BEING FORT2555
C USED, COMPUTE THE TOTAL ORBITAL OCCUPATIONS FOR THE FREE ATOMS. FORT2556
C FILL ATOMIC ORBITALS IN ORDER0 1S 2S 2P 3S 3P 3D 4S 4P 4D 5S 5P FORT2557
C 5D 6S 6P. FORT2558
C FORT2559
697 IF(METH.LT.3) GO TO 696 FORT2560
PRTCYC=MAXCYC
DO 700 I=1,NA FORT2561
KEYI=KEY(I) FORT2562
J=NP(KEYI) FORT2563
IF(J.EQ.0) GO TO 701 FORT2564
IL(2)=J*10+1 FORT2565
JL(2)=2 FORT2566
J=NS(KEYI) FORT2567
IL(1)=J*10 FORT2568
JL(1)=1 FORT2569
J=ND(KEYI) FORT2570
IF(J.EQ.0) GO TO 702 FORT2571
IL(3)=J*10+2 FORT2572
JL(3)=3 FORT2573
DO 710 J=2,3 FORT2574
IJ=4-J FORT2575
K=IJ+1 FORT2576
DO 710 L=1,IJ FORT2577
N=K-L FORT2578
IF(IL(K).GT.IL(N)) GO TO 710 FORT2579
I1=IL(K) FORT2580
J1=JL(K) FORT2581
IL(K)=IL(N) FORT2582
JL(K)=JL(N) FORT2583
IL(N)=I1 FORT2584
JL(N)=J1 FORT2585
710 CONTINUE FORT2586
GO TO 720 FORT2587
701 JL(1)=1 FORT2588
JL(2)=2 FORT2589
JL(3)=3 FORT2590
GO TO 720 FORT2591
702 JL(3)=3 FORT2592
IF(IL(1).LT.IL(2)) GO TO 720 FORT2593
JL(2)=1 FORT2594
JL(1)=2 FORT2595
720 L=VELEC(KEYI) FORT2596
DO 700 J=1,3 FORT2597
K=JL(J) FORT2598
J1=4*(K-1)+2 FORT2599
IF(L.LT.J1) J1=L FORT2600
L=L-J1 FORT2601
IF(K-2) 725,726,727 FORT2602
725 ZS(KEYI)=FLOAT(J1) FORT2603
GO TO 700 FORT2604
726 ZP(KEYI)=FLOAT(J1) FORT2605
GO TO 700 FORT2606
727 ZD(KEYI)=FLOAT(J1) FORT2607
700 CONTINUE FORT2608
696 IF(PRINTX) WRITE(6,751) FORT2609
751 FORMAT(////) FORT2610
PRINTX=((ICYCLE/PRTCYC)*PRTCYC.EQ.ICYCLE) FORT2611
IF(ICYCLE.EQ.MAXCYC) PRINTX=.TRUE. FORT2612
C FORT2613
C CALCULATE SUM OF ONE-ELECTRON ENERGIES. FORT2614
C FORT2615
753 SUM=0.0D0 FORT2616
DO 13 I=1,NDIM FORT2617
W(I)=DBLE(IOCC(I)) FORT2618
13 SUM=SUM+E(I)*W(I) FORT2619
C FORT2620
C COMPUTE ATOMIC ORBITAL OCCUPATIONS AND STORE IN E(I). FORT2621
C FORT2622
3004 IJ=1 FORT2623
DO 319 I=1,NDIM FORT2624
E(I)=0.0D0 FORT2625
DO 319 J=1,I FORT2626
UB=0.0D0 FORT2627
DO 320 K=NMIN,NDIM FORT2628
320 UB=UB+H(I,K)*H(J,K)*W(K) FORT2629
UB=UB*0.50D0 FORT2630
IF(I.EQ.J) GO TO 28 FORT2631
UB=(UB+UB)*C(IJ) FORT2632
IJ=IJ+1 FORT2633
28 E(I)=E(I)+UB FORT2634
E(J)=E(J)+UB FORT2635
319 CONTINUE FORT2636
C FORT2637
C COMPUTE THE ORBITAL OCCUPATION OF A GIVEN TYPE (S,P,D) WHICH FORT2638
C VARIES MOST FROM THE LAST CYCLE. FORT2639
C FORT2640
3005 DENOM=0.0D0 FORT2641
DENSE2=DENSE FORT2642
DCHG2=DCHG FORT2643
J=1 FORT2644
DO 5000 I=1,NA FORT2645
KEYI=KEY(I) FORT2646
SDENSE=E(J) FORT2647
SIGNS=SDENSE-X(I) FORT2648
C(J)=SIGNS FORT2649
DIFFS=DABS(SIGNS) FORT2650
N=J FORT2651
IF(DENOM.GT.DIFFS) GO TO 5001 FORT2652
DENSE=X(I) FORT2653
DCHG=SIGNS FORT2654
DENOM=DIFFS FORT2655
NLF=10*I FORT2656
5001 IF(NP(KEYI).EQ.0) GO TO 5000 FORT2657
PDENSE=E(J+1)+E(J+2)+E(J+3) FORT2658
SIGNP=PDENSE-Y(I) FORT2659
C(J+3)=SIGNP FORT2660
DIFFP=DABS(SIGNP) FORT2661
N=J+3 FORT2662
IF(DENOM.GT.DIFFP) GO TO 5002 FORT2663
DENSE=Y(I) FORT2664
DCHG=SIGNP FORT2665
DENOM=DIFFP FORT2666
NLF=10*I+1 FORT2667
5002 IF(ND(KEYI).EQ.0) GO TO 5000 FORT2668
DDENSE=E(J+4)+E(J+5)+E(J+6)+E(J+7)+E(J+8) FORT2669
SIGND=DDENSE-Z(I) FORT2670
C(J+8)=SIGND FORT2671
DIFFD=DABS(SIGND) FORT2672
N=J+8 FORT2673
IF(DENOM.GT.DIFFD) GO TO 5000 FORT2674
DENSE=Z(I) FORT2675
DCHG=SIGND FORT2676
DENOM=DIFFD FORT2677
NLF=10*I+2 FORT2678
5000 J=N+1 FORT2679
SIGNF=DCHG/DENOM FORT2680
C FORT2681
C DETERMINE LAMBDA, THE DAMPING FACTOR. FORT2682
C FORT2683
IF(ICYCLE.LE.2) SIGNN=SIGNF FORT2684
IF(ICYCLE.LE.2) NLN=NLF FORT2685
ISIGNN=SIGNN FORT2686
ISIGNF=SIGNF FORT2687
IF(ICYCLE.EQ.2) ADJUST=DAMP1*DENOM FORT2688
IF(ISIGNN.EQ.ISIGNF.OR.NLN.NE.NLF) GO TO 8001 FORT2689
IF(DAMP3.NE.0.0D0) GO TO 401 FORT2690
LAMBDA=(DENSE2-DENSE)/(DCHG-DCHG2) FORT2691
GO TO 402 FORT2692
401 ADJUST=DAMP3*ADJUST FORT2693
8001 IF((ADJUST/DENOM).GE.LAMPRI) ADJUST=DAMP2*DENOM FORT2694
IF(ICYCLE.EQ.1) ADJUST=DENOM FORT2695
LAMBDA=ADJUST/DENOM FORT2696
402 SIGNN=SIGNF FORT2697
NLN=NLF FORT2698
C FORT2699
C FORT2700
C IN THIS SECTION OF THE PROGRAM THREE CALCULATIONS ARE PERFORMED0 FORT2701
C FORT2702
C 1. THE TOTAL ORBITAL OCCUPATIONS OF A GIVEN TYPE (S,P,D) ARE FORT2703
C DAMPED AND STORED (IN X(I),Y(I),Z(I) RESPECTIVELY). FORT2704
C 2. THE NET CHARGES ARE CALCULATED AND STORED IN C(I). FORT2705
C 3. -VSIE'S ARE CALCULATED AND STORED IN W(I). FORT2706
C FORT2707
C FORT2708
J=1 FORT2709
DO 803 I=1,NA FORT2710
KEYI=KEY(I) FORT2711
SDENSE=X(I)+LAMBDA*C(J) FORT2712
X(I)=SDENSE FORT2713
UB=SDENSE FORT2714
N=J FORT2715
IF(NP(KEYI).EQ.0) GO TO 804 FORT2716
PDENSE=Y(I)+LAMBDA*C(J+3) FORT2717
Y(I)=PDENSE FORT2718
UB=UB+PDENSE FORT2719
N=J+3 FORT2720
IF(ND(KEYI).EQ.0) GO TO 805 FORT2721
N=J+8 FORT2722
DDENSE=Z(I)+LAMBDA*C(J+8) FORT2723
Z(I)=DDENSE FORT2724
UB=UB+DDENSE FORT2725
GO TO 806 FORT2726
804 Q=VELEC(KEYI)-UB FORT2727
IF(.NOT.PARTIT) GO TO 1111 FORT2728
IF(ITABLE(KEYI)) GO TO 1111 FORT2729
W(J)=COULS(KEYI) FORT2730
GO TO 807 FORT2731
1111 KEYI=MXUSR2-KEY(I)
W(J)=-((AS1(KEYI)*Q+BS1(KEYI))*Q+CS1(KEYI)) FORT2733
GO TO 807 FORT2734
805 Q=VELEC(KEYI)-UB FORT2735
IF(.NOT.PARTIT) GO TO 1113 FORT2736
IF(ITABLE(KEYI)) GO TO 1113 FORT2737
W(J)=COULS(KEYI) FORT2738
W(J+1)=COULP(KEYI) FORT2739
W(J+2)=COULP(KEYI) FORT2740
W(J+3)=COULP(KEYI) FORT2741
GO TO 807 FORT2742
1113 KEYI=MXUSR2-KEY(I)
VSIES1=(AS1(KEYI)*Q+BS1(KEYI))*Q+CS1(KEYI) FORT2744
W(J)=-VSIES1 FORT2745
VSIEP1=(AP1(KEYI)*Q+BP1(KEYI))*Q+CP1(KEYI) FORT2746
W(J+1)=-VSIEP1 FORT2747
W(J+2)=-VSIEP1 FORT2748
W(J+3)=-VSIEP1 FORT2749
GO TO 807 FORT2750
806 Q=VELEC(KEYI)-UB FORT2751
IF(.NOT.PARTIT) GO TO 1115 FORT2752
IF(ITABLE(KEYI)) GO TO 1115 FORT2753
W(J)=COULS(KEYI)
W(J+1)=COULP(KEYI) FORT2754
W(J+2)=COULP(KEYI) FORT2755
W(J+3)=COULP(KEYI) FORT2756
W(J+4)=COULD(KEYI) FORT2757
W(J+5)=COULD(KEYI) FORT2758
W(J+6)=COULD(KEYI) FORT2759
W(J+7)=COULD(KEYI) FORT2760
W(J+8)=COULD(KEYI) FORT2761
GO TO 807 FORT2762
1115 IF(ND(KEYI).EQ.NP(KEYI)) GO TO 1116 FORT2763
IF(NCON.NE.3) GO TO 1116 FORT2764
KEYI=MXUSR2-KEY(I)
VSIES1=(AS1(KEYI)*Q+BS1(KEYI))*Q+CS1(KEYI) FORT2766
VSIES2=(AS2(KEYI)*Q+BS2(KEYI))*Q+CS2(KEYI) FORT2767
VSIES3=(AS3(KEYI)*Q+BS3(KEYI))*Q+CS3(KEYI) FORT2768
W(J)=(SDENSE+PDENSE-2.0D0)*VSIES1+(1.0D0-SDENSE)*VSIES2- FORT2769
*PDENSE*VSIES3 FORT2770
VSIEP1=(AP1(KEYI)*Q+BP1(KEYI))*Q+CP1(KEYI) FORT2771
VSIEP2=(AP2(KEYI)*Q+BP2(KEYI))*Q+CP2(KEYI) FORT2772
VSIEP3=(AP3(KEYI)*Q+BP3(KEYI))*Q+CP3(KEYI) FORT2773
W(J+1)=(SDENSE+PDENSE-2.0D0)*VSIEP1+(1.0D0-PDENSE)*VSIEP2- FORT2774
*SDENSE*VSIEP3 FORT2775
W(J+2)=W(J+1) FORT2776
W(J+3)=W(J+1) FORT2777
VSIED1=(AD1(KEYI)*Q+BD1(KEYI))*Q+CD1(KEYI) FORT2778
VSIED2=(AD2(KEYI)*Q+BD2(KEYI))*Q+CD2(KEYI) FORT2779
VSIED3=(AD3(KEYI)*Q+BD3(KEYI))*Q+CD3(KEYI) FORT2780
W(J+4)=(SDENSE+PDENSE-1.0D0)*VSIED1-SDENSE*VSIED2-PDENSE*VSIED3 FORT2781
W(J+5)=W(J+4) FORT2782
W(J+6)=W(J+4) FORT2783
W(J+7)=W(J+4) FORT2784
W(J+8)=W(J+4) FORT2785
GO TO 807 FORT2786
1116 KEYI=MXUSR2-KEY(I)
VSIES1=(AS1(KEYI)*Q+BS1(KEYI))*Q+CS1(KEYI) FORT2788
W(J)=-VSIES1 FORT2789
VSIEP1=(AP1(KEYI)*Q+BP1(KEYI))*Q+CP1(KEYI) FORT2790
W(J+1)=-VSIEP1 FORT2791
W(J+2)=W(J+1) FORT2792
W(J+3)=W(J+1) FORT2793
VSIED1=(AD1(KEYI)*Q+BD1(KEYI))*Q+CD1(KEYI) FORT2794
W(J+4)=-VSIED1 FORT2795
W(J+5)=W(J+4) FORT2796
W(J+6)=W(J+4) FORT2797
W(J+7)=W(J+4) FORT2798
W(J+8)=W(J+4) FORT2799
807 J=N+1 FORT2800
C(I)=Q FORT2801
803 CONTINUE FORT2802
C FORT2803
C IF CHARGE ITERATION WITHOUT MADELUNG CORRECTION ( METH=2 ) IS FORT2804
C BEING USED, SETUP HUCKEL MATRIX. OTHERWISE SKIP THIS SECTION. FORT2805
C FORT2806
IF(METH.GT.2) GO TO 999 FORT2807
H(1,1)=W(1) FORT2808
CNST=CON FORT2809
DO 760 I=2,NDIM FORT2810
H(I,I)=W(I) FORT2811
J1=I-1 FORT2812
DO 760 J=1,J1 FORT2813
UB=W(I)+W(J) FORT2814
IF(.NOT.L5) GO TO 761 FORT2815
UC=(W(I)-W(J))/UB FORT2816
UC=UC*UC FORT2817
CNST=CON+UC/2.0D0+UC*UC*(0.5D0-CON) FORT2818
761 UB=CNST*U(I,J)*UB FORT2819
H(J,I)=UB FORT2820
760 H(I,J)=UB FORT2821
GO TO 850 FORT2822
C FORT2823
C IF CHARGE ITERATION WITH MADELUNG CORRECTION ( METH>2 ) IS BEING FORT2824
C USED, SETUP HUCKEL MATRIX. OTHERWISE SKIP THIS SECTION. FORT2825
C FORT2826
999 DGSUM=0.0D0 FORT2827
N=1 FORT2828
M=1 FORT2829
DO 880 I=1,NA FORT2830
KEYI=KEY(I) FORT2831
N1=N FORT2832
IF(NP(KEYI).NE.0) N1=N+1 FORT2833
IF(ND(KEYI).NE.0) N1=N+2 FORT2834
DO 881 J=N,N1 FORT2835
DG1=0.0D0 FORT2836
DG2=0.0D0 FORT2837
L=1 FORT2838
DO 882 K=1,NA FORT2839
KEYK=KEY(K) FORT2840
UB=(X(K)-DBLE(ZS(KEYK)))*MAD(J,L) FORT2841
UC=X(K)*MAD(J,L) FORT2842
L=L+1 FORT2843
IF(NP(KEYK).EQ.0) GO TO 883 FORT2844
UB=UB+(Y(K)-DBLE(ZP(KEYK)))*MAD(J,L) FORT2845
UC=UC+Y(K)*MAD(J,L) FORT2846
L=L+1 FORT2847
IF(ND(KEYK).EQ.0) GO TO 883 FORT2848
UB=UB+(Z(K)-DBLE(ZD(KEYK)))*MAD(J,L) FORT2849
UC=UC+Z(K)*MAD(J,L) FORT2850
L=L+1 FORT2851
883 IF(K.NE.I) DG1=DG1+UB FORT2852
882 DG2=DG2+UC FORT2853
J1=M+2*(J-N) FORT2854
DO 884 L=M,J1 FORT2855
H(L,L)=W(L)+DG1 FORT2856
IF(L5) HDG(L)=H(L,L)+QON*DG2 FORT2857
W(L)=DG2 FORT2858
884 DGSUM=DGSUM+DG2 FORT2859
881 M=J1+1 FORT2860
880 N=N1+1 FORT2861
CNST=CON FORT2862
DO 885 I=2,NDIM FORT2863
J1=I-1 FORT2864
DO 885 J=1,J1 FORT2865
IF(.NOT.L5) GO TO 886 FORT2866
UB=(HDG(I)-HDG(J))/(HDG(I)+HDG(J)) FORT2867
UB=UB*UB FORT2868
CNST=CON+UB/2.0D0+UB*UB*(0.5D0-CON) FORT2869
886 UB=U(I,J)*(CNST*(H(I,I)+H(J,J))-QON*(0.5D0-CNST)*(W(I)+W(J))) FORT2870
H(I,J)=UB FORT2871
885 H(J,I)=UB FORT2872
DGSUM=-(QON*DGSUM)/DFLOAT(NDIM) FORT2873
C FORT2874
C IF DOING LAST CYCLE CALCULATE ENERGY CORRECTIONS. FORT2875
C FORT2876
IF(ICYCLE.NE.MAXCYC) GO TO 850 FORT2877
N=1 FORT2878
UB=0.0D0 FORT2879
UC=0.0D0 FORT2880
DO 887 I=1,NA FORT2881
KEYI=KEY(I) FORT2882
K=N FORT2883
A1=X(I) FORT2884
E(1)=DBLE(ZS(KEYI)) FORT2885
UB=UB-0.5D0*(A1*A1-A1)*MAD(N,N) FORT2886
IF(NP(KEYI).EQ.0) GO TO 888 FORT2887
N=N+1 FORT2888
A2=Y(I) FORT2889
E(3)=DBLE(ZP(KEYI)) FORT2890
UB=UB-0.5D0*(A2*A2-A2)*MAD(N,N)-A1*A2*MAD(N-1,N) FORT2891
IF(ND(KEYI).EQ.0) GO TO 888 FORT2892
N=N+1 FORT2893
A3=Z(I) FORT2894
E(5)=DBLE(ZD(KEYI)) FORT2895
UB=UB-0.5D0*(A3*A3-A3)*MAD(N,N)-A1*A3*MAD(N-2,N)-A2*A3*MAD(N-1,N) FORT2896
888 M=1 FORT2897
I1=I-1 FORT2898
IF(I1.EQ.0) GO TO 887 FORT2899
DO 889 J=1,I1 FORT2900
KEYJ=KEY(J) FORT2901
L=M FORT2902
F(2)=X(J) FORT2903
E(2)=DBLE(ZS(KEYJ)) FORT2904
IF(NP(KEYJ).EQ.0) GO TO 890 FORT2905
M=M+1 FORT2906
F(4)=Y(J) FORT2907
E(4)=DBLE(ZP(KEYJ)) FORT2908
IF(ND(KEYJ).EQ.0) GO TO 890 FORT2909
M=M+1 FORT2910
F(6)=Z(J) FORT2911
E(6)=DBLE(ZD(KEYJ)) FORT2912
890 DO 891 IJ=K,N FORT2913
DO 891 JK=L,M FORT2914
N1=2*(IJ-K)+1 FORT2915
M1=2*(JK-L)+2 FORT2916
891 UC=UC-(F(N1)*F(M1)-E(N1)*E(M1))*MAD(IJ,JK) FORT2917
889 M=M+1 FORT2918
887 N=N+1 FORT2919
A1=UB FORT2920
A2=UC FORT2921
A3=UB+UC FORT2922
C FORT2923
C SAVE MADELUNG TERMS FOR USE IN SUBROUTINE OUTPUT IF DOING FORT2924
C THE LAST CYCLE. FORT2925
C FORT2926
K=1 FORT2927
L=1 FORT2928
DO 792 I=1,NA FORT2929
KEYI=KEY(I) FORT2930
MAD(K,K)=W(L) FORT2931
K=K+1 FORT2932
L=L+1 FORT2933
IF(NP(KEYI).EQ.0) GO TO 792 FORT2934
MAD(K,K)=W(L) FORT2935
K=K+1 FORT2936
L=L+3 FORT2937
IF(ND(KEYI).EQ.0) GO TO 792 FORT2938
MAD(K,K)=W(L) FORT2939
K=K+1 FORT2940
L=L+5 FORT2941
792 CONTINUE FORT2942
C FORT2943
C PRINT OUT RESULTS TO SHOW PROGRESS OF ITERATION PROCEDURE. FORT2944
C FORT2945
850 IF(PRTCYC.GT.0) WRITE(6,793) ICYCLE FORT2946
793 FORMAT('0CYCLE NO.',I3,'0') FORT2947
IF(PRTCYC.LT.0) WRITE(6,794) FORT2948
794 FORMAT('0CONVERGENCE REACHED - FINAL CYCLE FOLLOWS0',///) FORT2949
PRTCYC=IABS(PRTCYC) FORT2950
J=NLF/10 FORT2951
K=NLF-10*J FORT2952
WRITE(6,795) SUM,LAMBDA,J,ISIGNF,DENOM,ADJUST,K FORT2953
795 FORMAT('+',T25,'ENERGY =',F15.8,T52,'LAMBDA =',F8.5,T78,'ATOM =', FORT2954
1 I3,T92,'SIGN =',I3,/,T25,'DENOM =',D16.8,T52,'ADJUST =',D14.7, FORT2955
2 T78,'NL =',I3) FORT2956
C FORT2957
C PRINT OUT ATOMIC CHARGES, ORBITAL OCCUPATIONS, AND CORRECTED FORT2958
C H(I,I)'S IF PRINTX IS TRUE. FORT2959
C FORT2960
IF(.NOT.PRINTX) GO TO 500 FORT2961
WRITE(6,600) FORT2962
600 FORMAT(////,T16,'ATOM',T30,'NET CHG.-DAMPED',T60,'SUMMED ORBITAL O
1CCUPATIONS-DAMPED'/T65,'S',T75,'P',T85,'D'/) FORT2964
DO 650 I=1,NA FORT2965
KEYI=KEY(I) FORT2966
WRITE(6,601) X(I) FORT2967
601 FORMAT(T60,F10.5) FORT2968
IF(NP(KEYI).EQ.0) GO TO 625 FORT2969
WRITE(6,602) Y(I) FORT2970
602 FORMAT('+',T70,F10.5) FORT2971
IF(ND(KEYI).EQ.0) GO TO 625 FORT2972
WRITE(6,603) Z(I) FORT2973
603 FORMAT('+',T80,F10.5) FORT2974
625 UB=C(I) FORT2975
650 WRITE(6,604) SYMBOL(KEYI),I,UB FORT2976
604 FORMAT('+',T15,A2,I3,T30,F10.5) FORT2977
WRITE(6,809) FORT2978
809 FORMAT(///,T16,'ATOM',T70,'CORRECTED H(I,I)''S',/T40,'S',T50, FORT2979
1 'X',T60,'Y',T70,'Z',T80,'X2-Y2',T90,'Z2',T100,'XY',T110,'XZ', FORT2980
2 T120,'YZ'/) FORT2981
J=1 FORT2982
DO 810 I=1,NATM FORT2983
KEYI=KEY(I) FORT2984
N=J FORT2985
IF(NP(KEYI).NE.0) N=J+3 FORT2986
IF(ND(KEYI).NE.0) N=J+8 FORT2987
815 WRITE(6,816) SYMBOL(KEYI),I,(H(K,K),K=J,N) FORT2988
816 FORMAT(T15,A2,I3,T35,9F10.5) FORT2989
810 J=N+1 FORT2990
IF(METH.GE.3.AND.DABS(QON).GT.0.0001D0) WRITE(6,870) DGSUM FORT2991
870 FORMAT(///,T15,'AVERAGE SHIFT OF MO''S DUE TO NON-ZERO TOTAL CHARGFORT2992
1E =',F12.8,' EV.') FORT2993
IF(METH.GE.3.AND.ICYCLE.EQ.MAXCYC) WRITE(6,808) A1,A2,A3 FORT2994
808 FORMAT(///,T15,'ENERGY CORRECTIONS0',//,T20,'ONE-CENTER',T40, FORT2995
1F16.8,' EV.',//,T20,'TWO-CENTER',T40,F16.8,' EV.',//,T20,'TOTAL', FORT2996
2T40,F16.8,' EV.') FORT2997
WRITE(6,752) FORT2998
752 FORMAT(///) FORT2999
500 ICYCLE=ICYCLE+1 FORT3000
C FORT3001
C CHECK FOR CONVERGENCE ( IE. DENOM LESS THAN DELTAC ). FORT3002
C FORT3003
IF(ICYCLE.GE.MAXCYC) GO TO 433 FORT3004
IF(ICYCLE.LE.2) GO TO 433 FORT3005
IF(DENOM.GE.DELTAC) GO TO 433 FORT3006
ICYCLE=MAXCYC FORT3007
PRTCYC=-PRTCYC FORT3008
433 RETURN FORT3009
END FORT3010
SUBROUTINE OUTPUT(H,U,MAD,C,E,W,IOCC,HDG,NDIM,NTYPE,NC,NHDG) FORT3011
C FORT3012
C SUBROUTINE TO ANALYSE AND PRINT OUT RESULTS. FORT3013
C FORT3014
IMPLICIT REAL*8(A-H,O-Z) FORT3015
DIMENSION H(NDIM,NDIM),U(NDIM,NDIM),MAD(NTYPE,NTYPE),C(NC), FORT3016
1 E(NDIM),W(NDIM),IOCC(NDIM),HDG(NHDG) FORT3017
REAL*8 MAD FORT3018
REAL*4 IOCC FORT3019
PARAMETER (MAXATM=500)
PARAMETER (BB=250)
PARAMETER (MXUSER=230)
PARAMETER (MXUSR2=231)
COMMON/TITLE/AB(10) FORT3020
COMMON/CNTRL/CON,PEEP,COULH,NH,NA,NATM,KA,NELEC,METH,IPRINT, FORT3021
1 IPUNCH,L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT3022
LOGICAL*1 L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT3023
COMMON/OUT/PRT(20),PUN(20),IOVPOP(24),IENRGY(24) FORT3024
LOGICAL*1 PRT,PUN FORT3025
INTEGER*2 IOVPOP,IENRGY FORT3026
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 FORT3030
COMMON/ITPARM/DAMP1,DAMP2,DAMP3,LAMPRI,DELTAC,SENSE,MAXCYC,PRTCYC,
1 ICYCLE,NCON,PARTIT,PRINTX,ITABLE(20) FORT3032
REAL*8 LAMPRI FORT3033
INTEGER*4 PRTCYC FORT3034
LOGICAL*1 PARTIT,PRINTX,ITABLE FORT3035
INTEGER*2 SYMB,HYDROG FORT3036
DATA HYDROG/' H'/ FORT3037
DO 5 I=1,NDIM FORT3038
NMIN=I FORT3039
IF(IOCC(I).GT.0.0001) GO TO 7 FORT3040
5 CONTINUE FORT3041
7 IF(ITERAT) GO TO 38 FORT3042
C FORT3043
C CALCULATE AND PRINT OUT SUM OF ONE-ELECTRON ENERGIES. FORT3044
C FORT3045
10 SUM=0.0D0 FORT3046
DO 13 I=1,NDIM FORT3047
W(I)=DBLE(IOCC(I)) FORT3048
13 SUM=SUM+E(I)*W(I) FORT3049
IF(.NOT.PRT(9)) WRITE(6,2001) SUM FORT3050
2001 FORMAT(T10,'SUM OF ONE-ELECTRON ENERGIES =',F16.8,' EV.',///) FORT3051
IF(PUN(9)) WRITE(7,2002) SUM FORT3052
2002 FORMAT(F20.8) FORT3053
IF(ONEMAT) GO TO 9999 FORT3054
C FORT3055
C PRINT OUT WAVE FUNCTIONS. FORT3056
C FORT3057
IF(PRT(10)) GO TO 1003 FORT3058
WRITE(6,1002) FORT3059
1002 FORMAT('WAVE FUNCTIONS'/'MO''S IN COLUMNS, AO''S IN ROWS') FORT3060
CALL PEGLEG(H,NDIM,NDIM) FORT3061
C
C ** THIS CALL IS TO A SIMILAR ROUTINE TO WRITE THE MO COEFFICIENTS
C TO DISK FILE 13. THE FORMAT IS DIFFERENT.
C
CALL PEGLEG2(H,NDIM,NDIM)
C
C **
C
1003 IF(PUN(10)) WRITE(7,2003) H FORT3062
2003 FORMAT(8F9.6) FORT3063
C FORT3064
C CALCULATE AND PRINT OUT DENSITY MATRIX. FORT3065
C FORT3066
500 IF(PRT(11).AND..NOT.PUN(11)) GO TO 38 FORT3067
DO 300 I=1,NDIM FORT3068
DO 300 J=1,I FORT3069
U(I,J)=0.0D0 FORT3070
DO 310 K=1,NDIM FORT3071
310 U(I,J)=U(I,J)+H(I,K)*H(J,K)*W(K) FORT3072
300 U(J,I)=U(I,J) FORT3073
IF(PRT(11)) GO TO 360 FORT3074
WRITE(6,350) FORT3075
350 FORMAT('DENSITY MATRIX') FORT3076
CALL PEGLEG(U,NDIM,NDIM) FORT3077
360 IF(PUN(11)) WRITE(7,2003) U FORT3078
C FORT3079
C CALCULATE ATOMIC ORBITAL OCCUPATIONS AND STORE IN E(I). FORT3080
C CALCULATE OVERLAP POPULATION MATRIX. FORT3081
C FORT3082
38 IJ=1 FORT3083
DO 60 I=1,NDIM FORT3084
E(I)=0.0D0 FORT3085
DO 60 J=1,I FORT3086
UB=0.0D0 FORT3087
DO 41 K=NMIN,NDIM FORT3088
41 UB=UB+H(I,K)*H(J,K)*DBLE(IOCC(K)) FORT3089
UB=UB*0.5D0 FORT3090
IF(I.EQ.J) GO TO 50 FORT3091
UB=(UB+UB)*C(IJ) FORT3092
IJ=IJ+1 FORT3093
50 E(I)=E(I)+UB FORT3094
E(J)=E(J)+UB FORT3095
IF(ITERAT) GO TO 60 FORT3096
UB=UB+UB FORT3097
U(I,J)=UB FORT3098
U(J,I)=UB FORT3099
60 CONTINUE FORT3100
C FORT3101
C IF DOING CHARGE ITERATION ( METH=1 ) CALL LITER. FORT3102
C FORT3103
IF(.NOT.ITERAT) GO TO 80 FORT3104
K=ICYCLE FORT3105
PRINTX=((ICYCLE/PRTCYC)*PRTCYC.EQ.ICYCLE) FORT3106
IF(ICYCLE.EQ.MAXCYC) PRINTX=.TRUE. FORT3107
CALL LITER(NH,NA,E,W,J) FORT3108
IF(ICYCLE.EQ.15000) PRINTX=.TRUE. FORT3109
IF(K.EQ.1) WRITE(6,355) FORT3110
355 FORMAT(' ATOMIC CHARGES',//) FORT3111
GO TO (81,82,83,84),J FORT3112
81 WRITE(6,375) K,(X(I),I=1,NATM) FORT3113
375 FORMAT(/,T3,'CYCLE NO.',I3,(T20,10F10.5)) FORT3114
GO TO 9000 FORT3115
82 WRITE(6,375) K,(Y(I),I=1,NATM) FORT3116
GO TO 9000 FORT3117
83 WRITE(6,375) K,(Z(I),I=1,NATM) FORT3118
GO TO 9000 FORT3119
84 WRITE(6,375) K,(W(I),I=1,NATM) FORT3120
9000 RETURN FORT3121
C FORT3122
C PRINT OUT OVERLAP POPULATION MATRIX. FORT3123
C FORT3124
80 IF(PRT(12)) GO TO 1009 FORT3125
WRITE(6,1006) NELEC FORT3126
1006 FORMAT('OVERLAP POPULATION MATRIX FOR',I4,' ELECTRONS') FORT3127
CALL PEGLEG(U,NDIM,NDIM) FORT3128
1009 IF(PUN(12)) WRITE(7,2003) U FORT3129
C FORT3130
C CALL REDUCE TO CALCULATE REDUCED OVERLAP MATRIX. PRINT OUT FORT3131
C REDUCED OVERLAP MATRIX. FORT3132
C FORT3133
IF(PRT(13).AND..NOT.PUN(13)) GO TO 2005 FORT3134
CALL REDUCE(U,NDIM,NA,NH) FORT3135
DO 100 I=2,NATM FORT3136
K=I-1 FORT3137
DO 100 J=1,K FORT3138
100 U(J,I)=U(I,J) FORT3139
IF(PRT(13)) GO TO 1015 FORT3140
WRITE(6,1007) FORT3141
1007 FORMAT('0REDUCED OVERLAP POPULATION MATRIX, ATOM BY ATOM') FORT3142
CALL PEGLEG(U,NATM,NDIM) FORT3143
1015 IF(PUN(13)) WRITE(7,2003) ((U(I,J),I=1,NATM),J=1,NATM) FORT3144
C FORT3145
C IF L3 IS TRUE, CALCULATE AND PRINT OUT OVERLAP POPULATION ANALYSIS, FORT3146
C ORBITAL BY ORBITAL, FOR EACH MOLECULAR ORBITAL SPECIFIED. FORT3147
C FORT3148
2005 IF(.NOT.L3) GO TO 21 FORT3149
DO 600 N=1,23,2 FORT3150
IF(IOVPOP(N).EQ.0) GO TO 25 FORT3151
KMIN=IOVPOP(N) FORT3152
KMAX=IOVPOP(N+1) FORT3153
DO 600 K=KMIN,KMAX FORT3154
WRITE(6,2004) K,W(K) FORT3155
2004 FORMAT(///'0OVERLAP POPULATION MATRIX, ORBITAL BY ORBITAL, FOR MOLFORT3156
1ECULAR ORBITAL',I4,5X,'OCCUPATION IS',F7.4,' ELECTRONS') FORT3157
IF(IOCC(K).LT.0.0001) GO TO 600 FORT3158
SUM=0.5D0*W(K) FORT3159
IJ=1 FORT3160
DO 460 I=1,NDIM FORT3161
DO 460 J=1,I FORT3162
UB=H(I,K)*H(J,K) FORT3163
IF(I.EQ.J) GO TO 450 FORT3164
UB=(UB+UB)*C(IJ) FORT3165
IJ=IJ+1 FORT3166
450 UB=(UB+UB)*SUM FORT3167
U(J,I)=UB FORT3168
U(I,J)=UB FORT3169
460 CONTINUE FORT3170
CALL PEGLEG(U,NDIM,NDIM) FORT3171
600 CONTINUE FORT3172
25 WRITE(6,7003) FORT3173
7003 FORMAT(///) FORT3174
C FORT3175
C CALL FULCHM TO CALCULATE COMPLETE CHARGE MATRIX. PRINT OUT FORT3176
C COMPLETE CHARGE MATRIX. FORT3177
C FORT3178
21 L1=PRT(14).AND..NOT.PUN(14) FORT3179
L2=PRT(15).AND..NOT.PUN(15) FORT3180
IF(L1.AND.L2) GO TO 1020 FORT3181
CALL FULCHM(W,U,C,H,NDIM) FORT3182
IF(PRT(14)) GO TO 2021 FORT3183
WRITE(6,1008) FORT3184
1008 FORMAT('0COMPLETE CHARGE MATRIX FOR EACH MO, NORMALIZED TO TWO ELEFORT3185
1CTRONS REGARDLESS OF OCCUPATION') FORT3186
CALL PEGLEG(U,NDIM,NDIM) FORT3187
2021 IF(PUN(14)) WRITE(7,2003) U FORT3188
C FORT3189
C CALL REDCHM TO CALCULATE REDUCED CHARGE MATRIX. PRINT OUT FORT3190
C REDUCED CHARGE MATRIX. FORT3191
C FORT3192
IF(L2) GO TO 1020 FORT3193
CALL REDCHM(U,NDIM) FORT3194
IF(PRT(15)) GO TO 1022 FORT3195
WRITE(6,1019) FORT3196
1019 FORMAT('0REDUCED CHARGE MATRIX, MO''S IN COLUMNS, ATOMS IN ROWS') FORT3197
CALL OUTMAT(U,NDIM,NATM,NDIM) FORT3198
1022 IF(PUN(15)) WRITE(7,2003) ((U(I,J),I=1,NATM),J=1,NDIM) FORT3199
C FORT3200
C PRINT OUT ATOMIC CHARGES AND ORBITAL OCCUPATIONS. FORT3201
C FORT3202
1020 IF(PRT(16)) GO TO 40 FORT3203
WRITE(6,1010) FORT3204
1010 FORMAT(/'0ATOM',T12,'NET CHG.',T35,'ATOMIC ORBITAL OCCUPATION FOR FORT3205
*GIVEN MO OCCUPATION'/T35,'S',T45,'X',T55,'Y',T65,'Z',T75,'X2-Y2', FORT3206
*T85,'Z2',T95,'XY',T105,'XZ',T115,'YZ'/) FORT3207
SYMB=HYDROG FORT3208
J=1 FORT3209
DO 140 I=1,NATM FORT3210
IF(I.GT.NH) GO TO 120 FORT3211
UB=1.0-E(I) FORT3212
N=I FORT3213
GO TO 130 FORT3214
120 KITE=I-NH FORT3215
KEYI=KEY(KITE) FORT3216
SYMB=SYMBOL(KEYI) FORT3217
UB=VELEC(KEYI)-E(J) FORT3218
N=J FORT3219
IF(NP(KEYI).EQ.0) GO TO 130 FORT3220
UB=UB-E(J+1)-E(J+2)-E(J+3) FORT3221
N=J+3 FORT3222
IF(ND(KEYI).EQ.0) GO TO 130 FORT3223
UB=UB-E(J+4)-E(J+5)-E(J+6)-E(J+7)-E(J+8) FORT3224
N=J+8 FORT3225
130 WRITE(6,1011) SYMB,I,UB,(E(K),K=J,N) FORT3226
1011 FORMAT(1X,A2,I3,T10,F10.5,T30,9F10.5) FORT3227
140 J=N+1 FORT3228
C FORT3229
C IF CALCULATING ENERGY MATRIX, PUT DIAGONAL ELEMENTS OF HUCKEL FORT3230
C MATRIX IN W(I). FORT3231
C FORT3232
40 PRT(1)=PRT(17).AND..NOT.PUN(17) FORT3233
PRT(2)=PRT(18).AND..NOT.PUN(18) FORT3234
PRT(3)=PRT(19).AND..NOT.PUN(19) FORT3235
PRT(4)=PRT(20).AND..NOT.PUN(20) FORT3236
L1=PRT(1).AND.PRT(2) FORT3237
L2=PRT(3).AND.PRT(4) FORT3238
IF(L1.AND.L2.AND..NOT.L4) GO TO 9999 FORT3239
IF(NH.EQ.0) GO TO 23 FORT3240
DO 22 I=1,NH FORT3241
22 W(I)=X(I) FORT3242
23 J=NH+1 FORT3243
K=NH+1 FORT3244
DO 24 I=1,NA FORT3245
KEYI=KEY(I) FORT3246
W(J)=X(K) FORT3247
J=J+1 FORT3248
IF(NP(KEYI).EQ.0) GO TO 24 FORT3249
UB=Y(K) FORT3250
W(J)=UB FORT3251
W(J+1)=UB FORT3252
W(J+2)=UB FORT3253
J=J+3 FORT3254
IF(ND(KEYI).EQ.0) GO TO 24 FORT3255
UB=Z(K) FORT3256
W(J)=UB FORT3257
W(J+1)=UB FORT3258
W(J+2)=UB FORT3259
W(J+3)=UB FORT3260
W(J+4)=UB FORT3261
J=J+5 FORT3262
24 K=K+1 FORT3263
C FORT3264
C IF DOING CHARGE ITERATION WITH MADELUNG CORRECTION ON MOLECULE FORT3265
C WITH NON-ZERO CHARGE, PUT MADELUNG TERMS IN E(I). FORT3266
C FORT3267
QON= DFLOAT(KA)/DFLOAT(NELEC) FORT3268
IF(METH.LT.3.OR.DABS(QON).LT.0.0001D0) GO TO 180 FORT3269
IF(NH.EQ.0) GO TO 181 FORT3270
DO 182 I=1,NH FORT3271
182 E(I)=MAD(I,I) FORT3272
181 J=NH+1 FORT3273
K=NH+1 FORT3274
DO 183 I=1,NA FORT3275
KEYI=KEY(I) FORT3276
E(J)=MAD(K,K) FORT3277
J=J+1 FORT3278
K=K+1 FORT3279
IF(NP(KEYI).EQ.0) GO TO 183 FORT3280
UB=MAD(K,K) FORT3281
E(J)=UB FORT3282
E(J+1)=UB FORT3283
E(J+2)=UB FORT3284
J=J+3 FORT3285
K=K+1 FORT3286
IF(ND(KEYI).EQ.0) GO TO 183 FORT3287
UB=MAD(K,K) FORT3288
E(J)=UB FORT3289
E(J+1)=UB FORT3290
E(J+2)=UB FORT3291
E(J+3)=UB FORT3292
E(J+4)=UB FORT3293
J=J+5 FORT3294
K=K+1 FORT3295
183 CONTINUE FORT3296
C FORT3297
C CALCULATE AND PRINT OUT ENERGY MATRIX. FORT3298
C FORT3299
180 ONEMAT=.TRUE. FORT3300
IF(L1) GO TO 7000 FORT3301
170 SUM=1.0D0 FORT3302
IF(ONEMAT) SUM=0.0D0 FORT3303
IJ=1 FORT3304
CNST=2.0D0*CON FORT3305
CN2=CNST FORT3306
DO 28 I=1,NDIM FORT3307
U(I,I)=0.0D0 FORT3308
DO 28 J=1,I FORT3309
UB=0.0D0 FORT3310
DO 26 K=NMIN,NDIM FORT3311
26 UB=UB+H(I,K)*H(J,K)*DBLE(IOCC(K)) FORT3312
IF(I.EQ.J) GO TO 28 FORT3313
IF(.NOT.L5) GO TO 35 FORT3314
UC=W(I) FORT3315
ET=W(J) FORT3316
IF(NHDG.EQ.1) GO TO 36 FORT3317
UC=HDG(I) FORT3318
ET=HDG(J) FORT3319
36 UC=(UC-ET)/(UC+ET) FORT3320
UC=UC*UC FORT3321
CNST=CN2+UC+UC*UC*(1.0D0-CN2) FORT3322
35 IF(METH.LT.3.OR.DABS(QON).LT.0.0001D0) GO TO 31 FORT3323
UC=UB*C(IJ)*((CNST-SUM)*(W(I)+W(J))-QON*(1.0D0-CNST)*(E(I)+E(J))) FORT3324
GO TO 32 FORT3325
31 UC=UB*(CNST-SUM)*C(IJ)*(W(I)+W(J)) FORT3326
32 UB=SUM*UB*C(IJ) FORT3327
IJ=IJ+1 FORT3328
U(J,I)=UC FORT3329
U(I,J)=UC FORT3330
U(J,J)=U(J,J)+UB*W(J) FORT3331
28 U(I,I)=U(I,I)+UB*W(I) FORT3332
IF(PRT(17)) GO TO 3020 FORT3333
IF(ONEMAT) WRITE(6,3005) FORT3334
3005 FORMAT(///,'0ENERGY MATRIX') FORT3335
IF(.NOT.ONEMAT) WRITE(6,3006) FORT3336
3006 FORMAT('0ENERGY PARTITIONING') FORT3337
CALL PEGLEG(U,NDIM,NDIM) FORT3338
3020 IF(PUN(17)) WRITE(7,2800) U FORT3339
2800 FORMAT(8F9.5) FORT3340
C FORT3341
C CALL REDUCE TO CALCULATE REDUCED ENERGY MATRIX. PRINT OUT FORT3342
C REDUCED ENERGY MATRIX. FORT3343
C FORT3344
IF(PRT(2)) GO TO 7000 FORT3345
CALL REDUCE(U,NDIM,NA,NH) FORT3346
DO 700 I=2,NATM FORT3347
K=I-1 FORT3348
DO 700 J=1,K FORT3349
700 U(J,I)=U(I,J) FORT3350
IF(PRT(18)) GO TO 3021 FORT3351
IF(ONEMAT) WRITE(6,3007) FORT3352
3007 FORMAT('0REDUCED ENERGY MATRIX, ATOM BY ATOM') FORT3353
IF(.NOT.ONEMAT) WRITE(6,3008) FORT3354
3008 FORMAT('0REDUCED ENERGY PARTITIONING, ATOM BY ATOM') FORT3355
KITE=0 FORT3356
701 LOW=KITE+1 FORT3357
KITE=KITE+13 FORT3358
IF(KITE.GT.NATM) KITE=NATM FORT3359
WRITE(6,702) (I,I=LOW,KITE) FORT3360
702 FORMAT(/5X,13I9,//) FORT3361
DO 703 I=1,NATM FORT3362
703 WRITE(6,704) I,(U(I,J),J=LOW,KITE) FORT3363
704 FORMAT(I5,2X,13F9.4) FORT3364
IF(KITE.LT.NATM) GO TO 701 FORT3365
3021 IF(PUN(18)) WRITE(7,2111) ((U(I,J),I=1,NATM),J=1,NATM) FORT3366
2111 FORMAT(7F10.5) FORT3367
C FORT3368
C IF L4 IS TRUE, CALCULATE AND PRINT OUT ENERGY MATRIX ANALYSIS, FORT3369
C ORBITAL BY ORBITAL, FOR EACH MOLECULAR ORBITAL SPECIFIED. FORT3370
C FORT3371
7000 IF(.NOT.ONEMAT) GO TO 9999 FORT3372
IF(.NOT.L4) GO TO 71 FORT3373
DO 246 N=1,23,2 FORT3374
IF(IENRGY(N).EQ.0) GO TO 72 FORT3375
KMIN=IENRGY(N) FORT3376
KMAX=IENRGY(N+1) FORT3377
DO 246 K=KMIN,KMAX FORT3378
WRITE(6,3004) K,IOCC(K) FORT3379
3004 FORMAT(///'0ENERGY MATRIX, ORBITAL BY ORBITAL, FOR MOLECULAR ORBITFORT3380
1AL',I4,5X,'OCCUPATION IS',F7.4,' ELECTRONS') FORT3381
IF(IOCC(K).LT.0.0001) GO TO 246 FORT3382
IJ=1 FORT3383
CNST=CON FORT3384
EX=DBLE(IOCC(K)) FORT3385
DO 244 I=1,NDIM FORT3386
DO 244 J=1,I FORT3387
UB=H(I,K)*H(J,K) FORT3388
IF(I.EQ.J) GO TO 242 FORT3389
IF(.NOT.L5) GO TO 236 FORT3390
UC=W(I) FORT3391
ET=W(J) FORT3392
IF(NHDG.EQ.1) GO TO 237 FORT3393
UC=HDG(I) FORT3394
ET=HDG(J) FORT3395
237 UC=(UC-ET)/(UC+ET) FORT3396
UC=UC*UC FORT3397
CNST=CON+UC/2.0D0+UC*UC*(0.5D0-CON) FORT3398
236 IF(METH.LT.3.OR.DABS(QON).LT.0.0001D0) GO TO 247 FORT3399
UB=UB*C(IJ)*(CNST*(W(I)+W(J))-QON*(0.5D0-CNST)*(E(I)+E(J))) FORT3400
GO TO 248 FORT3401
247 UB=UB*CNST*C(IJ)*(W(I)+W(J)) FORT3402
248 IJ=IJ+1 FORT3403
UB=2.0D0*UB*EX FORT3404
U(I,J)=UB FORT3405
U(J,I)=UB FORT3406
GO TO 244 FORT3407
242 U(I,I)=UB*W(I)*EX FORT3408
244 CONTINUE FORT3409
CALL PEGLEG(U,NDIM,NDIM) FORT3410
246 CONTINUE FORT3411
72 WRITE(6,7003) FORT3412
C FORT3413
C CALCULATE AND PRINT OUT ENERGY PARTITIONING AND REDUCED FORT3414
C ENERGY PARTITIONING. FORT3415
C FORT3416
71 IF(L2) GO TO 9999 FORT3417
ONEMAT=.FALSE. FORT3418
PRT(17)=PRT(19) FORT3419
PUN(17)=PUN(19) FORT3420
PRT(18)=PRT(20) FORT3421
PUN(18)=PUN(20) FORT3422
PRT(2)=PRT(4) FORT3423
GO TO 170 FORT3424
9999 RETURN FORT3425
END FORT3426
SUBROUTINE LITER(NH,NA,E,W,NCYC) FORT3427
C FORT3428
C SUBROUTINE FOR CALCULATING Q*SENSE WHEN USING CHARGE ITERATION FORT3429
C OPTION ( METH = 1 ). FORT3430
C FORT3431
IMPLICIT REAL*8(A-H,O-Z) FORT3432
DIMENSION E(1),W(1) FORT3433
PARAMETER (MAXATM=500)
PARAMETER (BB=250)
PARAMETER (MXUSER=230)
PARAMETER (MXUSR2=231)
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 FORT3437
COMMON/ITPARM/DAMP1,DAMP2,DAMP3,LAMPRI,DELTAC,SENSE,MAXCYC,PRTCYC,
1 ICYCLE,NCON,PARTIT,PRINTX,ITABLE(20) FORT3439
REAL*8 LAMPRI FORT3440
INTEGER*4 PRTCYC FORT3441
LOGICAL*1 PARTIT,PRINTX,ITABLE FORT3442
NCYC=ICYCLE+1-(ICYCLE/4)*4 FORT3443
ICYCLE=ICYCLE+1 FORT3444
DELTA=0.D0 FORT3445
NATOM=NH+NA FORT3446
INDEX=1 FORT3447
DO 60 I=1,NATOM FORT3448
INCR=1 FORT3449
IF(I.GT.NH) GO TO 100 FORT3450
CHG=1.0D0-E(INDEX) FORT3451
GO TO 150 FORT3452
100 KEYI=KEY(I-NH) FORT3453
CHG=E(INDEX) FORT3454
IF(ND(KEYI).EQ.0) GO TO 120 FORT3455
INCR=9 FORT3456
CHG=CHG+E(INDEX+4)+E(INDEX+5)+E(INDEX+6)+E(INDEX+7)+E(INDEX+8) FORT3457
GO TO 130 FORT3458
120 IF(NP(KEYI).EQ.0) GO TO 140 FORT3459
INCR=4 FORT3460
130 CHG=CHG+E(INDEX+1)+E(INDEX+2)+E(INDEX+3) FORT3461
140 CHG=VELEC(KEYI)-CHG FORT3462
150 INDEX=INDEX+INCR FORT3463
GO TO (10,20,30,40),NCYC FORT3464
10 CHG=0.25D0*(CHG+Y(I)+Z(I)+W(I)) FORT3465
X(I)=CHG FORT3466
Y(I)=CHG FORT3467
Z(I)=CHG FORT3468
W(I)=CHG FORT3469
GO TO 60 FORT3470
20 DELTA=DELTA+DABS(CHG-X(I)) FORT3471
Y(I)=CHG FORT3472
GO TO 60 FORT3473
30 DELTA=DELTA+DABS(CHG-Y(I)) FORT3474
Z(I)=CHG FORT3475
GO TO 60 FORT3476
40 DELTA=DELTA+DABS(CHG-Z(I)) FORT3477
W(I)=CHG FORT3478
60 E(I)=CHG*SENSE FORT3479
IF(DELTAC.LT.DELTA.OR.NCYC.EQ.1) RETURN FORT3480
ICYCLE=15000 FORT3481
RETURN FORT3482
END FORT3483
SUBROUTINE MATRIX(N,NS1,NS2,NS3,NS4,NS5,NS6,NS7,NS8,NS9,NS10, FORT3484
X NS11,NS12,NS13,ND1,ND2,ND3,ND4,ND5,ND6,ND7,ND8,ND9,ND10, FORT3485
X ND11,ND12,ND13) FORT3486
IMPLICIT REAL*8(A-H,O-Z) FORT3487
DIMENSION A(2000000)
NSIZE=2000000
C FORT3490
C SUBROUTINE TO ALLOCATE STORAGE FOR MATRICES. FORT3491
C FORT3492
CALL ERRSET(74,0,0,0,1,0) FORT3493
I1=1 FORT3494
I2=I1+NS1 FORT3495
I3=I2+NS2 FORT3496
I4=I3+NS3 FORT3497
I5=I4+NS4 FORT3498
I6=I5+NS5 FORT3499
I7=I6+NS6 FORT3500
I8=I7+NS7 FORT3501
I9=I8+NS8 FORT3502
I10=I9+NS9 FORT3503
I11=I10+NS10 FORT3504
I12=I11+NS11 FORT3505
I13=I12+NS12 FORT3506
IF=I13+NS13 FORT3507
IF(IF.GT.NSIZE) GO TO 200 FORT3508
CALL MOVLAP (A(I1),A(I2),A(I3),A(I4),A(I5),A(I6),A(I7), FORT3509
X A(I8),A(I9),A(I10),A(I11),A(I12),A(I13),ND1,ND2,ND3,ND4, FORT3510
X ND5,ND6) FORT3511
CALL HUCKEL (A(I1),A(I2),A(I3),A(I4),A(I5),A(I6),A(I7), FORT3512
X A(I8),A(I9),A(I10),A(I11),A(I12),A(I13),ND1,ND2,ND3,ND4, FORT3513
X ND5,ND6) FORT3514
RETURN FORT3515
200 WRITE(6,1001) IF,NSIZE FORT3516
1001 FORMAT('0*** INSUFFICIENT SPACE FOR MATRICES'/'0PROGRAM REQUESTS' FORT3517
1,I7,' DOUBLE WORDS BUT ONLY',I7,' ARE AVAILABLE'/'0RECOMPILE WITH FORT3518
2LARGER DIMENSION FOR A AND INCREASED VALUE OF NSIZE') FORT3519
RETURN FORT3520
END FORT3521
SUBROUTINE ABFNS(A,B) FORT3522
IMPLICIT REAL*8(A-H,O-Z) FORT3523
PARAMETER (MXUSER=230)
PARAMETER (MXUSR2=231)
DIMENSION A(MXUSER),B(MXUSER)
COMMON/LOCLAP/SK1,SK2,RR,L1,L2,M,N1,N2,MAXCAL FORT3525
C FORT3526
C SUBROUTINE FOR CALCULATING A AND B FUNCTIONS FOR USE IN LOVLAP. FORT3527
C FORT3528
J=MAXCAL+1 FORT3529
RHO1=0.5D0*(SK1+SK2)*RR FORT3530
RHO2=0.5D0*(SK1-SK2)*RR FORT3531
IF(DABS(RHO1).GT.165.D0) GO TO 100 FORT3532
IF(DABS(RHO2).GT.165.D0) GO TO 100 FORT3533
C=DEXP(-RHO1) FORT3534
A(1)=C/RHO1 FORT3535
DO 15 I=2,J FORT3536
15 A(I)=(DFLOAT(I-1)*A(I-1)+C)/RHO1 FORT3537
IX=J FORT3538
IR=DABS(2.*RHO2) FORT3539
IS=MIN0(IR+1,19) FORT3540
IF(RHO2) 25,35,25 FORT3541
25 D=DEXP(RHO2) FORT3542
H=1.D0/D FORT3543
C FORT3544
C USE THE DSINH ROUTINE INSTEAD OF SUMMING THE INFINITE SERIES. FORT3545
C FORT3546
R=2.D0*DSINH(RHO2) FORT3547
C FORT3548
C AS MANY SUCCESSIVE B-FUNCTIONS ARE GENERATED FROM B-0 BY THE FORT3549
C RECURRENCE FORMULAE AS ACCURACY WILL PERMIT. FORT3550
C FORT3551
B(1)=R/RHO2 FORT3552
DO 51 I=2,IX,IS FORT3553
IF(IR.EQ.0) GO TO 40 FORT3554
IL=IS-1 FORT3555
C FORT3556
C MODIFICATION TO AVOID EXCEEDING STORAGE LIMITS. FORT3557
C D. WALLACE 04/14/71 FORT3558
C FORT3559
DO 31 K=I,IX FORT3560
IF((-1)**K) 29,29,30 FORT3561
29 B(K)=(R+DFLOAT(K-1)*B(K-1))/RHO2 FORT3562
GO TO 31 FORT3563
30 B(K)=-(D+H-DFLOAT(K-1)*B(K-1))/RHO2 FORT3564
31 CONTINUE FORT3565
40 IN=I+IS-1 FORT3566
C FORT3567
C AFTER THE RECURRENCE FORMULAE HAVE BEEN APPLIED AN APPROPRIATE FORT3568
C NUMBER OF TIMES THE NEXT B-FUNCTION IS OBTAINED BY SUMMATION FORT3569
C OF THE INFINITE SERIES. FORT3570
C FORT3571
IF(IN-IX) 39,39,38 FORT3572
39 IF((-1)**IN) 44,44,42 FORT3573
42 TR=RHO2 FORT3574
105 B(IN)=-2.*TR/DFLOAT(IN+1) FORT3575
DO 43 J=1,500 FORT3576
TR=TR*RHO2**2/DFLOAT((2*J)*(2*J+1)) FORT3577
IF(DABS(TR/B(IN))-1.0D-7 ) 51,51,43 FORT3578
43 B(IN)=B(IN)-2.*TR/DFLOAT(IN+1+2*J) FORT3579
GO TO 51
44 TR=1. FORT3580
107 B(IN)=2.*TR/DFLOAT(IN) FORT3581
DO 46 J=1,500 FORT3582
TR=TR*RHO2**2/DFLOAT((2*J)*(2*J-1)) FORT3583
IF(DABS(TR/B(IN))-1.0D-7 ) 51,51,46 FORT3584
46 B(IN)=B(IN)+2.*TR/DFLOAT(IN+2*J) FORT3585
51 CONTINUE FORT3586
C FORT3587
C IF THE ARGUMENT IS ZERO A SEPARATE FORMULA MUST BE USED. FORT3588
C FORT3589
GO TO 38 FORT3590
35 DO 36 I=1,IX,2 FORT3591
B(I)=2.D0/DFLOAT(I) FORT3592
36 B(I+1)=0.D0 FORT3593
38 RETURN FORT3594
100 DO 101 I=1,MXUSER
A(I)=0.D0 FORT3596
101 B(I)=0.D0 FORT3597
GO TO 38 FORT3598
END FORT3599
SUBROUTINE LOVLAP (STRAD,A,B) FORT3600
IMPLICIT REAL*8(A-H,O-Z) FORT3601
PARAMETER (MXUSER=230)
PARAMETER (MXUSR2=231)
DIMENSION A(MXUSER),B(MXUSER)
DIMENSION FACT(25) FORT3603
COMMON/LOCLAP/SK1,SK2,R,L1,L2,M1,N1,N2,MAX FORT3604
DIMENSION BINCOE(7,7) FORT3605
DATA BINCOE/7*1.D0, 0.D0,1.D0,2.D0,3.D0,4.D0,5.D0,6.D0, FORT3606
1 2*0.D0,1.D0,3.D0,6.D0,10.D0,15.D0, 3*0.D0,1.D0,4.D0,10.D0, FORT3607
2 20.D0, 4*0.D0,1.D0,5.D0,15.D0, 5*0.D0,1.D0,6.D0, 6*0.D0, FORT3608
3 1.D0/ FORT3609
LOGICAL*1 JGO FORT3610
DATA JGO/.FALSE./ FORT3611
C FORT3612
C SUBROUTINE TO CALCULATE OVERLAP INTEGRALS IN A LOCAL FORT3613
C COORDINATE SYSTEM. FORT3614
C FORT3615
C INTEGRALS ARE CALCULATED BY TRANSFORMATION TO ELLIPSOIDAL FORT3616
C COORDINATES AND THEREBY EXPRESSED IN TERMS OF C-FUNCTIONS. FORT3617
C SEE J.C.P.,24,201. ORIGINALLY WRITTEN BY R.M.STEVENS. FORT3618
C FORT3619
C FORT3620
C GENERATE FACTORIALS ONLY ONCE. FORT3621
C FORT3622
IF(JGO) GO TO 10 FORT3623
JGO=.TRUE. FORT3624
FACT(1)=1.D0 FORT3625
DO 5 I=2,25 FORT3626
5 FACT(I)=FACT(I-1)*DFLOAT(I-1) FORT3627
10 CONTINUE FORT3628
M2=M1 FORT3629
STRAD=0.D0 FORT3630
RHOA=R*SK1 FORT3631
RHOB=R*SK2 FORT3632
TERMA=0.5D0**(L1+L2+1) * DSQRT(DFLOAT((L1+L1+1)*(L2+L2+1))* FORT3633
1 FACT(L1-M1+1)*FACT(L2-M1+1)/(FACT(N1+N1+1)*FACT(N2+N2+1)* FORT3634
2 FACT(L1+M1+1)*FACT(L2+M1+1))*RHOA**(N1+N1+1)*RHOB**(N2+N2+1)) FORT3635
JEND=1+((L1-M1)/2) FORT3636
KEND=1+((L2-M2)/2) FORT3637
IEB=M1+1 FORT3638
DO 50 J=1,JEND FORT3639
JU=J-1 FORT3640
IAB=N1-L1+JU+JU+1 FORT3641
ICB=L1-M1-JU-JU+1 FORT3642
CON1=FACT(L1+L1-JU-JU+1)/(FACT(L1-M1-JU-JU+1)*FACT(JU+1)* FORT3643
1 FACT(L1-JU+1)) FORT3644
DO 50 K=1,KEND FORT3645
KU=K-1 FORT3646
CON12=CON1*FACT(L2+L2-KU-KU+1)/(FACT(L2-M2-KU-KU+1)*FACT(KU+1)* FORT3647
1 FACT(L2-KU+1)) FORT3648
IEV=JU+KU+L2 FORT3649
IF(2*(IEV/2).NE.IEV) CON12=-CON12 FORT3650
IBB=N2-L2+KU+KU+1 FORT3651
IDB=L2-M2-KU-KU+1 FORT3652
VALUE=0.D0 FORT3653
DO 90 I6=1,IEB FORT3654
DO 90 I5=1,IEB FORT3655
VALUE1=BINCOE(IEB,I6)*BINCOE(IEB,I5) FORT3656
IEV=I5+I6 FORT3657
IF(2*(IEV/2).NE.IEV) VALUE1=-VALUE1 FORT3658
DO 90 I4=1,IDB FORT3659
VALUE1=-VALUE1 FORT3660
VALUE2=BINCOE(IDB,I4)*VALUE1 FORT3661
DO 90 I3=1,ICB FORT3662
VALUE3=BINCOE(ICB,I3)*VALUE2 FORT3663
DO 90 I2=1,IBB FORT3664
VALUE3=-VALUE3 FORT3665
VALUE4=BINCOE(IBB,I2)*VALUE3 FORT3666
DO 90 I1=1,IAB FORT3667
TERM=VALUE4*BINCOE(IAB,I1) FORT3668
IR=I1+I2+IEB+IEB-I6-I6-I3+IDB-I4+ICB-1 FORT3669
IP=IAB-I1+IBB-I2+IEB+IEB-I5-I5+ICB-I3+IDB-I4+1 FORT3670
90 VALUE=VALUE+A(IP)*B(IR)*TERM FORT3671
50 STRAD=STRAD+VALUE*CON12 FORT3672
STRAD=STRAD*TERMA FORT3673
RETURN FORT3674
END FORT3675
SUBROUTINE GRMSCH(U,C,NDIM) FORT3676
IMPLICIT REAL*8(A-H,O-Z) FORT3677
DIMENSION U(NDIM,NDIM),C(NDIM) FORT3678
C FORT3679
C SUBROUTINE TO CARRY OUT A GRAM SCHMIDT ORTHOGONALIZATION ON A FORT3680
C SET OF VECTORS X(I). THEY ARE CONVERTED INTO A SET OF ORTHONORMAL FORT3681
C VECTORS Y(J). THE UNITARY TRANSFORMATION IS DEVELOPED. FORT3682
C FORT3683
C U CONTAINS THE OVERLAP MATRIX IN THE UPPER RIGHT TRIANGULAR FORT3684
C PART. THE LOWER TRIANGULAR PART IS NOT USED. THE UNITARY FORT3685
C TRANSFORMATION IS RETURNED IN THE UPPER RIGHT TRIANGLE, FORT3686
C INCLUDING THE DIAGONAL. FORT3687
C C IS A WORK VECTOR OF LENGTH NDIM. FORT3688
C FORT3689
C NDIM IS THE DIMENSION OF C AND U (IE. C(NDIM), U(NDIM,NDIM)). FORT3690
C FORT3691
C ACTUALLY THE FIRST COLUMN OF U MAY BE USED AS THE WORK AREA FORT3692
C PROVIDED THAT THE ELEMENT U(1,1) IS SET EQUAL TO 1.D0 AFTER FORT3693
C RETURN IS MADE TO THE CALLING PROGRAM. FORT3694
C FORT3695
U(1,1)=1.D0 FORT3696
DO 100 I=2,NDIM FORT3697
I1=I-1 FORT3698
XNORM=1.D0 FORT3699
DO 40 J=1,I1 FORT3700
SUM=0.D0 FORT3701
DO 30 K=1,J FORT3702
SUM=SUM+U(K,J)*U(K,I) FORT3703
30 CONTINUE FORT3704
C(J)=SUM FORT3705
XNORM=XNORM-SUM**2 FORT3706
40 CONTINUE FORT3707
XNORM=DSQRT(1.D0/XNORM) FORT3708
DO 70 J=1,I1 FORT3709
SUM=0.D0 FORT3710
DO 60 K=J,I1 FORT3711
SUM=SUM+C(K)*U(J,K) FORT3712
60 CONTINUE FORT3713
U(J,I)=-SUM*XNORM FORT3714
70 CONTINUE FORT3715
U(I,I)=XNORM FORT3716
100 CONTINUE FORT3717
RETURN FORT3718
END FORT3719
SUBROUTINE TRNFRM(S,H,C,COUL0,NDIM,SP,IEXIT) FORT3720
IMPLICIT REAL*8(A-H,O-Z) FORT3721
DIMENSION S(NDIM,NDIM),H(NDIM,NDIM),C(NDIM),COUL0(NDIM),SP(NDIM) FORT3722
LOGICAL*1 ONEMAT FORT3723
C FORT3724
C SUBROUTINE FOR CARRYING OUT A CHANGE IN BASIS SET ON A MATRIX H FORT3725
C BY MEANS OF A MATRIX U. FORT3726
C FORT3727
C C=U'*H*U FORT3728
C FORT3729
C H IS ASSUMED REAL SYMMETRIC AND ONLY THE LOWER LEFT TRIANGLE FORT3730
C IS REFERENCED. U IS ASSUMED TO BE ENTIRELY CONTAINED IN ITS FORT3731
C UPPER RIGHT TRIANGLE AND ONLY THIS PART IS REFERENCED. THE FORT3732
C RESULTING REAL SYMMETRIC MATRIX, C, IS STORED IN PACKED FORM. FORT3733
C FORT3734
C INDEX1 OF LOOP1 POINTS TO U(1,J), J=1,NDIM FORT3735
C INDEX2 OF LOOP2 POINTS TO H(I,1), I=1,J FORT3736
C INDEX3 OF LOOP3 POINTS TO U(K,J), K=1,I FORT3737
C INDEX4 OF LOOP4 POINTS TO U(K,J), K=I+1,J FORT3738
C INDEX5 OF LOOP5 POINTS TO U(1,I), I=1,J FORT3739
C INDEX6 OF LOOP6 POINTS TO U(K,I), K=1,I FORT3740
C FORT3741
FORT3742
ONEMAT=IEXIT.EQ.2.OR.IEXIT.EQ.3 FORT3743
ISUB=1 FORT3744
DO 100 J=1,NDIM FORT3745
DO 40 I=1,J FORT3746
SUM=0.D0 FORT3747
ILIM=I FORT3748
IF(.NOT.ONEMAT) GO TO 10 FORT3749
IF(I.LT.J) GO TO 10 FORT3750
ILIM=I-1 FORT3751
C FORT3752
C IN THE ONE MATRIX CASE, THE S AND H MATRICES OVERLAP ALONG THE FORT3753
C DIAGONAL. THE DIAGONAL OF S IS THEN PUT INTO SP. FORT3754
C FORT3755
SUM=SP(I)*H(I,I) FORT3756
IF(I.EQ.1) GO TO 25 FORT3757
10 DO 20 K=1,ILIM FORT3758
20 SUM=SUM+S(K,J)*H(I,K) FORT3759
25 IF(I.EQ.J) GO TO 40 FORT3760
I1=I+1 FORT3761
ILIM=J FORT3762
IF(.NOT.ONEMAT) GO TO 27 FORT3763
ILIM=J-1 FORT3764
SUM=SUM+SP(J)*H(J,I) FORT3765
IF(ILIM.LT.I1) GO TO 40 FORT3766
27 DO 30 K=I1,ILIM FORT3767
30 SUM=SUM+S(K,J)*H(K,I) FORT3768
40 COUL0(I)=SUM FORT3769
DO 60 I=1,J FORT3770
SUM=0.D0 FORT3771
ILIM=I FORT3772
IF(.NOT.ONEMAT) GO TO 45 FORT3773
ILIM=I-1 FORT3774
SUM=COUL0(I)*SP(I) FORT3775
IF(ILIM.EQ.0) GO TO 55 FORT3776
45 DO 50 K=1,ILIM FORT3777
50 SUM=SUM+COUL0(K)*S(K,I) FORT3778
55 C(ISUB)=SUM FORT3779
60 ISUB=ISUB+1 FORT3780
100 CONTINUE FORT3781
RETURN FORT3782
END FORT3783
DOUBLE PRECISION FUNCTION DSUM(B,A,IP1,LIMIT) FORT3784
IMPLICIT REAL*8(A-H,O-Z) FORT3785
DIMENSION B(1),A(1) FORT3786
C FORT3787
C FUNCTION FOR USE IN GIVENS DIAGNOLIZATION PACKAGE. FORT3788
C FORT3789
JJ=1 FORT3790
DSUM=0.D0 FORT3791
DO 180 II=IP1,LIMIT FORT3792
DSUM=DSUM+B(II+1)*A(JJ) FORT3793
180 JJ=JJ+II FORT3794
RETURN FORT3795
END FORT3796
SUBROUTINE ROTATE(V,C,S,NJX,JTOP) FORT3797
IMPLICIT REAL*8(A-H,O-Z) FORT3798
DIMENSION V(1) FORT3799
C FORT3800
C FUNCTION FOR USE IN GIVENS DIAGNOLIZATION PACKAGE. FORT3801
C FORT3802
JLIM=JTOP+1 FORT3803
DO 10 J=1,JLIM,NJX FORT3804
TA=V(J) FORT3805
TB=V(J+1) FORT3806
V(J)=TA*C+TB*S FORT3807
10 V(J+1)=TB*C-TA*S FORT3808
RETURN FORT3809
END FORT3810
DOUBLE PRECISION FUNCTION DOT(A,B) FORT3811
IMPLICIT REAL*8(A-H,O-Z) FORT3812
DIMENSION A(1),B(1) FORT3813
COMMON /VECTOR/ FACTOR,LIMIT FORT3814
C FORT3815
C FUNCTION FOR USE IN GIVENS DIAGNOLIZATION PACKAGE. FORT3816
C FORT3817
ITOP=LIMIT+1 FORT3818
DOT=0.D0 FORT3819
DO 10 I=1,ITOP FORT3820
10 DOT=DOT+A(I)*B(I) FORT3821
RETURN FORT3822
END FORT3823
SUBROUTINE VECSUM(A,B) FORT3824
IMPLICIT REAL*8(A-H,O-Z) FORT3825
DIMENSION A(1),B(1) FORT3826
COMMON /VECTOR/ FACTOR,LIMIT FORT3827
C FORT3828
C FUNCTION FOR USE IN GIVENS DIAGNOLIZATION PACKAGE. FORT3829
C FORT3830
ITOP=LIMIT+1 FORT3831
DO 10 I=1,ITOP FORT3832
10 A(I)=A(I)+FACTOR*B(I) FORT3833
RETURN FORT3834
END FORT3835
SUBROUTINE REDUCE(U,NDIM) FORT3836
IMPLICIT REAL*8(A-H,O-Z) FORT3837
DIMENSION U(NDIM,NDIM) FORT3838
PARAMETER (MAXATM=500)
PARAMETER (BB=250)
PARAMETER (MXUSER=230)
PARAMETER (MXUSR2=231)
COMMON/CNTRL/CON,PEEP,COULH,NH,NA,NATM,KA,NELEC,METH,IPRINT, FORT3839
1 IPUNCH,L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT3840
LOGICAL*1 L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT3841
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),
2 COULS(BB),COULP(BB),COULD(BB),X(MAXATM),Y(MAXATM),Z(MAXATM)
INTEGER*2 SYMBOL,KEY,VELEC FORT3845
C FORT3846
C SUBROUTINE TO CALCULATE REDUCED MATRIX. FORT3847
C FORT3848
CALL REDCHM(U,NDIM) FORT3849
ISUB=NH+1 FORT3850
IAO=ISUB FORT3851
DO 100 I=1,NA FORT3852
KEYI=KEY(I) FORT3853
IF(ND(KEYI)) 10,20,10 FORT3854
10 IATOP=IAO+8 FORT3855
GO TO 50 FORT3856
20 IF(NP(KEYI)) 30,40,30 FORT3857
30 IATOP=IAO+3 FORT3858
GO TO 50 FORT3859
40 IATOP=IAO FORT3860
50 DO 70 J=1,NATM FORT3861
SUM=0.D0 FORT3862
DO 60 K=IAO,IATOP FORT3863
SUM=SUM+U(J,K) FORT3864
60 CONTINUE FORT3865
U(J,ISUB)=SUM FORT3866
70 CONTINUE FORT3867
ISUB=ISUB+1 FORT3868
100 IAO=IATOP+1 FORT3869
RETURN FORT3870
END FORT3871
SUBROUTINE FULCHM(E,U,C,H,N5) FORT3872
REAL*8 E(N5),C(N5),U(N5,N5),H(N5,N5) FORT3873
LOGICAL*1 JGO FORT3874
C FORT3875
C SUBROUTINE TO CALCULATE COMPLETE CHARGE MATRIX. FORT3876
C FORT3877
KJ=1 FORT3878
DO 31 I=1,N5 FORT3879
JGO=.FALSE. FORT3880
IJ=KJ FORT3881
DO 32 K=1,N5 FORT3882
IF(JGO) GO TO 34 FORT3883
33 IF(I.NE.K) GO TO 35 FORT3884
JGO=.TRUE. FORT3885
E(K)=1.0D0 FORT3886
GO TO 32 FORT3887
35 E(K)=C(IJ) FORT3888
IJ=IJ+1 FORT3889
GO TO 32 FORT3890
34 IJ=IJ+K-2 FORT3891
E(K)=C(IJ) FORT3892
32 CONTINUE FORT3893
KJ=KJ+I-1 FORT3894
DO 31 J=1,N5 FORT3895
UB=0.0D0 FORT3896
DO 36 K=1,N5 FORT3897
36 UB=UB+H(K,J)*E(K) FORT3898
31 U(I,J)=2.0D0*H(I,J)*UB FORT3899
RETURN FORT3900
END FORT3901
SUBROUTINE REDCHM(U,N5) FORT3902
IMPLICIT REAL*8(A-H,O-Z) FORT3903
DIMENSION U(N5,N5) FORT3904
PARAMETER (MAXATM=500)
PARAMETER (BB=250)
PARAMETER (MXUSER=230)
PARAMETER (MXUSR2=231)
COMMON/CNTRL/CON,PEEP,COULH,NH,NA,NATM,KA,NELEC,METH,IPRINT, FORT3905
1 IPUNCH,L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT3906
LOGICAL*1 L1,L2,L3,L4,L5,ONEMAT,ITERAT FORT3907
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 FORT3911
C FORT3912
C SUBROUTINE TO CALCULATE REDUCED CHARGE MATRIX. FORT3913
C FORT3914
IF(NA.EQ.0) RETURN FORT3915
NH1=NH+1 FORT3916
IAO=NH1 FORT3917
ISUB=NH1 FORT3918
DO 128 I=1,NA FORT3919
KEYI=KEY(I) FORT3920
IF(NP(KEYI)) 122,121,122 FORT3921
121 IATOP=IAO FORT3922
GO TO 125 FORT3923
122 IF(ND(KEYI)) 124,123,124 FORT3924
123 IATOP=IAO+3 FORT3925
GO TO 125 FORT3926
124 IATOP=IAO+8 FORT3927
125 DO 127 J=1,N5 FORT3928
SUM=0.D0 FORT3929
DO 126 K=IAO,IATOP FORT3930
126 SUM=SUM+U(K,J) FORT3931
127 U(ISUB,J)=SUM FORT3932
ISUB=ISUB+1 FORT3933
128 IAO=IATOP+1 FORT3934
RETURN FORT3935
END FORT3936
SUBROUTINE CORECT(A) FORT3937
RETURN FORT3938
END
|
[ CCL Home Page ]
[ EHT ]
[ Raw Version of this page ]
| Modified: Wed Oct 13 16:00:00 1993 GMT |
| Page accessed 2154 times since Sat Apr 17 21:33:59 1999 GMT |