OSDN > ソフトウェアを探す > 科学/工学 > 分子科学 > Molby - 分子モデリングソフトウェア > SCM > git > Molby > コミット

Molby - 分子モデリングソフトウェア
Fork

R/O
HTTP
SSH
HTTPS

Molby: コミット

Molecular Modeling Software

コミットメタ情報

リビジョン	3d3ec3d7aa41664271bdfebe8fb8566f77e1e818 (tree)
日時	2022-10-16 16:27:02
作者	Toshi Nagata <alchemist.2005@nift...>
コミッター	Toshi Nagata

ログメッセージ

F/F7/G/G9 orbitals are implemented

変更サマリ

modified: MolLib/Molecule.c (diff)

差分

--- 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
3032	3032	return 3; /* Unexpected EOF */
3033	3033	}
3034	3034
	3035	+// Normalization constant for one gaussian component
	3036	+// 1/sqrt(Integral((Y(lm)(r^n)exp(-arr))^2, for all r = (x, y, z)))
	3037	+// where Y(lm) is a spherical harmonic function, r^n is an "additional exponent"
	3038	+// required in expanded Molden file generated by JANPA, and a is the exponent
	3039	+// of the gaussian component.
	3040	+// The function Y(lm) is assumed so that its norm equals sqrt(4*pi/(2l+1))
	3041	+// for each m in [-l..l].
	3042	+static double
	3043	+sGaussianNormalizationConstant(int l, double a, int n)
	3044	+{
	3045	+ 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))));
	3046	+}
	3047	+
