OSDN > Developer > ktns > Chambre > MolDS > Commit

MolDS
Fork

Dépôt original, Pas de origine de fork

R/O
HTTP
SSH
HTTPS

Commit

Commit MetaInfo

Révision	1c08913bd87dff3cf01e841552c14fece60cd54f (tree)
l'heure	2013-08-20 18:47:29
Auteur	Mikiya Fujii <mikiya.fujii@gmai...>
Commiter	Mikiya Fujii

Message de Log

Refactoring: pow is replaced. #31850

git-svn-id: https://svn.sourceforge.jp/svnroot/molds/trunk@1476 1136aad2-a195-0410-b898-f5ea1d11b9d8

Change Summary

modified: src/cndo/Cndo2.cpp (diff)
modified: src/mndo/Mndo.cpp (diff)

Modification

--- a/src/cndo/Cndo2.cpp

+++ b/src/cndo/Cndo2.cpp

		@@ -356,7 +356,7 @@ double Cndo2::GetDiatomCoreRepulsion1stDerivative(int indexAtomA, int indexAtomB
356	356	double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
357	357	value = atomA.GetCoreCharge()*atomB.GetCoreCharge();
358	358	value *= (atomA.GetXyz()[axisA] - atomB.GetXyz()[axisA])/distance;
359		- value *= -1.0/pow(distance,2.0);
	359	+ value = -1.0/(distancedistance);
360	360	return value;
361	361	}
362	362

		@@ -393,7 +393,8 @@ double Cndo2::GetVdwDampingValue1stDerivative(double vdWDistance, double distanc
393	393	double dampingFactor = Parameters::GetInstance()->GetVdWDampingFactorSCF();
394	394	return (dampingFactor/vdWDistance)
395	395	exp(-1.0dampingFactor*(distance/vdWDistance - 1.0))
396		- pow(1.0+exp(-1.0dampingFactor*(distance/vdWDistance - 1.0)),-2.0);
	396	+ /(1.0+exp(-1.0dampingFactor(distance/vdWDistance - 1.0)))
	397	+ /(1.0+exp(-1.0dampingFactor(distance/vdWDistance - 1.0)));
397	398	}
398	399
399	400	// See damping function in (2) in [G_2004] ((11) in [G_2006])

		@@ -402,8 +403,8 @@ double Cndo2::GetVdwDampingValue2ndDerivative(double vdWDistance, double distanc
402	403	double exponent = -1.0dampingFactor(distance/vdWDistance - 1.0);
403	404	double pre = dampingFactor/vdWDistance;
404	405	double dominator = 1.0+exp(exponent);
405		- return 2.0pow(dominator,-3.0)prepreexp(2.0*exponent)
406		- - pow(dominator,-2.0)prepre*exp( exponent);
	406	+ return 2.0prepreexp(2.0exponent)/(dominatordominatordominator)
	407	+ - prepreexp( exponent)/(dominator*dominator);
407	408	}
408	409
409	410	// See (2) in [G_2004] ((11) in [G_2006])

		@@ -416,7 +417,8 @@ double Cndo2::GetDiatomVdWCorrectionEnergy(int indexAtomA, int indexAtomB) const
416	417	/(atomA.GetVdWCoefficient()+atomB.GetVdWCoefficient());
417	418	double damping = this->GetVdwDampingValue(vdWDistance, distance);
418	419	double scalingFactor = Parameters::GetInstance()->GetVdWScalingFactorSCF();
419		- return -1.0scalingFactorvdWCoefficientspow(distance,-6.0)damping;
	420	+ return -1.0scalingFactorvdWCoefficients*damping
	421	+ /(distancedistancedistancedistancedistance*distance);
420	422	}
421	423
422	424	// First derivative of the vdW correction related to the coordinate of atom A.

		@@ -433,7 +435,9 @@ double Cndo2::GetDiatomVdWCorrection1stDerivative(int indexAtomA, int indexAtomB
433	435	double damping = this->GetVdwDampingValue(vdWDistance, distance);
434	436	double damping1stDerivative = this->GetVdwDampingValue1stDerivative(vdWDistance, distance);
435	437	double value=0.0;
436		- value += 6.0pow(distance,-7.0)damping - pow(distance,-6.0)*damping1stDerivative;
	438	+ double tmp = distancedistancedistancedistancedistance*distance;
	439	+ value += 6.0damping/(tmpdistance)
	440	+ -damping1stDerivative/tmp;
437	441	value *= vdWCoefficients;
438	442	value *= Parameters::GetInstance()->GetVdWScalingFactorSCF();
439	443	value *= (atomA.GetXyz()[axisA] - atomB.GetXyz()[axisA])/distance;

		@@ -461,24 +465,25 @@ double Cndo2::GetDiatomVdWCorrection2ndDerivative(int indexAtomA,
461	465	double damping1stDerivative = this->GetVdwDampingValue1stDerivative(vdWDistance, distance);
462	466	double damping2ndDerivative = this->GetVdwDampingValue2ndDerivative(vdWDistance, distance);
463	467
464		- double temp1 = -6.0pow(distance,-7.0)damping
465		- + pow(distance,-6.0)*damping1stDerivative;
466		- double temp2 = 42.0pow(distance,-8.0)damping
467		- -12.0pow(distance,-7.0)damping1stDerivative
468		- + pow(distance,-6.0)*damping2ndDerivative;
	468	+ double dis6 = distancedistancedistancedistancedistance*distance;
	469	+ double tmp1 = -6.0damping /(dis6distance)
	470	+ + damping1stDerivative/dis6;
	471	+ double tmp2 = 42.0damping /(dis6distance*distance)
	472	+ -12.0damping1stDerivative/(dis6distance)
	473	+ + damping2ndDerivative/dis6;
469	474
470	475	double pre1=0.0;
471	476	double pre2=0.0;
472	477	if(axisA1 != axisA2){
473		- pre1 = -dCartesian1*dCartesian2/pow(distance,3.0);
474		- pre2 = dCartesian1*dCartesian2/pow(distance,2.0);
	478	+ pre1 = -dCartesian1dCartesian2/(distancedistance*distance);
	479	+ pre2 = dCartesian1dCartesian2/(distancedistance);
475	480	}
476	481	else{
477		- pre1 = 1.0/distance - dCartesian1*dCartesian1/pow(distance,3.0);
478		- pre2 = pow(dCartesian1/distance,2.0);
	482	+ pre1 = 1.0/distance - dCartesian1dCartesian1/(distancedistance*distance);
	483	+ pre2 = (dCartesian1dCartesian1)/(distancedistance);
479	484	}
480	485
481		- double value= pre1temp1 + pre2temp2;
	486	+ double value= pre1tmp1 + pre2tmp2;
482	487	value = -1.0vdWScalingFacotor*vdWCoefficients;
483	488	return value;
484	489	}

		@@ -1330,7 +1335,8 @@ bool Cndo2::SatisfyConvergenceCriterion(double const* const * oldOrbitalElectron
1330	1335	for(int i=0; i<numberAOs; i++){
1331	1336	try{
1332	1337	for(int j=0; j<numberAOs; j++){
1333		- change += pow(oldOrbitalElectronPopulation[i][j] - orbitalElectronPopulation[i][j], 2.0);
	1338	+ change += (oldOrbitalElectronPopulation[i][j] - orbitalElectronPopulation[i][j])
	1339	+ *(oldOrbitalElectronPopulation[i][j] - orbitalElectronPopulation[i][j]);
1334	1340	}
1335	1341	}
1336	1342	catch(MolDSException ex){

		@@ -1949,9 +1955,10 @@ void Cndo2::CalcCartesianMatrixElementsByGTOExpansion(double& xComponent,
1949	1955	double gaussianExponentA = 0.0;
1950	1956	double gaussianExponentB = 0.0;
1951	1957	double overlapSASB = 0.0;
1952		- double rAB = sqrt( pow(atomA.GetXyz()[XAxis]-atomB.GetXyz()[XAxis], 2.0)
1953		- +pow(atomA.GetXyz()[YAxis]-atomB.GetXyz()[YAxis], 2.0)
1954		- +pow(atomA.GetXyz()[ZAxis]-atomB.GetXyz()[ZAxis], 2.0) );
	1958	+ double dX = atomA.GetXyz()[XAxis] - atomB.GetXyz()[XAxis];
	1959	+ double dY = atomA.GetXyz()[YAxis] - atomB.GetXyz()[YAxis];
	1960	+ double dZ = atomA.GetXyz()[ZAxis] - atomB.GetXyz()[ZAxis];
	1961	+ double rAB = sqrt(dXdX + dYdY + dZ*dZ);
1955	1962	double temp = 0.0;
1956	1963	double tempX = 0.0;
1957	1964	double tempY = 0.0;

		@@ -1966,12 +1973,12 @@ void Cndo2::CalcCartesianMatrixElementsByGTOExpansion(double& xComponent,
1966	1973	shellTypeB,
1967	1974	valenceOrbitalB,
1968	1975	j);
1969		- gaussianExponentA = pow(orbitalExponentA, 2.0) *
	1976	+ gaussianExponentA = (orbitalExponentAorbitalExponentA)
1970	1977	GTOExpansionSTO::GetInstance()->GetExponent(stonG,
1971	1978	shellTypeA,
1972	1979	valenceOrbitalA,
1973	1980	i);
1974		- gaussianExponentB = pow(orbitalExponentB, 2.0) *
	1981	+ gaussianExponentB = (orbitalExponentBorbitalExponentB)
1975	1982	GTOExpansionSTO::GetInstance()->GetExponent(stonG,
1976	1983	shellTypeB,
1977	1984	valenceOrbitalB,

		@@ -2162,8 +2169,8 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2162	2169	double temp2 = gaussianExponentAxyzA[axis] - gaussianExponentAxyzB[axis];
2163	2170	double temp3 = gaussianExponentBxyzA[axis] - gaussianExponentBxyzB[axis];
2164	2171	value = 0.5*(temp1+temp2-temp3);
2165		- value -= temp1temp2temp3*pow(beta,-1.0);
2166		- value = 4.0sqrt(gaussianExponentAgaussianExponentB)pow(beta,-2.0);
	2172	+ value -= temp1temp2temp3/beta;
	2173	+ value = 4.0sqrt(gaussianExponentAgaussianExponentB)/(betabeta);
2167	2174	value *= overlapSASB;
2168	2175	return value;
2169	2176	}

		@@ -2185,9 +2192,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2185	2192	}
2186	2193	double temp1 = gaussianExponentAxyzA[piDirection] - gaussianExponentAxyzB[piDirection];
2187	2194	double temp2 = gaussianExponentBxyzA[piDirection] - gaussianExponentBxyzB[piDirection];
2188		- value = 0.5 - temp1temp2pow(beta,-1.0);
	2195	+ value = 0.5 - temp1*temp2/beta;
2189	2196	value = gaussianExponentAxyzA[axis] + gaussianExponentB*xyzB[axis];
2190		- value = 4.0sqrt(gaussianExponentAgaussianExponentB)pow(beta,-2.0);
	2197	+ value = 4.0sqrt(gaussianExponentAgaussianExponentB)/(betabeta);
2191	2198	value *= overlapSASB;
2192	2199	return value;
2193	2200	}

		@@ -2209,9 +2216,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2209	2216	}
2210	2217	double temp1 = gaussianExponentAxyzA[axis] + gaussianExponentBxyzB[axis];
2211	2218	double temp2 = gaussianExponentAxyzA[axis] - gaussianExponentAxyzB[axis];
2212		- value = 0.5 + temp1temp2pow(beta,-1.0);
	2219	+ value = 0.5 + temp1*temp2/beta;
2213	2220	value = gaussianExponentBxyzA[piDirectionA] - gaussianExponentB*xyzB[piDirectionA];
2214		- value = -4.0sqrt(gaussianExponentAgaussianExponentB)pow(beta,-2.0);
	2221	+ value = -4.0sqrt(gaussianExponentAgaussianExponentB)/(betabeta);
2215	2222	value *= overlapSASB;
2216	2223	return value;
2217	2224	}

		@@ -2233,9 +2240,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2233	2240	}
2234	2241	double temp1 = gaussianExponentAxyzA[axis] + gaussianExponentBxyzB[axis];
2235	2242	double temp2 = gaussianExponentBxyzA[axis] - gaussianExponentBxyzB[axis];
2236		- value = 0.5 - temp1temp2pow(beta,-1.0);
	2243	+ value = 0.5 - temp1*temp2/beta;
2237	2244	value = gaussianExponentAxyzA[piDirectionB] - gaussianExponentA*xyzB[piDirectionB];
2238		- value = 4.0sqrt(gaussianExponentAgaussianExponentB)pow(beta,-2.0);
	2245	+ value = 4.0sqrt(gaussianExponentAgaussianExponentB)/(betabeta);
2239	2246	value *= overlapSASB;
2240	2247	return value;
2241	2248	}

		@@ -2268,7 +2275,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2268	2275	double temp1 = gaussianExponentBxyzA[piDirectionA] - gaussianExponentBxyzB[piDirectionA];
2269	2276	double temp2 = gaussianExponentAxyzA[axis] + gaussianExponentBxyzB[axis];
2270	2277	double temp3 = gaussianExponentAxyzA[piDirectionB] - gaussianExponentAxyzB[piDirectionB];
2271		- value = -4.0sqrt(gaussianExponentAgaussianExponentB)*pow(beta,-3.0);
	2278	+ value = -4.0sqrt(gaussianExponentAgaussianExponentB)/(betabetabeta);
2272	2279	value = temp1temp2*temp3;
2273	2280	value *= overlapSASB;
2274	2281	return value;

		@@ -2317,8 +2324,10 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2317	2324	overlapSASB);
2318	2325	value = 0.5gaussianExponentAdxyz[axis]
2319	2326	-gaussianExponentB*dxyz[axis]
2320		- +pow(gaussianExponentB,2.0)gaussianExponentApow(dxyz[axis],3.0)/beta;
2321		- value = 8.0pow(gaussianExponentA, 1.5)gaussianExponentBpow(beta,-3.0)dxyz[anotherAxis]overlapSASB;
	2327	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA
	2328	+ dxyz[axis]dxyz[axis]*dxyz[axis]/beta;
	2329	+ value = 8.0gaussianExponentA*sqrt(gaussianExponentA)
	2330	+ gaussianExponentBdxyz[anotherAxis]overlapSASB/(betabeta*beta);
2322	2331	value += xyzA[axis]*overlapAOs1;
2323	2332	return value;
2324	2333	}

		@@ -2366,8 +2375,10 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2366	2375	overlapSASB);
2367	2376	value = 0.5gaussianExponentBdxyz[axis]
2368	2377	-gaussianExponentA*dxyz[axis]
2369		- +pow(gaussianExponentA,2.0)gaussianExponentBpow(dxyz[axis],3.0)/beta;
2370		- value = 8.0pow(gaussianExponentB, 1.5)gaussianExponentApow(beta,-3.0)dxyz[anotherAxis]overlapSASB;
	2378	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB
	2379	+ dxyz[axis]dxyz[axis]*dxyz[axis]/beta;
	2380	+ value = 8.0gaussianExponentBsqrt(gaussianExponentB)gaussianExponentA
	2381	+ dxyz[anotherAxis]overlapSASB/(betabetabeta);
2371	2382	value += xyzB[axis]*overlapAOs1;
2372	2383	return value;
2373	2384	}

		@@ -2400,8 +2411,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2400	2411	dxyz[ZAxis],
2401	2412	rAB,
2402	2413	overlapSASB);
2403		- value = 0.5+pow(gaussianExponentB*dxyz[axis], 2.0)/beta;
2404		- value = 8.0pow(gaussianExponentA, 2.5)gaussianExponentBpow(beta, -3.0)*overlapSASB;
	2414	+ value = 0.5+gaussianExponentBgaussianExponentBdxyz[axis]*dxyz[axis]/beta;
	2415	+ value = 8.0gaussianExponentAgaussianExponentAsqrt(gaussianExponentA)
	2416	+ gaussianExponentBoverlapSASB/(betabetabeta);
2405	2417	value = dxyz[anotherAxis1]dxyz[anotherAxis2];
2406	2418	value += xyzA[axis]*overlapAOs1;
2407	2419	return value;

		@@ -2434,8 +2446,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2434	2446	dxyz[ZAxis],
2435	2447	rAB,
2436	2448	overlapSASB);
2437		- value = 0.5+pow(gaussianExponentA*dxyz[axis], 2.0)/beta;
2438		- value = 8.0pow(gaussianExponentB, 2.5)gaussianExponentApow(beta, -3.0)*overlapSASB;
	2449	+ value = 0.5+gaussianExponentAgaussianExponentAdxyz[axis]*dxyz[axis]/beta;
	2450	+ value = 8.0gaussianExponentBgaussianExponentBsqrt(gaussianExponentB)
	2451	+ gaussianExponentAoverlapSASB/(betabetabeta);
2439	2452	value = dxyz[anotherAxis1]dxyz[anotherAxis2];
2440	2453	value += xyzB[axis]*overlapAOs1;
2441	2454	return value;

		@@ -2588,11 +2601,12 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2588	2601	dxyz[ZAxis],
2589	2602	rAB,
2590	2603	overlapSASB);
2591		- value = 0.5-2.0gaussianExponentAgaussianExponentB*pow(dxyz[XAxis],2.0)/beta;
2592		- value += 0.5pow(gaussianExponentA,2.0)(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/beta;
2593		- value += pow(gaussianExponentAgaussianExponentBdxyz[XAxis]/beta,2.0)
2594		- *(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0));
2595		- value = 4.0pow(gaussianExponentA,0.5)gaussianExponentBpow(beta,-2.0);
	2604	+ value = 0.5-2.0gaussianExponentAgaussianExponentB(dxyz[XAxis]dxyz[XAxis])/beta;
	2605	+ value += 0.5(gaussianExponentAgaussianExponentA)((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))/beta;
	2606	+ value += (gaussianExponentAgaussianExponentBdxyz[XAxis]
	2607	+ gaussianExponentAgaussianExponentB*dxyz[XAxis])
	2608	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/(betabeta);
	2609	+ value = 4.0sqrt(gaussianExponentA)gaussianExponentB/(betabeta);
2596	2610	value *= overlapSASB;
2597	2611	value += xyzA[axis]*overlapAOs1;
2598	2612	return value;

		@@ -2609,11 +2623,12 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2609	2623	dxyz[ZAxis],
2610	2624	rAB,
2611	2625	overlapSASB);
2612		- value = 0.5-2.0gaussianExponentAgaussianExponentB*pow(dxyz[XAxis],2.0)/beta;
2613		- value += 0.5pow(gaussianExponentB,2.0)(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/beta;
2614		- value += pow(gaussianExponentAgaussianExponentBdxyz[XAxis]/beta,2.0)
2615		- *(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0));
2616		- value = 4.0pow(gaussianExponentB,0.5)gaussianExponentApow(beta,-2.0);
	2626	+ value = 0.5-2.0gaussianExponentAgaussianExponentB(dxyz[XAxis]dxyz[XAxis])/beta;
	2627	+ value += 0.5(gaussianExponentBgaussianExponentB)((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))/beta;
	2628	+ value += (gaussianExponentAgaussianExponentBdxyz[XAxis]
	2629	+ gaussianExponentAgaussianExponentB*dxyz[XAxis])
	2630	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/(betabeta);
	2631	+ value = 4.0sqrt(gaussianExponentB)gaussianExponentA/(betabeta);
2617	2632	value *= overlapSASB;
2618	2633	value += xyzB[axis]*overlapAOs1;
2619	2634	return value;

		@@ -2630,11 +2645,12 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2630	2645	dxyz[ZAxis],
2631	2646	rAB,
2632	2647	overlapSASB);
2633		- value = 0.5-2.0gaussianExponentAgaussianExponentB*pow(dxyz[YAxis],2.0)/beta;
2634		- value += 0.5pow(gaussianExponentA,2.0)(pow(dxyz[YAxis],2.0)-pow(dxyz[XAxis],2.0))/beta;
2635		- value += pow(gaussianExponentAgaussianExponentBdxyz[YAxis]/beta,2.0)
2636		- *(pow(dxyz[YAxis],2.0)-pow(dxyz[XAxis],2.0));
2637		- value = -4.0pow(gaussianExponentA,0.5)gaussianExponentBpow(beta,-2.0);
	2648	+ value = 0.5-2.0gaussianExponentAgaussianExponentB(dxyz[YAxis]dxyz[YAxis])/beta;
	2649	+ value += 0.5(gaussianExponentAgaussianExponentA)((dxyz[YAxis]dxyz[YAxis])-(dxyz[XAxis]*dxyz[XAxis]))/beta;
	2650	+ value += (gaussianExponentAgaussianExponentBdxyz[YAxis]
	2651	+ gaussianExponentAgaussianExponentB*dxyz[YAxis])
	2652	+ ((dxyz[YAxis]dxyz[YAxis])-(dxyz[XAxis]dxyz[XAxis]))/(betabeta);
	2653	+ value = -4.0sqrt(gaussianExponentA)gaussianExponentB/(betabeta);
2638	2654	value *= overlapSASB;
2639	2655	value += xyzA[axis]*overlapAOs1;
2640	2656	return value;

		@@ -2651,11 +2667,12 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2651	2667	dxyz[ZAxis],
2652	2668	rAB,
2653	2669	overlapSASB);
2654		- value = 0.5-2.0gaussianExponentAgaussianExponentB*pow(dxyz[YAxis],2.0)/beta;
2655		- value += 0.5pow(gaussianExponentB,2.0)(pow(dxyz[YAxis],2.0)-pow(dxyz[XAxis],2.0))/beta;
2656		- value += pow(gaussianExponentAgaussianExponentBdxyz[YAxis]/beta,2.0)
2657		- *(pow(dxyz[YAxis],2.0)-pow(dxyz[XAxis],2.0));
2658		- value = -4.0pow(gaussianExponentB,0.5)gaussianExponentApow(beta,-2.0);
	2670	+ value = 0.5-2.0gaussianExponentAgaussianExponentB(dxyz[YAxis]dxyz[YAxis])/beta;
	2671	+ value += 0.5(gaussianExponentBgaussianExponentB)((dxyz[YAxis]dxyz[YAxis])-(dxyz[XAxis]*dxyz[XAxis]))/beta;
	2672	+ value += (gaussianExponentAgaussianExponentBdxyz[YAxis]
	2673	+ gaussianExponentAgaussianExponentB*dxyz[YAxis])
	2674	+ ((dxyz[YAxis]dxyz[YAxis])-(dxyz[XAxis]dxyz[XAxis]))/(betabeta);
	2675	+ value = -4.0sqrt(gaussianExponentB)gaussianExponentA/(betabeta);
2659	2676	value *= overlapSASB;
2660	2677	value += xyzB[axis]*overlapAOs1;
2661	2678	return value;

		@@ -2672,11 +2689,12 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2672	2689	dxyz[ZAxis],
2673	2690	rAB,
2674	2691	overlapSASB);
2675		- value = 0.5*(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0));
2676		- value += pow(gaussianExponentB*dxyz[ZAxis],2.0)
2677		- *(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))
	2692	+ value = 0.5((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]));
	2693	+ value += gaussianExponentBgaussianExponentBdxyz[ZAxis]*dxyz[ZAxis]
	2694	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))
2678	2695	/beta;
2679		- value = 4.0pow(gaussianExponentA,2.5)gaussianExponentBpow(beta,-3.0);
	2696	+ value = 4.0gaussianExponentAgaussianExponentAsqrt(gaussianExponentA)
	2697	+ gaussianExponentB/(betabeta*beta);
2680	2698	value *= overlapSASB;
2681	2699	value += xyzA[axis]*overlapAOs1;
2682	2700	return value;

		@@ -2693,11 +2711,11 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2693	2711	dxyz[ZAxis],
2694	2712	rAB,
2695	2713	overlapSASB);
2696		- value = 0.5*(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0));
	2714	+ value = 0.5((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]));
2697	2715	value += pow(gaussianExponentA*dxyz[ZAxis],2.0)
2698		- *(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))
	2716	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))
2699	2717	/beta;
2700		- value = 4.0pow(gaussianExponentB,2.5)gaussianExponentApow(beta,-3.0);
	2718	+ value = 4.0pow(gaussianExponentB,2.5)gaussianExponentA/(betabeta*beta);
2701	2719	value *= overlapSASB;
2702	2720	value += xyzA[axis]*overlapAOs1;
2703	2721	return value;

		@@ -2716,12 +2734,12 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2716	2734	rAB,
2717	2735	overlapSASB);
2718	2736	value = -0.5
2719		- +2.0gaussianExponentAgaussianExponentB*pow(dxyz[axis],2.0)/beta
2720		- +0.5*pow(gaussianExponentA,2.0)
2721		- (2.0pow(dxyz[ZAxis],2.0)-pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/beta
	2737	+ +2.0gaussianExponentAgaussianExponentB(dxyz[axis]dxyz[axis])/beta
	2738	+ +0.5(gaussianExponentAgaussianExponentA)
	2739	+ (2.0(dxyz[ZAxis]dxyz[ZAxis])-(dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))/beta
2722	2740	+pow(gaussianExponentAgaussianExponentBdxyz[axis]/beta,2.0)
2723		- (2.0pow(dxyz[ZAxis],2.0)-pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0));
2724		- value = 4.0pow(gaussianExponentA,0.5)gaussianExponentBpow(beta,-2.0)/sqrt(3.0);
	2741	+ (2.0(dxyz[ZAxis]dxyz[ZAxis])-(dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]));
	2742	+ value = 4.0sqrt(gaussianExponentA)gaussianExponentB/(betabeta*sqrt(3.0));
2725	2743	value *= overlapSASB;
2726	2744	value += xyzA[axis]*overlapAOs1;
2727	2745	return value;

		@@ -2740,12 +2758,12 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2740	2758	rAB,
2741	2759	overlapSASB);
2742	2760	value = -0.5
2743		- +2.0gaussianExponentAgaussianExponentB*pow(dxyz[axis],2.0)/beta
2744		- +0.5*pow(gaussianExponentB,2.0)
2745		- (2.0pow(dxyz[ZAxis],2.0)-pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/beta
	2761	+ +2.0gaussianExponentAgaussianExponentB(dxyz[axis]dxyz[axis])/beta
	2762	+ +0.5(gaussianExponentBgaussianExponentB)
	2763	+ (2.0(dxyz[ZAxis]dxyz[ZAxis])-(dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))/beta
2746	2764	+pow(gaussianExponentAgaussianExponentBdxyz[axis]/beta,2.0)
2747		- (2.0pow(dxyz[ZAxis],2.0)-pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0));
2748		- value = 4.0pow(gaussianExponentB,0.5)gaussianExponentApow(beta,-2.0)/sqrt(3.0);
	2765	+ (2.0(dxyz[ZAxis]dxyz[ZAxis])-(dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]));
	2766	+ value = 4.0sqrt(gaussianExponentB)gaussianExponentA/(betabeta*sqrt(3.0));
2749	2767	value *= overlapSASB;
2750	2768	value += xyzB[axis]*overlapAOs1;
2751	2769	return value;

		@@ -2763,12 +2781,12 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2763	2781	rAB,
