C C $Id: overlp.F,v 1.3 1998/07/16 16:40:10 jjv5 Exp arjan $ C C------------------------------------------------------------------------ SUBROUTINE GOVERLP(IA,NORBSA,NPQA,XISA,XIPA,XIDA,XA,YA,ZA, . IB,NORBSB,NPQB,XISB,XIPB,XIDB,XB,YB,ZB, . RAB,SABDIA,DXSABDIA,DYSABDIA,DZSABDIA) C C ROUTINE TO CALCULATE DIATOMIC OVERLAPS FOR THE S-P-D SHELLS OF C OF ATOMS A AND B. RETURNS S, P, AND D OVERLAPS IN SABDIA. C C INPUT VARIABLES: C C IA,IB = ATOMIC NUMBERS FOR ATOMS A AND B. C C NORBSA, C NORBSB = NUMBERS OF ATOMIC ORBITALS CENTERED ON ATOMS A AND B. C C NPQA,NPQB = PRINCIPAL QUANTUM NUMBERS. C C XISA,XISB = ORBITAL EXPONENTS FOR S-TYPE ATOMIC ORBITALS. C C XIPA,XIPB = ORBITAL EXPONENTS FOR P-TYPE ATOMIC ORBITALS. C C XIDA,XIDB = ORBITAL EXPONENTS FOR D-TYPE ATOMIC ORBITALS. C C XA,YA,ZA, C XB,YB,ZB = CARTESIAN COORDINATE (ANGSTROMS). C C RAB = INTERATOMIC DISTANCE (ANGSTROMS). C C C C RETURNED: C C SABDIA = DIATOMIC OVERLAP MATRIX FOR ATOMS A AND B. THE FORMAT OF C SABDIA IS SHOWN BELOW. FOR EXAMPLE, SABDIA(6,4) WOULD C CONTAIN THE OVERLAP BETWEEN THE DXZ ORBITAL OF ATOM A AND C THE PZ ORBITAL OF ATOM B. C C C FORMAT OF SABDIA: C C S(B) PX(B) PY(B) PZ(B) DZZ(B) DXZ(B) DYZ(B) DXX-YY(B) DXY(B) C -------------------------------------------------------------- C S(A) | . C PX(A) | . C PY(A) | . C PZ(A) | . C DZZ(A) | . C DXZ(A) | . . . . . . . . .(6,4) = C DYZ(A) | C DXX-YY(A) | C DXY(A) | C C C C C UTILIZES THE FOLLOWING ADDITIONAL SUBPROGRAMS: SAB C SLIMIT C RTRANS C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" #include "constants.h" real*8 SABDIA(9,9) real*8 DXSABDIA(9,9),DYSABDIA(9,9),DZSABDIA(9,9) C C AUXILIARY INFORATION FOR EFFICIENT CALCULATION OF OVERLAPS: C C EXHALF(LA,IA) = SQRT(0.5)*XILA**(NPQA+0.5), C WHERE XILA = ORBITAL EXPONENT CORRESPONDING TO LA. C C C LOCAL ARRAYS: C real*8 SABLOC(9,9),SCTEMP(9,9),C(9,9),ICODE(18) real*8 DXSCTEMP(9,9),DYSCTEMP(9,9),DZSCTEMP(9,9) integer NIJSC(6),ISC(6,81),JSC(6,81),NKIJSC(6,81),KIJSC(6,81,9), . NKC(6,9),KC(6,9,9) real*8 DSABLOC(9,9),DXSABLOC(9,9),DYSABLOC(9,9),DZSABLOC(9,9) real*8 DXC(9,9),DYC(9,9),DZC(9,9) LOGICAL NEWEXP C LOGICAL FIRST DATA FIRST /.TRUE./ C C CODE FOR THE SIX POSSIBLE VALUES OF NORBSA+NORBSB: C DATA ICODE /0,1,0,0,2,0,0,3,0,4,0,0,5,0,0,0,0,6/ C SAVE FIRST,ICODE,NIJSC,ISC,JSC,NKIJSC,KIJSC,NKC,KC C C ON THE FIRST CALL, SET LIMITS TO TAKE ADVANTAGE OF THE SPARSITY C OF SABLOC AND C. C IF(FIRST)THEN FIRST = .FALSE. CALL GSLIMIT(NIJSC,ISC,JSC,NKIJSC,KIJSC,NKC,KC) ENDIF C C BEFORE DOING ANY OVERLAP CALCULATIONS, CHECK TO SEE IF THE ORBITAL C EXPONENTS AND INTERATOMIC DISTANCE ARE SUCH THAT THE OVERLAPS WOULD C BE NEGLIGIBLE. C XMIN = XISA + XISB IF(NORBSA.GT.1) XMIN = MIN(XMIN,XIPA+XISB) IF(NORBSA.GT.4) XMIN = MIN(XMIN,XIDA+XISB) IF(NORBSB.GT.1) XMIN = MIN(XMIN,XISA+XIPB) IF(NORBSB.GT.1) XMIN = MIN(XMIN,XISA+XIDB) IF(NORBSA.GT.1.AND.NORBSB.GT.1)THEN XMIN = MIN(XMIN,XIPA+XIPB) IF(NORBSA.GT.4) XMIN = MIN(XMIN,XIDA+XIPB) IF(NORBSB.GT.4) XMIN = MIN(XMIN,XIPA+XIDB) IF(NORBSA.GT.4.AND.NORBSB.GT.4) XMIN = MIN(XMIN,XIDA+XIDB) ENDIF XTEST = RAB*XMIN IF(XTEST.GT.25.0D0)THEN DO 20 I=1,NORBSA DO 10 J=1,NORBSB SABDIA(I,J) = 0.0D0 DXSABDIA(I,J) = 0.0D0 DYSABDIA(I,J) = 0.0D0 DZSABDIA(I,J) = 0.0D0 SABLOC(I,J) = 0.0D0 DSABLOC(I,J) = 0.0D0 DXSABLOC(I,J) = 0.0D0 DYSABLOC(I,J) = 0.0D0 DZSABLOC(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE RETURN ENDIF C C C ROTATION OF LOCAL OVERLAPS TO MOLECULAR FRAME REQUIRES THAT NORBSA C BE GREATER THAN OR EQUAL TO NORBSB. IF THIS ISN'T THE CASE, THEN C TEMPORARILY SWAP EVERYTHING. C ISWAP = 0 IF(NORBSB.GT.NORBSA)THEN ISWAP = 1 ITEMP = IA IA = IB IB = ITEMP XTEMP = XISA XISA = XISB XISB = XTEMP XTEMP = XIPA XIPA = XIPB XIPB = XTEMP XTEMP = XIDA XIDA = XIDB XIDB = XTEMP XTEMP = XA XA = XB XB = XTEMP YTEMP = YA YA = YB YB = YTEMP ZTEMP = ZA ZA = ZB ZB = ZTEMP NPTEMP = NPQA NPQA = NPQB NPQB = NPTEMP NTEMP = NORBSA NORBSA = NORBSB NORBSB = NTEMP ENDIF C C NOW COMPUTE THE NONZERO LOCAL OVERLAPS SABLOC. THESE ARE DEFINED C IN A FRAME WHEREIN ATOM B IS DIRECTLY ABOVE ATOM A (I.E., XA=XB, C YA=YB, AND ZB-ZA=RAB). THE LOCAL FRAME OVERLAPS ARE THEN ROTATED C TO YIELD THE MOLECULAR FRAME OVERLAPS SABDIA. C R = RAB/Bohr2Ang NCENTR = 2 EXHSA = EXHALF(0,IA) EXHPA = EXHALF(1,IA) EXHDA = EXHALF(2,IA) EXHSB = EXHALF(0,IB) EXHPB = EXHALF(1,IB) EXHDB = EXHALF(2,IB) RNAB = R**(NPQA+NPQB+1) C C : C LA = 0 LB = 0 M = 0 NEWEXP = .TRUE. call gfctsab(NCENTR,NPQA,LA,XISA,EXHSA, . NPQB,LB,XISB,EXHSB, . NEWEXP,M,R,RNAB,SABLOC(1,1),DSABLOC(1,1)) C C : C IF(NORBSA.GT.1)THEN LA = 1 LB = 0 M = 0 NEWEXP = .TRUE. call gfctsab(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XISB,EXHSB, . NEWEXP,M,R,RNAB,SABLOC(4,1),DSABLOC(4,1)) ENDIF C C : C IF(NORBSA.GT.4)THEN LA = 2 LB = 0 M = 0 NEWEXP = .TRUE. call gfctsab(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XISB,EXHSB, . NEWEXP,M,R,RNAB,SABLOC(5,1),DSABLOC(5,1)) ENDIF C C , , AND : C IF(NORBSA.GT.4.AND.NORBSB.GT.1)THEN LA = 2 LB = 1 M = 1 NEWEXP = .TRUE. call gfctsab(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB,SABLOC(6,2),DSABLOC(6,2)) SABLOC(7,3) = SABLOC(6,2) DSABLOC(7,3) = DSABLOC(6,2) M = 0 NEWEXP = .FALSE. call gfctsab(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB,SABLOC(5,4),DSABLOC(5,4)) SABLOC(5,4) = -SABLOC(5,4) DSABLOC(5,4) = -DSABLOC(5,4) ENDIF C C : C IF(NORBSB.GT.1)THEN IF(IA.EQ.IB)THEN SABLOC(1,4) = -SABLOC(4,1) DSABLOC(1,4) = -DSABLOC(4,1) ELSE LA = 0 LB = 1 M = 0 NEWEXP = .TRUE. call gfctsab(NCENTR,NPQA,LA,XISA,EXHSA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB,SABLOC(1,4),DSABLOC(1,4)) SABLOC(1,4) = -SABLOC(1,4) DSABLOC(1,4) = -DSABLOC(1,4) ENDIF ENDIF C C : C IF(NORBSB.GT.4)THEN IF(IA.EQ.IB)THEN SABLOC(1,5) = SABLOC(5,1) DSABLOC(1,5) = DSABLOC(5,1) ELSE LA = 0 LB = 2 M = 0 NEWEXP = .TRUE. call gfctsab(NCENTR,NPQA,LA,XISA,EXHSA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB,SABLOC(1,5),DSABLOC(1,5)) ENDIF ENDIF C C , , AND : C IF(NORBSA.GT.1.AND.NORBSB.GT.4)THEN IF(IA.EQ.IB)THEN SABLOC(2,6) = -SABLOC(6,2) DSABLOC(2,6) = -DSABLOC(6,2) SABLOC(4,5) = -SABLOC(5,4) DSABLOC(4,5) = -DSABLOC(5,4) ELSE LA = 1 LB = 2 M = 1 NEWEXP = .TRUE. call gfctsab(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB,SABLOC(2,6),DSABLOC(2,6)) SABLOC(2,6) = -SABLOC(2,6) DSABLOC(2,6) = -DSABLOC(2,6) M = 0 NEWEXP = .FALSE. call gfctsab(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB,SABLOC(4,5),DSABLOC(4,5)) ENDIF SABLOC(3,7) = SABLOC(2,6) DSABLOC(3,7) = DSABLOC(2,6) ENDIF C C , , : C IF(NORBSA.GT.1.AND.NORBSB.GT.1)THEN LA = 1 LB = 1 M = 1 NEWEXP = .TRUE. call gfctsab(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB,SABLOC(2,2),DSABLOC(2,2)) NEWEXP = .FALSE. SABLOC(3,3) = SABLOC(2,2) DSABLOC(3,3) = DSABLOC(2,2) M = 0 call gfctsab(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB,SABLOC(4,4),DSABLOC(4,4)) SABLOC(4,4) = -SABLOC(4,4) DSABLOC(4,4) = -DSABLOC(4,4) ENDIF C C THE FIVE OVERLAPS: C IF(NORBSA.GT.4.AND.NORBSB.GT.4)THEN LA = 2 LB = 2 M = 0 NEWEXP = .TRUE. call gfctsab(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB,SABLOC(5,5),DSABLOC(5,5)) M = 1 NEWEXP = .FALSE. call gfctsab(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB,SABLOC(6,6),DSABLOC(6,6)) SABLOC(6,6) = -SABLOC(6,6) DSABLOC(6,6) = -DSABLOC(6,6) SABLOC(7,7) = SABLOC(6,6) DSABLOC(7,7) = DSABLOC(6,6) M = 2 call gfctsab(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB,SABLOC(8,8),DSABLOC(8,8)) SABLOC(9,9) = SABLOC(8,8) DSABLOC(9,9) = DSABLOC(8,8) ENDIF c c D D c -- --> -- c DR Dx c IF(ISWAP.EQ.1)THEN do i=1,9 do j=1,9 DXSABLOC(j,i)=-DSABLOC(j,i)*(XA-XB)/RAB/Bohr2Ang DYSABLOC(j,i)=-DSABLOC(j,i)*(YA-YB)/RAB/Bohr2Ang DZSABLOC(j,i)=-DSABLOC(j,i)*(ZA-ZB)/RAB/Bohr2Ang enddo enddo ELSE do i=1,9 do j=1,9 DXSABLOC(j,i)=DSABLOC(j,i)*(XA-XB)/RAB/Bohr2Ang DYSABLOC(j,i)=DSABLOC(j,i)*(YA-YB)/RAB/Bohr2Ang DZSABLOC(j,i)=DSABLOC(j,i)*(ZA-ZB)/RAB/Bohr2Ang enddo enddo ENDIF C C NOW DETERMINE A TRANSFORMATION MATRIX C THAT RELATES THE LOCAL FRAME C OVERLAPS TO THE TRUE DIATOMIC OVERLAPS. THE MATRIX C IS COMPUTED C FROM A ROTATION MATRIX THAT RELATES THE TRUE COORDINATE SYSTEM TO C A COORDINATE SYSTEM WHICH IS CONSISTENT WITH THE LOCAL FRAME. C NORBS = MAX(NORBSA,NORBSB) IF(NORBS.EQ.1)THEN SABDIA(1,1) = SABLOC(1,1) DXSABLOC(1,1)=DSABLOC(1,1)*(XA-XB)/RAB/Bohr2Ang DYSABLOC(1,1)=DSABLOC(1,1)*(YA-YB)/RAB/Bohr2Ang DZSABLOC(1,1)=DSABLOC(1,1)*(ZA-ZB)/RAB/Bohr2Ang DXSABDIA(1,1)=DXSABLOC(1,1) DYSABDIA(1,1)=DYSABLOC(1,1) DZSABDIA(1,1)=DZSABLOC(1,1) ELSE CALL GRTRANS(XA,YA,ZA,XB,YB,ZB,RAB,NORBS,C,DXC,DYC,DZC,iswap) C C IF THERE ARE D ORBITALS THEN COMPUTE THE MOLECULAR FRAME OVERLAPS C USING THE GENERAL FORMULA, C C SABDIA = C'*SABLOC*C. C C IF WE ONLY HAVE S AND P ORBITALS, THEN DO ROTATIONS BY HAND. C IF(NORBSA.GT.4.OR.NORBSB.GT.4)THEN C C USE GENERAL FORMULA. C IAB = ICODE(NORBSA+NORBSB) DO 200 IJ=1,NIJSC(IAB) I = ISC(IAB,IJ) J = JSC(IAB,IJ) TMPIJ = 0.0D0 DXTMPIJ = 0.0D0 DYTMPIJ = 0.0D0 DZTMPIJ = 0.0D0 DO 160 K=1,NKIJSC(IAB,IJ) KINDX = KIJSC(IAB,IJ,K) TMPIJ = TMPIJ + SABLOC(I,KINDX)*C(KINDX,J) DXTMPIJ = DXTMPIJ + DXSABLOC(I,KINDX)*C(KINDX,J) . + SABLOC(I,KINDX)*DXC(KINDX,J) DYTMPIJ = DYTMPIJ + DYSABLOC(I,KINDX)*C(KINDX,J) . + SABLOC(I,KINDX)*DYC(KINDX,J) DZTMPIJ = DZTMPIJ + DZSABLOC(I,KINDX)*C(KINDX,J) . + SABLOC(I,KINDX)*DZC(KINDX,J) 160 CONTINUE SCTEMP(I,J) = TMPIJ DXSCTEMP(I,J) = DXTMPIJ DYSCTEMP(I,J) = DYTMPIJ DZSCTEMP(I,J) = DZTMPIJ 200 CONTINUE DO 260 I=1,NORBSA DO 240 J=1,NORBSB SIJ = 0.0D0 DXSIJ = 0.0D0 DYSIJ = 0.0D0 DZSIJ = 0.0D0 DO 220 K=1,NKC(IAB,I) KINDX = KC(IAB,I,K) SIJ = SIJ + C(KINDX,I)*SCTEMP(KINDX,J) DXSIJ = DXSIJ + DXC(KINDX,I)*SCTEMP(KINDX,J) . + C(KINDX,I)*DXSCTEMP(KINDX,J) DYSIJ = DYSIJ + DYC(KINDX,I)*SCTEMP(KINDX,J) . + C(KINDX,I)*DYSCTEMP(KINDX,J) DZSIJ = DZSIJ + DZC(KINDX,I)*SCTEMP(KINDX,J) . + C(KINDX,I)*DZSCTEMP(KINDX,J) 220 CONTINUE SABDIA(I,J) = SIJ DXSABDIA(I,J) = DXSIJ DYSABDIA(I,J) = DYSIJ DZSABDIA(I,J) = DZSIJ 240 CONTINUE 260 CONTINUE ELSE C C MANUAL ROTATIONS. C SABDIA(1,1) = SABLOC(1,1) DXSABDIA(1,1) = DXSABLOC(1,1) DYSABDIA(1,1) = DYSABLOC(1,1) DZSABDIA(1,1) = DZSABLOC(1,1) IF(NORBSB.GT.1)THEN SABDIA(1,2) = SABLOC(1,4)*C(4,2) DXSABDIA(1,2) = DXSABLOC(1,4)*C(4,2) + SABLOC(1,4)*DXC(4,2) DYSABDIA(1,2) = DYSABLOC(1,4)*C(4,2) + SABLOC(1,4)*DYC(4,2) DZSABDIA(1,2) = DZSABLOC(1,4)*C(4,2) + SABLOC(1,4)*DZC(4,2) SABDIA(1,3) = SABLOC(1,4)*C(4,3) DXSABDIA(1,3) = DXSABLOC(1,4)*C(4,3) + SABLOC(1,4)*DXC(4,3) DYSABDIA(1,3) = DYSABLOC(1,4)*C(4,3) + SABLOC(1,4)*DYC(4,3) DZSABDIA(1,3) = DZSABLOC(1,4)*C(4,3) + SABLOC(1,4)*DZC(4,3) SABDIA(1,4) = SABLOC(1,4)*C(4,4) DXSABDIA(1,4) = DXSABLOC(1,4)*C(4,4) + SABLOC(1,4)*DXC(4,4) DYSABDIA(1,4) = DYSABLOC(1,4)*C(4,4) + SABLOC(1,4)*DYC(4,4) DZSABDIA(1,4) = DZSABLOC(1,4)*C(4,4) + SABLOC(1,4)*DZC(4,4) ENDIF IF(NORBSA.GT.1)THEN SABDIA(2,1) = SABLOC(4,1)*C(4,2) DXSABDIA(2,1) = DXSABLOC(4,1)*C(4,2) + SABLOC(4,1)*DXC(4,2) DYSABDIA(2,1) = DYSABLOC(4,1)*C(4,2) + SABLOC(4,1)*DYC(4,2) DZSABDIA(2,1) = DZSABLOC(4,1)*C(4,2) + SABLOC(4,1)*DZC(4,2) SABDIA(3,1) = SABLOC(4,1)*C(4,3) DXSABDIA(3,1) = DXSABLOC(4,1)*C(4,3) + SABLOC(4,1)*DXC(4,3) DYSABDIA(3,1) = DYSABLOC(4,1)*C(4,3) + SABLOC(4,1)*DYC(4,3) DZSABDIA(3,1) = DZSABLOC(4,1)*C(4,3) + SABLOC(4,1)*DZC(4,3) SABDIA(4,1) = SABLOC(4,1)*C(4,4) DXSABDIA(4,1) = DXSABLOC(4,1)*C(4,4) + SABLOC(4,1)*DXC(4,4) DYSABDIA(4,1) = DYSABLOC(4,1)*C(4,4) + SABLOC(4,1)*DYC(4,4) DZSABDIA(4,1) = DZSABLOC(4,1)*C(4,4) + SABLOC(4,1)*DZC(4,4) ENDIF IF(NORBSA.GT.1.AND.NORBSB.GT.1)THEN SABDIA(2,2) = SABLOC(2,2)*C(2,2)*C(2,2) . + SABLOC(3,3)*C(3,2)*C(3,2) . + SABLOC(4,4)*C(4,2)*C(4,2) DXSABDIA(2,2) = DXSABLOC(2,2)*C(2,2)*C(2,2) . + SABLOC(2,2)*DXC(2,2)*C(2,2) . + SABLOC(2,2)*C(2,2)*DXC(2,2) . + DXSABLOC(3,3)*C(3,2)*C(3,2) . + SABLOC(3,3)*DXC(3,2)*C(3,2) . + SABLOC(3,3)*C(3,2)*DXC(3,2) . + DXSABLOC(4,4)*C(4,2)*C(4,2) . + SABLOC(4,4)*DXC(4,2)*C(4,2) . + SABLOC(4,4)*C(4,2)*DXC(4,2) DYSABDIA(2,2) = DYSABLOC(2,2)*C(2,2)*C(2,2) . + SABLOC(2,2)*DYC(2,2)*C(2,2) . + SABLOC(2,2)*C(2,2)*DYC(2,2) . + DYSABLOC(3,3)*C(3,2)*C(3,2) . + SABLOC(3,3)*DYC(3,2)*C(3,2) . + SABLOC(3,3)*C(3,2)*DYC(3,2) . + DYSABLOC(4,4)*C(4,2)*C(4,2) . + SABLOC(4,4)*DYC(4,2)*C(4,2) . + SABLOC(4,4)*C(4,2)*DYC(4,2) DZSABDIA(2,2) = DZSABLOC(2,2)*C(2,2)*C(2,2) . + SABLOC(2,2)*DZC(2,2)*C(2,2) . + SABLOC(2,2)*C(2,2)*DZC(2,2) . + DZSABLOC(3,3)*C(3,2)*C(3,2) . + SABLOC(3,3)*DZC(3,2)*C(3,2) . + SABLOC(3,3)*C(3,2)*DZC(3,2) . + DZSABLOC(4,4)*C(4,2)*C(4,2) . + SABLOC(4,4)*DZC(4,2)*C(4,2) . + SABLOC(4,4)*C(4,2)*DZC(4,2) C SABDIA(2,3) = SABLOC(2,2)*C(2,2)*C(2,3) . + SABLOC(3,3)*C(3,2)*C(3,3) . + SABLOC(4,4)*C(4,2)*C(4,3) DXSABDIA(2,3) = DXSABLOC(2,2)*C(2,2)*C(2,3) . + SABLOC(2,2)*DXC(2,2)*C(2,3) . + SABLOC(2,2)*C(2,2)*DXC(2,3) . + DXSABLOC(3,3)*C(3,2)*C(3,3) . + SABLOC(3,3)*DXC(3,2)*C(3,3) . + SABLOC(3,3)*C(3,2)*DXC(3,3) . + DXSABLOC(4,4)*C(4,2)*C(4,3) . + SABLOC(4,4)*DXC(4,2)*C(4,3) . + SABLOC(4,4)*C(4,2)*DXC(4,3) DYSABDIA(2,3) = DYSABLOC(2,2)*C(2,2)*C(2,3) . + SABLOC(2,2)*DYC(2,2)*C(2,3) . + SABLOC(2,2)*C(2,2)*DYC(2,3) . + DYSABLOC(3,3)*C(3,2)*C(3,3) . + SABLOC(3,3)*DYC(3,2)*C(3,3) . + SABLOC(3,3)*C(3,2)*DYC(3,3) . + DYSABLOC(4,4)*C(4,2)*C(4,3) . + SABLOC(4,4)*DYC(4,2)*C(4,3) . + SABLOC(4,4)*C(4,2)*DYC(4,3) DZSABDIA(2,3) = DZSABLOC(2,2)*C(2,2)*C(2,3) . + SABLOC(2,2)*DZC(2,2)*C(2,3) . + SABLOC(2,2)*C(2,2)*DZC(2,3) . + DZSABLOC(3,3)*C(3,2)*C(3,3) . + SABLOC(3,3)*DZC(3,2)*C(3,3) . + SABLOC(3,3)*C(3,2)*DZC(3,3) . + DZSABLOC(4,4)*C(4,2)*C(4,3) . + SABLOC(4,4)*DZC(4,2)*C(4,3) . + SABLOC(4,4)*C(4,2)*DZC(4,3) C SABDIA(2,4) = SABLOC(2,2)*C(2,2)*C(2,4) . + SABLOC(3,3)*C(3,2)*C(3,4) . + SABLOC(4,4)*C(4,2)*C(4,4) DXSABDIA(2,4) = DXSABLOC(2,2)*C(2,2)*C(2,4) . + SABLOC(2,2)*DXC(2,2)*C(2,4) . + SABLOC(2,2)*C(2,2)*DXC(2,4) . + DXSABLOC(3,3)*C(3,2)*C(3,4) . + SABLOC(3,3)*DXC(3,2)*C(3,4) . + SABLOC(3,3)*C(3,2)*DXC(3,4) . + DXSABLOC(4,4)*C(4,2)*C(4,4) . + SABLOC(4,4)*DXC(4,2)*C(4,4) . + SABLOC(4,4)*C(4,2)*DXC(4,4) DYSABDIA(2,4) = DYSABLOC(2,2)*C(2,2)*C(2,4) . + SABLOC(2,2)*DYC(2,2)*C(2,4) . + SABLOC(2,2)*C(2,2)*DYC(2,4) . + DYSABLOC(3,3)*C(3,2)*C(3,4) . + SABLOC(3,3)*DYC(3,2)*C(3,4) . + SABLOC(3,3)*C(3,2)*DYC(3,4) . + DYSABLOC(4,4)*C(4,2)*C(4,4) . + SABLOC(4,4)*DYC(4,2)*C(4,4) . + SABLOC(4,4)*C(4,2)*DYC(4,4) DZSABDIA(2,4) = DZSABLOC(2,2)*C(2,2)*C(2,4) . + SABLOC(2,2)*DZC(2,2)*C(2,4) . + SABLOC(2,2)*C(2,2)*DZC(2,4) . + DZSABLOC(3,3)*C(3,2)*C(3,4) . + SABLOC(3,3)*DZC(3,2)*C(3,4) . + SABLOC(3,3)*C(3,2)*DZC(3,4) . + DZSABLOC(4,4)*C(4,2)*C(4,4) . + SABLOC(4,4)*DZC(4,2)*C(4,4) . + SABLOC(4,4)*C(4,2)*DZC(4,4) C SABDIA(3,2) = SABDIA(2,3) DXSABDIA(3,2) = DXSABDIA(2,3) DYSABDIA(3,2) = DYSABDIA(2,3) DZSABDIA(3,2) = DZSABDIA(2,3) C SABDIA(3,3) = SABLOC(2,2)*C(2,3)*C(2,3) . + SABLOC(3,3)*C(3,3)*C(3,3) . + SABLOC(4,4)*C(4,3)*C(4,3) DXSABDIA(3,3) = DXSABLOC(2,2)*C(2,3)*C(2,3) . + SABLOC(2,2)*DXC(2,3)*C(2,3) . + SABLOC(2,2)*C(2,3)*DXC(2,3) . + DXSABLOC(3,3)*C(3,3)*C(3,3) . + SABLOC(3,3)*DXC(3,3)*C(3,3) . + SABLOC(3,3)*C(3,3)*DXC(3,3) . + DXSABLOC(4,4)*C(4,3)*C(4,3) . + SABLOC(4,4)*DXC(4,3)*C(4,3) . + SABLOC(4,4)*C(4,3)*DXC(4,3) DYSABDIA(3,3) = DYSABLOC(2,2)*C(2,3)*C(2,3) . + SABLOC(2,2)*DYC(2,3)*C(2,3) . + SABLOC(2,2)*C(2,3)*DYC(2,3) . + DYSABLOC(3,3)*C(3,3)*C(3,3) . + SABLOC(3,3)*DYC(3,3)*C(3,3) . + SABLOC(3,3)*C(3,3)*DYC(3,3) . + DYSABLOC(4,4)*C(4,3)*C(4,3) . + SABLOC(4,4)*DYC(4,3)*C(4,3) . + SABLOC(4,4)*C(4,3)*DYC(4,3) DZSABDIA(3,3) = DZSABLOC(2,2)*C(2,3)*C(2,3) . + SABLOC(2,2)*DZC(2,3)*C(2,3) . + SABLOC(2,2)*C(2,3)*DZC(2,3) . + DZSABLOC(3,3)*C(3,3)*C(3,3) . + SABLOC(3,3)*DZC(3,3)*C(3,3) . + SABLOC(3,3)*C(3,3)*DZC(3,3) . + DZSABLOC(4,4)*C(4,3)*C(4,3) . + SABLOC(4,4)*DZC(4,3)*C(4,3) . + SABLOC(4,4)*C(4,3)*DZC(4,3) C SABDIA(3,4) = SABLOC(2,2)*C(2,3)*C(2,4) . + SABLOC(3,3)*C(3,3)*C(3,4) . + SABLOC(4,4)*C(4,3)*C(4,4) DXSABDIA(3,4) = DXSABLOC(2,2)*C(2,3)*C(2,4) . + SABLOC(2,2)*DXC(2,3)*C(2,4) . + SABLOC(2,2)*C(2,3)*DXC(2,4) . + DXSABLOC(3,3)*C(3,3)*C(3,4) . + SABLOC(3,3)*DXC(3,3)*C(3,4) . + SABLOC(3,3)*C(3,3)*DXC(3,4) . + DXSABLOC(4,4)*C(4,3)*C(4,4) . + SABLOC(4,4)*DXC(4,3)*C(4,4) . + SABLOC(4,4)*C(4,3)*DXC(4,4) DYSABDIA(3,4) = DYSABLOC(2,2)*C(2,3)*C(2,4) . + SABLOC(2,2)*DYC(2,3)*C(2,4) . + SABLOC(2,2)*C(2,3)*DYC(2,4) . + DYSABLOC(3,3)*C(3,3)*C(3,4) . + SABLOC(3,3)*DYC(3,3)*C(3,4) . + SABLOC(3,3)*C(3,3)*DYC(3,4) . + DYSABLOC(4,4)*C(4,3)*C(4,4) . + SABLOC(4,4)*DYC(4,3)*C(4,4) . + SABLOC(4,4)*C(4,3)*DYC(4,4) DZSABDIA(3,4) = DZSABLOC(2,2)*C(2,3)*C(2,4) . + SABLOC(2,2)*DZC(2,3)*C(2,4) . + SABLOC(2,2)*C(2,3)*DZC(2,4) . + DZSABLOC(3,3)*C(3,3)*C(3,4) . + SABLOC(3,3)*DZC(3,3)*C(3,4) . + SABLOC(3,3)*C(3,3)*DZC(3,4) . + DZSABLOC(4,4)*C(4,3)*C(4,4) . + SABLOC(4,4)*DZC(4,3)*C(4,4) . + SABLOC(4,4)*C(4,3)*DZC(4,4) C SABDIA(4,2) = SABDIA(2,4) DXSABDIA(4,2) = DXSABDIA(2,4) DYSABDIA(4,2) = DYSABDIA(2,4) DZSABDIA(4,2) = DZSABDIA(2,4) SABDIA(4,3) = SABDIA(3,4) DXSABDIA(4,3) = DXSABDIA(3,4) DYSABDIA(4,3) = DYSABDIA(3,4) DZSABDIA(4,3) = DZSABDIA(3,4) C SABDIA(4,4) = SABLOC(2,2)*C(2,4)*C(2,4) . + SABLOC(3,3)*C(3,4)*C(3,4) . + SABLOC(4,4)*C(4,4)*C(4,4) DXSABDIA(4,4) = DXSABLOC(2,2)*C(2,4)*C(2,4) . + SABLOC(2,2)*DXC(2,4)*C(2,4) . + SABLOC(2,2)*C(2,4)*DXC(2,4) . + DXSABLOC(3,3)*C(3,4)*C(3,4) . + SABLOC(3,3)*DXC(3,4)*C(3,4) . + SABLOC(3,3)*C(3,4)*DXC(3,4) . + DXSABLOC(4,4)*C(4,4)*C(4,4) . + SABLOC(4,4)*DXC(4,4)*C(4,4) . + SABLOC(4,4)*C(4,4)*DXC(4,4) DYSABDIA(4,4) = DYSABLOC(2,2)*C(2,4)*C(2,4) . + SABLOC(2,2)*DYC(2,4)*C(2,4) . + SABLOC(2,2)*C(2,4)*DYC(2,4) . + DYSABLOC(3,3)*C(3,4)*C(3,4) . + SABLOC(3,3)*DYC(3,4)*C(3,4) . + SABLOC(3,3)*C(3,4)*DYC(3,4) . + DYSABLOC(4,4)*C(4,4)*C(4,4) . + SABLOC(4,4)*DYC(4,4)*C(4,4) . + SABLOC(4,4)*C(4,4)*DYC(4,4) DZSABDIA(4,4) = DZSABLOC(2,2)*C(2,4)*C(2,4) . + SABLOC(2,2)*DZC(2,4)*C(2,4) . + SABLOC(2,2)*C(2,4)*DZC(2,4) . + DZSABLOC(3,3)*C(3,4)*C(3,4) . + SABLOC(3,3)*DZC(3,4)*C(3,4) . + SABLOC(3,3)*C(3,4)*DZC(3,4) . + DZSABLOC(4,4)*C(4,4)*C(4,4) . + SABLOC(4,4)*DZC(4,4)*C(4,4) . + SABLOC(4,4)*C(4,4)*DZC(4,4) ENDIF ENDIF ENDIF C C IF SWAP WAS MADE PREVIOUSLY, RETURN INPUT VARIABLES TO ORIGINAL C VALUES AND TRANSPOSE THE DIATOMIC OVERLAP MATRIX. C IF(ISWAP.EQ.1)THEN ITEMP = IA IA = IB IB = ITEMP XTEMP = XISA XISA = XISB XISB = XTEMP XTEMP = XIPA XIPA = XIPB XIPB = XTEMP XTEMP = XIDA XIDA = XIDB XIDB = XTEMP XTEMP = XA XA = XB XB = XTEMP YTEMP = YA YA = YB YB = YTEMP ZTEMP = ZA ZA = ZB ZB = ZTEMP NPTEMP = NPQA NPQA = NPQB NPQB = NPTEMP NTEMP = NORBSA NORBSA = NORBSB NORBSB = NTEMP DO 500 I=2,NORBSB DO 480 J=1,I-1 STEMP = SABDIA(I,J) DXSTEMP = DXSABDIA(I,J) DYSTEMP = DYSABDIA(I,J) DZSTEMP = DZSABDIA(I,J) SABDIA(I,J) = SABDIA(J,I) DXSABDIA(I,J) = DXSABDIA(J,I) DYSABDIA(I,J) = DYSABDIA(J,I) DZSABDIA(I,J) = DZSABDIA(J,I) SABDIA(J,I) = STEMP DXSABDIA(J,I) = DXSTEMP DYSABDIA(J,I) = DYSTEMP DZSABDIA(J,I) = DZSTEMP 480 CONTINUE 500 CONTINUE ENDIF RETURN END C C C SUBROUTINE GFCTSAB(NCENTR,NA,LA,XIA,EXHLA, . NB,LB,XIB,EXHLB, . NEWEXP,M,RAB,RNAB,SAB,DSAB) C C PORTABLE OVERLAP FUNCTION. COMPUTES THE OVERLAP INTEGRAL BETWEEN C SLATER-TYPE ORBITALS A AND B BY TRANSFORMATION TO PROLATE SPHEROIDAL C COORDINATES. IN TERMS OF CARTESIAN COORDINATES, THIS MEANS THAT C THE XYZ COORDINATE SYSTEMS OF CENTERS A AND B ARE MIRROR IMAGES, C WITH THE Z-AXIS OF A POINTING TOWARD B AND VICE VERSA. C C HANDLES OVERLAP FOR ORBITALS OF TYPE S, P, D, F, ..., ETC. C C C INPUT: C C NCENTR = THE NUMBER OF CENTERS INVOLVED (NCENTR=2 MEANS ORBITALS C A AND B ARE CENTERED ON DIFFERENT ATOMS, NCENTR=1 MEANS C THEY ARE CENTERED ON THE SAME ATOM. DEFAULT=2). C NA,NB = PRINCIPAL QUANTUM NUMBERS FOR THE ORBITALS. C LA,LB = ANGULAR MOMENTUM QUANTUM NUMBERS FOR THE ORBITALS. C XIA,XIB = EXPONENTS FOR THE RADIAL TERMS IN THE STO. C EXHLA = SQRT(0.5)*XIA**(NA+0.5) C EXHLB = SQRT(0.5)*XIB**(NA+0.5) C NEWEXP = .TRUE. IF EITHER XIA OR XIB HAS CHANGED FROM PREVIOUS C CALL FOR A GIVEN ATOM, .FALSE. OTHERWISE. C M = MAGNETIC QUANTUM NUMBER (MUST BE THE SAME FOR BOTH ORBITALS). C RAB = DISTANCE BETWEEN THE ATOMIC CENTERS IN ATOMIC UNITS C (ATOMIC UNITS = ANGTROMS/Bohr2Ang) C Ken Changed back to original MOPAC value: C (ATOMIC UNITS = ANGTROMS/0.529167) C RNAB = RAB**(NA+NB+1) C C C RETURNED: C C SAB = OVERLAP BETWEEN ORBITALS A AND B. C C C C STO'S ARE OF THE FORM: C C STO(N,L,M,XI,R,THETA,PHI) = NORM*R**(N-1)*EXP(-XI*R)*Y(L,M,THETA,PHI); C C WHERE, C C N,L,M = QUANTUM NUMBERS; C XI = ORBITAL EXPONENT; C R,THETA,PHI = SPHERICAL POLAR COORDINATES; C NORM = RADIAL FUNCTION NORMALIZATION FACTOR; C Y(L,M,THETA,PHI) = NORMALIZED, COMPLEX SPHERICAL HARMONIC FUNCTION. C C C THE FORMULATION USED IS TAKEN DIRECTLY FROM: C J. P. P. STEWART, JOURNAL OF COMPUTER-AIDED MOLECULAR DESIGN, C VOL. 4, PP 1-105 (1990). THE FINAL FORMULA FOR THE OVERLAP C APPEARS ON PAGE 28 OF THIS ARTICLE. ONE CORRECTION WAS NECESSARY C FOR THE EXPONENT IEXP -- SEE BELOW IN CODE. C C PROGRAMMED BY S. L. DIXON, OCT., 1991. C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL NEWEXP real*8 FACT(0:15),CTERM(0:4,0:4,0:2),BINOM(0:10,0:10), . ATERM(0:15),BTERM(0:15),ABTERM(0:15,0:15),NEG1(-15:15), . DATERM(0:15),DBTERM(0:15),DABTERM(0:15,0:15) SAVE RPLUS,RMINUS,ABTERM,NEG1 C C INITIALIZE FACTORIAL VALUES (UP TO 15!). C DATA FACT /1.0D0,1.0D0,2.0D0,6.0D0,24.0D0,120.0D0,720.0D0, . 5040.0D0,40320.0D0,362880.0D0,36288.0D2,399168.0D2, . 4790016.0D2,62270208.0D2,871782912.0D2,1307674368.0D3/ C C INITIALIZE SPHERICAL HARMONIC NORMALIZATION FACTORS -- UP TO C AN ANGULAR MOMENTUM QUANTUM NUMBER OF 4. CORRESPONDS TO THE C TERM C(L,M,J) OF STEWART REFERENCE. C DATA CTERM(0,0,0) / 1.000000000000D0/ DATA CTERM(1,0,0) / 1.000000000000D0/ DATA CTERM(1,1,0) / 0.707106781187D0/ DATA CTERM(2,0,0) / 1.500000000000D0/ DATA CTERM(2,0,1) /-0.500000000000D0/ DATA CTERM(2,1,0) / 1.224744871392D0/ DATA CTERM(2,2,0) / 0.612372435696D0/ DATA CTERM(3,0,0) / 2.500000000000D0/ DATA CTERM(3,0,1) /-1.500000000000D0/ DATA CTERM(3,1,0) / 2.165063509461D0/ DATA CTERM(3,1,1) /-0.433012701892D0/ DATA CTERM(3,2,0) / 1.369306393763D0/ DATA CTERM(3,3,0) / 0.559016994375D0/ DATA CTERM(4,0,0) / 4.375000000000D0/ DATA CTERM(4,0,1) /-3.750000000000D0/ DATA CTERM(4,0,2) / 0.375000000000D0/ DATA CTERM(4,1,0) / 3.913118960625D0/ DATA CTERM(4,1,1) /-1.677050983125D0/ DATA CTERM(4,2,0) / 2.766992952647D0/ DATA CTERM(4,2,1) /-0.395284707521D0/ DATA CTERM(4,3,0) / 1.479019945775D0/ DATA CTERM(4,4,0) / 0.522912516584D0/ C C INITIALIZE BINOMIAL COEFFICIENTS: C C BINOM(N1,N2) = (N1!)/((N1-N2)!(N2!)) FOR N1.GE.N2 C C ASSIGN VALUES UP TO N1=10. C DATA BINOM /11*1.0D0,0.0D0,1.0D0,2.0D0,3.0D0,4.0D0,5.0D0,6.0D0, . 7.0D0,8.0D0,9.0D0,10.0D0,2*0.0D0,1.0D0,3.0D0,6.0D0, . 10.0D0,15.0D0,21.0D0,28.0D0,36.0D0,45.0D0,3*0.0D0, . 1.0D0,4.0D0,10.0D0,20.0D0,35.0D0,56.0D0,84.0D0, . 120.0D0,4*0.0D0,1.0D0,5.0D0,15.0D0,35.0D0,70.0D0, . 126.0D0,210.0D0,5*0.0D0,1.0D0,6.0D0,21.0D0,56.0D0, . 126.0D0,252.0D0,6*0.0D0,1.0D0,7.0D0,28.0D0,84.0D0, . 210.0D0,7*0.0D0,1.0D0,8.0D0,36.0D0,120.0D0,8*0.0D0, . 1.0D0,9.0D0,45.0D0,9*0.0D0,1.0D0,10.0D0,10*0.0D0, . 1.0D0/ C C INITIALIZE (-1)**I FOR I=-15,-14,...,0,...15. C DATA NEG1 /-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1, . -1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1/ C DIMENSION SQRTL(0:7,0:7),SQRTF(0:7,0:7) LOGICAL FIRST DATA FIRST /.TRUE./ SAVE FACT,CTERM,SQRTL,SQRTF,BINOM,FIRST IF(FIRST)THEN FIRST = .FALSE. DO 20 IA=0,7 DO 10 IB=0,7 TERML = (2*IA+1)*(2*IB+1) SQRTL(IA,IB) = DSQRT(TERML) TERMF = FACT(2*IA)*FACT(2*IB) SQRTF(IA,IB) = DSQRT(TERMF) 10 CONTINUE 20 CONTINUE ENDIF C C IF THE CALLING ROUTINE HAS SUPPLIED A NEGATIVE MAGNETIC QUANTUM C NUMBER, TEMPORARILY CHANGE IT TO THE EQUIVALENT POSITIVE NUMBER. C IF(NCENTR.NE.1.AND.NCENTR.NE.2) NCENTR = 2 MSAVE = M IF(M.LT.0) M = -M R = RAB C C SOME FREQUENTLY USED QUANTITIES: C IF(NEWEXP)THEN RPLUS = R*(XIA + XIB)*0.5D0 RMINUS = R*(XIA - XIB)*0.5D0 ENDIF C C IF THE ORBITALS ARE ON DIFFERENT CENTERS, COMPUTE THE OVERLAP C ACCORDING TO J.P.P. STEWART'S FORMULA. NOTE THAT TERM1, TERM2, C TERM3, TERM4, AND TERM5 CORRESPOND, IN ORDER, TO THE FIRST 5 C TERMS IN THE INNERMOST SUM ON P. 28 OF THE REFERENCE. THE TERM C (-1)**IEXP HAS BEEN CORRECTED AND DIFFERS FROM STEWART'S FORMULA. C ATERM AND BTERM CORRESPOND TO THE A AND B INTEGRALS ON P. 28. C IF(NCENTR.EQ.2)THEN IF(NA.EQ.1.AND.NB.EQ.1.AND.ABS(RMINUS).LT.1.0D-6)THEN C C USE SIMPLE FORMULA FOR HYDROGEN-HYDROGEN OVERLAPS. C EXPLUS = EXP(-RPLUS) DEXPLUS = -EXP(-RPLUS)*(XIA + XIB)*0.5D0 SUM1 = 1.0D0/(3.0D0*RPLUS) DSUM1 = -1.0D0/(3.0D0*RPLUS**2)*(XIA + XIB)*0.5D0 SUM2 = (RPLUS + 1.0D0)/RPLUS**3 DSUM2 = -(XIA + XIB)*0.5D0*(2.0d0/RPLUS**3+3.0d0/RPLUS**4) SAB = EXPLUS*(SUM1 + SUM2)*(XIA**3)*RNAB DSAB = DEXPLUS*(SUM1 + SUM2)*(XIA**3)*RNAB . + EXPLUS*(DSUM1 + DSUM2)*(XIA**3)*RNAB . + EXPLUS*(SUM1 + SUM2)*(XIA**3)*(NA+NB+1)*RNAB/R ELSE C C AT LEAST ONE NON-HYDROGEN INVOLVED. C C C ASSIGN A-TYPE AND B-TYPE INTEGRALS IF XIA OR XIB HAS CHANGED. C NAB = NA + NB IF(NEWEXP)THEN CALL GABINTG(NAB,RPLUS,RMINUS,ATERM,BTERM,DATERM,DBTERM) DO 120 I=0,NAB DO 110 J=0,NAB ABTERM(I,J) = ATERM(I)*BTERM(J) DABTERM(I,J) = DATERM(I)*BTERM(J)*(XIA + XIB)*0.5D0 . + ATERM(I)*DBTERM(J)*(XIA - XIB)*0.5D0 110 CONTINUE 120 CONTINUE ENDIF C C COMPUTE SUMS. C IF(LA.EQ.0.AND.LB.EQ.0)THEN C C S(A)-S(B) OVERLAP. C SUMJA = 0.0D0 DSUMJA = 0.0D0 DO 140 IPA=0,NA T45 = BINOM(NA,IPA) DO 130 IPB=0,NB T345 = BINOM(NB,IPB)*T45 NASUB = IPA + IPB NBSUB = NAB - NASUB IEXP = NB - IPB AB = ABTERM(NASUB,NBSUB) DAB = DABTERM(NASUB,NBSUB) SUMJA = SUMJA + T345*AB*NEG1(IEXP) DSUMJA = DSUMJA + T345*DAB*NEG1(IEXP) 130 CONTINUE 140 CONTINUE ELSEIF(LA.EQ.1.AND.LB.EQ.0)THEN C C P(A)-S(B) OVERLAP. C SUMJA = 0.0D0 DSUMJA = 0.0D0 NBSUB0 = NAB - 1 IEXP0 = NB DO 170 IPA=0,NA-1 T45 = BINOM(NA-1,IPA) NASUB5 = IPA NBSUB5 = NBSUB0 - IPA DO 160 IPB=0,NB T345 = BINOM(NB,IPB)*T45 NASUB6 = NASUB5 + IPB NBSUB6 = NBSUB5 - IPB IEXP6 = IEXP0 - IPB SUMQA = 0.0D0 DSUMQA = 0.0D0 DO 150 IQA=0,1 NASUB = NASUB6 + IQA NBSUB = NBSUB6 + IQA IEXP = IEXP6 AB = ABTERM(NASUB,NBSUB) DAB = DABTERM(NASUB,NBSUB) SUMQA = SUMQA + AB*NEG1(IEXP) DSUMQA = DSUMQA + DAB*NEG1(IEXP) 150 CONTINUE SUMJA = SUMJA + T345*SUMQA DSUMJA = DSUMJA + T345*DSUMQA 160 CONTINUE 170 CONTINUE ELSEIF(LA.EQ.0.AND.LB.EQ.1)THEN C C S(A)-P(B) OVERLAP. C SUMJA = 0.0D0 DSUMJA = 0.0D0 NBSUB0 = NAB - 1 IEXP0 = NB - 1 DO 200 IPA=0,NA T45 = BINOM(NA,IPA) NASUB5 = IPA NBSUB5 = NBSUB0 - IPA DO 190 IPB=0,NB-1 T345 = BINOM(NB-1,IPB)*T45 NASUB6 = NASUB5 + IPB NBSUB6 = NBSUB5 - IPB IEXP6 = IEXP0 - IPB SUMQB = 0.0D0 DSUMQB = 0.0D0 DO 180 IQB=0,1 NASUB = NASUB6 + IQB NBSUB = NBSUB6 + IQB IEXP = IEXP6 + IQB AB = ABTERM(NASUB,NBSUB) DAB = DABTERM(NASUB,NBSUB) SUMQB = SUMQB + AB*NEG1(IEXP) DSUMQB = DSUMQB + DAB*NEG1(IEXP) 180 CONTINUE SUMJA = SUMJA + SUMQB*T345 DSUMJA = DSUMJA + DSUMQB*T345 190 CONTINUE 200 CONTINUE ELSEIF(LA.EQ.1.AND.LB.EQ.1)THEN C C P(A)-P(B) OVERLAP. C NBSUB0 = NAB - 2 IEXP0 = M + NB - 1 SUMKA = 0.0D0 DSUMKA = 0.0D0 DO 260 KA=0,M NASUB3 = 2*KA IEXP3 = IEXP0 - KA DO 250 KB=0,M NBSUB4 = NBSUB0 + 2*KB IEXP4 = IEXP3 + KB DO 240 IPA=0,NA-1 T45 = BINOM(NA-1,IPA) NASUB5 = NASUB3 + IPA NBSUB5 = NBSUB4 - IPA DO 230 IPB=0,NB-1 T345 = BINOM(NB-1,IPB)*T45 NASUB6 = NASUB5 + IPB NBSUB6 = NBSUB5 - IPB IEXP6 = IEXP4 - IPB SUMQA = 0.0D0 DSUMQA = 0.0D0 DO 220 IQA=0,LA-M NASUB7 = NASUB6 + IQA NBSUB7 = NBSUB6 + IQA DO 210 IQB=0,LB-M NASUB = NASUB7 + IQB NBSUB = NBSUB7 + IQB IEXP = IEXP6 + IQB AB = ABTERM(NASUB,NBSUB) DAB = DABTERM(NASUB,NBSUB) SUMQA = SUMQA + AB*NEG1(IEXP) DSUMQA = DSUMQA + DAB*NEG1(IEXP) 210 CONTINUE 220 CONTINUE SUMKA = SUMKA + SUMQA*T345 DSUMKA = DSUMKA + DSUMQA*T345 230 CONTINUE 240 CONTINUE 250 CONTINUE 260 CONTINUE SUMJA = SUMKA*CTERM(1,M,0)**2 DSUMJA = DSUMKA*CTERM(1,M,0)**2 ELSE C C USE GENERAL FORMULA FROM STEWART. C JAMAX = (LA-M)/2 JBMAX = (LB-M)/2 SUMJA = 0.0D0 DSUMJA = 0.0D0 NBSUB0 = NA + NB - LA - LB IEXP0 = M + NB - LB IPAM0 = NA - LA IPBM0 = NB - LB IQAM0 = LA - M IQBM0 = LB - M DO 440 JA=0,JAMAX JA2 = 2*JA IPAMAX = IPAM0 + JA2 IQAMAX = IQAM0 - JA2 NBSUB1 = NBSUB0 + JA2 SUMJB = 0.0D0 DSUMJB = 0.0D0 DO 420 JB=0,JBMAX JB2 = 2*JB IPBMAX = IPBM0 + JB2 IQBMAX = IQBM0 - JB2 NBSUB2 = NBSUB1 + JB2 IEXP2 = IEXP0 + JB2 SUMKA = 0.0D0 DSUMKA = 0.0D0 DO 400 KA=0,M NASUB3 = 2*KA IEXP3 = IEXP2 - KA DO 380 KB=0,M TERM5 = BINOM(M,KA)*BINOM(M,KB) NBSUB4 = NBSUB2 + 2*KB IEXP4 = IEXP3 + KB DO 360 IPA=0,IPAMAX TERM4 = BINOM(IPAMAX,IPA) T45 = TERM4*TERM5 NASUB5 = NASUB3 + IPA NBSUB5 = NBSUB4 - IPA DO 340 IPB=0,IPBMAX TERM3 = BINOM(IPBMAX,IPB) T345 = TERM3*T45 NASUB6 = NASUB5 + IPB NBSUB6 = NBSUB5 - IPB IEXP6 = IEXP4 - IPB DO 320 IQA=0,IQAMAX TERM2 = BINOM(IQAMAX,IQA) T2345 = TERM2*T345 NASUB7 = NASUB6 + IQA NBSUB7 = NBSUB6 + IQA DO 300 IQB=0,IQBMAX TERM1 = BINOM(IQBMAX,IQB) T12345 = TERM1*T2345 NASUB = NASUB7 + IQB NBSUB = NBSUB7 + IQB IEXP = IEXP6 + IQB C C NOTE OVERALL QUANTITIES: C C NASUB = 2*KA + IPA + IPB + IQA + IQB C C NBSUB = 2*KB + NA - LA + 2*JA + NB - LB C + 2*JB - IPA - IPB + IQA + IQB C C IEXP = M - KA + KB + NB - LB + 2*JB - IPB + IQB C AB = ABTERM(NASUB,NBSUB) DAB = DABTERM(NASUB,NBSUB) SUMKA = SUMKA + T12345*AB*NEG1(IEXP) DSUMKA = DSUMKA + T12345*DAB*NEG1(IEXP) 300 CONTINUE 320 CONTINUE 340 CONTINUE 360 CONTINUE 380 CONTINUE 400 CONTINUE SUMJB = SUMJB + SUMKA*CTERM(LB,M,JB) DSUMJB = DSUMJB + DSUMKA*CTERM(LB,M,JB) 420 CONTINUE SUMJA = SUMJA + SUMJB*CTERM(LA,M,JA) DSUMJA = DSUMJA + DSUMJB*CTERM(LA,M,JA) 440 CONTINUE ENDIF TERMXI = EXHLA*EXHLB TERMNL = SQRTL(LA,LB)/SQRTF(NA,NB) SAB = TERMXI*TERMNL*RNAB*SUMJA DSAB = TERMXI*TERMNL*RNAB*SUMJA*(NA+NB+1)/R . + TERMXI*TERMNL*RNAB*DSUMJA ENDIF ELSE C C THE ORBITALS ARE CENTERED ON THE SAME ATOM. SINCE STO'S DO C NOT HAVE RADIAL NODES, THE RADIAL FUNCTIONS ARE NOT ORTHOGONAL, C AND THEREFORE THE ORBITALS MAY NOT BE ORTHOGONAL. C IF(LA.NE.LB)THEN C C ORBITALS ARE ORTHOGONAL DUE TO THEIR ORTHOGONAL SPHERICAL C HARMONICS. C SAB = 0.0D0 DSAB = 0.0D0 ELSE C C ORBITALS HAVE THE SAME TOTAL ANGULAR MOMENTUM. COMPUTE THE C OVERLAP BASED ON AN INTEGRATION OVER THE RADIAL PORTIONS OF C THE STO'S. C RADIAL = FACT(NA+NB)/((XIA+XIB)**(NA+NB+1)) TERMXI = (2**(NA+1))*EXHLA*(2**(NB+1))*EXHLB C C TAKE INTO ACCOUNT THE FACT THAT THE SPHERICAL HARMONICS COULD C BE 180 DEGREES OUT OF PHASE SINCE THE COORDINATE SYSTEMS ARE C MIRROR IMAGES. C SAB = NEG1(LA+M)*RADIAL*TERMXI/SQRTF(NA,NB) DSAB = 0.0d0 ENDIF ENDIF C C RETURN THE MAGNETIC QUANTUM NUMBER BACK TO ITS ORIGINAL VALUE C IF IT WAS CHANGED. C M = MSAVE RETURN END C C C SUBROUTINE GABINTG(N,XA,XB,ATERM,BTERM,DATERM,DBTERM) C C COMPUTES A-TYPE AND B-TYPE INTEGRALS. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) real*8 ATERM(0:15),BTERM(0:15) real*8 DATERM(0:15),DBTERM(0:15) C C LOCAL VARIABLES: C real*8 FACT(0:15) C C INITIALIZE FACTORIAL VALUES (UP TO 15!). C DATA FACT /1.0D0,1.0D0,2.0D0,6.0D0,24.0D0,120.0D0,720.0D0, . 5040.0D0,40320.0D0,362880.0D0,36288.0D2,399168.0D2, . 4790016.0D2,62270208.0D2,871782912.0D2,1307674368.0D3/ SAVE FACT C C ON FIRST CALL, INITIALIZE SOME QUANTITIES USED REPETITIVELY IN THE C B-TYPE INTEGRALS. C DIMENSION XBEXP(15),BFACTR(0:7,0:15),BTERM0(0:15),NEG1(0:16) DIMENSION DXBEXP(15) LOGICAL FIRST DATA FIRST /.TRUE./ C C INITIALIZE (-1)**I FOR I=0,1,...,15 C DATA NEG1 /1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1/ C SAVE FIRST,BFACTR,BTERM0,NEG1 IF(FIRST)THEN FIRST = .FALSE. DO 20 I=0,7 DO 10 J=0,15 IJ1 = I+J+1 BFACTR(I,J) = 2*MOD(IJ1,2)/(FACT(J)*IJ1) 10 CONTINUE 20 CONTINUE DO 30 I=0,15 BTERM0(I) = (1.0D0 - DFLOAT(NEG1(I+1)))/DFLOAT(I+1) 30 CONTINUE ENDIF C XAINV = 1.0D0/XA ATERM(0) = XAINV*EXP(-XA) DATERM(0) = -XAINV*EXP(-XA) - EXP(-XA)*XAINV**2 DO 130 I=1,N ATERM(I) = ATERM(0) + I*ATERM(I-1)*XAINV DATERM(I) = DATERM(0) + I*DATERM(I-1)*XAINV . - I*ATERM(I-1)*XAINV**2 130 CONTINUE C C THE DIRECT COMPUTATION OF BTERM CAN BE UNSTABLE, DEPENDING UPON C THE SIZE OF XB. IF ABS(XB) IS LARGE OR VERY SMALL, THEN THERE C IS NO PROBLEM. HOWEVER, FOR INTERMEDIATE SIZED VALUES OF ABS(XB), C A TAYLOR EXPANSION FOR THE INTEGRAL MUST BE USED. C AXB = ABS(XB) IF(AXB.LT.1.0D-6)THEN DO 140 I=0,N BTERM(I) = BTERM0(I) DBTERM(I) = 0.0d0 140 CONTINUE ELSEIF(AXB.GT.3.0D0)THEN XBINV = 1.0D0/XB DXBINV = -1.0D0/XB**2 EXB = EXP(XB) DEXB = EXP(XB) EXBINV = 1.0D0/EXB DEXBINV = -EXBINV BTERM(0) = XBINV*(EXB - EXBINV) DBTERM(0) = DXBINV*(EXB - EXBINV) + XBINV*(DEXB - DEXBINV) DO 160 I=1,N BTERM(I) = XBINV*(I*BTERM(I-1) + NEG1(I)*EXB - EXBINV) DBTERM(I) = DXBINV*(I*BTERM(I-1) + NEG1(I)*EXB - EXBINV) . + XBINV*(I*DBTERM(I-1) + NEG1(I)*DEXB - DEXBINV) 160 CONTINUE ELSE IF(AXB.GT.2.0D0)THEN JLAST = 15 ELSEIF(AXB.GT.1.0D0)THEN JLAST = 12 ELSEIF(AXB.GT.0.5D0)THEN JLAST = 7 ELSE JLAST = 6 ENDIF XBEXP(1) = -XB DXBEXP(1) = -1.0d0 DO 170 J=2,JLAST XBEXP(J) = -XB*XBEXP(J-1) DXBEXP(J) = -XB*DXBEXP(J-1) - XBEXP(J-1) 170 CONTINUE DO 200 I=0,N BI = BFACTR(I,0) DBI = 0.0d0 DO 180 J=1,JLAST BI = BI + XBEXP(J)*BFACTR(I,J) DBI = DBI + DXBEXP(J)*BFACTR(I,J) 180 CONTINUE BTERM(I) = BI DBTERM(I) = DBI 200 CONTINUE ENDIF RETURN END C C C SUBROUTINE GSLIMIT(NIJSC,ISC,JSC,NKIJSC,KIJSC,NKC,KC) C C SETS UP ARRAY LIMITS TO TAKE ADVANTAGE OF SPARSITY OF MATRICES C USED IN DIATOMIC OVERLAP CALCULATION. THE ULTIMATE GOAL IS C TO COMPUTE C'*SABLOC*C USING THE MINIMUM NUMBER OF FLOATING C POINT OPERATIONS. THE DEGREE OF SPARSITY DEPENDS DEPENDS ON C HOW MANY ATOMIC ORBITALS ARE ON EACH ATOM. FOR NORBSA.GE.NORBSB, C THE COMBINATIONS ARE CODED ACCORDING TO THE VARIABLE NAB: C C NORBSA NORBSB NAB C 1 1 1 C 4 1 2 C 4 4 3 C 9 1 4 C 9 4 5 C 9 9 6 C C ONCE THE VALUE OF NAB IS ESTABLISHED, THE SPARSITY VARIABLES ARE C DEFINED: C C NIJSC(NAB) = THE NUMBER OF NONZERO ELEMENTS IN SABLOC*C. C C ************************************************************************* C * * C * FOR IJ=1,NIJSC(NAB) * C * * C * ISC(NAB,IJ) = ROW NUMBER OF IJTH NONZERO ELEMENT OF SABLOC*C * C * * C * JSC(NAB,IJ) = COLUMN NUMBER OF IJTH NONZERO ELEMENT OF SABLOC*C * C * * C * NKIJSC(NAB,IJ) = NUMBER MULTIPLICATIONS REQUIRED TO COMPUTE THE * C * (ISC(NAB,IJ),JSC(NAB,IJ)) ENTRY OF SABLOC*C * C * * C * KIJSC(NAB,IJ,K) = INNER INDEX FOR ELEMENTS REQUIRED TO COMPUTE * C * THE (ISC(NAB,IJ),JSC(NAB,IJ)) ENTRY OF * C * SABLOC*C, I.E., * C * * C * I=ISC(NAB,IJ) * C * J=JSC(NAB,IJ) * C * SCIJ = 0.0 * C * DO K=1,NKIJSC(NAB,IJ) * C * KINDX = KIJSC(NAB,IJ,K) * C * SCIJ = SCIJ + SABLOC(I,KINDX)*C(KINDX,J) * C * CONTINUE * C * (I,J) ELMENT OF SABLOC*C = SCIJ * C * * C ************************************************************************* C C NKC(NAB,J) = THE NUMBER OF NONZERO ELEMENTS IN COLUMN J OF C. C (KC(NAB,J,K), K=1,NKC(NAB,J)) = LIST OF ROW NUMBERS IN COLUMN J C OF C THAT HAVE NONZERO ENTRIES. C integer NIJSC(6),ISC(6,81),JSC(6,81),NKIJSC(6,81),KIJSC(6,81,9), . NKC(6,9),KC(6,9,9) C C LOCAL ARRAYS: C integer IS(9,9),IC(9,9),NORBS(3) C C FLAGS FOR NON-ZERO ENTRIES OF LOCAL FRAME OVERLAP MATRIX SABLOC: C DATA IS /1,0,0,1,1,0,0,0,0, . 0,1,0,0,0,1,0,0,0, . 0,0,1,0,0,0,1,0,0, . 1,0,0,1,1,0,0,0,0, . 1,0,0,1,1,0,0,0,0, . 0,1,0,0,0,1,0,0,0, . 0,0,1,0,0,0,1,0,0, . 0,0,0,0,0,0,0,1,0, . 0,0,0,0,0,0,0,0,1/ C C FLAGS FOR NON-ZERO ENTRIES OF TRANSFORMATION MATRIX C: C DATA IC /1,0,0,0,0,0,0,0,0, . 0,1,1,1,0,0,0,0,0, . 0,1,1,1,0,0,0,0,0, . 0,1,1,1,0,0,0,0,0, . 0,0,0,0,1,1,1,1,1, . 0,0,0,0,1,1,1,1,1, . 0,0,0,0,1,1,1,1,1, . 0,0,0,0,1,1,1,1,1, . 0,0,0,0,1,1,1,1,1/ C C POSSIBLE NUMBERS OF ORBITALS ON AN ATOM: C DATA NORBS /1,4,9/ C SAVE IS,IC,NORBS C C DETERMINE WHICH ENTRIES OF SABLOC*C NEED TO BE COMPUTED BASED ON C THE NON-ZERO ENTRIES OF SABLOC AND C. CONSIDER EACH POSSIBLE C COMBINATION OF NORBSA AND NORBSB WITH NORBSA.GE.NORBSB. C IAB = 0 DO 500 NA=1,3 NORBSA = NORBS(NA) DO 480 NB=1,NA IAB = IAB + 1 IJ = 0 NORBSB = NORBS(NB) DO 160 I=1,NORBSA NKC(IAB,I) = 0 DO 140 J=1,NORBSA IF(IC(J,I).EQ.1)THEN NKC(IAB,I) = NKC(IAB,I) + 1 KC(IAB,I,NKC(IAB,I)) = J ENDIF IJFLAG = 0 DO 120 K=1,NORBSB IF(IS(I,K).NE.0.AND.IC(K,J).NE.0)THEN C C (I,J) ENTRY OF C*SABLOC IS NONZERO. C IF(IJFLAG.EQ.0)THEN IJFLAG = 1 IJ = IJ + 1 NKIJSC(IAB,IJ) = 0 ENDIF NKIJSC(IAB,IJ) = NKIJSC(IAB,IJ) + 1 KIJSC(IAB,IJ,NKIJSC(IAB,IJ)) = K ENDIF 120 CONTINUE IF(IJFLAG.EQ.1)THEN ISC(IAB,IJ) = I JSC(IAB,IJ) = J ENDIF 140 CONTINUE 160 CONTINUE NIJSC(IAB) = IJ 480 CONTINUE 500 CONTINUE RETURN END C C C SUBROUTINE GRTRANS(XA,YA,ZA,XB,YB,ZB,RAB,NORBS, . C,DXC,DYC,DZC,iswap) C C COMPUTES A MATRIX C THAT TRANSFORMS LOCAL FRAME OVERLAPS SABLOC C TO THE MOLECULAR FRAME OVERLAPS SABDIA: C C SABDIA = C'*SABLOC*C. C C C WILL BE EITHER 4 BY 4 OR 9 BY 9, DEPENDING ON WHETHER NORBS=4 C OR NORBS=9. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) real*8 C(9,9),DXC(9,9),DYC(9,9),DZC(9,9) C C LOCAL ARRAYS: C real*8 ROT(3,3),D(5,3,3),DXD(5,3,3),DYD(5,3,3),DZD(5,3,3) real*8 DXROT(3,3),DYROT(3,3),DZROT(3,3) C C SQRT(1/2) AND SQRT(1/6): C DATA SQ2,SQ6 /0.707106769084930420D0, 0.408248305320739746D0/ SAVE SQ2,SQ6 C C COEFFICIENT FOR S-TYPE ORBITALS NOT EFFECTED BY DIRECTIONALITY: C C(1,1) = 1.0D0 DXC(1,1) = 0.0D0 DYC(1,1) = 0.0D0 DZC(1,1) = 0.0D0 C C COMPUTE A ROTATION MATRIX ROT THAT RELATES THE X, Y, AND Z AXES C OF ATOMS A AND B TO A COORDINATE SYSTEM (X*,Y*,Z*), WHERE Z* POINTS C FROM ATOM A TO ATOM B. THE DIRECTIONS OF X* AND Y* CAN BE DEFINED C ARBITRARILY, AS LONG AS THE SYSTEM (X*,Y*,Z*) IS RIGHT-HANDED. C CALL GRMATRX(XA,YA,ZA,XB,YB,ZB,RAB,ROT,DXROT,DYROT,DZROT,iswap) C C NOW ASSIGN COEFFICIENTS FOR P-TYPE ORBITALS USING THE FACT THAT C THE PX, PY, AND PZ ORBITALS TRANSFORM JUST LIKE THE X, Y, AND Z C AXES. C DO 60 I=1,3 DO 40 J=1,3 C(I+1,J+1) = ROT(I,J) DXC(I+1,J+1) = DXROT(I,J) DYC(I+1,J+1) = DYROT(I,J) DZC(I+1,J+1) = DZROT(I,J) 40 CONTINUE 60 CONTINUE C C IF ATOMS A AND B HAVE NO D ORBITALS WE ARE DONE. C IF(NORBS.LE.4) GO TO 1000 C C NOW FORM FIVE SYMMETRIC TENSORS TO REPRESENT THE D(Z**2), D(XZ), C D(YZ), D(X**2-Y**2), AND D(XY) ORBITALS IN THE (X*,Y*,Z*) COORDINATE C SYSTEM. THIS WILL INVOLVE TAKING DYAD PRODUCTS OF COLUMN VECTORS. C FOR COLUMN VECTORS X AND Y, THE DYAD PRODUCT IS DEFINED AS: C C X DYAD Y = X*Y' = OUTER PRODUCT = TENSOR MATRIX. C DO 100 I=1,3 DO 80 J=1,3 C C D(Z**2) = (3*(Z DYAD Z) - I)/SQRT(6) (I = IDENTITY TENSOR): C D(1,I,J) = ROT(I,3)*ROT(J,3)*3.0D0*SQ6 DXD(1,I,J) = DXROT(I,3)*ROT(J,3)*3.0D0*SQ6 . + ROT(I,3)*DXROT(J,3)*3.0D0*SQ6 DYD(1,I,J) = DYROT(I,3)*ROT(J,3)*3.0D0*SQ6 . + ROT(I,3)*DYROT(J,3)*3.0D0*SQ6 DZD(1,I,J) = DZROT(I,3)*ROT(J,3)*3.0D0*SQ6 . + ROT(I,3)*DZROT(J,3)*3.0D0*SQ6 IF(I.EQ.J) then D(1,I,J) = D(1,I,J) - SQ6 DXD(1,I,J) = DXD(1,I,J) DYD(1,I,J) = DYD(1,I,J) DZD(1,I,J) = DZD(1,I,J) endif C C D(XZ) = (X DYAD Z + Z DYAD X)/SQRT(2): C D(2,I,J) = (ROT(I,1)*ROT(J,3) + ROT(I,3)*ROT(J,1))*SQ2 DXD(2,I,J) = (DXROT(I,1)*ROT(J,3) + DXROT(I,3)*ROT(J,1))*SQ2 . + (ROT(I,1)*DXROT(J,3) + ROT(I,3)*DXROT(J,1))*SQ2 DYD(2,I,J) = (DYROT(I,1)*ROT(J,3) + DYROT(I,3)*ROT(J,1))*SQ2 . + (ROT(I,1)*DYROT(J,3) + ROT(I,3)*DYROT(J,1))*SQ2 DZD(2,I,J) = (DZROT(I,1)*ROT(J,3) + DZROT(I,3)*ROT(J,1))*SQ2 . + (ROT(I,1)*DZROT(J,3) + ROT(I,3)*DZROT(J,1))*SQ2 C C D(YZ) = (Y DYAD Z + Z DYAD Y)/SQRT(2): C D(3,I,J) = (ROT(I,2)*ROT(J,3) + ROT(I,3)*ROT(J,2))*SQ2 DXD(3,I,J) = (DXROT(I,2)*ROT(J,3) + DXROT(I,3)*ROT(J,2))*SQ2 . + (ROT(I,2)*DXROT(J,3) + ROT(I,3)*DXROT(J,2))*SQ2 DYD(3,I,J) = (DYROT(I,2)*ROT(J,3) + DYROT(I,3)*ROT(J,2))*SQ2 . + (ROT(I,2)*DYROT(J,3) + ROT(I,3)*DYROT(J,2))*SQ2 DZD(3,I,J) = (DZROT(I,2)*ROT(J,3) + DZROT(I,3)*ROT(J,2))*SQ2 . + (ROT(I,2)*DZROT(J,3) + ROT(I,3)*DZROT(J,2))*SQ2 C C D(X**2-Y**2) = (X DYAD X - Y DYAD Y)/SQRT(2): C D(4,I,J) = (ROT(I,1)*ROT(J,1) - ROT(I,2)*ROT(J,2))*SQ2 DXD(4,I,J) = (DXROT(I,1)*ROT(J,1) - DXROT(I,2)*ROT(J,2))*SQ2 . + (ROT(I,1)*DXROT(J,1) - ROT(I,2)*DXROT(J,2))*SQ2 DYD(4,I,J) = (DYROT(I,1)*ROT(J,1) - DYROT(I,2)*ROT(J,2))*SQ2 . + (ROT(I,1)*DYROT(J,1) - ROT(I,2)*DYROT(J,2))*SQ2 DZD(4,I,J) = (DZROT(I,1)*ROT(J,1) - DZROT(I,2)*ROT(J,2))*SQ2 . + (ROT(I,1)*DZROT(J,1) - ROT(I,2)*DZROT(J,2))*SQ2 C C D(XY) = (X DYAD Y + Y DYAD X)/SQRT(2): C D(5,I,J) = (ROT(I,1)*ROT(J,2) + ROT(I,2)*ROT(J,1))*SQ2 DXD(5,I,J) = (DXROT(I,1)*ROT(J,2) + DXROT(I,2)*ROT(J,1))*SQ2 . + (ROT(I,1)*DXROT(J,2) + ROT(I,2)*DXROT(J,1))*SQ2 DYD(5,I,J) = (DYROT(I,1)*ROT(J,2) + DYROT(I,2)*ROT(J,1))*SQ2 . + (ROT(I,1)*DYROT(J,2) + ROT(I,2)*DYROT(J,1))*SQ2 DZD(5,I,J) = (DZROT(I,1)*ROT(J,2) + DZROT(I,2)*ROT(J,1))*SQ2 . + (ROT(I,1)*DZROT(J,2) + ROT(I,2)*DZROT(J,1))*SQ2 80 CONTINUE 100 CONTINUE C C NOW GET THE PROJECTIONS OF THE ABOVE TENSORS ON THE CORRESPONDING C TENSORS FORMED BY USING THE STANDARD BASIS: C C 1 0 0 C X* = 0 Y* = 1 Z* = 0 C 0 0 1 C C THIS REQUIRES A SCALAR PRODUCT OF TWO TENSORS. FOR TENSORS A AND C B, THE SCALAR PRODUCT (A|B) IS DEFINED AS: C C (A|B) = (B|A) = SUM(I=1,3;J=1,3: A(I,J)*B(I,J)) C DO 120 K=1,5 K4 = K+4 C(5,K4) = -D(K,1,1)*SQ6 - D(K,2,2)*SQ6 + D(K,3,3)*2.0*SQ6 DXC(5,K4) = -DXD(K,1,1)*SQ6 - DXD(K,2,2)*SQ6 + DXD(K,3,3)*2.0*SQ6 DYC(5,K4) = -DYD(K,1,1)*SQ6 - DYD(K,2,2)*SQ6 + DYD(K,3,3)*2.0*SQ6 DZC(5,K4) = -DZD(K,1,1)*SQ6 - DZD(K,2,2)*SQ6 + DZD(K,3,3)*2.0*SQ6 C(6,K4) = D(K,1,3)*SQ2 + D(K,3,1)*SQ2 DXC(6,K4) = DXD(K,1,3)*SQ2 + DXD(K,3,1)*SQ2 DYC(6,K4) = DYD(K,1,3)*SQ2 + DYD(K,3,1)*SQ2 DZC(6,K4) = DZD(K,1,3)*SQ2 + DZD(K,3,1)*SQ2 C(7,K4) = D(K,2,3)*SQ2 + D(K,3,2)*SQ2 DXC(7,K4) = DXD(K,2,3)*SQ2 + DXD(K,3,2)*SQ2 DYC(7,K4) = DYD(K,2,3)*SQ2 + DYD(K,3,2)*SQ2 DZC(7,K4) = DZD(K,2,3)*SQ2 + DZD(K,3,2)*SQ2 C(8,K4) = D(K,1,1)*SQ2 - D(K,2,2)*SQ2 DXC(8,K4) = DXD(K,1,1)*SQ2 - DXD(K,2,2)*SQ2 DYC(8,K4) = DYD(K,1,1)*SQ2 - DYD(K,2,2)*SQ2 DZC(8,K4) = DZD(K,1,1)*SQ2 - DZD(K,2,2)*SQ2 C(9,K4) = D(K,1,2)*SQ2 + D(K,2,1)*SQ2 DXC(9,K4) = DXD(K,1,2)*SQ2 + DXD(K,2,1)*SQ2 DYC(9,K4) = DYD(K,1,2)*SQ2 + DYD(K,2,1)*SQ2 DZC(9,K4) = DZD(K,1,2)*SQ2 + DZD(K,2,1)*SQ2 120 CONTINUE 1000 RETURN END C C C SUBROUTINE GRMATRX(XA,YA,ZA,XB,YB,ZB,RAB,ROT,DXROT,DYROT,DZROT & ,iswap) C C CONSTRUCTS A ROTATION MATRIX ROT THAT ALIGNS THE Z AXIS WITH THE C THE VECTOR (XB-XA,YB-YA,ZB-ZA). C IMPLICIT DOUBLE PRECISION (A-H,O-Z) real*8 ROT(3,3),X(3),Y(3),Z(3) real*8 DXROT(3,3),DYROT(3,3),DZROT(3,3) real*8 DXX(3),DXY(3),DXZ(3),DYX(3),DYY(3), & DYZ(3),DZX(3),DZY(3),DZZ(3) Z(1)= -(XA - XB) Z(2)= -(YA - YB) Z(3)= -(ZA - ZB) ZINV = 1/RAB if (iswap.ne.1) then DXZINV = -(XA-XB)/RAB**3 DYZINV = -(YA-YB)/RAB**3 DZZINV = -(ZA-ZB)/RAB**3 Z(1) = Z(1)*ZINV DXZ(1)= -(XA - XB)*DXZINV - ZINV DYZ(1)= -(XA - XB)*DYZINV DZZ(1)= -(XA - XB)*DZZINV Z(2) = Z(2)*ZINV DXZ(2)= -(YA - YB)*DXZINV DYZ(2)= -(YA - YB)*DYZINV - ZINV DZZ(2)= -(YA - YB)*DZZINV Z(3) = Z(3)*ZINV DXZ(3)= -(ZA - ZB)*DXZINV DYZ(3)= -(ZA - ZB)*DYZINV DZZ(3)= -(ZA - ZB)*DZZINV - ZINV else DXZINV = (XA-XB)/RAB**3 DYZINV = (YA-YB)/RAB**3 DZZINV = (ZA-ZB)/RAB**3 Z(1) = Z(1)*ZINV DXZ(1)= -(XA - XB)*DXZINV + ZINV DYZ(1)= -(XA - XB)*DYZINV DZZ(1)= -(XA - XB)*DZZINV Z(2) = Z(2)*ZINV DXZ(2)= -(YA - YB)*DXZINV DYZ(2)= -(YA - YB)*DYZINV + ZINV DZZ(2)= -(YA - YB)*DZZINV Z(3) = Z(3)*ZINV DXZ(3)= -(ZA - ZB)*DXZINV DYZ(3)= -(ZA - ZB)*DYZINV DZZ(3)= -(ZA - ZB)*DZZINV + ZINV endif IF(ABS(Z(3)).GT.0.99999999D0)THEN C C ATOMS A AND B LIE ALONG THE Z AXIS. CONSTRUCT A ZERO OR 180 DEGREE C ROTATION ABOUT THE X AXIS. C Z(1) = 0.0D0 DXZ(1) = 0.0D0 DYZ(1) = 0.0D0 DZZ(1) = 0.0D0 Z(2) = 0.0D0 DXZ(2) = 0.0D0 DYZ(2) = 0.0D0 DZZ(2) = 0.0D0 Z(3) = SIGN(1.0D0,Z(3)) DXZ(3) = 0.0D0 DYZ(3) = 0.0D0 DZZ(3) = 0.0D0 X(1) = 1.0D0 DXX(1) = 0.0D0 DYX(1) = 0.0D0 DZX(1) = 0.0D0 X(2) = 0.0D0 DXX(2) = 0.0D0 DYX(2) = 0.0D0 DZX(2) = 0.0D0 X(3) = 0.0D0 DXX(3) = 0.0D0 DYX(3) = 0.0D0 DZX(3) = 0.0D0 Y(1) = 0.0D0 DXY(1) = 0.0D0 DYY(1) = 0.0D0 DZY(1) = 0.0D0 Y(2) = Z(3) DXY(2) = 0.0D0 DYY(2) = 0.0D0 DZY(2) = 0.0D0 Y(3) = 0.0D0 DXY(3) = 0.0D0 DYY(3) = 0.0D0 DZY(3) = 0.0D0 ELSE Z12SQR = Z(1)**2 + Z(2)**2 DXZ12SQR = 2.0d0*(Z(1)*DXZ(1)+Z(2)*DXZ(2)) DYZ12SQR = 2.0d0*(Z(1)*DYZ(1)+Z(2)*DYZ(2)) DZZ12SQR = 2.0d0*(Z(1)*DZZ(1)+Z(2)*DZZ(2)) Z12 = DSQRT(Z12SQR) DXZ12 = DXZ12SQR/(2.0d0*Z12) DYZ12 = DYZ12SQR/(2.0d0*Z12) DZZ12 = DZZ12SQR/(2.0d0*Z12) X(1) = -Z(2)/Z12 DXX(1) = -DXZ(2)/Z12 + Z(2)*DXZ12/Z12**2 DYX(1) = -DYZ(2)/Z12 + Z(2)*DYZ12/Z12**2 DZX(1) = -DZZ(2)/Z12 + Z(2)*DZZ12/Z12**2 X(2) = Z(1)/Z12 DXX(2) = DXZ(1)/Z12 - Z(1)*DXZ12/Z12**2 DYX(2) = DYZ(1)/Z12 - Z(1)*DYZ12/Z12**2 DZX(2) = DZZ(1)/Z12 - Z(1)*DZZ12/Z12**2 X(3) = 0.0D0 DXX(3) = 0.0D0 DYX(3) = 0.0D0 DZX(3) = 0.0D0 IF(ABS(Z(3)).LT.1.0D-8)THEN Y(1) = 0.0D0 DXY(1) = 0.0D0 DYY(1) = 0.0D0 DZY(1) = 0.0D0 Y(2) = 0.0D0 DXY(2) = 0.0D0 DYY(2) = 0.0D0 DZY(2) = 0.0D0 Y(3) = 1.0D0 DXY(3) = 0.0D0 DYY(3) = 0.0D0 DZY(3) = 0.0D0 ELSE YINV = 1.0D0/(Z12*DSQRT(1.0D0 + Z12SQR/Z(3)**2)) DXYINV = -(DXZ12*DSQRT(1.0D0 + Z12SQR/Z(3)**2) & +Z12*0.5d0*(-2.0d0*Z12SQR*DXZ(3)/Z(3)**3 & + DXZ12SQR/Z(3)**2)/DSQRT(1.0D0 + Z12SQR/Z(3)**2)) & /(Z12*DSQRT(1.0D0 + Z12SQR/Z(3)**2))**2 DYYINV = -(DYZ12*DSQRT(1.0D0 + Z12SQR/Z(3)**2) & +Z12*0.5d0*(-2.0d0*Z12SQR*DYZ(3)/Z(3)**3 & + DYZ12SQR/Z(3)**2)/DSQRT(1.0D0 + Z12SQR/Z(3)**2)) & /(Z12*DSQRT(1.0D0 + Z12SQR/Z(3)**2))**2 DZYINV = -(DZZ12*DSQRT(1.0D0 + Z12SQR/Z(3)**2) & +Z12*0.5d0*(-2.0d0*Z12SQR*DZZ(3)/Z(3)**3 & + DZZ12SQR/Z(3)**2)/DSQRT(1.0D0 + Z12SQR/Z(3)**2)) & /(Z12*DSQRT(1.0D0 + Z12SQR/Z(3)**2))**2 Y(1) = Z(1)*YINV DXY(1) = DXZ(1)*YINV + Z(1)*DXYINV DYY(1) = DYZ(1)*YINV + Z(1)*DYYINV DZY(1) = DZZ(1)*YINV + Z(1)*DZYINV Y(2) = Z(2)*YINV DXY(2) = DXZ(2)*YINV + Z(2)*DXYINV DYY(2) = DYZ(2)*YINV + Z(2)*DYYINV DZY(2) = DZZ(2)*YINV + Z(2)*DZYINV Y(3) = -(Z12SQR/Z(3))*YINV DXY(3) = -(DXZ12SQR/Z(3))*YINV & +(Z12SQR*DXZ(3)/Z(3)**2)*YINV & - (Z12SQR/Z(3))*DXYINV DYY(3) = -(DYZ12SQR/Z(3))*YINV & +(Z12SQR*DYZ(3)/Z(3)**2)*YINV & - (Z12SQR/Z(3))*DYYINV DZY(3) = -(DZZ12SQR/Z(3))*YINV & +(Z12SQR*DZZ(3)/Z(3)**2)*YINV & - (Z12SQR/Z(3))*DZYINV ENDIF ENDIF DO 10 J=1,3 ROT(1,J) = X(J) ROT(2,J) = Y(J) ROT(3,J) = Z(J) DXROT(1,J) = DXX(J) DXROT(2,J) = DXY(J) DXROT(3,J) = DXZ(J) DYROT(1,J) = DYX(J) DYROT(2,J) = DYY(J) DYROT(3,J) = DYZ(J) DZROT(1,J) = DZX(J) DZROT(2,J) = DZY(J) DZROT(3,J) = DZZ(J) 10 CONTINUE RETURN END