mopac606_nbo/src/diat2.f

   1       SUBROUTINE DIAT2(NA,ESA,EPA,R12,NB,ESB,EPB,S)
   2       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   3       DIMENSION S(3,3,3)
   4 C***********************************************************************
   5 C
   6 C OVERLP CALCULATES OVERLAPS BETWEEN ATOMIC ORBITALS FOR PAIRS OF ATOMS
   7 C        IT CAN HANDLE THE ORBITALS 1S, 2S, 3S, 2P, AND 3P.
   8 C
   9 C***********************************************************************
  10       COMMON /SETC/ A(7),B(7),SA,SB,FACTOR,ISP,IPS
  11       DIMENSION INMB(17),III(78)
  12       SAVE INMB, III
  13       DATA INMB/1,0,2,2,3,4,5,6,7,0,8,8,8,9,10,11,12/
  14 C     NUMBERING CORRESPONDS TO BOND TYPE MATRIX GIVEN ABOVE
  15 C      THE CODE IS
  16 C
  17 C     III=1      FIRST - FIRST  ROW ELEMENTS
  18 C        =2      FIRST - SECOND
  19 C        =3      FIRST - THIRD
  20 C        =4      SECOND - SECOND
  21 C        =5      SECOND - THIRD
  22 C        =6      THIRD - THIRD
  23       DATA III/1,2,4,   2,4,4,   2,4,4,4,   2,4,4,4,4,
  24      1 2,4,4,4,4,4,   2,4,4,4,4,4,4,   3,5,5,5,5,5,5,6,
  25      2 3,5,5,5,5,5,5,6,6,   3,5,5,5,5,5,5,6,6,6,   3,5,5,5,5,5,5,6,6,6,6
  26      3, 3,5,5,5,5,5,5,6,6,6,6,6/
  27 C
  28 C      ASSIGNS BOND NUMBER
  29 C
  30       JMAX=MAX0(INMB(NA),INMB(NB))
  31       JMIN=MIN0(INMB(NA),INMB(NB))
  32       II=III((JMAX*(JMAX-1))/2+JMIN)
  33       DO 10 I=1,3
  34          DO 10 J=1,3
  35             DO 10 K=1,3
  36    10 S(I,J,K)=0.D0
  37       RAB=R12/0.529167D0
  38       GOTO (20,30,40,50,60,70), II
  39 C
  40 C     ------------------------------------------------------------------
  41 C *** THE ORDERING OF THE ELEMENTS WITHIN S IS
  42 C *** S(1,1,1)=(S(B)/S(A))
  43 C *** S(1,2,1)=(P-SIGMA(B)/S(A))
  44 C *** S(2,1,1)=(S(B)/P-SIGMA(A))
  45 C *** S(2,2,1)=(P-SIGMA(B)/P-SIGMA(A))
  46 C *** S(2,2,2)=(P-PI(B)/P-PI(A))
  47 C     ------------------------------------------------------------------
  48 C *** FIRST ROW - FIRST ROW OVERLAPS
  49 C
  50    20 CALL SET (ESA,ESB,NA,NB,RAB,II)
  51       S(1,1,1)=.25D00*SQRT((SA*SB*RAB*RAB)**3)*(A(3)*B(1)-B(3)*A(1))
  52       RETURN
  53 C
  54 C *** FIRST ROW - SECOND ROW OVERLAPS
  55 C
  56    30 CALL SET (ESA,ESB,NA,NB,RAB,II)
  57       W=SQRT((SA**3)*(SB**5))*(RAB**4)*0.125D00
  58       S(1,1,1) = SQRT(1.D00/3.D00)
  59       S(1,1,1)=W*S(1,1,1)*(A(4)*B(1)-B(4)*A(1)+A(3)*B(2)-B(3)*A(2))
  60       IF (NA.GT.1) CALL SET (EPA,ESB,NA,NB,RAB,II)
  61       IF (NB.GT.1) CALL SET (ESA,EPB,NA,NB,RAB,II)
  62       W=SQRT((SA**3)*(SB**5))*(RAB**4)*0.125D00
  63       S(ISP,IPS,1)=W*(A(3)*B(1)-B(3)*A(1)+A(4)*B(2)-B(4)*A(2))
  64       RETURN
  65 C
  66 C *** FIRST ROW - THIRD ROW OVERLAPS
  67 C
  68    40 CALL SET (ESA,ESB,NA,NB,RAB,II)
  69       W=SQRT((SA**3)*(SB**7)/7.5D00)*(RAB**5)*0.0625D00
  70       SROOT3 = SQRT(3.D00)
  71       S(1,1,1)=W*(A(5)*B(1)-B(5)*A(1)+
  72      12.D00*(A(4)*B(2)-B(4)*A(2)))/SROOT3
  73       IF (NA.GT.1) CALL SET (EPA,ESB,NA,NB,RAB,II)
  74       IF (NB.GT.1) CALL SET (ESA,EPB,NA,NB,RAB,II)
  75       W=SQRT((SA**3)*(SB**7)/7.5D00)*(RAB**5)*0.0625D00
  76       S(ISP,IPS,1)=W*(A(4)*(B(1)+B(3))-B(4)*(A(1)+A(3))+
  77      1B(2)*(A(3)+A(5))-A(2)*(B(3)+B(5)))
  78       RETURN
  79 C
  80 C *** SECOND ROW - SECOND ROW OVERLAPS
  81 C
  82    50 CALL SET (ESA,ESB,NA,NB,RAB,II)
  83       W=SQRT((SA*SB)**5)*(RAB**5)*0.0625D00
  84       RT3=1.D00/SQRT(3.D00)
  85       S(1,1,1)=W*(A(5)*B(1)+B(5)*A(1)-2.0D00*A(3)*B(3))/3.0D00
  86       CALL SET (ESA,EPB,NA,NB,RAB,II)
  87       IF (NA.GT.NB) CALL SET (EPA,ESB,NA,NB,RAB,II)
  88       W=SQRT((SA*SB)**5)*(RAB**5)*0.0625D00
  89       D=A(4)*(B(1)-B(3))-A(2)*(B(3)-B(5))
  90       E=B(4)*(A(1)-A(3))-B(2)*(A(3)-A(5))
  91       S(ISP,IPS,1)=W*RT3*(D+E)
  92       CALL SET (EPA,ESB,NA,NB,RAB,II)
  93       IF (NA.GT.NB) CALL SET (ESA,EPB,NA,NB,RAB,II)
  94       W=SQRT((SA*SB)**5)*(RAB**5)*0.0625D00
  95       D=A(4)*(B(1)-B(3))-A(2)*(B(3)-B(5))
  96       E=B(4)*(A(1)-A(3))-B(2)*(A(3)-A(5))
  97       S(IPS,ISP,1)=-W*RT3*(E-D)
  98       CALL SET (EPA,EPB,NA,NB,RAB,II)
  99       W=SQRT((SA*SB)**5)*(RAB**5)*0.0625D00
 100       S(2,2,1)=-W*(B(3)*(A(5)+A(1))-A(3)*(B(5)+B(1)))
 101       HD = .5D00
 102       S(2,2,2)=HD*W*(A(5)*(B(1)-B(3))-B(5)*(A(1)-A(3))
 103      1-A(3)*B(1)+B(3)*A(1))
 104       RETURN
 105 C
 106 C *** SECOND ROW - THIRD ROW OVERLAPS
 107 C
 108    60 CALL SET (ESA,ESB,NA,NB,RAB,II)
 109       W=SQRT((SA**5)*(SB**7)/7.5D00)*(RAB**6)*0.03125D00
 110       RT3 = 1.D00 / SQRT(3.D00)
 111       TD = 2.D00
 112       S(1,1,1)=W*(A(6)*B(1)+A(5)*B(2)-TD*(A(4)*B(3)+
 113      1A(3)*B(4))+A(2)*B(5)+A(
 114      21)*B(6))/3.D00
 115       CALL SET (ESA,EPB,NA,NB,RAB,II)
 116       IF (NA.GT.NB) CALL SET (EPA,ESB,NA,NB,RAB,II)
 117       W=SQRT((SA**5)*(SB**7)/7.5D00)*(RAB**6)*0.03125D00
 118       TD = 2.D00
 119       S(ISP,IPS,1)=W*RT3*(A(6)*B(2)+A(5)*B(1)-TD*(A(4)*B(4)+A(3)*B(3))
 120      1+A(2)*B(6)+A(1)*B(5))
 121       CALL SET (EPA,ESB,NA,NB,RAB,II)
 122       IF (NA.GT.NB) CALL SET (ESA,EPB,NA,NB,RAB,II)
 123       W=SQRT((SA**5)*SB**7/7.5D00)*(RAB**6)*0.03125D00
 124       TD = 2.D00
 125       S(IPS,ISP,1)=-W*RT3*(A(5)*(TD*B(3)-B(1))-B(5)*(TD*A(3)-A(1))-A(2
 126      1)*(B(6)-TD*B(4))+B(2)*(A(6)-TD*A(4)))
 127       CALL SET (EPA,EPB,NA,NB,RAB,II)
 128       W=SQRT((SA**5)*SB**7/7.5D00)*(RAB**6)*0.03125D00
 129       S(2,2,1)=-W*(B(4)*(A(1)+A(5))-A(4)*(B(1)+B(5))
 130      1+B(3)*(A(2)+A(6))-A(3)*(B(2)+B(6)))
 131       HD = .5D00
 132       S(2,2,2)=HD*W*(A(6)*(B(1)-B(3))-B(6)*(A(1)-
 133      1A(3))+A(5)*(B(2)-B(4))-B(5
 134      2)*(A(2)-A(4))-A(4)*B(1)+B(4)*A(1)-A(3)*B(2)+B(3)*A(2))
 135       RETURN
 136 C
 137 C *** THIRD ROW - THIRD ROW OVERLAPS
 138 C
 139    70 CALL SET (ESA,ESB,NA,NB,RAB,II)
 140       W=SQRT((SA*SB*RAB*RAB)**7)/480.D00
 141       RT3 = 1.D00 / SQRT(3.D00)
 142       S(1,1,1)=W*(A(7)*B(1)-3.D00*(A(5)*B(3)-A(3)*B(5))-A(1)*B(7))/3.D00
 143       CALL SET (ESA,EPB,NA,NB,RAB,II)
 144       IF (NA.GT.NB) CALL SET (EPA,ESB,NA,NB,RAB,II)
 145       W=SQRT((SA*SB*RAB*RAB)**7)/480.D00
 146       D=A(6)*(B(1)-B(3))-2.D00*A(4)*(B(3)-B(5))+A(2)*(B(5)-B(7))
 147       E=B(6)*(A(1)-A(3))-2.D00*B(4)*(A(3)-A(5))+B(2)*(A(5)-A(7))
 148       S(ISP,IPS,1)=W*RT3*(D-E)
 149       CALL SET (EPA,ESB,NA,NB,RAB,II)
 150       IF (NA.GT.NB) CALL SET (ESA,EPB,NA,NB,RAB,II)
 151       W=SQRT((SA*SB*RAB*RAB)**7)/480.D00
 152       D=A(6)*(B(1)-B(3))-2.D00*A(4)*(B(3)-B(5))+A(2)*(B(5)-B(7))
 153       E=B(6)*(A(1)-A(3))-2.D00*B(4)*(A(3)-A(5))+B(2)*(A(5)-A(7))
 154       S(IPS,ISP,1)=-W*RT3*(-D-E)
 155       CALL SET (EPA,EPB,NA,NB,RAB,II)
 156       W=SQRT((SA*SB*RAB*RAB)**7)/480.D00
 157       TD = 2.D00
 158       S(2,2,1)=-W*(A(3)*(B(7)+TD*B(3))-A(5)*(B(1)+
 159      1TD*B(5))-B(5)*A(1)+A(7)*B(3))
 160       HD = .5D00
 161       S(2,2,2)=HD*W*(A(7)*(B(1)-B(3))+B(7)*(A(1)-
 162      1A(3))+A(5)*(B(5)-B(3)-B(1)
 163      2)+B(5)*(A(5)-A(3)-A(1))+2.D00*A(3)*B(3))
 164       RETURN
 165 C
 166       END
 167       SUBROUTINE SET (S1,S2,NA,NB,RAB,II)
 168       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 169       COMMON /SETC/ A(7),B(7),SA,SB,FACTOR,ISP,IPS
 170 C***********************************************************************
 171 C
 172 C     SET IS PART OF THE OVERLAP CALCULATION, CALLED BY OVERLP.
 173 C         IT CALLS AINTGS AND BINTGS
 174 C
 175 C***********************************************************************
 176       IF (NA.GT.NB) GO TO 10
 177       ISP=1
 178       IPS=2
 179       SA=S1
 180       SB=S2
 181       GOTO 20
 182    10 ISP=2
 183       IPS=1
 184       SA=S2
 185       SB=S1
 186    20 J=II+2
 187       IF (II.GT.3) J=J-1
 188       ALPHA=0.5D00*RAB*(SA+SB)
 189       BETA=0.5D00*RAB*(SB-SA)
 190       JCALL=J-1
 191       CALL AINTGS (ALPHA,JCALL)
 192       CALL BINTGS (BETA,JCALL)
 193       RETURN
 194 C
 195       END
 196       SUBROUTINE AINTGS (X,K)
 197       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 198       COMMON /SETC/ A(7),B(7),SDUM(3),IDUM(2)
 199 C***********************************************************************
 200 C
 201 C    AINTGS FORMS THE "A" INTEGRALS FOR THE OVERLAP CALCULATION.
 202 C
 203 C***********************************************************************
 204       C=EXP(-X)
 205       A(1)=C/X
 206       DO 10 I=1,K
 207          A(I+1)=(A(I)*I+C)/X
 208    10 CONTINUE
 209       RETURN
 210 C
 211       END
 212       SUBROUTINE BINTGS (X,K)
 213       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 214       COMMON /SETC/ A(7),B(7),SDUM(3),IDUM(2)
 215       DIMENSION FACT(17)
 216 C**********************************************************************
 217 C
 218 C     BINTGS FORMS THE "B" INTEGRALS FOR THE OVERLAP CALCULATION.
 219 C
 220 C**********************************************************************
 221       SAVE FACT
 222       DATA FACT/1.D0,2.D0,6.D0,24.D0,120.D0,720.D0,5040.D0,40320.D0,
 223      1362880.D0,3628800.D0,39916800.D0,479001600.D0,6227020800.D0,
 224      28.71782912D10,1.307674368D12,2.092278989D13,3.556874281D14/
 225       IO=0
 226       ABSX = ABS(X)
 227       IF (ABSX.GT.3.D00) GO TO 40
 228       IF (ABSX.LE.2.D00) GO TO 10
 229       IF (K.LE.10) GO TO 40
 230       LAST=15
 231       GO TO 60
 232    10 IF (ABSX.LE.1.D00) GO TO 20
 233       IF (K.LE.7) GO TO 40
 234       LAST=12
 235       GO TO 60
 236    20 IF (ABSX.LE.0.5D00) GO TO 30
 237       IF (K.LE.5) GO TO 40
 238       LAST=7
 239       GO TO 60
 240    30 IF (ABSX.LE.1.D-6) GOTO 90
 241       LAST=6
 242       GO TO 60
 243    40 EXPX=EXP(X)
 244       EXPMX=1.D00/EXPX
 245       B(1)=(EXPX-EXPMX)/X
 246       DO 50 I=1,K
 247    50 B(I+1)=(I*B(I)+(-1.D00)**I*EXPX-EXPMX)/X
 248       GO TO 110
 249    60 DO 80 I=IO,K
 250          Y=0.0D00
 251          DO 70 M=IO,LAST
 252             XF=1.0D00
 253             IF(M.NE.0) XF=FACT(M)
 254    70    Y=Y+(-X)**M*(2*MOD(M+I+1,2))/(XF*(M+I+1))
 255    80 B(I+1)=Y
 256       GO TO 110
 257    90 DO 100 I=IO,K
 258   100 B(I+1)=(2*MOD(I+1,2))/(I+1.D0)
 259   110 CONTINUE
 260       RETURN
 261 C
 262       END