2764	2782	overlapSASB);
2765	2783	value = 1.0
2766		- -4.0gaussianExponentAgaussianExponentB*pow(dxyz[axis],2.0)/beta
2767		- +0.5*pow(gaussianExponentA,2.0)
2768		- (2.0pow(dxyz[ZAxis],2.0)-pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/beta
	2784	+ -4.0gaussianExponentAgaussianExponentB(dxyz[axis]dxyz[axis])/beta
	2785	+ +0.5(gaussianExponentAgaussianExponentA)
	2786	+ (2.0(dxyz[ZAxis]dxyz[ZAxis])-(dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))/beta
2769	2787	+pow(gaussianExponentAgaussianExponentBdxyz[axis]/beta,2.0)
2770		- (2.0pow(dxyz[ZAxis],2.0)-pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0));
2771		- value = 4.0pow(gaussianExponentA,0.5)gaussianExponentBpow(beta,-2.0)/sqrt(3.0);
	2788	+ (2.0(dxyz[ZAxis]dxyz[ZAxis])-(dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]));
	2789	+ value = 4.0sqrt(gaussianExponentA)gaussianExponentB/(betabeta*sqrt(3.0));
2772	2790	value *= overlapSASB;
2773	2791	value += xyzA[axis]*overlapAOs1;
2774	2792	return value;

		@@ -2786,12 +2804,12 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2786	2804	rAB,
2787	2805	overlapSASB);
2788	2806	value = 1.0
2789		- -4.0gaussianExponentAgaussianExponentB*pow(dxyz[axis],2.0)/beta
2790		- +0.5*pow(gaussianExponentB,2.0)
2791		- (2.0pow(dxyz[ZAxis],2.0)-pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/beta
	2807	+ -4.0gaussianExponentAgaussianExponentB(dxyz[axis]dxyz[axis])/beta
	2808	+ +0.5(gaussianExponentBgaussianExponentB)
	2809	+ (2.0(dxyz[ZAxis]dxyz[ZAxis])-(dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))/beta
2792	2810	+pow(gaussianExponentAgaussianExponentBdxyz[axis]/beta,2.0)
2793		- (2.0pow(dxyz[ZAxis],2.0)-pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0));
2794		- value = 4.0pow(gaussianExponentB,0.5)gaussianExponentApow(beta,-2.0)/sqrt(3.0);
	2811	+ (2.0(dxyz[ZAxis]dxyz[ZAxis])-(dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]));
	2812	+ value = 4.0sqrt(gaussianExponentB)gaussianExponentA/(betabeta*sqrt(3.0));
2795	2813	value *= overlapSASB;
2796	2814	value += xyzB[axis]*overlapAOs1;
2797	2815	return value;

		@@ -2842,7 +2860,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2842	2860	value = 0.5(gaussianExponentA-gaussianExponentB)dxyz[axis]
2843	2861	+(gaussianExponentA-gaussianExponentB)(gaussianExponentAgaussianExponentB)
2844	2862	dxyz[axis]pow(dxyz[anotherAxis],2.0)/beta;
2845		- value = 8.0(gaussianExponentAgaussianExponentB)pow(beta,-3.0);
	2863	+ value = 8.0(gaussianExponentAgaussianExponentB)/(betabeta*beta);
2846	2864	value *= overlapSASB;
2847	2865	value += axisAverage*overlapAOs1;
2848	2866	return value;

		@@ -2889,8 +2907,8 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2889	2907	overlapSASB);
2890	2908	value = (gaussianExponentA-gaussianExponentB)*dxyz[axis]
2891	2909	-(gaussianExponentA-gaussianExponentB)(gaussianExponentAgaussianExponentB)
2892		- (pow(dxyz[axis],2.0) - pow(dxyz[anotherAxis],2.0))dxyz[axis]/beta;
2893		- value = 4.0(gaussianExponentAgaussianExponentB)pow(beta,-3.0);
	2910	+ ((dxyz[axis]dxyz[axis]) - pow(dxyz[anotherAxis],2.0))*dxyz[axis]/beta;
	2911	+ value = 4.0(gaussianExponentAgaussianExponentB)/(betabeta*beta);
2894	2912	value *= overlapSASB;
2895	2913	value += axisAverage*overlapAOs1;
2896	2914	return value;

		@@ -2911,7 +2929,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2911	2929	overlapSASB);
2912	2930	value = (gaussianExponentA-gaussianExponentB)*dxyz[axis]
2913	2931	-(gaussianExponentA-gaussianExponentB)(gaussianExponentAgaussianExponentB)
2914		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))*dxyz[axis]/beta;
	2932	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]))dxyz[axis]/beta;
2915	2933	value = 4.0(gaussianExponentAgaussianExponentB)pow(beta,-3.0)/3.0;
2916	2934	value *= overlapSASB;
2917	2935	value += axisAverage*overlapAOs1;

		@@ -2932,7 +2950,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
2932	2950	overlapSASB);
2933	2951	value = 2.0(gaussianExponentA-gaussianExponentB)dxyz[axis]
2934	2952	-(gaussianExponentA-gaussianExponentB)(gaussianExponentAgaussianExponentB)
2935		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))*dxyz[axis]/beta;
	2953	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]))dxyz[axis]/beta;
2936	2954	value = 8.0(gaussianExponentAgaussianExponentB)pow(beta,-3.0)/3.0;
2937	2955	value *= overlapSASB;
2938	2956	value += axisAverage*overlapAOs1;

		@@ -3007,7 +3025,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3007	3025	rAB,
3008	3026	overlapSASB);
3009	3027	value = 0.5-gaussianExponentAgaussianExponentBpow(dxyz[anotherAxis1],2.0)/beta;
3010		- value = 8.0pow(gaussianExponentA,2.0)gaussianExponentBpow(beta,-3.0)*dxyz[anotherAxis2];
	3028	+ value = 8.0(gaussianExponentAgaussianExponentA)gaussianExponentBpow(beta,-3.0)dxyz[anotherAxis2];
3011	3029	value *= overlapSASB;
3012	3030	value += axisAverage*overlapAOs1;
3013	3031	return value;

		@@ -3058,7 +3076,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3058	3076	rAB,
3059	3077	overlapSASB);
3060	3078	value = 0.5-gaussianExponentAgaussianExponentBpow(dxyz[anotherAxis1],2.0)/beta;
3061		- value = -8.0pow(gaussianExponentB,2.0)gaussianExponentApow(beta,-3.0)*dxyz[anotherAxis2];
	3079	+ value = -8.0(gaussianExponentBgaussianExponentB)gaussianExponentApow(beta,-3.0)dxyz[anotherAxis2];
3062	3080	value *= overlapSASB;
3063	3081	value += axisAverage*overlapAOs1;
3064	3082	return value;

		@@ -3077,9 +3095,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3077	3095	rAB,
3078	3096	overlapSASB);
3079	3097	value = 0.5*beta
3080		- -pow(gaussianExponentA,2.0)gaussianExponentBpow(dxyz[XAxis],2.0)/beta
3081		- +pow(gaussianExponentB,2.0)*gaussianExponentA
3082		- (pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/(2.0beta);
	3098	+ -(gaussianExponentAgaussianExponentA)gaussianExponentB(dxyz[XAxis]dxyz[XAxis])/beta
	3099	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA
	3100	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/(2.0beta);
3083	3101	value = 8.0gaussianExponentAgaussianExponentBpow(beta,-3.0)*dxyz[YAxis];
3084	3102	value *= overlapSASB;
3085	3103	value += axisAverage*overlapAOs1;

		@@ -3099,9 +3117,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3099	3117	rAB,
3100	3118	overlapSASB);
3101	3119	value = 0.5*beta
3102		- -pow(gaussianExponentB,2.0)gaussianExponentApow(dxyz[XAxis],2.0)/beta
3103		- +pow(gaussianExponentA,2.0)*gaussianExponentB
3104		- (pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/(2.0beta);
	3120	+ -(gaussianExponentBgaussianExponentB)gaussianExponentA(dxyz[XAxis]dxyz[XAxis])/beta
	3121	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB
	3122	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/(2.0beta);
3105	3123	value = -8.0gaussianExponentAgaussianExponentBpow(beta,-3.0)*dxyz[YAxis];
3106	3124	value *= overlapSASB;
3107	3125	value += axisAverage*overlapAOs1;

		@@ -3121,9 +3139,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3121	3139	rAB,
3122	3140	overlapSASB);
3123	3141	value = -0.5*beta
3124		- +pow(gaussianExponentA,2.0)gaussianExponentBpow(dxyz[YAxis],2.0)/beta
3125		- +pow(gaussianExponentB,2.0)*gaussianExponentA
3126		- (pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/(2.0beta);
	3142	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB(dxyz[YAxis]dxyz[YAxis])/beta
	3143	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA
	3144	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/(2.0beta);
3127	3145	value = 8.0gaussianExponentAgaussianExponentBpow(beta,-3.0)*dxyz[XAxis];
3128	3146	value *= overlapSASB;
3129	3147	value += axisAverage*overlapAOs1;

		@@ -3143,9 +3161,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3143	3161	rAB,
3144	3162	overlapSASB);
3145	3163	value = -0.5*beta
3146		- +pow(gaussianExponentB,2.0)gaussianExponentApow(dxyz[YAxis],2.0)/beta
3147		- +pow(gaussianExponentA,2.0)*gaussianExponentB
3148		- (pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/(2.0beta);
	3164	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA(dxyz[YAxis]dxyz[YAxis])/beta
	3165	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB
	3166	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/(2.0beta);
3149	3167	value = -8.0gaussianExponentAgaussianExponentBpow(beta,-3.0)*dxyz[XAxis];
3150	3168	value *= overlapSASB;
3151	3169	value += axisAverage*overlapAOs1;

		@@ -3167,7 +3185,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3167	3185	overlapSASB);
3168	3186	value = 8.0pow(gaussianExponentA,4.0)pow(gaussianExponentB,3.0)
3169	3187	dxyz[XAxis]dxyz[YAxis]*dxyz[ZAxis]
3170		- *(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/pow(beta,5.0);
	3188	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))/pow(beta,5.0);
3171	3189	value *= overlapSASB;
3172	3190	value += xyzB[axis]*overlapAOs1;
3173	3191	return value;

		@@ -3188,7 +3206,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3188	3206	overlapSASB);
3189	3207	value = -8.0pow(gaussianExponentA,3.0)pow(gaussianExponentB,4.0)
3190	3208	dxyz[XAxis]dxyz[YAxis]*dxyz[ZAxis]
3191		- *(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/pow(beta,5.0);
	3209	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]*dxyz[YAxis]))/pow(beta,5.0);
3192	3210	value *= overlapSASB;
3193	3211	value += xyzA[axis]*overlapAOs1;
3194	3212	return value;

		@@ -3207,10 +3225,10 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3207	3225	rAB,
3208	3226	overlapSASB);
3209	3227	value = -0.5*gaussianExponentA
3210		- +pow(gaussianExponentA,2.0)gaussianExponentBpow(dxyz[YAxis],2.0)/beta
3211		- +pow(gaussianExponentB,2.0)*gaussianExponentA
3212		- (pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/(2.0beta);
3213		- value = 8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]*pow(beta,-3.0);
	3228	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB(dxyz[YAxis]dxyz[YAxis])/beta
	3229	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA
	3230	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/(2.0beta);
	3231	+ value = 8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]/(betabetabeta);
3214	3232	value *= overlapSASB;
3215	3233	value += axisAverage*overlapAOs1;
3216	3234	return value;

		@@ -3229,10 +3247,10 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3229	3247	rAB,
3230	3248	overlapSASB);
3231	3249	value = -0.5*gaussianExponentB
3232		- +pow(gaussianExponentB,2.0)gaussianExponentApow(dxyz[YAxis],2.0)/beta
3233		- +pow(gaussianExponentA,2.0)*gaussianExponentB
3234		- (pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/(2.0beta);
3235		- value = -8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]*pow(beta,-3.0);
	3250	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA(dxyz[YAxis]dxyz[YAxis])/beta
	3251	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB
	3252	+ ((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/(2.0beta);
	3253	+ value = -8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]/(betabetabeta);
3236	3254	value *= overlapSASB;
3237	3255	value += axisAverage*overlapAOs1;
3238	3256	return value;

		@@ -3251,10 +3269,10 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3251	3269	rAB,
3252	3270	overlapSASB);
3253	3271	value = -0.5*gaussianExponentA
3254		- +pow(gaussianExponentA,2.0)gaussianExponentBpow(dxyz[XAxis],2.0)/beta
3255		- +pow(gaussianExponentB,2.0)*gaussianExponentA
3256		- (pow(dxyz[YAxis],2.0)-pow(dxyz[XAxis],2.0))/(2.0beta);
3257		- value = -8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]*pow(beta,-3.0);
	3272	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB(dxyz[XAxis]dxyz[XAxis])/beta
	3273	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA
	3274	+ ((dxyz[YAxis]dxyz[YAxis])-(dxyz[XAxis]dxyz[XAxis]))/(2.0beta);
	3275	+ value = -8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]/(betabetabeta);
3258	3276	value *= overlapSASB;
3259	3277	value += axisAverage*overlapAOs1;
3260	3278	return value;

		@@ -3273,10 +3291,10 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3273	3291	rAB,
3274	3292	overlapSASB);
3275	3293	value = -0.5*gaussianExponentB
3276		- +pow(gaussianExponentB,2.0)gaussianExponentApow(dxyz[YAxis],2.0)/beta
3277		- +pow(gaussianExponentA,2.0)*gaussianExponentB
3278		- (pow(dxyz[YAxis],2.0)-pow(dxyz[XAxis],2.0))/(2.0beta);
3279		- value = 8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]*pow(beta,-3.0);
	3294	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA(dxyz[YAxis]dxyz[YAxis])/beta
	3295	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB
	3296	+ ((dxyz[YAxis]dxyz[YAxis])-(dxyz[XAxis]dxyz[XAxis]))/(2.0beta);
	3297	+ value = 8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]/(betabetabeta);
3280	3298	value *= overlapSASB;
3281	3299	value += axisAverage*overlapAOs1;
3282	3300	return value;

		@@ -3294,8 +3312,8 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3294	3312	dxyz[ZAxis],
3295	3313	rAB,
3296	3314	overlapSASB);
3297		- value = gaussianExponentAgaussianExponentB(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/beta+1.0;
3298		- value = 4.0gaussianExponentApow(gaussianExponentB,2.0)dxyz[YAxis]*pow(beta,-3.0);
	3315	+ value = gaussianExponentAgaussianExponentB((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/beta+1.0;
	3316	+ value = 4.0gaussianExponentA(gaussianExponentBgaussianExponentB)dxyz[YAxis]/(betabeta*beta);
3299	3317	value *= overlapSASB;
3300	3318	value += axisAverage*overlapAOs1;
3301	3319	return value;

		@@ -3313,8 +3331,8 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3313	3331	dxyz[ZAxis],
3314	3332	rAB,
3315	3333	overlapSASB);
3316		- value = gaussianExponentAgaussianExponentB(pow(dxyz[XAxis],2.0)-pow(dxyz[YAxis],2.0))/beta+1.0;
3317		- value = -4.0gaussianExponentBpow(gaussianExponentA,2.0)dxyz[YAxis]*pow(beta,-3.0);
	3334	+ value = gaussianExponentAgaussianExponentB((dxyz[XAxis]dxyz[XAxis])-(dxyz[YAxis]dxyz[YAxis]))/beta+1.0;
	3335	+ value = -4.0gaussianExponentB(gaussianExponentAgaussianExponentA)dxyz[YAxis]/(betabeta*beta);
3318	3336	value *= overlapSASB;
3319	3337	value += axisAverage*overlapAOs1;
3320	3338	return value;

		@@ -3332,8 +3350,8 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3332	3350	dxyz[ZAxis],
3333	3351	rAB,
3334	3352	overlapSASB);
3335		- value = gaussianExponentAgaussianExponentB(pow(dxyz[YAxis],2.0)-pow(dxyz[XAxis],2.0))/beta+1.0;
3336		- value = -4.0gaussianExponentApow(gaussianExponentB,2.0)dxyz[XAxis]*pow(beta,-3.0);
	3353	+ value = gaussianExponentAgaussianExponentB((dxyz[YAxis]dxyz[YAxis])-(dxyz[XAxis]dxyz[XAxis]))/beta+1.0;
	3354	+ value = -4.0gaussianExponentA(gaussianExponentBgaussianExponentB)dxyz[XAxis]/(betabeta*beta);
3337	3355	value *= overlapSASB;
3338	3356	value += axisAverage*overlapAOs1;
3339	3357	return value;

		@@ -3351,8 +3369,8 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3351	3369	dxyz[ZAxis],
3352	3370	rAB,
3353	3371	overlapSASB);
3354		- value = gaussianExponentAgaussianExponentB(pow(dxyz[YAxis],2.0)-pow(dxyz[XAxis],2.0))/beta+1.0;
3355		- value = 4.0gaussianExponentBpow(gaussianExponentA,2.0)dxyz[XAxis]*pow(beta,-3.0);
	3372	+ value = gaussianExponentAgaussianExponentB((dxyz[YAxis]dxyz[YAxis])-(dxyz[XAxis]dxyz[XAxis]))/beta+1.0;
	3373	+ value = 4.0gaussianExponentB(gaussianExponentAgaussianExponentA)dxyz[XAxis]/(betabeta*beta);
3356	3374	value *= overlapSASB;
3357	3375	value += axisAverage*overlapAOs1;
3358	3376	return value;

		@@ -3371,7 +3389,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3371	3389	dxyz[ZAxis],
3372	3390	rAB,
3373	3391	overlapSASB);
3374		- value = 2.0*pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0);
	3392	+ value = 2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]);
3375	3393	value = dxyz[XAxis]dxyz[YAxis]*dxyz[ZAxis];
3376	3394	value = 8.0pow(gaussianExponentA,4.0)*pow(gaussianExponentB,3.0);
3377	3395	value /= sqrt(3.0)*pow(beta,5.0);

		@@ -3393,7 +3411,7 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3393	3411	dxyz[ZAxis],
3394	3412	rAB,
3395	3413	overlapSASB);
3396		- value = 2.0*pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0);
	3414	+ value = 2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]);
3397	3415	value = dxyz[XAxis]dxyz[YAxis]*dxyz[ZAxis];
3398	3416	value = -8.0pow(gaussianExponentB,4.0)*pow(gaussianExponentA,3.0);
3399	3417	value /= sqrt(3.0)*pow(beta,5.0);

		@@ -3422,12 +3440,14 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3422	3440	rAB,
3423	3441	overlapSASB);
3424	3442	value = 0.5*(gaussianExponentB-gaussianExponentA)
3425		- +3.0pow(gaussianExponentA,2.0)gaussianExponentBpow(dxyz[axis],2.0)pow(beta,-2.0)
3426		- +gaussianExponentA*pow(gaussianExponentB,2.0)
3427		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/(2.0*beta)
3428		- +pow(gaussianExponentA,3.0)pow(gaussianExponentB,2.0)pow(dxyz[axis],2.0)
3429		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))*pow(beta,-2.0);
3430		- value = 8.0gaussianExponentAgaussianExponentBdxyz[anotherAxis]/(sqrt(3.0)*pow(beta,3.0));
	3443	+ +3.0(gaussianExponentAgaussianExponentA)*gaussianExponentB
	3444	+ (dxyz[axis]dxyz[axis])/(beta*beta)
	3445	+ +gaussianExponentA(gaussianExponentBgaussianExponentB)
	3446	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]))/(2.0beta)
	3447	+ +pow(gaussianExponentA,3.0)(gaussianExponentBgaussianExponentB)(dxyz[axis]dxyz[axis])
	3448	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))
	3449	+ /(beta*beta);
	3450	+ value = 8.0gaussianExponentAgaussianExponentBdxyz[anotherAxis]/(sqrt(3.0)betabeta*beta);
3431	3451	value *= overlapSASB;
3432	3452	value += xyzB[axis]*overlapAOs1;
3433	3453	return value;

		@@ -3453,12 +3473,13 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3453	3473	rAB,
3454	3474	overlapSASB);
3455	3475	value = 0.5*(gaussianExponentA-gaussianExponentB)
3456		- +3.0pow(gaussianExponentB,2.0)gaussianExponentApow(dxyz[axis],2.0)pow(beta,-2.0)
3457		- +gaussianExponentB*pow(gaussianExponentA,2.0)
3458		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/(2.0*beta)
3459		- +pow(gaussianExponentB,3.0)pow(gaussianExponentA,2.0)pow(dxyz[axis],2.0)
3460		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))*pow(beta,-2.0);
3461		- value = -8.0gaussianExponentAgaussianExponentBdxyz[anotherAxis]/(sqrt(3.0)*pow(beta,3.0));
	3476	+ +3.0(gaussianExponentBgaussianExponentB)*gaussianExponentA
	3477	+ (dxyz[axis]dxyz[axis])/(beta*beta)
	3478	+ +gaussianExponentB(gaussianExponentAgaussianExponentA)
	3479	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]))/(2.0beta)
	3480	+ +pow(gaussianExponentB,3.0)(gaussianExponentAgaussianExponentA)(dxyz[axis]dxyz[axis])
	3481	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]))/(betabeta);
	3482	+ value = -8.0gaussianExponentAgaussianExponentBdxyz[anotherAxis]/(sqrt(3.0)(betabeta*beta));
3462	3483	value *= overlapSASB;
3463	3484	value += xyzA[axis]*overlapAOs1;
3464	3485	return value;

		@@ -3477,13 +3498,13 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3477	3498	rAB,
3478	3499	overlapSASB);
3479	3500	value = gaussianExponentB-0.5*gaussianExponentA
3480		- +gaussianExponentA*pow(gaussianExponentB,2.0)
3481		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/(2.0*beta)
	3501	+ +gaussianExponentA(gaussianExponentBgaussianExponentB)
	3502	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]))/(2.0beta)
3482	3503	+pow(gaussianExponentA,3.0)pow(gaussianExponentBdxyz[axis],2.0)
3483		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))
3484		- *pow(beta,-2.0);
	3504	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))
	3505	+ /(beta*beta);
3485	3506	value = 8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]
3486		- /(sqrt(3.0)*pow(beta,3.0));
	3507	+ /(sqrt(3.0)(betabeta*beta));
3487	3508	value *= overlapSASB;
3488	3509	value += xyzB[axis]*overlapAOs1;
3489	3510	return value;

		@@ -3502,13 +3523,13 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3502	3523	rAB,
3503	3524	overlapSASB);
3504	3525	value = gaussianExponentA-0.5*gaussianExponentB
3505		- +gaussianExponentB*pow(gaussianExponentA,2.0)
3506		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/(2.0*beta)
3507		- +pow(gaussianExponentB,3.0)pow(gaussianExponentAdxyz[axis],2.0)
3508		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))
3509		- *pow(beta,-2.0);
	3526	+ +gaussianExponentB(gaussianExponentAgaussianExponentA)
	3527	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]))/(2.0beta)
	3528	+ +pow(gaussianExponentB,3.0)(gaussianExponentAgaussianExponentAdxyz[axis]dxyz[axis])
	3529	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))
	3530	+ /(beta*beta);
3510	3531	value = -8.0gaussianExponentAgaussianExponentBdxyz[ZAxis]
3511		- /(sqrt(3.0)*pow(beta,3.0));
	3532	+ /(sqrt(3.0)(betabeta*beta));
3512	3533	value *= overlapSASB;
3513	3534	value += xyzA[axis]*overlapAOs1;
3514	3535	return value;

		@@ -3534,14 +3555,14 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3534	3555	rAB,
3535	3556	overlapSASB);
3536	3557	value = gaussianExponentA-0.5*gaussianExponentB
3537		- -3.0pow(gaussianExponentA,2.0)gaussianExponentB*pow(dxyz[axis]/beta,2.0)
3538		- +gaussianExponentA*pow(gaussianExponentB,2.0)
3539		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/(2.0*beta)
	3558	+ -3.0(gaussianExponentAgaussianExponentA)gaussianExponentBpow(dxyz[axis]/beta,2.0)
	3559	+ +gaussianExponentA(gaussianExponentBgaussianExponentB)
	3560	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]))/(2.0beta)
3540	3561	+pow(gaussianExponentA,3.0)pow(gaussianExponentBdxyz[axis],2.0)
3541		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))
3542		- *pow(beta,-2.0);
	3562	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))
	3563	+ /(beta*beta);
3543	3564	value = 8.0gaussianExponentAgaussianExponentBdxyz[anotherAxis]
3544		- /(sqrt(3.0)*pow(beta,3.0));
	3565	+ /(sqrt(3.0)(betabeta*beta));
3545	3566	value *= overlapSASB;
3546	3567	value += xyzB[axis]*overlapAOs1;
3547	3568	return value;

		@@ -3567,14 +3588,15 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3567	3588	rAB,
3568	3589	overlapSASB);
3569	3590	value = gaussianExponentB-0.5*gaussianExponentA
3570		- -3.0pow(gaussianExponentB,2.0)gaussianExponentA*pow(dxyz[axis]/beta,2.0)
3571		- +gaussianExponentB*pow(gaussianExponentA,2.0)
3572		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/(2.0*beta)
3573		- +pow(gaussianExponentB,3.0)pow(gaussianExponentAdxyz[axis],2.0)
3574		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))
3575		- *pow(beta,-2.0);
	3591	+ -3.0(gaussianExponentBgaussianExponentB)gaussianExponentAdxyz[axis]dxyz[axis]/(betabeta)
	3592	+ +gaussianExponentB(gaussianExponentAgaussianExponentA)
	3593	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]))/(2.0beta)
	3594	+ +pow(gaussianExponentB,3.0)
	3595	+ (gaussianExponentAgaussianExponentAdxyz[axis]dxyz[axis])
	3596	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))
	3597	+ /(beta*beta);
3576	3598	value = -8.0gaussianExponentAgaussianExponentBdxyz[anotherAxis]
3577		- /(sqrt(3.0)*pow(beta,3.0));
	3599	+ /(sqrt(3.0)(betabeta*beta));
3578	3600	value *= overlapSASB;
3579	3601	value += xyzA[axis]*overlapAOs1;
3580	3602	return value;

		@@ -3593,11 +3615,11 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3593	3615	rAB,
3594	3616	overlapSASB);
3595	3617	value = gaussianExponentB - gaussianExponentA
3596		- +gaussianExponentA*pow(gaussianExponentB,2.0)
3597		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/beta
3598		- +pow(gaussianExponentA,2.0)*gaussianExponentB
3599		- *(pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/beta;
3600		- value = 4.0gaussianExponentAgaussianExponentBdxyz[XAxis]*pow(beta,-3.0)/sqrt(3.0);
	3618	+ +gaussianExponentA(gaussianExponentBgaussianExponentB)
	3619	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))/beta
	3620	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB
	3621	+ ((dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))/beta;
	3622	+ value = 4.0gaussianExponentAgaussianExponentBdxyz[XAxis]/(betabetabeta*sqrt(3.0));
3601	3623	value *= overlapSASB;
3602	3624	value += axisAverage*overlapAOs1;
3603	3625	return value;

		@@ -3616,11 +3638,11 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3616	3638	rAB,
3617	3639	overlapSASB);
3618	3640	value = gaussianExponentA - gaussianExponentB
3619		- +gaussianExponentB*pow(gaussianExponentA,2.0)
3620		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/beta
3621		- +pow(gaussianExponentB,2.0)*gaussianExponentA
3622		- *(pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/beta;
3623		- value = -4.0gaussianExponentAgaussianExponentBdxyz[XAxis]*pow(beta,-3.0)/sqrt(3.0);
	3641	+ +gaussianExponentB(gaussianExponentAgaussianExponentA)
	3642	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))/beta
	3643	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA
	3644	+ ((dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))/beta;
	3645	+ value = -4.0gaussianExponentAgaussianExponentBdxyz[XAxis]/(betabetabeta*sqrt(3.0));
