OSDN -- 오픈 소스 소프트웨어 개발 및 다운로드
Translation Status of 한국말 translation status for ko
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

Commit

Tags
No Tags

Frequently used words (click to add to your profile)

javac++androidlinuxc#windowsobjective-ccocoa誰得qtpythonphprubygameguibathyscaphec計画中(planning stage)翻訳omegatframeworktwitterdomtestvb.netdirectxゲームエンジンbtronarduinopreviewer

Commit MetaInfo

Revision527797991a235e911aeebed374c4c729fc0900ae (tree)
Time2011-10-20 01:14:16
AuthorMikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Log Message

Mndo::GetSemiEmpiricalMultipoleInteractionFirstDerivative and Mndo::GetNddoRepulsionIntegralFirstDerivative are implemented. #26395

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

Change Summary

Incremental Difference

--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -68,9 +68,11 @@ protected:
6868 virtual void CalcForce(int electronicStateIndex);
6969 private:
7070 string errorMessageGetSemiEmpiricalMultipoleInteractionBadMultipoles;
71+ string errorMessageGetSemiEmpiricalMultipoleInteractionFirstDeriBadMultipoles;
7172 string errorMessageMultipoleA;
7273 string errorMessageMultipoleB;
7374 string errorMessageGetNddoRepulsionIntegral;
75+ string errorMessageGetNddoRepulsionIntegralFirstDerivative;
7476 string errorMessageCalcTwoElecTwoCoreDiatomicNullMatrix;
7577 string errorMessageCalcTwoElecTwoCoreNullMatrix;
7678 string errorMessageCalcTwoElecTwoCoreDiatomicSameAtoms;
@@ -84,6 +86,10 @@ private:
8486 void RotateTwoElecTwoCoreDiatomicToSpaceFramegc(double**** matrix, double** rotatingMatrix);
8587 double GetNddoRepulsionIntegral(Atom* atomA, OrbitalType mu, OrbitalType nu,
8688 Atom* atomB, OrbitalType lambda, OrbitalType sigma);
89+ double GetNddoRepulsionIntegralFirstDerivative(
90+ Atom* atomA, OrbitalType mu, OrbitalType nu,
91+ Atom* atomB, OrbitalType lambda, OrbitalType sigma,
92+ CartesianType axisA);
8793 double GetSemiEmpiricalMultipoleInteraction(MultipoleType multipoleA,
8894 MultipoleType multipoleB,
8995 double rhoA,
@@ -91,6 +97,14 @@ private:
9197 double DA,
9298 double DB,
9399 double Rab);
100+ double GetSemiEmpiricalMultipoleInteractionFirstDerivative(
101+ MultipoleType multipoleA,
102+ MultipoleType multipoleB,
103+ double rhoA,
104+ double rhoB,
105+ double DA,
106+ double DB,
107+ double Rab);
94108 };
95109
96110 Mndo::Mndo() : MolDS_zindo::ZindoS(){
@@ -138,9 +152,14 @@ void Mndo::SetMessages(){
138152 this->errorMessageDavidsonNotConverged = "Error in mndo::Mndo::DoesCISDavidson: Davidson did not met convergence criterion. \n";
139153 this->errorMessageGetSemiEmpiricalMultipoleInteractionBadMultipoles
140154 = "Error in mndo:: Mndo::GetSemiEmpiricalMultipoleInteraction: Bad multipole combintaion is set\n";
155+ this->errorMessageGetSemiEmpiricalMultipoleInteractionFirstDeriBadMultipoles
156+ = "Error in mndo:: Mndo::GetSemiEmpiricalMultipoleInteractionFirstDerivative: Bad multipole combintaion is set\n";
141157 this->errorMessageMultipoleA = "Multipole A is: ";
142158 this->errorMessageMultipoleB = "Multipole B is: ";
143- this->errorMessageGetNddoRepulsionIntegral = "Error in mndo::Mndo::GetNddoRepulsionIntegral: Bad orbital is set.\n";
159+ this->errorMessageGetNddoRepulsionIntegral
160+ = "Error in mndo::Mndo::GetNddoRepulsionIntegral: Bad orbital is set.\n";
161+ this->errorMessageGetNddoRepulsionIntegralFirstDerivative
162+ = "Error in mndo::Mndo::GetNddoRepulsionIntegralFirstDerivative: Bad orbital is set.\n";
144163 this->errorMessageCalcTwoElecTwoCoreNullMatrix
145164 = "Error in mndo::Mndo::CalcTwoElecTwoCore: The two elec two core matrix is NULL.\n";
146165 this->errorMessageCalcTwoElecTwoCoreDiatomicNullMatrix
@@ -1397,7 +1416,6 @@ double Mndo::GetNddoRepulsionIntegral(Atom* atomA, OrbitalType mu, OrbitalType n
13971416 lambda == dxy || lambda == dyz || lambda == dzz || lambda == dzx || lambda == dxxyy ||
13981417 sigma == dxy || sigma == dyz || sigma == dzz || sigma == dzx || sigma == dxxyy){
13991418
1400- // ToDo: error log.
14011419 stringstream ss;
14021420 ss << this->errorMessageGetNddoRepulsionIntegral;
14031421 ss << this->errorMessageAtomA << AtomTypeStr(atomA->GetAtomType()) << endl;
@@ -1414,98 +1432,881 @@ double Mndo::GetNddoRepulsionIntegral(Atom* atomA, OrbitalType mu, OrbitalType n
14141432 return value;
14151433 }
14161434
1435+// First derivative of NDDO repulsion integral.
1436+// This derivation is related to the coordinate of atomA
14171437 // See Apendix in [DT_1977]
1418-double Mndo::GetSemiEmpiricalMultipoleInteraction(MultipoleType multipoleA,
1419- MultipoleType multipoleB,
1420- double rhoA,
1421- double rhoB,
1422- double DA,
1423- double DB,
1424- double Rab){
1438+// Orbital mu and nu belong atom A,
1439+// orbital lambda and sigma belong atomB.
1440+double Mndo::GetNddoRepulsionIntegralFirstDerivative(
1441+ Atom* atomA, OrbitalType mu, OrbitalType nu,
1442+ Atom* atomB, OrbitalType lambda, OrbitalType sigma,
1443+ CartesianType axisA){
14251444 double value = 0.0;
1426- double a = rhoA + rhoB;
1427-
1428- if(multipoleA == sQ && multipoleB == sQ){
1429- value = pow(pow(Rab,2.0) + pow(a,2.0), -0.5);
1430- }
1431- else if(multipoleA == sQ && multipoleB == muz){
1432- double temp1 = pow(Rab+DB,2.0) + pow(a,2.0);
1433- double temp2 = pow(Rab-DB,2.0) + pow(a,2.0);
1434- value = pow(temp1,-0.5)/2.0 - pow(temp2,-0.5)/2.0;
1435- }
1436- else if(multipoleA == muz && multipoleB == sQ){
1437- value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
1438- rhoB, rhoA, DB, DA, Rab);
1439- value *= pow(-1.0,1.0);
1440- }
1441- else if(multipoleA == sQ && multipoleB == Qxx){
1442- double temp1 = pow(Rab,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
1443- double temp2 = pow(Rab,2.0) + pow(a,2.0);
1444- value = pow(temp1,-0.5)/2.0 - pow(temp2,-0.5)/2.0;
1445+ double DA=0.0;
1446+ double DB=0.0;
1447+ double rhoA = 0.0;
1448+ double rhoB = 0.0;
1449+ double Rab = this->molecule->GetDistanceAtoms(atomA, atomB);
1450+ double dRabDa = (atomA->GetXyz()[axisA] - atomB->GetXyz()[axisA])/Rab;
1451+ int lA = 0;
1452+ int lB = 0;
1453+ // (28) in [DT_1977]
1454+ if(mu == s && nu == s && lambda == s && sigma == s){
1455+ DA = atomA->GetMndoDerivedParameterD(0);
1456+ DB = atomB->GetMndoDerivedParameterD(0);
1457+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1458+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1459+ value = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1460+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1461+ value *= dRabDa;
14451462 }
1446- else if(multipoleA == Qxx && multipoleB == sQ){
1447- value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
1448- rhoB, rhoA, DB, DA, Rab);
1449- value *= pow(-1.0,2.0);
1463+ // (29) in [DT_1977]
1464+ else if(mu == s && nu == s && lambda == px && sigma == px){
1465+ DA = atomA->GetMndoDerivedParameterD(0);
1466+ DB = atomB->GetMndoDerivedParameterD(0);
1467+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1468+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1469+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1470+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1471+ DA = atomA->GetMndoDerivedParameterD(0);
1472+ DB = atomB->GetMndoDerivedParameterD(2);
1473+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1474+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1475+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1476+ sQ, Qxx, rhoA, rhoB, DA, DB, Rab);
1477+ value = temp1 + temp2;
1478+ value *= dRabDa;
14501479 }
1451- else if(multipoleA == sQ && multipoleB == Qyy){
1452- value = this->GetSemiEmpiricalMultipoleInteraction(multipoleA, Qxx,
1453- rhoA, rhoB, DA, DB, Rab);
1480+ else if(mu == s && nu == s && lambda == py && sigma == py){
1481+ DA = atomA->GetMndoDerivedParameterD(0);
1482+ DB = atomB->GetMndoDerivedParameterD(0);
1483+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1484+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1485+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1486+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1487+ DA = atomA->GetMndoDerivedParameterD(0);
1488+ DB = atomB->GetMndoDerivedParameterD(2);
1489+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1490+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1491+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1492+ sQ, Qyy, rhoA, rhoB, DA, DB, Rab);
1493+ value = temp1 + temp2;
1494+ value *= dRabDa;
14541495 }
1455- else if(multipoleA == Qyy && multipoleB == sQ){
1456- value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
1457- rhoB, rhoA, DB, DA, Rab);
1458- value *= pow(-1.0,2.0);
1496+ // (30) in [DT_1977]
1497+ else if(mu == s && nu == s && lambda == pz && sigma == pz){
1498+ DA = atomA->GetMndoDerivedParameterD(0);
1499+ DB = atomB->GetMndoDerivedParameterD(0);
1500+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1501+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1502+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1503+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1504+ DA = atomA->GetMndoDerivedParameterD(0);
1505+ DB = atomB->GetMndoDerivedParameterD(2);
1506+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1507+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1508+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1509+ sQ, Qzz, rhoA, rhoB, DA, DB, Rab);
1510+ value = temp1 + temp2;
1511+ value *= dRabDa;
14591512 }
1460- else if(multipoleA == sQ && multipoleB == Qzz){
1461- double temp1 = pow(Rab+2.0*DB,2.0) + pow(a,2.0);
1462- double temp2 = pow(Rab,2.0) + pow(a,2.0);
1463- double temp3 = pow(Rab-2.0*DB,2.0) + pow(a,2.0);
1464- value = pow(temp1,-0.5)/4.0 - pow(temp2,-0.5)/2.0 + pow(temp3,-0.5)/4.0;
1513+ // (31) in [DT_1977]
1514+ else if(mu == px && nu == px && lambda == s && sigma == s){
1515+ DA = atomA->GetMndoDerivedParameterD(0);
1516+ DB = atomB->GetMndoDerivedParameterD(0);
1517+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1518+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1519+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1520+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1521+ DA = atomA->GetMndoDerivedParameterD(2);
1522+ DB = atomB->GetMndoDerivedParameterD(0);
1523+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1524+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1525+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1526+ Qxx, sQ, rhoA, rhoB, DA, DB, Rab);
1527+ value = temp1 + temp2;
1528+ value *= dRabDa;
14651529 }
1466- else if(multipoleA == Qzz && multipoleB == sQ){
1467- value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
1468- rhoB, rhoA, DB, DA, Rab);
1469- value *= pow(-1.0,2.0);
1530+ else if(mu == py && nu == py && lambda == s && sigma == s){
1531+ DA = atomA->GetMndoDerivedParameterD(0);
1532+ DB = atomB->GetMndoDerivedParameterD(0);
1533+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1534+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1535+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1536+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1537+ DA = atomA->GetMndoDerivedParameterD(2);
1538+ DB = atomB->GetMndoDerivedParameterD(0);
1539+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1540+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1541+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1542+ Qyy, sQ, rhoA, rhoB, DA, DB, Rab);
1543+ value = temp1 + temp2;
1544+ value *= dRabDa;
14701545 }
1471- else if(multipoleA == mux && multipoleB == mux){
1472- double temp1 = pow(Rab,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
1473- double temp2 = pow(Rab,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
1474- value = pow(temp1,-0.5)/2.0 - pow(temp2,-0.5)/2.0;
1546+ // (32) in [DT_1977]
1547+ else if(mu == pz && nu == pz && lambda == s && sigma == s){
1548+ DA = atomA->GetMndoDerivedParameterD(0);
1549+ DB = atomB->GetMndoDerivedParameterD(0);
1550+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1551+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1552+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1553+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1554+ DA = atomA->GetMndoDerivedParameterD(2);
1555+ DB = atomB->GetMndoDerivedParameterD(0);
1556+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1557+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1558+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1559+ Qzz, sQ, rhoA, rhoB, DA, DB, Rab);
1560+ value = temp1 + temp2;
1561+ value *= dRabDa;
14751562 }
1476- else if(multipoleA == muy && multipoleB == muy){
1477- value = this->GetSemiEmpiricalMultipoleInteraction(mux, mux,
1478- rhoA, rhoB, DA, DB, Rab);
1563+ // (33) in [DT_1977]
1564+ else if(mu == px && nu == px && lambda == px && sigma == px){
1565+ DA = atomA->GetMndoDerivedParameterD(0);
1566+ DB = atomB->GetMndoDerivedParameterD(0);
1567+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1568+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1569+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1570+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1571+ DA = atomA->GetMndoDerivedParameterD(0);
1572+ DB = atomB->GetMndoDerivedParameterD(2);
1573+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1574+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1575+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1576+ sQ, Qxx, rhoA, rhoB, DA, DB, Rab);
1577+ DA = atomA->GetMndoDerivedParameterD(2);
1578+ DB = atomB->GetMndoDerivedParameterD(0);
1579+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1580+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1581+ double temp3 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1582+ Qxx, sQ, rhoA, rhoB, DA, DB, Rab);
1583+ DA = atomA->GetMndoDerivedParameterD(2);
1584+ DB = atomB->GetMndoDerivedParameterD(2);
1585+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1586+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1587+ double temp4 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1588+ Qxx, Qxx, rhoA, rhoB, DA, DB, Rab);
1589+ value = temp1 + temp2 + temp3 + temp4;
1590+ value *= dRabDa;
14791591 }
1480- else if(multipoleA == muz && multipoleB == muz){
1481- double temp1 = pow(Rab+DA-DB,2.0) + pow(a,2.0);
1482- double temp2 = pow(Rab+DA+DB,2.0) + pow(a,2.0);
1483- double temp3 = pow(Rab-DA-DB,2.0) + pow(a,2.0);
1484- double temp4 = pow(Rab-DA+DB,2.0) + pow(a,2.0);
1485- value = pow(temp1,-0.5)/4.0 - pow(temp2,-0.5)/4.0
1486- -pow(temp3,-0.5)/4.0 + pow(temp4,-0.5)/4.0;
1592+ else if(mu == py && nu == py && lambda == py && sigma == py){
1593+ DA = atomA->GetMndoDerivedParameterD(0);
1594+ DB = atomB->GetMndoDerivedParameterD(0);
1595+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1596+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1597+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1598+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1599+ DA = atomA->GetMndoDerivedParameterD(0);
1600+ DB = atomB->GetMndoDerivedParameterD(2);
1601+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1602+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1603+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1604+ sQ, Qyy, rhoA, rhoB, DA, DB, Rab);
1605+ DA = atomA->GetMndoDerivedParameterD(2);
1606+ DB = atomB->GetMndoDerivedParameterD(0);
1607+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1608+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1609+ double temp3 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1610+ Qyy, sQ, rhoA, rhoB, DA, DB, Rab);
1611+ DA = atomA->GetMndoDerivedParameterD(2);
1612+ DB = atomB->GetMndoDerivedParameterD(2);
1613+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1614+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1615+ double temp4 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1616+ Qyy, Qyy, rhoA, rhoB, DA, DB, Rab);
1617+ value = temp1 + temp2 + temp3 + temp4;
1618+ value *= dRabDa;
14871619 }
1488- else if(multipoleA == mux && multipoleB == Qxz){
1489- double temp1 = pow(Rab-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
1490- double temp2 = pow(Rab-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
1491- double temp3 = pow(Rab+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
1492- double temp4 = pow(Rab+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
1493- value =-pow(temp1,-0.5)/4.0 + pow(temp2,-0.5)/4.0
1494- +pow(temp3,-0.5)/4.0 - pow(temp4,-0.5)/4.0;
1620+ // (34) in [DT_1977]
1621+ else if(mu == px && nu == px && lambda == py && sigma == py){
1622+ DA = atomA->GetMndoDerivedParameterD(0);
1623+ DB = atomB->GetMndoDerivedParameterD(0);
1624+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1625+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1626+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1627+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1628+ DA = atomA->GetMndoDerivedParameterD(0);
1629+ DB = atomB->GetMndoDerivedParameterD(2);
1630+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1631+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1632+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1633+ sQ, Qyy, rhoA, rhoB, DA, DB, Rab);
1634+ DA = atomA->GetMndoDerivedParameterD(2);
1635+ DB = atomB->GetMndoDerivedParameterD(0);
1636+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1637+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1638+ double temp3 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1639+ Qxx, sQ, rhoA, rhoB, DA, DB, Rab);
1640+ DA = atomA->GetMndoDerivedParameterD(2);
1641+ DB = atomB->GetMndoDerivedParameterD(2);
1642+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1643+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1644+ double temp4 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1645+ Qxx, Qyy, rhoA, rhoB, DA, DB, Rab);
1646+ value = temp1 + temp2 + temp3 + temp4;
1647+ value *= dRabDa;
14951648 }
1496- else if(multipoleA == Qxz && multipoleB == mux){
1497- value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
1498- rhoB, rhoA, DB, DA, Rab);
1499- value *= pow(-1.0,3.0);
1649+ else if(mu == py && nu == py && lambda == px && sigma == px){
1650+ DA = atomA->GetMndoDerivedParameterD(0);
1651+ DB = atomB->GetMndoDerivedParameterD(0);
1652+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1653+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1654+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1655+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1656+ DA = atomA->GetMndoDerivedParameterD(0);
1657+ DB = atomB->GetMndoDerivedParameterD(2);
1658+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1659+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1660+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1661+ sQ, Qxx, rhoA, rhoB, DA, DB, Rab);
1662+ DA = atomA->GetMndoDerivedParameterD(2);
1663+ DB = atomB->GetMndoDerivedParameterD(0);
1664+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1665+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1666+ double temp3 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1667+ Qyy, sQ, rhoA, rhoB, DA, DB, Rab);
1668+ DA = atomA->GetMndoDerivedParameterD(2);
1669+ DB = atomB->GetMndoDerivedParameterD(2);
1670+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1671+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1672+ double temp4 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1673+ Qyy, Qxx, rhoA, rhoB, DA, DB, Rab);
1674+ value = temp1 + temp2 + temp3 + temp4;
1675+ value *= dRabDa;
15001676 }
1501- else if(multipoleA == muy && multipoleB == Qyz){
1502- value = this->GetSemiEmpiricalMultipoleInteraction(mux, Qxz,
1503- rhoA, rhoB, DA, DB, Rab);
1677+ // (35) in [DT_1977]
1678+ else if(mu == px && nu == px && lambda == pz && sigma == pz){
1679+ DA = atomA->GetMndoDerivedParameterD(0);
1680+ DB = atomB->GetMndoDerivedParameterD(0);
1681+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1682+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1683+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1684+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1685+ DA = atomA->GetMndoDerivedParameterD(0);
1686+ DB = atomB->GetMndoDerivedParameterD(2);
1687+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1688+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1689+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1690+ sQ, Qzz, rhoA, rhoB, DA, DB, Rab);
1691+ DA = atomA->GetMndoDerivedParameterD(2);
1692+ DB = atomB->GetMndoDerivedParameterD(0);
1693+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1694+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1695+ double temp3 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1696+ Qxx, sQ, rhoA, rhoB, DA, DB, Rab);
1697+ DA = atomA->GetMndoDerivedParameterD(2);
1698+ DB = atomB->GetMndoDerivedParameterD(2);
1699+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1700+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1701+ double temp4 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1702+ Qxx, Qzz, rhoA, rhoB, DA, DB, Rab);
1703+ value = temp1 + temp2 + temp3 + temp4;
1704+ value *= dRabDa;
15041705 }
1505- else if(multipoleA == Qyz && multipoleB == muy){
1506- value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
1507- rhoB, rhoA, DB, DA, Rab);
1508- value *= pow(-1.0,3.0);
1706+ else if(mu == py && nu == py && lambda == pz && sigma == pz){
1707+ DA = atomA->GetMndoDerivedParameterD(0);
1708+ DB = atomB->GetMndoDerivedParameterD(0);
1709+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1710+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1711+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1712+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1713+ DA = atomA->GetMndoDerivedParameterD(0);
1714+ DB = atomB->GetMndoDerivedParameterD(2);
1715+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1716+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1717+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1718+ sQ, Qzz, rhoA, rhoB, DA, DB, Rab);
1719+ DA = atomA->GetMndoDerivedParameterD(2);
1720+ DB = atomB->GetMndoDerivedParameterD(0);
1721+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1722+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1723+ double temp3 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1724+ Qyy, sQ, rhoA, rhoB, DA, DB, Rab);
1725+ DA = atomA->GetMndoDerivedParameterD(2);
1726+ DB = atomB->GetMndoDerivedParameterD(2);
1727+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1728+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1729+ double temp4 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1730+ Qyy, Qzz, rhoA, rhoB, DA, DB, Rab);
1731+ value = temp1 + temp2 + temp3 + temp4;
1732+ value *= dRabDa;
1733+ }
1734+ // (36) in [DT_1977]
1735+ else if(mu == pz && nu == pz && lambda == px && sigma == px){
1736+ DA = atomA->GetMndoDerivedParameterD(0);
1737+ DB = atomB->GetMndoDerivedParameterD(0);
1738+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1739+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1740+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1741+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1742+ DA = atomA->GetMndoDerivedParameterD(0);
1743+ DB = atomB->GetMndoDerivedParameterD(2);
1744+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1745+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1746+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1747+ sQ, Qxx, rhoA, rhoB, DA, DB, Rab);
1748+ DA = atomA->GetMndoDerivedParameterD(2);
1749+ DB = atomB->GetMndoDerivedParameterD(0);
1750+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1751+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1752+ double temp3 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1753+ Qzz, sQ, rhoA, rhoB, DA, DB, Rab);
1754+ DA = atomA->GetMndoDerivedParameterD(2);
1755+ DB = atomB->GetMndoDerivedParameterD(2);
1756+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1757+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1758+ double temp4 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1759+ Qzz, Qxx, rhoA, rhoB, DA, DB, Rab);
1760+ value = temp1 + temp2 + temp3 + temp4;
1761+ value *= dRabDa;
1762+ }
1763+ else if(mu == pz && nu == pz && lambda == py && sigma == py){
1764+ DA = atomA->GetMndoDerivedParameterD(0);
1765+ DB = atomB->GetMndoDerivedParameterD(0);
1766+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1767+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1768+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1769+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1770+ DA = atomA->GetMndoDerivedParameterD(0);
1771+ DB = atomB->GetMndoDerivedParameterD(2);
1772+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1773+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1774+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1775+ sQ, Qyy, rhoA, rhoB, DA, DB, Rab);
1776+ DA = atomA->GetMndoDerivedParameterD(2);
1777+ DB = atomB->GetMndoDerivedParameterD(0);
1778+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1779+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1780+ double temp3 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1781+ Qzz, sQ, rhoA, rhoB, DA, DB, Rab);
1782+ DA = atomA->GetMndoDerivedParameterD(2);
1783+ DB = atomB->GetMndoDerivedParameterD(2);
1784+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1785+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1786+ double temp4 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1787+ Qzz, Qyy, rhoA, rhoB, DA, DB, Rab);
1788+ value = temp1 + temp2 + temp3 + temp4;
1789+ value *= dRabDa;
1790+ }
1791+ // (37) in [DT_1977]
1792+ else if(mu == pz && nu == pz && lambda == pz && sigma == pz){
1793+ DA = atomA->GetMndoDerivedParameterD(0);
1794+ DB = atomB->GetMndoDerivedParameterD(0);
1795+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1796+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1797+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1798+ sQ, sQ, rhoA, rhoB, DA, DB, Rab);
1799+ DA = atomA->GetMndoDerivedParameterD(0);
1800+ DB = atomB->GetMndoDerivedParameterD(2);
1801+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1802+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1803+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1804+ sQ, Qzz, rhoA, rhoB, DA, DB, Rab);
1805+ DA = atomA->GetMndoDerivedParameterD(2);
1806+ DB = atomB->GetMndoDerivedParameterD(0);
1807+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1808+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1809+ double temp3 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1810+ Qzz, sQ, rhoA, rhoB, DA, DB, Rab);
1811+ DA = atomA->GetMndoDerivedParameterD(2);
1812+ DB = atomB->GetMndoDerivedParameterD(2);
1813+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1814+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1815+ double temp4 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1816+ Qzz, Qzz, rhoA, rhoB, DA, DB, Rab);
1817+ value = temp1 + temp2 + temp3 + temp4;
1818+ value *= dRabDa;
1819+ }
1820+ // (38) in [DT_1977]
1821+ else if(mu == s && nu == pz && lambda == s && sigma == s){
1822+ DA = atomA->GetMndoDerivedParameterD(1);
1823+ DB = atomB->GetMndoDerivedParameterD(0);
1824+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1825+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1826+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1827+ muz, sQ, rhoA, rhoB, DA, DB, Rab);
1828+ value = temp1;
1829+ value *= dRabDa;
1830+ }
1831+ else if(mu == pz && nu == s && lambda == s && sigma == s){
1832+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1833+ atomA, nu, mu, atomB, lambda, sigma, axisA);
1834+ }
1835+ // (39) in [DT_1977]
1836+ else if(mu == s && nu == pz && lambda == px && sigma == px){
1837+ DA = atomA->GetMndoDerivedParameterD(1);
1838+ DB = atomB->GetMndoDerivedParameterD(0);
1839+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1840+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1841+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1842+ muz, sQ, rhoA, rhoB, DA, DB, Rab);
1843+ DA = atomA->GetMndoDerivedParameterD(1);
1844+ DB = atomB->GetMndoDerivedParameterD(2);
1845+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1846+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1847+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1848+ muz, Qxx, rhoA, rhoB, DA, DB, Rab);
1849+ value = temp1 + temp2;
1850+ value *= dRabDa;
1851+ }
1852+ else if(mu == pz && nu == s && lambda == px && sigma == px){
1853+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1854+ atomA, nu, mu, atomB, lambda, sigma, axisA);
1855+ }
1856+ else if(mu == s && nu == pz && lambda == py && sigma == py){
1857+ DA = atomA->GetMndoDerivedParameterD(1);
1858+ DB = atomB->GetMndoDerivedParameterD(0);
1859+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1860+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1861+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1862+ muz, sQ, rhoA, rhoB, DA, DB, Rab);
1863+ DA = atomA->GetMndoDerivedParameterD(1);
1864+ DB = atomB->GetMndoDerivedParameterD(2);
1865+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1866+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1867+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1868+ muz, Qyy, rhoA, rhoB, DA, DB, Rab);
1869+ value = temp1 + temp2;
1870+ value *= dRabDa;
1871+ }
1872+ else if(mu == pz && nu == s && lambda == py && sigma == py){
1873+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1874+ atomA, nu, mu, atomB, lambda, sigma, axisA);
1875+ }
1876+ // (40) in [DT_1977]
1877+ else if(mu == s && nu == pz && lambda == pz && sigma == pz){
1878+ DA = atomA->GetMndoDerivedParameterD(1);
1879+ DB = atomB->GetMndoDerivedParameterD(0);
1880+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1881+ rhoB = atomB->GetMndoDerivedParameterRho(0);
1882+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1883+ muz, sQ, rhoA, rhoB, DA, DB, Rab);
1884+ DA = atomA->GetMndoDerivedParameterD(1);
1885+ DB = atomB->GetMndoDerivedParameterD(2);
1886+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1887+ rhoB = atomB->GetMndoDerivedParameterRho(2);
1888+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1889+ muz, Qzz, rhoA, rhoB, DA, DB, Rab);
1890+ value = temp1 + temp2;
1891+ value *= dRabDa;
1892+ }
1893+ else if(mu == pz && nu == s && lambda == pz && sigma == pz){
1894+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1895+ atomA, nu, mu, atomB, lambda, sigma, axisA);
1896+ }
1897+ // (41) in [DT_1977]
1898+ else if(mu == s && nu == s && lambda == s && sigma == pz){
1899+ DA = atomA->GetMndoDerivedParameterD(0);
1900+ DB = atomB->GetMndoDerivedParameterD(1);
1901+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1902+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1903+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1904+ sQ, muz, rhoA, rhoB, DA, DB, Rab);
1905+ value = temp1;
1906+ value *= dRabDa;
1907+ }
1908+ else if(mu == s && nu == s && lambda == pz && sigma == s){
1909+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1910+ atomA, mu, nu, atomB, sigma, lambda, axisA);
1911+ }
1912+ // (42) in [DT_1977]
1913+ else if(mu == px && nu == px && lambda == s && sigma == pz){
1914+ DA = atomA->GetMndoDerivedParameterD(0);
1915+ DB = atomB->GetMndoDerivedParameterD(1);
1916+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1917+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1918+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1919+ sQ, muz, rhoA, rhoB, DA, DB, Rab);
1920+ DA = atomA->GetMndoDerivedParameterD(2);
1921+ DB = atomB->GetMndoDerivedParameterD(1);
1922+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1923+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1924+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1925+ Qxx, muz, rhoA, rhoB, DA, DB, Rab);
1926+ value = temp1 + temp2;
1927+ value *= dRabDa;
1928+ }
1929+ else if(mu == px && nu == px && lambda == pz && sigma == s){
1930+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1931+ atomA, mu, nu, atomB, sigma, lambda, axisA);
1932+ }
1933+ else if(mu == py && nu == py && lambda == s && sigma == pz){
1934+ DA = atomA->GetMndoDerivedParameterD(0);
1935+ DB = atomB->GetMndoDerivedParameterD(1);
1936+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1937+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1938+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1939+ sQ, muz, rhoA, rhoB, DA, DB, Rab);
1940+ DA = atomA->GetMndoDerivedParameterD(2);
1941+ DB = atomB->GetMndoDerivedParameterD(1);
1942+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1943+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1944+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1945+ Qyy, muz, rhoA, rhoB, DA, DB, Rab);
1946+ value = temp1 + temp2;
1947+ value *= dRabDa;
1948+ }
1949+ else if(mu == py && nu == py && lambda == pz && sigma == s){
1950+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1951+ atomA, mu, nu, atomB, sigma, lambda, axisA);
1952+ }
1953+ // (43) in [DT_1977]
1954+ else if(mu == pz && nu == pz && lambda == s && sigma == pz){
1955+ DA = atomA->GetMndoDerivedParameterD(0);
1956+ DB = atomB->GetMndoDerivedParameterD(1);
1957+ rhoA = atomA->GetMndoDerivedParameterRho(0);
1958+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1959+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1960+ sQ, muz, rhoA, rhoB, DA, DB, Rab);
1961+ DA = atomA->GetMndoDerivedParameterD(2);
1962+ DB = atomB->GetMndoDerivedParameterD(1);
1963+ rhoA = atomA->GetMndoDerivedParameterRho(2);
1964+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1965+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1966+ Qzz, muz, rhoA, rhoB, DA, DB, Rab);
1967+ value = temp1 + temp2;
1968+ value *= dRabDa;
1969+ }
1970+ else if(mu == pz && nu == pz && lambda == pz && sigma == s){
1971+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1972+ atomA, mu, nu, atomB, sigma, lambda, axisA);
1973+ }
1974+ // (44) in [DT_1977]
1975+ else if(mu == s && nu == px && lambda == s && sigma == px){
1976+ DA = atomA->GetMndoDerivedParameterD(1);
1977+ DB = atomB->GetMndoDerivedParameterD(1);
1978+ rhoA = atomA->GetMndoDerivedParameterRho(1);
1979+ rhoB = atomB->GetMndoDerivedParameterRho(1);
1980+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
1981+ mux, mux, rhoA, rhoB, DA, DB, Rab);
1982+ value = temp1;
1983+ value *= dRabDa;
1984+ }
1985+ else if(mu == px && nu == s && lambda == s && sigma == px){
1986+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1987+ atomA, nu, mu, atomB, lambda, sigma, axisA);
1988+ }
1989+ else if(mu == s && nu == px && lambda == px && sigma == s){
1990+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1991+ atomA, mu, nu, atomB, sigma, lambda, axisA);
1992+ }
1993+ else if(mu == px && nu == s && lambda == px && sigma == s){
1994+ value = this->GetNddoRepulsionIntegralFirstDerivative(
1995+ atomA, nu, mu, atomB, sigma, lambda, axisA);
1996+ }
1997+ else if(mu == s && nu == py && lambda == s && sigma == py){
1998+ DA = atomA->GetMndoDerivedParameterD(1);
1999+ DB = atomB->GetMndoDerivedParameterD(1);
2000+ rhoA = atomA->GetMndoDerivedParameterRho(1);
2001+ rhoB = atomB->GetMndoDerivedParameterRho(1);
2002+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
2003+ muy, muy, rhoA, rhoB, DA, DB, Rab);
2004+ value = temp1;
2005+ value *= dRabDa;
2006+ }
2007+ else if(mu == py && nu == s && lambda == s && sigma == py){
2008+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2009+ atomA, nu, mu, atomB, lambda, sigma, axisA);
2010+ }
2011+ else if(mu == s && nu == py && lambda == py && sigma == s){
2012+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2013+ atomA, mu, nu, atomB, sigma, lambda, axisA);
2014+ }
2015+ else if(mu == py && nu == s && lambda == py && sigma == s){
2016+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2017+ atomA, nu, mu, atomB, sigma, lambda, axisA);
2018+ }
2019+ // (45) in [DT_1977]
2020+ else if(mu == s && nu == pz && lambda == s && sigma == pz){
2021+ DA = atomA->GetMndoDerivedParameterD(1);
2022+ DB = atomB->GetMndoDerivedParameterD(1);
2023+ rhoA = atomA->GetMndoDerivedParameterRho(1);
2024+ rhoB = atomB->GetMndoDerivedParameterRho(1);
2025+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
2026+ muz, muz, rhoA, rhoB, DA, DB, Rab);
2027+ value = temp1;
2028+ value *= dRabDa;
2029+ }
2030+ else if(mu == pz && nu == s && lambda == s && sigma == pz){
2031+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2032+ atomA, nu, mu, atomB, lambda, sigma, axisA);
2033+ }
2034+ else if(mu == s && nu == pz && lambda == pz && sigma == s){
2035+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2036+ atomA, mu, nu, atomB, sigma, lambda, axisA);
2037+ }
2038+ else if(mu == pz && nu == s && lambda == pz && sigma == s){
2039+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2040+ atomA, nu, mu, atomB, sigma, lambda, axisA);
2041+ }
2042+ // (46) in [DT_1977]
2043+ else if(mu == s && nu == px && lambda == px && sigma == pz){
2044+ DA = atomA->GetMndoDerivedParameterD(1);
2045+ DB = atomB->GetMndoDerivedParameterD(2);
2046+ rhoA = atomA->GetMndoDerivedParameterRho(1);
2047+ rhoB = atomB->GetMndoDerivedParameterRho(2);
2048+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
2049+ mux, Qxz, rhoA, rhoB, DA, DB, Rab);
2050+ value = temp1;
2051+ value *= dRabDa;
2052+ }
2053+ else if(mu == px && nu == s && lambda == px && sigma == pz){
2054+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2055+ atomA, nu, mu, atomB, lambda, sigma, axisA);
2056+ }
2057+ else if(mu == s && nu == px && lambda == pz && sigma == px){
2058+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2059+ atomA, mu, nu, atomB, sigma, lambda, axisA);
2060+ }
2061+ else if(mu == px && nu == s && lambda == pz && sigma == px){
2062+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2063+ atomA, nu, mu, atomB, sigma, lambda, axisA);
2064+ }
2065+ else if(mu == s && nu == py && lambda == py && sigma == pz){
2066+ DA = atomA->GetMndoDerivedParameterD(1);
2067+ DB = atomB->GetMndoDerivedParameterD(2);
2068+ rhoA = atomA->GetMndoDerivedParameterRho(1);
2069+ rhoB = atomB->GetMndoDerivedParameterRho(2);
2070+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
2071+ muy, Qyz, rhoA, rhoB, DA, DB, Rab);
2072+ value = temp1;
2073+ value *= dRabDa;
2074+ }
2075+ else if(mu == py && nu == s && lambda == py && sigma == pz){
2076+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2077+ atomA, nu, mu, atomB, lambda, sigma, axisA);
2078+ }
2079+ else if(mu == s && nu == py && lambda == pz && sigma == py){
2080+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2081+ atomA, mu, nu, atomB, sigma, lambda, axisA);
2082+ }
2083+ else if(mu == py && nu == s && lambda == pz && sigma == py){
2084+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2085+ atomA, nu, mu, atomB, sigma, lambda, axisA);
2086+ }
2087+ // (47) in [DT_1977]
2088+ else if(mu == px && nu == pz && lambda == s && sigma == px){
2089+ DA = atomA->GetMndoDerivedParameterD(2);
2090+ DB = atomB->GetMndoDerivedParameterD(1);
2091+ rhoA = atomA->GetMndoDerivedParameterRho(2);
2092+ rhoB = atomB->GetMndoDerivedParameterRho(1);
2093+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
2094+ Qxz, mux, rhoA, rhoB, DA, DB, Rab);
2095+ value = temp1;
2096+ value *= dRabDa;
2097+ }
2098+ else if(mu == pz && nu == px && lambda == s && sigma == px){
2099+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2100+ atomA, nu, mu, atomB, lambda, sigma, axisA);
2101+ }
2102+ else if(mu == px && nu == pz && lambda == px && sigma == s){
2103+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2104+ atomA, mu, nu, atomB, sigma, lambda, axisA);
2105+ }
2106+ else if(mu == pz && nu == px && lambda == px && sigma == s){
2107+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2108+ atomA, nu, mu, atomB, sigma, lambda, axisA);
2109+ }
2110+ else if(mu == py && nu == pz && lambda == s && sigma == py){
2111+ DA = atomA->GetMndoDerivedParameterD(2);
2112+ DB = atomB->GetMndoDerivedParameterD(1);
2113+ rhoA = atomA->GetMndoDerivedParameterRho(2);
2114+ rhoB = atomB->GetMndoDerivedParameterRho(1);
2115+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
2116+ Qyz, muy, rhoA, rhoB, DA, DB, Rab);
2117+ value = temp1;
2118+ value *= dRabDa;
2119+ }
2120+ else if(mu == pz && nu == py && lambda == s && sigma == py){
2121+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2122+ atomA, nu, mu, atomB, lambda, sigma, axisA);
2123+ }
2124+ else if(mu == py && nu == pz && lambda == py && sigma == s){
2125+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2126+ atomA, mu, nu, atomB, sigma, lambda, axisA);
2127+ }
2128+ else if(mu == pz && nu == py && lambda == py && sigma == s){
2129+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2130+ atomA, nu, mu, atomB, sigma, lambda, axisA);
2131+ }
2132+ // (48) in [DT_1977]
2133+ else if(mu == px && nu == pz && lambda == px && sigma == pz){
2134+ DA = atomA->GetMndoDerivedParameterD(2);
2135+ DB = atomB->GetMndoDerivedParameterD(2);
2136+ rhoA = atomA->GetMndoDerivedParameterRho(2);
2137+ rhoB = atomB->GetMndoDerivedParameterRho(2);
2138+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
2139+ Qxz, Qxz, rhoA, rhoB, DA, DB, Rab);
2140+ value = temp1;
2141+ value *= dRabDa;
2142+ }
2143+ else if(mu == pz && nu == px && lambda == px && sigma == pz){
2144+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2145+ atomA, nu, mu, atomB, lambda, sigma, axisA);
2146+ }
2147+ else if(mu == px && nu == pz && lambda == pz && sigma == px){
2148+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2149+ atomA, mu, nu, atomB, sigma, lambda, axisA);
2150+ }
2151+ else if(mu == pz && nu == px && lambda == pz && sigma == px){
2152+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2153+ atomA, nu, mu, atomB, sigma, lambda, axisA);
2154+ }
2155+ else if(mu == py && nu == pz && lambda == py && sigma == pz){
2156+ DA = atomA->GetMndoDerivedParameterD(2);
2157+ DB = atomB->GetMndoDerivedParameterD(2);
2158+ rhoA = atomA->GetMndoDerivedParameterRho(2);
2159+ rhoB = atomB->GetMndoDerivedParameterRho(2);
2160+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionFirstDerivative(
2161+ Qyz, Qyz, rhoA, rhoB, DA, DB, Rab);
2162+ value = temp1;
2163+ value *= dRabDa;
2164+ }
2165+ else if(mu == pz && nu == py && lambda == py && sigma == pz){
2166+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2167+ atomA, nu, mu, atomB, lambda, sigma, axisA);
2168+ }
2169+ else if(mu == py && nu == pz && lambda == pz && sigma == py){
2170+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2171+ atomA, mu, nu, atomB, sigma, lambda, axisA);
2172+ }
2173+ else if(mu == pz && nu == py && lambda == pz && sigma == py){
2174+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2175+ atomA, nu, mu, atomB, sigma, lambda, axisA);
2176+ }
2177+ // (49) in [DT_1977] and p19 in [MOPAC_1990]
2178+ else if(mu == px && nu == py && lambda == px && sigma == py){
2179+ value = 0.5*(this->GetNddoRepulsionIntegralFirstDerivative(
2180+ atomA, mu, mu, atomB, mu, mu, axisA)
2181+ -this->GetNddoRepulsionIntegralFirstDerivative(
2182+ atomA, mu, mu, atomB, nu, nu, axisA));
2183+ }
2184+ else if(mu == py && nu == px && lambda == px && sigma == py){
2185+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2186+ atomA, nu, mu, atomB, lambda, sigma, axisA);
2187+ }
2188+ else if(mu == px && nu == py && lambda == py && sigma == px){
2189+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2190+ atomA, mu, nu, atomB, sigma, lambda, axisA);
2191+ }
2192+ else if(mu == py && nu == px && lambda == py && sigma == px){
2193+ value = this->GetNddoRepulsionIntegralFirstDerivative(
2194+ atomA, nu, mu, atomB, sigma, lambda, axisA);
2195+ }
2196+ // d-orbitals
2197+ else if(mu == dxy || mu == dyz || mu == dzz || mu == dzx || mu == dxxyy ||
2198+ nu == dxy || nu == dyz || nu == dzz || nu == dzx || nu == dxxyy ||
2199+ lambda == dxy || lambda == dyz || lambda == dzz || lambda == dzx || lambda == dxxyy ||
2200+ sigma == dxy || sigma == dyz || sigma == dzz || sigma == dzx || sigma == dxxyy){
2201+
2202+ stringstream ss;
2203+ ss << this->errorMessageGetNddoRepulsionIntegralFirstDerivative;
2204+ ss << this->errorMessageAtomA << AtomTypeStr(atomA->GetAtomType()) << endl;
2205+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(mu) << endl;
2206+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(nu) << endl;
2207+ ss << this->errorMessageAtomB << AtomTypeStr(atomB->GetAtomType()) << endl;
2208+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(lambda) << endl;
2209+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(sigma) << endl;
2210+ throw MolDSException(ss.str());
2211+ }
2212+ else{
2213+ value = 0.0;
2214+ }
2215+ return value;
2216+}
2217+
2218+// See Apendix in [DT_1977]
2219+double Mndo::GetSemiEmpiricalMultipoleInteraction(MultipoleType multipoleA,
2220+ MultipoleType multipoleB,
2221+ double rhoA,
2222+ double rhoB,
2223+ double DA,
2224+ double DB,
2225+ double Rab){
2226+ double value = 0.0;
2227+ double a = rhoA + rhoB;
2228+
2229+ if(multipoleA == sQ && multipoleB == sQ){
2230+ value = pow(pow(Rab,2.0) + pow(a,2.0), -0.5);
2231+ }
2232+ else if(multipoleA == sQ && multipoleB == muz){
2233+ double temp1 = pow(Rab+DB,2.0) + pow(a,2.0);
2234+ double temp2 = pow(Rab-DB,2.0) + pow(a,2.0);
2235+ value = pow(temp1,-0.5)/2.0 - pow(temp2,-0.5)/2.0;
2236+ }
2237+ else if(multipoleA == muz && multipoleB == sQ){
2238+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2239+ rhoB, rhoA, DB, DA, Rab);
2240+ value *= pow(-1.0,1.0);
2241+ }
2242+ else if(multipoleA == sQ && multipoleB == Qxx){
2243+ double temp1 = pow(Rab,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
2244+ double temp2 = pow(Rab,2.0) + pow(a,2.0);
2245+ value = pow(temp1,-0.5)/2.0 - pow(temp2,-0.5)/2.0;
2246+ }
2247+ else if(multipoleA == Qxx && multipoleB == sQ){
2248+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2249+ rhoB, rhoA, DB, DA, Rab);
2250+ value *= pow(-1.0,2.0);
2251+ }
2252+ else if(multipoleA == sQ && multipoleB == Qyy){
2253+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleA, Qxx,
2254+ rhoA, rhoB, DA, DB, Rab);
2255+ }
2256+ else if(multipoleA == Qyy && multipoleB == sQ){
2257+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2258+ rhoB, rhoA, DB, DA, Rab);
2259+ value *= pow(-1.0,2.0);
2260+ }
2261+ else if(multipoleA == sQ && multipoleB == Qzz){
2262+ double temp1 = pow(Rab+2.0*DB,2.0) + pow(a,2.0);
2263+ double temp2 = pow(Rab,2.0) + pow(a,2.0);
2264+ double temp3 = pow(Rab-2.0*DB,2.0) + pow(a,2.0);
2265+ value = pow(temp1,-0.5)/4.0 - pow(temp2,-0.5)/2.0 + pow(temp3,-0.5)/4.0;
2266+ }
2267+ else if(multipoleA == Qzz && multipoleB == sQ){
2268+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2269+ rhoB, rhoA, DB, DA, Rab);
2270+ value *= pow(-1.0,2.0);
2271+ }
2272+ else if(multipoleA == mux && multipoleB == mux){
2273+ double temp1 = pow(Rab,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2274+ double temp2 = pow(Rab,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2275+ value = pow(temp1,-0.5)/2.0 - pow(temp2,-0.5)/2.0;
2276+ }
2277+ else if(multipoleA == muy && multipoleB == muy){
2278+ value = this->GetSemiEmpiricalMultipoleInteraction(mux, mux,
2279+ rhoA, rhoB, DA, DB, Rab);
2280+ }
2281+ else if(multipoleA == muz && multipoleB == muz){
2282+ double temp1 = pow(Rab+DA-DB,2.0) + pow(a,2.0);
2283+ double temp2 = pow(Rab+DA+DB,2.0) + pow(a,2.0);
2284+ double temp3 = pow(Rab-DA-DB,2.0) + pow(a,2.0);
2285+ double temp4 = pow(Rab-DA+DB,2.0) + pow(a,2.0);
2286+ value = pow(temp1,-0.5)/4.0 - pow(temp2,-0.5)/4.0
2287+ -pow(temp3,-0.5)/4.0 + pow(temp4,-0.5)/4.0;
2288+ }
2289+ else if(multipoleA == mux && multipoleB == Qxz){
2290+ double temp1 = pow(Rab-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2291+ double temp2 = pow(Rab-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2292+ double temp3 = pow(Rab+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2293+ double temp4 = pow(Rab+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2294+ value =-pow(temp1,-0.5)/4.0 + pow(temp2,-0.5)/4.0
2295+ +pow(temp3,-0.5)/4.0 - pow(temp4,-0.5)/4.0;
2296+ }
2297+ else if(multipoleA == Qxz && multipoleB == mux){
2298+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2299+ rhoB, rhoA, DB, DA, Rab);
2300+ value *= pow(-1.0,3.0);
2301+ }
2302+ else if(multipoleA == muy && multipoleB == Qyz){
2303+ value = this->GetSemiEmpiricalMultipoleInteraction(mux, Qxz,
2304+ rhoA, rhoB, DA, DB, Rab);
2305+ }
2306+ else if(multipoleA == Qyz && multipoleB == muy){
2307+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2308+ rhoB, rhoA, DB, DA, Rab);
2309+ value *= pow(-1.0,3.0);
15092310 }
15102311 else if(multipoleA == muz && multipoleB == Qxx){
15112312 double temp1 = pow(Rab+DA,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
@@ -1648,6 +2449,287 @@ double Mndo::GetSemiEmpiricalMultipoleInteraction(MultipoleType multipoleA,
16482449 return value;
16492450 }
16502451
2452+// First derivative of semiempirical multipole-multipole interactions.
2453+// This derivativ is related to the nuclear distance Rab.
2454+// See Apendix in [DT_1977]
2455+double Mndo::GetSemiEmpiricalMultipoleInteractionFirstDerivative(
2456+ MultipoleType multipoleA,
2457+ MultipoleType multipoleB,
2458+ double rhoA,
2459+ double rhoB,
2460+ double DA,
2461+ double DB,
2462+ double Rab){
2463+ double value = 0.0;
2464+ double a = rhoA + rhoB;
2465+
2466+ if(multipoleA == sQ && multipoleB == sQ){
2467+ value = -1.0*Rab*pow(pow(Rab,2.0) + pow(a,2.0), -1.5);
2468+ }
2469+ else if(multipoleA == sQ && multipoleB == muz){
2470+ double temp1 = pow(Rab+DB,2.0) + pow(a,2.0);
2471+ double temp2 = pow(Rab-DB,2.0) + pow(a,2.0);
2472+ value = (Rab+DB)*pow(temp1,-1.5)/2.0
2473+ -(Rab-DB)*pow(temp2,-1.5)/2.0;
2474+ value *= -1.0;
2475+ }
2476+ else if(multipoleA == muz && multipoleB == sQ){
2477+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2478+ rhoB, rhoA, DB, DA, Rab);
2479+ value *= pow(-1.0,1.0);
2480+ }
2481+ else if(multipoleA == sQ && multipoleB == Qxx){
2482+ double temp1 = pow(Rab,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
2483+ double temp2 = pow(Rab,2.0) + pow(a,2.0);
2484+ value = Rab*pow(temp1,-1.5)/2.0
2485+ -Rab*pow(temp2,-1.5)/2.0;
2486+ value *= -1.0;
2487+ }
2488+ else if(multipoleA == Qxx && multipoleB == sQ){
2489+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2490+ rhoB, rhoA, DB, DA, Rab);
2491+ value *= pow(-1.0,2.0);
2492+ }
2493+ else if(multipoleA == sQ && multipoleB == Qyy){
2494+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleA, Qxx,
2495+ rhoA, rhoB, DA, DB, Rab);
2496+ }
2497+ else if(multipoleA == Qyy && multipoleB == sQ){
2498+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2499+ rhoB, rhoA, DB, DA, Rab);
2500+ value *= pow(-1.0,2.0);
2501+ }
2502+ else if(multipoleA == sQ && multipoleB == Qzz){
2503+ double temp1 = pow(Rab+2.0*DB,2.0) + pow(a,2.0);
2504+ double temp2 = pow(Rab,2.0) + pow(a,2.0);
2505+ double temp3 = pow(Rab-2.0*DB,2.0) + pow(a,2.0);
2506+ value = (Rab+2.0*DB)*pow(temp1,-1.5)/4.0
2507+ -(Rab)*pow(temp2,-1.5)/2.0
2508+ +(Rab-2.0*DB)*pow(temp3,-1.5)/4.0;
2509+ value *= -1.0;
2510+ }
2511+ else if(multipoleA == Qzz && multipoleB == sQ){
2512+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2513+ rhoB, rhoA, DB, DA, Rab);
2514+ value *= pow(-1.0,2.0);
2515+ }
2516+ else if(multipoleA == mux && multipoleB == mux){
2517+ double temp1 = pow(Rab,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2518+ double temp2 = pow(Rab,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2519+ value = (Rab)*pow(temp1,-1.5)/2.0
2520+ -(Rab)*pow(temp2,-1.5)/2.0;
2521+ value *= -1.0;
2522+ }
2523+ else if(multipoleA == muy && multipoleB == muy){
2524+ value = this->GetSemiEmpiricalMultipoleInteraction(mux, mux,
2525+ rhoA, rhoB, DA, DB, Rab);
2526+ }
2527+ else if(multipoleA == muz && multipoleB == muz){
2528+ double temp1 = pow(Rab+DA-DB,2.0) + pow(a,2.0);
2529+ double temp2 = pow(Rab+DA+DB,2.0) + pow(a,2.0);
2530+ double temp3 = pow(Rab-DA-DB,2.0) + pow(a,2.0);
2531+ double temp4 = pow(Rab-DA+DB,2.0) + pow(a,2.0);
2532+ value = (Rab+DA-DB)*pow(temp1,-1.5)/4.0
2533+ -(Rab+DA+DB)*pow(temp2,-1.5)/4.0
2534+ -(Rab-DA-DB)*pow(temp3,-1.5)/4.0
2535+ +(Rab-DA+DB)*pow(temp4,-1.5)/4.0;
2536+ value *= -1.0;
2537+ }
2538+ else if(multipoleA == mux && multipoleB == Qxz){
2539+ double temp1 = pow(Rab-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2540+ double temp2 = pow(Rab-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2541+ double temp3 = pow(Rab+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2542+ double temp4 = pow(Rab+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2543+ value =-(Rab-DB)*pow(temp1,-1.5)/4.0
2544+ +(Rab-DB)*pow(temp2,-1.5)/4.0
2545+ +(Rab+DB)*pow(temp3,-1.5)/4.0
2546+ -(Rab+DB)*pow(temp4,-1.5)/4.0;
2547+ value *= -1.0;
2548+ }
2549+ else if(multipoleA == Qxz && multipoleB == mux){
2550+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2551+ rhoB, rhoA, DB, DA, Rab);
2552+ value *= pow(-1.0,3.0);
2553+ }
2554+ else if(multipoleA == muy && multipoleB == Qyz){
2555+ value = this->GetSemiEmpiricalMultipoleInteraction(mux, Qxz,
2556+ rhoA, rhoB, DA, DB, Rab);
2557+ }
2558+ else if(multipoleA == Qyz && multipoleB == muy){
2559+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2560+ rhoB, rhoA, DB, DA, Rab);
2561+ value *= pow(-1.0,3.0);
2562+ }
2563+ else if(multipoleA == muz && multipoleB == Qxx){
2564+ double temp1 = pow(Rab+DA,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
2565+ double temp2 = pow(Rab-DA,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
2566+ double temp3 = pow(Rab+DA,2.0) + pow(a,2.0);
2567+ double temp4 = pow(Rab-DA,2.0) + pow(a,2.0);
2568+ value =-(Rab+DA)*pow(temp1,-1.5)/4.0
2569+ +(Rab-DA)*pow(temp2,-1.5)/4.0
2570+ +(Rab+DA)*pow(temp3,-1.5)/4.0
2571+ -(Rab-DA)*pow(temp4,-1.5)/4.0;
2572+ value *= -1.0;
2573+ }
2574+ else if(multipoleA == Qxx && multipoleB == muz){
2575+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2576+ rhoB, rhoA, DB, DA, Rab);
2577+ value *= pow(-1.0,3.0);
2578+ }
2579+ else if(multipoleA == muz && multipoleB == Qyy){
2580+ value = this->GetSemiEmpiricalMultipoleInteraction(muz, Qxx,
2581+ rhoA, rhoB, DA, DB, Rab);
2582+ }
2583+ else if(multipoleA == Qyy && multipoleB == muz){
2584+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2585+ rhoB, rhoA, DB, DA, Rab);
2586+ value *= pow(-1.0,3.0);
2587+ }
2588+ else if(multipoleA == muz && multipoleB == Qzz){
2589+ double temp1 = pow(Rab+DA-2.0*DB,2.0) + pow(a,2.0);
2590+ double temp2 = pow(Rab-DA-2.0*DB,2.0) + pow(a,2.0);
2591+ double temp3 = pow(Rab+DA+2.0*DB,2.0) + pow(a,2.0);
2592+ double temp4 = pow(Rab-DA+2.0*DB,2.0) + pow(a,2.0);
2593+ double temp5 = pow(Rab+DA,2.0) + pow(a,2.0);
2594+ double temp6 = pow(Rab-DA,2.0) + pow(a,2.0);
2595+ value =-(Rab+DA-2.0*DB)*pow(temp1,-1.5)/8.0
2596+ +(Rab-DA-2.0*DB)*pow(temp2,-1.5)/8.0
2597+ -(Rab+DA+2.0*DB)*pow(temp3,-1.5)/8.0
2598+ +(Rab-DA+2.0*DB)*pow(temp4,-1.5)/8.0
2599+ +(Rab+DA )*pow(temp5,-1.5)/4.0
2600+ -(Rab-DA )*pow(temp6,-1.5)/4.0;
2601+ value *= -1.0;
2602+ }
2603+ else if(multipoleA == Qzz && multipoleB == muz){
2604+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2605+ rhoB, rhoA, DB, DA, Rab);
2606+ value *= pow(-1.0,3.0);
2607+ }
2608+ else if(multipoleA == Qxx && multipoleB == Qxx){
2609+ double temp1 = pow(Rab,2.0) + 4.0*pow(DA-DB,2.0) + pow(a,2.0);
2610+ double temp2 = pow(Rab,2.0) + 4.0*pow(DA+DB,2.0) + pow(a,2.0);
2611+ double temp3 = pow(Rab,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
2612+ double temp4 = pow(Rab,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
2613+ double temp5 = pow(Rab,2.0) + pow(a,2.0);
2614+ value = (Rab)*pow(temp1,-1.5)/8.0
2615+ +(Rab)*pow(temp2,-1.5)/8.0
2616+ -(Rab)*pow(temp3,-1.5)/4.0
2617+ -(Rab)*pow(temp4,-1.5)/4.0
2618+ +(Rab)*pow(temp5,-1.5)/4.0;
2619+ value *= -1.0;
2620+ }
2621+ else if(multipoleA == Qyy && multipoleB == Qyy){
2622+ value = this->GetSemiEmpiricalMultipoleInteraction(Qxx, Qxx,
2623+ rhoA, rhoB, DA, DB, Rab);
2624+ }
2625+ else if(multipoleA == Qxx && multipoleB == Qyy){
2626+ double temp1 = pow(Rab,2.0) + pow(2.0*DA,2.0) + pow(2.0*DB,2.0)+ pow(a,2.0);
2627+ double temp2 = pow(Rab,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
2628+ double temp3 = pow(Rab,2.0) + pow(2.0*DB,2.0) + pow(a,2.0);
2629+ double temp4 = pow(Rab,2.0) + pow(a,2.0);
2630+ value = (Rab)*pow(temp1,-1.5)/4.0
2631+ -(Rab)*pow(temp2,-1.5)/4.0
2632+ -(Rab)*pow(temp3,-1.5)/4.0
2633+ +(Rab)*pow(temp4,-1.5)/4.0;
2634+ value *= -1.0;
2635+ }
2636+ else if(multipoleA == Qyy && multipoleB == Qxx){
2637+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2638+ rhoB, rhoA, DB, DA, Rab);
2639+ value *= pow(-1.0,4.0);
2640+ }
2641+ else if(multipoleA == Qxx && multipoleB == Qzz){
2642+ double temp1 = pow(Rab-2.0*DB,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
2643+ double temp2 = pow(Rab+2.0*DB,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
2644+ double temp3 = pow(Rab-2.0*DB,2.0) + pow(a,2.0);
2645+ double temp4 = pow(Rab+2.0*DB,2.0) + pow(a,2.0);
2646+ double temp5 = pow(Rab,2.0) + pow(2.0*DA,2.0) + pow(a,2.0);
2647+ double temp6 = pow(Rab,2.0) + pow(a,2.0);
2648+ value = (Rab-2.0*DB)*pow(temp1,-1.5)/8.0
2649+ +(Rab+2.0*DB)*pow(temp2,-1.5)/8.0
2650+ -(Rab-2.0*DB)*pow(temp3,-1.5)/8.0
2651+ -(Rab+2.0*DB)*pow(temp4,-1.5)/8.0
2652+ -(Rab )*pow(temp5,-1.5)/4.0
2653+ +(Rab )*pow(temp6,-1.5)/4.0;
2654+ value *= -1.0;
2655+ }
2656+ else if(multipoleA == Qzz && multipoleB == Qxx){
2657+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2658+ rhoB, rhoA, DB, DA, Rab);
2659+ value *= pow(-1.0,4.0);
2660+ }
2661+ else if(multipoleA == Qyy && multipoleB == Qzz){
2662+ value = this->GetSemiEmpiricalMultipoleInteraction(Qxx, multipoleB,
2663+ rhoA, rhoB, DA, DB, Rab);
2664+ }
2665+ else if(multipoleA == Qzz && multipoleB == Qyy){
2666+ value = this->GetSemiEmpiricalMultipoleInteraction(multipoleB, multipoleA,
2667+ rhoB, rhoA, DB, DA, Rab);
2668+ value *= pow(-1.0,4.0);
2669+ }
2670+ else if(multipoleA == Qzz && multipoleB == Qzz){
2671+ double temp1 = pow(Rab+2.0*DA-2.0*DB,2.0) + pow(a,2.0);
2672+ double temp2 = pow(Rab+2.0*DA+2.0*DB,2.0) + pow(a,2.0);
2673+ double temp3 = pow(Rab-2.0*DA-2.0*DB,2.0) + pow(a,2.0);
2674+ double temp4 = pow(Rab-2.0*DA+2.0*DB,2.0) + pow(a,2.0);
2675+ double temp5 = pow(Rab+2.0*DA,2.0) + pow(a,2.0);
2676+ double temp6 = pow(Rab-2.0*DA,2.0) + pow(a,2.0);
2677+ double temp7 = pow(Rab+2.0*DB,2.0) + pow(a,2.0);
2678+ double temp8 = pow(Rab-2.0*DB,2.0) + pow(a,2.0);
2679+ double temp9 = pow(Rab,2.0) + pow(a,2.0);
2680+ value = (Rab+2.0*DA-2.0*DB)*pow(temp1,-1.5)/16.0
2681+ +(Rab+2.0*DA+2.0*DB)*pow(temp2,-1.5)/16.0
2682+ +(Rab-2.0*DA-2.0*DB)*pow(temp3,-1.5)/16.0
2683+ +(Rab-2.0*DA+2.0*DB)*pow(temp4,-1.5)/16.0
2684+ -(Rab+2.0*DA)*pow(temp5,-1.5)/8.0
2685+ -(Rab-2.0*DA)*pow(temp6,-1.5)/8.0
2686+ -(Rab+2.0*DB)*pow(temp7,-1.5)/8.0
2687+ -(Rab-2.0*DB)*pow(temp8,-1.5)/8.0
2688+ +(Rab)*pow(temp9,-1.5)/4.0;
2689+ value *= -1.0;
2690+ }
2691+ else if(multipoleA == Qxz && multipoleB == Qxz){
2692+ double temp1 = pow(Rab+DA-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2693+ double temp2 = pow(Rab+DA-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2694+ double temp3 = pow(Rab+DA+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2695+ double temp4 = pow(Rab+DA+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2696+ double temp5 = pow(Rab-DA-DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2697+ double temp6 = pow(Rab-DA-DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2698+ double temp7 = pow(Rab-DA+DB,2.0) + pow(DA-DB,2.0) + pow(a,2.0);
2699+ double temp8 = pow(Rab-DA+DB,2.0) + pow(DA+DB,2.0) + pow(a,2.0);
2700+ value = (Rab+DA-DB)*pow(temp1,-1.5)/8.0
2701+ -(Rab+DA-DB)*pow(temp2,-1.5)/8.0
2702+ -(Rab+DA+DB)*pow(temp3,-1.5)/8.0
2703+ +(Rab+DA+DB)*pow(temp4,-1.5)/8.0
2704+ -(Rab-DA-DB)*pow(temp5,-1.5)/8.0
2705+ +(Rab-DA-DB)*pow(temp6,-1.5)/8.0
2706+ +(Rab-DA+DB)*pow(temp7,-1.5)/8.0
2707+ -(Rab-DA+DB)*pow(temp8,-1.5)/8.0;
2708+ value *= -1.0;
2709+ }
2710+ else if(multipoleA == Qyz && multipoleB == Qyz){
2711+ value = this->GetSemiEmpiricalMultipoleInteraction(Qxz, Qxz,
2712+ rhoA, rhoB, DA, DB, Rab);
2713+ }
2714+ else if(multipoleA == Qxy && multipoleB == Qxy){
2715+ double temp1 = pow(Rab,2.0) + 2.0*pow(DA-DB,2.0) + pow(a,2.0);
2716+ double temp2 = pow(Rab,2.0) + 2.0*pow(DA+DB,2.0) + pow(a,2.0);
2717+ double temp3 = pow(Rab,2.0) + 2.0*pow(DA,2.0) + 2.0*pow(DB,2.0) + pow(a,2.0);
2718+ value = (Rab)*pow(temp1,-1.5)/4.0
2719+ +(Rab)*pow(temp2,-1.5)/4.0
2720+ -(Rab)*pow(temp3,-1.5)/2.0;
2721+ value *= -1.0;
2722+ }
2723+ else{
2724+ stringstream ss;
2725+ ss << this->errorMessageGetSemiEmpiricalMultipoleInteractionFirstDeriBadMultipoles;
2726+ ss << this->errorMessageMultipoleA << MultipoleTypeStr(multipoleA) << endl;
2727+ ss << this->errorMessageMultipoleB << MultipoleTypeStr(multipoleB) << endl;
2728+ throw MolDSException(ss.str());
2729+ }
2730+ return value;
2731+}
2732+
16512733 }
16522734 #endif
16532735