From 3d3ec3d7aa41664271bdfebe8fb8566f77e1e818 Mon Sep 17 00:00:00 2001 From: Toshi Nagata Date: Sun, 16 Oct 2022 16:27:02 +0900 Subject: [PATCH] F/F7/G/G9 orbitals are implemented --- MolLib/Molecule.c | 291 ++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 270 insertions(+), 21 deletions(-) diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c index 4672586..ac8f541 100755 --- a/MolLib/Molecule.c +++ b/MolLib/Molecule.c @@ -3032,12 +3032,25 @@ sReadNumberArray(void *basep, Int *countp, Int size, Int num, FILE *fp, int *lnp return 3; /* Unexpected EOF */ } +// Normalization constant for one gaussian component +// 1/sqrt(Integral((Y(lm)*(r^n)*exp(-a*r*r))^2, for all r = (x, y, z))) +// where Y(lm) is a spherical harmonic function, r^n is an "additional exponent" +// required in expanded Molden file generated by JANPA, and a is the exponent +// of the gaussian component. +// The function Y(lm) is assumed so that its norm equals sqrt(4*pi/(2l+1)) +// for each m in [-l..l]. +static double +sGaussianNormalizationConstant(int l, double a, int n) +{ + return 1.0/(sqrt(4 * PI / (2 * l + 1.0)) * sqrt(tgamma(l + n + 1.5) / (2.0 * pow(2.0 * a, l + n + 1.5)))); +} + static int sSetupGaussianCoefficients(BasisSet *bset) { ShellInfo *sp; PrimInfo *pp; - int i, j, k; + int i, j, k, n; Double *dp, d; /* Cache the contraction coefficients for efficient calculation */ @@ -3053,44 +3066,134 @@ sSetupGaussianCoefficients(BasisSet *bset) dp = bset->cns; for (i = 0, sp = bset->shells; i < bset->nshells; i++, sp++) { for (j = 0, pp = bset->priminfos + sp->p_idx; j < sp->nprim; j++, pp++) { + n = sp->add_exp; switch (sp->sym) { case kGTOType_S: + // GNC(0,a,n) * r^n * exp(-a*r^2) + d = pp->C * sGaussianNormalizationConstant(0, pp->A, n); + *dp++ = d; + //{ printf("type_S: %g %g\n", d, pp->C * pow(pp->A, 0.75) * 0.71270547); } // (8 alpha^3/pi^3)^0.25 exp(-alpha r^2) - *dp++ = pp->C * pow(pp->A, 0.75) * 0.71270547; + //*dp++ = pp->C * pow(pp->A, 0.75) * 0.71270547; break; case kGTOType_P: - // (128 alpha^5/pi^3)^0.25 [x|y|z]exp(-alpha r^2) - d = pp->C * pow(pp->A, 1.25) * 1.425410941; + // GNC(1,a,n) * [x|y|z] * r^n * exp(-a*r^2) + d = pp->C * sGaussianNormalizationConstant(1, pp->A, n); + //{ printf("type_P: %g %g\n", d, pp->C * pow(pp->A, 1.25) * 1.425410941); } + // (128 alpha^5/pi^3)^0.25 [x|y|z]exp(-alpha r^2) + // d = pp->C * pow(pp->A, 1.25) * 1.425410941; *dp++ = d; *dp++ = d; *dp++ = d; break; case kGTOType_SP: - *dp++ = pp->C * pow(pp->A, 0.75) * 0.71270547; - d = pp->Csp * pow(pp->A, 1.25) * 1.425410941; + // GNC(0,a,n) * r^n * exp(-a*r^2) + *dp++ = d = pp->C * sGaussianNormalizationConstant(0, pp->A, n); + //{ printf("type_SP(s): %g %g\n", d, pp->C * pow(pp->A, 0.75) * 0.71270547); } + // GNC(1,a,n) * [x|y|z] * r^n * exp(-a*r^2) + d = pp->Csp * sGaussianNormalizationConstant(1, pp->A, n); + //{ printf("type_SP(p): %g %g\n", d, pp->Csp * pow(pp->A, 1.25) * 1.425410941); } + //*dp++ = pp->C * pow(pp->A, 0.75) * 0.71270547; + //d = pp->Csp * pow(pp->A, 1.25) * 1.425410941; *dp++ = d; *dp++ = d; *dp++ = d; break; case kGTOType_D: + // GNC(2,a,n) * [xx|yy|zz] * r^n * exp(-a*r^2) + // GNC(2,a,n) * sqrt(3) * [xy|yz|zx] * r^n * exp(-a*r^2) + d = pp->C * sGaussianNormalizationConstant(2, pp->A, n); + //{ printf("type_D[0-2]: %g %g\n", d, pp->C * pow(pp->A, 1.75) * 1.645922781); } + //{ printf("type_D[3-5]: %g %g\n", d * sqrt(3), pp->C * pow(pp->A, 1.75) * 2.850821881); } + dp[0] = dp[1] = dp[2] = d; + dp[3] = dp[4] = dp[5] = d * sqrt(3); // xx|yy|zz: (2048 alpha^7/9pi^3)^0.25 [xx|yy|zz]exp(-alpha r^2) // xy|yz|zx: (2048 alpha^7/pi^3)^0.25 [xy|xz|yz]exp(-alpha r^2) - d = pp->C * pow(pp->A, 1.75); - dp[0] = dp[1] = dp[2] = d * 1.645922781; - dp[3] = dp[4] = dp[5] = d * 2.850821881; + // d = pp->C * pow(pp->A, 1.75); + //dp[0] = dp[1] = dp[2] = d * 1.645922781; + //dp[3] = dp[4] = dp[5] = d * 2.850821881; dp += 6; break; case kGTOType_D5: + // D(0): GNC(2,a,n) * (1/2) * (3zz-rr) * r^n * exp(-a*r^2) + // D(+1): GNC(2,a,n) * sqrt(3) * xz * r^n * exp(-a*r^2) + // D(-1): GNC(2,a,n) * sqrt(3) * yz * r^n * exp(-a*r^2) + // D(+2): GNC(2,a,n) * (sqrt(3)/2) * (xx-yy) * r^n * exp(-a*r^2) + // D(-2): GNC(2,a,n) * sqrt(3) * xy * r^n * exp(-a*r^2) + d = pp->C * sGaussianNormalizationConstant(2, pp->A, n); + //{ printf("type_D5[0]: %g %g\n", d * 0.5, pp->C * pow(pp->A, 1.75) * 0.822961390); } + //{ printf("type_D5[1,2,4]: %g %g\n", d * sqrt(3), pp->C * pow(pp->A, 1.75) * 2.850821881); } + //{ printf("type_D5[3]: %g %g\n", d * sqrt(3) * 0.5, pp->C * pow(pp->A, 1.75) * 1.425410941); } + dp[0] = d * 0.5; + dp[1] = dp[2] = dp[4] = d * sqrt(3); + dp[3] = d * sqrt(3) * 0.5; // 3zz-rr: (128 alpha^7/9pi^3)^0.25 (3zz-rr)exp(-alpha r^2) // xy|yz|zx: (2048 alpha^7/pi^3)^0.25 [xy|xz|yz]exp(-alpha r^2) // xx-yy: (128 alpha^7/pi^3)^0.25 (xx-yy)exp(-alpha r^2) - d = pp->C * pow(pp->A, 1.75); - dp[0] = d * 0.822961390; - dp[1] = dp[2] = dp[4] = d * 2.850821881; - dp[3] = d * 1.425410941; + //d = pp->C * pow(pp->A, 1.75); + //dp[0] = d * 0.822961390; + //dp[1] = dp[2] = dp[4] = d * 2.850821881; + //dp[3] = d * 1.425410941; dp += 5; break; - /* TODO: Support F/F7 and G/G9 type orbitals */ + case kGTOType_F: + // GNC(3,a,n) * [xxx|yyy|zzz] * r^n * exp(-a*r^2) + // GNC(3,a,n) * sqrt(5) * [xyy|xxy|xxz|xzz|yzz|yyz] * r^n * exp(-a*r^2) + // GNC(3,a,n) * sqrt(15) * xyz * r^n * exp(-a*r^2) + d = pp->C * sGaussianNormalizationConstant(3, pp->A, n); + dp[0] = dp[1] = dp[2] = d; + dp[3] = dp[4] = dp[5] = dp[6] = dp[7] = dp[8] = d * sqrt(5); + dp[9] = d * sqrt(15); + dp += 10; + break; + case kGTOType_F7: + // F(0): GNC(3,a,n) * (1/2) * (5zzz-3zrr) * r^n * exp(-a*r^2) + // F(+1): GNC(3,a,n) * sqrt(3/8) * (5xzz-xrr) * r^n * exp(-a*r^2) + // F(-1): GNC(3,a,n) * sqrt(3/8) * (5yzz-yrr) * r^n * exp(-a*r^2) + // F(+2): GNC(3,a,n) * sqrt(15/4) * (xxz-yyz) * r^n * exp(-a*r^2) + // F(-2): GNC(3,a,n) * sqrt(15) * xyz * r^n * exp(-a*r^2) + // F(+3): GNC(3,a,n) * sqrt(5/8) * (xxx-3xyy) * r^n * exp(-a*r^2) + // F(-3): GNC(3,a,n) * sqrt(5/8) * (3xxy-yyy) * r^n * exp(-a*r^2) + d = pp->C * sGaussianNormalizationConstant(3, pp->A, n); + dp[0] = d * 0.5; + dp[1] = dp[2] = d * sqrt(3/8.0); + dp[3] = d * sqrt(15/4.0); + dp[4] = d * sqrt(15); + dp[5] = dp[6] = d * sqrt(5/8.0); + dp += 7; + break; + case kGTOType_G: + // GNC(4,a,n) * [xxxx|yyyy|zzzz] * exp(-a*r^2) + // GNC(4,a,n) * sqrt(7) * [xxxy|xxxz|yyyx|yyyz|zzzx|zzzy] * exp(-a*r^2) + // GNC(4,a,n) * sqrt(35/3) * [xxyy|xxzz|yyzz] * exp(-a*r^2) + // GNC(4,a,n) * sqrt(35) * [xxyz|yyzx|zzxy] * exp(-a*r^2) + d = pp->C * sGaussianNormalizationConstant(4, pp->A, n); + dp[0] = dp[1] = dp[2] = d; + dp[3] = dp[4] = dp[5] = dp[6] = dp[7] = dp[8] = d * sqrt(7); + dp[9] = dp[10] = dp[11] = d * sqrt(35/3.0); + dp[12] = dp[13] = dp[14] = d * sqrt(35); + dp += 15; + break; + case kGTOType_G9: + // G(0): GNC(4,a,n) * (1/8) * (35zzzz-30zzrr+3rrrr) * exp(-a*r^2) + // G(+1): GNC(4,a,n) * sqrt(5/8) * (7xzzz-3xzrr) * exp(-a*r^2) + // G(-1): GNC(4,a,n) * sqrt(5/8) * (7yzzz-3yzrr) * exp(-a*r^2) + // G(+2): GNC(4,a,n) * sqrt(5/16) * (xx-yy)(7zz-rr) * exp(-a*r^2) + // G(-2): GNC(4,a,n) * sqrt(5/4) * (7xyzz-xyrr) * exp(-a*r^2) + // G(+3): GNC(4,a,n) * sqrt(35/8) * (xxxz-3xyyz) * exp(-a*r^2) + // G(-3): GNC(4,a,n) * sqrt(35/8) * (3xxyz-yyyz) * exp(-a*r^2) + // G(+4): GNC(4,a,n) * sqrt(35/64) * (xxxx-6xxyy+yyyy) * exp(-a*r^2) + // G(-4): GNC(4,a,n) * sqrt(35/4) * (xxxy-xyyy) * exp(-a*r^2) + d = pp->C * sGaussianNormalizationConstant(4, pp->A, n); + dp[0] = d * 0.125; + dp[1] = dp[2] = d * sqrt(5/8.0); + dp[3] = d * sqrt(5/16.0); + dp[4] = d * sqrt(5/4.0); + dp[5] = dp[6] = d * sqrt(35/8.0); + dp[7] = d * sqrt(35/64.0); + dp[8] = d * sqrt(35/4.0); + dp += 9; + break; } } } @@ -11360,17 +11463,22 @@ sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Do val = 0.0; mobasep = bset->mo + (index - 1) * bset->ncomps; for (i = 0, sp = bset->shells; i < bset->nshells; i++, sp++) { - pp = bset->priminfos + sp->p_idx; + Double rn; + pp = bset->priminfos + sp->p_idx; cnp = bset->cns + sp->cn_idx; if (sp->a_idx >= mp->natoms) return 0.0; /* This may happen when molecule is edited after setting up MO info */ tmpp = tmp + sp->a_idx * 4; mop = mobasep + sp->m_idx; + if (sp->add_exp == 0) + rn = 1.0; + else + rn = pow(tmpp[3], sp->add_exp * 0.5); switch (sp->sym) { case kGTOType_S: { tval = 0; for (j = 0; j < sp->nprim; j++) { - tval += *cnp++ * exp(-pp->A * tmpp[3]); + tval += *cnp++ * rn * exp(-pp->A * tmpp[3]); pp++; } val += mop[0] * tval; @@ -11380,7 +11488,7 @@ sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Do Double x, y, z; x = y = z = 0; for (j = 0; j < sp->nprim; j++) { - tval = exp(-pp->A * tmpp[3]); + tval = rn * exp(-pp->A * tmpp[3]); x += *cnp++ * tval; y += *cnp++ * tval; z += *cnp++ * tval; @@ -11396,7 +11504,7 @@ sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Do Double t, x, y, z; t = x = y = z = 0; for (j = 0; j < sp->nprim; j++) { - tval = exp(-pp->A * tmpp[3]); + tval = rn * exp(-pp->A * tmpp[3]); t += *cnp++ * tval; x += *cnp++ * tval; y += *cnp++ * tval; @@ -11414,7 +11522,7 @@ sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Do Double xx, yy, zz, xy, xz, yz; xx = yy = zz = xy = xz = yz = 0; for (j = 0; j < sp->nprim; j++) { - tval = exp(-pp->A * tmpp[3]); + tval = rn * exp(-pp->A * tmpp[3]); xx += *cnp++ * tval; yy += *cnp++ * tval; zz += *cnp++ * tval; @@ -11436,7 +11544,7 @@ sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Do Double d0, d1p, d1n, d2p, d2n; d0 = d1p = d1n = d2p = d2n = 0; for (j = 0; j < sp->nprim; j++) { - tval = exp(-pp->A * tmpp[3]); + tval = rn * exp(-pp->A * tmpp[3]); d0 += *cnp++ * tval; d1p += *cnp++ * tval; d1n += *cnp++ * tval; @@ -11452,7 +11560,148 @@ sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Do val += d0 + d1p + d1n + d2p + d2n; break; } - /* TODO: Support F/F7 and G/G9 type orbitals */ + case kGTOType_F: { + Double xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz; + xxx = yyy = zzz = xyy = xxy = xxz = xzz = yzz = yyz = xyz = 0; + for (j = 0; j < sp->nprim; j++) { + tval = rn * exp(-pp->A * tmpp[3]); + xxx += *cnp++ * tval; + yyy += *cnp++ * tval; + zzz += *cnp++ * tval; + xyy += *cnp++ * tval; + xxy += *cnp++ * tval; + xxz += *cnp++ * tval; + xzz += *cnp++ * tval; + yzz += *cnp++ * tval; + yyz += *cnp++ * tval; + xyz += *cnp++ * tval; + pp++; + } + xxx *= mop[0] * tmpp[0] * tmpp[0] * tmpp[0]; + yyy *= mop[1] * tmpp[1] * tmpp[1] * tmpp[1]; + zzz *= mop[2] * tmpp[2] * tmpp[2] * tmpp[2]; + xyy *= mop[3] * tmpp[0] * tmpp[1] * tmpp[1]; + xxy *= mop[4] * tmpp[0] * tmpp[0] * tmpp[1]; + xxz *= mop[5] * tmpp[0] * tmpp[0] * tmpp[2]; + xzz *= mop[6] * tmpp[0] * tmpp[2] * tmpp[2]; + yzz *= mop[7] * tmpp[1] * tmpp[2] * tmpp[2]; + yyz *= mop[8] * tmpp[1] * tmpp[1] * tmpp[2]; + xyz *= mop[9] * tmpp[0] * tmpp[1] * tmpp[2]; + val += xxx + yyy + zzz + xyy + xxy + xxz + xzz + yzz + yyz + xyz; + break; + } + case kGTOType_F7: { + Double f0, f1p, f1n, f2p, f2n, f3p, f3n; + f0 = f1p = f1n = f2p = f2n = f3p = f3n = 0; + for (j = 0; j < sp->nprim; j++) { + tval = rn * exp(-pp->A * tmpp[3]); + f0 += *cnp++ * tval; + f1p += *cnp++ * tval; + f1n += *cnp++ * tval; + f2p += *cnp++ * tval; + f2n += *cnp++ * tval; + f3p += *cnp++ * tval; + f3n += *cnp++ * tval; + pp++; + } + // F(0): GNC(3,a,n) * (1/2) * (5zzz-3zrr) * r^n * exp(-a*r^2) + // F(+1): GNC(3,a,n) * sqrt(3/8) * (5xzz-xrr) * r^n * exp(-a*r^2) + // F(-1): GNC(3,a,n) * sqrt(3/8) * (5yzz-yrr) * r^n * exp(-a*r^2) + // F(+2): GNC(3,a,n) * sqrt(15/4) * (xxz-yyz) * r^n * exp(-a*r^2) + // F(-2): GNC(3,a,n) * sqrt(15) * xyz * r^n * exp(-a*r^2) + // F(+3): GNC(3,a,n) * sqrt(5/8) * (xxx-3xyy) * r^n * exp(-a*r^2) + // F(-3): GNC(3,a,n) * sqrt(5/8) * (3xxy-yyy) * r^n * exp(-a*r^2) + f0 *= mop[0] * tmpp[2] * (5 * tmpp[2] * tmpp[2] - 3 * tmpp[3]); + f1p *= mop[1] * tmpp[0] * (5 * tmpp[2] * tmpp[2] - tmpp[3]); + f1n *= mop[2] * tmpp[1] * (5 * tmpp[2] * tmpp[2] - tmpp[3]); + f2p *= mop[3] * tmpp[2] * (tmpp[0] * tmpp[0] - tmpp[1] * tmpp[1]); + f2n *= mop[4] * tmpp[0] * tmpp[1] * tmpp[2]; + f3p *= mop[5] * tmpp[0] * (tmpp[0] * tmpp[0] - 3 * tmpp[1] * tmpp[1]); + f3n *= mop[6] * tmpp[1] * (3 * tmpp[0] * tmpp[0] - tmpp[2] * tmpp[2]); + val += f0 + f1p + f1n + f2p + f2n + f3p + f3n; + break; + } + case kGTOType_G: { + Double xxxx, yyyy, zzzz, xxxy, xxxz, yyyx, yyyz, zzzx, zzzy, xxyy, xxzz, yyzz, xxyz, yyxz, zzxy; + xxxx = yyyy = zzzz = xxxy = xxxz = yyyx = yyyz = zzzx = zzzy = xxyy = xxzz = yyzz = xxyz = yyxz = zzxy = 0; + for (j = 0; j < sp->nprim; j++) { + tval = rn * exp(-pp->A * tmpp[3]); + xxxx += *cnp++ * tval; + yyyy += *cnp++ * tval; + zzzz += *cnp++ * tval; + xxxy += *cnp++ * tval; + xxxz += *cnp++ * tval; + yyyx += *cnp++ * tval; + yyyz += *cnp++ * tval; + zzzx += *cnp++ * tval; + zzzy += *cnp++ * tval; + xxyy += *cnp++ * tval; + xxzz += *cnp++ * tval; + yyzz += *cnp++ * tval; + xxyz += *cnp++ * tval; + yyxz += *cnp++ * tval; + zzxy += *cnp++ * tval; + pp++; + } + xxxx *= mop[0] * tmpp[0] * tmpp[0] * tmpp[0] * tmpp[0]; + yyyy *= mop[1] * tmpp[1] * tmpp[1] * tmpp[1] * tmpp[1]; + zzzz *= mop[2] * tmpp[2] * tmpp[2] * tmpp[2] * tmpp[2]; + xxxy *= mop[3] * tmpp[0] * tmpp[0] * tmpp[0] * tmpp[1]; + xxxz *= mop[4] * tmpp[0] * tmpp[0] * tmpp[0] * tmpp[2]; + yyyx *= mop[5] * tmpp[1] * tmpp[1] * tmpp[1] * tmpp[0]; + yyyz *= mop[6] * tmpp[1] * tmpp[1] * tmpp[1] * tmpp[2]; + zzzx *= mop[7] * tmpp[2] * tmpp[2] * tmpp[2] * tmpp[0]; + zzzy *= mop[8] * tmpp[2] * tmpp[2] * tmpp[2] * tmpp[1]; + xxyy *= mop[9] * tmpp[0] * tmpp[0] * tmpp[1] * tmpp[1]; + xxzz *= mop[10] * tmpp[0] * tmpp[0] * tmpp[2] * tmpp[2]; + yyzz *= mop[11] * tmpp[1] * tmpp[1] * tmpp[2] * tmpp[2]; + xxyz *= mop[12] * tmpp[0] * tmpp[0] * tmpp[1] * tmpp[2]; + yyxz *= mop[13] * tmpp[1] * tmpp[1] * tmpp[0] * tmpp[2]; + zzxy *= mop[14] * tmpp[2] * tmpp[2] * tmpp[0] * tmpp[1]; + val += xxxx + yyyy + zzzz + xxxy + xxxz + yyyx + yyyz + zzzx + zzzy + xxyy + xxzz + yyzz + xxyz + yyxz + zzxy; + break; + } + case kGTOType_G9: { + Double g0, g1p, g1n, g2p, g2n, g3p, g3n, g4p, g4n; + Double xx = tmpp[0] * tmpp[0]; + Double yy = tmpp[1] * tmpp[1]; + Double zz = tmpp[2] * tmpp[2]; + Double rr = tmpp[3]; + g0 = g1p = g1n = g2p = g2n = g3p = g3n = g4p = g4n = 0; + for (j = 0; j < sp->nprim; j++) { + tval = rn * exp(-pp->A * tmpp[3]); + g0 += *cnp++ * tval; + g1p += *cnp++ * tval; + g1n += *cnp++ * tval; + g2p += *cnp++ * tval; + g2n += *cnp++ * tval; + g3p += *cnp++ * tval; + g3n += *cnp++ * tval; + g4p += *cnp++ * tval; + g4n += *cnp++ * tval; + pp++; + } + // G(0): GNC(4,a,n) * (1/8) * (35zzzz-30zzrr+3rrrr) * r^n * exp(-a*r^2) + // G(+1): GNC(4,a,n) * sqrt(5/8) * (7xzzz-3xzrr) * r^n * exp(-a*r^2) + // G(-1): GNC(4,a,n) * sqrt(5/8) * (7yzzz-3yzrr) * r^n * exp(-a*r^2) + // G(+2): GNC(4,a,n) * sqrt(5/16) * (xx-yy)(7zz-rr) * r^n * exp(-a*r^2) + // G(-2): GNC(4,a,n) * sqrt(5/4) * (7xyzz-xyrr) * r^n * exp(-a*r^2) + // G(+3): GNC(4,a,n) * sqrt(35/8) * (xxxz-3xyyz) * r^n * exp(-a*r^2) + // G(-3): GNC(4,a,n) * sqrt(35/8) * (3xxyz-yyyz) * r^n * exp(-a*r^2) + // G(+4): GNC(4,a,n) * sqrt(35/64) * (xxxx-6xxyy+yyyy) * r^n * exp(-a*r^2) + // G(-4): GNC(4,a,n) * sqrt(35/4) * (xxxy-xyyy) * r^n * exp(-a*r^2) + g0 *= mop[0] * (35 * zz * zz - 30 * zz * rr + 3 * rr * rr); + g1p *= mop[1] * tmpp[0] * tmpp[2] * (7 * zz - 3 * rr); + g1n *= mop[2] * tmpp[1] * tmpp[2] * (7 * zz - 3 * rr); + g2p *= mop[3] * (xx - yy) * (7 * zz - rr); + g2n *= mop[4] * tmpp[0] * tmpp[1] * (7 * zz - rr); + g3p *= mop[5] * tmpp[0] * tmpp[2] * (xx - 3 * yy); + g3n *= mop[6] * tmpp[1] * tmpp[2] * (3 * xx - yy); + g4p *= mop[7] * (xx * xx - 6 * xx * yy + yy * yy); + g4n *= mop[8] * tmpp[0] * tmpp[1] * (xx - yy); + val += g0 + g1p + g1n + g2p + g2n + g3p + g3n + g4p + g4n; + break; + } } } return val; -- 2.11.0