3624	3646	value *= overlapSASB;
3625	3647	value += axisAverage*overlapAOs1;
3626	3648	return value;

		@@ -3639,11 +3661,11 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3639	3661	rAB,
3640	3662	overlapSASB);
3641	3663	value = gaussianExponentB - gaussianExponentA
3642		- +gaussianExponentA*pow(gaussianExponentB,2.0)
3643		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/beta
3644		- +pow(gaussianExponentA,2.0)*gaussianExponentB
3645		- *(pow(dxyz[YAxis],2.0) - pow(dxyz[XAxis],2.0))/beta;
3646		- value = -4.0gaussianExponentAgaussianExponentBdxyz[YAxis]*pow(beta,-3.0)/sqrt(3.0);
	3664	+ +gaussianExponentA(gaussianExponentBgaussianExponentB)
	3665	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))/beta
	3666	+ +(gaussianExponentAgaussianExponentA)gaussianExponentB
	3667	+ ((dxyz[YAxis]dxyz[YAxis]) - (dxyz[XAxis]*dxyz[XAxis]))/beta;
	3668	+ value = -4.0gaussianExponentAgaussianExponentBdxyz[YAxis]/(betabetabeta*sqrt(3.0));
3647	3669	value *= overlapSASB;
3648	3670	value += axisAverage*overlapAOs1;
3649	3671	return value;

		@@ -3662,11 +3684,11 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3662	3684	rAB,
3663	3685	overlapSASB);
3664	3686	value = gaussianExponentA - gaussianExponentB
3665		- +gaussianExponentB*pow(gaussianExponentA,2.0)
3666		- (2.0pow(dxyz[ZAxis],2.0) - pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0))/beta
3667		- +pow(gaussianExponentB,2.0)*gaussianExponentA
3668		- *(pow(dxyz[YAxis],2.0) - pow(dxyz[XAxis],2.0))/beta;
3669		- value = 4.0gaussianExponentAgaussianExponentBdxyz[YAxis]*pow(beta,-3.0)/sqrt(3.0);
	3687	+ +gaussianExponentB(gaussianExponentAgaussianExponentA)
	3688	+ (2.0(dxyz[ZAxis]dxyz[ZAxis]) - (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]*dxyz[YAxis]))/beta
	3689	+ +(gaussianExponentBgaussianExponentB)gaussianExponentA
	3690	+ ((dxyz[YAxis]dxyz[YAxis]) - (dxyz[XAxis]*dxyz[XAxis]))/beta;
	3691	+ value = 4.0gaussianExponentAgaussianExponentBdxyz[YAxis]/(betabetabeta*sqrt(3.0));
3670	3692	value *= overlapSASB;
3671	3693	value += axisAverage*overlapAOs1;
3672	3694	return value;

		@@ -3684,9 +3706,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3684	3706	dxyz[ZAxis],
3685	3707	rAB,
3686	3708	overlapSASB);
3687		- value = pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0);
3688		- value = -8.0pow(gaussianExponentA,3.0)pow(gaussianExponentB,2.0)dxyz[ZAxis];
3689		- value /= sqrt(3.0)*pow(beta,4.0);
	3709	+ value = (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]);
	3710	+ value = -8.0pow(gaussianExponentA,3.0)(gaussianExponentBgaussianExponentB)*dxyz[ZAxis];
	3711	+ value /= sqrt(3.0)betabetabetabeta;
3690	3712	value *= overlapSASB;
3691	3713	value += axisAverage*overlapAOs1;
3692	3714	return value;

		@@ -3704,9 +3726,9 @@ double Cndo2::GetGaussianCartesianMatrix(AtomType atomTypeA,
3704	3726	dxyz[ZAxis],
3705	3727	rAB,
3706	3728	overlapSASB);
3707		- value = pow(dxyz[XAxis],2.0) - pow(dxyz[YAxis],2.0);
3708		- value = 8.0pow(gaussianExponentA,2.0)pow(gaussianExponentB,3.0)dxyz[ZAxis];
3709		- value /= sqrt(3.0)*pow(beta,4.0);
	3729	+ value = (dxyz[XAxis]dxyz[XAxis]) - (dxyz[YAxis]dxyz[YAxis]);
	3730	+ value = 8.0(gaussianExponentAgaussianExponentA)pow(gaussianExponentB,3.0)*dxyz[ZAxis];
	3731	+ value /= sqrt(3.0)betabetabetabeta;
3710	3732	value *= overlapSASB;
3711	3733	value += axisAverage*overlapAOs1;
3712	3734	return value;

		@@ -4001,9 +4023,7 @@ void Cndo2::CalcDiatomicOverlapAOs1stDerivatives(double*** diatomicOverlapAOs1st
4001	4023	double cartesian[CartesianType_end] = {atomA.GetXyz()[XAxis] - atomB.GetXyz()[XAxis],
4002	4024	atomA.GetXyz()[YAxis] - atomB.GetXyz()[YAxis],
4003	4025	atomA.GetXyz()[ZAxis] - atomB.GetXyz()[ZAxis]};
4004		- double R = sqrt( pow(cartesian[XAxis],2.0) +
4005		- pow(cartesian[YAxis],2.0) +
4006		- pow(cartesian[ZAxis],2.0) );
	4026	+ double R = this->molecule->GetDistanceAtoms(atomA, atomB);
4007	4027
4008	4028	this->CalcDiatomicOverlapAOsInDiatomicFrame(tmpDiaOverlapAOsInDiaFrame, atomA, atomB);
4009	4029	this->CalcDiatomicOverlapAOs1stDerivativeInDiatomicFrame(tmpDiaOverlapAOs1stDerivInDiaFrame, atomA, atomB);

		@@ -4133,9 +4153,7 @@ void Cndo2::CalcDiatomicOverlapAOs2ndDerivatives(double**** diatomicOverlapAOs2n
4133	4153	double cartesian[CartesianType_end] = {atomA.GetXyz()[XAxis] - atomB.GetXyz()[XAxis],
4134	4154	atomA.GetXyz()[YAxis] - atomB.GetXyz()[YAxis],
4135	4155	atomA.GetXyz()[ZAxis] - atomB.GetXyz()[ZAxis]};
4136		- double R = sqrt( pow(cartesian[XAxis],2.0) +
4137		- pow(cartesian[YAxis],2.0) +
4138		- pow(cartesian[ZAxis],2.0) );
	4156	+ double R = this->molecule->GetDistanceAtoms(atomA, atomB);
4139	4157
4140	4158	this->CalcDiatomicOverlapAOsInDiatomicFrame(tmpDiaOverlapAOsInDiaFrame, atomA, atomB);
4141	4159	this->CalcDiatomicOverlapAOs1stDerivativeInDiatomicFrame(tmpDiaOverlapAOs1stDerivInDiaFrame, atomA, atomB);

		@@ -4266,13 +4284,13 @@ double Cndo2::Get2ndDerivativeElementFromDistanceDerivatives(double firstDistanc
4266	4284	double rAB) const{
4267	4285	double value=0.0;
4268	4286	if(axisA1 != axisA2){
4269		- value = -1.0pow(rAB, -3.0)firstDistanceDeri;
4270		- value += pow(rAB, -2.0)*secondDistanceDeri;
	4287	+ value = -1.0firstDistanceDeri/(rABrAB*rAB);
	4288	+ value += secondDistanceDeri/(rAB*rAB);
4271	4289	value = cartesian[axisA1]cartesian[axisA2];
4272	4290	}
4273	4291	else{
4274		- value = (pow(rAB,2.0) - pow(cartesian[axisA1],2.0))pow(rAB, -3.0)firstDistanceDeri;
4275		- value += pow(cartesian[axisA1]/rAB, 2.0)*secondDistanceDeri;
	4292	+ value = (rABrAB - cartesian[axisA1]cartesian[axisA1])firstDistanceDeri/(rABrAB*rAB);
	4293	+ value += cartesian[axisA1]cartesian[axisA1]secondDistanceDeri/(rAB*rAB);
4276	4294	}
4277	4295	return value;
4278	4296	}

		@@ -4341,7 +4359,7 @@ double Cndo2::GetOverlapAOsElementByGTOExpansion(const Atom& atomA, int valenceI
4341	4359	double dx = atomA.GetXyz()[XAxis] - atomB.GetXyz()[XAxis];
4342	4360	double dy = atomA.GetXyz()[YAxis] - atomB.GetXyz()[YAxis];
4343	4361	double dz = atomA.GetXyz()[ZAxis] - atomB.GetXyz()[ZAxis];
4344		- double rAB = sqrt( pow(dx, 2.0) + pow(dy, 2.0) + pow(dz,2.0) );
	4362	+ double rAB = sqrt( dxdx + dydy + dz*dz );
4345	4363	ShellType shellTypeA = atomA.GetValenceShellType();
4346	4364	ShellType shellTypeB = atomB.GetValenceShellType();
4347	4365	OrbitalType valenceOrbitalA = atomA.GetValence(valenceIndexA);

		@@ -4366,12 +4384,12 @@ double Cndo2::GetOverlapAOsElementByGTOExpansion(const Atom& atomA, int valenceI
4366	4384	shellTypeB,
4367	4385	valenceOrbitalB,
4368	4386	j);
4369		- gaussianExponentA = pow(orbitalExponentA, 2.0) *
	4387	+ gaussianExponentA = (orbitalExponentAorbitalExponentA)
4370	4388	GTOExpansionSTO::GetInstance()->GetExponent(stonG,
4371	4389	shellTypeA,
4372	4390	valenceOrbitalA,
4373	4391	i);
4374		- gaussianExponentB = pow(orbitalExponentB, 2.0) *
	4392	+ gaussianExponentB = (orbitalExponentBorbitalExponentB)
4375	4393	GTOExpansionSTO::GetInstance()->GetExponent(stonG,
4376	4394	shellTypeB,
4377	4395	valenceOrbitalB,

		@@ -4397,7 +4415,7 @@ double Cndo2::GetGaussianOverlapAOsSASB(double gaussianExponentA,
4397	4415	double value;
4398	4416	double temp1 = 0.0;
4399	4417	double temp2 = 0.0;
4400		- temp1 = 2.0pow(gaussianExponentAgaussianExponentB, 0.5)
	4418	+ temp1 = 2.0sqrt(gaussianExponentAgaussianExponentB)
4401	4419	/(gaussianExponentA+gaussianExponentB);
4402	4420	temp2 = -1.0* gaussianExponentA*gaussianExponentB
4403	4421	/(gaussianExponentA+gaussianExponentB);

		@@ -4445,125 +4463,125 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4445	4463	}
4446	4464
4447	4465	else if(valenceOrbitalA == s && valenceOrbitalB == px){
4448		- value = 2.0gaussianExponentApow(gaussianExponentB, 0.5)*dx;
	4466	+ value = 2.0gaussianExponentAsqrt(gaussianExponentB)*dx;
4449	4467	value /= (gaussianExponentA+gaussianExponentB);
4450	4468	}
4451	4469	else if(valenceOrbitalA == s && valenceOrbitalB == py){
4452		- value = 2.0gaussianExponentApow(gaussianExponentB, 0.5)*dy;
	4470	+ value = 2.0gaussianExponentAsqrt(gaussianExponentB)*dy;
4453	4471	value /= (gaussianExponentA+gaussianExponentB);
4454	4472	}
4455	4473	else if(valenceOrbitalA == s && valenceOrbitalB == pz){
4456		- value = 2.0gaussianExponentApow(gaussianExponentB, 0.5)*dz;
	4474	+ value = 2.0gaussianExponentAsqrt(gaussianExponentB)*dz;
4457	4475	value /= (gaussianExponentA+gaussianExponentB);
4458	4476	}
4459	4477
4460	4478	else if(valenceOrbitalA == px && valenceOrbitalB == s){
4461		- value = -2.0pow(gaussianExponentA, 0.5)gaussianExponentB*dx;
	4479	+ value = -2.0sqrt(gaussianExponentA)gaussianExponentB*dx;
4462	4480	value /= (gaussianExponentA+gaussianExponentB);
4463	4481	}
4464	4482	else if(valenceOrbitalA == py && valenceOrbitalB == s){
4465		- value = -2.0pow(gaussianExponentA, 0.5)gaussianExponentB*dy;
	4483	+ value = -2.0sqrt(gaussianExponentA)gaussianExponentB*dy;
4466	4484	value /= (gaussianExponentA+gaussianExponentB);
4467	4485	}
4468	4486	else if(valenceOrbitalA == pz && valenceOrbitalB == s){
4469		- value = -2.0pow(gaussianExponentA, 0.5)gaussianExponentB*dz;
	4487	+ value = -2.0sqrt(gaussianExponentA)gaussianExponentB*dz;
4470	4488	value /= (gaussianExponentA+gaussianExponentB);
4471	4489	}
4472	4490
4473	4491	else if(valenceOrbitalA == px && valenceOrbitalB == px){
4474	4492	double temp = 0.0;
4475		- temp = -1.0pow(dx,2.0)gaussianExponentA*gaussianExponentB;
	4493	+ temp = -1.0(dxdx)gaussianExponentAgaussianExponentB;
4476	4494	temp /= (gaussianExponentA+gaussianExponentB);
4477	4495	temp += 0.5;
4478		- value = 4.0pow(gaussianExponentAgaussianExponentB, 0.5);
	4496	+ value = 4.0sqrt(gaussianExponentAgaussianExponentB);
4479	4497	value /= (gaussianExponentA+gaussianExponentB);
4480	4498	value *= temp;
4481	4499	}
4482	4500	else if(valenceOrbitalA == px && valenceOrbitalB == py){
4483		- value = -4.0pow(gaussianExponentAgaussianExponentB, 1.5);
4484		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4501	+ value = -4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB);
	4502	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4485	4503	value = dxdy;
4486	4504	}
4487	4505	else if(valenceOrbitalA == px && valenceOrbitalB == pz){
4488		- value = -4.0pow(gaussianExponentAgaussianExponentB, 1.5);
4489		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4506	+ value = -4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB);
	4507	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4490	4508	value = dxdz;
4491	4509	}
4492	4510
4493	4511	else if(valenceOrbitalA == py && valenceOrbitalB == px){
4494		- value = -4.0pow(gaussianExponentAgaussianExponentB, 1.5);
4495		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4512	+ value = -4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB);
	4513	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4496	4514	value = dydx;
4497	4515	}
4498	4516	else if(valenceOrbitalA == py && valenceOrbitalB == py){
4499	4517	double temp = 0.0;
4500		- temp = -1.0pow(dy,2.0)gaussianExponentA*gaussianExponentB;
	4518	+ temp = -1.0(dydy)gaussianExponentAgaussianExponentB;
4501	4519	temp /= (gaussianExponentA+gaussianExponentB);
4502	4520	temp += 0.5;
4503		- value = 4.0pow(gaussianExponentAgaussianExponentB, 0.5);
	4521	+ value = 4.0sqrt(gaussianExponentAgaussianExponentB);
4504	4522	value /= (gaussianExponentA+gaussianExponentB);
4505	4523	value *= temp;
4506	4524	}
4507	4525	else if(valenceOrbitalA == py && valenceOrbitalB == pz){
4508		- value = -4.0pow(gaussianExponentAgaussianExponentB, 1.5);
4509		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4526	+ value = -4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB);
	4527	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4510	4528	value = dydz;
4511	4529	}
4512	4530
4513	4531	else if(valenceOrbitalA == pz && valenceOrbitalB == px){
4514		- value = -4.0pow(gaussianExponentAgaussianExponentB, 1.5);
4515		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4532	+ value = -4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB);
	4533	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4516	4534	value = dzdx;
4517	4535	}
4518	4536	else if(valenceOrbitalA == pz && valenceOrbitalB == py){
4519		- value = -4.0pow(gaussianExponentAgaussianExponentB, 1.5);
4520		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4537	+ value = -4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB);
	4538	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4521	4539	value = dzdy;
4522	4540	}
4523	4541	else if(valenceOrbitalA == pz && valenceOrbitalB == pz){
4524	4542	double temp = 0.0;
4525		- temp = -1.0pow(dz,2.0)gaussianExponentA*gaussianExponentB;
	4543	+ temp = -1.0(dzdz)gaussianExponentAgaussianExponentB;
4526	4544	temp /= (gaussianExponentA+gaussianExponentB);
4527	4545	temp += 0.5;
4528		- value = 4.0pow(gaussianExponentAgaussianExponentB, 0.5);
	4546	+ value = 4.0sqrt(gaussianExponentAgaussianExponentB);
4529	4547	value /= (gaussianExponentA+gaussianExponentB);
4530	4548	value *= temp;
4531	4549	}
4532	4550
4533	4551	else if(valenceOrbitalA == dxy && valenceOrbitalB == s){
4534	4552	value = 4.0*gaussianExponentA;
4535		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4553	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4536	4554	value = gaussianExponentBdx;
4537	4555	value = gaussianExponentBdy;
4538	4556	}
4539	4557	else if(valenceOrbitalA == dyz && valenceOrbitalB == s){
4540	4558	value = 4.0*gaussianExponentA;
4541		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4559	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4542	4560	value = gaussianExponentBdy;
4543	4561	value = gaussianExponentBdz;
4544	4562	}
4545	4563	else if(valenceOrbitalA == dzx && valenceOrbitalB == s){
4546	4564	value = 4.0*gaussianExponentA;
4547		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4565	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4548	4566	value = gaussianExponentBdz;
4549	4567	value = gaussianExponentBdx;
4550	4568	}
4551	4569
4552	4570	else if(valenceOrbitalA == s && valenceOrbitalB == dxy){
4553	4571	value = 4.0*gaussianExponentB;
4554		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4572	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4555	4573	value = gaussianExponentAdx;
4556	4574	value = gaussianExponentAdy;
4557	4575	}
4558	4576	else if(valenceOrbitalA == s && valenceOrbitalB == dyz){
4559	4577	value = 4.0*gaussianExponentB;
4560		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4578	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4561	4579	value = gaussianExponentAdy;
4562	4580	value = gaussianExponentAdz;
4563	4581	}
4564	4582	else if(valenceOrbitalA == s && valenceOrbitalB == dzx){
4565	4583	value = 4.0*gaussianExponentB;
4566		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4584	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4567	4585	value = gaussianExponentAdz;
4568	4586	value = gaussianExponentAdx;
4569	4587	}

		@@ -4573,7 +4591,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4573	4591	double temp2 = gaussianExponentBdxgaussianExponentAdxgaussianExponentB*dy;
4574	4592	temp2 /= gaussianExponentA+gaussianExponentB;
4575	4593	value = temp1 + temp2;
4576		- value = 8.0gaussianExponentA*pow(gaussianExponentB, 0.5);
	4594	+ value = 8.0gaussianExponentA*sqrt(gaussianExponentB);
4577	4595	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4578	4596	}
4579	4597	else if(valenceOrbitalA == dxy && valenceOrbitalB == py){

		@@ -4581,7 +4599,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4581	4599	double temp2 = gaussianExponentBdygaussianExponentAdygaussianExponentB*dx;
4582	4600	temp2 /= gaussianExponentA+gaussianExponentB;
4583	4601	value = temp1 + temp2;
4584		- value = 8.0gaussianExponentA*pow(gaussianExponentB, 0.5);
	4602	+ value = 8.0gaussianExponentA*sqrt(gaussianExponentB);
4585	4603	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4586	4604	}
4587	4605	else if(valenceOrbitalA == dyz && valenceOrbitalB == py){

		@@ -4589,7 +4607,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4589	4607	double temp2 = gaussianExponentBdygaussianExponentAdygaussianExponentB*dz;
4590	4608	temp2 /= gaussianExponentA+gaussianExponentB;
4591	4609	value = temp1 + temp2;
4592		- value = 8.0gaussianExponentA*pow(gaussianExponentB, 0.5);
	4610	+ value = 8.0gaussianExponentA*sqrt(gaussianExponentB);
4593	4611	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4594	4612	}
4595	4613	else if(valenceOrbitalA == dyz && valenceOrbitalB == pz){

		@@ -4597,7 +4615,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4597	4615	double temp2 = gaussianExponentBdzgaussianExponentAdzgaussianExponentB*dy;
4598	4616	temp2 /= gaussianExponentA+gaussianExponentB;
4599	4617	value = temp1 + temp2;
4600		- value = 8.0gaussianExponentA*pow(gaussianExponentB, 0.5);
	4618	+ value = 8.0gaussianExponentA*sqrt(gaussianExponentB);
4601	4619	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4602	4620	}
4603	4621	else if(valenceOrbitalA == dzx && valenceOrbitalB == pz){

		@@ -4605,7 +4623,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4605	4623	double temp2 = gaussianExponentBdzgaussianExponentAdzgaussianExponentB*dx;
4606	4624	temp2 /= gaussianExponentA+gaussianExponentB;
4607	4625	value = temp1 + temp2;
4608		- value = 8.0gaussianExponentA*pow(gaussianExponentB, 0.5);
	4626	+ value = 8.0gaussianExponentA*sqrt(gaussianExponentB);
4609	4627	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4610	4628	}
4611	4629	else if(valenceOrbitalA == dzx && valenceOrbitalB == px){

		@@ -4613,7 +4631,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4613	4631	double temp2 = gaussianExponentBdxgaussianExponentAdxgaussianExponentB*dz;
4614	4632	temp2 /= gaussianExponentA+gaussianExponentB;
4615	4633	value = temp1 + temp2;
4616		- value = 8.0gaussianExponentA*pow(gaussianExponentB, 0.5);
	4634	+ value = 8.0gaussianExponentA*sqrt(gaussianExponentB);
4617	4635	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4618	4636	}
4619	4637

		@@ -4622,7 +4640,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4622	4640	double temp2 = -1.0gaussianExponentAdxgaussianExponentBdxgaussianExponentAdy;
4623	4641	temp2 /= gaussianExponentA+gaussianExponentB;
4624	4642	value = temp1 + temp2;
4625		- value = 8.0gaussianExponentB*pow(gaussianExponentA, 0.5);
	4643	+ value = 8.0gaussianExponentB*sqrt(gaussianExponentA);
