795 SUBROUTINE coulomb3(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
796 lc_max, zetc, rpgfc, lc_min, gccc, rab, rab2, rac, rac2, rbc2, vabc, int_abc, &
797 v, f, maxder, vabc_plus)
798 INTEGER,
INTENT(IN) :: la_max, npgfa
799 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: zeta, rpgfa
800 INTEGER,
INTENT(IN) :: la_min, lb_max, npgfb
801 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: zetb, rpgfb
802 INTEGER,
INTENT(IN) :: lb_min, lc_max
803 REAL(kind=
dp),
INTENT(IN) :: zetc, rpgfc
804 INTEGER,
INTENT(IN) :: lc_min
805 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: gccc
806 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: rab
807 REAL(kind=
dp),
INTENT(IN) :: rab2
808 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: rac
809 REAL(kind=
dp),
INTENT(IN) :: rac2, rbc2
810 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: vabc
811 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(OUT) :: int_abc
812 REAL(kind=
dp),
DIMENSION(:, :, :, :) :: v
813 REAL(kind=
dp),
DIMENSION(0:) :: f
814 INTEGER,
INTENT(IN),
OPTIONAL :: maxder
815 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL :: vabc_plus
817 INTEGER :: ax, ay, az, bx, by, bz, coc, cocx, cocy, &
818 cocz, cx, cy, cz, i, ipgf, j, jpgf, k, &
819 kk, la, la_start, lb, lc, &
820 maxder_local, n, na, nap, nb, nmax
821 REAL(kind=
dp) :: dab, dac, dbc, f0, f1, f2, f3, f4, f5, &
822 f6, f7, fcx, fcy, fcz, fx, fy, fz, t, &
824 REAL(kind=
dp),
DIMENSION(3) :: rap, rbp, rcp, rcw, rpw
829 IF (
PRESENT(maxder))
THEN
830 maxder_local = maxder
833 nmax = la_max + lb_max + lc_max + 1
852 IF (rpgfa(ipgf) + rpgfc < dac)
THEN
853 na = na +
ncoset(la_max - maxder_local)
854 nap = nap +
ncoset(la_max)
864 (rpgfb(jpgf) + rpgfc < dbc) .OR. &
865 (rpgfa(ipgf) + rpgfb(jpgf) < dab))
THEN
872 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf))
874 zetw = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc)
876 f0 = 2.0_dp*sqrt(
pi**5*zetw)*zetp*zetq
881 f0 = f0*exp(-zeta(ipgf)*f1*rab2)
884 rcp(:) = rap(:) - rac(:)
889 t = -f4*(rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3))/zetp
891 CALL fgamma(nmax - 1, t, f)
896 v(1, 1, 1, n) = f0*f(n - 1)
909 v(2, 1, 1, n) = rap(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
910 v(3, 1, 1, n) = rap(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
911 v(4, 1, 1, n) = rap(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
925 v(
coset(0, 0, la), 1, 1, n) = &
926 rap(3)*v(
coset(0, 0, la - 1), 1, 1, n) + &
927 rpw(3)*v(
coset(0, 0, la - 1), 1, 1, n + 1) + &
928 f2*real(la - 1,
dp)*(v(
coset(0, 0, la - 2), 1, 1, n) + &
929 f4*v(
coset(0, 0, la - 2), 1, 1, n + 1))
934 v(
coset(0, 1, az), 1, 1, n) = &
935 rap(2)*v(
coset(0, 0, az), 1, 1, n) + &
936 rpw(2)*v(
coset(0, 0, az), 1, 1, n + 1)
940 v(
coset(0, ay, az), 1, 1, n) = &
941 rap(2)*v(
coset(0, ay - 1, az), 1, 1, n) + &
942 rpw(2)*v(
coset(0, ay - 1, az), 1, 1, n + 1) + &
943 f2*real(ay - 1,
dp)*(v(
coset(0, ay - 2, az), 1, 1, n) + &
944 f4*v(
coset(0, ay - 2, az), 1, 1, n + 1))
951 v(
coset(1, ay, az), 1, 1, n) = &
952 rap(1)*v(
coset(0, ay, az), 1, 1, n) + &
953 rpw(1)*v(
coset(0, ay, az), 1, 1, n + 1)
957 f3 = f2*real(ax - 1,
dp)
960 v(
coset(ax, ay, az), 1, 1, n) = &
961 rap(1)*v(
coset(ax - 1, ay, az), 1, 1, n) + &
962 rpw(1)*v(
coset(ax - 1, ay, az), 1, 1, n + 1) + &
963 f3*(v(
coset(ax - 2, ay, az), 1, 1, n) + &
964 f4*v(
coset(ax - 2, ay, az), 1, 1, n + 1))
978 rbp(:) = rap(:) - rab(:)
982 la_start = max(0, la_min - 1)
984 DO la = la_start, la_max - 1
985 DO n = 1, nmax - la - 1
989 v(
coset(ax, ay, az), 2, 1, n) = &
990 v(
coset(ax + 1, ay, az), 1, 1, n) - &
991 rab(1)*v(
coset(ax, ay, az), 1, 1, n)
992 v(
coset(ax, ay, az), 3, 1, n) = &
993 v(
coset(ax, ay + 1, az), 1, 1, n) - &
994 rab(2)*v(
coset(ax, ay, az), 1, 1, n)
995 v(
coset(ax, ay, az), 4, 1, n) = &
996 v(
coset(ax, ay, az + 1), 1, 1, n) - &
997 rab(3)*v(
coset(ax, ay, az), 1, 1, n)
1010 DO n = 1, nmax - la_max - 1
1012 fx = f2*real(ax,
dp)
1013 DO ay = 0, la_max - ax
1014 fy = f2*real(ay,
dp)
1015 az = la_max - ax - ay
1016 fz = f2*real(az,
dp)
1019 v(
coset(ax, ay, az), 2, 1, n) = &
1020 rbp(1)*v(
coset(ax, ay, az), 1, 1, n) + &
1021 rpw(1)*v(
coset(ax, ay, az), 1, 1, n + 1)
1023 v(
coset(ax, ay, az), 2, 1, n) = &
1024 rbp(1)*v(
coset(ax, ay, az), 1, 1, n) + &
1025 rpw(1)*v(
coset(ax, ay, az), 1, 1, n + 1) + &
1026 fx*(v(
coset(ax - 1, ay, az), 1, 1, n) + &
1027 f4*v(
coset(ax - 1, ay, az), 1, 1, n + 1))
1031 v(
coset(ax, ay, az), 3, 1, n) = &
1032 rbp(2)*v(
coset(ax, ay, az), 1, 1, n) + &
1033 rpw(2)*v(
coset(ax, ay, az), 1, 1, n + 1)
1035 v(
coset(ax, ay, az), 3, 1, n) = &
1036 rbp(2)*v(
coset(ax, ay, az), 1, 1, n) + &
1037 rpw(2)*v(
coset(ax, ay, az), 1, 1, n + 1) + &
1038 fy*(v(
coset(ax, ay - 1, az), 1, 1, n) + &
1039 f4*v(
coset(ax, ay - 1, az), 1, 1, n + 1))
1043 v(
coset(ax, ay, az), 4, 1, n) = &
1044 rbp(3)*v(
coset(ax, ay, az), 1, 1, n) + &
1045 rpw(3)*v(
coset(ax, ay, az), 1, 1, n + 1)
1047 v(
coset(ax, ay, az), 4, 1, n) = &
1048 rbp(3)*v(
coset(ax, ay, az), 1, 1, n) + &
1049 rpw(3)*v(
coset(ax, ay, az), 1, 1, n + 1) + &
1050 fz*(v(
coset(ax, ay, az - 1), 1, 1, n) + &
1051 f4*v(
coset(ax, ay, az - 1), 1, 1, n + 1))
1067 la_start = max(0, la_min - 1)
1069 DO la = la_start, la_max - 1
1070 DO n = 1, nmax - la - lb
1077 v(
coset(ax, ay, az),
coset(0, 0, lb), 1, n) = &
1078 v(
coset(ax, ay, az + 1),
coset(0, 0, lb - 1), 1, n) - &
1079 rab(3)*v(
coset(ax, ay, az),
coset(0, 0, lb - 1), 1, n)
1085 v(
coset(ax, ay, az),
coset(0, by, bz), 1, n) = &
1086 v(
coset(ax, ay + 1, az),
coset(0, by - 1, bz), 1, n) - &
1087 rab(2)*v(
coset(ax, ay, az),
coset(0, by - 1, bz), 1, n)
1095 v(
coset(ax, ay, az),
coset(bx, by, bz), 1, n) = &
1096 v(
coset(ax + 1, ay, az),
coset(bx - 1, by, bz), 1, n) - &
1097 rab(1)*v(
coset(ax, ay, az),
coset(bx - 1, by, bz), 1, n)
1115 DO n = 1, nmax - la_max - lb
1117 fx = f2*real(ax,
dp)
1118 DO ay = 0, la_max - ax
1119 fy = f2*real(ay,
dp)
1120 az = la_max - ax - ay
1121 fz = f2*real(az,
dp)
1125 f3 = f2*real(lb - 1,
dp)
1128 v(
coset(ax, ay, az),
coset(0, 0, lb), 1, n) = &
1129 rbp(3)*v(
coset(ax, ay, az),
coset(0, 0, lb - 1), 1, n) + &
1130 rpw(3)*v(
coset(ax, ay, az),
coset(0, 0, lb - 1), 1, n + 1) + &
1131 f3*(v(
coset(ax, ay, az),
coset(0, 0, lb - 2), 1, n) + &
1132 f4*v(
coset(ax, ay, az),
coset(0, 0, lb - 2), 1, n + 1))
1134 v(
coset(ax, ay, az),
coset(0, 0, lb), 1, n) = &
1135 rbp(3)*v(
coset(ax, ay, az),
coset(0, 0, lb - 1), 1, n) + &
1136 rpw(3)*v(
coset(ax, ay, az),
coset(0, 0, lb - 1), 1, n + 1) + &
1137 fz*(v(
coset(ax, ay, az - 1),
coset(0, 0, lb - 1), 1, n) + &
1138 f4*v(
coset(ax, ay, az - 1),
coset(0, 0, lb - 1), 1, n + 1)) + &
1139 f3*(v(
coset(ax, ay, az),
coset(0, 0, lb - 2), 1, n) + &
1140 f4*v(
coset(ax, ay, az),
coset(0, 0, lb - 2), 1, n + 1))
1147 v(
coset(ax, ay, az),
coset(0, 1, bz), 1, n) = &
1148 rbp(2)*v(
coset(ax, ay, az),
coset(0, 0, bz), 1, n) + &
1149 rpw(2)*v(
coset(ax, ay, az),
coset(0, 0, bz), 1, n + 1)
1152 f3 = f2*real(by - 1,
dp)
1153 v(
coset(ax, ay, az),
coset(0, by, bz), 1, n) = &
1154 rbp(2)*v(
coset(ax, ay, az),
coset(0, by - 1, bz), 1, n) + &
1155 rpw(2)*v(
coset(ax, ay, az),
coset(0, by - 1, bz), 1, n + 1) + &
1156 f3*(v(
coset(ax, ay, az),
coset(0, by - 2, bz), 1, n) + &
1157 f4*v(
coset(ax, ay, az),
coset(0, by - 2, bz), 1, n + 1))
1161 v(
coset(ax, ay, az),
coset(0, 1, bz), 1, n) = &
1162 rbp(2)*v(
coset(ax, ay, az),
coset(0, 0, bz), 1, n) + &
1163 rpw(2)*v(
coset(ax, ay, az),
coset(0, 0, bz), 1, n + 1) + &
1164 fy*(v(
coset(ax, ay - 1, az),
coset(0, 0, bz), 1, n) + &
1165 f4*v(
coset(ax, ay - 1, az),
coset(0, 0, bz), 1, n + 1))
1168 f3 = f2*real(by - 1,
dp)
1169 v(
coset(ax, ay, az),
coset(0, by, bz), 1, n) = &
1170 rbp(2)*v(
coset(ax, ay, az),
coset(0, by - 1, bz), 1, n) + &
1171 rpw(2)*v(
coset(ax, ay, az),
coset(0, by - 1, bz), 1, n + 1) + &
1172 fy*(v(
coset(ax, ay - 1, az),
coset(0, by - 1, bz), 1, n) + &
1173 f4*v(
coset(ax, ay - 1, az), &
1174 coset(0, by - 1, bz), 1, n + 1)) + &
1175 f3*(v(
coset(ax, ay, az),
coset(0, by - 2, bz), 1, n) + &
1176 f4*v(
coset(ax, ay, az),
coset(0, by - 2, bz), 1, n + 1))
1185 v(
coset(ax, ay, az),
coset(1, by, bz), 1, n) = &
1186 rbp(1)*v(
coset(ax, ay, az),
coset(0, by, bz), 1, n) + &
1187 rpw(1)*v(
coset(ax, ay, az),
coset(0, by, bz), 1, n + 1)
1190 f3 = f2*real(bx - 1,
dp)
1193 v(
coset(ax, ay, az),
coset(bx, by, bz), 1, n) = &
1194 rbp(1)*v(
coset(ax, ay, az),
coset(bx - 1, by, bz), 1, n) + &
1195 rpw(1)*v(
coset(ax, ay, az), &
1196 coset(bx - 1, by, bz), 1, n + 1) + &
1197 f3*(v(
coset(ax, ay, az),
coset(bx - 2, by, bz), 1, n) + &
1198 f4*v(
coset(ax, ay, az),
coset(bx - 2, by, bz), 1, n + 1))
1204 v(
coset(ax, ay, az),
coset(1, by, bz), 1, n) = &
1205 rbp(1)*v(
coset(ax, ay, az),
coset(0, by, bz), 1, n) + &
1206 rpw(1)*v(
coset(ax, ay, az),
coset(0, by, bz), 1, n + 1) + &
1207 fx*(v(
coset(ax - 1, ay, az),
coset(0, by, bz), 1, n) + &
1208 f4*v(
coset(ax - 1, ay, az),
coset(0, by, bz), 1, n + 1))
1211 f3 = f2*real(bx - 1,
dp)
1214 v(
coset(ax, ay, az),
coset(bx, by, bz), 1, n) = &
1215 rbp(1)*v(
coset(ax, ay, az),
coset(bx - 1, by, bz), 1, n) + &
1216 rpw(1)*v(
coset(ax, ay, az), &
1217 coset(bx - 1, by, bz), 1, n + 1) + &
1218 fx*(v(
coset(ax - 1, ay, az), &
1219 coset(bx - 1, by, bz), 1, n) + &
1220 f4*v(
coset(ax - 1, ay, az), &
1221 coset(bx - 1, by, bz), 1, n + 1)) + &
1222 f3*(v(
coset(ax, ay, az),
coset(bx - 2, by, bz), 1, n) + &
1223 f4*v(
coset(ax, ay, az),
coset(bx - 2, by, bz), 1, n + 1))
1238 IF (lb_max > 0)
THEN
1242 rbp(:) = rap(:) - rab(:)
1248 v(1, 2, 1, n) = rbp(1)*v(1, 1, 1, n) + rpw(1)*v(1, 1, 1, n + 1)
1249 v(1, 3, 1, n) = rbp(2)*v(1, 1, 1, n) + rpw(2)*v(1, 1, 1, n + 1)
1250 v(1, 4, 1, n) = rbp(3)*v(1, 1, 1, n) + rpw(3)*v(1, 1, 1, n + 1)
1264 v(1,
coset(0, 0, lb), 1, n) = &
1265 rbp(3)*v(1,
coset(0, 0, lb - 1), 1, n) + &
1266 rpw(3)*v(1,
coset(0, 0, lb - 1), 1, n + 1) + &
1267 f2*real(lb - 1,
dp)*(v(1,
coset(0, 0, lb - 2), 1, n) + &
1268 f4*v(1,
coset(0, 0, lb - 2), 1, n + 1))
1273 v(1,
coset(0, 1, bz), 1, n) = &
1274 rbp(2)*v(1,
coset(0, 0, bz), 1, n) + &
1275 rpw(2)*v(1,
coset(0, 0, bz), 1, n + 1)
1279 v(1,
coset(0, by, bz), 1, n) = &
1280 rbp(2)*v(1,
coset(0, by - 1, bz), 1, n) + &
1281 rpw(2)*v(1,
coset(0, by - 1, bz), 1, n + 1) + &
1282 f2*real(by - 1,
dp)*(v(1,
coset(0, by - 2, bz), 1, n) + &
1283 f4*v(1,
coset(0, by - 2, bz), 1, n + 1))
1290 v(1,
coset(1, by, bz), 1, n) = &
1291 rbp(1)*v(1,
coset(0, by, bz), 1, n) + &
1292 rpw(1)*v(1,
coset(0, by, bz), 1, n + 1)
1296 f3 = f2*real(bx - 1,
dp)
1299 v(1,
coset(bx, by, bz), 1, n) = &
1300 rbp(1)*v(1,
coset(bx - 1, by, bz), 1, n) + &
1301 rpw(1)*v(1,
coset(bx - 1, by, bz), 1, n + 1) + &
1302 f3*(v(1,
coset(bx - 2, by, bz), 1, n) + &
1303 f4*v(1,
coset(bx - 2, by, bz), 1, n + 1))
1317 IF (lc_max > 0)
THEN
1325 rcw(:) = rcp(:) + rpw(:)
1330 v(1, 1, 2, n) = rcw(1)*v(1, 1, 1, n + 1)
1331 v(1, 1, 3, n) = rcw(2)*v(1, 1, 1, n + 1)
1332 v(1, 1, 4, n) = rcw(3)*v(1, 1, 1, n + 1)
1345 v(1, 1,
coset(0, 0, lc), n) = &
1346 rcw(3)*v(1, 1,
coset(0, 0, lc - 1), n + 1) + &
1347 f7*real(lc - 1,
dp)*(v(1, 1,
coset(0, 0, lc - 2), n) + &
1348 f5*v(1, 1,
coset(0, 0, lc - 2), n + 1))
1353 v(1, 1,
coset(0, 1, cz), n) = rcw(2)*v(1, 1,
coset(0, 0, cz), n + 1)
1357 v(1, 1,
coset(0, cy, cz), n) = &
1358 rcw(2)*v(1, 1,
coset(0, cy - 1, cz), n + 1) + &
1359 f7*real(cy - 1,
dp)*(v(1, 1,
coset(0, cy - 2, cz), n) + &
1360 f5*v(1, 1,
coset(0, cy - 2, cz), n + 1))
1367 v(1, 1,
coset(1, cy, cz), n) = rcw(1)*v(1, 1,
coset(0, cy, cz), n + 1)
1373 v(1, 1,
coset(cx, cy, cz), n) = &
1374 rcw(1)*v(1, 1,
coset(cx - 1, cy, cz), n + 1) + &
1375 f7*real(cx - 1,
dp)*(v(1, 1,
coset(cx - 2, cy, cz), n) + &
1376 f5*v(1, 1,
coset(cx - 2, cy, cz), n + 1))
1392 coc =
coset(cx, cy, cz)
1393 cocx =
coset(max(0, cx - 1), cy, cz)
1394 cocy =
coset(cx, max(0, cy - 1), cz)
1395 cocz =
coset(cx, cy, max(0, cz - 1))
1397 fcx = f6*real(cx,
dp)
1398 fcy = f6*real(cy,
dp)
1399 fcz = f6*real(cz,
dp)
1403 IF (la_max > 0)
THEN
1411 DO n = 1, nmax - 1 - lc
1412 v(2, 1, coc, n) = rap(1)*v(1, 1, coc, n) + &
1413 rpw(1)*v(1, 1, coc, n + 1) + &
1414 fcx*v(1, 1, cocx, n + 1)
1415 v(3, 1, coc, n) = rap(2)*v(1, 1, coc, n) + &
1416 rpw(2)*v(1, 1, coc, n + 1) + &
1417 fcy*v(1, 1, cocy, n + 1)
1418 v(4, 1, coc, n) = rap(3)*v(1, 1, coc, n) + &
1419 rpw(3)*v(1, 1, coc, n + 1) + &
1420 fcz*v(1, 1, cocz, n + 1)
1431 DO n = 1, nmax - la - lc
1435 v(
coset(0, 0, la), 1, coc, n) = &
1436 rap(3)*v(
coset(0, 0, la - 1), 1, coc, n) + &
1437 rpw(3)*v(
coset(0, 0, la - 1), 1, coc, n + 1) + &
1438 f2*real(la - 1,
dp)*(v(
coset(0, 0, la - 2), 1, coc, n) + &
1439 f4*v(
coset(0, 0, la - 2), 1, coc, n + 1)) + &
1440 fcz*v(
coset(0, 0, la - 1), 1, cocz, n + 1)
1445 v(
coset(0, 1, az), 1, coc, n) = &
1446 rap(2)*v(
coset(0, 0, az), 1, coc, n) + &
1447 rpw(2)*v(
coset(0, 0, az), 1, coc, n + 1) + &
1448 fcy*v(
coset(0, 0, az), 1, cocy, n + 1)
1451 f3 = f2*real(ay - 1,
dp)
1453 v(
coset(0, ay, az), 1, coc, n) = &
1454 rap(2)*v(
coset(0, ay - 1, az), 1, coc, n) + &
1455 rpw(2)*v(
coset(0, ay - 1, az), 1, coc, n + 1) + &
1456 f3*(v(
coset(0, ay - 2, az), 1, coc, n) + &
1457 f4*v(
coset(0, ay - 2, az), 1, coc, n + 1)) + &
1458 fcy*v(
coset(0, ay - 1, az), 1, cocy, n + 1)
1465 v(
coset(1, ay, az), 1, coc, n) = &
1466 rap(1)*v(
coset(0, ay, az), 1, coc, n) + &
1467 rpw(1)*v(
coset(0, ay, az), 1, coc, n + 1) + &
1468 fcx*v(
coset(0, ay, az), 1, cocx, n + 1)
1472 f3 = f2*real(ax - 1,
dp)
1475 v(
coset(ax, ay, az), 1, coc, n) = &
1476 rap(1)*v(
coset(ax - 1, ay, az), 1, coc, n) + &
1477 rpw(1)*v(
coset(ax - 1, ay, az), 1, coc, n + 1) + &
1478 f3*(v(
coset(ax - 2, ay, az), 1, coc, n) + &
1479 f4*v(
coset(ax - 2, ay, az), 1, coc, n + 1)) + &
1480 fcx*v(
coset(ax - 1, ay, az), 1, cocx, n + 1)
1490 IF (lb_max > 0)
THEN
1496 la_start = max(0, la_min - 1)
1498 DO la = la_start, la_max - 1
1499 DO n = 1, nmax - la - 1 - lc
1503 v(
coset(ax, ay, az), 2, coc, n) = &
1504 v(
coset(ax + 1, ay, az), 1, coc, n) - &
1505 rab(1)*v(
coset(ax, ay, az), 1, coc, n)
1506 v(
coset(ax, ay, az), 3, coc, n) = &
1507 v(
coset(ax, ay + 1, az), 1, coc, n) - &
1508 rab(2)*v(
coset(ax, ay, az), 1, coc, n)
1509 v(
coset(ax, ay, az), 4, coc, n) = &
1510 v(
coset(ax, ay, az + 1), 1, coc, n) - &
1511 rab(3)*v(
coset(ax, ay, az), 1, coc, n)
1525 DO n = 1, nmax - la_max - 1 - lc
1527 fx = f2*real(ax,
dp)
1528 DO ay = 0, la_max - ax
1529 fy = f2*real(ay,
dp)
1530 az = la_max - ax - ay
1531 fz = f2*real(az,
dp)
1534 v(
coset(ax, ay, az), 2, coc, n) = &
1535 rbp(1)*v(
coset(ax, ay, az), 1, coc, n) + &
1536 rpw(1)*v(
coset(ax, ay, az), 1, coc, n + 1) + &
1537 fcx*v(
coset(ax, ay, az), 1, cocx, n + 1)
1539 v(
coset(ax, ay, az), 2, coc, n) = &
1540 rbp(1)*v(
coset(ax, ay, az), 1, coc, n) + &
1541 rpw(1)*v(
coset(ax, ay, az), 1, coc, n + 1) + &
1542 fx*(v(
coset(ax - 1, ay, az), 1, coc, n) + &
1543 f4*v(
coset(ax - 1, ay, az), 1, coc, n + 1)) + &
1544 fcx*v(
coset(ax, ay, az), 1, cocx, n + 1)
1548 v(
coset(ax, ay, az), 3, coc, n) = &
1549 rbp(2)*v(
coset(ax, ay, az), 1, coc, n) + &
1550 rpw(2)*v(
coset(ax, ay, az), 1, coc, n + 1) + &
1551 fcy*v(
coset(ax, ay, az), 1, cocy, n + 1)
1553 v(
coset(ax, ay, az), 3, coc, n) = &
1554 rbp(2)*v(
coset(ax, ay, az), 1, coc, n) + &
1555 rpw(2)*v(
coset(ax, ay, az), 1, coc, n + 1) + &
1556 fy*(v(
coset(ax, ay - 1, az), 1, coc, n) + &
1557 f4*v(
coset(ax, ay - 1, az), 1, coc, n + 1)) + &
1558 fcy*v(
coset(ax, ay, az), 1, cocy, n + 1)
1562 v(
coset(ax, ay, az), 4, coc, n) = &
1563 rbp(3)*v(
coset(ax, ay, az), 1, coc, n) + &
1564 rpw(3)*v(
coset(ax, ay, az), 1, coc, n + 1) + &
1565 fcz*v(
coset(ax, ay, az), 1, cocz, n + 1)
1567 v(
coset(ax, ay, az), 4, coc, n) = &
1568 rbp(3)*v(
coset(ax, ay, az), 1, coc, n) + &
1569 rpw(3)*v(
coset(ax, ay, az), 1, coc, n + 1) + &
1570 fz*(v(
coset(ax, ay, az - 1), 1, coc, n) + &
1571 f4*v(
coset(ax, ay, az - 1), 1, coc, n + 1)) + &
1572 fcz*v(
coset(ax, ay, az), 1, cocz, n + 1)
1588 la_start = max(0, la_min - 1)
1590 DO la = la_start, la_max - 1
1591 DO n = 1, nmax - la - lb - lc
1598 v(
coset(ax, ay, az),
coset(0, 0, lb), coc, n) = &
1599 v(
coset(ax, ay, az + 1), &
1600 coset(0, 0, lb - 1), coc, n) - &
1601 rab(3)*v(
coset(ax, ay, az), &
1602 coset(0, 0, lb - 1), coc, n)
1608 v(
coset(ax, ay, az),
coset(0, by, bz), coc, n) = &
1609 v(
coset(ax, ay + 1, az), &
1610 coset(0, by - 1, bz), coc, n) - &
1611 rab(2)*v(
coset(ax, ay, az), &
1612 coset(0, by - 1, bz), coc, n)
1620 v(
coset(ax, ay, az),
coset(bx, by, bz), coc, n) = &
1621 v(
coset(ax + 1, ay, az), &
1622 coset(bx - 1, by, bz), coc, n) - &
1623 rab(1)*v(
coset(ax, ay, az), &
1624 coset(bx - 1, by, bz), coc, n)
1643 DO n = 1, nmax - la_max - lb - lc
1645 fx = f2*real(ax,
dp)
1646 DO ay = 0, la_max - ax
1647 fy = f2*real(ay,
dp)
1648 az = la_max - ax - ay
1649 fz = f2*real(az,
dp)
1653 f3 = f2*real(lb - 1,
dp)
1656 v(
coset(ax, ay, az),
coset(0, 0, lb), coc, n) = &
1657 rbp(3)*v(
coset(ax, ay, az), &
1658 coset(0, 0, lb - 1), coc, n) + &
1659 rpw(3)*v(
coset(ax, ay, az), &
1660 coset(0, 0, lb - 1), coc, n + 1) + &
1661 f3*(v(
coset(ax, ay, az), &
1662 coset(0, 0, lb - 2), coc, n) + &
1663 f4*v(
coset(ax, ay, az), &
1664 coset(0, 0, lb - 2), coc, n + 1)) + &
1665 fcz*v(
coset(ax, ay, az), &
1666 coset(0, 0, lb - 1), cocz, n + 1)
1668 v(
coset(ax, ay, az),
coset(0, 0, lb), coc, n) = &
1669 rbp(3)*v(
coset(ax, ay, az), &
1670 coset(0, 0, lb - 1), coc, n) + &
1671 rpw(3)*v(
coset(ax, ay, az), &
1672 coset(0, 0, lb - 1), coc, n + 1) + &
1673 fz*(v(
coset(ax, ay, az - 1), &
1674 coset(0, 0, lb - 1), coc, n) + &
1675 f4*v(
coset(ax, ay, az - 1), &
1676 coset(0, 0, lb - 1), coc, n + 1)) + &
1677 f3*(v(
coset(ax, ay, az), &
1678 coset(0, 0, lb - 2), coc, n) + &
1679 f4*v(
coset(ax, ay, az), &
1680 coset(0, 0, lb - 2), coc, n + 1)) + &
1681 fcz*v(
coset(ax, ay, az), &
1682 coset(0, 0, lb - 1), cocz, n + 1)
1689 v(
coset(ax, ay, az),
coset(0, 1, bz), coc, n) = &
1690 rbp(2)*v(
coset(ax, ay, az), &
1691 coset(0, 0, bz), coc, n) + &
1692 rpw(2)*v(
coset(ax, ay, az), &
1693 coset(0, 0, bz), coc, n + 1) + &
1694 fcy*v(
coset(ax, ay, az), &
1695 coset(0, 0, bz), cocy, n + 1)
1698 f3 = f2*real(by - 1,
dp)
1699 v(
coset(ax, ay, az),
coset(0, by, bz), coc, n) = &
1700 rbp(2)*v(
coset(ax, ay, az), &
1701 coset(0, by - 1, bz), coc, n) + &
1702 rpw(2)*v(
coset(ax, ay, az), &
1703 coset(0, by - 1, bz), coc, n + 1) + &
1704 f3*(v(
coset(ax, ay, az), &
1705 coset(0, by - 2, bz), coc, n) + &
1706 f4*v(
coset(ax, ay, az), &
1707 coset(0, by - 2, bz), coc, n + 1)) + &
1708 fcy*v(
coset(ax, ay, az), &
1709 coset(0, by - 1, bz), cocy, n + 1)
1713 v(
coset(ax, ay, az),
coset(0, 1, bz), coc, n) = &
1714 rbp(2)*v(
coset(ax, ay, az), &
1715 coset(0, 0, bz), coc, n) + &
1716 rpw(2)*v(
coset(ax, ay, az), &
1717 coset(0, 0, bz), coc, n + 1) + &
1718 fy*(v(
coset(ax, ay - 1, az), &
1719 coset(0, 0, bz), coc, n) + &
1720 f4*v(
coset(ax, ay - 1, az), &
1721 coset(0, 0, bz), coc, n + 1)) + &
1722 fcy*v(
coset(ax, ay, az), &
1723 coset(0, 0, bz), cocy, n + 1)
1726 f3 = f2*real(by - 1,
dp)
1727 v(
coset(ax, ay, az),
coset(0, by, bz), coc, n) = &
1728 rbp(2)*v(
coset(ax, ay, az), &
1729 coset(0, by - 1, bz), coc, n) + &
1730 rpw(2)*v(
coset(ax, ay, az), &
1731 coset(0, by - 1, bz), coc, n + 1) + &
1732 fy*(v(
coset(ax, ay - 1, az), &
1733 coset(0, by - 1, bz), coc, n) + &
1734 f4*v(
coset(ax, ay - 1, az), &
1735 coset(0, by - 1, bz), coc, n + 1)) + &
1736 f3*(v(
coset(ax, ay, az), &
1737 coset(0, by - 2, bz), coc, n) + &
1738 f4*v(
coset(ax, ay, az), &
1739 coset(0, by - 2, bz), coc, n + 1)) + &
1740 fcy*v(
coset(ax, ay, az), &
1741 coset(0, by - 1, bz), cocy, n + 1)
1750 v(
coset(ax, ay, az),
coset(1, by, bz), coc, n) = &
1751 rbp(1)*v(
coset(ax, ay, az), &
1752 coset(0, by, bz), coc, n) + &
1753 rpw(1)*v(
coset(ax, ay, az), &
1754 coset(0, by, bz), coc, n + 1) + &
1755 fcx*v(
coset(ax, ay, az), &
1756 coset(0, by, bz), cocx, n + 1)
1759 f3 = f2*real(bx - 1,
dp)
1762 v(
coset(ax, ay, az),
coset(bx, by, bz), coc, n) = &
1763 rbp(1)*v(
coset(ax, ay, az), &
1764 coset(bx - 1, by, bz), coc, n) + &
1765 rpw(1)*v(
coset(ax, ay, az), &
1766 coset(bx - 1, by, bz), coc, n + 1) + &
1767 f3*(v(
coset(ax, ay, az), &
1768 coset(bx - 2, by, bz), coc, n) + &
1769 f4*v(
coset(ax, ay, az), &
1770 coset(bx - 2, by, bz), coc, n + 1)) + &
1771 fcx*v(
coset(ax, ay, az), &
1772 coset(bx - 1, by, bz), cocx, n + 1)
1778 v(
coset(ax, ay, az),
coset(1, by, bz), coc, n) = &
1779 rbp(1)*v(
coset(ax, ay, az), &
1780 coset(0, by, bz), coc, n) + &
1781 rpw(1)*v(
coset(ax, ay, az), &
1782 coset(0, by, bz), coc, n + 1) + &
1783 fx*(v(
coset(ax - 1, ay, az), &
1784 coset(0, by, bz), coc, n) + &
1785 f4*v(
coset(ax - 1, ay, az), &
1786 coset(0, by, bz), coc, n + 1)) + &
1787 fcx*v(
coset(ax, ay, az), &
1788 coset(0, by, bz), cocx, n + 1)
1791 f3 = f2*real(bx - 1,
dp)
1794 v(
coset(ax, ay, az),
coset(bx, by, bz), coc, n) = &
1795 rbp(1)*v(
coset(ax, ay, az), &
1796 coset(bx - 1, by, bz), coc, n) + &
1797 rpw(1)*v(
coset(ax, ay, az), &
1798 coset(bx - 1, by, bz), coc, n + 1) + &
1799 fx*(v(
coset(ax - 1, ay, az), &
1800 coset(bx - 1, by, bz), coc, n) + &
1801 f4*v(
coset(ax - 1, ay, az), &
1802 coset(bx - 1, by, bz), coc, n + 1)) + &
1803 f3*(v(
coset(ax, ay, az), &
1804 coset(bx - 2, by, bz), coc, n) + &
1805 f4*v(
coset(ax, ay, az), &
1806 coset(bx - 2, by, bz), coc, n + 1)) + &
1807 fcx*v(
coset(ax, ay, az), &
1808 coset(bx - 1, by, bz), cocx, n + 1)
1822 IF (lb_max > 0)
THEN
1830 DO n = 1, nmax - 1 - lc
1831 v(1, 2, coc, n) = rbp(1)*v(1, 1, coc, n) + &
1832 rpw(1)*v(1, 1, coc, n + 1) + &
1833 fcx*v(1, 1, cocx, n + 1)
1834 v(1, 3, coc, n) = rbp(2)*v(1, 1, coc, n) + &
1835 rpw(2)*v(1, 1, coc, n + 1) + &
1836 fcy*v(1, 1, cocy, n + 1)
1837 v(1, 4, coc, n) = rbp(3)*v(1, 1, coc, n) + &
1838 rpw(3)*v(1, 1, coc, n + 1) + &
1839 fcz*v(1, 1, cocz, n + 1)
1850 DO n = 1, nmax - lb - lc
1854 v(1,
coset(0, 0, lb), coc, n) = &
1855 rbp(3)*v(1,
coset(0, 0, lb - 1), coc, n) + &
1856 rpw(3)*v(1,
coset(0, 0, lb - 1), coc, n + 1) + &
1857 f2*real(lb - 1,
dp)*(v(1,
coset(0, 0, lb - 2), coc, n) + &
1858 f4*v(1,
coset(0, 0, lb - 2), coc, n + 1)) + &
1859 fcz*v(1,
coset(0, 0, lb - 1), cocz, n + 1)
1864 v(1,
coset(0, 1, bz), coc, n) = &
1865 rbp(2)*v(1,
coset(0, 0, bz), coc, n) + &
1866 rpw(2)*v(1,
coset(0, 0, bz), coc, n + 1) + &
1867 fcy*v(1,
coset(0, 0, bz), cocy, n + 1)
1870 f3 = f2*real(by - 1,
dp)
1872 v(1,
coset(0, by, bz), coc, n) = &
1873 rbp(2)*v(1,
coset(0, by - 1, bz), coc, n) + &
1874 rpw(2)*v(1,
coset(0, by - 1, bz), coc, n + 1) + &
1875 f3*(v(1,
coset(0, by - 2, bz), coc, n) + &
1876 f4*v(1,
coset(0, by - 2, bz), coc, n + 1)) + &
1877 fcy*v(1,
coset(0, by - 1, bz), cocy, n + 1)
1884 v(1,
coset(1, by, bz), coc, n) = &
1885 rbp(1)*v(1,
coset(0, by, bz), coc, n) + &
1886 rpw(1)*v(1,
coset(0, by, bz), coc, n + 1) + &
1887 fcx*v(1,
coset(0, by, bz), cocx, n + 1)
1891 f3 = f2*real(bx - 1,
dp)
1894 v(1,
coset(bx, by, bz), coc, n) = &
1895 rbp(1)*v(1,
coset(bx - 1, by, bz), coc, n) + &
1896 rpw(1)*v(1,
coset(bx - 1, by, bz), coc, n + 1) + &
1897 f3*(v(1,
coset(bx - 2, by, bz), coc, n) + &
1898 f4*v(1,
coset(bx - 2, by, bz), coc, n + 1)) + &
1899 fcx*v(1,
coset(bx - 1, by, bz), cocx, n + 1)
1922 kk = k -
ncoset(lc_min - 1)
1924 DO i =
ncoset(la_min - 1) + 1,
ncoset(la_max - maxder_local)
1925 vabc(na + i, nb + j) = vabc(na + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1926 int_abc(na + i, nb + j, kk) = v(i, j, k, 1)
1931 IF (
PRESENT(maxder))
THEN
1933 kk = k -
ncoset(lc_min - 1)
1936 vabc_plus(nap + i, nb + j) = vabc_plus(nap + i, nb + j) + gccc(kk)*v(i, j, k, 1)
1946 na = na +
ncoset(la_max - maxder_local)
1947 nap = nap +
ncoset(la_max)