3035	3048	static int
3036	3049	sSetupGaussianCoefficients(BasisSet *bset)
3037	3050	{
3038	3051	ShellInfo *sp;
3039	3052	PrimInfo *pp;
3040		- int i, j, k;
	3053	+ int i, j, k, n;
3041	3054	Double *dp, d;
3042	3055
3043	3056	/* Cache the contraction coefficients for efficient calculation */

		@@ -3053,44 +3066,134 @@ sSetupGaussianCoefficients(BasisSet *bset)
3053	3066	dp = bset->cns;
3054	3067	for (i = 0, sp = bset->shells; i < bset->nshells; i++, sp++) {
3055	3068	for (j = 0, pp = bset->priminfos + sp->p_idx; j < sp->nprim; j++, pp++) {
	3069	+ n = sp->add_exp;
3056	3070	switch (sp->sym) {
3057	3071	case kGTOType_S:
	3072	+ // GNC(0,a,n) * r^n * exp(-a*r^2)
	3073	+ d = pp->C * sGaussianNormalizationConstant(0, pp->A, n);
	3074	+ *dp++ = d;
	3075	+ //{ printf("type_S: %g %g\n", d, pp->C * pow(pp->A, 0.75) * 0.71270547); }
3058	3076	// (8 alpha^3/pi^3)^0.25 exp(-alpha r^2)
3059		- dp++ = pp->C pow(pp->A, 0.75) * 0.71270547;
	3077	+ //dp++ = pp->C pow(pp->A, 0.75) * 0.71270547;
3060	3078	break;
3061	3079	case kGTOType_P:
3062		- // (128 alpha^5/pi^3)^0.25 [x\|y\|z]exp(-alpha r^2)
3063		- d = pp->C * pow(pp->A, 1.25) * 1.425410941;
	3080	+ // GNC(1,a,n) * [x\|y\|z] * r^n * exp(-a*r^2)
	3081	+ d = pp->C * sGaussianNormalizationConstant(1, pp->A, n);
	3082	+ //{ printf("type_P: %g %g\n", d, pp->C * pow(pp->A, 1.25) * 1.425410941); }
	3083	+ // (128 alpha^5/pi^3)^0.25 [x\|y\|z]exp(-alpha r^2)
	3084	+ // d = pp->C * pow(pp->A, 1.25) * 1.425410941;
3064	3085	*dp++ = d;
3065	3086	*dp++ = d;
3066	3087	*dp++ = d;
3067	3088	break;
3068	3089	case kGTOType_SP:
3069		- dp++ = pp->C pow(pp->A, 0.75) * 0.71270547;
3070		- d = pp->Csp * pow(pp->A, 1.25) * 1.425410941;
	3090	+ // GNC(0,a,n) * r^n * exp(-a*r^2)
	3091	+ dp++ = d = pp->C sGaussianNormalizationConstant(0, pp->A, n);
	3092	+ //{ printf("type_SP(s): %g %g\n", d, pp->C * pow(pp->A, 0.75) * 0.71270547); }
	3093	+ // GNC(1,a,n) * [x\|y\|z] * r^n * exp(-a*r^2)
	3094	+ d = pp->Csp * sGaussianNormalizationConstant(1, pp->A, n);
	3095	+ //{ printf("type_SP(p): %g %g\n", d, pp->Csp * pow(pp->A, 1.25) * 1.425410941); }
	3096	+ //dp++ = pp->C pow(pp->A, 0.75) * 0.71270547;
	3097	+ //d = pp->Csp * pow(pp->A, 1.25) * 1.425410941;
3071	3098	*dp++ = d;
3072	3099	*dp++ = d;
3073	3100	*dp++ = d;
3074	3101	break;
3075	3102	case kGTOType_D:
	3103	+ // GNC(2,a,n) * [xx\|yy\|zz] * r^n * exp(-a*r^2)
	3104	+ // GNC(2,a,n) * sqrt(3) * [xy\|yz\|zx] * r^n * exp(-a*r^2)
	3105	+ d = pp->C * sGaussianNormalizationConstant(2, pp->A, n);
	3106	+ //{ printf("type_D[0-2]: %g %g\n", d, pp->C * pow(pp->A, 1.75) * 1.645922781); }
	3107	+ //{ printf("type_D[3-5]: %g %g\n", d * sqrt(3), pp->C * pow(pp->A, 1.75) * 2.850821881); }
	3108	+ dp[0] = dp[1] = dp[2] = d;
	3109	+ dp[3] = dp[4] = dp[5] = d * sqrt(3);
3076	3110	// xx\|yy\|zz: (2048 alpha^7/9pi^3)^0.25 [xx\|yy\|zz]exp(-alpha r^2)
3077	3111	// xy\|yz\|zx: (2048 alpha^7/pi^3)^0.25 [xy\|xz\|yz]exp(-alpha r^2)
3078		- d = pp->C * pow(pp->A, 1.75);
3079		- dp[0] = dp[1] = dp[2] = d * 1.645922781;
3080		- dp[3] = dp[4] = dp[5] = d * 2.850821881;
	3112	+ // d = pp->C * pow(pp->A, 1.75);
	3113	+ //dp[0] = dp[1] = dp[2] = d * 1.645922781;
	3114	+ //dp[3] = dp[4] = dp[5] = d * 2.850821881;
3081	3115	dp += 6;
3082	3116	break;
3083	3117	case kGTOType_D5:
	3118	+ // D(0): GNC(2,a,n) * (1/2) * (3zz-rr) * r^n * exp(-a*r^2)
	3119	+ // D(+1): GNC(2,a,n) * sqrt(3) * xz * r^n * exp(-a*r^2)
	3120	+ // D(-1): GNC(2,a,n) * sqrt(3) * yz * r^n * exp(-a*r^2)
	3121	+ // D(+2): GNC(2,a,n) * (sqrt(3)/2) * (xx-yy) * r^n * exp(-a*r^2)
	3122	+ // D(-2): GNC(2,a,n) * sqrt(3) * xy * r^n * exp(-a*r^2)
	3123	+ d = pp->C * sGaussianNormalizationConstant(2, pp->A, n);
	3124	+ //{ printf("type_D5[0]: %g %g\n", d * 0.5, pp->C * pow(pp->A, 1.75) * 0.822961390); }
	3125	+ //{ printf("type_D5[1,2,4]: %g %g\n", d * sqrt(3), pp->C * pow(pp->A, 1.75) * 2.850821881); }
	3126	+ //{ printf("type_D5[3]: %g %g\n", d * sqrt(3) * 0.5, pp->C * pow(pp->A, 1.75) * 1.425410941); }
	3127	+ dp[0] = d * 0.5;
	3128	+ dp[1] = dp[2] = dp[4] = d * sqrt(3);
	3129	+ dp[3] = d * sqrt(3) * 0.5;
3084	3130	// 3zz-rr: (128 alpha^7/9pi^3)^0.25 (3zz-rr)exp(-alpha r^2)
3085	3131	// xy\|yz\|zx: (2048 alpha^7/pi^3)^0.25 [xy\|xz\|yz]exp(-alpha r^2)
3086	3132	// xx-yy: (128 alpha^7/pi^3)^0.25 (xx-yy)exp(-alpha r^2)
3087		- d = pp->C * pow(pp->A, 1.75);
3088		- dp[0] = d * 0.822961390;
3089		- dp[1] = dp[2] = dp[4] = d * 2.850821881;
3090		- dp[3] = d * 1.425410941;
	3133	+ //d = pp->C * pow(pp->A, 1.75);
	3134	+ //dp[0] = d * 0.822961390;
	3135	+ //dp[1] = dp[2] = dp[4] = d * 2.850821881;
	3136	+ //dp[3] = d * 1.425410941;
3091	3137	dp += 5;
3092	3138	break;
3093		- /* TODO: Support F/F7 and G/G9 type orbitals */
	3139	+ case kGTOType_F:
	3140	+ // GNC(3,a,n) * [xxx\|yyy\|zzz] * r^n * exp(-a*r^2)
	3141	+ // GNC(3,a,n) * sqrt(5) * [xyy\|xxy\|xxz\|xzz\|yzz\|yyz] * r^n * exp(-a*r^2)
	3142	+ // GNC(3,a,n) * sqrt(15) * xyz * r^n * exp(-a*r^2)
	3143	+ d = pp->C * sGaussianNormalizationConstant(3, pp->A, n);
	3144	+ dp[0] = dp[1] = dp[2] = d;
	3145	+ dp[3] = dp[4] = dp[5] = dp[6] = dp[7] = dp[8] = d * sqrt(5);
	3146	+ dp[9] = d * sqrt(15);
	3147	+ dp += 10;
	3148	+ break;
	3149	+ case kGTOType_F7:
	3150	+ // F(0): GNC(3,a,n) * (1/2) * (5zzz-3zrr) * r^n * exp(-a*r^2)
	3151	+ // F(+1): GNC(3,a,n) * sqrt(3/8) * (5xzz-xrr) * r^n * exp(-a*r^2)
	3152	+ // F(-1): GNC(3,a,n) * sqrt(3/8) * (5yzz-yrr) * r^n * exp(-a*r^2)
	3153	+ // F(+2): GNC(3,a,n) * sqrt(15/4) * (xxz-yyz) * r^n * exp(-a*r^2)
	3154	+ // F(-2): GNC(3,a,n) * sqrt(15) * xyz * r^n * exp(-a*r^2)
	3155	+ // F(+3): GNC(3,a,n) * sqrt(5/8) * (xxx-3xyy) * r^n * exp(-a*r^2)
	3156	+ // F(-3): GNC(3,a,n) * sqrt(5/8) * (3xxy-yyy) * r^n * exp(-a*r^2)
	3157	+ d = pp->C * sGaussianNormalizationConstant(3, pp->A, n);
	3158	+ dp[0] = d * 0.5;
	3159	+ dp[1] = dp[2] = d * sqrt(3/8.0);
	3160	+ dp[3] = d * sqrt(15/4.0);
	3161	+ dp[4] = d * sqrt(15);
	3162	+ dp[5] = dp[6] = d * sqrt(5/8.0);
	3163	+ dp += 7;
	3164	+ break;
	3165	+ case kGTOType_G:
	3166	+ // GNC(4,a,n) * [xxxx\|yyyy\|zzzz] * exp(-a*r^2)
	3167	+ // GNC(4,a,n) * sqrt(7) * [xxxy\|xxxz\|yyyx\|yyyz\|zzzx\|zzzy] * exp(-a*r^2)
	3168	+ // GNC(4,a,n) * sqrt(35/3) * [xxyy\|xxzz\|yyzz] * exp(-a*r^2)
	3169	+ // GNC(4,a,n) * sqrt(35) * [xxyz\|yyzx\|zzxy] * exp(-a*r^2)
	3170	+ d = pp->C * sGaussianNormalizationConstant(4, pp->A, n);
	3171	+ dp[0] = dp[1] = dp[2] = d;
	3172	+ dp[3] = dp[4] = dp[5] = dp[6] = dp[7] = dp[8] = d * sqrt(7);
	3173	+ dp[9] = dp[10] = dp[11] = d * sqrt(35/3.0);
	3174	+ dp[12] = dp[13] = dp[14] = d * sqrt(35);
	3175	+ dp += 15;
	3176	+ break;
	3177	+ case kGTOType_G9:
	3178	+ // G(0): GNC(4,a,n) * (1/8) * (35zzzz-30zzrr+3rrrr) * exp(-a*r^2)
	3179	+ // G(+1): GNC(4,a,n) * sqrt(5/8) * (7xzzz-3xzrr) * exp(-a*r^2)
	3180	+ // G(-1): GNC(4,a,n) * sqrt(5/8) * (7yzzz-3yzrr) * exp(-a*r^2)
	3181	+ // G(+2): GNC(4,a,n) * sqrt(5/16) * (xx-yy)(7zz-rr) * exp(-a*r^2)
	3182	+ // G(-2): GNC(4,a,n) * sqrt(5/4) * (7xyzz-xyrr) * exp(-a*r^2)
	3183	+ // G(+3): GNC(4,a,n) * sqrt(35/8) * (xxxz-3xyyz) * exp(-a*r^2)
	3184	+ // G(-3): GNC(4,a,n) * sqrt(35/8) * (3xxyz-yyyz) * exp(-a*r^2)
	3185	+ // G(+4): GNC(4,a,n) * sqrt(35/64) * (xxxx-6xxyy+yyyy) * exp(-a*r^2)
	3186	+ // G(-4): GNC(4,a,n) * sqrt(35/4) * (xxxy-xyyy) * exp(-a*r^2)
	3187	+ d = pp->C * sGaussianNormalizationConstant(4, pp->A, n);
	3188	+ dp[0] = d * 0.125;
	3189	+ dp[1] = dp[2] = d * sqrt(5/8.0);
	3190	+ dp[3] = d * sqrt(5/16.0);
	3191	+ dp[4] = d * sqrt(5/4.0);
	3192	+ dp[5] = dp[6] = d * sqrt(35/8.0);
	3193	+ dp[7] = d * sqrt(35/64.0);
	3194	+ dp[8] = d * sqrt(35/4.0);
	3195	+ dp += 9;
	3196	+ break;
3094	3197	}
3095	3198	}
3096	3199	}

		@@ -11360,17 +11463,22 @@ sCalcMOPoint(Molecule mp, const BasisSet bset, Int index, const Vector *vp, Do
11360	11463	val = 0.0;
11361	11464	mobasep = bset->mo + (index - 1) * bset->ncomps;
11362	11465	for (i = 0, sp = bset->shells; i < bset->nshells; i++, sp++) {
11363		- pp = bset->priminfos + sp->p_idx;
	11466	+ Double rn;
	11467	+ pp = bset->priminfos + sp->p_idx;
11364	11468	cnp = bset->cns + sp->cn_idx;
11365	11469	if (sp->a_idx >= mp->natoms)
11366	11470	return 0.0; /* This may happen when molecule is edited after setting up MO info */
11367	11471	tmpp = tmp + sp->a_idx * 4;
11368	11472	mop = mobasep + sp->m_idx;
	11473	+ if (sp->add_exp == 0)
	11474	+ rn = 1.0;
	11475	+ else
	11476	+ rn = pow(tmpp[3], sp->add_exp * 0.5);
11369	11477	switch (sp->sym) {
11370	11478	case kGTOType_S: {
11371	11479	tval = 0;
11372	11480	for (j = 0; j < sp->nprim; j++) {
11373		- tval += cnp++ exp(-pp->A * tmpp[3]);
	11481	+ tval += cnp++ rn * exp(-pp->A * tmpp[3]);
11374	11482	pp++;
11375	11483	}
11376	11484	val += mop[0] * tval;

		@@ -11380,7 +11488,7 @@ sCalcMOPoint(Molecule mp, const BasisSet bset, Int index, const Vector *vp, Do
11380	11488	Double x, y, z;
11381	11489	x = y = z = 0;
11382	11490	for (j = 0; j < sp->nprim; j++) {
11383		- tval = exp(-pp->A * tmpp[3]);
	11491	+ tval = rn * exp(-pp->A * tmpp[3]);
11384	11492	x += cnp++ tval;
11385	11493	y += cnp++ tval;
11386	11494	z += cnp++ tval;

		@@ -11396,7 +11504,7 @@ sCalcMOPoint(Molecule mp, const BasisSet bset, Int index, const Vector *vp, Do
11396	11504	Double t, x, y, z;
11397	11505	t = x = y = z = 0;
11398	11506	for (j = 0; j < sp->nprim; j++) {
11399		- tval = exp(-pp->A * tmpp[3]);
	11507	+ tval = rn * exp(-pp->A * tmpp[3]);
11400	11508	t += cnp++ tval;
11401	11509	x += cnp++ tval;
11402	11510	y += cnp++ tval;

		@@ -11414,7 +11522,7 @@ sCalcMOPoint(Molecule mp, const BasisSet bset, Int index, const Vector *vp, Do
11414	11522	Double xx, yy, zz, xy, xz, yz;
11415	11523	xx = yy = zz = xy = xz = yz = 0;
11416	11524	for (j = 0; j < sp->nprim; j++) {
11417		- tval = exp(-pp->A * tmpp[3]);
	11525	+ tval = rn * exp(-pp->A * tmpp[3]);
11418	11526	xx += cnp++ tval;
11419	11527	yy += cnp++ tval;
11420	11528	zz += cnp++ tval;

		@@ -11436,7 +11544,7 @@ sCalcMOPoint(Molecule mp, const BasisSet bset, Int index, const Vector *vp, Do
11436	11544	Double d0, d1p, d1n, d2p, d2n;
11437	11545	d0 = d1p = d1n = d2p = d2n = 0;
11438	11546	for (j = 0; j < sp->nprim; j++) {
11439		- tval = exp(-pp->A * tmpp[3]);
	11547	+ tval = rn * exp(-pp->A * tmpp[3]);
11440	11548	d0 += cnp++ tval;
11441	11549	d1p += cnp++ tval;
11442	11550	d1n += cnp++ tval;

		@@ -11452,7 +11560,148 @@ sCalcMOPoint(Molecule mp, const BasisSet bset, Int index, const Vector *vp, Do
11452	11560	val += d0 + d1p + d1n + d2p + d2n;
11453	11561	break;
11454	11562	}
11455		- /* TODO: Support F/F7 and G/G9 type orbitals */
	11563	+ case kGTOType_F: {
	11564	+ Double xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz;
	11565	+ xxx = yyy = zzz = xyy = xxy = xxz = xzz = yzz = yyz = xyz = 0;
	11566	+ for (j = 0; j < sp->nprim; j++) {
	11567	+ tval = rn * exp(-pp->A * tmpp[3]);
	11568	+ xxx += cnp++ tval;
	11569	+ yyy += cnp++ tval;
	11570	+ zzz += cnp++ tval;
	11571	+ xyy += cnp++ tval;
	11572	+ xxy += cnp++ tval;
	11573	+ xxz += cnp++ tval;
	11574	+ xzz += cnp++ tval;
	11575	+ yzz += cnp++ tval;
	11576	+ yyz += cnp++ tval;
	11577	+ xyz += cnp++ tval;
	11578	+ pp++;
	11579	+ }
	11580	+ xxx = mop[0] tmpp[0] * tmpp[0] * tmpp[0];
	11581	+ yyy = mop[1] tmpp[1] * tmpp[1] * tmpp[1];
	11582	+ zzz = mop[2] tmpp[2] * tmpp[2] * tmpp[2];
	11583	+ xyy = mop[3] tmpp[0] * tmpp[1] * tmpp[1];
	11584	+ xxy = mop[4] tmpp[0] * tmpp[0] * tmpp[1];
	11585	+ xxz = mop[5] tmpp[0] * tmpp[0] * tmpp[2];
	11586	+ xzz = mop[6] tmpp[0] * tmpp[2] * tmpp[2];
	11587	+ yzz = mop[7] tmpp[1] * tmpp[2] * tmpp[2];
	11588	+ yyz = mop[8] tmpp[1] * tmpp[1] * tmpp[2];
	11589	+ xyz = mop[9] tmpp[0] * tmpp[1] * tmpp[2];
	11590	+ val += xxx + yyy + zzz + xyy + xxy + xxz + xzz + yzz + yyz + xyz;
	11591	+ break;
	11592	+ }
	11593	+ case kGTOType_F7: {
	11594	+ Double f0, f1p, f1n, f2p, f2n, f3p, f3n;
	11595	+ f0 = f1p = f1n = f2p = f2n = f3p = f3n = 0;
	11596	+ for (j = 0; j < sp->nprim; j++) {
	11597	+ tval = rn * exp(-pp->A * tmpp[3]);
	11598	+ f0 += cnp++ tval;
	11599	+ f1p += cnp++ tval;
	11600	+ f1n += cnp++ tval;
	11601	+ f2p += cnp++ tval;
	11602	+ f2n += cnp++ tval;
	11603	+ f3p += cnp++ tval;
	11604	+ f3n += cnp++ tval;
	11605	+ pp++;
	11606	+ }
	11607	+ // F(0): GNC(3,a,n) * (1/2) * (5zzz-3zrr) * r^n * exp(-a*r^2)
	11608	+ // F(+1): GNC(3,a,n) * sqrt(3/8) * (5xzz-xrr) * r^n * exp(-a*r^2)
	11609	+ // F(-1): GNC(3,a,n) * sqrt(3/8) * (5yzz-yrr) * r^n * exp(-a*r^2)
	11610	+ // F(+2): GNC(3,a,n) * sqrt(15/4) * (xxz-yyz) * r^n * exp(-a*r^2)
	11611	+ // F(-2): GNC(3,a,n) * sqrt(15) * xyz * r^n * exp(-a*r^2)
	11612	+ // F(+3): GNC(3,a,n) * sqrt(5/8) * (xxx-3xyy) * r^n * exp(-a*r^2)
	11613	+ // F(-3): GNC(3,a,n) * sqrt(5/8) * (3xxy-yyy) * r^n * exp(-a*r^2)
	11614	+ f0 = mop[0] tmpp[2] * (5 * tmpp[2] * tmpp[2] - 3 * tmpp[3]);
	11615	+ f1p = mop[1] tmpp[0] * (5 * tmpp[2] * tmpp[2] - tmpp[3]);
	11616	+ f1n = mop[2] tmpp[1] * (5 * tmpp[2] * tmpp[2] - tmpp[3]);
	11617	+ f2p = mop[3] tmpp[2] * (tmpp[0] * tmpp[0] - tmpp[1] * tmpp[1]);
	11618	+ f2n = mop[4] tmpp[0] * tmpp[1] * tmpp[2];
	11619	+ f3p = mop[5] tmpp[0] * (tmpp[0] * tmpp[0] - 3 * tmpp[1] * tmpp[1]);
	11620	+ f3n = mop[6] tmpp[1] * (3 * tmpp[0] * tmpp[0] - tmpp[2] * tmpp[2]);
	11621	+ val += f0 + f1p + f1n + f2p + f2n + f3p + f3n;
	11622	+ break;
	11623	+ }
	11624	+ case kGTOType_G: {
	11625	+ Double xxxx, yyyy, zzzz, xxxy, xxxz, yyyx, yyyz, zzzx, zzzy, xxyy, xxzz, yyzz, xxyz, yyxz, zzxy;
	11626	+ xxxx = yyyy = zzzz = xxxy = xxxz = yyyx = yyyz = zzzx = zzzy = xxyy = xxzz = yyzz = xxyz = yyxz = zzxy = 0;
	11627	+ for (j = 0; j < sp->nprim; j++) {
	11628	+ tval = rn * exp(-pp->A * tmpp[3]);
	11629	+ xxxx += cnp++ tval;
	11630	+ yyyy += cnp++ tval;
	11631	+ zzzz += cnp++ tval;
	11632	+ xxxy += cnp++ tval;
	11633	+ xxxz += cnp++ tval;
	11634	+ yyyx += cnp++ tval;
	11635	+ yyyz += cnp++ tval;
	11636	+ zzzx += cnp++ tval;
	11637	+ zzzy += cnp++ tval;
	11638	+ xxyy += cnp++ tval;
	11639	+ xxzz += cnp++ tval;
	11640	+ yyzz += cnp++ tval;
	11641	+ xxyz += cnp++ tval;
	11642	+ yyxz += cnp++ tval;
	11643	+ zzxy += cnp++ tval;
	11644	+ pp++;
	11645	+ }
	11646	+ xxxx = mop[0] tmpp[0] * tmpp[0] * tmpp[0] * tmpp[0];
	11647	+ yyyy = mop[1] tmpp[1] * tmpp[1] * tmpp[1] * tmpp[1];
	11648	+ zzzz = mop[2] tmpp[2] * tmpp[2] * tmpp[2] * tmpp[2];
	11649	+ xxxy = mop[3] tmpp[0] * tmpp[0] * tmpp[0] * tmpp[1];
	11650	+ xxxz = mop[4] tmpp[0] * tmpp[0] * tmpp[0] * tmpp[2];
	11651	+ yyyx = mop[5] tmpp[1] * tmpp[1] * tmpp[1] * tmpp[0];
	11652	+ yyyz = mop[6] tmpp[1] * tmpp[1] * tmpp[1] * tmpp[2];
	11653	+ zzzx = mop[7] tmpp[2] * tmpp[2] * tmpp[2] * tmpp[0];
	11654	+ zzzy = mop[8] tmpp[2] * tmpp[2] * tmpp[2] * tmpp[1];
	11655	+ xxyy = mop[9] tmpp[0] * tmpp[0] * tmpp[1] * tmpp[1];
	11656	+ xxzz = mop[10] tmpp[0] * tmpp[0] * tmpp[2] * tmpp[2];
	11657	+ yyzz = mop[11] tmpp[1] * tmpp[1] * tmpp[2] * tmpp[2];
	11658	+ xxyz = mop[12] tmpp[0] * tmpp[0] * tmpp[1] * tmpp[2];
	11659	+ yyxz = mop[13] tmpp[1] * tmpp[1] * tmpp[0] * tmpp[2];
	11660	+ zzxy = mop[14] tmpp[2] * tmpp[2] * tmpp[0] * tmpp[1];
	11661	+ val += xxxx + yyyy + zzzz + xxxy + xxxz + yyyx + yyyz + zzzx + zzzy + xxyy + xxzz + yyzz + xxyz + yyxz + zzxy;
	11662	+ break;
	11663	+ }
	11664	+ case kGTOType_G9: {
	11665	+ Double g0, g1p, g1n, g2p, g2n, g3p, g3n, g4p, g4n;
	11666	+ Double xx = tmpp[0] * tmpp[0];
	11667	+ Double yy = tmpp[1] * tmpp[1];
	11668	+ Double zz = tmpp[2] * tmpp[2];
	11669	+ Double rr = tmpp[3];
	11670	+ g0 = g1p = g1n = g2p = g2n = g3p = g3n = g4p = g4n = 0;
	11671	+ for (j = 0; j < sp->nprim; j++) {
	11672	+ tval = rn * exp(-pp->A * tmpp[3]);
	11673	+ g0 += cnp++ tval;
	11674	+ g1p += cnp++ tval;
	11675	+ g1n += cnp++ tval;
	11676	+ g2p += cnp++ tval;
	11677	+ g2n += cnp++ tval;
	11678	+ g3p += cnp++ tval;
	11679	+ g3n += cnp++ tval;
	11680	+ g4p += cnp++ tval;
	11681	+ g4n += cnp++ tval;
	11682	+ pp++;
	11683	+ }
	11684	+ // G(0): GNC(4,a,n) * (1/8) * (35zzzz-30zzrr+3rrrr) * r^n * exp(-a*r^2)
	11685	+ // G(+1): GNC(4,a,n) * sqrt(5/8) * (7xzzz-3xzrr) * r^n * exp(-a*r^2)
	11686	+ // G(-1): GNC(4,a,n) * sqrt(5/8) * (7yzzz-3yzrr) * r^n * exp(-a*r^2)
	11687	+ // G(+2): GNC(4,a,n) * sqrt(5/16) * (xx-yy)(7zz-rr) * r^n * exp(-a*r^2)
	11688	+ // G(-2): GNC(4,a,n) * sqrt(5/4) * (7xyzz-xyrr) * r^n * exp(-a*r^2)
	11689	+ // G(+3): GNC(4,a,n) * sqrt(35/8) * (xxxz-3xyyz) * r^n * exp(-a*r^2)
	11690	+ // G(-3): GNC(4,a,n) * sqrt(35/8) * (3xxyz-yyyz) * r^n * exp(-a*r^2)
	11691	+ // G(+4): GNC(4,a,n) * sqrt(35/64) * (xxxx-6xxyy+yyyy) * r^n * exp(-a*r^2)
	11692	+ // G(-4): GNC(4,a,n) * sqrt(35/4) * (xxxy-xyyy) * r^n * exp(-a*r^2)
	11693	+ g0 = mop[0] (35 * zz * zz - 30 * zz * rr + 3 * rr * rr);
	11694	+ g1p = mop[1] tmpp[0] * tmpp[2] * (7 * zz - 3 * rr);
	11695	+ g1n = mop[2] tmpp[1] * tmpp[2] * (7 * zz - 3 * rr);
	11696	+ g2p = mop[3] (xx - yy) * (7 * zz - rr);
	11697	+ g2n = mop[4] tmpp[0] * tmpp[1] * (7 * zz - rr);
	11698	+ g3p = mop[5] tmpp[0] * tmpp[2] * (xx - 3 * yy);
	11699	+ g3n = mop[6] tmpp[1] * tmpp[2] * (3 * xx - yy);
	11700	+ g4p = mop[7] (xx * xx - 6 * xx * yy + yy * yy);
	11701	+ g4n = mop[8] tmpp[0] * tmpp[1] * (xx - yy);
	11702	+ val += g0 + g1p + g1n + g2p + g2n + g3p + g3n + g4p + g4n;
	11703	+ break;
	11704	+ }
11456	11705	}
11457	11706	}
11458	11707	return val;

旧リポジトリブラウザで表示

Molby - 分子モデリングソフトウェア Fork

Molby: コミット

コミットメタ情報

ログメッセージ

変更サマリ

差分

Molby - 分子モデリングソフトウェア
Fork