C C $Id: overlp.F,v 1.3 1998/07/16 16:40:10 jjv5 Exp arjan $ C C------------------------------------------------------------------------ SUBROUTINE OVERLP(IA,NORBSA,NPQA,XISA,XIPA,XIDA,XA,YA,ZA, . IB,NORBSB,NPQB,XISB,XIPB,XIDB,XB,YB,ZB, . RAB,SABDIA) 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" DIMENSION SABDIA(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 DIMENSION SABLOC(9,9),SCTEMP(9,9),C(9,9),ICODE(18) DIMENSION NIJSC(6),ISC(6,81),JSC(6,81),NKIJSC(6,81),KIJSC(6,81,9), . NKC(6,9),KC(6,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 SLIMIT(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 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. SABLOC(1,1) = SAB(NCENTR,NPQA,LA,XISA,EXHSA, . NPQB,LB,XISB,EXHSB, . NEWEXP,M,R,RNAB) C C : C IF(NORBSA.GT.1)THEN LA = 1 LB = 0 M = 0 NEWEXP = .TRUE. SABLOC(4,1) = SAB(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XISB,EXHSB, . NEWEXP,M,R,RNAB) ENDIF C C : C IF(NORBSA.GT.4)THEN LA = 2 LB = 0 M = 0 NEWEXP = .TRUE. SABLOC(5,1) = SAB(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XISB,EXHSB, . NEWEXP,M,R,RNAB) ENDIF C C , , AND : C IF(NORBSA.GT.4.AND.NORBSB.GT.1)THEN LA = 2 LB = 1 M = 1 NEWEXP = .TRUE. SABLOC(6,2) = SAB(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB) SABLOC(7,3) = SABLOC(6,2) M = 0 NEWEXP = .FALSE. SABLOC(5,4) = -SAB(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB) ENDIF C C : C IF(NORBSB.GT.1)THEN IF(IA.EQ.IB)THEN SABLOC(1,4) = -SABLOC(4,1) ELSE LA = 0 LB = 1 M = 0 NEWEXP = .TRUE. SABLOC(1,4) = -SAB(NCENTR,NPQA,LA,XISA,EXHSA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB) ENDIF ENDIF C C : C IF(NORBSB.GT.4)THEN IF(IA.EQ.IB)THEN SABLOC(1,5) = SABLOC(5,1) ELSE LA = 0 LB = 2 M = 0 NEWEXP = .TRUE. SABLOC(1,5) = SAB(NCENTR,NPQA,LA,XISA,EXHSA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB) 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) SABLOC(4,5) = -SABLOC(5,4) ELSE LA = 1 LB = 2 M = 1 NEWEXP = .TRUE. SABLOC(2,6) = -SAB(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB) M = 0 NEWEXP = .FALSE. SABLOC(4,5) = SAB(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB) ENDIF SABLOC(3,7) = SABLOC(2,6) ENDIF C C , , : C IF(NORBSA.GT.1.AND.NORBSB.GT.1)THEN LA = 1 LB = 1 M = 1 NEWEXP = .TRUE. SABLOC(2,2) = SAB(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB) NEWEXP = .FALSE. SABLOC(3,3) = SABLOC(2,2) M = 0 SABLOC(4,4) = -SAB(NCENTR,NPQA,LA,XIPA,EXHPA, . NPQB,LB,XIPB,EXHPB, . NEWEXP,M,R,RNAB) ENDIF C C THE FIVE OVERLAPS: C IF(NORBSA.GT.4.AND.NORBSB.GT.4)THEN LA = 2 LB = 2 M = 0 NEWEXP = .TRUE. SABLOC(5,5) = SAB(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB) M = 1 NEWEXP = .FALSE. SABLOC(6,6) = -SAB(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB) SABLOC(7,7) = SABLOC(6,6) M = 2 SABLOC(8,8) = SAB(NCENTR,NPQA,LA,XIDA,EXHDA, . NPQB,LB,XIDB,EXHDB, . NEWEXP,M,R,RNAB) SABLOC(9,9) = SABLOC(8,8) 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) ELSE CALL RTRANS(XA,YA,ZA,XB,YB,ZB,RAB,NORBS,C) 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 DO 160 K=1,NKIJSC(IAB,IJ) KINDX = KIJSC(IAB,IJ,K) TMPIJ = TMPIJ + SABLOC(I,KINDX)*C(KINDX,J) 160 CONTINUE SCTEMP(I,J) = TMPIJ 200 CONTINUE DO 260 I=1,NORBSA DO 240 J=1,NORBSB SIJ = 0.0D0 DO 220 K=1,NKC(IAB,I) KINDX = KC(IAB,I,K) SIJ = SIJ + C(KINDX,I)*SCTEMP(KINDX,J) 220 CONTINUE SABDIA(I,J) = SIJ 240 CONTINUE 260 CONTINUE ELSE C C MANUAL ROTATIONS. C SABDIA(1,1) = SABLOC(1,1) IF(NORBSB.GT.1)THEN SABDIA(1,2) = SABLOC(1,4)*C(4,2) SABDIA(1,3) = SABLOC(1,4)*C(4,3) SABDIA(1,4) = SABLOC(1,4)*C(4,4) ENDIF IF(NORBSA.GT.1)THEN SABDIA(2,1) = SABLOC(4,1)*C(4,2) SABDIA(3,1) = SABLOC(4,1)*C(4,3) SABDIA(4,1) = SABLOC(4,1)*C(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) 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) 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) C SABDIA(3,2) = SABDIA(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) 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) C SABDIA(4,2) = SABDIA(2,4) SABDIA(4,3) = SABDIA(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) 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) SABDIA(I,J) = SABDIA(J,I) SABDIA(J,I) = STEMP 480 CONTINUE 500 CONTINUE ENDIF RETURN END C C C DOUBLE PRECISION FUNCTION SAB(NCENTR,NA,LA,XIA,EXHLA, . NB,LB,XIB,EXHLB, . NEWEXP,M,RAB,RNAB) 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 DIMENSION 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) 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) SUM1 = 1.0D0/(3.0D0*RPLUS) SUM2 = (RPLUS + 1.0D0)/RPLUS**3 SAB = EXPLUS*(SUM1 + SUM2)*(XIA**3)*RNAB 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 ABINTG(NAB,RPLUS,RMINUS,ATERM,BTERM) DO 120 I=0,NAB DO 110 J=0,NAB ABTERM(I,J) = ATERM(I)*BTERM(J) 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 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) SUMJA = SUMJA + T345*AB*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 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 DO 150 IQA=0,1 NASUB = NASUB6 + IQA NBSUB = NBSUB6 + IQA IEXP = IEXP6 AB = ABTERM(NASUB,NBSUB) SUMQA = SUMQA + AB*NEG1(IEXP) 150 CONTINUE SUMJA = SUMJA + T345*SUMQA 160 CONTINUE 170 CONTINUE ELSEIF(LA.EQ.0.AND.LB.EQ.1)THEN C C S(A)-P(B) OVERLAP. C SUMJA = 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 DO 180 IQB=0,1 NASUB = NASUB6 + IQB NBSUB = NBSUB6 + IQB IEXP = IEXP6 + IQB AB = ABTERM(NASUB,NBSUB) SUMQB = SUMQB + AB*NEG1(IEXP) 180 CONTINUE SUMJA = SUMJA + SUMQB*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 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 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) SUMQA = SUMQA + AB*NEG1(IEXP) 210 CONTINUE 220 CONTINUE SUMKA = SUMKA + SUMQA*T345 230 CONTINUE 240 CONTINUE 250 CONTINUE 260 CONTINUE SUMJA = SUMKA*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 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 DO 420 JB=0,JBMAX JB2 = 2*JB IPBMAX = IPBM0 + JB2 IQBMAX = IQBM0 - JB2 NBSUB2 = NBSUB1 + JB2 IEXP2 = IEXP0 + JB2 SUMKA = 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) SUMKA = SUMKA + T12345*AB*NEG1(IEXP) 300 CONTINUE 320 CONTINUE 340 CONTINUE 360 CONTINUE 380 CONTINUE 400 CONTINUE SUMJB = SUMJB + SUMKA*CTERM(LB,M,JB) 420 CONTINUE SUMJA = SUMJA + SUMJB*CTERM(LA,M,JA) 440 CONTINUE ENDIF TERMXI = EXHLA*EXHLB TERMNL = SQRTL(LA,LB)/SQRTF(NA,NB) SAB = TERMXI*TERMNL*RNAB*SUMJA 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 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) 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 ABINTG(N,XA,XB,ATERM,BTERM) C C COMPUTES A-TYPE AND B-TYPE INTEGRALS. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ATERM(0:15),BTERM(0:15) C C LOCAL VARIABLES: C DIMENSION 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) 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) DO 130 I=1,N ATERM(I) = ATERM(0) + I*ATERM(I-1)*XAINV 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) 140 CONTINUE ELSEIF(AXB.GT.3.0D0)THEN XBINV = 1.0D0/XB EXB = EXP(XB) EXBINV = 1.0D0/EXB BTERM(0) = XBINV*(EXB - EXBINV) DO 160 I=1,N BTERM(I) = XBINV*(I*BTERM(I-1) + NEG1(I)*EXB - EXBINV) 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 DO 170 J=2,JLAST XBEXP(J) = -XB*XBEXP(J-1) 170 CONTINUE DO 200 I=0,N BI = BFACTR(I,0) DO 180 J=1,JLAST BI = BI + XBEXP(J)*BFACTR(I,J) 180 CONTINUE BTERM(I) = BI 200 CONTINUE ENDIF RETURN END C C C SUBROUTINE SLIMIT(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 DIMENSION 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 DIMENSION 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 RTRANS(XA,YA,ZA,XB,YB,ZB,RAB,NORBS,C) 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) DIMENSION C(9,9) C C LOCAL ARRAYS: C DIMENSION ROT(3,3),D(5,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 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 RMATRX(XA,YA,ZA,XB,YB,ZB,RAB,ROT) 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) 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 IF(I.EQ.J) D(1,I,J) = D(1,I,J) - SQ6 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 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 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 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 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 C(6,K4) = D(K,1,3)*SQ2 + D(K,3,1)*SQ2 C(7,K4) = D(K,2,3)*SQ2 + D(K,3,2)*SQ2 C(8,K4) = D(K,1,1)*SQ2 - D(K,2,2)*SQ2 C(9,K4) = D(K,1,2)*SQ2 + D(K,2,1)*SQ2 120 CONTINUE 1000 RETURN END C C C SUBROUTINE RMATRX(XA,YA,ZA,XB,YB,ZB,RAB,ROT) 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) DIMENSION ROT(3,3),X(3),Y(3),Z(3) Z(1) = (XB - XA) Z(2) = (YB - YA) Z(3) = (ZB - ZA) ZINV = 1.0D0/RAB Z(1) = Z(1)*ZINV Z(2) = Z(2)*ZINV Z(3) = Z(3)*ZINV 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 Z(2) = 0.0D0 Z(3) = SIGN(1.0D0,Z(3)) X(1) = 1.0D0 X(2) = 0.0D0 X(3) = 0.0D0 Y(1) = 0.0D0 Y(2) = Z(3) Y(3) = 0.0D0 ELSE Z12SQR = Z(1)**2 + Z(2)**2 Z12 = DSQRT(Z12SQR) X(1) = -Z(2)/Z12 X(2) = Z(1)/Z12 X(3) = 0.0D0 IF(ABS(Z(3)).LT.1.0D-8)THEN Y(1) = 0.0D0 Y(2) = 0.0D0 Y(3) = 1.0D0 ELSE YINV = 1.0D0/(Z12*DSQRT(1.0D0 + Z12SQR/Z(3)**2)) Y(1) = Z(1)*YINV Y(2) = Z(2)*YINV Y(3) = -(Z12SQR/Z(3))*YINV ENDIF ENDIF DO 10 J=1,3 ROT(1,J) = X(J) ROT(2,J) = Y(J) ROT(3,J) = Z(J) 10 CONTINUE RETURN END