OSDN -- Desenvolver e Download de Software Open Source
Estado tradução de Português translation status for pt
  • 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

Revisão8d6fa7e9025c204e893bf1f0c657b6114e3999a6 (tree)
Hora2013-08-06 06:24:12
AutorMikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Mensagem de Log

Refactoring CalcDiatomicOverlapAOs1stDerivatives in ZindoS and Mndo. #31814

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

Mudança Sumário

Diff

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -3781,8 +3781,15 @@ void Cndo2::CalcOverlapAOs(double** overlapAOs, const Molecule& molecule) const{
37813781 // Note that this method can not treat d-obitals
37823782 // because CalcRotatingMatrix1stDerivatives can not treat d-orbitals.
37833783 void Cndo2::CalcDiatomicOverlapAOs1stDerivatives(double*** diatomicOverlapAOs1stDerivs,
3784- const Atom& atomA,
3785- const Atom& atomB) const{
3784+ double** tmpDiaOverlapAOsInDiaFrame,
3785+ double** tmpDiaOverlapAOs1stDerivInDiaFrame,
3786+ double** tmpRotMat,
3787+ double** tmpRotMat1stDeriv,
3788+ double*** tmpRotMat1stDerivs,
3789+ double** tmpRotatedDiatomicOverlap,
3790+ double** tmpMatrix,
3791+ const Atom& atomA,
3792+ const Atom& atomB) const{
37863793 double cartesian[CartesianType_end] = {atomA.GetXyz()[XAxis] - atomB.GetXyz()[XAxis],
37873794 atomA.GetXyz()[YAxis] - atomB.GetXyz()[YAxis],
37883795 atomA.GetXyz()[ZAxis] - atomB.GetXyz()[ZAxis]};
@@ -3790,134 +3797,112 @@ void Cndo2::CalcDiatomicOverlapAOs1stDerivatives(double*** diatomicOverlapAOs1st
37903797 pow(cartesian[YAxis],2.0) +
37913798 pow(cartesian[ZAxis],2.0) );
37923799
3793- double** diaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
3794- double** diaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
3795- double** rotMat = NULL; // rotating Matrix from the diatomic frame to space fixed frame.
3796- double*** rotMat1stDerivs = NULL; // first derivatives of the rotMat.
3797- double** rotatedDiatomicOverlap = NULL;
3798- double** tmpRotMat1stDeriv = NULL;
3799- double** tmpMatrix = NULL;
3800-
3801- try{
3802- this->MallocDiatomicOverlapAOs1stDeriTemps(&diaOverlapAOsInDiaFrame,
3803- &diaOverlapAOs1stDerivInDiaFrame,
3804- &rotMat,
3805- &rotMat1stDerivs,
3806- &rotatedDiatomicOverlap,
3807- &tmpRotMat1stDeriv,
3808- &tmpMatrix);
3809- this->CalcDiatomicOverlapAOsInDiatomicFrame(diaOverlapAOsInDiaFrame, atomA, atomB);
3810- this->CalcDiatomicOverlapAOs1stDerivativeInDiatomicFrame(diaOverlapAOs1stDerivInDiaFrame, atomA, atomB);
3811- this->CalcRotatingMatrix(rotMat, atomA, atomB);
3812- this->CalcRotatingMatrix1stDerivatives(rotMat1stDerivs, atomA, atomB);
3800+ this->CalcDiatomicOverlapAOsInDiatomicFrame(tmpDiaOverlapAOsInDiaFrame, atomA, atomB);
3801+ this->CalcDiatomicOverlapAOs1stDerivativeInDiatomicFrame(tmpDiaOverlapAOs1stDerivInDiaFrame, atomA, atomB);
3802+ this->CalcRotatingMatrix(tmpRotMat, atomA, atomB);
3803+ this->CalcRotatingMatrix1stDerivatives(tmpRotMat1stDerivs, atomA, atomB);
3804+
3805+ // rotate (fast algorithm, see also slow algorithm shown later)
3806+ int incrementOne = 1;
3807+ bool isColumnMajorRotatingMatrix = false;
3808+ bool isColumnMajorDiaOverlapAOs = false;
3809+ double alpha = 0.0;
3810+ double beta = 0.0;
3811+ for(int c=0; c<CartesianType_end; c++){
3812+ MolDS_wrappers::Blas::GetInstance()->Dcopy(OrbitalType_end*OrbitalType_end,
3813+ &tmpRotMat1stDerivs[0][0][c], CartesianType_end,
3814+ &tmpRotMat1stDeriv[0][0], incrementOne);
3815+ alpha = cartesian[c]/R;
3816+ beta = 0.0;
3817+ MolDS_wrappers::Blas::GetInstance()->Dgemmm(isColumnMajorRotatingMatrix,
3818+ isColumnMajorDiaOverlapAOs,
3819+ !isColumnMajorRotatingMatrix,
3820+ OrbitalType_end, OrbitalType_end, OrbitalType_end, OrbitalType_end,
3821+ alpha,
3822+ tmpRotMat,
3823+ tmpDiaOverlapAOs1stDerivInDiaFrame,
3824+ tmpRotMat,
3825+ beta,
3826+ tmpRotatedDiatomicOverlap,
3827+ tmpMatrix);
3828+ alpha = 1.0;
3829+ beta = 1.0;
3830+ MolDS_wrappers::Blas::GetInstance()->Dgemmm(isColumnMajorRotatingMatrix,
3831+ isColumnMajorDiaOverlapAOs,
3832+ !isColumnMajorRotatingMatrix,
3833+ OrbitalType_end, OrbitalType_end, OrbitalType_end, OrbitalType_end,
3834+ alpha,
3835+ tmpRotMat1stDeriv,
3836+ tmpDiaOverlapAOsInDiaFrame,
3837+ tmpRotMat,
3838+ beta,
3839+ tmpRotatedDiatomicOverlap,
3840+ tmpMatrix);
3841+ MolDS_wrappers::Blas::GetInstance()->Dgemmm(isColumnMajorRotatingMatrix,
3842+ isColumnMajorDiaOverlapAOs,
3843+ !isColumnMajorRotatingMatrix,
3844+ OrbitalType_end, OrbitalType_end, OrbitalType_end, OrbitalType_end,
3845+ alpha,
3846+ tmpRotMat,
3847+ tmpDiaOverlapAOsInDiaFrame,
3848+ tmpRotMat1stDeriv,
3849+ beta,
3850+ tmpRotatedDiatomicOverlap,
3851+ tmpMatrix);
3852+ MolDS_wrappers::Blas::GetInstance()->Dcopy(OrbitalType_end*OrbitalType_end,
3853+ &tmpRotatedDiatomicOverlap[0][0], incrementOne,
3854+ &diatomicOverlapAOs1stDerivs[0][0][c], CartesianType_end);
3855+ }
38133856
3814- // rotate (fast algorithm, see also slow algorithm shown later)
3815- int incrementOne = 1;
3816- bool isColumnMajorRotatingMatrix = false;
3817- bool isColumnMajorDiaOverlapAOs = false;
3818- double alpha = 0.0;
3819- double beta = 0.0;
3820- for(int c=0; c<CartesianType_end; c++){
3821- MolDS_wrappers::Blas::GetInstance()->Dcopy(OrbitalType_end*OrbitalType_end,
3822- &rotMat1stDerivs[0][0][c], CartesianType_end,
3823- &tmpRotMat1stDeriv[0][0], incrementOne);
3824- alpha = cartesian[c]/R;
3825- beta = 0.0;
3826- MolDS_wrappers::Blas::GetInstance()->Dgemmm(isColumnMajorRotatingMatrix,
3827- isColumnMajorDiaOverlapAOs,
3828- !isColumnMajorRotatingMatrix,
3829- OrbitalType_end, OrbitalType_end, OrbitalType_end, OrbitalType_end,
3830- alpha,
3831- rotMat,
3832- diaOverlapAOs1stDerivInDiaFrame,
3833- rotMat,
3834- beta,
3835- rotatedDiatomicOverlap,
3836- tmpMatrix);
3837- alpha = 1.0;
3838- beta = 1.0;
3839- MolDS_wrappers::Blas::GetInstance()->Dgemmm(isColumnMajorRotatingMatrix,
3840- isColumnMajorDiaOverlapAOs,
3841- !isColumnMajorRotatingMatrix,
3842- OrbitalType_end, OrbitalType_end, OrbitalType_end, OrbitalType_end,
3843- alpha,
3844- tmpRotMat1stDeriv,
3845- diaOverlapAOsInDiaFrame,
3846- rotMat,
3847- beta,
3848- rotatedDiatomicOverlap,
3849- tmpMatrix);
3850- MolDS_wrappers::Blas::GetInstance()->Dgemmm(isColumnMajorRotatingMatrix,
3851- isColumnMajorDiaOverlapAOs,
3852- !isColumnMajorRotatingMatrix,
3853- OrbitalType_end, OrbitalType_end, OrbitalType_end, OrbitalType_end,
3854- alpha,
3855- rotMat,
3856- diaOverlapAOsInDiaFrame,
3857- tmpRotMat1stDeriv,
3858- beta,
3859- rotatedDiatomicOverlap,
3860- tmpMatrix);
3861- MolDS_wrappers::Blas::GetInstance()->Dcopy(OrbitalType_end*OrbitalType_end,
3862- &rotatedDiatomicOverlap[0][0], incrementOne,
3863- &diatomicOverlapAOs1stDerivs[0][0][c], CartesianType_end);
3864- }
3865-
3866- /*
3867- // rotate (slow)
3868- for(int i=0; i<OrbitalType_end; i++){
3869- for(int j=0; j<OrbitalType_end; j++){
3870- for(int c=0; c<CartesianType_end; c++){
3871- diatomicOverlapAOs1stDerivs[i][j][c] = 0.0;
3872-
3873- double temp1 = 0.0;
3874- double temp2 = 0.0;
3875- double temp3 = 0.0;
3876- for(int k=0; k<OrbitalType_end; k++){
3877- for(int l=0; l<OrbitalType_end; l++){
3878- temp1 += rotMat[i][k]
3879- *rotMat[j][l]
3880- *(cartesian[c]/R)
3881- *diaOverlapAOs1stDerivInDiaFrame[k][l];
3882- temp2 += rotMat1stDerivs[i][k][c]
3883- *rotMat[j][l]
3884- *diaOverlapAOsInDiaFrame[k][l];
3885- temp3 += rotMat[i][k]
3886- *rotMat1stDerivs[j][l][c]
3887- *diaOverlapAOsInDiaFrame[k][l];
3888- }
3857+ /*
3858+ // rotate (slow)
3859+ for(int i=0; i<OrbitalType_end; i++){
3860+ for(int j=0; j<OrbitalType_end; j++){
3861+ for(int c=0; c<CartesianType_end; c++){
3862+ diatomicOverlapAOs1stDerivs[i][j][c] = 0.0;
3863+
3864+ double temp1 = 0.0;
3865+ double temp2 = 0.0;
3866+ double temp3 = 0.0;
3867+ for(int k=0; k<OrbitalType_end; k++){
3868+ for(int l=0; l<OrbitalType_end; l++){
3869+ temp1 += tmpRotMat[i][k]
3870+ *tmpRotMat[j][l]
3871+ *(cartesian[c]/R)
3872+ *tmpDiaOverlapAOs1stDerivInDiaFrame[k][l];
3873+ temp2 += tmpRotMat1stDerivs[i][k][c]
3874+ *tmpRotMat[j][l]
3875+ *tmpDiaOverlapAOsInDiaFrame[k][l];
3876+ temp3 += tmpRotMat[i][k]
3877+ *tmpRotMat1stDerivs[j][l][c]
3878+ *tmpDiaOverlapAOsInDiaFrame[k][l];
38893879 }
3890- diatomicOverlapAOs1stDerivs[i][j][c] = temp1 + temp2 + temp3;
38913880 }
3881+ diatomicOverlapAOs1stDerivs[i][j][c] = temp1 + temp2 + temp3;
38923882 }
38933883 }
3894- */
3895-
3896- }
3897- catch(MolDSException ex){
3898- this->FreeDiatomicOverlapAOs1stDeriTemps(&diaOverlapAOsInDiaFrame,
3899- &diaOverlapAOs1stDerivInDiaFrame,
3900- &rotMat,
3901- &rotMat1stDerivs,
3902- &rotatedDiatomicOverlap,
3903- &tmpRotMat1stDeriv,
3904- &tmpMatrix);
3905- throw ex;
39063884 }
3907- // free
3908- this->FreeDiatomicOverlapAOs1stDeriTemps(&diaOverlapAOsInDiaFrame,
3909- &diaOverlapAOs1stDerivInDiaFrame,
3910- &rotMat,
3911- &rotMat1stDerivs,
3912- &rotatedDiatomicOverlap,
3913- &tmpRotMat1stDeriv,
3914- &tmpMatrix);
3885+ */
39153886 }
39163887
39173888 void Cndo2::CalcDiatomicOverlapAOs1stDerivatives(double*** diatomicOverlapAOs1stDerivs,
3918- int indexAtomA,
3919- int indexAtomB) const{
3889+ double** tmpDiaOverlapAOsInDiaFrame,
3890+ double** tmpDiaOverlapAOs1stDerivInDiaFrame,
3891+ double** tmpRotMat,
3892+ double** tmpRotMat1stDeriv,
3893+ double*** tmpRotMat1stDerivs,
3894+ double** tmpRotatedDiatomicOverlap,
3895+ double** tmpMatrix,
3896+ int indexAtomA,
3897+ int indexAtomB) const{
39203898 this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs,
3899+ tmpDiaOverlapAOsInDiaFrame,
3900+ tmpDiaOverlapAOs1stDerivInDiaFrame,
3901+ tmpRotMat,
3902+ tmpRotMat1stDeriv,
3903+ tmpRotMat1stDerivs,
3904+ tmpRotatedDiatomicOverlap,
3905+ tmpMatrix,
39213906 *this->molecule->GetAtom(indexAtomA),
39223907 *this->molecule->GetAtom(indexAtomB));
39233908 }
@@ -4099,22 +4084,6 @@ double Cndo2::Get2ndDerivativeElementFromDistanceDerivatives(double firstDistanc
40994084 return value;
41004085 }
41014086
4102-void Cndo2::MallocDiatomicOverlapAOs1stDeriTemps(double*** diaOverlapAOsInDiaFrame,
4103- double*** diaOverlapAOs1stDerivInDiaFrame,
4104- double*** rotMat,
4105- double**** rotMat1stDerivs,
4106- double*** rotatedDiatomicOverlap,
4107- double*** tmpRotMat1stDeriv,
4108- double*** tmpMatrix) const{
4109- MallocerFreer::GetInstance()->Malloc<double>(diaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
4110- MallocerFreer::GetInstance()->Malloc<double>(diaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
4111- MallocerFreer::GetInstance()->Malloc<double>(rotMat, OrbitalType_end, OrbitalType_end);
4112- MallocerFreer::GetInstance()->Malloc<double>(rotMat1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
4113- MallocerFreer::GetInstance()->Malloc<double>(rotatedDiatomicOverlap, OrbitalType_end, OrbitalType_end);
4114- MallocerFreer::GetInstance()->Malloc<double>(tmpMatrix, OrbitalType_end, OrbitalType_end);
4115- MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat1stDeriv, OrbitalType_end, OrbitalType_end);
4116-}
4117-
41184087 void Cndo2::MallocDiatomicOverlapAOs2ndDeriTemps(double*** diaOverlapAOsInDiaFrame,
41194088 double*** diaOverlapAOs1stDerivInDiaFrame,
41204089 double*** diaOverlapAOs2ndDerivInDiaFrame,
@@ -4133,22 +4102,6 @@ void Cndo2::MallocDiatomicOverlapAOs2ndDeriTemps(double*** diaOverlapAOsInDiaFra
41334102 MallocerFreer::GetInstance()->Malloc<double>(tempDiaOverlapAOs2ndDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end, CartesianType_end);
41344103 }
41354104
4136-void Cndo2::FreeDiatomicOverlapAOs1stDeriTemps(double*** diaOverlapAOsInDiaFrame,
4137- double*** diaOverlapAOs1stDerivInDiaFrame,
4138- double*** rotMat,
4139- double**** rotMat1stDerivs,
4140- double*** rotatedDiatomicOverlap,
4141- double*** tmpRotMat1stDeriv,
4142- double*** tmpMatrix) const{
4143- MallocerFreer::GetInstance()->Free<double>(diaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
4144- MallocerFreer::GetInstance()->Free<double>(diaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
4145- MallocerFreer::GetInstance()->Free<double>(rotMat, OrbitalType_end, OrbitalType_end);
4146- MallocerFreer::GetInstance()->Free<double>(rotMat1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
4147- MallocerFreer::GetInstance()->Free<double>(rotatedDiatomicOverlap, OrbitalType_end, OrbitalType_end);
4148- MallocerFreer::GetInstance()->Free<double>(tmpMatrix, OrbitalType_end, OrbitalType_end);
4149- MallocerFreer::GetInstance()->Free<double>(tmpRotMat1stDeriv, OrbitalType_end, OrbitalType_end);
4150-}
4151-
41524105 void Cndo2::FreeDiatomicOverlapAOs2ndDeriTemps(double*** diaOverlapAOsInDiaFrame,
41534106 double*** diaOverlapAOs1stDerivInDiaFrame,
41544107 double*** diaOverlapAOs2ndDerivInDiaFrame,
@@ -5250,54 +5203,46 @@ void Cndo2::CalcRotatingMatrix(double** rotatingMatrix,
52505203
52515204 // rotating matrix for d-function
52525205 // dMatrix is (37) in J. Mol. Strct. 419, 19(1997) (ref. [BFB_1997])
5253- double** dMatrix = NULL;
5254- try{
5255- MallocerFreer::GetInstance()->Malloc<double>(&dMatrix, OrbitalType_end, OrbitalType_end);
5256- dMatrix[dzz][dzz] = 0.5*(3.0*pow(cos(beta),2.0) - 1.0);
5257- dMatrix[dxxyy][dxxyy] = pow(cos(0.5*beta),4.0);
5258- dMatrix[dzx][dzx] = (2.0*cos(beta)-1.0)*pow(cos(0.5*beta),2.0);
5259- dMatrix[dxxyy][dzx] = -2.0*sin(0.5*beta)*pow(cos(0.5*beta),3.0);
5260- dMatrix[dxxyy][dzz] = sqrt(6.0)*pow(sin(0.5*beta),2.0)*pow(cos(0.5*beta),2.0);
5261- dMatrix[dxxyy][dyz] = -2.0*pow(sin(0.5*beta),3.0)*pow(cos(0.5*beta),1.0);
5262- dMatrix[dxxyy][dxy] = pow(sin(0.5*beta),4.0);
5263- dMatrix[dzx][dzz] = -sqrt(6.0)*cos(beta)*cos(0.5*beta)*sin(0.5*beta);
5264- dMatrix[dzx][dyz] = (2.0*cos(beta)+1.0)*pow(sin(0.5*beta),2.0);
5265-
5266- rotatingMatrix[dxy][dxy] = cos(2.0*alpha)* (dMatrix[dxxyy][dxxyy] - dMatrix[dxxyy][dxy]);
5267- rotatingMatrix[dxy][dyz] = cos(2.0*alpha)* (-1.0*dMatrix[dxxyy][dzx] - dMatrix[dxxyy][dyz]);
5268- rotatingMatrix[dxy][dzz] = sqrt(2.0)*sin(2.0*alpha)* dMatrix[dxxyy][dzz];
5269- rotatingMatrix[dxy][dzx] = sin(2.0*alpha)* (-1.0*dMatrix[dxxyy][dzx] + dMatrix[dxxyy][dyz]);
5270- rotatingMatrix[dxy][dxxyy] = sin(2.0*alpha)* (dMatrix[dxxyy][dxxyy] + dMatrix[dxxyy][dxy]);
5271-
5272- rotatingMatrix[dyz][dxy] = cos(alpha)* (dMatrix[dxxyy][dzx] + dMatrix[dxxyy][dyz]);
5273- rotatingMatrix[dyz][dyz] = cos(alpha)* (dMatrix[dzx][dzx] + dMatrix[dzx][dyz]);
5274- rotatingMatrix[dyz][dzz] = -1.0*sqrt(2.0)*sin(alpha)* dMatrix[dzx][dzz];
5275- rotatingMatrix[dyz][dzx] = sin(alpha)* (dMatrix[dzx][dzx] - dMatrix[dzx][dyz]);
5276- rotatingMatrix[dyz][dxxyy] = sin(alpha)* (dMatrix[dxxyy][dzx] - dMatrix[dxxyy][dyz]);
5277-
5278- rotatingMatrix[dzz][dxy] = 0.0;
5279- rotatingMatrix[dzz][dyz] = 0.0;
5280- rotatingMatrix[dzz][dzz] = dMatrix[dzz][dzz];
5281- rotatingMatrix[dzz][dzx] = sqrt(2.0)*dMatrix[dzx][dzz];
5282- rotatingMatrix[dzz][dxxyy] = sqrt(2.0)*dMatrix[dxxyy][dzz];
5283-
5284- rotatingMatrix[dzx][dxy] = -1.0*sin(alpha)* (dMatrix[dxxyy][dzx] + dMatrix[dxxyy][dyz]);
5285- rotatingMatrix[dzx][dyz] = -1.0*sin(alpha)* (dMatrix[dzx][dzx] + dMatrix[dzx][dyz]);
5286- rotatingMatrix[dzx][dzz] = -1.0*sqrt(2.0)*cos(alpha)* dMatrix[dzx][dzz];
5287- rotatingMatrix[dzx][dzx] = cos(alpha)* (dMatrix[dzx][dzx] - dMatrix[dzx][dyz]);
5288- rotatingMatrix[dzx][dxxyy] = cos(alpha)* (dMatrix[dxxyy][dzx] - dMatrix[dxxyy][dyz]);
5289-
5290- rotatingMatrix[dxxyy][dxy] = -1.0*sin(2.0*alpha)* (dMatrix[dxxyy][dxxyy] - dMatrix[dxxyy][dxy]);
5291- rotatingMatrix[dxxyy][dyz] = -1.0*sin(2.0*alpha)* (-1.0*dMatrix[dxxyy][dzx] - dMatrix[dxxyy][dyz]);
5292- rotatingMatrix[dxxyy][dzz] = sqrt(2.0)*cos(2.0*alpha)*dMatrix[dxxyy][dzz];
5293- rotatingMatrix[dxxyy][dzx] = cos(2.0*alpha)* (-1.0*dMatrix[dxxyy][dzx] + dMatrix[dxxyy][dyz]);
5294- rotatingMatrix[dxxyy][dxxyy] = cos(2.0*alpha)* (dMatrix[dxxyy][dxxyy] + dMatrix[dxxyy][dxy]);
5295- }
5296- catch(MolDSException ex){
5297- MallocerFreer::GetInstance()->Free<double>(&dMatrix, OrbitalType_end, OrbitalType_end);
5298- throw ex;
5299- }
5300- MallocerFreer::GetInstance()->Free<double>(&dMatrix, OrbitalType_end, OrbitalType_end);
5206+ double dMatrix[OrbitalType_end][OrbitalType_end];
5207+ dMatrix[dzz][dzz] = 0.5*(3.0*pow(cos(beta),2.0) - 1.0);
5208+ dMatrix[dxxyy][dxxyy] = pow(cos(0.5*beta),4.0);
5209+ dMatrix[dzx][dzx] = (2.0*cos(beta)-1.0)*pow(cos(0.5*beta),2.0);
5210+ dMatrix[dxxyy][dzx] = -2.0*sin(0.5*beta)*pow(cos(0.5*beta),3.0);
5211+ dMatrix[dxxyy][dzz] = sqrt(6.0)*pow(sin(0.5*beta),2.0)*pow(cos(0.5*beta),2.0);
5212+ dMatrix[dxxyy][dyz] = -2.0*pow(sin(0.5*beta),3.0)*pow(cos(0.5*beta),1.0);
5213+ dMatrix[dxxyy][dxy] = pow(sin(0.5*beta),4.0);
5214+ dMatrix[dzx][dzz] = -sqrt(6.0)*cos(beta)*cos(0.5*beta)*sin(0.5*beta);
5215+ dMatrix[dzx][dyz] = (2.0*cos(beta)+1.0)*pow(sin(0.5*beta),2.0);
5216+
5217+ rotatingMatrix[dxy][dxy] = cos(2.0*alpha)* (dMatrix[dxxyy][dxxyy] - dMatrix[dxxyy][dxy]);
5218+ rotatingMatrix[dxy][dyz] = cos(2.0*alpha)* (-1.0*dMatrix[dxxyy][dzx] - dMatrix[dxxyy][dyz]);
5219+ rotatingMatrix[dxy][dzz] = sqrt(2.0)*sin(2.0*alpha)* dMatrix[dxxyy][dzz];
5220+ rotatingMatrix[dxy][dzx] = sin(2.0*alpha)* (-1.0*dMatrix[dxxyy][dzx] + dMatrix[dxxyy][dyz]);
5221+ rotatingMatrix[dxy][dxxyy] = sin(2.0*alpha)* (dMatrix[dxxyy][dxxyy] + dMatrix[dxxyy][dxy]);
5222+
5223+ rotatingMatrix[dyz][dxy] = cos(alpha)* (dMatrix[dxxyy][dzx] + dMatrix[dxxyy][dyz]);
5224+ rotatingMatrix[dyz][dyz] = cos(alpha)* (dMatrix[dzx][dzx] + dMatrix[dzx][dyz]);
5225+ rotatingMatrix[dyz][dzz] = -1.0*sqrt(2.0)*sin(alpha)* dMatrix[dzx][dzz];
5226+ rotatingMatrix[dyz][dzx] = sin(alpha)* (dMatrix[dzx][dzx] - dMatrix[dzx][dyz]);
5227+ rotatingMatrix[dyz][dxxyy] = sin(alpha)* (dMatrix[dxxyy][dzx] - dMatrix[dxxyy][dyz]);
5228+
5229+ rotatingMatrix[dzz][dxy] = 0.0;
5230+ rotatingMatrix[dzz][dyz] = 0.0;
5231+ rotatingMatrix[dzz][dzz] = dMatrix[dzz][dzz];
5232+ rotatingMatrix[dzz][dzx] = sqrt(2.0)*dMatrix[dzx][dzz];
5233+ rotatingMatrix[dzz][dxxyy] = sqrt(2.0)*dMatrix[dxxyy][dzz];
5234+
5235+ rotatingMatrix[dzx][dxy] = -1.0*sin(alpha)* (dMatrix[dxxyy][dzx] + dMatrix[dxxyy][dyz]);
5236+ rotatingMatrix[dzx][dyz] = -1.0*sin(alpha)* (dMatrix[dzx][dzx] + dMatrix[dzx][dyz]);
5237+ rotatingMatrix[dzx][dzz] = -1.0*sqrt(2.0)*cos(alpha)* dMatrix[dzx][dzz];
5238+ rotatingMatrix[dzx][dzx] = cos(alpha)* (dMatrix[dzx][dzx] - dMatrix[dzx][dyz]);
5239+ rotatingMatrix[dzx][dxxyy] = cos(alpha)* (dMatrix[dxxyy][dzx] - dMatrix[dxxyy][dyz]);
5240+
5241+ rotatingMatrix[dxxyy][dxy] = -1.0*sin(2.0*alpha)* (dMatrix[dxxyy][dxxyy] - dMatrix[dxxyy][dxy]);
5242+ rotatingMatrix[dxxyy][dyz] = -1.0*sin(2.0*alpha)* (-1.0*dMatrix[dxxyy][dzx] - dMatrix[dxxyy][dyz]);
5243+ rotatingMatrix[dxxyy][dzz] = sqrt(2.0)*cos(2.0*alpha)*dMatrix[dxxyy][dzz];
5244+ rotatingMatrix[dxxyy][dzx] = cos(2.0*alpha)* (-1.0*dMatrix[dxxyy][dzx] + dMatrix[dxxyy][dyz]);
5245+ rotatingMatrix[dxxyy][dxxyy] = cos(2.0*alpha)* (dMatrix[dxxyy][dxxyy] + dMatrix[dxxyy][dxy]);
53015246 }
53025247
53035248 // First derivative of rotating matirx.
--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -195,9 +195,23 @@ protected:
195195 const MolDS_base_atoms::Atom& atomA,
196196 const MolDS_base_atoms::Atom& atomB) const;
197197 void CalcDiatomicOverlapAOs1stDerivatives(double*** diatomicOverlapAOs1stDerivs,
198+ double** tmpDiaOverlapAOsInDiaFrame,
199+ double** tmpDiaOverlapAOs1stDerivInDiaFrame,
200+ double** tmpRotMat,
201+ double** tmpRotMat1stDeriv,
202+ double*** tmpRotMat1stDerivs,
203+ double** tmpRotatedDiatomicOverlap,
204+ double** tmpMatrix,
198205 const MolDS_base_atoms::Atom& atomA,
199206 const MolDS_base_atoms::Atom& atomB) const;
200207 void CalcDiatomicOverlapAOs1stDerivatives(double*** diatomicOverlapAOs1stDerivs,
208+ double** tmpDiaOverlapAOsInDiaFrame,
209+ double** tmpDiaOverlapAOs1stDerivInDiaFrame,
210+ double** tmpRotMat,
211+ double** tmpRotMat1stDeriv,
212+ double*** tmpRotMat1stDerivs,
213+ double** tmpRotatedDiatomicOverlap,
214+ double** tmpMatrix,
201215 int indexAtomA,
202216 int indexAtomB) const;
203217 void CalcDiatomicOverlapAOs2ndDerivatives(double**** overlapAOs2ndDeri,
@@ -482,13 +496,6 @@ private:
482496 double**** diisStoredErrorVect,
483497 double*** diisErrorProducts,
484498 double** diisErrorCoefficients);
485- void MallocDiatomicOverlapAOs1stDeriTemps(double*** diaOverlapAOsInDiaFrame,
486- double*** diaOverlapAOs1stDerivInDiaFrame,
487- double*** rotMat,
488- double**** rotMat1stDerivs,
489- double*** rotatedDiatomicOverlap,
490- double*** tmpRotMat1stDeriv,
491- double*** tmpMatrix) const;
492499 void MallocDiatomicOverlapAOs2ndDeriTemps(double*** diaOverlapAOsInDiaFrame,
493500 double*** diaOverlapAOs1stDerivInDiaFrame,
494501 double*** diaOverlapAOs2ndDerivInDiaFrame,
@@ -497,13 +504,6 @@ private:
497504 double***** rotMat2ndDerivs,
498505 double**** tempDiaOverlapAOs1stDerivs,
499506 double***** tempDiaOverlapAOs2ndDerivs) const;
500- void FreeDiatomicOverlapAOs1stDeriTemps(double*** diaOverlapAOsInDiaFrame,
501- double*** diaOverlapAOs1stDerivInDiaFrame,
502- double*** rotMat,
503- double**** rotMat1stDerivs,
504- double*** rotatedDiatomicOverlap,
505- double*** tmpRotMat1stDeriv,
506- double*** tmpMatrix) const;
507507 void FreeDiatomicOverlapAOs2ndDeriTemps(double*** diaOverlapAOsInDiaFrame,
508508 double*** diaOverlapAOs1stDerivInDiaFrame,
509509 double*** diaOverlapAOs2ndDerivInDiaFrame,
--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -974,10 +974,15 @@ void Mndo::MallocTempMatricesEachThreadCalcHessianSCF(double***** diatomicOve
974974 double******* diatomicTwoElecTwoCore1stDerivs,
975975 double******** diatomicTwoElecTwoCore2ndDerivs,
976976 double*** tmpRotMat,
977+ double*** tmpRotMat1stDeriv,
977978 double**** tmpRotMat1stDerivs,
978979 double***** tmpRotMat2ndDerivs,
979980 double***** tmpDiatomicTwoElecTwoCore,
980- double****** tmpDiatomicTwoElecTwoCore1stDerivs) const{
981+ double****** tmpDiatomicTwoElecTwoCore1stDerivs,
982+ double*** tmpDiaOverlapAOsInDiaFrame,
983+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
984+ double*** tmpRotatedDiatomicOverlap,
985+ double*** tmpMatrix) const{
981986 MallocerFreer::GetInstance()->Malloc<double>(diatomicOverlapAOs1stDerivs,
982987 this->molecule->GetNumberAtoms(),
983988 OrbitalType_end,
@@ -1007,6 +1012,9 @@ void Mndo::MallocTempMatricesEachThreadCalcHessianSCF(double***** diatomicOve
10071012 MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat,
10081013 OrbitalType_end,
10091014 OrbitalType_end);
1015+ MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat1stDeriv,
1016+ OrbitalType_end,
1017+ OrbitalType_end);
10101018 MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat1stDerivs,
10111019 OrbitalType_end,
10121020 OrbitalType_end,
@@ -1027,6 +1035,17 @@ void Mndo::MallocTempMatricesEachThreadCalcHessianSCF(double***** diatomicOve
10271035 dxy,
10281036 dxy,
10291037 CartesianType_end);
1038+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiaOverlapAOsInDiaFrame,
1039+ OrbitalType_end,
1040+ OrbitalType_end);
1041+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiaOverlapAOs1stDerivInDiaFrame,
1042+ OrbitalType_end,
1043+ OrbitalType_end);
1044+ MallocerFreer::GetInstance()->Malloc<double>(tmpRotatedDiatomicOverlap,
1045+ OrbitalType_end, OrbitalType_end);
1046+ MallocerFreer::GetInstance()->Malloc<double>(tmpMatrix,
1047+ OrbitalType_end,
1048+ OrbitalType_end);
10301049 }
10311050
10321051 void Mndo::FreeTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverlapAOs1stDerivs,
@@ -1034,10 +1053,15 @@ void Mndo::FreeTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverl
10341053 double******* diatomicTwoElecTwoCore1stDerivs,
10351054 double******** diatomicTwoElecTwoCore2ndDerivs,
10361055 double*** tmpRotMat,
1056+ double*** tmpRotMat1stDeriv,
10371057 double**** tmpRotMat1stDerivs,
10381058 double***** tmpRotMat2ndDerivs,
10391059 double***** tmpDiatomicTwoElecTwoCore,
1040- double****** tmpDiatomicTwoElecTwoCore1stDerivs) const{
1060+ double****** tmpDiatomicTwoElecTwoCore1stDerivs,
1061+ double*** tmpDiaOverlapAOsInDiaFrame,
1062+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
1063+ double*** tmpRotatedDiatomicOverlap,
1064+ double*** tmpMatrix) const{
10411065 MallocerFreer::GetInstance()->Free<double>(diatomicOverlapAOs1stDerivs,
10421066 this->molecule->GetNumberAtoms(),
10431067 OrbitalType_end,
@@ -1067,6 +1091,9 @@ void Mndo::FreeTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverl
10671091 MallocerFreer::GetInstance()->Free<double>(tmpRotMat,
10681092 OrbitalType_end,
10691093 OrbitalType_end);
1094+ MallocerFreer::GetInstance()->Free<double>(tmpRotMat1stDeriv,
1095+ OrbitalType_end,
1096+ OrbitalType_end);
10701097 MallocerFreer::GetInstance()->Free<double>(tmpRotMat1stDerivs,
10711098 OrbitalType_end,
10721099 OrbitalType_end,
@@ -1087,6 +1114,18 @@ void Mndo::FreeTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverl
10871114 dxy,
10881115 dxy,
10891116 CartesianType_end);
1117+ MallocerFreer::GetInstance()->Free<double>(tmpDiaOverlapAOsInDiaFrame,
1118+ OrbitalType_end,
1119+ OrbitalType_end);
1120+ MallocerFreer::GetInstance()->Free<double>(tmpDiaOverlapAOs1stDerivInDiaFrame,
1121+ OrbitalType_end,
1122+ OrbitalType_end);
1123+ MallocerFreer::GetInstance()->Free<double>(tmpRotatedDiatomicOverlap,
1124+ OrbitalType_end,
1125+ OrbitalType_end);
1126+ MallocerFreer::GetInstance()->Free<double>(tmpMatrix,
1127+ OrbitalType_end,
1128+ OrbitalType_end);
10901129 }
10911130
10921131 // mu and nu is included in atomA' AO.
@@ -1637,16 +1676,26 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
16371676 double**** tmpRotMat2ndDerivs = NULL;
16381677 double**** tmpDiatomicTwoElecTwoCore = NULL;
16391678 double***** tmpDiatomicTwoElecTwoCore1stDerivs = NULL;
1679+ double** tmpDiaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
1680+ double** tmpDiaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
1681+ double** tmpRotMat1stDeriv = NULL;
1682+ double** tmpRotatedDiatomicOverlap = NULL;
1683+ double** tmpMatrix = NULL;
16401684 try{
16411685 this->MallocTempMatricesEachThreadCalcHessianSCF(&diatomicOverlapAOs1stDerivs,
16421686 &diatomicOverlapAOs2ndDerivs,
16431687 &diatomicTwoElecTwoCore1stDerivs,
16441688 &diatomicTwoElecTwoCore2ndDerivs,
16451689 &tmpRotMat,
1690+ &tmpRotMat1stDeriv,
16461691 &tmpRotMat1stDerivs,
16471692 &tmpRotMat2ndDerivs,
16481693 &tmpDiatomicTwoElecTwoCore,
1649- &tmpDiatomicTwoElecTwoCore1stDerivs);
1694+ &tmpDiatomicTwoElecTwoCore1stDerivs,
1695+ &tmpDiaOverlapAOsInDiaFrame,
1696+ &tmpDiaOverlapAOs1stDerivInDiaFrame,
1697+ &tmpRotatedDiatomicOverlap,
1698+ &tmpMatrix);
16501699 #pragma omp for schedule(auto)
16511700 for(int indexAtomA=0; indexAtomA<this->molecule->GetNumberAtoms(); indexAtomA++){
16521701 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
@@ -1658,6 +1707,13 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
16581707 for(int indexAtomB=0; indexAtomB<this->molecule->GetNumberAtoms(); indexAtomB++){
16591708 if(indexAtomA != indexAtomB){
16601709 this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs[indexAtomB],
1710+ tmpDiaOverlapAOsInDiaFrame,
1711+ tmpDiaOverlapAOs1stDerivInDiaFrame,
1712+ tmpRotMat,
1713+ tmpRotMat1stDeriv,
1714+ tmpRotMat1stDerivs,
1715+ tmpRotatedDiatomicOverlap,
1716+ tmpMatrix,
16611717 indexAtomA,
16621718 indexAtomB);
16631719 this->CalcDiatomicOverlapAOs2ndDerivatives(diatomicOverlapAOs2ndDerivs[indexAtomB],
@@ -1735,10 +1791,15 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
17351791 &diatomicTwoElecTwoCore1stDerivs,
17361792 &diatomicTwoElecTwoCore2ndDerivs,
17371793 &tmpRotMat,
1794+ &tmpRotMat1stDeriv,
17381795 &tmpRotMat1stDerivs,
17391796 &tmpRotMat2ndDerivs,
17401797 &tmpDiatomicTwoElecTwoCore,
1741- &tmpDiatomicTwoElecTwoCore1stDerivs);
1798+ &tmpDiatomicTwoElecTwoCore1stDerivs,
1799+ &tmpDiaOverlapAOsInDiaFrame,
1800+ &tmpDiaOverlapAOs1stDerivInDiaFrame,
1801+ &tmpRotatedDiatomicOverlap,
1802+ &tmpMatrix);
17421803 }// end of omp-region
17431804 // Exception throwing for omp-region
17441805 if(!ompErrors.str().empty()){
@@ -1925,12 +1986,22 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
19251986 double*** tmpRotMat1stDerivs = NULL;
19261987 double**** tmpDiatomicTwoElecTwoCore = NULL;
19271988
1989+ double** tmpDiaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
1990+ double** tmpDiaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
1991+ double** tmpRotMat1stDeriv = NULL;
1992+ double** tmpRotatedDiatomicOverlap = NULL;
1993+ double** tmpMatrix = NULL;
19281994 try{
19291995 this->MallocTempMatricesStaticFirstOrderFock(&diatomicTwoElecTwoCore1stDerivs,
19301996 &diatomicOverlapAOs1stDerivs,
19311997 &tmpRotMat,
19321998 &tmpRotMat1stDerivs,
19331999 &tmpDiatomicTwoElecTwoCore);
2000+ MallocerFreer::GetInstance()->Malloc<double>(&tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
2001+ MallocerFreer::GetInstance()->Malloc<double>(&tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
2002+ MallocerFreer::GetInstance()->Malloc<double>(&tmpRotMat1stDeriv, OrbitalType_end, OrbitalType_end);
2003+ MallocerFreer::GetInstance()->Malloc<double>(&tmpRotatedDiatomicOverlap, OrbitalType_end, OrbitalType_end);
2004+ MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix, OrbitalType_end, OrbitalType_end);
19342005 const Atom& atomA = *molecule->GetAtom(indexAtomA);
19352006 int firstAOIndexA = atomA.GetFirstAOIndex();
19362007 int lastAOIndexA = atomA.GetLastAOIndex();
@@ -1949,7 +2020,16 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
19492020 tmpDiatomicTwoElecTwoCore,
19502021 indexAtomA, indexAtomB);
19512022 // calc. first derivative of overlapAOs.
1952- this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs, atomA, atomB);
2023+ this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs,
2024+ tmpDiaOverlapAOsInDiaFrame,
2025+ tmpDiaOverlapAOs1stDerivInDiaFrame,
2026+ tmpRotMat,
2027+ tmpRotMat1stDeriv,
2028+ tmpRotMat1stDerivs,
2029+ tmpRotatedDiatomicOverlap,
2030+ tmpMatrix,
2031+ atomA,
2032+ atomB);
19532033
19542034 for(int i=0; i<nonRedundantQIndeces.size()+redundantQIndeces.size();i++){
19552035 int moI=0, moJ=0;
@@ -2045,6 +2125,13 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
20452125 &tmpRotMat,
20462126 &tmpRotMat1stDerivs,
20472127 &tmpDiatomicTwoElecTwoCore);
2128+ MallocerFreer::GetInstance()->Free<double>(&tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
2129+ MallocerFreer::GetInstance()->Free<double>(&tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
2130+ //MallocerFreer::GetInstance()->Free<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
2131+ MallocerFreer::GetInstance()->Free<double>(&tmpRotMat1stDeriv, OrbitalType_end, OrbitalType_end);
2132+ //MallocerFreer::GetInstance()->Free<double>(&tmpRotMat1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
2133+ MallocerFreer::GetInstance()->Free<double>(&tmpRotatedDiatomicOverlap, OrbitalType_end, OrbitalType_end);
2134+ MallocerFreer::GetInstance()->Free<double>(&tmpMatrix, OrbitalType_end, OrbitalType_end);
20482135 throw ex;
20492136 }
20502137 this->FreeTempMatricesStaticFirstOrderFock(&diatomicTwoElecTwoCore1stDerivs,
@@ -2052,6 +2139,13 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
20522139 &tmpRotMat,
20532140 &tmpRotMat1stDerivs,
20542141 &tmpDiatomicTwoElecTwoCore);
2142+ MallocerFreer::GetInstance()->Free<double>(&tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
2143+ MallocerFreer::GetInstance()->Free<double>(&tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
2144+ //MallocerFreer::GetInstance()->Free<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
2145+ MallocerFreer::GetInstance()->Free<double>(&tmpRotMat1stDeriv, OrbitalType_end, OrbitalType_end);
2146+ //MallocerFreer::GetInstance()->Free<double>(&tmpRotMat1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
2147+ MallocerFreer::GetInstance()->Free<double>(&tmpRotatedDiatomicOverlap, OrbitalType_end, OrbitalType_end);
2148+ MallocerFreer::GetInstance()->Free<double>(&tmpMatrix, OrbitalType_end, OrbitalType_end);
20552149
20562150 /*
20572151 printf("staticFirstOrderFock(atomA:%d axis:%s)\n",indexAtomA,CartesianTypeStr(axisA));
@@ -2413,16 +2507,27 @@ void Mndo::CalcForce(const vector<int>& elecStates){
24132507 stringstream ompErrors;
24142508 #pragma omp parallel
24152509 {
2416- double***** diatomicTwoElecTwoCore1stDerivs = NULL;
24172510 double*** diatomicOverlapAOs1stDerivs = NULL;
2511+ double***** diatomicTwoElecTwoCore1stDerivs = NULL;
24182512 double** tmpRotMat = NULL;
24192513 double*** tmpRotMat1stDerivs = NULL;
24202514 double**** tmpDiatomicTwoElecTwoCore = NULL;
2515+
2516+ double** tmpDiaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
2517+ double** tmpDiaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
2518+ double** tmpRotMat1stDeriv = NULL;
2519+ double** tmpRotatedDiatomicOverlap = NULL;
2520+ double** tmpMatrix = NULL;
24212521 try{
24222522 this->MallocTempMatricesCalcForce(&diatomicOverlapAOs1stDerivs,
24232523 &diatomicTwoElecTwoCore1stDerivs,
2524+ &tmpDiaOverlapAOsInDiaFrame,
2525+ &tmpDiaOverlapAOs1stDerivInDiaFrame,
24242526 &tmpRotMat,
2527+ &tmpRotMat1stDeriv,
24252528 &tmpRotMat1stDerivs,
2529+ &tmpRotatedDiatomicOverlap,
2530+ &tmpMatrix,
24262531 &tmpDiatomicTwoElecTwoCore);
24272532
24282533 #pragma omp for schedule(auto)
@@ -2433,7 +2538,16 @@ void Mndo::CalcForce(const vector<int>& elecStates){
24332538 int lastAOIndexB = atomB.GetLastAOIndex();
24342539
24352540 // calc. first derivative of overlapAOs.
2436- this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs, atomA, atomB);
2541+ this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs,
2542+ tmpDiaOverlapAOsInDiaFrame,
2543+ tmpDiaOverlapAOs1stDerivInDiaFrame,
2544+ tmpRotMat,
2545+ tmpRotMat1stDeriv,
2546+ tmpRotMat1stDerivs,
2547+ tmpRotatedDiatomicOverlap,
2548+ tmpMatrix,
2549+ atomA,
2550+ atomB);
24372551 // calc. first derivative of two elec two core interaction
24382552 this->CalcDiatomicTwoElecTwoCore1stDerivatives(diatomicTwoElecTwoCore1stDerivs,
24392553 tmpRotMat,
@@ -2548,8 +2662,13 @@ void Mndo::CalcForce(const vector<int>& elecStates){
25482662 }
25492663 this->FreeTempMatricesCalcForce(&diatomicOverlapAOs1stDerivs,
25502664 &diatomicTwoElecTwoCore1stDerivs,
2665+ &tmpDiaOverlapAOsInDiaFrame,
2666+ &tmpDiaOverlapAOs1stDerivInDiaFrame,
25512667 &tmpRotMat,
2668+ &tmpRotMat1stDeriv,
25522669 &tmpRotMat1stDerivs,
2670+ &tmpRotatedDiatomicOverlap,
2671+ &tmpMatrix,
25532672 &tmpDiatomicTwoElecTwoCore);
25542673 } // end of omp-parallelized region
25552674 // Exception throwing for omp-region
@@ -2565,8 +2684,13 @@ void Mndo::CalcForce(const vector<int>& elecStates){
25652684
25662685 void Mndo::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
25672686 double****** diatomicTwoElecTwoCore1stDerivs,
2687+ double*** tmpDiaOverlapAOsInDiaFrame,
2688+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
25682689 double*** tmpRotMat,
2690+ double*** tmpRotMat1stDeriv,
25692691 double**** tmpRotMat1stDerivs,
2692+ double*** tmpRotatedDiatomicOverlap,
2693+ double*** tmpMatrix,
25702694 double***** tmpDiatomicTwoElecTwoCore) const{
25712695 MallocerFreer::GetInstance()->Malloc<double>(diatomicOverlapAOs1stDerivs,
25722696 OrbitalType_end,
@@ -2578,13 +2702,28 @@ void Mndo::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
25782702 dxy,
25792703 dxy,
25802704 CartesianType_end);
2705+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiaOverlapAOsInDiaFrame,
2706+ OrbitalType_end,
2707+ OrbitalType_end);
2708+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiaOverlapAOs1stDerivInDiaFrame,
2709+ OrbitalType_end,
2710+ OrbitalType_end);
25812711 MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat,
25822712 OrbitalType_end,
25832713 OrbitalType_end);
2714+ MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat1stDeriv,
2715+ OrbitalType_end,
2716+ OrbitalType_end);
25842717 MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat1stDerivs,
25852718 OrbitalType_end,
25862719 OrbitalType_end,
25872720 CartesianType_end);
2721+ MallocerFreer::GetInstance()->Malloc<double>(tmpRotatedDiatomicOverlap,
2722+ OrbitalType_end,
2723+ OrbitalType_end);
2724+ MallocerFreer::GetInstance()->Malloc<double>(tmpMatrix,
2725+ OrbitalType_end,
2726+ OrbitalType_end);
25882727 MallocerFreer::GetInstance()->Malloc<double>(tmpDiatomicTwoElecTwoCore,
25892728 dxy,
25902729 dxy,
@@ -2592,10 +2731,15 @@ void Mndo::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
25922731 dxy);
25932732 }
25942733
2595-void Mndo::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
2734+void Mndo::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
25962735 double****** diatomicTwoElecTwoCore1stDerivs,
2736+ double*** tmpDiaOverlapAOsInDiaFrame,
2737+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
25972738 double*** tmpRotMat,
2739+ double*** tmpRotMat1stDeriv,
25982740 double**** tmpRotMat1stDerivs,
2741+ double*** tmpRotatedDiatomicOverlap,
2742+ double*** tmpMatrix,
25992743 double***** tmpDiatomicTwoElecTwoCore) const{
26002744 MallocerFreer::GetInstance()->Free<double>(diatomicOverlapAOs1stDerivs,
26012745 OrbitalType_end,
@@ -2607,13 +2751,28 @@ void Mndo::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
26072751 dxy,
26082752 dxy,
26092753 CartesianType_end);
2754+ MallocerFreer::GetInstance()->Free<double>(tmpDiaOverlapAOsInDiaFrame,
2755+ OrbitalType_end,
2756+ OrbitalType_end);
2757+ MallocerFreer::GetInstance()->Free<double>(tmpDiaOverlapAOs1stDerivInDiaFrame,
2758+ OrbitalType_end,
2759+ OrbitalType_end);
26102760 MallocerFreer::GetInstance()->Free<double>(tmpRotMat,
26112761 OrbitalType_end,
26122762 OrbitalType_end);
2763+ MallocerFreer::GetInstance()->Free<double>(tmpRotMat1stDeriv,
2764+ OrbitalType_end,
2765+ OrbitalType_end);
26132766 MallocerFreer::GetInstance()->Free<double>(tmpRotMat1stDerivs,
26142767 OrbitalType_end,
26152768 OrbitalType_end,
26162769 CartesianType_end);
2770+ MallocerFreer::GetInstance()->Free<double>(tmpRotatedDiatomicOverlap,
2771+ OrbitalType_end,
2772+ OrbitalType_end);
2773+ MallocerFreer::GetInstance()->Free<double>(tmpMatrix,
2774+ OrbitalType_end,
2775+ OrbitalType_end);
26172776 MallocerFreer::GetInstance()->Free<double>(tmpDiatomicTwoElecTwoCore,
26182777 dxy,
26192778 dxy,
--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -152,19 +152,29 @@ private:
152152 double******* diatomicTwoElecTwoCore1stDerivs,
153153 double******** diatomicTwoElecTwoCore2ndDerivs,
154154 double*** tmpRotMat,
155+ double*** tmpRotMat1stDeriv,
155156 double**** tmpRotMat1stDerivs,
156157 double***** tmpRotMat2ndDerivs,
157- double***** tmpDiatomicTwoElecTwo,
158- double****** tmpDiatomicTwoElecTwo1stDerivs) const;
158+ double***** tmpDiatomicTwoElecTwoCore,
159+ double****** tmpDiatomicTwoElecTwoCore1stDerivs,
160+ double*** tmpDiaOverlapAOsInDiaFrame,
161+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
162+ double*** tmpRotatedDiatomicOverlap,
163+ double*** tmpMatrix) const;
159164 void FreeTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverlapAOs1stDerivs,
160165 double****** diatomicOverlapAOs2ndDerivs,
161166 double******* diatomicTwoElecTwoCore1stDerivs,
162167 double******** diatomicTwoElecTwoCore2ndDerivs,
163168 double*** tmpRotMat,
169+ double*** tmpRotMat1stDeriv,
164170 double**** tmpRotMat1stDerivs,
165171 double***** tmpRotMat2ndDerivs,
166- double***** tmpDiatomicTwoElecTwo,
167- double****** tmpDiatomicTwoElecTwo1stDerivs) const;
172+ double***** tmpDiatomicTwoElecTwoCore,
173+ double****** tmpDiatomicTwoElecTwoCore1stDerivs,
174+ double*** tmpDiaOverlapAOsInDiaFrame,
175+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
176+ double*** tmpRotatedDiatomicOverlap,
177+ double*** tmpMatrix) const;
168178 double GetAuxiliaryHessianElement1(int mu,
169179 int nu,
170180 int indexAtomA,
@@ -364,14 +374,24 @@ private:
364374 MolDS_base::MultipoleType multipoleB,
365375 double rAB) const;
366376 void MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
367- double****** diatomiTwoElecTwoCore1stDerivs,
377+ double****** diatomicTwoElecTwoCore1stDerivs,
378+ double*** tmpDiaOverlapAOsInDiaFrame,
379+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
368380 double*** tmpRotMat,
381+ double*** tmpRotMat1stDeriv,
369382 double**** tmpRotMat1stDerivs,
383+ double*** tmpRotatedDiatomicOverlap,
384+ double*** tmpMatrix,
370385 double***** tmpDiatomicTwoElecTwoCore) const;
371386 void FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
372- double****** diatomiTwoElecTwoCore1stDerivs,
387+ double****** diatomicTwoElecTwoCore1stDerivs,
388+ double*** tmpDiaOverlapAOsInDiaFrame,
389+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
373390 double*** tmpRotMat,
391+ double*** tmpRotMat1stDeriv,
374392 double**** tmpRotMat1stDerivs,
393+ double*** tmpRotatedDiatomicOverlap,
394+ double*** tmpMatrix,
375395 double***** tmpDiatomicTwoElecTwoCore) const;
376396 void CalcForceSCFElecCoreAttractionPart(double* force,
377397 int indexAtomA,
--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -3593,17 +3593,25 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
35933593 stringstream ompErrors;
35943594 #pragma omp parallel
35953595 {
3596- double*** diatomicTwoElecTwoCore1stDerivs = NULL;
3597- double*** diatomicOverlapAOs1stDerivs = NULL;
3596+ double*** diatomicOverlapAOs1stDerivs = NULL;
3597+ double*** diatomicTwoElecTwoCore1stDerivs = NULL;
3598+ double** tmpDiaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
3599+ double** tmpDiaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
3600+ double** tmpRotMat = NULL; // rotating Matrix from the diatomic frame to space fixed frame.
3601+ double** tmpRotMat1stDeriv = NULL;
3602+ double*** tmpRotMat1stDerivs = NULL; // first derivatives of the rotMat.
3603+ double** tmpRotatedDiatomicOverlap = NULL;
3604+ double** tmpMatrix = NULL;
35983605 try{
3599- MallocerFreer::GetInstance()->Malloc<double>(&diatomicTwoElecTwoCore1stDerivs,
3600- OrbitalType_end,
3601- OrbitalType_end,
3602- CartesianType_end);
3603- MallocerFreer::GetInstance()->Malloc<double>(&diatomicOverlapAOs1stDerivs,
3604- OrbitalType_end,
3605- OrbitalType_end,
3606- CartesianType_end);
3606+ MallocTempMatricesCalcForce(&diatomicOverlapAOs1stDerivs,
3607+ &diatomicTwoElecTwoCore1stDerivs,
3608+ &tmpDiaOverlapAOsInDiaFrame,
3609+ &tmpDiaOverlapAOs1stDerivInDiaFrame,
3610+ &tmpRotMat,
3611+ &tmpRotMat1stDeriv,
3612+ &tmpRotMat1stDerivs,
3613+ &tmpRotatedDiatomicOverlap,
3614+ &tmpMatrix);
36073615 #pragma omp for schedule(auto)
36083616 for(int b=0; b<this->molecule->GetNumberAtoms(); b++){
36093617 if(a == b){continue;}
@@ -3613,7 +3621,16 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
36133621 double rAB = this->molecule->GetDistanceAtoms(atomA, atomB);
36143622
36153623 // calc. first derivative of overlapAOs.
3616- this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs, atomA, atomB);
3624+ this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs,
3625+ tmpDiaOverlapAOsInDiaFrame,
3626+ tmpDiaOverlapAOs1stDerivInDiaFrame,
3627+ tmpRotMat,
3628+ tmpRotMat1stDeriv,
3629+ tmpRotMat1stDerivs,
3630+ tmpRotatedDiatomicOverlap,
3631+ tmpMatrix,
3632+ atomA,
3633+ atomB);
36173634
36183635 // calc. first derivative of two elec two core interaction by Nishimoto-Mataga
36193636 this->CalcDiatomicTwoElecTwoCore1stDerivatives(diatomicTwoElecTwoCore1stDerivs, a, b);
@@ -3725,14 +3742,15 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
37253742 #pragma omp critical
37263743 ex.Serialize(ompErrors);
37273744 }
3728- MallocerFreer::GetInstance()->Free<double>(&diatomicTwoElecTwoCore1stDerivs,
3729- OrbitalType_end,
3730- OrbitalType_end,
3731- CartesianType_end);
3732- MallocerFreer::GetInstance()->Free<double>(&diatomicOverlapAOs1stDerivs,
3733- OrbitalType_end,
3734- OrbitalType_end,
3735- CartesianType_end);
3745+ FreeTempMatricesCalcForce(&diatomicOverlapAOs1stDerivs,
3746+ &diatomicTwoElecTwoCore1stDerivs,
3747+ &tmpDiaOverlapAOsInDiaFrame,
3748+ &tmpDiaOverlapAOs1stDerivInDiaFrame,
3749+ &tmpRotMat,
3750+ &tmpRotMat1stDeriv,
3751+ &tmpRotMat1stDerivs,
3752+ &tmpRotatedDiatomicOverlap,
3753+ &tmpMatrix);
37363754 } //end of omp-parallelized region
37373755 // Exception throwing for omp-region
37383756 if(!ompErrors.str().empty()){
@@ -3743,35 +3761,6 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
37433761 // communication to reduce thsi->matrixForce on all node (namely, all_reduce)
37443762 int numTransported = elecStates.size()*this->molecule->GetNumberAtoms()*CartesianType_end;
37453763 MolDS_mpi::MpiProcess::GetInstance()->AllReduce(&this->matrixForce[0][0][0], numTransported, std::plus<double>());
3746- /*
3747- double*** tmp=NULL;
3748- try{
3749- MallocerFreer::GetInstance()->Malloc<double>(&tmp,
3750- elecStates.size(),
3751- this->molecule->GetNumberAtoms(),
3752- CartesianType_end);
3753- MolDS_mpi::MpiProcess::GetInstance()->AllReduce(&this->matrixForce[0][0][0], numTransported, &tmp[0][0][0], std::plus<double>());
3754- for(int n=0; n<elecStates.size(); n++){
3755- for(int a=0; a<this->molecule->GetNumberAtoms(); a++){
3756- for(int i=0; i<CartesianType_end; i++){
3757- this->matrixForce[n][a][i] = tmp[n][a][i];
3758- }
3759- }
3760- }
3761- }
3762- catch(MolDSException ex){
3763- MallocerFreer::GetInstance()->Free<double>(&tmp,
3764- elecStates.size(),
3765- this->molecule->GetNumberAtoms(),
3766- CartesianType_end);
3767- throw ex;
3768- }
3769- MallocerFreer::GetInstance()->Free<double>(&tmp,
3770- elecStates.size(),
3771- this->molecule->GetNumberAtoms(),
3772- CartesianType_end);
3773- */
3774- // end of communication
37753764
37763765 /*
37773766 // Calculate force (on the ground state only).
@@ -3856,6 +3845,46 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
38563845 */
38573846 }
38583847
3848+void ZindoS::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
3849+ double**** diatomicTwoElecTwoCore1stDerivs,
3850+ double*** tmpDiaOverlapAOsInDiaFrame,
3851+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
3852+ double*** tmpRotMat,
3853+ double*** tmpRotMat1stDeriv,
3854+ double**** tmpRotMat1stDerivs,
3855+ double*** tmpRotatedDiatomicOverlap,
3856+ double*** tmpMatrix) const{
3857+ MallocerFreer::GetInstance()->Malloc<double>(diatomicOverlapAOs1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3858+ MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecTwoCore1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3859+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
3860+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
3861+ MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat, OrbitalType_end, OrbitalType_end);
3862+ MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat1stDeriv, OrbitalType_end, OrbitalType_end);
3863+ MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3864+ MallocerFreer::GetInstance()->Malloc<double>(tmpRotatedDiatomicOverlap, OrbitalType_end, OrbitalType_end);
3865+ MallocerFreer::GetInstance()->Malloc<double>(tmpMatrix, OrbitalType_end, OrbitalType_end);
3866+}
3867+
3868+void ZindoS::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
3869+ double**** diatomicTwoElecTwoCore1stDerivs,
3870+ double*** tmpDiaOverlapAOsInDiaFrame,
3871+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
3872+ double*** tmpRotMat,
3873+ double*** tmpRotMat1stDeriv,
3874+ double**** tmpRotMat1stDerivs,
3875+ double*** tmpRotatedDiatomicOverlap,
3876+ double*** tmpMatrix) const{
3877+ MallocerFreer::GetInstance()->Free<double>(diatomicOverlapAOs1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3878+ MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecTwoCore1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3879+ MallocerFreer::GetInstance()->Free<double>(tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
3880+ MallocerFreer::GetInstance()->Free<double>(tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
3881+ MallocerFreer::GetInstance()->Free<double>(tmpRotMat, OrbitalType_end, OrbitalType_end);
3882+ MallocerFreer::GetInstance()->Free<double>(tmpRotMat1stDeriv, OrbitalType_end, OrbitalType_end);
3883+ MallocerFreer::GetInstance()->Free<double>(tmpRotMat1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3884+ MallocerFreer::GetInstance()->Free<double>(tmpRotatedDiatomicOverlap, OrbitalType_end, OrbitalType_end);
3885+ MallocerFreer::GetInstance()->Free<double>(tmpMatrix, OrbitalType_end, OrbitalType_end);
3886+}
3887+
38593888 void ZindoS::CalcForceExcitedStaticPart(double* force,
38603889 int elecStateIndex,
38613890 int indexAtomA,
--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -271,6 +271,24 @@ private:
271271 void CalcDiatomicTwoElecTwoCore1stDerivatives(double*** matrix,
272272 int indexAtomA,
273273 int indexAtomB) const;
274+ void MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
275+ double**** diatomicTwoElecTwoCore1stDerivs,
276+ double*** tmpDiaOverlapAOsInDiaFrame,
277+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
278+ double*** tmpRotMat,
279+ double*** tmpRotMat1stDeriv,
280+ double**** tmpRotMat1stDerivs,
281+ double*** tmpRotatedDiatomicOverlap,
282+ double*** tmpMatrix) const;
283+ void FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
284+ double**** diatomicTwoElecTwoCore1stDerivs,
285+ double*** tmpDiaOverlapAOsInDiaFrame,
286+ double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
287+ double*** tmpRotMat,
288+ double*** tmpRotMat1stDeriv,
289+ double**** tmpRotMat1stDerivs,
290+ double*** tmpRotatedDiatomicOverlap,
291+ double*** tmpMatrix) const;
274292 void CalcForceExcitedStaticPart(double* force,
275293 int elecStateIndex,
276294 int indexAtomA,