4626	4644	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4627	4645	}
4628	4646	else if(valenceOrbitalA == py && valenceOrbitalB == dxy){

		@@ -4630,7 +4648,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4630	4648	double temp2 = -1.0gaussianExponentAdygaussianExponentBdygaussianExponentAdx;
4631	4649	temp2 /= gaussianExponentA+gaussianExponentB;
4632	4650	value = temp1 + temp2;
4633		- value = 8.0gaussianExponentB*pow(gaussianExponentA, 0.5);
	4651	+ value = 8.0gaussianExponentB*sqrt(gaussianExponentA);
4634	4652	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4635	4653	}
4636	4654	else if(valenceOrbitalA == py && valenceOrbitalB == dyz){

		@@ -4638,7 +4656,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4638	4656	double temp2 = -1.0gaussianExponentAdygaussianExponentBdygaussianExponentAdz;
4639	4657	temp2 /= gaussianExponentA+gaussianExponentB;
4640	4658	value = temp1 + temp2;
4641		- value = 8.0gaussianExponentB*pow(gaussianExponentA, 0.5);
	4659	+ value = 8.0gaussianExponentB*sqrt(gaussianExponentA);
4642	4660	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4643	4661	}
4644	4662	else if(valenceOrbitalA == pz && valenceOrbitalB == dyz){

		@@ -4646,7 +4664,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4646	4664	double temp2 = -1.0gaussianExponentAdzgaussianExponentBdzgaussianExponentAdy;
4647	4665	temp2 /= gaussianExponentA+gaussianExponentB;
4648	4666	value = temp1 + temp2;
4649		- value = 8.0gaussianExponentB*pow(gaussianExponentA, 0.5);
	4667	+ value = 8.0gaussianExponentB*sqrt(gaussianExponentA);
4650	4668	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4651	4669	}
4652	4670	else if(valenceOrbitalA == pz && valenceOrbitalB == dzx){

		@@ -4654,7 +4672,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4654	4672	double temp2 = -1.0gaussianExponentAdzgaussianExponentBdzgaussianExponentAdx;
4655	4673	temp2 /= gaussianExponentA+gaussianExponentB;
4656	4674	value = temp1 + temp2;
4657		- value = 8.0gaussianExponentB*pow(gaussianExponentA, 0.5);
	4675	+ value = 8.0gaussianExponentB*sqrt(gaussianExponentA);
4658	4676	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4659	4677	}
4660	4678	else if(valenceOrbitalA == px && valenceOrbitalB == dzx){

		@@ -4662,170 +4680,170 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4662	4680	double temp2 = -1.0gaussianExponentAdxgaussianExponentBdxgaussianExponentAdz;
4663	4681	temp2 /= gaussianExponentA+gaussianExponentB;
4664	4682	value = temp1 + temp2;
4665		- value = 8.0gaussianExponentB*pow(gaussianExponentA, 0.5);
	4683	+ value = 8.0gaussianExponentB*sqrt(gaussianExponentA);
4666	4684	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4667	4685	}
4668	4686
4669	4687	else if(valenceOrbitalA == dxy && valenceOrbitalB == pz){
4670		- value = 8.0gaussianExponentApow(gaussianExponentB, 0.5);
	4688	+ value = 8.0gaussianExponentAsqrt(gaussianExponentB);
4671	4689	value /= pow(gaussianExponentA+gaussianExponentB, 3.0);
4672	4690	value = gaussianExponentBdxgaussianExponentBdygaussianExponentAdz;
4673	4691	}
4674	4692	else if(valenceOrbitalA == dyz && valenceOrbitalB == px){
4675		- value = 8.0gaussianExponentApow(gaussianExponentB, 0.5);
	4693	+ value = 8.0gaussianExponentAsqrt(gaussianExponentB);
4676	4694	value /= pow(gaussianExponentA+gaussianExponentB, 3.0);
4677	4695	value = gaussianExponentBdygaussianExponentBdzgaussianExponentAdx;
4678	4696	}
4679	4697	else if(valenceOrbitalA == dzx && valenceOrbitalB == py){
4680		- value = 8.0gaussianExponentApow(gaussianExponentB, 0.5);
	4698	+ value = 8.0gaussianExponentAsqrt(gaussianExponentB);
4681	4699	value /= pow(gaussianExponentA+gaussianExponentB, 3.0);
4682	4700	value = gaussianExponentBdzgaussianExponentBdxgaussianExponentAdy;
4683	4701	}
4684	4702
4685	4703	else if(valenceOrbitalA == pz && valenceOrbitalB == dxy){
4686		- value = -8.0gaussianExponentBpow(gaussianExponentA, 0.5);
	4704	+ value = -8.0gaussianExponentBsqrt(gaussianExponentA);
4687	4705	value /= pow(gaussianExponentA+gaussianExponentB, 3.0);
4688	4706	value = gaussianExponentAdxgaussianExponentAdygaussianExponentBdz;
4689	4707	}
4690	4708	else if(valenceOrbitalA == px && valenceOrbitalB == dyz){
4691		- value = -8.0gaussianExponentBpow(gaussianExponentA, 0.5);
	4709	+ value = -8.0gaussianExponentBsqrt(gaussianExponentA);
4692	4710	value /= pow(gaussianExponentA+gaussianExponentB, 3.0);
4693	4711	value = gaussianExponentAdygaussianExponentAdzgaussianExponentBdx;
4694	4712	}
4695	4713	else if(valenceOrbitalA == py && valenceOrbitalB == dzx){
4696		- value = -8.0gaussianExponentBpow(gaussianExponentA, 0.5);
	4714	+ value = -8.0gaussianExponentBsqrt(gaussianExponentA);
4697	4715	value /= pow(gaussianExponentA+gaussianExponentB, 3.0);
4698	4716	value = gaussianExponentAdzgaussianExponentAdxgaussianExponentBdy;
4699	4717	}
4700	4718
4701	4719	else if(valenceOrbitalA == dxxyy && valenceOrbitalB == s){
4702	4720	value = 2.0*gaussianExponentA;
4703		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
4704		- value = pow(gaussianExponentBdx, 2.0) - pow(gaussianExponentB*dy, 2.0);
	4721	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
	4722	+ value = (gaussianExponentBgaussianExponentBdxdx) - (gaussianExponentBgaussianExponentBdy*dy);
4705	4723	}
4706	4724	else if(valenceOrbitalA == s && valenceOrbitalB == dxxyy){
4707	4725	value = 2.0*gaussianExponentB;
4708		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
4709		- value = pow(gaussianExponentAdx, 2.0) - pow(gaussianExponentA*dy, 2.0);
	4726	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
	4727	+ value = (gaussianExponentAgaussianExponentAdxdx) - (gaussianExponentAgaussianExponentAdy*dy);
4710	4728	}
4711	4729	else if(valenceOrbitalA == dxxyy && valenceOrbitalB == px){
4712	4730	value = gaussianExponentB*dx;
4713		- value -= pow(gaussianExponentBdx, 2.0)gaussianExponentA*dx/(gaussianExponentA+gaussianExponentB);
4714		- value += pow(gaussianExponentBdy, 2.0)gaussianExponentA*dx/(gaussianExponentA+gaussianExponentB);
4715		- value = -4.0gaussianExponentA*pow(gaussianExponentB, 0.5);
	4731	+ value -= (gaussianExponentBgaussianExponentBdxdx)gaussianExponentA*dx/(gaussianExponentA+gaussianExponentB);
	4732	+ value += (gaussianExponentBgaussianExponentBdydy)gaussianExponentA*dx/(gaussianExponentA+gaussianExponentB);
	4733	+ value = -4.0gaussianExponentA*sqrt(gaussianExponentB);
4716	4734	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4717	4735	}
4718	4736	else if(valenceOrbitalA == px && valenceOrbitalB == dxxyy){
4719	4737	value = gaussianExponentA*dx;
4720		- value -= pow(gaussianExponentAdx, 2.0)gaussianExponentB*dx/(gaussianExponentA+gaussianExponentB);
4721		- value += pow(gaussianExponentAdy, 2.0)gaussianExponentB*dx/(gaussianExponentA+gaussianExponentB);
4722		- value = 4.0gaussianExponentB*pow(gaussianExponentA, 0.5);
	4738	+ value -= (gaussianExponentAgaussianExponentAdxdx)gaussianExponentB*dx/(gaussianExponentA+gaussianExponentB);
	4739	+ value += (gaussianExponentAgaussianExponentAdydy)gaussianExponentB*dx/(gaussianExponentA+gaussianExponentB);
	4740	+ value = 4.0gaussianExponentB*sqrt(gaussianExponentA);
4723	4741	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4724	4742	}
4725	4743	else if(valenceOrbitalA == dxxyy && valenceOrbitalB == py){
4726	4744	value = gaussianExponentB*dy;
4727		- value += pow(gaussianExponentBdx, 2.0)gaussianExponentA*dy/(gaussianExponentA+gaussianExponentB);
4728		- value -= pow(gaussianExponentBdy, 2.0)gaussianExponentA*dy/(gaussianExponentA+gaussianExponentB);
4729		- value = 4.0gaussianExponentA*pow(gaussianExponentB, 0.5);
	4745	+ value += (gaussianExponentBgaussianExponentBdxdx)gaussianExponentA*dy/(gaussianExponentA+gaussianExponentB);
	4746	+ value -= (gaussianExponentBgaussianExponentBdydy)gaussianExponentA*dy/(gaussianExponentA+gaussianExponentB);
	4747	+ value = 4.0gaussianExponentA*sqrt(gaussianExponentB);
4730	4748	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4731	4749	}
4732	4750	else if(valenceOrbitalA == py && valenceOrbitalB == dxxyy){
4733	4751	value = gaussianExponentA*dy;
4734		- value += pow(gaussianExponentAdx, 2.0)gaussianExponentB*dy/(gaussianExponentA+gaussianExponentB);
4735		- value -= pow(gaussianExponentAdy, 2.0)gaussianExponentB*dy/(gaussianExponentA+gaussianExponentB);
4736		- value = -4.0gaussianExponentB*pow(gaussianExponentA, 0.5);
	4752	+ value += (gaussianExponentAgaussianExponentAdxdx)gaussianExponentB*dy/(gaussianExponentA+gaussianExponentB);
	4753	+ value -= (gaussianExponentAgaussianExponentAdydy)gaussianExponentB*dy/(gaussianExponentA+gaussianExponentB);
	4754	+ value = -4.0gaussianExponentB*sqrt(gaussianExponentA);
4737	4755	value /= pow(gaussianExponentA+gaussianExponentB, 2.0);
4738	4756	}
4739	4757	else if(valenceOrbitalA == dxxyy && valenceOrbitalB == pz){
4740		- value = pow(gaussianExponentBdx, 2.0) - pow(gaussianExponentBdy, 2.0);
	4758	+ value = (gaussianExponentBgaussianExponentBdxdx) - (gaussianExponentBgaussianExponentBdydy);
4741	4759	value = gaussianExponentAdz;
4742		- value = 4.0gaussianExponentA*pow(gaussianExponentB, 0.5);
	4760	+ value = 4.0gaussianExponentA*sqrt(gaussianExponentB);
4743	4761	value /= pow(gaussianExponentA+gaussianExponentB, 3.0);
4744	4762	}
4745	4763	else if(valenceOrbitalA == pz && valenceOrbitalB == dxxyy){
4746		- value = pow(gaussianExponentAdx, 2.0) - pow(gaussianExponentAdy, 2.0);
	4764	+ value = (gaussianExponentAgaussianExponentAdxdx) - (gaussianExponentAgaussianExponentAdydy);
4747	4765	value = gaussianExponentBdz;
4748		- value = -4.0gaussianExponentB*pow(gaussianExponentA, 0.5);
	4766	+ value = -4.0gaussianExponentB*sqrt(gaussianExponentA);
4749	4767	value /= pow(gaussianExponentA+gaussianExponentB, 3.0);
4750	4768	}
4751	4769
4752	4770	else if(valenceOrbitalA == dzz && valenceOrbitalB == s){
4753	4771	double temp = 0.0;
4754		- temp = 2.0pow(gaussianExponentBdz, 2.0)
4755		- - pow(gaussianExponentB*dx, 2.0)
4756		- - pow(gaussianExponentB*dy, 2.0);
	4772	+ temp = 2.0(gaussianExponentBgaussianExponentBdzdz)
	4773	+ - (gaussianExponentBgaussianExponentBdx*dx)
	4774	+ - (gaussianExponentBgaussianExponentBdy*dy);
4757	4775	value = 2.0*gaussianExponentA/sqrt(3.0);
4758		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4776	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4759	4777	value *= temp;
4760	4778	}
4761	4779	else if(valenceOrbitalA == s && valenceOrbitalB == dzz){
4762	4780	double temp = 0.0;
4763		- temp = 2.0pow(gaussianExponentAdz, 2.0)
4764		- - pow(gaussianExponentA*dx, 2.0)
4765		- - pow(gaussianExponentA*dy, 2.0);
	4781	+ temp = 2.0(gaussianExponentAgaussianExponentAdzdz)
	4782	+ - (gaussianExponentAgaussianExponentAdx*dx)
	4783	+ - (gaussianExponentAgaussianExponentAdy*dy);
4766	4784	value = 2.0*gaussianExponentB/sqrt(3.0);
4767		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4785	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4768	4786	value *= temp;
4769	4787	}
4770	4788	else if(valenceOrbitalA == dzz && valenceOrbitalB == px){
4771	4789	double temp = 0.0;
4772	4790	temp = gaussianExponentB*dx;
4773		- temp += 2.0pow(gaussianExponentBdz, 2.0)gaussianExponentAdx/(gaussianExponentA+gaussianExponentB);
4774		- temp -= pow(gaussianExponentBdx, 2.0)gaussianExponentA*dx/(gaussianExponentA+gaussianExponentB);
4775		- temp -= pow(gaussianExponentBdy, 2.0)gaussianExponentA*dx/(gaussianExponentA+gaussianExponentB);
	4791	+ temp += 2.0(gaussianExponentBgaussianExponentBdzdz)gaussianExponentAdx/(gaussianExponentA+gaussianExponentB);
	4792	+ temp -= (gaussianExponentBgaussianExponentBdxdx)gaussianExponentA*dx/(gaussianExponentA+gaussianExponentB);
	4793	+ temp -= (gaussianExponentBgaussianExponentBdydy)gaussianExponentA*dx/(gaussianExponentA+gaussianExponentB);
4776	4794	value = temp;
4777		- value = 4.0gaussianExponentA*pow(gaussianExponentB, 0.5)/sqrt(3.0);
4778		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4795	+ value = 4.0gaussianExponentA*sqrt(gaussianExponentB)/sqrt(3.0);
	4796	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4779	4797	}
4780	4798	else if(valenceOrbitalA == px && valenceOrbitalB == dzz){
4781	4799	double temp = 0.0;
4782	4800	temp = gaussianExponentA*dx;
4783		- temp += 2.0pow(gaussianExponentAdz, 2.0)gaussianExponentBdx/(gaussianExponentA+gaussianExponentB);
4784		- temp -= pow(gaussianExponentAdx, 2.0)gaussianExponentB*dx/(gaussianExponentA+gaussianExponentB);
4785		- temp -= pow(gaussianExponentAdy, 2.0)gaussianExponentB*dx/(gaussianExponentA+gaussianExponentB);
	4801	+ temp += 2.0(gaussianExponentAgaussianExponentAdzdz)gaussianExponentBdx/(gaussianExponentA+gaussianExponentB);
	4802	+ temp -= (gaussianExponentAgaussianExponentAdxdx)gaussianExponentB*dx/(gaussianExponentA+gaussianExponentB);
	4803	+ temp -= (gaussianExponentAgaussianExponentAdydy)gaussianExponentB*dx/(gaussianExponentA+gaussianExponentB);
4786	4804	value = temp;
4787		- value = -4.0gaussianExponentB*pow(gaussianExponentA, 0.5)/sqrt(3.0);
4788		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4805	+ value = -4.0gaussianExponentB*sqrt(gaussianExponentA)/sqrt(3.0);
	4806	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4789	4807	}
4790	4808	else if(valenceOrbitalA == dzz && valenceOrbitalB == py){
4791	4809	double temp = 0.0;
4792	4810	temp = gaussianExponentB*dy;
4793		- temp += 2.0pow(gaussianExponentBdz, 2.0)gaussianExponentAdy/(gaussianExponentA+gaussianExponentB);
4794		- temp -= pow(gaussianExponentBdx, 2.0)gaussianExponentA*dy/(gaussianExponentA+gaussianExponentB);
4795		- temp -= pow(gaussianExponentBdy, 2.0)gaussianExponentA*dy/(gaussianExponentA+gaussianExponentB);
	4811	+ temp += 2.0(gaussianExponentBgaussianExponentBdzdz)gaussianExponentAdy/(gaussianExponentA+gaussianExponentB);
	4812	+ temp -= (gaussianExponentBgaussianExponentBdxdx)gaussianExponentA*dy/(gaussianExponentA+gaussianExponentB);
	4813	+ temp -= (gaussianExponentBgaussianExponentBdydy)gaussianExponentA*dy/(gaussianExponentA+gaussianExponentB);
4796	4814	value = temp;
4797		- value = 4.0gaussianExponentA*pow(gaussianExponentB, 0.5)/sqrt(3.0);
4798		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4815	+ value = 4.0gaussianExponentA*sqrt(gaussianExponentB)/sqrt(3.0);
	4816	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4799	4817	}
4800	4818	else if(valenceOrbitalA == py && valenceOrbitalB == dzz){
4801	4819	double temp = 0.0;
4802	4820	temp = gaussianExponentA*dy;
4803		- temp += 2.0pow(gaussianExponentAdz, 2.0)gaussianExponentBdy/(gaussianExponentA+gaussianExponentB);
4804		- temp -= pow(gaussianExponentAdx, 2.0)gaussianExponentB*dy/(gaussianExponentA+gaussianExponentB);
4805		- temp -= pow(gaussianExponentAdy, 2.0)gaussianExponentB*dy/(gaussianExponentA+gaussianExponentB);
	4821	+ temp += 2.0(gaussianExponentAgaussianExponentAdzdz)gaussianExponentBdy/(gaussianExponentA+gaussianExponentB);
	4822	+ temp -= (gaussianExponentAgaussianExponentAdxdx)gaussianExponentB*dy/(gaussianExponentA+gaussianExponentB);
	4823	+ temp -= (gaussianExponentAgaussianExponentAdydy)gaussianExponentB*dy/(gaussianExponentA+gaussianExponentB);
4806	4824	value = temp;
4807		- value = -4.0gaussianExponentB*pow(gaussianExponentA, 0.5)/sqrt(3.0);
4808		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4825	+ value = -4.0gaussianExponentB*sqrt(gaussianExponentA)/sqrt(3.0);
	4826	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4809	4827	}
4810	4828	else if(valenceOrbitalA == dzz && valenceOrbitalB == pz){
4811	4829	double temp = 0.0;
4812	4830	temp = -2.0gaussianExponentBdz;
4813		- temp += 2.0pow(gaussianExponentBdz, 2.0)gaussianExponentAdz/(gaussianExponentA+gaussianExponentB);
4814		- temp -= pow(gaussianExponentBdx, 2.0)gaussianExponentA*dz/(gaussianExponentA+gaussianExponentB);
4815		- temp -= pow(gaussianExponentBdy, 2.0)gaussianExponentA*dz/(gaussianExponentA+gaussianExponentB);
	4831	+ temp += 2.0(gaussianExponentBgaussianExponentBdzdz)gaussianExponentAdz/(gaussianExponentA+gaussianExponentB);
	4832	+ temp -= (gaussianExponentBgaussianExponentBdxdx)gaussianExponentA*dz/(gaussianExponentA+gaussianExponentB);
	4833	+ temp -= (gaussianExponentBgaussianExponentBdydy)gaussianExponentA*dz/(gaussianExponentA+gaussianExponentB);
4816	4834	value = temp;
4817		- value = 4.0gaussianExponentA*pow(gaussianExponentB, 0.5)/sqrt(3.0);
4818		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4835	+ value = 4.0gaussianExponentA*sqrt(gaussianExponentB)/sqrt(3.0);
	4836	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4819	4837	}
4820	4838	else if(valenceOrbitalA == pz && valenceOrbitalB == dzz){
4821	4839	double temp = 0.0;
4822	4840	temp = -2.0gaussianExponentAdz;
4823		- temp += 2.0pow(gaussianExponentAdz, 2.0)gaussianExponentBdz/(gaussianExponentA+gaussianExponentB);
4824		- temp -= pow(gaussianExponentAdx, 2.0)gaussianExponentB*dz/(gaussianExponentA+gaussianExponentB);
4825		- temp -= pow(gaussianExponentAdy, 2.0)gaussianExponentB*dz/(gaussianExponentA+gaussianExponentB);
	4841	+ temp += 2.0(gaussianExponentAgaussianExponentAdzdz)gaussianExponentBdz/(gaussianExponentA+gaussianExponentB);
	4842	+ temp -= (gaussianExponentAgaussianExponentAdxdx)gaussianExponentB*dz/(gaussianExponentA+gaussianExponentB);
	4843	+ temp -= (gaussianExponentAgaussianExponentAdydy)gaussianExponentB*dz/(gaussianExponentA+gaussianExponentB);
4826	4844	value = temp;
4827		- value = -4.0gaussianExponentB*pow(gaussianExponentA, 0.5)/sqrt(3.0);
4828		- value *= pow(gaussianExponentA+gaussianExponentB, -2.0);
	4845	+ value = -4.0gaussianExponentB*sqrt(gaussianExponentA)/sqrt(3.0);
	4846	+ value /= (gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB);
4829	4847	}
4830	4848
4831	4849	else if(valenceOrbitalA == dxy && valenceOrbitalB == dxy){

		@@ -4834,9 +4852,9 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4834	4852	temp -= 0.5gaussianExponentBdygaussianExponentAdy/(gaussianExponentA+gaussianExponentB);
4835	4853	temp += gaussianExponentBdxgaussianExponentA*dx
4836	4854	gaussianExponentBdygaussianExponentAdy
4837		- *pow(gaussianExponentA+gaussianExponentB, -2.0);
	4855	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
4838	4856	value = 16.0tempgaussianExponentA*gaussianExponentB
4839		- *pow(gaussianExponentA+gaussianExponentB, -2.0);
	4857	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
4840	4858	}
4841	4859	else if(valenceOrbitalA == dyz && valenceOrbitalB == dyz){
4842	4860	double temp = 0.25;

		@@ -4844,9 +4862,9 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4844	4862	temp -= 0.5gaussianExponentBdzgaussianExponentAdz/(gaussianExponentA+gaussianExponentB);
4845	4863	temp += gaussianExponentBdygaussianExponentA*dy
4846	4864	gaussianExponentBdzgaussianExponentAdz
4847		- *pow(gaussianExponentA+gaussianExponentB, -2.0);
	4865	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
4848	4866	value = 16.0tempgaussianExponentA*gaussianExponentB
4849		- *pow(gaussianExponentA+gaussianExponentB, -2.0);
	4867	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
4850	4868	}
4851	4869	else if(valenceOrbitalA == dzx && valenceOrbitalB == dzx){
4852	4870	double temp = 0.25;

		@@ -4854,42 +4872,41 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4854	4872	temp -= 0.5gaussianExponentBdxgaussianExponentAdx/(gaussianExponentA+gaussianExponentB);
4855	4873	temp += gaussianExponentBdzgaussianExponentA*dz
4856	4874	gaussianExponentBdxgaussianExponentAdx
4857		- *pow(gaussianExponentA+gaussianExponentB, -2.0);
	4875	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
4858	4876	value = 16.0tempgaussianExponentA*gaussianExponentB
4859		- *pow(gaussianExponentA+gaussianExponentB, -2.0);
	4877	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
4860	4878	}
4861	4879	else if(valenceOrbitalA == dxxyy && valenceOrbitalB == dxxyy){
4862	4880	double temp1 = 1.0;
4863	4881	temp1 -= 2.0gaussianExponentBdxgaussianExponentAdx/(gaussianExponentA+gaussianExponentB);
4864	4882	temp1 -= 2.0gaussianExponentBdygaussianExponentAdy/(gaussianExponentA+gaussianExponentB);
4865		- double temp2 = gaussianExponentAgaussianExponentB(pow(dx,2.0)-pow(dy,2.0))
	4883	+ double temp2 = gaussianExponentAgaussianExponentB((dxdx)-(dydy))
4866	4884	/(gaussianExponentA+gaussianExponentB);
4867	4885	temp1 += pow(temp2,2.0);
4868	4886	value = 4.0temp1gaussianExponentA*gaussianExponentB
4869		- *pow(gaussianExponentA+gaussianExponentB, -2.0);
	4887	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
4870	4888	}
4871	4889	else if(valenceOrbitalA == dzz && valenceOrbitalB == dzz){
4872	4890	double temp = 3.0;
4873	4891	temp -= gaussianExponentA*gaussianExponentB
4874	4892	/(gaussianExponentA+gaussianExponentB)
4875		- (8.0pow(dz,2.0)+2.0pow(dx,2.0)+2.0pow(dy,2.0));
	4893	+ (8.0(dzdz)+2.0(dxdx)+2.0(dy*dy));
4876	4894	temp += pow(gaussianExponentA*gaussianExponentB,2.0)
4877		- *pow(gaussianExponentA+gaussianExponentB, -2.0)
4878		- (4.0pow(dz,4.0)
4879		- +pow(dx,4.0)
4880		- +pow(dy,4.0)
4881		- -4.0pow(dxdz,2.0)
4882		- -4.0pow(dydz,2.0)
4883		- +2.0pow(dxdy,2.0));
	4895	+ (4.0(dzdzdz*dz)
	4896	+ +(dxdxdx*dx)
	4897	+ +(dydydy*dy)
	4898	+ -4.0(dxdxdzdz)
	4899	+ -4.0(dydydzdz)
	4900	+ +2.0(dxdxdydy))
	4901	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
4884	4902	value = 4.0tempgaussianExponentA*gaussianExponentB
4885		- *pow(gaussianExponentA+gaussianExponentB, -2.0)
4886		- /3.0;
	4903	+ /(3.0(gaussianExponentA+gaussianExponentB)(gaussianExponentA+gaussianExponentB));
4887	4904	}
4888	4905
4889	4906	else if((valenceOrbitalA == dxy && valenceOrbitalB == dyz) \|\|
4890	4907	(valenceOrbitalA == dyz && valenceOrbitalB == dxy)){
4891	4908	double temp = 0.5;
4892		- temp -= gaussianExponentAgaussianExponentBpow(dy,2.0)/(gaussianExponentA+gaussianExponentB);
	4909	+ temp -= gaussianExponentAgaussianExponentB(dy*dy)/(gaussianExponentA+gaussianExponentB);
4893	4910	value = -16.0pow(gaussianExponentAgaussianExponentB,2.0)
4894	4911	*pow(gaussianExponentA+gaussianExponentB,-3.0)
4895	4912	dxdz*temp;

		@@ -4897,7 +4914,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4897	4914	else if((valenceOrbitalA == dyz && valenceOrbitalB == dzx) \|\|
4898	4915	(valenceOrbitalA == dzx && valenceOrbitalB == dyz)){
4899	4916	double temp = 0.5;
4900		- temp -= gaussianExponentAgaussianExponentBpow(dz,2.0)/(gaussianExponentA+gaussianExponentB);
	4917	+ temp -= gaussianExponentAgaussianExponentB(dz*dz)/(gaussianExponentA+gaussianExponentB);
4901	4918	value = -16.0pow(gaussianExponentAgaussianExponentB,2.0)
4902	4919	*pow(gaussianExponentA+gaussianExponentB,-3.0)
4903	4920	dydx*temp;

		@@ -4905,7 +4922,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4905	4922	else if((valenceOrbitalA == dzx && valenceOrbitalB == dxy) \|\|
4906	4923	(valenceOrbitalA == dxy && valenceOrbitalB == dzx)){
4907	4924	double temp = 0.5;
4908		- temp -= gaussianExponentAgaussianExponentBpow(dx,2.0)/(gaussianExponentA+gaussianExponentB);
	4925	+ temp -= gaussianExponentAgaussianExponentB(dx*dx)/(gaussianExponentA+gaussianExponentB);
4909	4926	value = -16.0pow(gaussianExponentAgaussianExponentB,2.0)
4910	4927	*pow(gaussianExponentA+gaussianExponentB,-3.0)
4911	4928	dzdy*temp;

		@@ -4922,7 +4939,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4922	4939	double temp = 2.0gaussianExponentAgaussianExponentB;
4923	4940	value = pow(temp,3.0)(dydz*(gaussianExponentA+gaussianExponentB)
4924	4941	/(gaussianExponentA*gaussianExponentB)
4925		- +(pow(dx,2.0)dydz - pow(dy,3.0)*dz))
	4942	+ +((dxdx)dydz - pow(dy,3.0)dz))
4926	4943	*pow(gaussianExponentA+gaussianExponentB,-4.0);
4927	4944	}
4928	4945	else if((valenceOrbitalA == dxxyy && valenceOrbitalB == dzx) \|\|

		@@ -4930,13 +4947,13 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4930	4947	double temp = 2.0gaussianExponentAgaussianExponentB;
4931	4948	value = -1.0pow(temp,3.0)(dxdz(gaussianExponentA+gaussianExponentB)
4932	4949	/(gaussianExponentA*gaussianExponentB)
4933		- +(pow(dy,2.0)dxdz - pow(dx,3.0)*dz))
	4950	+ +((dydy)dxdz - pow(dx,3.0)dz))
4934	4951	*pow(gaussianExponentA+gaussianExponentB,-4.0);
4935	4952	}
4936	4953
4937	4954	else if((valenceOrbitalA == dzz && valenceOrbitalB == dxy) \|\|
4938	4955	(valenceOrbitalA == dxy && valenceOrbitalB == dzz)){
4939		- double temp = 2.0dxdypow(dz,2.0) - pow(dx,3.0)dy - dx*pow(dy,3.0);
	4956	+ double temp = 2.0dxdy(dzdz) - pow(dx,3.0)dy - dxpow(dy,3.0);
4940	4957	temp = gaussianExponentAgaussianExponentB/(gaussianExponentA+gaussianExponentB);
4941	4958	temp += 2.0dxdy;
4942	4959	value = 8.0pow(gaussianExponentAgaussianExponentB,2.0)*temp;

		@@ -4945,7 +4962,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4945	4962	else if((valenceOrbitalA == dzz && valenceOrbitalB == dyz) \|\|
4946	4963	(valenceOrbitalA == dyz && valenceOrbitalB == dzz)){
4947	4964	double temp1 = -1.0dydz;
4948		- double temp2 = 2.0dypow(dz,3.0) - pow(dy,3.0)dz - pow(dx,2.0)dy*dz;
	4965	+ double temp2 = 2.0dypow(dz,3.0) - pow(dy,3.0)dz - (dxdx)dydz;
4949	4966	temp2 = gaussianExponentAgaussianExponentB/(gaussianExponentA+gaussianExponentB);
4950	4967	temp1 += temp2;
4951	4968	value = 8.0pow(gaussianExponentAgaussianExponentB,2.0)*temp1;

		@@ -4954,7 +4971,7 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4954	4971	else if((valenceOrbitalA == dzz && valenceOrbitalB == dzx) \|\|
4955	4972	(valenceOrbitalA == dzx && valenceOrbitalB == dzz)){
4956	4973	double temp1 = -1.0dxdz;
4957		- double temp2 = 2.0dxpow(dz,3.0) - pow(dx,3.0)dz - pow(dy,2.0)dx*dz;
	4974	+ double temp2 = 2.0dxpow(dz,3.0) - pow(dx,3.0)dz - (dydy)dxdz;
4958	4975	temp2 = gaussianExponentAgaussianExponentB/(gaussianExponentA+gaussianExponentB);
4959	4976	temp1 += temp2;
4960	4977	value = 8.0pow(gaussianExponentAgaussianExponentB,2.0)*temp1;

		@@ -4962,12 +4979,12 @@ double Cndo2::GetGaussianOverlapAOs(AtomType atomTypeA,
4962	4979	}
4963	4980	else if((valenceOrbitalA == dxxyy && valenceOrbitalB == dzz) \|\|
4964	4981	(valenceOrbitalA == dzz && valenceOrbitalB == dxxyy)){
4965		- double temp = 2.0*pow(dz,2.0)-pow(dx,2.0)-pow(dy,2.0);
	4982	+ double temp = 2.0(dzdz)-(dxdx)-(dydy);
4966	4983	temp = gaussianExponentAgaussianExponentB/(gaussianExponentA+gaussianExponentB);
4967	4984	temp += 2.0;
4968	4985	value = 4.0pow(gaussianExponentAgaussianExponentB,2.0);
4969	4986	value /= sqrt(3.0)*pow(gaussianExponentA+gaussianExponentB,3.0);
4970		- value = (pow(dx,2.0)-pow(dy,2.0))temp;
	4987	+ value = ((dxdx)-(dydy))temp;
4971	4988	}
4972	4989
4973	4990	else{

		@@ -5000,7 +5017,7 @@ double Cndo2::GetOverlapAOsElement1stDerivativeByGTOExpansion(const Atom& atomA,
5000	5017	double dx = atomA.GetXyz()[XAxis] - atomB.GetXyz()[XAxis];
5001	5018	double dy = atomA.GetXyz()[YAxis] - atomB.GetXyz()[YAxis];
5002	5019	double dz = atomA.GetXyz()[ZAxis] - atomB.GetXyz()[ZAxis];
5003		- double rAB = sqrt( pow(dx, 2.0) + pow(dy, 2.0) + pow(dz,2.0) );
	5020	+ double rAB = sqrt( dxdx + dydy + dz*dz );
5004	5021	ShellType shellTypeA = atomA.GetValenceShellType();
5005	5022	ShellType shellTypeB = atomB.GetValenceShellType();
5006	5023	OrbitalType valenceOrbitalA = atomA.GetValence(valenceIndexA);

		@@ -5025,12 +5042,12 @@ double Cndo2::GetOverlapAOsElement1stDerivativeByGTOExpansion(const Atom& atomA,
5025	5042	shellTypeB,
5026	5043	valenceOrbitalB,
5027	5044	j);
5028		- gaussianExponentA = pow(orbitalExponentA, 2.0)
	5045	+ gaussianExponentA = (orbitalExponentA*orbitalExponentA)
5029	5046	*GTOExpansionSTO::GetInstance()->GetExponent(stonG,
5030	5047	shellTypeA,
5031	5048	valenceOrbitalA,
5032	5049	i);
5033		- gaussianExponentB = pow(orbitalExponentB, 2.0)
	5050	+ gaussianExponentB = (orbitalExponentB*orbitalExponentB)
5034	5051	*GTOExpansionSTO::GetInstance()->GetExponent(stonG,
5035	5052	shellTypeB,
5036	5053	valenceOrbitalB,

		@@ -5079,10 +5096,10 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5079	5096	}
5080	5097	}
5081	5098	else if(valenceOrbitalA == s && valenceOrbitalB == px){
5082		- double temp1 = 4.0pow(gaussianExponentA,2.0)pow(gaussianExponentB, 1.5)
5083		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5099	+ double temp1 = 4.0(gaussianExponentAgaussianExponentA)gaussianExponentBsqrt(gaussianExponentB)
	5100	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5084	5101	if(axisA == XAxis){
5085		- double temp2 = 2.0gaussianExponentApow(gaussianExponentB, 0.5)
	5102	+ double temp2 = 2.0gaussianExponentAsqrt(gaussianExponentB)
5086	5103	/(gaussianExponentA+gaussianExponentB);
5087	5104	value = temp2-temp1dxdx;
5088	5105	}

		@@ -5094,13 +5111,13 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5094	5111	}
5095	5112	}
5096	5113	else if(valenceOrbitalA == s && valenceOrbitalB == py){
5097		- double temp1 = 4.0pow(gaussianExponentA,2.0)pow(gaussianExponentB, 1.5)
5098		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5114	+ double temp1 = 4.0(gaussianExponentAgaussianExponentA)gaussianExponentBsqrt(gaussianExponentB)
	5115	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5099	5116	if(axisA == XAxis){
5100	5117	value = -1.0temp1dx*dy;
5101	5118	}
5102	5119	else if(axisA == YAxis){
5103		- double temp2 = 2.0gaussianExponentApow(gaussianExponentB, 0.5)
	5120	+ double temp2 = 2.0gaussianExponentAsqrt(gaussianExponentB)
5104	5121	/(gaussianExponentA+gaussianExponentB);
5105	5122	value = temp2-temp1dydy;
5106	5123	}

		@@ -5109,8 +5126,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5109	5126	}
5110	5127	}
5111	5128	else if(valenceOrbitalA == s && valenceOrbitalB == pz){
5112		- double temp1 = 4.0pow(gaussianExponentA,2.0)pow(gaussianExponentB, 1.5)
5113		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5129	+ double temp1 = 4.0(gaussianExponentAgaussianExponentA)gaussianExponentBsqrt(gaussianExponentB)
	5130	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5114	5131	if(axisA == XAxis){
5115	5132	value = -1.0temp1dx*dz;
5116	5133	}

		@@ -5118,16 +5135,17 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5118	5135	value = -1.0temp1dy*dz;
5119	5136	}
5120	5137	else if(axisA == ZAxis){
5121		- double temp2 = 2.0gaussianExponentApow(gaussianExponentB, 0.5)
	5138	+ double temp2 = 2.0gaussianExponentAsqrt(gaussianExponentB)
5122	5139	/(gaussianExponentA+gaussianExponentB);
5123	5140	value = temp2-temp1dzdz;
5124	5141	}
5125	5142	}
5126	5143	else if(valenceOrbitalA == px && valenceOrbitalB == s){
5127		- double temp1 = 4.0pow(gaussianExponentA,1.5)pow(gaussianExponentB, 2.0)
5128		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5144	+ double temp1 = 4.0gaussianExponentAsqrt(gaussianExponentA)
	5145	+ gaussianExponentBgaussianExponentB
	5146	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5129	5147	if(axisA == XAxis){
5130		- double temp2 = 2.0pow(gaussianExponentA,0.5)gaussianExponentB
	5148	+ double temp2 = 2.0sqrt(gaussianExponentA)gaussianExponentB
5131	5149	/(gaussianExponentA+gaussianExponentB);
5132	5150	value = -1.0temp2+temp1dx*dx;
5133	5151	}

		@@ -5139,13 +5157,14 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5139	5157	}
5140	5158	}
5141	5159	else if(valenceOrbitalA == py && valenceOrbitalB == s){
5142		- double temp1 = 4.0pow(gaussianExponentA,1.5)pow(gaussianExponentB, 2.0)
5143		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5160	+ double temp1 = 4.0gaussianExponentAsqrt(gaussianExponentA)
	5161	+ gaussianExponentBgaussianExponentB
	5162	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5144	5163	if(axisA == XAxis){
5145	5164	value = temp1dxdy;
5146	5165	}
5147	5166	else if(axisA == YAxis){
5148		- double temp2 = 2.0pow(gaussianExponentA,0.5)gaussianExponentB
	5167	+ double temp2 = 2.0sqrt(gaussianExponentA)gaussianExponentB
5149	5168	/(gaussianExponentA+gaussianExponentB);
5150	5169	value = -1.0temp2+temp1dy*dy;
5151	5170	}

		@@ -5154,8 +5173,9 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5154	5173	}
