1C Copyright (c) 2003-2010 University of Florida 2C 3C This program is free software; you can redistribute it and/or modify 4C it under the terms of the GNU General Public License as published by 5C the Free Software Foundation; either version 2 of the License, or 6C (at your option) any later version. 7 8C This program is distributed in the hope that it will be useful, 9C but WITHOUT ANY WARRANTY; without even the implied warranty of 10C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11C GNU General Public License for more details. 12 13C The GNU General Public License is included in this distribution 14C in the file COPYRIGHT. 15 SUBROUTINE ERD__2D_ATOM_COEFFICIENTS 16 + 17 + ( MIJ,MKL,MIJKL, 18 + NGQP,MGQIJKL, 19 + P, 20 + PINVHF,QINVHF,PQPINV, 21 + RTS, 22 + CASE2D, 23 + 24 + B00,B01,B10 ) 25 + 26C------------------------------------------------------------------------ 27C OPERATION : ERD__2D_ATOM_COEFFICIENTS 28C MODULE : ELECTRON REPULSION INTEGRALS DIRECT 29C MODULE-ID : ERD 30C SUBROUTINES : none 31C DESCRIPTION : This operation evaluates the atomic VRR coefficients 32C for the 2D integrals for the present set of NGQP roots 33C corresponding to all i,j,k,l exponent quadruplets. 34C 35C 36C Input: 37C 38C MIJ(KL) = current # of ij (kl) primitive 39C index pairs corresponding to 40C the csh pairs A,B (C,D) 41C MIJKL = current # of ijkl primitive 42C index quadruplets (= MIJ*MKL) 43C NGQP = # of gaussian quadrature points 44C (roots) 45C MGQIJKL = # of roots times # of ijkl 46C quadruplets (= NGQP*MIJKL) 47C P = current MIJ exponent sums for 48C csh A and B 49C P(Q)INVHF = current MIJ (MKL) values of 50C 1/(2*P(Q)), where P and Q are 51C the corresponding exponent sums 52C for csh A and B (C and D) 53C PQPINV = current MIJKL values of 1/(P+Q) 54C RTS = current MGQIJKL values of all 55C quadrature roots 56C CASE2D = integer value within the range 57C from 1 to 9, indicating which 58C 2D coefficient evaluation case 59C is present to trigger specific 60C simplified sections of the code 61C 62C Output: 63C 64C Bxx = the coordinate independent atomic 65C B-coefficients (xx=00,01,10) 66C 67C 68C AUTHOR : Norbert Flocke 69C------------------------------------------------------------------------ 70C 71C 72C ...include files and declare variables. 73C 74C 75 IMPLICIT NONE 76 77 INTEGER CASE2D 78 INTEGER IJ,KL 79 INTEGER M,N 80 INTEGER MIJ,MKL,MIJKL 81 INTEGER NG,NGQP,MGQIJKL 82 83 DOUBLE PRECISION PIJ 84 DOUBLE PRECISION PSCALE,QSCALE 85 DOUBLE PRECISION ROOT 86 DOUBLE PRECISION TWOP,TWOQ,TWOPQ 87 DOUBLE PRECISION HALF,ONE 88 89 DOUBLE PRECISION P (1:MIJ) 90 DOUBLE PRECISION PINVHF (1:MIJ) 91 DOUBLE PRECISION QINVHF (1:MKL) 92 93 DOUBLE PRECISION PQPINV (1:MIJKL) 94 95 DOUBLE PRECISION B00 (1:MGQIJKL) 96 DOUBLE PRECISION B01 (1:MGQIJKL) 97 DOUBLE PRECISION B10 (1:MGQIJKL) 98 DOUBLE PRECISION RTS (1:MGQIJKL) 99 100 PARAMETER (HALF = 0.5D0) 101 PARAMETER (ONE = 1.D0) 102C 103C 104C------------------------------------------------------------------------ 105C 106C 107C ...jump according to the 9 different cases that can arise: 108C 109C P-shell = s- ,p- or higher angular momentum 110C Q-shell = s- ,p- or higher angular momentum 111C 112C each leading to simplifications in the VRR formulas. 113C The case present has been evaluated outside this 114C routine and is transmitted via argument. 115C 116C 117 GOTO (1,1,5,1,4,7,6,8,9) CASE2D 118C 119C 120C ...the cases: P = s-shell and Q = s-shell 121C P = p-shell and Q = s-shell 122C P = s-shell and Q = p-shell 123C (no coefficients here) 124C 125C 126 1 RETURN 127C 128C 129C ...the case P = p-shell and Q = p-shell. 130C (here we know that NGQP = 2) 131C 132C 133 4 N = 1 134 DO 400 M = 1,MIJKL 135 TWOPQ = HALF * PQPINV (M) 136 B00 (N) = RTS (N) * TWOPQ 137 B00 (N+1) = RTS (N+1) * TWOPQ 138 N = N + 2 139 400 CONTINUE 140 141 RETURN 142C 143C 144C ...the case P > p-shell and Q = s-shell. 145C 146C 147 5 M = 0 148 N = 0 149 DO 500 IJ = 1,MIJ 150 PIJ = P (IJ) 151 TWOP = PINVHF (IJ) 152 DO 510 KL = 1,MKL 153 M = M + 1 154 QSCALE = ONE - PIJ * PQPINV (M) 155 156 DO 520 NG = 1,NGQP 157 N = N + 1 158 B10 (N)= (ONE - QSCALE * RTS (N)) * TWOP 159 520 CONTINUE 160 161 510 CONTINUE 162 500 CONTINUE 163 164 RETURN 165C 166C 167C ...the case P = s-shell and Q > p-shell. 168C 169C 170 6 M = 0 171 N = 0 172 DO 600 IJ = 1,MIJ 173 PIJ = P (IJ) 174 DO 610 KL = 1,MKL 175 TWOQ = QINVHF (KL) 176 M = M + 1 177 PSCALE = PIJ * PQPINV (M) 178 179 DO 620 NG = 1,NGQP 180 N = N + 1 181 B01 (N) = (ONE - PSCALE * RTS (N)) * TWOQ 182 620 CONTINUE 183 184 610 CONTINUE 185 600 CONTINUE 186 187 RETURN 188C 189C 190C ...the case P > p-shell and Q = p-shell. 191C 192C 193 7 M = 0 194 N = 0 195 DO 700 IJ = 1,MIJ 196 PIJ = P (IJ) 197 TWOP = PINVHF (IJ) 198 DO 710 KL = 1,MKL 199 M = M + 1 200 TWOPQ = HALF * PQPINV (M) 201 QSCALE = ONE - PIJ * PQPINV (M) 202 203 DO 720 NG = 1,NGQP 204 N = N + 1 205 ROOT = RTS (N) 206 B00 (N) = ROOT * TWOPQ 207 B10 (N) = (ONE - QSCALE * ROOT) * TWOP 208 720 CONTINUE 209 210 710 CONTINUE 211 700 CONTINUE 212 213 RETURN 214C 215C 216C ...the case P = p-shell and Q > p-shell. 217C 218C 219 8 M = 0 220 N = 0 221 DO 800 IJ = 1,MIJ 222 PIJ = P (IJ) 223 DO 810 KL = 1,MKL 224 TWOQ = QINVHF (KL) 225 M = M + 1 226 TWOPQ = HALF * PQPINV (M) 227 PSCALE = PIJ * PQPINV (M) 228 229 DO 820 NG = 1,NGQP 230 N = N + 1 231 ROOT = RTS (N) 232 B00 (N) = ROOT * TWOPQ 233 B01 (N) = (ONE - PSCALE * ROOT) * TWOQ 234 820 CONTINUE 235 236 810 CONTINUE 237 800 CONTINUE 238 239 RETURN 240C 241C 242C ...the case P > p-shell and Q > p-shell. 243C 244C 245 9 M = 0 246 N = 0 247 DO 900 IJ = 1,MIJ 248 PIJ = P (IJ) 249 TWOP = PINVHF (IJ) 250 DO 910 KL = 1,MKL 251 TWOQ = QINVHF (KL) 252 M = M + 1 253 TWOPQ = HALF * PQPINV (M) 254 PSCALE = PIJ * PQPINV (M) 255 QSCALE = ONE - PSCALE 256 257 DO 920 NG = 1,NGQP 258 N = N + 1 259 ROOT = RTS (N) 260 B00 (N) = ROOT * TWOPQ 261 B01 (N) = (ONE - PSCALE * ROOT) * TWOQ 262 B10 (N) = (ONE - QSCALE * ROOT) * TWOP 263 920 CONTINUE 264 265 910 CONTINUE 266 900 CONTINUE 267C 268C 269C ...ready! 270C 271C 272 RETURN 273 END 274