C C $Id: repul.F,v 1.4 1998/07/16 16:40:37 jjv5 Exp arjan $ C C------------------------------------------------------------------------ SUBROUTINE REPUL(IAI,NORBSI,AI0,AI1,AI2,DI1,DI2,XI,YI,ZI, . IAJ,NORBSJ,AJ0,AJ1,AJ2,DJ1,DJ2,XJ,YJ,ZJ, . RIJ,REPIJ) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" #include "constants.h" C C ROUTINE TO CALCULATE REPULSION INTEGRALS FOR A PAIR OF ATOMS I,J. C C INPUT: C C IAI,IAJ = ATOMIC NUMBERS FOR ATOMS I AND J. C C NORBSI, C NORBSJ, = NUMBERS OF ATOMIC ORBITALS CENTERED ON ATOMS I AND J. C C AI0,AI1,AI2, C AJ0,AJ1,AJ2 = MONOPOLE, DIPOLE, AND QUADRUPOLE ADDITIVE TERMS. C C DI1,DI2, C DJ1,DJ2, = DIPOLE AND QUADRUPOLE CHARGE SEPARATIONS. C C XI,YI,ZI, C XJ,YJ,ZJ, = CARTESIAN COORDINATES (ANGSTROMS). C C RIJ = INTERATOMIC DISTANCE (ANGSTROMS). C C C RETURNED: C C REPIJ = 10X10 ARRAY OF REPULSION INTEGRALS ORGANIZED ACCORDING TO: C C S_S S_PX S_PY S_PZ PX_PX PX_PY PX_PZ PY_PY PY_PZ PZ_PZ C S_S ---------------------------------------------------------------- C S_PX | . C S_PY | . C S_PZ | . C PX_PX | . C PX_PY |................(6,4) = (PX(I)*PY(I)|S(J)*PZ(J)) C PX_PZ | C PY_PY | C PY_PZ | C PZ_PZ | C DIMENSION REPIJ(10,10) C C LOCAL VARIABLES: C CHARACTER QCODE*6 DIMENSION REPLOC(22),SQROOT(100) C C EE = CONVERSION FACTOR FROM ATOMIC UNITS TO eV. C DATA EE /27.210000D0/ DATA EEHALF /13.605000D0/ DATA EE4TH / 6.802500D0/ DATA EE8TH / 3.401250D0/ DATA EE16TH / 1.700625D0/ SAVE EE,EEHALF,EE4TH,EE8TH,EE16TH C DO 20 I=1,10 DO 10 J=1,10 REPIJ(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE C C CONVERT RIJ TO ATOMIC UNITS. C A0 = Bohr2Ang R = RIJ/A0 RSQR = R*R C C ASSIGN SEMIEMPIRICAL "ADDITIVE" TERMS THAT DETERMINE ONE-CENTER C LIMITS FOR TWO-CENTER INTEGRALS. C RI0 = 0.5D0/AI0 RJ0 = 0.5D0/AJ0 IF(NORBSI.GT.1)THEN RI1 = 0.5D0/AI1 RI2 = 0.5D0/AI2 ENDIF IF(NORBSJ.GT.1)THEN RJ1 = 0.5D0/AJ1 RJ2 = 0.5D0/AJ2 ENDIF C C COMPUTE THE NONZERO MULTIPOLE-MULTIPOLE INTERACTIONS FOR A LOCAL C FRAME WITH XI=XJ, YI=YJ, AND ZJ-ZI=R, I.E., C C z C | C | C | C ATOM I |_______ y C /. C / . C x . C . C . C z C | C | C | C ATOM J |_______ y C / C / C x C C C [Q(I),Q(I)]: C A00 = (RI0 + RJ0)**2 SQROOT(1) = DSQRT(RSQR + A00) QQ = EE/SQROOT(1) IF(NORBSI.GT.1)THEN C C [DZ(I),Q(J)]: C A10 = (RI1 + RJ0)**2 SQROOT(2) = DSQRT((R+DI1)**2 + A10) SQROOT(3) = DSQRT((R-DI1)**2 + A10) DZQ = EEHALF/SQROOT(2) - EEHALF/SQROOT(3) C C [QZZ(I),Q(J)]: C A20 = (RI2 + RJ0)**2 DI22 = 2.0D0*DI2 SQROOT(4) = DSQRT((R+DI22)**2 + A20) SQROOT(5) = DSQRT((R-DI22)**2 + A20) SQROOT(6) = DSQRT(RSQR + A20) QZZQ = EE4TH/SQROOT(4) + EE4TH/SQROOT(5) - EEHALF/SQROOT(6) C C [QXX(I),Q(J)] OR (QYY(I),Q(J)]: C DI22SQ = DI22*DI22 SQROOT(7) = DSQRT(RSQR + A20 + DI22SQ) SQROOT(8) = SQROOT(6) QXXQ = EEHALF/SQROOT(7) - EEHALF/SQROOT(8) C IF(NORBSJ.GT.1)THEN C C [QXZ(I),DX(J)] OR [QYZ(I),DY(J)]: C A21 = (RI2 + RJ1)**2 RP = (R + DI2)**2 RM = (R - DI2)**2 RPA = RP + A21 RMA = RM + A21 DP = (DI2 + DJ1)**2 DM = (DI2 - DJ1)**2 SQROOT(9) = DSQRT(RPA + DM) SQROOT(10) = DSQRT(RMA + DM) SQROOT(11) = DSQRT(RMA + DP) SQROOT(12) = DSQRT(RPA + DP) QXZDX = EE4TH/SQROOT(9) - EE4TH/SQROOT(10) . + EE4TH/SQROOT(11) - EE4TH/SQROOT(12) C C [QZZ(I),DZ(J)]: C RP = R + DJ1 RPSQR = RP*RP RM = R - DJ1 RMSQR = RM*RM RPP = (RP + DI22)**2 RMM = (RM - DI22)**2 RPM = (RP - DI22)**2 RMP = (RM + DI22)**2 SQROOT(13) = DSQRT(RMM + A21) SQROOT(14) = DSQRT(RPP + A21) SQROOT(15) = DSQRT(RMP + A21) SQROOT(16) = DSQRT(RPM + A21) SQROOT(17) = DSQRT(RMSQR + A21) SQROOT(18) = DSQRT(RPSQR + A21) QZZDZ = EE8TH/SQROOT(13) - EE8TH/SQROOT(14) + EE8TH/SQROOT(15) . - EE8TH/SQROOT(16) - EE4TH/SQROOT(17) + EE4TH/SQROOT(18) C C [QXX(I),DZ(J)] OR [QYY(I),DZ(J)]: C DI2A = DI22SQ + A21 SQROOT(19) = DSQRT(RMSQR + DI2A) SQROOT(20) = DSQRT(RPSQR + DI2A) SQROOT(21) = SQROOT(17) SQROOT(22) = SQROOT(18) QXXDZ = EE4TH/SQROOT(19) - EE4TH/SQROOT(20) . - EE4TH/SQROOT(21) + EE4TH/SQROOT(22) ENDIF ENDIF IF(NORBSJ.GT.1)THEN C C [Q(I),DZ(J)]: C IF(IAI.EQ.IAJ)THEN QDZ = -DZQ ELSE A01 = (RI0 + RJ1)**2 SQROOT(23) = DSQRT((R+DJ1)**2 + A01) SQROOT(24) = DSQRT((R-DJ1)**2 + A01) QDZ = -EEHALF/SQROOT(23) + EEHALF/SQROOT(24) ENDIF C C [Q(I),QZZ(J)]: C DJ22 = 2.0D0*DJ2 IF(IAI.EQ.IAJ)THEN QQZZ = QZZQ ELSE A02 = (RI0 + RJ2)**2 SQROOT(25) = DSQRT((R+DJ22)**2 + A02) SQROOT(26) = DSQRT((R-DJ22)**2 + A02) SQROOT(27) = DSQRT(RSQR + A02) QQZZ = EE4TH/SQROOT(25) + EE4TH/SQROOT(26) - EEHALF/SQROOT(27) ENDIF C C [Q(I),QXX(J)]: C DJ22SQ = DJ22*DJ22 IF(IAI.EQ.IAJ)THEN QQXX = QXXQ ELSE SQROOT(28) = DSQRT(RSQR + A02 + DJ22SQ) SQROOT(29) = SQROOT(27) QQXX = EEHALF/SQROOT(28) - EEHALF/SQROOT(29) ENDIF IF(NORBSI.GT.1)THEN C C [DX(I),QXZ(J)] OR [DY(I),QYZ(J)]: C IF(IAI.EQ.IAJ)THEN DXQXZ = -QXZDX ELSE A12 = (RI1 + RJ2)**2 RP = (R + DJ2)**2 RM = (R - DJ2)**2 RPA = RP + A12 RMA = RM + A12 DP = (DJ2+DI1)**2 DM = (DJ2-DI1)**2 SQROOT(30) = DSQRT(RPA + DM) SQROOT(31) = DSQRT(RMA + DM) SQROOT(32) = DSQRT(RMA + DP) SQROOT(33) = DSQRT(RPA + DP) DXQXZ = -EE4TH/SQROOT(30) + EE4TH/SQROOT(31) . -EE4TH/SQROOT(32) + EE4TH/SQROOT(33) ENDIF C C [DZ(I),QZZ(J)]: C IF(IAI.EQ.IAJ)THEN DZQZZ = -QZZDZ ELSE RP = R + DI1 RPSQR = RP*RP RM = R - DI1 RMSQR = RM*RM RPP = (RP + DJ22)**2 RMM = (RM - DJ22)**2 RPM = (RP - DJ22)**2 RMP = (RM + DJ22)**2 SQROOT(34) = DSQRT(RMM + A12) SQROOT(35) = DSQRT(RPP + A12) SQROOT(36) = DSQRT(RMP + A12) SQROOT(37) = DSQRT(RPM + A12) SQROOT(38) = DSQRT(RMSQR + A12) SQROOT(39) = DSQRT(RPSQR + A12) DZQZZ = -EE8TH/SQROOT(34) + EE8TH/SQROOT(35) . -EE8TH/SQROOT(36) + EE8TH/SQROOT(37) . +EE4TH/SQROOT(38) - EE4TH/SQROOT(39) ENDIF C C [DZ(I),QXX(J)] OR [DZ(I),QYY(J)]: C IF(IAI.EQ.IAJ)THEN DZQXX = -QXXDZ ELSE DJ2A = DJ22SQ + A12 SQROOT(40) = DSQRT(RMSQR + DJ2A) SQROOT(41) = DSQRT(RPSQR + DJ2A) SQROOT(42) = SQROOT(38) SQROOT(43) = SQROOT(39) DZQXX = -EE4TH/SQROOT(40) + EE4TH/SQROOT(41) . +EE4TH/SQROOT(42) - EE4TH/SQROOT(43) ENDIF C C [DX(I),DX(J)] OR [DY(I),DY(J)]: C A11 = (RI1 + RJ1)**2 RA = RSQR + A11 DIPJ = DI1 + DJ1 DIMJ = DI1 - DJ1 SQROOT(44) = DSQRT(RA + DIMJ*DIMJ) SQROOT(45) = DSQRT(RA + DIPJ*DIPJ) DXDX = EEHALF/SQROOT(44) - EEHALF/SQROOT(45) C C [DZ(I),DZ(J)]: C SQROOT(46) = DSQRT((R+DIMJ)**2 + A11) SQROOT(47) = DSQRT((R+DIPJ)**2 + A11) SQROOT(48) = DSQRT((R-DIPJ)**2 + A11) SQROOT(49) = DSQRT((R-DIMJ)**2 + A11) DZDZ = EE4TH/SQROOT(46) - EE4TH/SQROOT(47) . - EE4TH/SQROOT(48) + EE4TH/SQROOT(49) #ifdef CUTREPUL_IS_ON if (rij.le.cutrepul) then #endif C C [QZZ(I),QZZ(J)]: C A22 = (RI2 + RJ2)**2 RA = RSQR + A22 RPI = R + DI22 RMI = R - DI22 RPJ = R + DJ22 RMJ = R - DJ22 SQROOT(50) = DSQRT((RPI - DJ22)**2 + A22) SQROOT(51) = DSQRT((RPI + DJ22)**2 + A22) SQROOT(52) = DSQRT((RMI - DJ22)**2 + A22) SQROOT(53) = DSQRT((RMI + DJ22)**2 + A22) SQROOT(54) = DSQRT(RPI*RPI + A22) SQROOT(55) = DSQRT(RMI*RMI + A22) SQROOT(56) = DSQRT(RPJ*RPJ + A22) SQROOT(57) = DSQRT(RMJ*RMJ + A22) SQROOT(58) = DSQRT(RA) QZZQZZ = EE16TH/SQROOT(50) + EE16TH/SQROOT(51) . + EE16TH/SQROOT(52) + EE16TH/SQROOT(53) . - EE8TH/SQROOT(54) - EE8TH/SQROOT(55) . - EE8TH/SQROOT(56) - EE8TH/SQROOT(57) . + EE4TH/SQROOT(58) C C [QXX(I),QXX(J)] OR [QYY(I),QYY(I)]: C SQROOT(59) = DSQRT(RA + (DI22-DJ22)**2) SQROOT(60) = DSQRT(RA + (DI22+DJ22)**2) SQROOT(61) = DSQRT(RA + DI22SQ) SQROOT(62) = DSQRT(RA + DJ22SQ) SQROOT(63) = SQROOT(58) QXXQXX = EE8TH/SQROOT(59) + EE8TH/SQROOT(60) . - EE4TH/SQROOT(61) - EE4TH/SQROOT(62) . + EE4TH/SQROOT(63) C C [QZZ(I),QXX(J)] OR [QZZ(I),QYY(I)]: C AJJ2SQ = A22 + DJ22SQ RPISQ = RPI*RPI RMISQ = RMI*RMI SQROOT(64) = DSQRT(RMISQ + AJJ2SQ) SQROOT(65) = DSQRT(RPISQ + AJJ2SQ) SQROOT(66) = SQROOT(55) SQROOT(67) = SQROOT(54) SQROOT(68) = SQROOT(62) SQROOT(69) = SQROOT(58) QZZQXX = EE8TH/SQROOT(64) + EE8TH/SQROOT(65) . - EE8TH/SQROOT(66) - EE8TH/SQROOT(67) . - EE4TH/SQROOT(68) + EE4TH/SQROOT(69) C C [QXX(I),QZZ(J)] OR [QYY(I),QZZ(J)] C IF(IAI.EQ.IAJ)THEN QXXQZZ = QZZQXX ELSE AII2SQ = A22 + DI22SQ RPJSQ = RPJ*RPJ RMJSQ = RMJ*RMJ SQROOT(70) = DSQRT(RMJSQ + AII2SQ) SQROOT(71) = DSQRT(RPJSQ + AII2SQ) SQROOT(72) = SQROOT(57) SQROOT(73) = SQROOT(56) SQROOT(74) = SQROOT(61) SQROOT(75) = SQROOT(58) QXXQZZ = EE8TH/SQROOT(70) + EE8TH/SQROOT(71) . - EE8TH/SQROOT(72) - EE8TH/SQROOT(73) . - EE4TH/SQROOT(74) + EE4TH/SQROOT(75) ENDIF C C [QXX(I),QYY(J)] OR [QYY(I),QXX(J)]: C SQROOT(76) = DSQRT(RA + DI22SQ + DJ22SQ) SQROOT(77) = SQROOT(61) SQROOT(78) = SQROOT(62) SQROOT(79) = SQROOT(58) QXXQYY = EE4TH/SQROOT(76) - EE4TH/SQROOT(77) . - EE4TH/SQROOT(78) + EE4TH/SQROOT(79) C C [QXZ(I),QXZ(J)] OR [QYZ(I),QYZ(J)]: C DIPJ = DI2 + DJ2 DIMJ = DI2 - DJ2 DIPJSQ = DIPJ*DIPJ DIMJSQ = DIMJ*DIMJ RPP = (R+DIPJ)**2 RMP = (R-DIPJ)**2 RPM = (R+DIMJ)**2 RMM = (R-DIMJ)**2 DPSQA = DIPJSQ + A22 DMSQA = DIMJSQ + A22 SQROOT(80) = DSQRT(RPM + DMSQA) SQROOT(81) = DSQRT(RPM + DPSQA) SQROOT(82) = DSQRT(RPP + DMSQA) SQROOT(83) = DSQRT(RPP + DPSQA) SQROOT(84) = DSQRT(RMP + DMSQA) SQROOT(85) = DSQRT(RMP + DPSQA) SQROOT(86) = DSQRT(RMM + DMSQA) SQROOT(87) = DSQRT(RMM + DPSQA) QXZQXZ = EE8TH/SQROOT(80) - EE8TH/SQROOT(81) . - EE8TH/SQROOT(82) + EE8TH/SQROOT(83) . - EE8TH/SQROOT(84) + EE8TH/SQROOT(85) . + EE8TH/SQROOT(86) - EE8TH/SQROOT(87) C C [QXY(I),QXY(J)] OR [QXX-YY(I),QXX-YY(J)]: C QXYQXY = 0.5D0*(QXXQXX - QXXQYY) #ifdef CUTREPUL_IS_ON else qzzqzz = 0.0 qxxqxx = 0.0 qzzqxx = 0.0 qxxqzz = 0.0 qxxqyy = 0.0 qxzqxz = 0.0 qxyqxy = 0.0 endif #endif ENDIF ENDIF C C NOW USE THE LOCAL MULTIPOLE-MULTIPOLE INTERACTIONS TO CONSTRUCT C THE 22 NONZERO LOCAL REPULSION INTEGRALS. C C (S S | S S): C REPLOC(1) = QQ IF(NORBSI.GT.1)THEN C C (S PZ | S S): C REPLOC(2) = DZQ C C (PZ PZ | S S): C REPLOC(3) = QZZQ + QQ C C (PX PX | S S) OR (PY PY | S S): C REPLOC(4) = QXXQ + QQ ENDIF IF(NORBSJ.GT.1)THEN C C (S S | S PZ): C REPLOC(5) = QDZ IF(NORBSI.GT.1)THEN C C (S PZ | S PZ): C REPLOC(6) = DZDZ C C (S PX | S PX) OR (S PY | S PY): C REPLOC(7) = DXDX C C (PZ PZ | S PZ): C REPLOC(8) = QZZDZ + QDZ C C (PX PX | S PZ) OR (PY PY | S PZ): C REPLOC(9) = QXXDZ + QDZ C C (PX PZ | S PX): C REPLOC(10) = QXZDX ENDIF C C (S S | PZ PZ): C REPLOC(11) = QQZZ + QQ C C (S S | PX PX) OR (S S | PY PY): C REPLOC(12) = QQXX + QQ IF(NORBSI.GT.1)THEN C C (S PZ | PZ PZ): C REPLOC(13) = DZQZZ + DZQ C C (S PZ | PX PX) OR (S PZ | PY PY): C REPLOC(14) = DZQXX + DZQ C C (S PX | PX PZ) OR (S PY | PY PZ): C REPLOC(15) = DXQXZ QQQQZZ = QQ + QQZZ C C (PZ PZ | PZ PZ): C REPLOC(16) = QQQQZZ + QZZQ + QZZQZZ C C (PX PX | PZ PZ) OR (PY PY | PZ PZ): C REPLOC(17) = QQQQZZ + QXXQ + QXXQZZ QQQQXX = QQ + QQXX C C (PZ PZ | PX PX) OR (PZ PZ | PY PY): C REPLOC(18) = QQQQXX + QZZQ + QZZQXX C C (PX PX | PX PX) OR (PY PY | PY PY): C REPLOC(19) = QQQQXX + QXXQ + QXXQXX C C (PX PZ | PX PZ) OR (PY PZ | PY PZ): C REPLOC(20) = QXZQXZ C C (PX PX | PY PY) OR (PY PY | PX PX): C REPLOC(21) = QQQQXX + QXXQ + QXXQYY C C (PX PY | PX PY): C REPLOC(22) = QXYQXY ENDIF ENDIF C C ROTATE LOCAL REPULSION INTEGRALS TO MOLECULAR FRAME. C IF(NORBSI.EQ.1.AND.NORBSJ.EQ.1)THEN REPIJ(1,1) = REPLOC(1) ELSE CALL RREP(NORBSI,NORBSJ,XI,YI,ZI,XJ,YJ,ZJ,RIJ,REPLOC,REPIJ) ENDIF RETURN END C C C SUBROUTINE RREP(NORBSI,NORBSJ,XI,YI,ZI,XJ,YJ,ZJ,RIJ,REPLOC,REPIJ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION REPLOC(22),REPIJ(10,10) C C ROTATES THE LOCAL FRAME REPULSION INTEGRALS TO THE MOLECULAR FRAME. C DIMENSION R(3,3) C C GENERATE ROTATION MATRIX. C R(3,1) = (XI-XJ)/RIJ R(3,2) = (YI-YJ)/RIJ R(3,3) = (ZI-ZJ)/RIJ IF(ABS(R(3,3)).GT.0.99999999D0)THEN R(3,1) = 0.0D0 R(3,2) = 0.0D0 R(3,3) = SIGN(1.0D0,R(3,3)) R(1,1) = 1.0D0 R(1,2) = 0.0D0 R(1,3) = 0.0D0 R(2,1) = 0.0D0 R(2,2) = R(3,3) R(2,3) = 0.0D0 ELSE Z12SQR = R(3,1)**2 + R(3,2)**2 Z12 = DSQRT(Z12SQR) R(1,1) = -R(3,2)/Z12 R(1,2) = R(3,1)/Z12 R(1,3) = 0.0D0 IF(ABS(R(3,3)).LT.1.0D-8)THEN R(2,1) = 0.0D0 R(2,2) = 0.0D0 R(2,3) = 1.0D0 ELSE YINV = 1.0D0/(Z12*DSQRT(1.0D0 + Z12SQR/R(3,3)**2)) R(2,1) = R(3,1)*YINV R(2,2) = R(3,2)*YINV R(2,3) = -(Z12SQR/R(3,3))*YINV ENDIF ENDIF R11R11 = R(1,1)*R(1,1) R11R12 = R(1,1)*R(1,2) R11R13 = R(1,1)*R(1,3) R12R12 = R(1,2)*R(1,2) R12R13 = R(1,2)*R(1,3) R13R13 = R(1,3)*R(1,3) R21R21 = R(2,1)*R(2,1) R21R22 = R(2,1)*R(2,2) R21R23 = R(2,1)*R(2,3) R22R22 = R(2,2)*R(2,2) R22R23 = R(2,2)*R(2,3) R23R23 = R(2,3)*R(2,3) R31R31 = R(3,1)*R(3,1) R31R32 = R(3,1)*R(3,2) R31R33 = R(3,1)*R(3,3) R32R32 = R(3,2)*R(3,2) R32R33 = R(3,2)*R(3,3) R33R33 = R(3,3)*R(3,3) IF(NORBSI.GT.1.OR.NORBSJ.GT.1)THEN RADD01 = R11R11+R21R21 RADD02 = R11R12+R21R22 RADD03 = R11R13+R21R23 RADD04 = R12R12+R22R22 RADD05 = R12R13+R22R23 RADD06 = R13R13+R23R23 ENDIF IF(NORBSI.GT.1.AND.NORBSJ.GT.1)THEN R11R21 = R(1,1)*R(2,1)*2.0D0 R11R22 = R(1,1)*R(2,2) R11R23 = R(1,1)*R(2,3) R11R31 = R(1,1)*R(3,1)*2.0D0 R11R32 = R(1,1)*R(3,2) R11R33 = R(1,1)*R(3,3) R12R22 = R(1,2)*R(2,2)*2.0D0 R12R23 = R(1,2)*R(2,3) R12R31 = R(1,2)*R(3,1) R12R32 = R(1,2)*R(3,2)*2.0D0 R12R33 = R(1,2)*R(3,3) R13R23 = R(1,3)*R(2,3)*2.0D0 R13R33 = R(1,3)*R(3,3)*2.0D0 R21R11 = R11R21 R21R12 = R(2,1)*R(1,2) R21R13 = R(2,1)*R(1,3) R21R21 = R(2,1)*R(2,1) R21R31 = R(2,1)*R(3,1)*2.0D0 R21R32 = R(2,1)*R(3,2) R21R33 = R(2,1)*R(3,3) R22R12 = R12R22 R22R13 = R(2,2)*R(1,3) R22R31 = R(2,2)*R(3,1) R22R32 = R(2,2)*R(3,2)*2.0D0 R22R33 = R(2,2)*R(3,3) R23R13 = R13R23 R23R33 = R(2,3)*R(3,3)*2.0D0 R31R11 = R11R31 R31R12 = R(3,1)*R(1,2) R31R13 = R(3,1)*R(1,3) R31R21 = R21R31 R31R22 = R22R31 R31R23 = R(3,1)*R(2,3) R32R12 = R12R32 R32R13 = R(3,2)*R(1,3) R32R22 = R22R32 R32R23 = R(3,2)*R(2,3) R33R13 = R13R33 R33R23 = R23R33 RADD07 = R11R32+R12R31 RADD08 = R21R32+R22R31 RADD09 = R11R33+R31R13 RADD10 = R21R33+R31R23 RADD11 = R11R32+R31R12 RADD12 = R21R32+R31R22 RADD13 = R12R33+R32R13 RADD14 = R22R33+R32R23 RADD15 = R11R23+R21R13 RADD16 = R12R23+R22R13 RADD17 = R11R22+R21R12 ENDIF C C (S S | S S): C REPIJ(1,1) = REPLOC(1) IF(NORBSJ.GT.1)THEN C C (S S | S PX): C REPIJ(1,2) = REPLOC(5)*R(3,1) C C (S S | S PY): C REPIJ(1,3) = REPLOC(5)*R(3,2) C C (S S | S PZ): C REPIJ(1,4) = REPLOC(5)*R(3,3) C C (S S | PX PX): C REPIJ(1,5) = REPLOC(11)*R31R31 . + REPLOC(12)*RADD01 C C (S S | PX PY): C REPIJ(1,6) = REPLOC(11)*R31R32 . + REPLOC(12)*RADD02 C C (S S | PX PZ): C REPIJ(1,7) = REPLOC(11)*R31R33 . + REPLOC(12)*RADD03 C C (S S | PY PY): C REPIJ(1,8) = REPLOC(11)*R32R32 . + REPLOC(12)*RADD04 C C (S S | PY PZ): C REPIJ(1,9) = REPLOC(11)*R32R33 . + REPLOC(12)*RADD05 C C (S S | PZ PZ): C REPIJ(1,10) = REPLOC(11)*R33R33 . + REPLOC(12)*RADD06 ENDIF IF(NORBSI.GT.1)THEN C C (S PX | S S): C REPIJ(2,1) = REPLOC(2)*R(3,1) C C (S PY | S S): C REPIJ(3,1) = REPLOC(2)*R(3,2) C C (S PZ | S S): C REPIJ(4,1) = REPLOC(2)*R(3,3) C C (PX PX | S S): C REPIJ(5,1) = REPLOC(3)*R31R31 . + REPLOC(4)*RADD01 C C (PX PY | S S): C REPIJ(6,1) = REPLOC(3)*R31R32 . + REPLOC(4)*RADD02 C C (PX PZ | S S): C REPIJ(7,1) = REPLOC(3)*R31R33 . + REPLOC(4)*RADD03 C C (PY PY | S S): C REPIJ(8,1) = REPLOC(3)*R32R32 . + REPLOC(4)*RADD04 C C (PY PZ | S S); C REPIJ(9,1) = REPLOC(3)*R32R33 . + REPLOC(4)*RADD05 C C (PZ PZ | S S): C REPIJ(10,1) = REPLOC(3)*R33R33 . + REPLOC(4)*RADD06 ENDIF IF(NORBSI.GT.1.AND.NORBSJ.GT.1)THEN C C (S PX | S PX): C REPIJ(2,2) = REPLOC(6)*R31R31 . + REPLOC(7)*RADD01 C C (S PX | S PY): C REPIJ(2,3) = REPLOC(6)*R31R32 . + REPLOC(7)*RADD02 C C (S PY | S PX): C REPIJ(3,2) = REPIJ(2,3) C C (S PX | S PZ): C REPIJ(2,4) = REPLOC(6)*R31R33 . + REPLOC(7)*RADD03 C C (S PZ | S PX): C REPIJ(4,2) = REPIJ(2,4) C C (S PX | PX PX): C REPIJ(2,5) = REPLOC(13)*R(3,1)*R31R31 . + REPLOC(14)*R(3,1)*RADD01 . + REPLOC(15)*(R(1,1)*R11R31 . + R(2,1)*R21R31) C C (PX PX | S PX): C REPIJ(5,2) = REPLOC(8)*R(3,1)*R31R31 . + REPLOC(9)*R(3,1)*RADD01 . + REPLOC(10)*(R(1,1)*R11R31 . + R(2,1)*R21R31) C C (S PX | PX PY): C REPIJ(2,6) = REPLOC(13)*R(3,1)*R31R32 . + REPLOC(14)*R(3,1)*RADD02 . + REPLOC(15)*(R(1,1)*RADD07 . + R(2,1)*RADD08) C C (PX PY | S PX): C REPIJ(6,2) = REPLOC(8)*R(3,1)*R31R32 . + REPLOC(9)*R(3,1)*RADD02 . + REPLOC(10)*(R(1,1)*RADD07 . + R(2,1)*RADD08) C C (S PX | PX PZ): C REPIJ(2,7) = REPLOC(13)*R(3,1)*R31R33 . + REPLOC(14)*R(3,1)*RADD03 . + REPLOC(15)*(R(1,1)*RADD09 . + R(2,1)*RADD10) C C (PX PZ | S PX): C REPIJ(7,2) = REPLOC(8)*R(3,1)*R31R33 . + REPLOC(9)*(R(3,1)*R11R13 . + R(3,1)*R21R23) . + REPLOC(10)*(R(1,1)*RADD09 . + R(2,1)*RADD10) C C (S PX | PY PY): C REPIJ(2,8) = REPLOC(13)*R(3,1)*R32R32 . + REPLOC(14)*R(3,1)*RADD04 . + REPLOC(15)*(R(1,1)*R12R32 . + R(2,1)*R22R32) C C (PY PY | S PX): C REPIJ(8,2) = REPLOC(8)*R(3,1)*R32R32 . + REPLOC(9)*R(3,1)*RADD04 . + REPLOC(10)*(R(1,1)*R12R32 . + R(2,1)*R22R32) C C (S PX | PY PZ): C C REPIJ(2,9) = REPLOC(13)*R(3,1)*R32R33 . + REPLOC(14)*R(3,1)*RADD05 . + REPLOC(15)*(R(1,1)*RADD13 . + R(2,1)*RADD14) C C (PY PZ | S PX): C REPIJ(9,2) = REPLOC(8)*R(3,1)*R32R33 . + REPLOC(9)*R(3,1)*RADD05 . + REPLOC(10)*(R(1,1)*RADD13 . + R(2,1)*RADD14) C C (S PX | PZ PZ): C REPIJ(2,10) = REPLOC(13)*R(3,1)*R33R33 . + REPLOC(14)*R(3,1)*RADD06 . + REPLOC(15)*(R(1,1)*R13R33 . + R(2,1)*R23R33) C C (PZ PZ | S PX): C REPIJ(10,2) = REPLOC(8)*R(3,1)*R33R33 . + REPLOC(9)*R(3,1)*RADD06 . + REPLOC(10)*(R(1,1)*R13R33 . + R(2,1)*R23R33) C C (S PY | S PY): C REPIJ(3,3) = REPLOC(6)*R32R32 . + REPLOC(7)*RADD04 C C (S PY | S PZ): C REPIJ(3,4) = REPLOC(6)*R32R33 . + REPLOC(7)*RADD05 C C (S PZ | S PY): C REPIJ(4,3) = REPIJ(3,4) C C (S PY | PX PX): C REPIJ(3,5) = REPLOC(13)*R(3,2)*R31R31 . + REPLOC(14)*R(3,2)*RADD01 . + REPLOC(15)*(R(1,2)*R11R31 . + R(2,2)*R21R31) C C (PX PX | S PY): C REPIJ(5,3) = REPLOC(8)*R(3,2)*R31R31 . + REPLOC(9)*R(3,2)*RADD01 . + REPLOC(10)*(R(1,2)*R11R31 . + R(2,2)*R21R31) C C (S PY | PX PY): C REPIJ(3,6) = REPLOC(13)*R(3,2)*R31R32 . + REPLOC(14)*R(3,2)*RADD02 . + REPLOC(15)*(R(1,2)*RADD11 . + R(2,2)*RADD12) C C (PX PY | S PY): C REPIJ(6,3) = REPLOC(8)*R(3,2)*R31R32 . + REPLOC(9)*R(3,2)*RADD02 . + REPLOC(10)*(R(1,2)*RADD11 . + R(2,2)*RADD12) C C (S PY | PX PZ): C REPIJ(3,7) = REPLOC(13)*R(3,2)*R31R33 . + REPLOC(14)*R(3,2)*RADD03 . + REPLOC(15)*(R(1,2)*RADD09 . + R(2,2)*RADD10) C C (PX PZ | S PY): C REPIJ(7,3) = REPLOC(8)*R(3,2)*R31R33 . + REPLOC(9)*R(3,2)*RADD03 . + REPLOC(10)*(R(1,2)*RADD09 . + R(2,2)*RADD10) C C (S PY | PY PY): C REPIJ(3,8) = REPLOC(13)*R(3,2)*R32R32 . + REPLOC(14)*R(3,2)*RADD04 . + REPLOC(15)*(R(1,2)*R12R32 . + R(2,2)*R22R32) C C (PY PY | S PY): C REPIJ(8,3) = REPLOC(8)*R(3,2)*R32R32 . + REPLOC(9)*R(3,2)*RADD04 . + REPLOC(10)*(R(1,2)*R12R32 . + R(2,2)*R22R32) C C (S PY | PY PZ): C REPIJ(3,9) = REPLOC(13)*R(3,2)*R32R33 . + REPLOC(14)*R(3,2)*RADD05 . + REPLOC(15)*(R(1,2)*RADD13 . + R(2,2)*RADD14) C C (PY PZ | S PY): C REPIJ(9,3) = REPLOC(8)*R(3,2)*R32R33 . + REPLOC(9)*R(3,2)*RADD05 . + REPLOC(10)*(R(1,2)*RADD13 . + R(2,2)*RADD14) C C (S PY | PZ PZ): C REPIJ(3,10) = REPLOC(13)*R(3,2)*R33R33 . + REPLOC(14)*R(3,2)*RADD06 . + REPLOC(15)*(R(1,2)*R13R33 . + R(2,2)*R23R33) C C (PZ PZ | S PY): C REPIJ(10,3) = REPLOC(8)*R(3,2)*R33R33 . + REPLOC(9)*R(3,2)*RADD06 . + REPLOC(10)*(R(1,2)*R13R33 . + R(2,2)*R23R33) C C (S PZ | S PZ): C REPIJ(4,4) = REPLOC(6)*R33R33 . + REPLOC(7)*RADD06 C C (S PZ | PX PX): C REPIJ(4,5) = REPLOC(13)*R(3,3)*R31R31 . + REPLOC(14)*R(3,3)*RADD01 . + REPLOC(15)*(R(1,3)*R11R31 . + R(2,3)*R21R31) C C (PX PX | S PZ): C REPIJ(5,4) = REPLOC(8)*R(3,3)*R31R31 . + REPLOC(9)*R(3,3)*RADD01 . + REPLOC(10)*(R(1,3)*R11R31 . + R(2,3)*R21R31) C C (S PZ | PX PY): C REPIJ(4,6) = REPLOC(13)*R(3,3)*R31R32 . + REPLOC(14)*R(3,3)*RADD02 . + REPLOC(15)*(R(1,3)*RADD11 . + R(2,3)*RADD12) C C (PX PY | S PZ): C REPIJ(6,4) = REPLOC(8)*R(3,3)*R31R32 . + REPLOC(9)*R(3,3)*RADD02 . + REPLOC(10)*(R(1,3)*RADD11 . + R(2,3)*RADD12) C C (S PZ | PX PZ): C REPIJ(4,7) = REPLOC(13)*R(3,3)*R31R33 . + REPLOC(14)*R(3,3)*RADD03 . + REPLOC(15)*(R(1,3)*RADD09 . + R(2,3)*RADD10) C C (PX PZ | S PZ): C REPIJ(7,4) = REPLOC(8)*R(3,3)*R31R33 . + REPLOC(9)*R(3,3)*RADD03 . + REPLOC(10)*(R(1,3)*RADD09 . + R(2,3)*RADD10) C C (S PZ | PY PY): C REPIJ(4,8) = REPLOC(13)*R(3,3)*R32R32 . + REPLOC(14)*R(3,3)*RADD04 . + REPLOC(15)*(R(1,3)*R12R32 . + R(2,3)*R22R32) C C (PY PY | S PZ): C REPIJ(8,4) = REPLOC(8)*R(3,3)*R32R32 . + REPLOC(9)*R(3,3)*RADD04 . + REPLOC(10)*(R(1,3)*R12R32 . + R(2,3)*R22R32) C C (S PZ | PY PZ): C REPIJ(4,9) = REPLOC(13)*R(3,3)*R32R33 . + REPLOC(14)*R(3,3)*RADD05 . + REPLOC(15)*(R(1,3)*RADD13 . + R(2,3)*RADD14) C C (PY PZ | S PZ): C REPIJ(9,4) = REPLOC(8)*R(3,3)*R32R33 . + REPLOC(9)*R(3,3)*RADD05 . + REPLOC(10)*(R(1,3)*RADD13 . + R(2,3)*RADD14) C C (S PZ | PZ PZ): C REPIJ(4,10) = REPLOC(13)*R(3,3)*R33R33 . + REPLOC(14)*R(3,3)*RADD06 . + REPLOC(15)*(R(1,3)*R13R33 . + R(2,3)*R23R33) C C (PZ PZ | S PZ): C REPIJ(10,4) = REPLOC(8)*R(3,3)*R33R33 . + REPLOC(9)*R(3,3)*RADD06 . + REPLOC(10)*(R(1,3)*R13R33 . + R(2,3)*R23R33) C C (PX PX | PX PX): C REPIJ(5,5) = REPLOC(16)*R31R31*R31R31 . + REPLOC(17)*RADD01*R31R31 . + REPLOC(18)*RADD01*R31R31 . + REPLOC(19)*(R11R11*R11R11 . + R21R21*R21R21) . + REPLOC(20)*(R11R31*R11R31 . + R21R31*R21R31) . + REPLOC(21)*R11R11*R21R21*2.0D0 . + REPLOC(22)*R11R21*R11R21 C C (PX PX | PX PY): C REPIJ(5,6) = REPLOC(16)*R31R31*R31R32 . + REPLOC(17)*RADD01*R31R32 . + REPLOC(18)*RADD02*R31R31 . + REPLOC(19)*(R11R11*R11R12 . + R21R21*R21R22) . + REPLOC(20)*(R11R31*R11R32 . + R31R11*R31R12 . + R21R31*R21R32 . + R31R21*R31R22) . + REPLOC(21)*(R11R11*R21R22 . + R21R21*R11R12) . + REPLOC(22)*(R11R21*R11R22 . + R21R11*R21R12) C C (PX PY | PX PX): C REPIJ(6,5) = REPLOC(16)*R31R32*R31R31 . + REPLOC(17)*RADD02*R31R31 . + REPLOC(18)*RADD01*R31R32 . + REPLOC(19)*(R11R12*R11R11 . + R21R22*R21R21) . + REPLOC(20)*(RADD11*R11R31 . + RADD12*R21R31) . + REPLOC(21)*(R11R12*R21R21 . + R21R22*R11R11) . + REPLOC(22)*RADD17*R11R21 C C (PX PX | PX PZ): C REPIJ(5,7) = REPLOC(16)*R31R31*R31R33 . + REPLOC(17)*RADD01*R31R33 . + REPLOC(18)*RADD03*R31R31 . + REPLOC(19)*(R11R11*R11R13 . + R21R21*R21R23) . + REPLOC(20)*(R11R31*R11R33 . + R31R11*R31R13 . + R21R31*R21R33 . + R31R21*R31R23) . + REPLOC(21)*(R11R11*R21R23 . + R21R21*R11R13) . + REPLOC(22)*(R11R21*R11R23 . + R21R11*R21R13) C C (PX PZ | PX PX): C REPIJ(7,5) = REPLOC(16)*R31R33*R31R31 . + REPLOC(17)*RADD03*R31R31 . + REPLOC(18)*RADD01*R31R33 . + REPLOC(19)*(R11R13*R11R11 . + R21R23*R21R21) . + REPLOC(20)*(RADD09*R11R31 . +RADD10*R21R31) . + REPLOC(21)*(R11R13*R21R21 . + R21R23*R11R11) . + REPLOC(22)*RADD15*R11R21 C C (PX PX | PY PY): C REPIJ(5,8) = REPLOC(16)*R31R31*R32R32 . + REPLOC(17)*RADD01*R32R32 . + REPLOC(18)*RADD04*R31R31 . + REPLOC(19)*(R11R11*R12R12 . + R21R21*R22R22) . + REPLOC(20)*(R11R31*R12R32 . + R21R31*R22R32) . + REPLOC(21)*(R11R11*R22R22 . + R21R21*R12R12) . + REPLOC(22)*R11R21*R12R22 C C (PY PY | PX PX): C REPIJ(8,5) = REPLOC(16)*R32R32*R31R31 . + REPLOC(17)*RADD04*R31R31 . + REPLOC(18)*RADD01*R32R32 . + REPLOC(19)*(R12R12*R11R11 . + R22R22*R21R21) . + REPLOC(20)*(R12R32*R11R31 . + R22R32*R21R31) . + REPLOC(21)*(R12R12*R21R21 . + R22R22*R11R11) . + REPLOC(22)*R12R22*R11R21 C C (PX PX | PY PZ): C REPIJ(5,9) = REPLOC(16)*R31R31*R32R33 . + REPLOC(17)*RADD01*R32R33 . + REPLOC(18)*RADD05*R31R31 . + REPLOC(19)*(R11R11*R12R13 . + R21R21*R22R23) . + REPLOC(20)*(R11R31*R12R33 . + R31R11*R32R13 . + R21R31*R22R33 . + R31R21*R32R23) . + REPLOC(21)*(R11R11*R22R23 . + R21R21*R12R13) . + REPLOC(22)*(R11R21*R12R23 . + R21R11*R22R13) C C (PY PZ | PX PX): C REPIJ(9,5) = REPLOC(16)*R32R33*R31R31 . + REPLOC(17)*RADD05*R31R31 . + REPLOC(18)*RADD01*R32R33 . + REPLOC(19)*(R12R13*R11R11 . + R22R23*R21R21) . + REPLOC(20)*(RADD13*R11R31 . + RADD14*R21R31) . + REPLOC(21)*(R12R13*R21R21 . + R22R23*R11R11) . + REPLOC(22)*RADD16*R11R21 C C (PX PX | PZ PZ): C REPIJ(5,10) = REPLOC(16)*R31R31*R33R33 . + REPLOC(17)*RADD01*R33R33 . + REPLOC(18)*RADD06*R31R31 . + REPLOC(19)*(R11R11*R13R13 . + R21R21*R23R23) . + REPLOC(20)*(R11R31*R13R33 . + R21R31*R23R33) . + REPLOC(21)*(R11R11*R23R23 . + R21R21*R13R13) . + REPLOC(22)*R11R21*R13R23 C C (PZ PZ | PX PX): C REPIJ(10,5) = REPLOC(16)*R33R33*R31R31 . + REPLOC(17)*RADD06*R31R31 . + REPLOC(18)*RADD01*R33R33 . + REPLOC(19)*(R13R13*R11R11 . + R23R23*R21R21) . + REPLOC(20)*(R13R33*R11R31 . + R23R33*R21R31) . + REPLOC(21)*(R13R13*R21R21 . + R23R23*R11R11) . + REPLOC(22)*R13R23*R11R21 C C (PX PY | PX PY): C REPIJ(6,6) = REPLOC(16)*R31R32*R31R32 . + REPLOC(17)*RADD02*R31R32 . + REPLOC(18)*RADD02*R31R32 . + REPLOC(19)*(R11R12*R11R12 . + R21R22*R21R22) . + REPLOC(20)*(RADD11**2 . + RADD12**2) . + REPLOC(21)*2.0D0*R11R12*R21R22 . + REPLOC(22)*RADD17**2 C C (PX PY | PX PZ): C REPIJ(6,7) = REPLOC(16)*R31R32*R31R33 . + REPLOC(17)*RADD02*R31R33 . + REPLOC(18)*R31R32*RADD03 . + REPLOC(19)*(R11R12*R11R13 . + R21R22*R21R23) . + REPLOC(20)*(RADD11*RADD09 . + RADD12*RADD10) . + REPLOC(21)*(R11R12*R21R23 . + R21R22*R11R13) . + REPLOC(22)*RADD17*(R11R23+R21R13) C C (PX PZ | PX PY): C REPIJ(7,6) = REPIJ(6,7) C C (PX PY | PY PY): C REPIJ(6,8) = REPLOC(16)*R31R32*R32R32 . + REPLOC(17)*RADD02*R32R32 . + REPLOC(18)*RADD04*R31R32 . + REPLOC(19)*(R11R12*R12R12 . + R21R22*R22R22) . + REPLOC(20)*(RADD11*R12R32 . + RADD12*R22R32) . + REPLOC(21)*(R11R12*R22R22 . + R21R22*R12R12) . + REPLOC(22)*RADD17*R12R22 C C (PY PY | PX PY): C REPIJ(8,6) = REPLOC(16)*R32R32*R31R32 . + REPLOC(17)*RADD04*R31R32 . + REPLOC(18)*RADD02*R32R32 . + REPLOC(19)*(R12R12*R11R12 . + R22R22*R21R22) . + REPLOC(20)*(R12R32*R11R32 . + R32R12*R31R12 . + R22R32*R21R32 . + R32R22*R31R22) . + REPLOC(21)*(R12R12*R21R22 . + R22R22*R11R12) . + REPLOC(22)*(R12R22*R11R22 . + R22R12*R21R12) C C (PX PY | PY PZ): C REPIJ(6,9) = REPLOC(16)*R31R32*R32R33 . + REPLOC(17)*RADD02*R32R33 . + REPLOC(18)*R31R32*RADD05 . + REPLOC(19)*(R11R12*R12R13 . + R21R22*R22R23) . + REPLOC(20)*(RADD11*RADD13 . + RADD12*RADD14) . + REPLOC(21)*(R11R12*R22R23 . + R21R22*R12R13) . + REPLOC(22)*RADD17*(R12R23+R22R13) C C (PY PZ | PX PY): C REPIJ(9,6) = REPIJ(6,9) C C (PX PY | PZ PZ): C REPIJ(6,10) = REPLOC(16)*R31R32*R33R33 . + REPLOC(17)*RADD02*R33R33 . + REPLOC(18)*RADD06*R31R32 . + REPLOC(19)*(R11R12*R13R13 . + R21R22*R23R23) . + REPLOC(20)*(RADD11*R13R33 . + RADD12*R23R33) . + REPLOC(21)*(R11R12*R23R23 . + R21R22*R13R13) . + REPLOC(22)*RADD17*R13R23 C C (PZ PZ | PX PY): C REPIJ(10,6) = REPLOC(16)*R33R33*R31R32 . + REPLOC(17)*RADD06*R31R32 . + REPLOC(18)*RADD02*R33R33 . + REPLOC(19)*(R13R13*R11R12 . + R23R23*R21R22) . + REPLOC(20)*(R13R33*R11R32 . + R33R13*R31R12 . + R23R33*R21R32 . + R33R23*R31R22) . + REPLOC(21)*(R13R13*R21R22 . + R23R23*R11R12) . + REPLOC(22)*(R13R23*R11R22 . + R23R13*R21R12) C C (PX PZ | PX PZ): C REPIJ(7,7) = REPLOC(16)*R31R33*R31R33 . + REPLOC(17)*RADD03*R31R33 . + REPLOC(18)*RADD03*R31R33 . + REPLOC(19)*(R11R13*R11R13 . + R21R23*R21R23) . + REPLOC(20)*(RADD09**2 . + RADD10**2) . + REPLOC(21)*R11R13*R21R23*2.0D0 . + REPLOC(22)*RADD15**2 C C (PX PZ | PY PY): C REPIJ(7,8) = REPLOC(16)*R31R33*R32R32 . + REPLOC(17)*RADD03*R32R32 . + REPLOC(18)*RADD04*R31R33 . + REPLOC(19)*(R11R13*R12R12 . + R21R23*R22R22) . + REPLOC(20)*(RADD09*R12R32 . + RADD10*R22R32) . + REPLOC(21)*(R11R13*R22R22 . + R21R23*R12R12) . + REPLOC(22)*RADD15*R12R22 C C (PY PY | PX PZ): C REPIJ(8,7) = REPLOC(16)*R32R32*R31R33 . + REPLOC(17)*RADD04*R31R33 . + REPLOC(18)*RADD03*R32R32 . + REPLOC(19)*(R12R12*R11R13 . + R22R22*R21R23) . + REPLOC(20)*(R12R32*R11R33 . + R32R12*R31R13 . + R22R32*R21R33 . + R32R22*R31R23) . + REPLOC(21)*(R12R12*R21R23 . + R22R22*R11R13) . + REPLOC(22)*(R12R22*R11R23 . + R22R12*R21R13) C C (PX PZ | PY PZ): C REPIJ(7,9) = REPLOC(16)*R31R33*R32R33 . + REPLOC(17)*RADD03*R32R33 . + REPLOC(18)*R31R33*RADD05 . + REPLOC(19)*(R11R13*R12R13 . + R21R23*R22R23) . + REPLOC(20)*(RADD09*RADD13 . + RADD10*R22R33 . + RADD10*R32R23) . + REPLOC(21)*(R11R13*R22R23 . + R21R23*R12R13) . + REPLOC(22)*RADD15*RADD16 C C (PY PZ | PX PZ): C REPIJ(9,7) = REPIJ(7,9) C C (PX PZ | PZ PZ): C REPIJ(7,10) = REPLOC(16)*R31R33*R33R33 . + REPLOC(17)*RADD03*R33R33 . + REPLOC(18)*RADD06*R31R33 . + REPLOC(19)*(R11R13*R13R13 . + R21R23*R23R23) . + REPLOC(20)*(RADD09*R13R33 . + RADD10*R23R33) . + REPLOC(21)*(R11R13*R23R23 . + R21R23*R13R13) . + REPLOC(22)*RADD15*R13R23 C C (PZ PZ | PX PZ): C REPIJ(10,7) = REPLOC(16)*R33R33*R31R33 . + REPLOC(17)*RADD06*R31R33 . + REPLOC(18)*RADD03*R33R33 . + REPLOC(19)*(R13R13*R11R13 . + R23R23*R21R23) . + REPLOC(20)*(R13R33*R11R33 . + R33R13*R31R13 . + R23R33*R21R33 . + R33R23*R31R23) . + REPLOC(21)*(R13R13*R21R23 . + R23R23*R11R13) . + REPLOC(22)*(R13R23*R11R23 . + R23R13*R21R13) C C (PY PY | PY PY): C REPIJ(8,8) = REPLOC(16)*R32R32*R32R32 . + REPLOC(17)*RADD04*R32R32 . + REPLOC(18)*RADD04*R32R32 . + REPLOC(19)*(R12R12*R12R12 . + R22R22*R22R22) . + REPLOC(20)*(R12R32*R12R32 . + R22R32*R22R32) . + REPLOC(21)*R12R12*R22R22*2.0D0 . + REPLOC(22)*R12R22*R12R22 C C (PY PY | PY PZ): C REPIJ(8,9) = REPLOC(16)*R32R32*R32R33 . + REPLOC(17)*RADD04*R32R33 . + REPLOC(18)*RADD05*R32R32 . + REPLOC(19)*(R12R12*R12R13 . + R22R22*R22R23) . + REPLOC(20)*(R12R32*R12R33 . + R32R12*R32R13 . + R22R32*R22R33 . + R32R22*R32R23) . + REPLOC(21)*(R12R12*R22R23 . + R22R22*R12R13) . + REPLOC(22)*(R12R22*R12R23 . + R22R12*R22R13) C C (PY PZ | PY PY): C REPIJ(9,8) = REPLOC(16)*R32R33*R32R32 . + REPLOC(17)*RADD05*R32R32 . + REPLOC(18)*RADD04*R32R33 . + REPLOC(19)*(R12R13*R12R12 . + R22R23*R22R22) . + REPLOC(20)*(RADD13*R12R32 . + RADD14*R22R32) . + REPLOC(21)*(R12R13*R22R22 . + R22R23*R12R12) . + REPLOC(22)*RADD16*R12R22 C C (PY PY | PZ PZ): C REPIJ(8,10) = REPLOC(16)*R32R32*R33R33 . + REPLOC(17)*RADD04*R33R33 . + REPLOC(18)*RADD06*R32R32 . + REPLOC(19)*(R12R12*R13R13 . + R22R22*R23R23) . + REPLOC(20)*(R12R32*R13R33 . + R22R32*R23R33) . + REPLOC(21)*(R12R12*R23R23 . + R22R22*R13R13) . + REPLOC(22)*R12R22*R13R23 C C (PZ PZ | PY PY): C REPIJ(10,8) = REPLOC(16)*R33R33*R32R32 . + REPLOC(17)*RADD06*R32R32 . + REPLOC(18)*RADD04*R33R33 . + REPLOC(19)*(R13R13*R12R12 . + R23R23*R22R22) . + REPLOC(20)*(R13R33*R12R32 . + R23R33*R22R32) . + REPLOC(21)*(R13R13*R22R22 . + R23R23*R12R12) . + REPLOC(22)*R13R23*R12R22 C C (PY PZ | PY PZ): C REPIJ(9,9) = REPLOC(16)*R32R33*R32R33 . + REPLOC(17)*RADD05*R32R33 . + REPLOC(18)*RADD05*R32R33 . + REPLOC(19)*(R12R13*R12R13 . + R22R23*R22R23) . + REPLOC(20)*(RADD13**2 . + RADD14**2) . + REPLOC(21)*(R12R13*R22R23)*2.0D0 . + REPLOC(22)*RADD16**2 C C (PY PZ | PZ PZ): C REPIJ(9,10) = REPLOC(16)*R32R33*R33R33 . + REPLOC(17)*RADD05*R33R33 . + REPLOC(18)*RADD06*R32R33 . + REPLOC(19)*(R12R13*R13R13 . + R22R23*R23R23) . + REPLOC(20)*(RADD13*R13R33 . + RADD14*R23R33) . + REPLOC(21)*(R12R13*R23R23 . + R22R23*R13R13) . + REPLOC(22)*RADD16*R13R23 C C (PZ PZ | PY PZ): C REPIJ(10,9) = REPLOC(16)*R33R33*R32R33 . + REPLOC(17)*RADD06*R32R33 . + REPLOC(18)*RADD05*R33R33 . + REPLOC(19)*(R13R13*R12R13 . + R23R23*R22R23) . + REPLOC(20)*(R13R33*R12R33 . + R33R13*R32R13 . + R23R33*R22R33 . + R33R23*R32R23) . + REPLOC(21)*(R13R13*R22R23 . + R23R23*R12R13) . + REPLOC(22)*(R13R23*R12R23 . + R23R13*R22R13) C C (PZ PZ | PZ PZ): C REPIJ(10,10) = REPLOC(16)*R33R33*R33R33 . + REPLOC(17)*RADD06*R33R33 . + REPLOC(18)*RADD06*R33R33 . + REPLOC(19)*(R13R13*R13R13 . + R23R23*R23R23) . + REPLOC(20)*(R13R33*R13R33 . + R23R33*R23R33) . + REPLOC(21)*R13R13*R23R23*2.0D0 . + REPLOC(22)*R13R23*R13R23 ENDIF RETURN END