5155	5174	}
5156	5175	else if(valenceOrbitalA == pz && valenceOrbitalB == s){
5157		- double temp1 = 4.0pow(gaussianExponentA,1.5)pow(gaussianExponentB, 2.0)
5158		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5176	+ double temp1 = 4.0gaussianExponentAsqrt(gaussianExponentA)
	5177	+ gaussianExponentBgaussianExponentB
	5178	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5159	5179	if(axisA == XAxis){
5160	5180	value = temp1dxdz;
5161	5181	}

		@@ -5163,7 +5183,7 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5163	5183	value = temp1dydz;
5164	5184	}
5165	5185	else if(axisA == ZAxis){
5166		- double temp2 = 2.0pow(gaussianExponentA,0.5)gaussianExponentB
	5186	+ double temp2 = 2.0sqrt(gaussianExponentA)gaussianExponentB
5167	5187	/(gaussianExponentA+gaussianExponentB);
5168	5188	value = -1.0temp2+temp1dz*dz;
5169	5189	}

		@@ -5171,8 +5191,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5171	5191	else if(valenceOrbitalA == px && valenceOrbitalB == py){
5172	5192	double temp1 = 8.0pow(gaussianExponentAgaussianExponentB,2.5)
5173	5193	/pow(gaussianExponentA+gaussianExponentB,3.0);
5174		- double temp2 = 4.0pow(gaussianExponentAgaussianExponentB,1.5)
5175		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5194	+ double temp2 = 4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB)
	5195	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5176	5196	if(axisA == XAxis){
5177	5197	value = -1.0temp2dy+temp1dxdx*dy;
5178	5198	}

		@@ -5186,8 +5206,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5186	5206	else if(valenceOrbitalA == py && valenceOrbitalB == px){
5187	5207	double temp1 = 8.0pow(gaussianExponentAgaussianExponentB,2.5)
5188	5208	/pow(gaussianExponentA+gaussianExponentB,3.0);
5189		- double temp2 = 4.0pow(gaussianExponentAgaussianExponentB,1.5)
5190		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5209	+ double temp2 = 4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB)
	5210	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5191	5211	if(axisA == XAxis){
5192	5212	value = -1.0temp2dy+temp1dydx*dx;
5193	5213	}

		@@ -5201,8 +5221,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5201	5221	else if(valenceOrbitalA == px && valenceOrbitalB == pz){
5202	5222	double temp1 = 8.0pow(gaussianExponentAgaussianExponentB,2.5)
5203	5223	/pow(gaussianExponentA+gaussianExponentB,3.0);
5204		- double temp2 = 4.0pow(gaussianExponentAgaussianExponentB,1.5)
5205		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5224	+ double temp2 = 4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB)
	5225	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5206	5226	if(axisA == XAxis){
5207	5227	value = -1.0temp2dz+temp1dxdx*dz;
5208	5228	}

		@@ -5216,8 +5236,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5216	5236	else if(valenceOrbitalA == pz && valenceOrbitalB == px){
5217	5237	double temp1 = 8.0pow(gaussianExponentAgaussianExponentB,2.5)
5218	5238	/pow(gaussianExponentA+gaussianExponentB,3.0);
5219		- double temp2 = 4.0pow(gaussianExponentAgaussianExponentB,1.5)
5220		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5239	+ double temp2 = 4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB)
	5240	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5221	5241	if(axisA == XAxis){
5222	5242	value = -1.0temp2dz+temp1dzdx*dx;
5223	5243	}

		@@ -5231,8 +5251,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5231	5251	else if(valenceOrbitalA == py && valenceOrbitalB == pz){
5232	5252	double temp1 = 8.0pow(gaussianExponentAgaussianExponentB,2.5)
5233	5253	/pow(gaussianExponentA+gaussianExponentB,3.0);
5234		- double temp2 = 4.0pow(gaussianExponentAgaussianExponentB,1.5)
5235		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5254	+ double temp2 = 4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB)
	5255	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5236	5256	if(axisA == XAxis){
5237	5257	value = temp1dxdy*dz;
5238	5258	}

		@@ -5246,8 +5266,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5246	5266	else if(valenceOrbitalA == pz && valenceOrbitalB == py){
5247	5267	double temp1 = 8.0pow(gaussianExponentAgaussianExponentB,2.5)
5248	5268	/pow(gaussianExponentA+gaussianExponentB,3.0);
5249		- double temp2 = 4.0pow(gaussianExponentAgaussianExponentB,1.5)
5250		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5269	+ double temp2 = 4.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB)
	5270	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5251	5271	if(axisA == XAxis){
5252	5272	value = temp1dxdy*dz;
5253	5273	}

		@@ -5259,8 +5279,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5259	5279	}
5260	5280	}
5261	5281	else if(valenceOrbitalA == px && valenceOrbitalB == px){
5262		- double temp1 = 8.0pow(gaussianExponentAgaussianExponentB,1.5)
5263		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5282	+ double temp1 = 8.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB)
	5283	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5264	5284	double temp2 = gaussianExponentA*gaussianExponentB
5265	5285	/(gaussianExponentA+gaussianExponentB);
5266	5286	if(axisA == XAxis){

		@@ -5274,8 +5294,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5274	5294	}
5275	5295	}
5276	5296	else if(valenceOrbitalA == py && valenceOrbitalB == py){
5277		- double temp1 = 8.0pow(gaussianExponentAgaussianExponentB,1.5)
5278		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5297	+ double temp1 = 8.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB)
	5298	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5279	5299	double temp2 = gaussianExponentA*gaussianExponentB
