OSDN -- Desarrollar y descargar software de código abierto
Translation Status of Español translation status for es
  • 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

Revisión3f3330ea0ab811198347c254ba0d7688a5f316c3 (tree)
Tiempo2012-07-06 03:53:02
AutorMikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Log Message

Some methods to calculate Hessian elements are added. #28554

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

Cambiar Resumen

Diferencia incremental

--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -1807,9 +1807,106 @@ void Mndo::CalcXiMatrices(double** xiOcc,
18071807 }
18081808 }
18091809
1810+// mu and nu is included in atomA' AO.
1811+// s is included in atomB's AO.
1812+double Mndo::GetAuxiliaryHessianElement1(int mu,
1813+ int nu,
1814+ int firstAOIndexA,
1815+ MolDS_base::CartesianType axisA1,
1816+ MolDS_base::CartesianType axisA2,
1817+ int atomBIndex,
1818+ double const* const* orbitalElectronPopulation,
1819+ double const* const* const* const* const* const* const* diatomicTwoElecTwoCoreSecondDerivs) const{
1820+ double value = orbitalElectronPopulation[mu]
1821+ [nu]
1822+ *diatomicTwoElecTwoCoreSecondDerivs[atomBIndex]
1823+ [mu-firstAOIndexA]
1824+ [nu-firstAOIndexA]
1825+ [s]
1826+ [s]
1827+ [axisA1]
1828+ [axisA2];
1829+ return value*this->molecule->GetAtom(atomBIndex)->GetCoreCharge();
1830+}
1831+
1832+// mu and nu is included in atomA' AO.
1833+// s is included in atomB's AO.
1834+double Mndo::GetAuxiliaryHessianElement2(int mu,
1835+ int nu,
1836+ int firstAOIndexA,
1837+ MolDS_base::CartesianType axisA1,
1838+ MolDS_base::CartesianType axisA2,
1839+ int atomAIndex,
1840+ int atomBIndex,
1841+ double const* const* const* const* orbitalElectronPopulationFirstDerivs,
1842+ double const* const* const* const* const* const* diatomicTwoElecTwoCoreFirstDerivs) const{
1843+ double value = orbitalElectronPopulationFirstDerivs[mu]
1844+ [nu]
1845+ [atomAIndex]
1846+ [axisA2]
1847+ *diatomicTwoElecTwoCoreFirstDerivs[atomBIndex]
1848+ [mu-firstAOIndexA]
1849+ [nu-firstAOIndexA]
1850+ [s]
1851+ [s]
1852+ [axisA1];
1853+ return value*this->molecule->GetAtom(atomBIndex)->GetCoreCharge();
1854+}
1855+
1856+// lambda and sigma is included in atomB' AO.
1857+// s is included in atomA's AO.
1858+double Mndo::GetAuxiliaryHessianElement3(int lambda,
1859+ int sigma,
1860+ int firstAOIndexB,
1861+ MolDS_base::CartesianType axisA1,
1862+ MolDS_base::CartesianType axisA2,
1863+ int atomAIndex,
1864+ int atomBIndex,
1865+ double const* const* orbitalElectronPopulation,
1866+ double const* const* const* const* const* const* const* diatomicTwoElecTwoCoreSecondDerivs) const{
1867+ double value = orbitalElectronPopulation[lambda]
1868+ [sigma]
1869+ *diatomicTwoElecTwoCoreSecondDerivs[atomBIndex]
1870+ [s]
1871+ [s]
1872+ [lambda-firstAOIndexB]
1873+ [sigma-firstAOIndexB]
1874+ [axisA1]
1875+ [axisA2];
1876+ return value*this->molecule->GetAtom(atomAIndex)->GetCoreCharge();
1877+}
1878+
1879+// lambda and sigma is included in atomB' AO.
1880+// s is included in atomA's AO.
1881+double Mndo::GetAuxiliaryHessianElement4(int lambda,
1882+ int sigma,
1883+ int firstAOIndexB,
1884+ MolDS_base::CartesianType axisA1,
1885+ MolDS_base::CartesianType axisA2,
1886+ int atomAIndex,
1887+ int atomBIndex,
1888+ double const* const* const* const* orbitalElectronPopulationFirstDerivs,
1889+ double const* const* const* const* const* const* diatomicTwoElecTwoCoreFirstDerivs) const{
1890+ double value = orbitalElectronPopulationFirstDerivs[lambda]
1891+ [sigma]
1892+ [atomAIndex]
1893+ [axisA2]
1894+ *diatomicTwoElecTwoCoreFirstDerivs[atomBIndex]
1895+ [s]
1896+ [s]
1897+ [lambda-firstAOIndexB]
1898+ [sigma-firstAOIndexB]
1899+ [axisA1];
1900+ return value*this->molecule->GetAtom(atomAIndex)->GetCoreCharge();
1901+}
1902+
18101903 void Mndo::CalcHessianSCF(double** hessianSCF) const{
18111904 int totalNumberAOs = this->molecule->GetTotalNumberAOs();
18121905 double**** orbitalElectronPopulationFirstDerivatives = NULL;
1906+ double****** diatomicTwoElecTwoCoreFirstDerivs = NULL;
1907+ double******* diatomicTwoElecTwoCoreSecondDerivs = NULL;
1908+ double**** diatomicOverlapFirstDerivs = NULL;
1909+ double***** diatomicOverlapSecondDerivs = NULL;
18131910
18141911 try{
18151912 MallocerFreer::GetInstance()->Malloc<double>(&orbitalElectronPopulationFirstDerivatives,
@@ -1822,13 +1919,80 @@ void Mndo::CalcHessianSCF(double** hessianSCF) const{
18221919
18231920 for(int atomAIndex=0; atomAIndex<this->molecule->GetNumberAtoms(); atomAIndex++){
18241921 for(int axisA = XAxis; axisA<CartesianType_end; axisA++){
1825- int k = atomAIndex*CartesianType_end + axisA;
1922+ int firstAOIndexA = this->molecule->GetAtom(atomAIndex)->GetFirstAOIndex();
1923+ int numberAOsA = this->molecule->GetAtom(atomAIndex)->GetValenceSize();
1924+ int k = atomAIndex*CartesianType_end + axisA; // hessian index, i.e. hessian[k][l]
1925+
1926+ // calculation of derivatives of the overlaps and two electron integrals
1927+ for(int atomBIndex=0; atomBIndex<this->molecule->GetNumberAtoms(); atomBIndex++){
1928+ if(atomAIndex != atomBIndex){
1929+ this->CalcDiatomicOverlapFirstDerivatives(diatomicOverlapFirstDerivs[atomBIndex],
1930+ atomAIndex,
1931+ atomBIndex);
1932+ this->CalcDiatomicOverlapSecondDerivatives(diatomicOverlapSecondDerivs[atomBIndex],
1933+ atomAIndex,
1934+ atomBIndex);
1935+ this->CalcDiatomicTwoElecTwoCoreFirstDerivatives(diatomicTwoElecTwoCoreFirstDerivs[atomBIndex],
1936+ atomAIndex,
1937+ atomBIndex);
1938+ this->CalcDiatomicTwoElecTwoCoreSecondDerivatives(diatomicTwoElecTwoCoreSecondDerivs[atomBIndex],
1939+ atomAIndex,
1940+ atomBIndex);
1941+ }
1942+ }
18261943
1827- // atomA == atomB
1944+ // hessian element (atomA == atomB)
18281945 for(int axisA2 = axisA; axisA2<CartesianType_end; axisA2++){
1829- int l = atomAIndex*CartesianType_end + axisA2;
1946+ int l = atomAIndex*CartesianType_end + axisA2; // hessian index, i.e. hessian[k][l]
18301947 hessianSCF[k][l] = 0.0;
18311948 for(int atomCIndex=0; atomCIndex<this->molecule->GetNumberAtoms(); atomCIndex++){
1949+ int firstAOIndexC = this->molecule->GetAtom(atomCIndex)->GetFirstAOIndex();
1950+ int numberAOsC = this->molecule->GetAtom(atomCIndex)->GetValenceSize();
1951+ // second derivative of electronic part
1952+ for(int mu=firstAOIndexA; mu<firstAOIndexA+numberAOsA; mu++){
1953+ for(int nu=firstAOIndexA; nu<firstAOIndexA+numberAOsA; nu++){
1954+ hessianSCF[k][l] -= this->GetAuxiliaryHessianElement1(mu,
1955+ nu,
1956+ firstAOIndexA,
1957+ static_cast<CartesianType>(axisA),
1958+ static_cast<CartesianType>(axisA2),
1959+ atomCIndex,
1960+ this->orbitalElectronPopulation,
1961+ diatomicTwoElecTwoCoreSecondDerivs);
1962+ hessianSCF[k][l] -= this->GetAuxiliaryHessianElement2(mu,
1963+ nu,
1964+ firstAOIndexA,
1965+ static_cast<CartesianType>(axisA),
1966+ static_cast<CartesianType>(axisA2),
1967+ atomAIndex,
1968+ atomCIndex,
1969+ orbitalElectronPopulationFirstDerivatives,
1970+ diatomicTwoElecTwoCoreFirstDerivs);
1971+ }
1972+ }
1973+ for(int lambda=firstAOIndexC; lambda<firstAOIndexC+numberAOsC; lambda++){
1974+ for(int sigma=firstAOIndexC; sigma<firstAOIndexC+numberAOsC; sigma++){
1975+ hessianSCF[k][l] -= this->GetAuxiliaryHessianElement3(lambda,
1976+ sigma,
1977+ firstAOIndexC,
1978+ static_cast<CartesianType>(axisA),
1979+ static_cast<CartesianType>(axisA2),
1980+ atomAIndex,
1981+ atomCIndex,
1982+ this->orbitalElectronPopulation,
1983+ diatomicTwoElecTwoCoreSecondDerivs);
1984+ hessianSCF[k][l] -= this->GetAuxiliaryHessianElement4(lambda,
1985+ sigma,
1986+ firstAOIndexC,
1987+ static_cast<CartesianType>(axisA),
1988+ static_cast<CartesianType>(axisA2),
1989+ atomAIndex,
1990+ atomCIndex,
1991+ orbitalElectronPopulationFirstDerivatives,
1992+ diatomicTwoElecTwoCoreFirstDerivs);
1993+ }
1994+ }
1995+
18321996 // second derivatives of the nuclear repulsions
18331997 hessianSCF[k][l] += this->GetDiatomCoreRepulsionSecondDerivative(atomAIndex,
18341998 atomCIndex,
@@ -1845,7 +2009,7 @@ void Mndo::CalcHessianSCF(double** hessianSCF) const{
18452009 }
18462010 }
18472011
1848- // atomA < atomB
2012+ // hessian element atomA < atomB
18492013 for(int atomBIndex=atomAIndex+1; atomBIndex<this->molecule->GetNumberAtoms(); atomBIndex++){
18502014 for(int axisB = XAxis; axisB<CartesianType_end; axisB++){
18512015 int l = atomAIndex*CartesianType_end + axisB;
--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -208,6 +208,41 @@ private:
208208 int numberActiveOcc,
209209 int numberActiveVir) const;
210210 void CalcHessianSCF(double** hessianSCF) const;
211+ double GetAuxiliaryHessianElement1(int mu,
212+ int nu,
213+ int firstAOIndexA,
214+ MolDS_base::CartesianType axisA1,
215+ MolDS_base::CartesianType axisA2,
216+ int atomBIndex,
217+ double const* const* orbitalElectronPopulation,
218+ double const* const* const* const* const* const* const* diatomicTwoElecTwoCoreSecondDerivs) const;
219+ double GetAuxiliaryHessianElement2(int mu,
220+ int nu,
221+ int firstAOIndexA,
222+ MolDS_base::CartesianType axisA1,
223+ MolDS_base::CartesianType axisA2,
224+ int atomAIndex,
225+ int atomBIndex,
226+ double const* const* const* const* orbitalElectronPopulationFirstDerivs,
227+ double const* const* const* const* const* const* diatomicTwoElecTwoCoreFirstDerivs) const;
228+ double GetAuxiliaryHessianElement3(int lambda,
229+ int sigma,
230+ int firstAOIndexB,
231+ MolDS_base::CartesianType axisA1,
232+ MolDS_base::CartesianType axisA2,
233+ int atomAIndex,
234+ int atomBIndex,
235+ double const* const* orbitalElectronPopulation,
236+ double const* const* const* const* const* const* const* diatomicTwoElecTwoCoreSecondDerivs) const;
237+ double GetAuxiliaryHessianElement4(int lambda,
238+ int sigma,
239+ int firstAOIndexB,
240+ MolDS_base::CartesianType axisA1,
241+ MolDS_base::CartesianType axisA2,
242+ int atomAIndex,
243+ int atomBIndex,
244+ double const* const* const* const* orbitalElectronPopulationFirstDerivs,
245+ double const* const* const* const* const* const* diatomicTwoElecTwoCoreFirstDerivs) const;
211246 void CalcOrbitalElectronPopulationFirstDerivatives(double**** orbitalElectronPopulationFirstDerivatives) const;
212247 void SolveCPHF(double** solutionsCPHF,
213248 const std::vector<MoIndexPair>& nonRedundantQIndeces,