PageRenderTime 47ms CodeModel.GetById 18ms RepoModel.GetById 1ms app.codeStats 0ms

/TwoE/INT/IntB10100101.f90

#
FORTRAN Modern | 243 lines | 189 code | 0 blank | 54 comment | 0 complexity | 9e7f744c2bc121b5217203f51f5b6c2f MD5 | raw file
Possible License(s): GPL-2.0, GPL-3.0, BSD-3-Clause 
    

  1. !------------------------------------------------------------------------------
    2. 

  2. !    This code is part of the MondoSCF suite of programs for linear scaling
    3. 

  3. !    electronic structure theory and ab initio molecular dynamics.
    4. 

  4. !
    5. 

  5. !    Copyright (2004). The Regents of the University of California. This
    6. 

  6. !    material was produced under U.S. Government contract W-7405-ENG-36
    7. 

  7. !    for Los Alamos National Laboratory, which is operated by the University
    8. 

  8. !    of California for the U.S. Department of Energy. The U.S. Government has
    9. 

  9. !    rights to use, reproduce, and distribute this software.  NEITHER THE
    10. 

  10. !    GOVERNMENT NOR THE UNIVERSITY MAKES ANY WARRANTY, EXPRESS OR IMPLIED,
    11. 

  11. !    OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
    12. 

  12. !
    13. 

  13. !    This program is free software; you can redistribute it and/or modify
    14. 

  14. !    it under the terms of the GNU General Public License as published by the
    15. 

  15. !    Free Software Foundation; either version 2 of the License, or (at your
    16. 

  16. !    option) any later version. Accordingly, this program is distributed in
    17. 

  17. !    the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
    18. 

  18. !    the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
    19. 

  19. !    PURPOSE. See the GNU General Public License at www.gnu.org for details.
    20. 

  20. !
    21. 

  21. !    While you may do as you like with this software, the GNU license requires
    22. 

  22. !    that you clearly mark derivative software.  In addition, you are encouraged
    23. 

  23. !    to return derivative works to the MondoSCF group for review, and possible
    24. 

  24. !    disemination in future releases.
    25. 

  25. !------------------------------------------------------------------------------
    26. 

  26. ! ---------------------------------------------------------- 
    27. 

  27. ! COMPUTES THE INTEGRAL CLASS (f f|s s) 
    28. 

  28. ! ---------------------------------------------------------- 
    29. 

  29.    SUBROUTINE IntB10100101(PrmBufB,LBra,PrmBufK,LKet,ACInfo,BDInfo, & 
    30. 

  30.                               OA,LDA,OB,LDB,OC,LDC,OD,LDD,PBC,INTGRL) 
    31. 

  31.       USE DerivedTypes
    32. 

  32.       USE VScratchB
    33. 

  33.       USE GlobalScalars
    34. 

  34.       USE ShellPairStruct
    35. 

  35.       USE GammaF6
    36. 

  36.       IMPLICIT REAL(DOUBLE) (W)
    37. 

  37.       INTEGER        :: LBra,LKet,CDOffSet
    38. 

  38.       REAL(DOUBLE)   :: PrmBufB(10,LBra),PrmBufK(10,LKet)
    39. 

  39.       TYPE(SmallAtomInfo) :: ACInfo,BDInfo
    40. 

  40.       TYPE(PBCInfo) :: PBC
    41. 

  41.       REAL(DOUBLE)  :: INTGRL(*)
    42. 

  42.       REAL(DOUBLE)  :: Ax,Ay,Az,Bx,By,Bz,Cx,Cy,Cz
    43. 

  43.       REAL(DOUBLE)  :: Dx,Dy,Dz,Qx,Qy,Qz,Px,Py,Pz
    44. 

  44.       REAL(DOUBLE)  :: PQx,PQy,PQz,FPQx,FPQy,FPQz
    45. 

  45.       REAL(DOUBLE)  :: Zeta,Eta,Omega,Up,Uq,Upq
    46. 

  46.       REAL(DOUBLE)  :: T,ET,TwoT,InvT,SqInvT
    47. 

  47.       REAL(DOUBLE)  :: VRR(84,1,0:6)
    48. 

  48.       REAL(DOUBLE)  :: HRR(84,1,1)
    49. 

  49.       INTEGER       :: OffSet,OA,LDA,OB,LDB,OC,LDC,OD,LDD,I,J,K,L
    50. 

  50.       EXTERNAL InitDbl
    51. 

  51.       CALL InitDbl(84*1,HRR(1,1,1))
    52. 

  52.       Ax=ACInfo%Atm1X
    53. 

  53.       Ay=ACInfo%Atm1Y
    54. 

  54.       Az=ACInfo%Atm1Z
    55. 

  55.       Bx=ACInfo%Atm2X
    56. 

  56.       By=ACInfo%Atm2Y
    57. 

  57.       Bz=ACInfo%Atm2Z
    58. 

  58.       Cx=BDInfo%Atm1X
    59. 

  59.       Cy=BDInfo%Atm1Y
    60. 

  60.       Cz=BDInfo%Atm1Z
    61. 

  61.       Dx=BDInfo%Atm2X
    62. 

  62.       Dy=BDInfo%Atm2Y
    63. 

  63.       Dz=BDInfo%Atm2Z
    64. 

  64.       ABx=Ax-Bx
    65. 

  65.       ABy=Ay-By
    66. 

  66.       ABz=Az-Bz
    67. 

  67.       CDx=Cx-Dx
    68. 

  68.       CDy=Cy-Dy
    69. 

  69.       CDz=Cz-Dz
    70. 

  70.       DO J=1,LKet ! K^2 VRR |N0) loop 
    71. 

  71.          Eta=PrmBufK(1,J)
    72. 

  72.          Qx=PrmBufK(2,J)
    73. 

  73.          Qy=PrmBufK(3,J)
    74. 

  74.          Qz=PrmBufK(4,J)
    75. 

  75.          Uq=PrmBufK(5,J)
    76. 

  76.          QCx=Qx-Cx
    77. 

  77.          QCy=Qy-Cy
    78. 

  78.          QCz=Qz-Cz
    79. 

  79.          DO K=1,LBra ! K^2 VRR (M0| loop 
    80. 

  80.             Zeta=PrmBufB(1,K)
    81. 

  81.             Px=PrmBufB(2,K)
    82. 

  82.             Py=PrmBufB(3,K)
    83. 

  83.             Pz=PrmBufB(4,K)
    84. 

  84.             Up=PrmBufB(5,K)
    85. 

  85.             r1xZpE=One/(Zeta+Eta)
    86. 

  86.             Upq=SQRT(r1xZpE)*Up*Uq
    87. 

  87.             HfxZpE=Half/(Zeta+Eta)
    88. 

  88.             r1x2E=Half/Eta
    89. 

  89.             r1x2Z=Half/Zeta
    90. 

  90.             ExZpE=Eta*r1xZpE
    91. 

  91.             ZxZpE=Zeta*r1xZpE
    92. 

  92.             Omega=Eta*Zeta*r1xZpE
    93. 

  93.             PAx=Px-Ax
    94. 

  94.             PAy=Py-Ay
    95. 

  95.             PAz=Pz-Az
    96. 

  96.             PQx=Px-Qx
    97. 

  97.             PQy=Py-Qy
    98. 

  98.             PQz=Pz-Qz
    99. 

  99.             ! Begin Minimum Image Convention 
    100. 

  100.             FPQx = PQx*PBC%InvBoxSh%D(1,1)+PQy*PBC%InvBoxSh%D(1,2)+PQz*PBC%InvBoxSh%D(1,3)
    101. 

  101.             FPQy = PQy*PBC%InvBoxSh%D(2,2)+PQz*PBC%InvBoxSh%D(2,3)
    102. 

  102.             FPQz = PQz*PBC%InvBoxSh%D(3,3)
    103. 

  103.             IF(PBC%AutoW%I(1)==1)FPQx=FPQx-ANINT(FPQx-SIGN(1D-15,FPQx))
    104. 

  104.             IF(PBC%AutoW%I(2)==1)FPQy=FPQy-ANINT(FPQy-SIGN(1D-15,FPQy))
    105. 

  105.             IF(PBC%AutoW%I(3)==1)FPQz=FPQz-ANINT(FPQz-SIGN(1D-15,FPQz))
    106. 

  106.             PQx=FPQx*PBC%BoxShape%D(1,1)+FPQy*PBC%BoxShape%D(1,2)+FPQz*PBC%BoxShape%D(1,3)
    107. 

  107.             PQy=FPQy*PBC%BoxShape%D(2,2)+FPQz*PBC%BoxShape%D(2,3)
    108. 

  108.             PQz=FPQz*PBC%BoxShape%D(3,3)
    109. 

  109.             ! End MIC
    110. 

  110.             WPx = -Eta*PQx*r1xZpE
    111. 

  111.             WPy = -Eta*PQy*r1xZpE
    112. 

  112.             WPz = -Eta*PQz*r1xZpE
    113. 

  113.             WQx = Zeta*PQx*r1xZpE
    114. 

  114.             WQy = Zeta*PQy*r1xZpE
    115. 

  115.             WQz = Zeta*PQz*r1xZpE
    116. 

  116.             T=Omega*(PQx*PQx+PQy*PQy+PQz*PQz)
    117. 

  117.             IF(T<Gamma_Switch)THEN
    118. 

  118.               L=AINT(T*Gamma_Grid)
    119. 

  119.               ET=EXP(-T)
    120. 

  120.               TwoT=Two*T
    121. 

  121.               W6=(F6_0(L)+T*(F6_1(L)+T*(F6_2(L)+T*(F6_3(L)+T*F6_4(L)))))
    122. 

  122.               W5=+9.090909090909090D-02*(TwoT*W6+ET)
    123. 

  123.               W4=+1.111111111111111D-01*(TwoT*W5+ET)
    124. 

  124.               W3=+1.428571428571428D-01*(TwoT*W4+ET)
    125. 

  125.               W2=+2.000000000000000D-01*(TwoT*W3+ET)
    126. 

  126.               W1=+3.333333333333333D-01*(TwoT*W2+ET)
    127. 

  127.               W0=TwoT*W1+ET
    128. 

  128.               VRR(1,1,0)=Upq*W0
    129. 

  129.               VRR(1,1,1)=Upq*W1
    130. 

  130.               VRR(1,1,2)=Upq*W2
    131. 

  131.               VRR(1,1,3)=Upq*W3
    132. 

  132.               VRR(1,1,4)=Upq*W4
    133. 

  133.               VRR(1,1,5)=Upq*W5
    134. 

  134.               VRR(1,1,6)=Upq*W6
    135. 

  135.             ELSE
    136. 

  136.               InvT=One/T
    137. 

  137.               SqInvT=DSQRT(InvT)
    138. 

  138.               VRR(1,1,0)=+8.862269254527580D-01*Upq*SqInvT
    139. 

  139.               SqInvT=SqInvT*InvT
    140. 

  140.               VRR(1,1,1)=+4.431134627263790D-01*Upq*SqInvT
    141. 

  141.               SqInvT=SqInvT*InvT
    142. 

  142.               VRR(1,1,2)=+6.646701940895685D-01*Upq*SqInvT
    143. 

  143.               SqInvT=SqInvT*InvT
    144. 

  144.               VRR(1,1,3)=+1.661675485223921D+00*Upq*SqInvT
    145. 

  145.               SqInvT=SqInvT*InvT
    146. 

  146.               VRR(1,1,4)=+5.815864198283724D+00*Upq*SqInvT
    147. 

  147.               SqInvT=SqInvT*InvT
    148. 

  148.               VRR(1,1,5)=+2.617138889227676D+01*Upq*SqInvT
    149. 

  149.               SqInvT=SqInvT*InvT
    150. 

  150.               VRR(1,1,6)=+1.439426389075222D+02*Upq*SqInvT
    151. 

  151.             ENDIF
    152. 

  152.             ! Generating (p0|s0)^(5)
    153. 

  153.             VRR(2,1,5)=PAx*VRR(1,1,5)+WPx*VRR(1,1,6) 
    154. 

  154.             VRR(3,1,5)=PAy*VRR(1,1,5)+WPy*VRR(1,1,6) 
    155. 

  155.             VRR(4,1,5)=PAz*VRR(1,1,5)+WPz*VRR(1,1,6) 
    156. 

  156.             ! Generating (p0|s0)^(4)
    157. 

  157.             VRR(2,1,4)=PAx*VRR(1,1,4)+WPx*VRR(1,1,5) 
    158. 

  158.             VRR(3,1,4)=PAy*VRR(1,1,4)+WPy*VRR(1,1,5) 
    159. 

  159.             VRR(4,1,4)=PAz*VRR(1,1,4)+WPz*VRR(1,1,5) 
    160. 

  160.             ! Generating (p0|s0)^(3)
    161. 

  161.             VRR(2,1,3)=PAx*VRR(1,1,3)+WPx*VRR(1,1,4) 
    162. 

  162.             VRR(3,1,3)=PAy*VRR(1,1,3)+WPy*VRR(1,1,4) 
    163. 

  163.             VRR(4,1,3)=PAz*VRR(1,1,3)+WPz*VRR(1,1,4) 
    164. 

  164.             ! Generating (p0|s0)^(2)
    165. 

  165.             VRR(2,1,2)=PAx*VRR(1,1,2)+WPx*VRR(1,1,3) 
    166. 

  166.             VRR(3,1,2)=PAy*VRR(1,1,2)+WPy*VRR(1,1,3) 
    167. 

  167.             VRR(4,1,2)=PAz*VRR(1,1,2)+WPz*VRR(1,1,3) 
    168. 

  168.             ! Generating (p0|s0)^(1)
    169. 

  169.             VRR(2,1,1)=PAx*VRR(1,1,1)+WPx*VRR(1,1,2) 
    170. 

  170.             VRR(3,1,1)=PAy*VRR(1,1,1)+WPy*VRR(1,1,2) 
    171. 

  171.             VRR(4,1,1)=PAz*VRR(1,1,1)+WPz*VRR(1,1,2) 
    172. 

  172.             ! Generating (p0|s0)^(0)
    173. 

  173.             VRR(2,1,0)=PAx*VRR(1,1,0)+WPx*VRR(1,1,1) 
    174. 

  174.             VRR(3,1,0)=PAy*VRR(1,1,0)+WPy*VRR(1,1,1) 
    175. 

  175.             VRR(4,1,0)=PAz*VRR(1,1,0)+WPz*VRR(1,1,1) 
    176. 

  176.             ! Generating (d0|s0)^(4)
    177. 

  177.             VRR(5,1,4)=PAx*VRR(2,1,4)+r1x2Z*(VRR(1,1,4)-ExZpE*VRR(1,1,5))+WPx*VRR(2,1,5)
    178. 

  178.             VRR(6,1,4)=PAx*VRR(3,1,4)+WPx*VRR(3,1,5)
    179. 

  179.             VRR(7,1,4)=PAy*VRR(3,1,4)+r1x2Z*(VRR(1,1,4)-ExZpE*VRR(1,1,5))+WPy*VRR(3,1,5)
    180. 

  180.             VRR(8,1,4)=PAx*VRR(4,1,4)+WPx*VRR(4,1,5)
    181. 

  181.             VRR(9,1,4)=PAy*VRR(4,1,4)+WPy*VRR(4,1,5)
    182. 

  182.             VRR(10,1,4)=PAz*VRR(4,1,4)+r1x2Z*(VRR(1,1,4)-ExZpE*VRR(1,1,5))+WPz*VRR(4,1,5)
    183. 

  183.             ! Generating (d0|s0)^(3)
    184. 

  184.             VRR(5,1,3)=PAx*VRR(2,1,3)+r1x2Z*(VRR(1,1,3)-ExZpE*VRR(1,1,4))+WPx*VRR(2,1,4)
    185. 

  185.             VRR(6,1,3)=PAx*VRR(3,1,3)+WPx*VRR(3,1,4)
    186. 

  186.             VRR(7,1,3)=PAy*VRR(3,1,3)+r1x2Z*(VRR(1,1,3)-ExZpE*VRR(1,1,4))+WPy*VRR(3,1,4)
    187. 

  187.             VRR(8,1,3)=PAx*VRR(4,1,3)+WPx*VRR(4,1,4)
    188. 

  188.             VRR(9,1,3)=PAy*VRR(4,1,3)+WPy*VRR(4,1,4)
    189. 

  189.             VRR(10,1,3)=PAz*VRR(4,1,3)+r1x2Z*(VRR(1,1,3)-ExZpE*VRR(1,1,4))+WPz*VRR(4,1,4)
    190. 

  190.             ! Generating (d0|s0)^(2)
    191. 

  191.             VRR(5,1,2)=PAx*VRR(2,1,2)+r1x2Z*(VRR(1,1,2)-ExZpE*VRR(1,1,3))+WPx*VRR(2,1,3)
    192. 

  192.             VRR(6,1,2)=PAx*VRR(3,1,2)+WPx*VRR(3,1,3)
    193. 

  193.             VRR(7,1,2)=PAy*VRR(3,1,2)+r1x2Z*(VRR(1,1,2)-ExZpE*VRR(1,1,3))+WPy*VRR(3,1,3)
    194. 

  194.             VRR(8,1,2)=PAx*VRR(4,1,2)+WPx*VRR(4,1,3)
    195. 

  195.             VRR(9,1,2)=PAy*VRR(4,1,2)+WPy*VRR(4,1,3)
    196. 

  196.             VRR(10,1,2)=PAz*VRR(4,1,2)+r1x2Z*(VRR(1,1,2)-ExZpE*VRR(1,1,3))+WPz*VRR(4,1,3)
    197. 

  197.             ! Generating (d0|s0)^(1)
    198. 

  198.             VRR(5,1,1)=PAx*VRR(2,1,1)+r1x2Z*(VRR(1,1,1)-ExZpE*VRR(1,1,2))+WPx*VRR(2,1,2)
    199. 

  199.             VRR(6,1,1)=PAx*VRR(3,1,1)+WPx*VRR(3,1,2)
    200. 

  200.             VRR(7,1,1)=PAy*VRR(3,1,1)+r1x2Z*(VRR(1,1,1)-ExZpE*VRR(1,1,2))+WPy*VRR(3,1,2)
    201. 

  201.             VRR(8,1,1)=PAx*VRR(4,1,1)+WPx*VRR(4,1,2)
    202. 

  202.             VRR(9,1,1)=PAy*VRR(4,1,1)+WPy*VRR(4,1,2)
    203. 

  203.             VRR(10,1,1)=PAz*VRR(4,1,1)+r1x2Z*(VRR(1,1,1)-ExZpE*VRR(1,1,2))+WPz*VRR(4,1,2)
    204. 

  204.             ! Generating (d0|s0)^(0)
    205. 

  205.             VRR(5,1,0)=PAx*VRR(2,1,0)+r1x2Z*(VRR(1,1,0)-ExZpE*VRR(1,1,1))+WPx*VRR(2,1,1)
    206. 

  206.             VRR(6,1,0)=PAx*VRR(3,1,0)+WPx*VRR(3,1,1)
    207. 

  207.             VRR(7,1,0)=PAy*VRR(3,1,0)+r1x2Z*(VRR(1,1,0)-ExZpE*VRR(1,1,1))+WPy*VRR(3,1,1)
    208. 

  208.             VRR(8,1,0)=PAx*VRR(4,1,0)+WPx*VRR(4,1,1)
    209. 

  209.             VRR(9,1,0)=PAy*VRR(4,1,0)+WPy*VRR(4,1,1)
    210. 

  210.             VRR(10,1,0)=PAz*VRR(4,1,0)+r1x2Z*(VRR(1,1,0)-ExZpE*VRR(1,1,1))+WPz*VRR(4,1,1)
    211. 

  211.             ! Generating (f0|s0)^(3)
    212. 

  212.             CALL VRRf0s0(84,1,VRR(1,1,3),VRR(1,1,4))
    213. 

  213.             ! Generating (f0|s0)^(2)
    214. 

  214.             CALL VRRf0s0(84,1,VRR(1,1,2),VRR(1,1,3))
    215. 

  215.             ! Generating (f0|s0)^(1)
    216. 

  216.             CALL VRRf0s0(84,1,VRR(1,1,1),VRR(1,1,2))
    217. 

  217.             ! Generating (f0|s0)^(0)
    218. 

  218.             CALL VRRf0s0(84,1,VRR(1,1,0),VRR(1,1,1))
    219. 

  219.             ! Generating (g0|s0)^(2)
    220. 

  220.             CALL VRRg0s0(84,1,VRR(1,1,2),VRR(1,1,3))
    221. 

  221.             ! Generating (g0|s0)^(1)
    222. 

  222.             CALL VRRg0s0(84,1,VRR(1,1,1),VRR(1,1,2))
    223. 

  223.             ! Generating (g0|s0)^(0)
    224. 

  224.             CALL VRRg0s0(84,1,VRR(1,1,0),VRR(1,1,1))
    225. 

  225.             ! Generating (h0|s0)^(1)
    226. 

  226.             CALL VRRh0s0(84,1,VRR(1,1,1),VRR(1,1,2))
    227. 

  227.             ! Generating (h0|s0)^(0)
    228. 

  228.             CALL VRRh0s0(84,1,VRR(1,1,0),VRR(1,1,1))
    229. 

  229.             ! Generating (i0|s0)^(0)
    230. 

  230.             CALL VRRi0s0(84,1,VRR(1,1,0),VRR(1,1,1))
    231. 

  231.             ! Contracting ... 
    232. 

  232.             CALL DBLAXPY(84,HRR(1,1,1),VRR(1,1,0)) 
    233. 

  233.          ENDDO ! (M0| loop
    234. 

  234.       ENDDO ! |N0) loop
    235. 

  235.       ! No need to generate (f,0|s,s)^(0) 
    236. 

  236.       ! Generating (f,f|s,s)^(0) 
    237. 

  237.       DO L=1,1
    238. 

  238.          DO K=1,1
    239. 

  239.             CDOffSet=(OC+K-1)*LDC+(OD+L-1)*LDD 
    240. 

  240.             CALL BraHRR1010(OA,OB,LDA,LDB,CDOffSet,HRR(1,K,L),INTGRL) 
    241. 

  241.           ENDDO 
    242. 

  242.       ENDDO 
    243. 

  243.     END SUBROUTINE IntB10100101
    244. 

  244.