5280	5300	/(gaussianExponentA+gaussianExponentB);
5281	5301	if(axisA == XAxis){

		@@ -5289,8 +5309,8 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
5289	5309	}
5290	5310	}
5291	5311	else if(valenceOrbitalA == pz && valenceOrbitalB == pz){
5292		- double temp1 = 8.0pow(gaussianExponentAgaussianExponentB,1.5)
5293		- /pow(gaussianExponentA+gaussianExponentB,2.0);
	5312	+ double temp1 = 8.0gaussianExponentAgaussianExponentBsqrt(gaussianExponentAgaussianExponentB)
	5313	+ /((gaussianExponentA+gaussianExponentB)*(gaussianExponentA+gaussianExponentB));
5294	5314	double temp2 = gaussianExponentA*gaussianExponentB
5295	5315	/(gaussianExponentA+gaussianExponentB);
5296	5316	if(axisA == XAxis){

		@@ -5361,15 +5381,15 @@ void Cndo2::CalcRotatingMatrix(double** rotatingMatrix,
5361	5381	// rotating matrix for d-function
5362	5382	// dMatrix is (37) in J. Mol. Strct. 419, 19(1997) (ref. [BFB_1997])
5363	5383	double dMatrix[OrbitalType_end][OrbitalType_end];
5364		- dMatrix[dzz][dzz] = 0.5(3.0pow(cos(beta),2.0) - 1.0);
5365		- dMatrix[dxxyy][dxxyy] = pow(cos(0.5*beta),4.0);
5366		- dMatrix[dzx][dzx] = (2.0cos(beta)-1.0)pow(cos(0.5*beta),2.0);
5367		- dMatrix[dxxyy][dzx] = -2.0sin(0.5beta)pow(cos(0.5beta),3.0);
5368		- dMatrix[dxxyy][dzz] = sqrt(6.0)pow(sin(0.5beta),2.0)pow(cos(0.5beta),2.0);
5369		- dMatrix[dxxyy][dyz] = -2.0pow(sin(0.5beta),3.0)pow(cos(0.5beta),1.0);
5370		- dMatrix[dxxyy][dxy] = pow(sin(0.5*beta),4.0);
	5384	+ dMatrix[dzz][dzz] = 0.5(3.0(cos(beta)*cos(beta)) - 1.0);
	5385	+ dMatrix[dxxyy][dxxyy] = cos(0.5beta)cos(0.5beta)cos(0.5beta)cos(0.5*beta);
	5386	+ dMatrix[dzx][dzx] = (2.0cos(beta)-1.0)cos(0.5beta)cos(0.5*beta);
	5387	+ dMatrix[dxxyy][dzx] = -2.0sin(0.5beta)cos(0.5beta)cos(0.5beta)cos(0.5beta);
	5388	+ dMatrix[dxxyy][dzz] = sqrt(6.0)(sin(0.5beta)sin(0.5beta))((cos(0.5beta))(cos(0.5beta)));
	5389	+ dMatrix[dxxyy][dyz] = -2.0sin(0.5beta)sin(0.5beta)sin(0.5beta)cos(0.5beta);
	5390	+ dMatrix[dxxyy][dxy] = sin(0.5beta)sin(0.5beta)sin(0.5beta)sin(0.5*beta);
5371	5391	dMatrix[dzx][dzz] = -sqrt(6.0)cos(beta)cos(0.5beta)sin(0.5*beta);
5372		- dMatrix[dzx][dyz] = (2.0cos(beta)+1.0)pow(sin(0.5*beta),2.0);
	5392	+ dMatrix[dzx][dyz] = (2.0cos(beta)+1.0)(sin(0.5beta)sin(0.5*beta));
5373	5393
5374	5394	rotatingMatrix[dxy][dxy] = cos(2.0alpha) (dMatrix[dxxyy][dxxyy] - dMatrix[dxxyy][dxy]);
5375	5395	rotatingMatrix[dxy][dyz] = cos(2.0alpha) (-1.0*dMatrix[dxxyy][dzx] - dMatrix[dxxyy][dyz]);

		@@ -5420,8 +5440,8 @@ void Cndo2::CalcRotatingMatrix1stDerivatives(double*** rotMat1stDerivatives,
5420	5440	double x = atomB.GetXyz()[0] - atomA.GetXyz()[0];
5421	5441	double y = atomB.GetXyz()[1] - atomA.GetXyz()[1];
5422	5442	double z = atomB.GetXyz()[2] - atomA.GetXyz()[2];
5423		- double r = sqrt( pow(x,2.0) + pow(y,2.0) );
5424		- double R = sqrt( pow(x,2.0) + pow(y,2.0) + pow(z,2.0) );
	5443	+ double r = sqrt( xx + yy );
	5444	+ double R = sqrt( xx + yy + z*z );
5425	5445
5426	5446	if(r==0e0){
5427	5447	return;

		@@ -5433,43 +5453,43 @@ void Cndo2::CalcRotatingMatrix1stDerivatives(double*** rotMat1stDerivatives,
5433	5453	rotMat1stDerivatives[s][s][ZAxis] = 0.0;
5434	5454
5435	5455	// for p-function
5436		- rotMat1stDerivatives[py][py][XAxis] = -1.0/r + pow(x,2.0)/pow(r,3.0);
5437		- rotMat1stDerivatives[py][pz][XAxis] = x*y/pow(R,3.0);
5438		- rotMat1stDerivatives[py][px][XAxis] = (1.0/(pow(r,3.0)R) + 1.0/(pow(R,3.0)r))xy*z;
	5456	+ rotMat1stDerivatives[py][py][XAxis] = -1.0/r + xx/(rr*r);
	5457	+ rotMat1stDerivatives[py][pz][XAxis] = xy/(RR*R);
	5458	+ rotMat1stDerivatives[py][px][XAxis] = (1.0/(rrrR) + 1.0/(RRRr))xy*z;
5439	5459
5440	5460	rotMat1stDerivatives[pz][py][XAxis] = 0.0;
5441		- rotMat1stDerivatives[pz][pz][XAxis] = x*z/pow(R,3.0);
5442		- rotMat1stDerivatives[pz][px][XAxis] = x/(rR) - xr/pow(R,3.0);
	5461	+ rotMat1stDerivatives[pz][pz][XAxis] = xz/(RR*R);
	5462	+ rotMat1stDerivatives[pz][px][XAxis] = x/(rR) - xr/(RRR);
5443	5463
5444		- rotMat1stDerivatives[px][py][XAxis] = -1.0xy/pow(r,3.0);
5445		- rotMat1stDerivatives[px][pz][XAxis] = -1.0/R + x*x/pow(R,3.0);
	5464	+ rotMat1stDerivatives[px][py][XAxis] = -1.0xy/(rrr);
	5465	+ rotMat1stDerivatives[px][pz][XAxis] = -1.0/R + xx/(RR*R);
5446	5466	rotMat1stDerivatives[px][px][XAxis] = -1.0z/(rR) +
5447		- (1.0/(pow(r,3.0)R) + 1.0/(pow(R,3.0)r))xx*z;
	5467	+ (1.0/(rrrR) + 1.0/(RRRr))xx*z;
5448	5468
5449		- rotMat1stDerivatives[py][py][YAxis] = x*y/pow(r,3.0);
5450		- rotMat1stDerivatives[py][pz][YAxis] = -1.0/R + y*y/pow(R,3.0);
	5469	+ rotMat1stDerivatives[py][py][YAxis] = xy/(rr*r);
	5470	+ rotMat1stDerivatives[py][pz][YAxis] = -1.0/R + yy/(RR*R);
5451	5471	rotMat1stDerivatives[py][px][YAxis] = -1.0z/(rR) +
5452		- (1.0/(pow(r,3.0)R) + 1.0/(pow(R,3.0)r))yy*z;
	5472	+ (1.0/(rrrR) + 1.0/(RRRr))yy*z;
5453	5473
5454	5474	rotMat1stDerivatives[pz][py][YAxis] = 0.0;
5455		- rotMat1stDerivatives[pz][pz][YAxis] = y*z/pow(R,3.0);
5456		- rotMat1stDerivatives[pz][px][YAxis] = y/(rR) - yr/pow(R,3.0);
	5475	+ rotMat1stDerivatives[pz][pz][YAxis] = yz/(RR*R);
	5476	+ rotMat1stDerivatives[pz][px][YAxis] = y/(rR) - yr/(RRR);
5457	5477
5458		- rotMat1stDerivatives[px][py][YAxis] = 1.0/r - y*y/pow(r,3.0);
5459		- rotMat1stDerivatives[px][pz][YAxis] = x*y/pow(R,3.0);
5460		- rotMat1stDerivatives[px][px][YAxis] = (1.0/(pow(r,3.0)R) + 1.0/(pow(R,3.0)r))xy*z;
	5478	+ rotMat1stDerivatives[px][py][YAxis] = 1.0/r - yy/(rr*r);
	5479	+ rotMat1stDerivatives[px][pz][YAxis] = xy/(RR*R);
	5480	+ rotMat1stDerivatives[px][px][YAxis] = (1.0/(rrrR) + 1.0/(RRRr))xy*z;
5461	5481
5462	5482	rotMat1stDerivatives[py][py][ZAxis] = 0.0;
5463		- rotMat1stDerivatives[py][pz][ZAxis] = y*z/pow(R,3.0);
5464		- rotMat1stDerivatives[py][px][ZAxis] = -1.0y/(rR) + yzz/(r*pow(R,3.0));
	5483	+ rotMat1stDerivatives[py][pz][ZAxis] = yz/(RR*R);
	5484	+ rotMat1stDerivatives[py][px][ZAxis] = -1.0y/(rR) + yzz/(rRR*R);
5465	5485
5466	5486	rotMat1stDerivatives[pz][py][ZAxis] = 0.0;
5467		- rotMat1stDerivatives[pz][pz][ZAxis] = -1.0/R + z*z/pow(R,3.0);
5468		- rotMat1stDerivatives[pz][px][ZAxis] = -1.0zr/pow(R,3.0);
	5487	+ rotMat1stDerivatives[pz][pz][ZAxis] = -1.0/R + zz/(RR*R);
	5488	+ rotMat1stDerivatives[pz][px][ZAxis] = -1.0zr/(RRR);
5469	5489
5470	5490	rotMat1stDerivatives[px][py][ZAxis] = 0.0;
5471		- rotMat1stDerivatives[px][pz][ZAxis] = x*z/pow(R,3.0);
5472		- rotMat1stDerivatives[px][px][ZAxis] = -1.0x/(rR) + xzz/(r*pow(R,3.0));
	5491	+ rotMat1stDerivatives[px][pz][ZAxis] = xz/(RR*R);
	5492	+ rotMat1stDerivatives[px][px][ZAxis] = -1.0x/(rR) + xzz/(rRR*R);
5473	5493
5474	5494	// for d-function
5475	5495	// ToDo: First derivative of rotating matrix for d-orbital...

		@@ -5495,16 +5515,16 @@ void Cndo2::CalcRotatingMatrix2ndDerivatives(double**** rotMat2ndDerivatives,
5495	5515	double x = atomB.GetXyz()[0] - atomA.GetXyz()[0];
5496	5516	double y = atomB.GetXyz()[1] - atomA.GetXyz()[1];
5497	5517	double z = atomB.GetXyz()[2] - atomA.GetXyz()[2];
5498		- double r = sqrt( pow(x,2.0) + pow(y,2.0) );
5499		- double R = sqrt( pow(x,2.0) + pow(y,2.0) + pow(z,2.0) );
	5518	+ double r = sqrt( xx + yy );
	5519	+ double R = sqrt( xx + yy + z*z );
5500	5520
5501	5521	if(r==0e0){
5502	5522	return;
5503	5523	}
5504	5524
5505		- double temp1 = 1.0/(pow(r,3.0)R) + 1.0/(rpow(R,3.0));
5506		- double temp2 = 2.0pow(rR,-3.0) + 3.0/(pow(r,5.0)R) + 3.0/(rpow(R,5.0));
5507		- double temp3 = pow(rR,-3.0) + 3.0/(rpow(R,5.0));
	5525	+ double temp1 = 1.0/(rrrR) + 1.0/(rRRR);
	5526	+ double temp2 = 2.0/(rrrRRR) + 3.0/(rrrrrR) + 3.0/(rRRRR*R);
	5527	+ double temp3 = 1.0/(rrrRRR) + 3.0/(rRRRRR);
5508	5528
5509	5529	// for s-function
5510	5530	rotMat2ndDerivatives[s][s][XAxis][XAxis] = 0.0;

		@@ -5518,29 +5538,29 @@ void Cndo2::CalcRotatingMatrix2ndDerivatives(double**** rotMat2ndDerivatives,
5518	5538	rotMat2ndDerivatives[s][s][ZAxis][ZAxis] = 0.0;
5519	5539
5520	5540	// for p-function, xx-derivatives
5521		- rotMat2ndDerivatives[py][py][XAxis][XAxis] = -3.0xpow(r,-3.0) + 3.0pow(x,3.0)pow(r,-5.0);
5522		- rotMat2ndDerivatives[py][pz][XAxis][XAxis] = -1.0ypow(R,-3.0) + 3.0pow(x,2.0)y*pow(R,-5.0);
5523		- rotMat2ndDerivatives[py][px][XAxis][XAxis] = -1.0temp1yz+temp2pow(x,2.0)yz;
	5541	+ rotMat2ndDerivatives[py][py][XAxis][XAxis] = -3.0x/(rrr) + 3.0xxx/(rrrrr);
	5542	+ rotMat2ndDerivatives[py][pz][XAxis][XAxis] = -1.0y/(RRR) + 3.0xxy/(RRRRR);
	5543	+ rotMat2ndDerivatives[py][px][XAxis][XAxis] = -1.0temp1yz+temp2xxy*z;
5524	5544
5525	5545	rotMat2ndDerivatives[pz][py][XAxis][XAxis] = 0.0;
5526		- rotMat2ndDerivatives[pz][pz][XAxis][XAxis] = -1.0zpow(R,-3.0) + 3.0pow(x,2.0)z*pow(R,-5.0);
5527		- rotMat2ndDerivatives[pz][px][XAxis][XAxis] = -1.0pow(rR,-1.0) + temp1*pow(x,2.0)
5528		- +rpow(R,-3.0) - 3.0pow(x,2.0)rpow(R,-5.0) + pow(x,2.0)pow(r,-1.0)pow(R,-3.0);
	5546	+ rotMat2ndDerivatives[pz][pz][XAxis][XAxis] = -1.0z/(RRR) + 3.0xxz/(RRRRR);
	5547	+ rotMat2ndDerivatives[pz][px][XAxis][XAxis] = -1.0/(rR) + temp1x*x
	5548	+ +r/(RRR) - 3.0xxr/(RRRRR) + xx/(rRR*R);
5529	5549
5530		- rotMat2ndDerivatives[px][py][XAxis][XAxis] = ypow(r,-3.0) - 3.0pow(x,2.0)ypow(r,-5.0);
5531		- rotMat2ndDerivatives[px][pz][XAxis][XAxis] = -3.0xpow(R,-3.0) + 3.0pow(x,3.0)pow(R,-5.0);
5532		- rotMat2ndDerivatives[px][px][XAxis][XAxis] = -3.0temp1xz+temp2pow(x,3.0)*z;
	5550	+ rotMat2ndDerivatives[px][py][XAxis][XAxis] = y/(rrr) - 3.0xxy/(rrrr*r);
	5551	+ rotMat2ndDerivatives[px][pz][XAxis][XAxis] = -3.0x/(RRR) + 3.0xxx/(RRRRR);
	5552	+ rotMat2ndDerivatives[px][px][XAxis][XAxis] = -3.0temp1xz+temp2xxx*z;
5533	5553
5534	5554	// for p-function, xy-derivatives
5535		- rotMat2ndDerivatives[py][py][XAxis][YAxis] = -1.0ypow(r,-3.0) + 3.0pow(x,2.0)y*pow(r,-5.0);
5536		- rotMat2ndDerivatives[py][pz][XAxis][YAxis] = -1.0xpow(R,-3.0) + 3.0xpow(y,2.0)*pow(R,-5.0);
5537		- rotMat2ndDerivatives[py][px][XAxis][YAxis] = -1.0temp1xz+temp2xpow(y,2.0)z;
	5555	+ rotMat2ndDerivatives[py][py][XAxis][YAxis] = -1.0y/(rrr) + 3.0xxy/(rrrrr);
	5556	+ rotMat2ndDerivatives[py][pz][XAxis][YAxis] = -1.0x/(RRR) + 3.0xyy/(RRRRR);
	5557	+ rotMat2ndDerivatives[py][px][XAxis][YAxis] = -1.0temp1xz+temp2xyy*z;
5538	5558
5539	5559	rotMat2ndDerivatives[pz][py][XAxis][YAxis] = 0.0;
5540		- rotMat2ndDerivatives[pz][pz][XAxis][YAxis] = 3.0xyzpow(R,-5.0);
5541		- rotMat2ndDerivatives[pz][px][XAxis][YAxis] = temp1xy + xypow(r,-1.0)pow(R,-3.0) - 3.0xyr*pow(R,-5.0);
	5560	+ rotMat2ndDerivatives[pz][pz][XAxis][YAxis] = 3.0xyz/(RRRR*R);
	5561	+ rotMat2ndDerivatives[pz][px][XAxis][YAxis] = temp1xy + xy/(rRRR) - 3.0xyr/(RRRR*R);
5542	5562
5543		- rotMat2ndDerivatives[px][py][XAxis][YAxis] = xpow(r,-3.0) - 3.0xpow(y,2.0)pow(r,-5.0);
	5563	+ rotMat2ndDerivatives[px][py][XAxis][YAxis] = x/(rrr) - 3.0xyy/(rrrr*r);
5544	5564	rotMat2ndDerivatives[px][pz][XAxis][YAxis] = rotMat2ndDerivatives[py][pz][XAxis][XAxis];
5545	5565	rotMat2ndDerivatives[px][px][XAxis][YAxis] = rotMat2ndDerivatives[py][px][XAxis][XAxis];
5546	5566

		@@ -5553,16 +5573,16 @@ void Cndo2::CalcRotatingMatrix2ndDerivatives(double**** rotMat2ndDerivatives,
5553	5573	// for p-function, xz-derivatives
5554	5574	rotMat2ndDerivatives[py][py][XAxis][ZAxis] = 0.0;
5555	5575	rotMat2ndDerivatives[py][pz][XAxis][ZAxis] = rotMat2ndDerivatives[pz][pz][XAxis][YAxis];
5556		- rotMat2ndDerivatives[py][px][XAxis][ZAxis] = -1.0temp1xy +temp3xypow(z,2.0);
	5576	+ rotMat2ndDerivatives[py][px][XAxis][ZAxis] = -1.0temp1xy +temp3xyz*z;
5557	5577
5558	5578	rotMat2ndDerivatives[pz][py][XAxis][ZAxis] = 0.0;
5559		- rotMat2ndDerivatives[pz][pz][XAxis][ZAxis] = -1.0xpow(R,-3.0) + 3.0xpow(z,2.0)*pow(R,-5.0);
5560		- rotMat2ndDerivatives[pz][px][XAxis][ZAxis] = xzpow(r,-1.0)pow(R,-3.0) - 3.0xzr*pow(R,-5.0);
	5579	+ rotMat2ndDerivatives[pz][pz][XAxis][ZAxis] = -1.0x/(RRR) + 3.0xzz/(RRRRR);
	5580	+ rotMat2ndDerivatives[pz][px][XAxis][ZAxis] = xz/(rRRR) - 3.0xzr/(RRRR*R);
5561	5581
5562	5582	rotMat2ndDerivatives[px][py][XAxis][ZAxis] = 0.0;
5563	5583	rotMat2ndDerivatives[px][pz][XAxis][ZAxis] = rotMat2ndDerivatives[pz][pz][XAxis][XAxis];
5564		- rotMat2ndDerivatives[px][px][XAxis][ZAxis] = pow(rR,-1.0) - pow(z,2.0)pow(r,-1.0)*pow(R,-3.0)
5565		- -1.0temp1pow(x,2.0)+temp3pow(xz,2.0);
	5584	+ rotMat2ndDerivatives[px][px][XAxis][ZAxis] = 1.0/(rR) - zz/(rRR*R)
	5585	+ -1.0temp1xx+temp3xxz*z;
5566	5586
5567	5587
5568	5588	// for p-function, zx-derivatives

		@@ -5573,32 +5593,32 @@ void Cndo2::CalcRotatingMatrix2ndDerivatives(double**** rotMat2ndDerivatives,
5573	5593	}
5574	5594
5575	5595	// for p-function, yy-derivatives
5576		- rotMat2ndDerivatives[py][py][YAxis][YAxis] = -1.0xpow(r,-3.0) + 3.0xpow(y,2.0)*pow(r,-5.0);
5577		- rotMat2ndDerivatives[py][pz][YAxis][YAxis] = -3.0ypow(R,-3.0) + 3.0pow(y,3.0)pow(R,-5.0);
5578		- rotMat2ndDerivatives[py][px][YAxis][YAxis] = -3.0temp1yz+temp2pow(y,3.0)*z;
	5596	+ rotMat2ndDerivatives[py][py][YAxis][YAxis] = -1.0x/(rrr) + 3.0xyy/(rrrrr);
	5597	+ rotMat2ndDerivatives[py][pz][YAxis][YAxis] = -3.0y/(RRR) + 3.0yyy/(RRRRR);
	5598	+ rotMat2ndDerivatives[py][px][YAxis][YAxis] = -3.0temp1yz+temp2yyy*z;
5579	5599
5580	5600	rotMat2ndDerivatives[pz][py][YAxis][YAxis] = 0.0;
5581		- rotMat2ndDerivatives[pz][pz][YAxis][YAxis] = -1.0zpow(R,-3.0) + 3.0pow(y,2.0)z*pow(R,-5.0);
5582		- rotMat2ndDerivatives[pz][px][YAxis][YAxis] = -1.0pow(rR,-1.0) + temp1*pow(y,2.0)
5583		- +rpow(R,-3.0) - 3.0pow(y,2.0)rpow(R,-5.0) + pow(y,2.0)pow(r,-1.0)pow(R,-3.0);
	5601	+ rotMat2ndDerivatives[pz][pz][YAxis][YAxis] = -1.0z/(RRR) + 3.0yyz/(RRRRR);
	5602	+ rotMat2ndDerivatives[pz][px][YAxis][YAxis] = -1.0/(rR) + temp1y*y
	5603	+ +r/(RRR) - 3.0yyr/(RRRRR) + yy/(rRR*R);
5584	5604
5585		- rotMat2ndDerivatives[px][py][YAxis][YAxis] = 3.0ypow(r,-3.0) - 3.0pow(y,3.0)pow(r,-5.0);
	5605	+ rotMat2ndDerivatives[px][py][YAxis][YAxis] = 3.0y/(rrr) - 3.0yyy/(rrrrr);
5586	5606	rotMat2ndDerivatives[px][pz][YAxis][YAxis] = rotMat2ndDerivatives[py][pz][XAxis][YAxis];
5587		- rotMat2ndDerivatives[px][px][YAxis][YAxis] = -1.0temp1xz+temp2xpow(y,2.0)z;
	5607	+ rotMat2ndDerivatives[px][px][YAxis][YAxis] = -1.0temp1xz+temp2xyy*z;
5588	5608
5589	5609	// for p-function, yz-derivatives
5590	5610	rotMat2ndDerivatives[py][py][YAxis][ZAxis] = 0.0;
5591	5611	rotMat2ndDerivatives[py][pz][YAxis][ZAxis] = rotMat2ndDerivatives[pz][pz][YAxis][YAxis];
5592		- rotMat2ndDerivatives[py][px][YAxis][ZAxis] = pow(rR,-1.0) - pow(z,2.0)pow(r,-1.0)*pow(R,-3.0)
5593		- -1.0temp1pow(y,2.0)+temp3pow(yz,2.0);
	5612	+ rotMat2ndDerivatives[py][px][YAxis][ZAxis] = 1.0/(rR) - zz/(rRR*R)
	5613	+ -1.0temp1yy+temp3yyz*z;
5594	5614
5595	5615	rotMat2ndDerivatives[pz][py][YAxis][ZAxis] = 0.0;
5596		- rotMat2ndDerivatives[pz][pz][YAxis][ZAxis] = -1.0ypow(R,-3.0) + 3.0ypow(z,2.0)*pow(R,-5.0);
5597		- rotMat2ndDerivatives[pz][px][YAxis][ZAxis] = yzpow(r,-1.0)pow(R,-3.0) - 3.0yzr*pow(R,-5.0);
	5616	+ rotMat2ndDerivatives[pz][pz][YAxis][ZAxis] = -1.0y/(RRR) + 3.0yzz/(RRRRR);
	5617	+ rotMat2ndDerivatives[pz][px][YAxis][ZAxis] = yz/(rRRR) - 3.0yzr/(RRRR*R);
5598	5618
5599	5619	rotMat2ndDerivatives[px][py][YAxis][ZAxis] = 0.0;
5600	5620	rotMat2ndDerivatives[px][pz][YAxis][ZAxis] = rotMat2ndDerivatives[pz][pz][XAxis][YAxis];
5601		- rotMat2ndDerivatives[px][px][YAxis][ZAxis] = -1.0temp1xy+temp3xypow(z,2.0);
	5621	+ rotMat2ndDerivatives[px][px][YAxis][ZAxis] = -1.0temp1xy+temp3xyz*z;
5602	5622
5603	5623
5604	5624	// for p-function, zy-derivatives

		@@ -5611,15 +5631,15 @@ void Cndo2::CalcRotatingMatrix2ndDerivatives(double**** rotMat2ndDerivatives,
5611	5631	// for p-function, zz-derivatives
5612	5632	rotMat2ndDerivatives[py][py][ZAxis][ZAxis] = 0.0;
5613	5633	rotMat2ndDerivatives[py][pz][ZAxis][ZAxis] = rotMat2ndDerivatives[pz][pz][YAxis][ZAxis];
5614		- rotMat2ndDerivatives[py][px][ZAxis][ZAxis] = -3.0yzpow(r,-1.0)pow(R,-3.0) + 3.0ypow(z,3.0)pow(r,-1.0)pow(R,-5.0);
	5634	+ rotMat2ndDerivatives[py][px][ZAxis][ZAxis] = -3.0yz/(rRRR) + 3.0yzzz/(rRRRRR);
5615	5635
5616	5636	rotMat2ndDerivatives[pz][py][ZAxis][ZAxis] = 0.0;
5617		- rotMat2ndDerivatives[pz][pz][ZAxis][ZAxis] = -3.0zpow(R,-3.0) + 3.0pow(z,3.0)pow(R,-5.0);
5618		- rotMat2ndDerivatives[pz][px][ZAxis][ZAxis] = -3.0pow(z,2.0)rpow(R,-5.0) + rpow(R,-3.0);
	5637	+ rotMat2ndDerivatives[pz][pz][ZAxis][ZAxis] = -3.0z/(RRR) + 3.0zzz/(RRRRR);
	5638	+ rotMat2ndDerivatives[pz][px][ZAxis][ZAxis] = -3.0zzr/(RRRRR) + r/(RR*R);
5619	5639
5620	5640	rotMat2ndDerivatives[px][py][ZAxis][ZAxis] = 0.0;
5621	5641	rotMat2ndDerivatives[px][pz][ZAxis][ZAxis] = rotMat2ndDerivatives[pz][pz][XAxis][ZAxis];
5622		- rotMat2ndDerivatives[px][px][ZAxis][ZAxis] = -3.0xzpow(r,-1.0)pow(R,-3.0) + 3.0xpow(z,3.0)pow(r,-1.0)pow(R,-5.0);
	5642	+ rotMat2ndDerivatives[px][px][ZAxis][ZAxis] = -3.0xz/(rRRR) + 3.0xzzz/(rRRRRR);
5623	5643
5624	5644	// for d-function
5625	5645	// ToDo: Second derivative of rotating matrix for d-orbital...

		@@ -6275,7 +6295,8 @@ double Cndo2::GetAuxiliaryD(int la, int lb, int m) const{
6275	6295	throw MolDSException(ss.str());
6276	6296	}
6277	6297
6278		- double pre = pow(Factorial(m+1)/8.0, 2.0);
	6298	+ double tmp = Factorial(m+1)/8.0;
	6299	+ double pre = tmp*tmp;
6279	6300	double termA = ( (2.0la+1.0)Factorial(la-m) ) / ( 2.0*Factorial(la+m) );
6280	6301	double termB = ( (2.0lb+1.0)Factorial(lb-m) ) / ( 2.0*Factorial(lb+m) );
6281	6302	value = presqrt(termA)sqrt(termB);

--- a/src/mndo/Mndo.cpp

+++ b/src/mndo/Mndo.cpp

		@@ -226,12 +226,12 @@ double Mndo::GetAuxiliaryDiatomCoreRepulsionEnergy2ndDerivative(const Atom& atom
226	226	double pre1=0.0;
227	227	double pre2=0.0;
228	228	if(axisA1 == axisA2){
229		- pre1 = 1.0/distanceAB - pow(dCartesian1,2.0)/pow(distanceAB,3.0);
230		- pre2 = pow(dCartesian1/distanceAB,2.0);
	229	+ pre1 = 1.0/distanceAB - dCartesian1dCartesian1/(distanceABdistanceAB*distanceAB);
	230	+ pre2 = (dCartesian1dCartesian1)/(distanceABdistanceAB);
231	231	}
232	232	else{
233		- pre1 = -dCartesian1*dCartesian2/pow(distanceAB,3.0);
234		- pre2 = dCartesian1*dCartesian2/pow(distanceAB,2.0);
	233	+ pre1 = -dCartesian1dCartesian2/(distanceABdistanceAB*distanceAB);
	234	+ pre2 = dCartesian1dCartesian2/(distanceABdistanceAB);
235	235	}
236	236
237	237	double ang2AU = Parameters::GetInstance()->GetAngstrom2AU();

		@@ -6127,50 +6127,51 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction(const Atom& atomA,
6127	6127
6128	6128	// Eq. (52) in [DT_1977]
6129	6129	if(multipoleA == sQ && multipoleB == sQ){
6130		- value = pow(pow(rAB,2.0) + pow(a,2.0), -0.5);
	6130	+ value = pow((rABrAB) + (aa), -0.5);
6131	6131	}
6132	6132	// Eq. (53) in [DT_1977]
6133	6133	else if(multipoleA == sQ && multipoleB == muz){
6134		- double temp1 = pow(rAB+DB,2.0) + pow(a,2.0);
6135		- double temp2 = pow(rAB-DB,2.0) + pow(a,2.0);
	6134	+ double temp1 = ((rAB+DB)(rAB+DB)) + (aa);
	6135	+ double temp2 = ((rAB-DB)(rAB-DB)) + (aa);
6136	6136	value = pow(temp1,-0.5)/2.0 - pow(temp2,-0.5)/2.0;
6137	6137	}
6138	6138	else if(multipoleA == muz && multipoleB == sQ){
6139	6139	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6140		- value *= pow(-1.0,1.0);
	6140	+ //value *= pow(-1.0,1.0);
	6141	+ value *= -1.0;
6141	6142	}
6142	6143	// Eq. (54) in [DT_1977]
6143	6144	else if(multipoleA == sQ && multipoleB == Qxx){
6144		- double temp1 = pow(rAB,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6145		- double temp2 = pow(rAB,2.0) + pow(a,2.0);
	6145	+ double temp1 = (rABrAB) + (4.0DBDB) + (aa);
	6146	+ double temp2 = (rABrAB) + (aa);
6146	6147	value = pow(temp1,-0.5)/2.0 - pow(temp2,-0.5)/2.0;
6147	6148	}
6148	6149	else if(multipoleA == Qxx && multipoleB == sQ){
6149	6150	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6150		- value *= pow(-1.0,2.0);
	6151	+ //value *= pow(-1.0,2.0);
6151	6152	}
6152	6153	else if(multipoleA == sQ && multipoleB == Qyy){
6153	6154	value = this->GetSemiEmpiricalMultipoleInteraction(atomA, atomB, multipoleA, Qxx, rAB);
6154	6155	}
6155	6156	else if(multipoleA == Qyy && multipoleB == sQ){
6156	6157	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6157		- value *= pow(-1.0,2.0);
	6158	+ //value *= pow(-1.0,2.0);
6158	6159	}
6159	6160	// Eq. (55) in [DT_1977]
6160	6161	else if(multipoleA == sQ && multipoleB == Qzz){
6161		- double temp1 = pow(rAB+2.0*DB,2.0) + pow(a,2.0);
6162		- double temp2 = pow(rAB,2.0) + pow(a,2.0);
6163		- double temp3 = pow(rAB-2.0*DB,2.0) + pow(a,2.0);
	6162	+ double temp1 = ((rAB+2.0DB)(rAB+2.0DB)) + (aa);
	6163	+ double temp2 = (rABrAB) + (aa);
	6164	+ double temp3 = ((rAB-2.0DB)(rAB-2.0DB)) + (aa);
6164	6165	value = pow(temp1,-0.5)/4.0 - pow(temp2,-0.5)/2.0 + pow(temp3,-0.5)/4.0;
6165	6166	}
6166	6167	else if(multipoleA == Qzz && multipoleB == sQ){
6167	6168	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6168		- value *= pow(-1.0,2.0);
	6169	+ //value *= pow(-1.0,2.0);
6169	6170	}
6170	6171	// Eq. (56) in [DT_1977]
6171	6172	else if(multipoleA == mux && multipoleB == mux){
6172		- double temp1 = pow(rAB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6173		- double temp2 = pow(rAB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
	6173	+ double temp1 = (rABrAB) + ((DA-DB)(DA-DB)) + (a*a);
	6174	+ double temp2 = (rABrAB) + ((DA+DB)(DA+DB)) + (a*a);
6174	6175	value = pow(temp1,-0.5)/2.0 - pow(temp2,-0.5)/2.0;
6175	6176	}
6176	6177	else if(multipoleA == muy && multipoleB == muy){

		@@ -6178,76 +6179,81 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction(const Atom& atomA,
6178	6179	}
6179	6180	// Eq. (57) in [DT_1977]
6180	6181	else if(multipoleA == muz && multipoleB == muz){
6181		- double temp1 = pow(rAB+DA-DB,2.0) + pow(a,2.0);
6182		- double temp2 = pow(rAB+DA+DB,2.0) + pow(a,2.0);
6183		- double temp3 = pow(rAB-DA-DB,2.0) + pow(a,2.0);
6184		- double temp4 = pow(rAB-DA+DB,2.0) + pow(a,2.0);
	6182	+ double temp1 = ((rAB+DA-DB)(rAB+DA-DB)) + (aa);
	6183	+ double temp2 = ((rAB+DA+DB)(rAB+DA+DB)) + (aa);
	6184	+ double temp3 = ((rAB-DA-DB)(rAB-DA-DB)) + (aa);
	6185	+ double temp4 = ((rAB-DA+DB)(rAB-DA+DB)) + (aa);
6185	6186	value = pow(temp1,-0.5)/4.0 - pow(temp2,-0.5)/4.0
6186	6187	-pow(temp3,-0.5)/4.0 + pow(temp4,-0.5)/4.0;
6187	6188	}
6188	6189	// Eq. (58) in [DT_1977]
6189	6190	else if(multipoleA == mux && multipoleB == Qxz){
6190		- double temp1 = pow(rAB-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6191		- double temp2 = pow(rAB-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
6192		- double temp3 = pow(rAB+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6193		- double temp4 = pow(rAB+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
	6191	+ double temp1 = ((rAB-DB)(rAB-DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6192	+ double temp2 = ((rAB-DB)(rAB-DB)) + ((DA+DB)(DA+DB)) + (a*a);
	6193	+ double temp3 = ((rAB+DB)(rAB+DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6194	+ double temp4 = ((rAB+DB)(rAB+DB)) + ((DA+DB)(DA+DB)) + (a*a);
6194	6195	value =-pow(temp1,-0.5)/4.0 + pow(temp2,-0.5)/4.0
6195	6196	+pow(temp3,-0.5)/4.0 - pow(temp4,-0.5)/4.0;
6196	6197	}
6197	6198	else if(multipoleA == Qxz && multipoleB == mux){
6198	6199	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6199		- value *= pow(-1.0,3.0);
	6200	+ //value *= pow(-1.0,3.0);
	6201	+ value *= -1.0;
6200	6202	}
6201	6203	else if(multipoleA == muy && multipoleB == Qyz){
6202	6204	value = this->GetSemiEmpiricalMultipoleInteraction(atomA, atomB, mux, Qxz, rAB);
6203	6205	}
6204	6206	else if(multipoleA == Qyz && multipoleB == muy){
6205	6207	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6206		- value *= pow(-1.0,3.0);
	6208	+ //value *= pow(-1.0,3.0);
	6209	+ value *= -1.0;
6207	6210	}
6208	6211	// Eq. (59) in [DT_1977]
6209	6212	else if(multipoleA == muz && multipoleB == Qxx){
6210		- double temp1 = pow(rAB+DA,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6211		- double temp2 = pow(rAB-DA,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6212		- double temp3 = pow(rAB+DA,2.0) + pow(a,2.0);
6213		- double temp4 = pow(rAB-DA,2.0) + pow(a,2.0);
	6213	+ double temp1 = ((rAB+DA)(rAB+DA)) + (4.0DBDB) + (aa);
	6214	+ double temp2 = ((rAB-DA)(rAB-DA)) + (4.0DBDB) + (aa);
	6215	+ double temp3 = ((rAB+DA)(rAB+DA)) + (aa);
	6216	+ double temp4 = ((rAB-DA)(rAB-DA)) + (aa);
6214	6217	value =-pow(temp1,-0.5)/4.0 + pow(temp2,-0.5)/4.0
6215	6218	+pow(temp3,-0.5)/4.0 - pow(temp4,-0.5)/4.0;
6216	6219	}
6217	6220	else if(multipoleA == Qxx && multipoleB == muz){
6218	6221	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6219		- value *= pow(-1.0,3.0);
	6222	+ //value *= pow(-1.0,3.0);
	6223	+ value *= -1.0;
6220	6224	}
6221	6225	else if(multipoleA == muz && multipoleB == Qyy){
6222	6226	value = this->GetSemiEmpiricalMultipoleInteraction(atomA, atomB, muz, Qxx, rAB);
6223	6227	}
6224	6228	else if(multipoleA == Qyy && multipoleB == muz){
6225	6229	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6226		- value *= pow(-1.0,3.0);
	6230	+ //value *= pow(-1.0,3.0);
	6231	+ value *= -1.0;
6227	6232	}
6228	6233	// Eq. (60) in [DT_1977]
6229	6234	else if(multipoleA == muz && multipoleB == Qzz){
6230		- double temp1 = pow(rAB+DA-2.0*DB,2.0) + pow(a,2.0);
6231		- double temp2 = pow(rAB-DA-2.0*DB,2.0) + pow(a,2.0);
6232		- double temp3 = pow(rAB+DA+2.0*DB,2.0) + pow(a,2.0);
6233		- double temp4 = pow(rAB-DA+2.0*DB,2.0) + pow(a,2.0);
6234		- double temp5 = pow(rAB+DA,2.0) + pow(a,2.0);
6235		- double temp6 = pow(rAB-DA,2.0) + pow(a,2.0);
	6235	+ double temp1 = ((rAB+DA-2.0DB)(rAB+DA-2.0DB)) + (aa);
	6236	+ double temp2 = ((rAB-DA-2.0DB)(rAB-DA-2.0DB)) + (aa);
	6237	+ double temp3 = ((rAB+DA+2.0DB)(rAB+DA+2.0DB)) + (aa);
	6238	+ double temp4 = ((rAB-DA+2.0DB)(rAB-DA+2.0DB)) + (aa);
	6239	+ double temp5 = ((rAB+DA)(rAB+DA)) + (aa);
	6240	+ double temp6 = ((rAB-DA)(rAB-DA)) + (aa);
6236	6241	value =-pow(temp1,-0.5)/8.0 + pow(temp2,-0.5)/8.0
6237	6242	-pow(temp3,-0.5)/8.0 + pow(temp4,-0.5)/8.0
6238	6243	+pow(temp5,-0.5)/4.0 - pow(temp6,-0.5)/4.0;
6239	6244	}
6240	6245	else if(multipoleA == Qzz && multipoleB == muz){
6241	6246	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6242		- value *= pow(-1.0,3.0);
	6247	+ //value *= pow(-1.0,3.0);
	6248	+ value *= -1.0;
6243	6249	}
6244	6250	// Eq. (61) in [DT_1977]
6245	6251	else if(multipoleA == Qxx && multipoleB == Qxx){
6246		- double temp1 = pow(rAB,2.0) + 4.0*pow(DA-DB,2.0) + pow(a,2.0);
6247		- double temp2 = pow(rAB,2.0) + 4.0*pow(DA+DB,2.0) + pow(a,2.0);
6248		- double temp3 = pow(rAB,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
6249		- double temp4 = pow(rAB,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6250		- double temp5 = pow(rAB,2.0) + pow(a,2.0);
	6252	+ double temp1 = (rABrAB) + 4.0((DA-DB)(DA-DB)) + (aa);
	6253	+ double temp2 = (rABrAB) + 4.0((DA+DB)(DA+DB)) + (aa);
	6254	+ double temp3 = (rABrAB) + (4.0DADA) + (aa);
	6255	+ double temp4 = (rABrAB) + (4.0DBDB) + (aa);
	6256	+ double temp5 = (rABrAB) + (aa);
6251	6257	value = pow(temp1,-0.5)/8.0 + pow(temp2,-0.5)/8.0
6252	6258	-pow(temp3,-0.5)/4.0 - pow(temp4,-0.5)/4.0
6253	6259	+pow(temp5,-0.5)/4.0;

		@@ -6257,51 +6263,51 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction(const Atom& atomA,
6257	6263	}
6258	6264	// Eq. (62) in [DT_1977]
6259	6265	else if(multipoleA == Qxx && multipoleB == Qyy){
6260		- double temp1 = pow(rAB,2.0) + pow(2.0DA,2.0) + pow(2.0DB,2.0)+ pow(a,2.0);
6261		- double temp2 = pow(rAB,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
6262		- double temp3 = pow(rAB,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6263		- double temp4 = pow(rAB,2.0) + pow(a,2.0);
	6266	+ double temp1 = (rABrAB) + (4.0DADA) + (4.0DBDB)+ (aa);
	6267	+ double temp2 = (rABrAB) + (4.0DADA) + (aa);
	6268	+ double temp3 = (rABrAB) + (4.0DBDB) + (aa);
	6269	+ double temp4 = (rABrAB) + (aa);
6264	6270	value = pow(temp1,-0.5)/4.0 - pow(temp2,-0.5)/4.0
6265	6271	-pow(temp3,-0.5)/4.0 + pow(temp4,-0.5)/4.0;
6266	6272	}
6267	6273	else if(multipoleA == Qyy && multipoleB == Qxx){
6268	6274	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6269		- value *= pow(-1.0,4.0);
	6275	+ //value *= pow(-1.0,4.0);
6270	6276	}
6271	6277	// Eq. (63) in [DT_1977]
6272	6278	else if(multipoleA == Qxx && multipoleB == Qzz){
6273		- double temp1 = pow(rAB-2.0DB,2.0) + pow(2.0DA,2.0) + pow(a,2.0);
6274		- double temp2 = pow(rAB+2.0DB,2.0) + pow(2.0DA,2.0) + pow(a,2.0);
6275		- double temp3 = pow(rAB-2.0*DB,2.0) + pow(a,2.0);
6276		- double temp4 = pow(rAB+2.0*DB,2.0) + pow(a,2.0);
6277		- double temp5 = pow(rAB,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
6278		- double temp6 = pow(rAB,2.0) + pow(a,2.0);
	6279	+ double temp1 = ((rAB-2.0DB)(rAB-2.0DB)) + (4.0DADA) + (aa);
	6280	+ double temp2 = ((rAB+2.0DB)(rAB+2.0DB)) + (4.0DADA) + (aa);
	6281	+ double temp3 = ((rAB-2.0DB)(rAB-2.0DB)) + (aa);
	6282	+ double temp4 = ((rAB+2.0DB)(rAB+2.0DB)) + (aa);
	6283	+ double temp5 = (rABrAB) + (4.0DADA) + (aa);
	6284	+ double temp6 = (rABrAB) + (aa);
6279	6285	value = pow(temp1,-0.5)/8.0 + pow(temp2,-0.5)/8.0
6280	6286	-pow(temp3,-0.5)/8.0 - pow(temp4,-0.5)/8.0
6281	6287	-pow(temp5,-0.5)/4.0 + pow(temp6,-0.5)/4.0;
6282	6288	}
6283	6289	else if(multipoleA == Qzz && multipoleB == Qxx){
6284	6290	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6285		- value *= pow(-1.0,4.0);
	6291	+ //value *= pow(-1.0,4.0);
6286	6292	}
6287	6293	else if(multipoleA == Qyy && multipoleB == Qzz){
6288	6294	value = this->GetSemiEmpiricalMultipoleInteraction(atomA, atomB, Qxx, multipoleB, rAB);
6289	6295	}
6290	6296	else if(multipoleA == Qzz && multipoleB == Qyy){
6291	6297	value = this->GetSemiEmpiricalMultipoleInteraction(atomB, atomA, multipoleB, multipoleA, rAB);
6292		- value *= pow(-1.0,4.0);
	6298	+ //value *= pow(-1.0,4.0);
6293	6299	}
6294	6300	// Eq. (64) in [DT_1977]
6295	6301	else if(multipoleA == Qzz && multipoleB == Qzz){
6296		- double temp1 = pow(rAB+2.0DA-2.0DB,2.0) + pow(a,2.0);
6297		- double temp2 = pow(rAB+2.0DA+2.0DB,2.0) + pow(a,2.0);
6298		- double temp3 = pow(rAB-2.0DA-2.0DB,2.0) + pow(a,2.0);
6299		- double temp4 = pow(rAB-2.0DA+2.0DB,2.0) + pow(a,2.0);
6300		- double temp5 = pow(rAB+2.0*DA,2.0) + pow(a,2.0);
6301		- double temp6 = pow(rAB-2.0*DA,2.0) + pow(a,2.0);
6302		- double temp7 = pow(rAB+2.0*DB,2.0) + pow(a,2.0);
6303		- double temp8 = pow(rAB-2.0*DB,2.0) + pow(a,2.0);
6304		- double temp9 = pow(rAB,2.0) + pow(a,2.0);
	6302	+ double temp1 = ((rAB+2.0DA-2.0DB)(rAB+2.0DA-2.0DB)) + (aa);
	6303	+ double temp2 = ((rAB+2.0DA+2.0DB)(rAB+2.0DA+2.0DB)) + (aa);
	6304	+ double temp3 = ((rAB-2.0DA-2.0DB)(rAB-2.0DA-2.0DB)) + (aa);
	6305	+ double temp4 = ((rAB-2.0DA+2.0DB)(rAB-2.0DA+2.0DB)) + (aa);
	6306	+ double temp5 = pow(rAB+2.0DA,2.0) + (aa);
	6307	+ double temp6 = pow(rAB-2.0DA,2.0) + (aa);
	6308	+ double temp7 = ((rAB+2.0DB)(rAB+2.0DB)) + (aa);
	6309	+ double temp8 = ((rAB-2.0DB)(rAB-2.0DB)) + (aa);
	6310	+ double temp9 = (rABrAB) + (aa);
6305	6311	value = pow(temp1,-0.5)/16.0 + pow(temp2,-0.5)/16.0
6306	6312	+pow(temp3,-0.5)/16.0 + pow(temp4,-0.5)/16.0
6307	6313	-pow(temp5,-0.5)/8.0 - pow(temp6,-0.5)/8.0

		@@ -6310,14 +6316,14 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction(const Atom& atomA,
6310	6316	}
6311	6317	// Eq. (65) in [DT_1977]
6312	6318	else if(multipoleA == Qxz && multipoleB == Qxz){
6313		- double temp1 = pow(rAB+DA-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6314		- double temp2 = pow(rAB+DA-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
6315		- double temp3 = pow(rAB+DA+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6316		- double temp4 = pow(rAB+DA+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
6317		- double temp5 = pow(rAB-DA-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6318		- double temp6 = pow(rAB-DA-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
6319		- double temp7 = pow(rAB-DA+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6320		- double temp8 = pow(rAB-DA+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
	6319	+ double temp1 = ((rAB+DA-DB)(rAB+DA-DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6320	+ double temp2 = ((rAB+DA-DB)(rAB+DA-DB)) + ((DA+DB)(DA+DB)) + (a*a);
	6321	+ double temp3 = ((rAB+DA+DB)(rAB+DA+DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6322	+ double temp4 = ((rAB+DA+DB)(rAB+DA+DB)) + ((DA+DB)(DA+DB)) + (a*a);
	6323	+ double temp5 = ((rAB-DA-DB)(rAB-DA-DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6324	+ double temp6 = ((rAB-DA-DB)(rAB-DA-DB)) + ((DA+DB)(DA+DB)) + (a*a);
	6325	+ double temp7 = ((rAB-DA+DB)(rAB-DA+DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6326	+ double temp8 = ((rAB-DA+DB)(rAB-DA+DB)) + ((DA+DB)(DA+DB)) + (a*a);
6321	6327	value = pow(temp1,-0.5)/8.0 - pow(temp2,-0.5)/8.0
6322	6328	-pow(temp3,-0.5)/8.0 + pow(temp4,-0.5)/8.0
6323	6329	-pow(temp5,-0.5)/8.0 + pow(temp6,-0.5)/8.0

		@@ -6328,9 +6334,9 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction(const Atom& atomA,
6328	6334	}
6329	6335	// Eq. (66) in [DT_1977]
6330	6336	else if(multipoleA == Qxy && multipoleB == Qxy){
6331		- double temp1 = pow(rAB,2.0) + 2.0*pow(DA-DB,2.0) + pow(a,2.0);
6332		- double temp2 = pow(rAB,2.0) + 2.0*pow(DA+DB,2.0) + pow(a,2.0);
6333		- double temp3 = pow(rAB,2.0) + 2.0pow(DA,2.0) + 2.0pow(DB,2.0) + pow(a,2.0);
	6337	+ double temp1 = (rABrAB) + 2.0((DA-DB)(DA-DB)) + (aa);
	6338	+ double temp2 = (rABrAB) + 2.0((DA+DB)(DA+DB)) + (aa);
	6339	+ double temp3 = (rABrAB) + 2.0(DADA) + 2.0(DBDB) + (aa);
6334	6340	value = pow(temp1,-0.5)/4.0 + pow(temp2,-0.5)/4.0
6335	6341	-pow(temp3,-0.5)/2.0;
6336	6342	}

		@@ -6361,44 +6367,45 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6361	6367
6362	6368	// Eq. (52) in [DT_1977]
6363	6369	if(multipoleA == sQ && multipoleB == sQ){
6364		- value = -1.0rABpow(pow(rAB,2.0) + pow(a,2.0), -1.5);
	6370	+ value = -1.0rABpow((rABrAB) + (aa), -1.5);
6365	6371	}
6366	6372	// Eq. (53) in [DT_1977]
6367	6373	else if(multipoleA == sQ && multipoleB == muz){
6368		- double temp1 = pow(rAB+DB,2.0) + pow(a,2.0);
6369		- double temp2 = pow(rAB-DB,2.0) + pow(a,2.0);
	6374	+ double temp1 = ((rAB+DB)(rAB+DB)) + (aa);
	6375	+ double temp2 = ((rAB-DB)(rAB-DB)) + (aa);
6370	6376	value = (rAB+DB)*pow(temp1,-1.5)/2.0
6371	6377	-(rAB-DB)*pow(temp2,-1.5)/2.0;
6372	6378	value *= -1.0;
6373	6379	}
6374	6380	else if(multipoleA == muz && multipoleB == sQ){
6375	6381	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6376		- value *= pow(-1.0,1.0);
	6382	+ //value *= pow(-1.0,1.0);
	6383	+ value *= -1.0;
6377	6384	}
6378	6385	// Eq. (54) in [DT_1977]
6379	6386	else if(multipoleA == sQ && multipoleB == Qxx){
6380		- double temp1 = pow(rAB,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6381		- double temp2 = pow(rAB,2.0) + pow(a,2.0);
	6387	+ double temp1 = (rABrAB) + (4.0DBDB) + (aa);
	6388	+ double temp2 = (rABrAB) + (aa);
6382	6389	value = rAB*pow(temp1,-1.5)/2.0
6383	6390	-rAB*pow(temp2,-1.5)/2.0;
6384	6391	value *= -1.0;
6385	6392	}
6386	6393	else if(multipoleA == Qxx && multipoleB == sQ){
6387	6394	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6388		- value *= pow(-1.0,2.0);
	6395	+ //value *= pow(-1.0,2.0);
6389	6396	}
6390	6397	else if(multipoleA == sQ && multipoleB == Qyy){
6391	6398	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomA, atomB, multipoleA, Qxx, rAB);
6392	6399	}
6393	6400	else if(multipoleA == Qyy && multipoleB == sQ){
6394	6401	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6395		- value *= pow(-1.0,2.0);
	6402	+ //value *= pow(-1.0,2.0);
6396	6403	}
6397	6404	// Eq. (55) in [DT_1977]
6398	6405	else if(multipoleA == sQ && multipoleB == Qzz){
6399		- double temp1 = pow(rAB+2.0*DB,2.0) + pow(a,2.0);
6400		- double temp2 = pow(rAB,2.0) + pow(a,2.0);
6401		- double temp3 = pow(rAB-2.0*DB,2.0) + pow(a,2.0);
	6406	+ double temp1 = ((rAB+2.0DB)(rAB+2.0DB)) + (aa);
	6407	+ double temp2 = (rABrAB) + (aa);
	6408	+ double temp3 = ((rAB-2.0DB)(rAB-2.0DB)) + (aa);
6402	6409	value = (rAB+2.0DB)pow(temp1,-1.5)/4.0
6403	6410	-(rAB)*pow(temp2,-1.5)/2.0
6404	6411	+(rAB-2.0DB)pow(temp3,-1.5)/4.0;

		@@ -6406,12 +6413,12 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6406	6413	}
6407	6414	else if(multipoleA == Qzz && multipoleB == sQ){
6408	6415	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6409		- value *= pow(-1.0,2.0);
	6416	+ //value *= pow(-1.0,2.0);
6410	6417	}
6411	6418	// Eq. (56) in [DT_1977]
6412	6419	else if(multipoleA == mux && multipoleB == mux){
6413		- double temp1 = pow(rAB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6414		- double temp2 = pow(rAB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
	6420	+ double temp1 = (rABrAB) + ((DA-DB)(DA-DB)) + (a*a);
	6421	+ double temp2 = (rABrAB) + ((DA+DB)(DA+DB)) + (a*a);
6415	6422	value = (rAB)*pow(temp1,-1.5)/2.0
6416	6423	-(rAB)*pow(temp2,-1.5)/2.0;
6417	6424	value *= -1.0;

		@@ -6421,10 +6428,10 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6421	6428	}
6422	6429	// Eq. (57) in [DT_1977]
6423	6430	else if(multipoleA == muz && multipoleB == muz){
6424		- double temp1 = pow(rAB+DA-DB,2.0) + pow(a,2.0);
6425		- double temp2 = pow(rAB+DA+DB,2.0) + pow(a,2.0);
6426		- double temp3 = pow(rAB-DA-DB,2.0) + pow(a,2.0);
6427		- double temp4 = pow(rAB-DA+DB,2.0) + pow(a,2.0);
	6431	+ double temp1 = ((rAB+DA-DB)(rAB+DA-DB)) + (aa);
	6432	+ double temp2 = ((rAB+DA+DB)(rAB+DA+DB)) + (aa);
	6433	+ double temp3 = ((rAB-DA-DB)(rAB-DA-DB)) + (aa);
	6434	+ double temp4 = ((rAB-DA+DB)(rAB-DA+DB)) + (aa);
6428	6435	value = (rAB+DA-DB)*pow(temp1,-1.5)/4.0
6429	6436	-(rAB+DA+DB)*pow(temp2,-1.5)/4.0
6430	6437	-(rAB-DA-DB)*pow(temp3,-1.5)/4.0

		@@ -6433,10 +6440,10 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6433	6440	}
6434	6441	// Eq. (58) in [DT_1977]
6435	6442	else if(multipoleA == mux && multipoleB == Qxz){
6436		- double temp1 = pow(rAB-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6437		- double temp2 = pow(rAB-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
6438		- double temp3 = pow(rAB+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6439		- double temp4 = pow(rAB+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
	6443	+ double temp1 = ((rAB-DB)(rAB-DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6444	+ double temp2 = ((rAB-DB)(rAB-DB)) + ((DA+DB)(DA+DB)) + (a*a);
	6445	+ double temp3 = ((rAB+DB)(rAB+DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6446	+ double temp4 = ((rAB+DB)(rAB+DB)) + ((DA+DB)(DA+DB)) + (a*a);
6440	6447	value =-(rAB-DB)*pow(temp1,-1.5)/4.0
6441	6448	+(rAB-DB)*pow(temp2,-1.5)/4.0
6442	6449	+(rAB+DB)*pow(temp3,-1.5)/4.0

		@@ -6445,21 +6452,23 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6445	6452	}
6446	6453	else if(multipoleA == Qxz && multipoleB == mux){
6447	6454	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6448		- value *= pow(-1.0,3.0);
	6455	+ //value *= pow(-1.0,3.0);
	6456	+ value *= -1.0;
6449	6457	}
6450	6458	else if(multipoleA == muy && multipoleB == Qyz){
6451	6459	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomA, atomB, mux, Qxz, rAB);
6452	6460	}
6453	6461	else if(multipoleA == Qyz && multipoleB == muy){
6454	6462	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6455		- value *= pow(-1.0,3.0);
	6463	+ //value *= pow(-1.0,3.0);
	6464	+ value *= -1.0;
6456	6465	}
6457	6466	// Eq. (59) in [DT_1977]
6458	6467	else if(multipoleA == muz && multipoleB == Qxx){
6459		- double temp1 = pow(rAB+DA,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6460		- double temp2 = pow(rAB-DA,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6461		- double temp3 = pow(rAB+DA,2.0) + pow(a,2.0);
6462		- double temp4 = pow(rAB-DA,2.0) + pow(a,2.0);
	6468	+ double temp1 = ((rAB+DA)(rAB+DA)) + (4.0DBDB) + (aa);
	6469	+ double temp2 = ((rAB-DA)(rAB-DA)) + (4.0DBDB) + (aa);
	6470	+ double temp3 = ((rAB+DA)(rAB+DA)) + (aa);
	6471	+ double temp4 = ((rAB-DA)(rAB-DA)) + (aa);
6463	6472	value =-(rAB+DA)*pow(temp1,-1.5)/4.0
6464	6473	+(rAB-DA)*pow(temp2,-1.5)/4.0
6465	6474	+(rAB+DA)*pow(temp3,-1.5)/4.0

		@@ -6468,23 +6477,25 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6468	6477	}
6469	6478	else if(multipoleA == Qxx && multipoleB == muz){
6470	6479	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6471		- value *= pow(-1.0,3.0);
	6480	+ //value *= pow(-1.0,3.0);
	6481	+ value *= -1.0;
6472	6482	}
6473	6483	else if(multipoleA == muz && multipoleB == Qyy){
6474	6484	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomA, atomB, muz, Qxx, rAB);
6475	6485	}
6476	6486	else if(multipoleA == Qyy && multipoleB == muz){
6477	6487	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6478		- value *= pow(-1.0,3.0);
	6488	+ //value *= pow(-1.0,3.0);
	6489	+ value *= -1.0;
6479	6490	}
6480	6491	// Eq. (60) in [DT_1977]
6481	6492	else if(multipoleA == muz && multipoleB == Qzz){
6482		- double temp1 = pow(rAB+DA-2.0*DB,2.0) + pow(a,2.0);
6483		- double temp2 = pow(rAB-DA-2.0*DB,2.0) + pow(a,2.0);
6484		- double temp3 = pow(rAB+DA+2.0*DB,2.0) + pow(a,2.0);
6485		- double temp4 = pow(rAB-DA+2.0*DB,2.0) + pow(a,2.0);
6486		- double temp5 = pow(rAB+DA,2.0) + pow(a,2.0);
6487		- double temp6 = pow(rAB-DA,2.0) + pow(a,2.0);
	6493	+ double temp1 = ((rAB+DA-2.0DB)(rAB+DA-2.0DB)) + (aa);
	6494	+ double temp2 = ((rAB-DA-2.0DB)(rAB-DA-2.0DB)) + (aa);
	6495	+ double temp3 = ((rAB+DA+2.0DB)(rAB+DA+2.0DB)) + (aa);
	6496	+ double temp4 = ((rAB-DA+2.0DB)(rAB-DA+2.0DB)) + (aa);
	6497	+ double temp5 = ((rAB+DA)(rAB+DA)) + (aa);
	6498	+ double temp6 = ((rAB-DA)(rAB-DA)) + (aa);
6488	6499	value =-(rAB+DA-2.0DB)pow(temp1,-1.5)/8.0
6489	6500	+(rAB-DA-2.0DB)pow(temp2,-1.5)/8.0
6490	6501	-(rAB+DA+2.0DB)pow(temp3,-1.5)/8.0

		@@ -6495,15 +6506,16 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6495	6506	}
6496	6507	else if(multipoleA == Qzz && multipoleB == muz){
6497	6508	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6498		- value *= pow(-1.0,3.0);
	6509	+ //value *= pow(-1.0,3.0);
	6510	+ value *= -1.0;
6499	6511	}
6500	6512	// Eq. (61) in [DT_1977]
6501	6513	else if(multipoleA == Qxx && multipoleB == Qxx){
6502		- double temp1 = pow(rAB,2.0) + 4.0*pow(DA-DB,2.0) + pow(a,2.0);
6503		- double temp2 = pow(rAB,2.0) + 4.0*pow(DA+DB,2.0) + pow(a,2.0);
6504		- double temp3 = pow(rAB,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
6505		- double temp4 = pow(rAB,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6506		- double temp5 = pow(rAB,2.0) + pow(a,2.0);
	6514	+ double temp1 = (rABrAB) + 4.0((DA-DB)(DA-DB)) + (aa);
	6515	+ double temp2 = (rABrAB) + 4.0((DA+DB)(DA+DB)) + (aa);
	6516	+ double temp3 = (rABrAB) + (4.0DADA) + (aa);
	6517	+ double temp4 = (rABrAB) + (4.0DBDB) + (aa);
	6518	+ double temp5 = (rABrAB) + (aa);
6507	6519	value = (rAB)*pow(temp1,-1.5)/8.0
6508	6520	+(rAB)*pow(temp2,-1.5)/8.0
6509	6521	-(rAB)*pow(temp3,-1.5)/4.0

		@@ -6516,10 +6528,10 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6516	6528	}
6517	6529	// Eq. (62) in [DT_1977]
6518	6530	else if(multipoleA == Qxx && multipoleB == Qyy){
6519		- double temp1 = pow(rAB,2.0) + pow(2.0DA,2.0) + pow(2.0DB,2.0)+ pow(a,2.0);
6520		- double temp2 = pow(rAB,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
6521		- double temp3 = pow(rAB,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
6522		- double temp4 = pow(rAB,2.0) + pow(a,2.0);
	6531	+ double temp1 = (rABrAB) + (4.0DADA) + (4.0DBDB)+ (aa);
	6532	+ double temp2 = (rABrAB) + (4.0DADA) + (aa);
	6533	+ double temp3 = (rABrAB) + (4.0DBDB) + (aa);
	6534	+ double temp4 = (rABrAB) + (aa);
6523	6535	value = (rAB)*pow(temp1,-1.5)/4.0
6524	6536	-(rAB)*pow(temp2,-1.5)/4.0
6525	6537	-(rAB)*pow(temp3,-1.5)/4.0

		@@ -6528,16 +6540,16 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6528	6540	}
6529	6541	else if(multipoleA == Qyy && multipoleB == Qxx){
6530	6542	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6531		- value *= pow(-1.0,4.0);
	6543	+ //value *= pow(-1.0,4.0);
6532	6544	}
6533	6545	// Eq. (63) in [DT_1977]
6534	6546	else if(multipoleA == Qxx && multipoleB == Qzz){
6535		- double temp1 = pow(rAB-2.0DB,2.0) + pow(2.0DA,2.0) + pow(a,2.0);
6536		- double temp2 = pow(rAB+2.0DB,2.0) + pow(2.0DA,2.0) + pow(a,2.0);
6537		- double temp3 = pow(rAB-2.0*DB,2.0) + pow(a,2.0);
6538		- double temp4 = pow(rAB+2.0*DB,2.0) + pow(a,2.0);
6539		- double temp5 = pow(rAB,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
6540		- double temp6 = pow(rAB,2.0) + pow(a,2.0);
	6547	+ double temp1 = ((rAB-2.0DB)(rAB-2.0DB)) + (4.0DADA) + (aa);
	6548	+ double temp2 = ((rAB+2.0DB)(rAB+2.0DB)) + (4.0DADA) + (aa);
	6549	+ double temp3 = ((rAB-2.0DB)(rAB-2.0DB)) + (aa);
	6550	+ double temp4 = ((rAB+2.0DB)(rAB+2.0DB)) + (aa);
	6551	+ double temp5 = (rABrAB) + (4.0DADA) + (aa);
	6552	+ double temp6 = (rABrAB) + (aa);
6541	6553	value = (rAB-2.0DB)pow(temp1,-1.5)/8.0
6542	6554	+(rAB+2.0DB)pow(temp2,-1.5)/8.0
6543	6555	-(rAB-2.0DB)pow(temp3,-1.5)/8.0

		@@ -6548,26 +6560,26 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6548	6560	}
6549	6561	else if(multipoleA == Qzz && multipoleB == Qxx){
6550	6562	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6551		- value *= pow(-1.0,4.0);
	6563	+ //value *= pow(-1.0,4.0);
6552	6564	}
6553	6565	else if(multipoleA == Qyy && multipoleB == Qzz){
6554	6566	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomA, atomB, Qxx, multipoleB, rAB);
6555	6567	}
6556	6568	else if(multipoleA == Qzz && multipoleB == Qyy){
6557	6569	value = this->GetSemiEmpiricalMultipoleInteraction1stDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6558		- value *= pow(-1.0,4.0);
	6570	+ //value *= pow(-1.0,4.0);
6559	6571	}
6560	6572	// Eq. (64) in [DT_1977]
6561	6573	else if(multipoleA == Qzz && multipoleB == Qzz){
6562		- double temp1 = pow(rAB+2.0DA-2.0DB,2.0) + pow(a,2.0);
6563		- double temp2 = pow(rAB+2.0DA+2.0DB,2.0) + pow(a,2.0);
6564		- double temp3 = pow(rAB-2.0DA-2.0DB,2.0) + pow(a,2.0);
6565		- double temp4 = pow(rAB-2.0DA+2.0DB,2.0) + pow(a,2.0);
6566		- double temp5 = pow(rAB+2.0*DA,2.0) + pow(a,2.0);
6567		- double temp6 = pow(rAB-2.0*DA,2.0) + pow(a,2.0);
6568		- double temp7 = pow(rAB+2.0*DB,2.0) + pow(a,2.0);
6569		- double temp8 = pow(rAB-2.0*DB,2.0) + pow(a,2.0);
6570		- double temp9 = pow(rAB,2.0) + pow(a,2.0);
	6574	+ double temp1 = ((rAB+2.0DA-2.0DB)(rAB+2.0DA-2.0DB)) + (aa);
	6575	+ double temp2 = ((rAB+2.0DA+2.0DB)(rAB+2.0DA+2.0DB)) + (aa);
	6576	+ double temp3 = ((rAB-2.0DA-2.0DB)(rAB-2.0DA-2.0DB)) + (aa);
	6577	+ double temp4 = ((rAB-2.0DA+2.0DB)(rAB-2.0DA+2.0DB)) + (aa);
	6578	+ double temp5 = pow(rAB+2.0DA,2.0) + (aa);
	6579	+ double temp6 = pow(rAB-2.0DA,2.0) + (aa);
	6580	+ double temp7 = ((rAB+2.0DB)(rAB+2.0DB)) + (aa);
	6581	+ double temp8 = ((rAB-2.0DB)(rAB-2.0DB)) + (aa);
	6582	+ double temp9 = (rABrAB) + (aa);
6571	6583	value = (rAB+2.0DA-2.0DB)*pow(temp1,-1.5)/16.0
6572	6584	+(rAB+2.0DA+2.0DB)*pow(temp2,-1.5)/16.0
6573	6585	+(rAB-2.0DA-2.0DB)*pow(temp3,-1.5)/16.0

		@@ -6581,14 +6593,14 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6581	6593	}
6582	6594	// Eq. (65) in [DT_1977]
6583	6595	else if(multipoleA == Qxz && multipoleB == Qxz){
6584		- double temp1 = pow(rAB+DA-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6585		- double temp2 = pow(rAB+DA-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
6586		- double temp3 = pow(rAB+DA+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6587		- double temp4 = pow(rAB+DA+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
6588		- double temp5 = pow(rAB-DA-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6589		- double temp6 = pow(rAB-DA-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
6590		- double temp7 = pow(rAB-DA+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
6591		- double temp8 = pow(rAB-DA+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
	6596	+ double temp1 = ((rAB+DA-DB)(rAB+DA-DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6597	+ double temp2 = ((rAB+DA-DB)(rAB+DA-DB)) + ((DA+DB)(DA+DB)) + (a*a);
	6598	+ double temp3 = ((rAB+DA+DB)(rAB+DA+DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6599	+ double temp4 = ((rAB+DA+DB)(rAB+DA+DB)) + ((DA+DB)(DA+DB)) + (a*a);
	6600	+ double temp5 = ((rAB-DA-DB)(rAB-DA-DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6601	+ double temp6 = ((rAB-DA-DB)(rAB-DA-DB)) + ((DA+DB)(DA+DB)) + (a*a);
	6602	+ double temp7 = ((rAB-DA+DB)(rAB-DA+DB)) + ((DA-DB)(DA-DB)) + (a*a);
	6603	+ double temp8 = ((rAB-DA+DB)(rAB-DA+DB)) + ((DA+DB)(DA+DB)) + (a*a);
6592	6604	value = (rAB+DA-DB)*pow(temp1,-1.5)/8.0
6593	6605	-(rAB+DA-DB)*pow(temp2,-1.5)/8.0
6594	6606	-(rAB+DA+DB)*pow(temp3,-1.5)/8.0

		@@ -6604,9 +6616,9 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction1stDerivative(const Atom& atomA
6604	6616	}
6605	6617	// Eq. (66) in [DT_1977]
6606	6618	else if(multipoleA == Qxy && multipoleB == Qxy){
6607		- double temp1 = pow(rAB,2.0) + 2.0*pow(DA-DB,2.0) + pow(a,2.0);
6608		- double temp2 = pow(rAB,2.0) + 2.0*pow(DA+DB,2.0) + pow(a,2.0);
6609		- double temp3 = pow(rAB,2.0) + 2.0pow(DA,2.0) + 2.0pow(DB,2.0) + pow(a,2.0);
	6619	+ double temp1 = (rABrAB) + 2.0((DA-DB)(DA-DB)) + (aa);
	6620	+ double temp2 = (rABrAB) + 2.0((DA+DB)(DA+DB)) + (aa);
	6621	+ double temp3 = (rABrAB) + 2.0(DADA) + 2.0(DBDB) + (aa);
6610	6622	value = (rAB)*pow(temp1,-1.5)/4.0
6611	6623	+(rAB)*pow(temp2,-1.5)/4.0
6612	6624	-(rAB)*pow(temp3,-1.5)/2.0;

		@@ -6640,74 +6652,75 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6640	6652	// Eq. (52) in [DT_1977]
6641	6653	if(multipoleA == sQ && multipoleB == sQ){
6642	6654	double c1 = 1.0;
6643		- double f1 = pow(rAB,2.0);
6644		- double a1 = pow(a,2.0);
	6655	+ double f1 = (rAB*rAB);
	6656	+ double a1 = (a*a);
6645	6657	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6646	6658	}
6647	6659	// Eq. (53) in [DT_1977]
6648	6660	else if(multipoleA == sQ && multipoleB == muz){
6649	6661	double c1 = 0.5;
6650	6662	double c2 = -0.5;
6651		- double f1 = pow(rAB+DB,2.0);
6652		- double f2 = pow(rAB-DB,2.0);
6653		- double a1 = pow(a,2.0);
6654		- double a2 = pow(a,2.0);
	6663	+ double f1 = ((rAB+DB)*(rAB+DB));
	6664	+ double f2 = ((rAB-DB)*(rAB-DB));
	6665	+ double a1 = (a*a);
	6666	+ double a2 = (a*a);
6655	6667	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6656	6668	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6657	6669	}
6658	6670	else if(multipoleA == muz && multipoleB == sQ){
6659	6671	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6660		- value *= pow(-1.0,1.0);
	6672	+ //value *= pow(-1.0,1.0);
	6673	+ value *= -1.0;
6661	6674	}
6662	6675	// Eq. (54) in [DT_1977]
6663	6676	else if(multipoleA == sQ && multipoleB == Qxx){
6664	6677	double c1 = 0.5;
6665	6678	double c2 = -0.5;
6666		- double f1 = pow(rAB,2.0);
6667		- double f2 = pow(rAB,2.0);
6668		- double a1 = pow(2.0*DB,2.0) + pow(a,2.0);
6669		- double a2 = pow(a,2.0);
	6679	+ double f1 = (rAB*rAB);
	6680	+ double f2 = (rAB*rAB);
	6681	+ double a1 = (4.0DBDB) + (a*a);
	6682	+ double a2 = (a*a);
6670	6683	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6671	6684	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6672	6685	}
6673	6686	else if(multipoleA == Qxx && multipoleB == sQ){
6674	6687	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6675		- value *= pow(-1.0,2.0);
	6688	+ //value *= pow(-1.0,2.0);
6676	6689	}
6677	6690	else if(multipoleA == sQ && multipoleB == Qyy){
6678	6691	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomA, atomB, multipoleA, Qxx, rAB);
6679	6692	}
6680	6693	else if(multipoleA == Qyy && multipoleB == sQ){
6681	6694	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6682		- value *= pow(-1.0,2.0);
	6695	+ //value *= pow(-1.0,2.0);
6683	6696	}
6684	6697	// Eq. (55) in [DT_1977]
6685	6698	else if(multipoleA == sQ && multipoleB == Qzz){
6686	6699	double c1 = 0.25;
6687	6700	double c2 = -0.50;
6688	6701	double c3 = 0.25;
6689		- double f1 = pow(rAB+2.0*DB,2.0);
6690		- double f2 = pow(rAB,2.0);
6691		- double f3 = pow(rAB-2.0*DB,2.0);
6692		- double a1 = pow(a,2.0);
6693		- double a2 = pow(a,2.0);
6694		- double a3 = pow(a,2.0);
	6702	+ double f1 = ((rAB+2.0DB)(rAB+2.0*DB));
	6703	+ double f2 = (rAB*rAB);
	6704	+ double f3 = ((rAB-2.0DB)(rAB-2.0*DB));
	6705	+ double a1 = (a*a);
	6706	+ double a2 = (a*a);
	6707	+ double a3 = (a*a);
6695	6708	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6696	6709	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6697	6710	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));
6698	6711	}
6699	6712	else if(multipoleA == Qzz && multipoleB == sQ){
6700	6713	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6701		- value *= pow(-1.0,2.0);
	6714	+ //value *= pow(-1.0,2.0);
6702	6715	}
6703	6716	// Eq. (56) in [DT_1977]
6704	6717	else if(multipoleA == mux && multipoleB == mux){
6705	6718	double c1 = 0.50;
6706	6719	double c2 = -0.50;
6707		- double f1 = pow(rAB,2.0);
6708		- double f2 = pow(rAB,2.0);
6709		- double a1 = pow(DA-DB,2.0) + pow(a,2.0);
6710		- double a2 = pow(DA+DB,2.0) + pow(a,2.0);
	6720	+ double f1 = (rAB*rAB);
	6721	+ double f2 = (rAB*rAB);
	6722	+ double a1 = ((DA-DB)(DA-DB)) + (aa);
	6723	+ double a2 = ((DA+DB)(DA+DB)) + (aa);
6711	6724	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6712	6725	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6713	6726	}

		@@ -6720,14 +6733,14 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6720	6733	double c2 = -0.25;
6721	6734	double c3 = -0.25;
6722	6735	double c4 = 0.25;
6723		- double f1 = pow(rAB+DA-DB,2.0);
6724		- double f2 = pow(rAB+DA+DB,2.0);
6725		- double f3 = pow(rAB-DA-DB,2.0);
6726		- double f4 = pow(rAB-DA+DB,2.0);
6727		- double a1 = pow(a,2.0);
6728		- double a2 = pow(a,2.0);
6729		- double a3 = pow(a,2.0);
6730		- double a4 = pow(a,2.0);
	6736	+ double f1 = ((rAB+DA-DB)*(rAB+DA-DB));
	6737	+ double f2 = ((rAB+DA+DB)*(rAB+DA+DB));
	6738	+ double f3 = ((rAB-DA-DB)*(rAB-DA-DB));
	6739	+ double f4 = ((rAB-DA+DB)*(rAB-DA+DB));
	6740	+ double a1 = (a*a);
	6741	+ double a2 = (a*a);
	6742	+ double a3 = (a*a);
	6743	+ double a4 = (a*a);
6731	6744	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6732	6745	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6733	6746	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

		@@ -6739,14 +6752,14 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6739	6752	double c2 = 0.25;
6740	6753	double c3 = 0.25;
6741	6754	double c4 = -0.25;
6742		- double f1 = pow(rAB-DB,2.0);
6743		- double f2 = pow(rAB-DB,2.0);
6744		- double f3 = pow(rAB+DB,2.0);
6745		- double f4 = pow(rAB+DB,2.0);
6746		- double a1 = pow(DA-DB,2.0) + pow(a,2.0);
6747		- double a2 = pow(DA+DB,2.0) + pow(a,2.0);
6748		- double a3 = pow(DA-DB,2.0) + pow(a,2.0);
6749		- double a4 = pow(DA+DB,2.0) + pow(a,2.0);
	6755	+ double f1 = ((rAB-DB)*(rAB-DB));
	6756	+ double f2 = ((rAB-DB)*(rAB-DB));
	6757	+ double f3 = ((rAB+DB)*(rAB+DB));
	6758	+ double f4 = ((rAB+DB)*(rAB+DB));
	6759	+ double a1 = ((DA-DB)(DA-DB)) + (aa);
	6760	+ double a2 = ((DA+DB)(DA+DB)) + (aa);
	6761	+ double a3 = ((DA-DB)(DA-DB)) + (aa);
	6762	+ double a4 = ((DA+DB)(DA+DB)) + (aa);
6750	6763	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6751	6764	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6752	6765	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

		@@ -6754,14 +6767,16 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6754	6767	}
6755	6768	else if(multipoleA == Qxz && multipoleB == mux){
6756	6769	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6757		- value *= pow(-1.0,3.0);
	6770	+ //value *= pow(-1.0,3.0);
	6771	+ value *= -1.0;
6758	6772	}
6759	6773	else if(multipoleA == muy && multipoleB == Qyz){
6760	6774	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomA, atomB, mux, Qxz, rAB);
6761	6775	}
6762	6776	else if(multipoleA == Qyz && multipoleB == muy){
6763	6777	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6764		- value *= pow(-1.0,3.0);
	6778	+ //value *= pow(-1.0,3.0);
	6779	+ value *= -1.0;
6765	6780	}
6766	6781	// Eq. (59) in [DT_1977]
6767	6782	else if(multipoleA == muz && multipoleB == Qxx){

		@@ -6769,14 +6784,14 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6769	6784	double c2 = 0.25;
6770	6785	double c3 = 0.25;
6771	6786	double c4 = -0.25;
6772		- double f1 = pow(rAB+DA,2.0);
6773		- double f2 = pow(rAB-DA,2.0);
6774		- double f3 = pow(rAB+DA,2.0);
6775		- double f4 = pow(rAB-DA,2.0);
6776		- double a1 = pow(2.0*DB,2.0) + pow(a,2.0);
6777		- double a2 = pow(2.0*DB,2.0) + pow(a,2.0);
6778		- double a3 = pow(a,2.0);
6779		- double a4 = pow(a,2.0);
	6787	+ double f1 = ((rAB+DA)*(rAB+DA));
	6788	+ double f2 = ((rAB-DA)*(rAB-DA));
	6789	+ double f3 = ((rAB+DA)*(rAB+DA));
	6790	+ double f4 = ((rAB-DA)*(rAB-DA));
	6791	+ double a1 = (4.0DBDB) + (a*a);
	6792	+ double a2 = (4.0DBDB) + (a*a);
	6793	+ double a3 = (a*a);
	6794	+ double a4 = (a*a);
6780	6795	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6781	6796	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6782	6797	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

		@@ -6784,14 +6799,16 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6784	6799	}
6785	6800	else if(multipoleA == Qxx && multipoleB == muz){
6786	6801	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6787		- value *= pow(-1.0,3.0);
	6802	+ //value *= pow(-1.0,3.0);
	6803	+ value *= -1.0;
6788	6804	}
6789	6805	else if(multipoleA == muz && multipoleB == Qyy){
6790	6806	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomA, atomB, muz, Qxx, rAB);
6791	6807	}
6792	6808	else if(multipoleA == Qyy && multipoleB == muz){
6793	6809	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6794		- value *= pow(-1.0,3.0);
	6810	+ //value *= pow(-1.0,3.0);
	6811	+ value *= -1.0;
6795	6812	}
6796	6813	// Eq. (60) in [DT_1977]
6797	6814	else if(multipoleA == muz && multipoleB == Qzz){

		@@ -6801,18 +6818,18 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6801	6818	double c4 = 0.125;
6802	6819	double c5 = 0.25;
6803	6820	double c6 = -0.25;
6804		- double f1 = pow(rAB+DA-2.0*DB,2.0);
6805		- double f2 = pow(rAB-DA-2.0*DB,2.0);
6806		- double f3 = pow(rAB+DA+2.0*DB,2.0);
6807		- double f4 = pow(rAB-DA+2.0*DB,2.0);
6808		- double f5 = pow(rAB+DA,2.0);
6809		- double f6 = pow(rAB-DA,2.0);
6810		- double a1 = pow(a,2.0);
6811		- double a2 = pow(a,2.0);
6812		- double a3 = pow(a,2.0);
6813		- double a4 = pow(a,2.0);
6814		- double a5 = pow(a,2.0);
6815		- double a6 = pow(a,2.0);
	6821	+ double f1 = ((rAB+DA-2.0DB)(rAB+DA-2.0*DB));
	6822	+ double f2 = ((rAB-DA-2.0DB)(rAB-DA-2.0*DB));
	6823	+ double f3 = ((rAB+DA+2.0DB)(rAB+DA+2.0*DB));
	6824	+ double f4 = ((rAB-DA+2.0DB)(rAB-DA+2.0*DB));
	6825	+ double f5 = ((rAB+DA)*(rAB+DA));
	6826	+ double f6 = ((rAB-DA)*(rAB-DA));
	6827	+ double a1 = (a*a);
	6828	+ double a2 = (a*a);
	6829	+ double a3 = (a*a);
	6830	+ double a4 = (a*a);
	6831	+ double a5 = (a*a);
	6832	+ double a6 = (a*a);
6816	6833	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6817	6834	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6818	6835	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

		@@ -6822,7 +6839,8 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6822	6839	}
6823	6840	else if(multipoleA == Qzz && multipoleB == muz){
6824	6841	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6825		- value *= pow(-1.0,3.0);
	6842	+ //value *= pow(-1.0,3.0);
	6843	+ value *= -1.0;
6826	6844	}
6827	6845	// Eq. (61) in [DT_1977]
6828	6846	else if(multipoleA == Qxx && multipoleB == Qxx){

		@@ -6831,16 +6849,16 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6831	6849	double c3 = -0.25;
6832	6850	double c4 = -0.25;
6833	6851	double c5 = 0.25;
6834		- double f1 = pow(rAB,2.0);
6835		- double f2 = pow(rAB,2.0);
6836		- double f3 = pow(rAB,2.0);
6837		- double f4 = pow(rAB,2.0);
6838		- double f5 = pow(rAB,2.0);
6839		- double a1 = 4.0*pow(DA-DB,2.0) + pow(a,2.0);
6840		- double a2 = 4.0*pow(DA+DB,2.0) + pow(a,2.0);
6841		- double a3 = pow(2.0*DA,2.0) + pow(a,2.0);
6842		- double a4 = pow(2.0*DB,2.0) + pow(a,2.0);
6843		- double a5 = pow(a,2.0);
	6852	+ double f1 = (rAB*rAB);
	6853	+ double f2 = (rAB*rAB);
	6854	+ double f3 = (rAB*rAB);
	6855	+ double f4 = (rAB*rAB);
	6856	+ double f5 = (rAB*rAB);
	6857	+ double a1 = 4.0((DA-DB)(DA-DB)) + (a*a);
	6858	+ double a2 = 4.0((DA+DB)(DA+DB)) + (a*a);
	6859	+ double a3 = (4.0DADA) + (a*a);
	6860	+ double a4 = (4.0DBDB) + (a*a);
	6861	+ double a5 = (a*a);
6844	6862	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6845	6863	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6846	6864	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

		@@ -6856,14 +6874,14 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6856	6874	double c2 = -0.25;
6857	6875	double c3 = -0.25;
6858	6876	double c4 = 0.25;
6859		- double f1 = pow(rAB,2.0);
6860		- double f2 = pow(rAB,2.0);
6861		- double f3 = pow(rAB,2.0);
6862		- double f4 = pow(rAB,2.0);
6863		- double a1 = pow(2.0DA,2.0) + pow(2.0DB,2.0) + pow(a,2.0);
6864		- double a2 = pow(2.0*DA,2.0) + pow(a,2.0);
6865		- double a3 = pow(2.0*DB,2.0) + pow(a,2.0);
6866		- double a4 = pow(a,2.0);
	6877	+ double f1 = (rAB*rAB);
	6878	+ double f2 = (rAB*rAB);
	6879	+ double f3 = (rAB*rAB);
	6880	+ double f4 = (rAB*rAB);
	6881	+ double a1 = (4.0DADA) + (4.0DBDB) + (a*a);
	6882	+ double a2 = (4.0DADA) + (a*a);
	6883	+ double a3 = (4.0DBDB) + (a*a);
	6884	+ double a4 = (a*a);
6867	6885	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6868	6886	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6869	6887	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

		@@ -6871,7 +6889,7 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6871	6889	}
6872	6890	else if(multipoleA == Qyy && multipoleB == Qxx){
6873	6891	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6874		- value *= pow(-1.0,4.0);
	6892	+ //value *= pow(-1.0,4.0);
6875	6893	}
6876	6894	// Eq. (63) in [DT_1977]
6877	6895	else if(multipoleA == Qxx && multipoleB == Qzz){

		@@ -6881,18 +6899,18 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6881	6899	double c4 = -0.125;
6882	6900	double c5 = -0.25;
6883	6901	double c6 = 0.25;
6884		- double f1 = pow(rAB-2.0*DB,2.0);
6885		- double f2 = pow(rAB+2.0*DB,2.0);
6886		- double f3 = pow(rAB-2.0*DB,2.0);
6887		- double f4 = pow(rAB+2.0*DB,2.0);
6888		- double f5 = pow(rAB ,2.0);
6889		- double f6 = pow(rAB ,2.0);
6890		- double a1 = pow(2.0*DA,2.0) + pow(a,2.0);
6891		- double a2 = pow(2.0*DA,2.0) + pow(a,2.0);
6892		- double a3 = pow(a,2.0);
6893		- double a4 = pow(a,2.0);
6894		- double a5 = pow(2.0*DA,2.0) + pow(a,2.0);
6895		- double a6 = pow(a,2.0);
	6902	+ double f1 = ((rAB-2.0DB)(rAB-2.0*DB));
	6903	+ double f2 = ((rAB+2.0DB)(rAB+2.0*DB));
	6904	+ double f3 = ((rAB-2.0DB)(rAB-2.0*DB));
	6905	+ double f4 = ((rAB+2.0DB)(rAB+2.0*DB));
	6906	+ double f5 = rAB*rAB;
	6907	+ double f6 = rAB*rAB;
	6908	+ double a1 = (4.0DADA) + (a*a);
	6909	+ double a2 = (4.0DADA) + (a*a);
	6910	+ double a3 = (a*a);
	6911	+ double a4 = (a*a);
	6912	+ double a5 = (4.0DADA) + (a*a);
	6913	+ double a6 = (a*a);
6896	6914	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6897	6915	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6898	6916	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

		@@ -6902,14 +6920,14 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6902	6920	}
6903	6921	else if(multipoleA == Qzz && multipoleB == Qxx){
6904	6922	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6905		- value *= pow(-1.0,4.0);
	6923	+ //value *= pow(-1.0,4.0);
6906	6924	}
6907	6925	else if(multipoleA == Qyy && multipoleB == Qzz){
6908	6926	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomA, atomB, Qxx, multipoleB, rAB);
6909	6927	}
6910	6928	else if(multipoleA == Qzz && multipoleB == Qyy){
6911	6929	value = this->GetSemiEmpiricalMultipoleInteraction2ndDerivative(atomB, atomA, multipoleB, multipoleA, rAB);
6912		- value *= pow(-1.0,4.0);
	6930	+ //value *= pow(-1.0,4.0);
6913	6931	}
6914	6932	// Eq. (64) in [DT_1977]
6915	6933	else if(multipoleA == Qzz && multipoleB == Qzz){

		@@ -6922,24 +6940,24 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6922	6940	double c7 = -0.125;
6923	6941	double c8 = -0.125;
6924	6942	double c9 = 0.25;
6925		- double f1 = pow(rAB+2.0DA-2.0DB,2.0);
6926		- double f2 = pow(rAB+2.0DA+2.0DB,2.0);
6927		- double f3 = pow(rAB-2.0DA-2.0DB,2.0);
6928		- double f4 = pow(rAB-2.0DA+2.0DB,2.0);
6929		- double f5 = pow(rAB+2.0*DA ,2.0);
6930		- double f6 = pow(rAB-2.0*DA ,2.0);
6931		- double f7 = pow(rAB+2.0*DB ,2.0);
6932		- double f8 = pow(rAB-2.0*DB ,2.0);
6933		- double f9 = pow(rAB ,2.0);
6934		- double a1 = pow(a,2.0);
6935		- double a2 = pow(a,2.0);
6936		- double a3 = pow(a,2.0);
6937		- double a4 = pow(a,2.0);
6938		- double a5 = pow(a,2.0);
6939		- double a6 = pow(a,2.0);
6940		- double a7 = pow(a,2.0);
6941		- double a8 = pow(a,2.0);
6942		- double a9 = pow(a,2.0);
	6943	+ double f1 = ((rAB+2.0DA-2.0DB)(rAB+2.0DA-2.0*DB));
	6944	+ double f2 = ((rAB+2.0DA+2.0DB)(rAB+2.0DA+2.0*DB));
	6945	+ double f3 = ((rAB-2.0DA-2.0DB)(rAB-2.0DA-2.0*DB));
	6946	+ double f4 = ((rAB-2.0DA+2.0DB)(rAB-2.0DA+2.0*DB));
	6947	+ double f5 = ((rAB+2.0DA)(rAB+2.0*DA));
	6948	+ double f6 = ((rAB-2.0DA)(rAB-2.0*DA));
	6949	+ double f7 = ((rAB+2.0DB)(rAB+2.0*DB));
	6950	+ double f8 = ((rAB-2.0DB)(rAB-2.0*DB));
	6951	+ double f9 = (rAB*rAB);
	6952	+ double a1 = (a*a);
	6953	+ double a2 = (a*a);
	6954	+ double a3 = (a*a);
	6955	+ double a4 = (a*a);
	6956	+ double a5 = (a*a);
	6957	+ double a6 = (a*a);
	6958	+ double a7 = (a*a);
	6959	+ double a8 = (a*a);
	6960	+ double a9 = (a*a);
6943	6961	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6944	6962	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6945	6963	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

		@@ -6960,22 +6978,22 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6960	6978	double c6 = 0.125;
6961	6979	double c7 = 0.125;
6962	6980	double c8 = -0.125;
6963		- double f1 = pow(rAB+DA-DB,2.0);
6964		- double f2 = pow(rAB+DA-DB,2.0);
6965		- double f3 = pow(rAB+DA+DB,2.0);
6966		- double f4 = pow(rAB+DA+DB,2.0);
6967		- double f5 = pow(rAB-DA-DB,2.0);
6968		- double f6 = pow(rAB-DA-DB,2.0);
6969		- double f7 = pow(rAB-DA+DB,2.0);
6970		- double f8 = pow(rAB-DA+DB,2.0);
6971		- double a1 = pow(DA-DB,2.0) + pow(a,2.0);
6972		- double a2 = pow(DA+DB,2.0) + pow(a,2.0);
6973		- double a3 = pow(DA-DB,2.0) + pow(a,2.0);
6974		- double a4 = pow(DA+DB,2.0) + pow(a,2.0);
6975		- double a5 = pow(DA-DB,2.0) + pow(a,2.0);
6976		- double a6 = pow(DA+DB,2.0) + pow(a,2.0);
6977		- double a7 = pow(DA-DB,2.0) + pow(a,2.0);
6978		- double a8 = pow(DA+DB,2.0) + pow(a,2.0);
	6981	+ double f1 = ((rAB+DA-DB)*(rAB+DA-DB));
	6982	+ double f2 = ((rAB+DA-DB)*(rAB+DA-DB));
	6983	+ double f3 = ((rAB+DA+DB)*(rAB+DA+DB));
	6984	+ double f4 = ((rAB+DA+DB)*(rAB+DA+DB));
	6985	+ double f5 = ((rAB-DA-DB)*(rAB-DA-DB));
	6986	+ double f6 = ((rAB-DA-DB)*(rAB-DA-DB));
	6987	+ double f7 = ((rAB-DA+DB)*(rAB-DA+DB));
	6988	+ double f8 = ((rAB-DA+DB)*(rAB-DA+DB));
	6989	+ double a1 = ((DA-DB)(DA-DB)) + (aa);
	6990	+ double a2 = ((DA+DB)(DA+DB)) + (aa);
	6991	+ double a3 = ((DA-DB)(DA-DB)) + (aa);
	6992	+ double a4 = ((DA+DB)(DA+DB)) + (aa);
	6993	+ double a5 = ((DA-DB)(DA-DB)) + (aa);
	6994	+ double a6 = ((DA+DB)(DA+DB)) + (aa);
	6995	+ double a7 = ((DA-DB)(DA-DB)) + (aa);
	6996	+ double a8 = ((DA+DB)(DA+DB)) + (aa);
6979	6997	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
6980	6998	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
6981	6999	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

		@@ -6993,12 +7011,12 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction2ndDerivative(const Atom& atomA
6993	7011	double c1 = 0.25;
6994	7012	double c2 = 0.25;
6995	7013	double c3 = -0.50;
6996		- double f1 = pow(rAB,2.0);
6997		- double f2 = pow(rAB,2.0);
6998		- double f3 = pow(rAB,2.0);
6999		- double a1 = 2.0*pow(DA-DB,2.0) + pow(a,2.0);
7000		- double a2 = 2.0*pow(DA+DB,2.0) + pow(a,2.0);
7001		- double a3 = 2.0pow(DA,2.0) + 2.0pow(DB,2.0) + pow(a,2.0);
	7014	+ double f1 = (rAB*rAB);
	7015	+ double f2 = (rAB*rAB);
	7016	+ double f3 = (rAB*rAB);
	7017	+ double a1 = 2.0((DA-DB)(DA-DB)) + (a*a);
	7018	+ double a2 = 2.0((DA+DB)(DA+DB)) + (a*a);
	7019	+ double a3 = 2.0(DADA) + 2.0(DBDB) + (a*a);
7002	7020	value = c1(3.0f1*pow(f1+a1,-2.5) - pow(f1+a1,-1.5));
7003	7021	value += c2(3.0f2*pow(f2+a2,-2.5) - pow(f2+a2,-1.5));
7004	7022	value += c3(3.0f3*pow(f3+a3,-2.5) - pow(f3+a3,-1.5));

MolDS Fork