提交: 1dcb9a22eaf6cdaab95f10fac5d962eca59e03ad - MolDS (git)

Commit MetaInfo

• 修订版: 1dcb9a22eaf6cdaab95f10fac5d962eca59e03ad (tree)
• 时间: 2013-10-18 10:54:22
• 作者: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Bug in handling of MolDSException in openMP/MP region is fixed. #32309

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

更改概述

• modified: src/cndo/Cndo2.cpp
• modified: src/mndo/Mndo.cpp
• modified: src/zindo/ZindoS.cpp

差异

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -1389,24 +1389,23 @@ void Cndo2::CalcFockMatrix(double** fockMatrix,
                            bool isGuess) const{
    int totalNumberAOs   = molecule.GetTotalNumberAOs();
    int totalNumberAtoms = molecule.GetNumberAtoms();
+   MallocerFreer::GetInstance()->Initialize<double>(fockMatrix, totalNumberAOs, totalNumberAOs);
 
    // MPI setting of each rank
    int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
    int mpiSize       = MolDS_mpi::MpiProcess::GetInstance()->GetSize();
    int mpiHeadRank   = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
+   stringstream errorStream;
    MolDS_mpi::AsyncCommunicator asyncCommunicator;
-   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, 
-                                                  &asyncCommunicator) );
+   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
 
-   MallocerFreer::GetInstance()->Initialize<double>(fockMatrix, totalNumberAOs, totalNumberAOs);
-   for(int A=totalNumberAtoms-1; 0<=A; A--){
+   for(int A=0; A<totalNumberAtoms; A++){
       const Atom& atomA = *molecule.GetAtom(A);
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
       for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
          int calcRank = mu%mpiSize;
          if(mpiRank == calcRank){
-            stringstream ompErrors;
 #pragma omp parallel for schedule(auto)
             for(int B=A; B<totalNumberAtoms; B++){
                try{
@@ -1444,36 +1443,34 @@ void Cndo2::CalcFockMatrix(double** fockMatrix,
                      else{
                         // lower left part (not calculated)
                      }
-                  }  // end of loop nu
-               }  // end of try
+                  }
+               }
                catch(MolDSException ex){
-#pragma omp critical
-                  ex.Serialize(ompErrors);
+#pragma omp critical 
+                  ex.Serialize(errorStream);
                }
-            }  // end of loop B parallelized with openMP
-            // Exception throwing for omp-region
-            if(!ompErrors.str().empty()){
-               throw MolDSException::Deserialize(ompErrors);
             }
-         }  // end of if(mpiRank == calcRank) 
-
-         // set data to gather in mpiHeadRank with asynchronous MPI 
-         int tag                      = mu;
-         int source                   = calcRank;
-         int dest                     = mpiHeadRank;
-         double* buff                 = &fockMatrix[mu][mu];
-         MolDS_mpi::molds_mpi_int num = totalNumberAOs-mu;
-         if(mpiRank == mpiHeadRank && mpiRank != calcRank){
-            asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
          }
-         if(mpiRank != mpiHeadRank && mpiRank == calcRank){
-            asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+         if(errorStream.str().empty()){
+            int tag                      = mu;
+            int source                   = calcRank;
+            int dest                     = mpiHeadRank;
+            double* buff                 = &fockMatrix[mu][mu];
+            MolDS_mpi::molds_mpi_int num = totalNumberAOs-mu;
+            if(mpiRank == mpiHeadRank && mpiRank != calcRank){
+               asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
+            }
+            if(mpiRank != mpiHeadRank && mpiRank == calcRank){
+               asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+            }
          }
-      }  // end of loop mu parallelized with MPI
-   }  // end of loop A 
-   // Delete the communication thread.
+      }
+   }
    asyncCommunicator.Finalize();
    communicationThread.join();
+   if(!errorStream.str().empty()){
+      throw MolDSException::Deserialize(errorStream);
+   }
    double* buff                 = &fockMatrix[0][0];
    MolDS_mpi::molds_mpi_int num = totalNumberAOs*totalNumberAOs;
    MolDS_mpi::MpiProcess::GetInstance()->Broadcast(buff, num, mpiHeadRank);
@@ -1609,11 +1606,10 @@ void Cndo2::CalcGammaAB(double** gammaAB, const Molecule& molecule) const{
    int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
    int mpiSize       = MolDS_mpi::MpiProcess::GetInstance()->GetSize();
    int mpiHeadRank   = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
+   stringstream errorStream;
    MolDS_mpi::AsyncCommunicator asyncCommunicator;
-   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, 
-                                                  &asyncCommunicator) );
+   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
 
-   // This loop (A) is parallelized by MPI
    for(int A=0; A<totalAtomNumber; A++){
       int calcRank = A%mpiSize;
       if(mpiRank == calcRank){
@@ -1621,7 +1617,6 @@ void Cndo2::CalcGammaAB(double** gammaAB, const Molecule& molecule) const{
          int na = atomA.GetValenceShellType() + 1;
          double orbitalExponentA = atomA.GetOrbitalExponent(
                                          atomA.GetValenceShellType(), s, this->theory);
-         stringstream ompErrors;
 #pragma omp parallel for schedule(auto)
          for(int B=A; B<totalAtomNumber; B++){
             try{
@@ -1675,30 +1670,29 @@ void Cndo2::CalcGammaAB(double** gammaAB, const Molecule& molecule) const{
             }
             catch(MolDSException ex){
             #pragma omp critical
-               ex.Serialize(ompErrors);
+               ex.Serialize(errorStream);
             }
-         }  // end of loop B parallelized by openMP
-         // Exception throwing for omp-region
-         if(!ompErrors.str().empty()){
-            throw MolDSException::Deserialize(ompErrors);
          }
-      }  //  end of if(mpiRank==calcRank)
-
-      // set data to gater in mpiHeadRank with asynchronous MPI 
-      int tag                      = A;
-      int source                   = calcRank;
-      int dest                     = mpiHeadRank;
-      double* buff                 = &gammaAB[A][A];
-      MolDS_mpi::molds_mpi_int num = totalAtomNumber-A;
-      if(mpiRank == mpiHeadRank && mpiRank != calcRank){
-         asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
       }
-      if(mpiRank != mpiHeadRank && mpiRank == calcRank){
-         asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+      if(errorStream.str().empty()){
+         int tag                      = A;
+         int source                   = calcRank;
+         int dest                     = mpiHeadRank;
+         double* buff                 = &gammaAB[A][A];
+         MolDS_mpi::molds_mpi_int num = totalAtomNumber-A;
+         if(mpiRank == mpiHeadRank && mpiRank != calcRank){
+            asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
+         }
+         if(mpiRank != mpiHeadRank && mpiRank == calcRank){
+            asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+         }
       }
-   }  // end of loop A prallelized by MPI
+   }
    asyncCommunicator.Finalize();
    communicationThread.join();
+   if(!errorStream.str().empty()){
+      throw MolDSException::Deserialize(errorStream);
+   }
    double* buff                 = &gammaAB[0][0];
    MolDS_mpi::molds_mpi_int num = totalAtomNumber*totalAtomNumber;
    MolDS_mpi::MpiProcess::GetInstance()->Broadcast(buff, num, mpiHeadRank);
@@ -1807,11 +1801,10 @@ void Cndo2::CalcCartesianMatrixByGTOExpansion(double*** cartesianMatrix,
    int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
    int mpiSize       = MolDS_mpi::MpiProcess::GetInstance()->GetSize();
    int mpiHeadRank   = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
+   stringstream errorStream;
    MolDS_mpi::AsyncCommunicator asyncCommunicator;
-   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, 
-                                                  &asyncCommunicator) );
+   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
 
-   // This loop (A and mu) is parallelized by MPI
    for(int A=0; A<totalAtomNumber; A++){
       const Atom& atomA  = *molecule.GetAtom(A);
       int firstAOIndexA  = atomA.GetFirstAOIndex();
@@ -1820,7 +1813,6 @@ void Cndo2::CalcCartesianMatrixByGTOExpansion(double*** cartesianMatrix,
       if(mpiRank == calcRank){
          for(int a=0; a<numValenceAOsA; a++){
             int mu = firstAOIndexA + a;      
-            stringstream ompErrors;
 #pragma omp parallel for schedule(auto)
             for(int B=0; B<totalAtomNumber; B++){
                try{
@@ -1837,40 +1829,38 @@ void Cndo2::CalcCartesianMatrixByGTOExpansion(double*** cartesianMatrix,
                }
                catch(MolDSException ex){
 #pragma omp critical
-                  ex.Serialize(ompErrors);
+                  ex.Serialize(errorStream);
                }  
-            }// end of loop for int B with openMP
-            // Exception throwing for omp-region
-            if(!ompErrors.str().empty()){
-               throw MolDSException::Deserialize(ompErrors);
-            }
-         } 
-      } // end lof if(mpiRank == calcRank)
-
-      // set data to gater in mpiHeadRank with asynchronous MPI 
-      int tagX                     = A* CartesianType_end + XAxis;
-      int tagY                     = A* CartesianType_end + YAxis;
-      int tagZ                     = A* CartesianType_end + ZAxis;
-      int source                   = calcRank;
-      int dest                     = mpiHeadRank;
-      double* buffX                = &cartesianMatrix[XAxis][firstAOIndexA][0];
-      double* buffY                = &cartesianMatrix[YAxis][firstAOIndexA][0];
-      double* buffZ                = &cartesianMatrix[ZAxis][firstAOIndexA][0];
-      MolDS_mpi::molds_mpi_int num = numValenceAOsA*totalAONumber;
-      if(mpiRank == mpiHeadRank && mpiRank != calcRank){
-         asyncCommunicator.SetRecvedMessage(buffX, num, source, tagX);
-         asyncCommunicator.SetRecvedMessage(buffY, num, source, tagY);
-         asyncCommunicator.SetRecvedMessage(buffZ, num, source, tagZ);
-      }
-      if(mpiRank != mpiHeadRank && mpiRank == calcRank){
-         asyncCommunicator.SetSentMessage(buffX, num, dest, tagX);
-         asyncCommunicator.SetSentMessage(buffY, num, dest, tagY);
-         asyncCommunicator.SetSentMessage(buffZ, num, dest, tagZ);
-      }
-   } // end of loop for int A with MPI
-   // Delete the communication thread.
+            } 
+         }
+      }
+      if(errorStream.str().empty()){
+         int tagX                     = A* CartesianType_end + XAxis;
+         int tagY                     = A* CartesianType_end + YAxis;
+         int tagZ                     = A* CartesianType_end + ZAxis;
+         int source                   = calcRank;
+         int dest                     = mpiHeadRank;
+         double* buffX                = &cartesianMatrix[XAxis][firstAOIndexA][0];
+         double* buffY                = &cartesianMatrix[YAxis][firstAOIndexA][0];
+         double* buffZ                = &cartesianMatrix[ZAxis][firstAOIndexA][0];
+         MolDS_mpi::molds_mpi_int num = numValenceAOsA*totalAONumber;
+         if(mpiRank == mpiHeadRank && mpiRank != calcRank){
+            asyncCommunicator.SetRecvedMessage(buffX, num, source, tagX);
+            asyncCommunicator.SetRecvedMessage(buffY, num, source, tagY);
+            asyncCommunicator.SetRecvedMessage(buffZ, num, source, tagZ);
+         }
+         if(mpiRank != mpiHeadRank && mpiRank == calcRank){
+            asyncCommunicator.SetSentMessage(buffX, num, dest, tagX);
+            asyncCommunicator.SetSentMessage(buffY, num, dest, tagY);
+            asyncCommunicator.SetSentMessage(buffZ, num, dest, tagZ);
+         }
+      }
+   }
    asyncCommunicator.Finalize();
    communicationThread.join();
+   if(!errorStream.str().empty()){
+      throw MolDSException::Deserialize(errorStream);
+   }
    double* buff                 = &cartesianMatrix[0][0][0];
    MolDS_mpi::molds_mpi_int num = CartesianType_end*totalAONumber*totalAONumber;
    MolDS_mpi::MpiProcess::GetInstance()->Broadcast(buff, num, mpiHeadRank);
@@ -3904,27 +3894,22 @@ void Cndo2::CalcOverlapESsWithAnotherElectronicStructure(double** overlapESs,
 void Cndo2::CalcOverlapAOs(double** overlapAOs, const Molecule& molecule) const{
    int totalAONumber = molecule.GetTotalNumberAOs();
    int totalAtomNumber = molecule.GetNumberAtoms();
+   MallocerFreer::GetInstance()->Initialize<double>(overlapAOs, totalAONumber, totalAONumber);
 
    // MPI setting of each rank
    int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
    int mpiSize       = MolDS_mpi::MpiProcess::GetInstance()->GetSize();
    int mpiHeadRank   = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
+   stringstream errorStream;
    MolDS_mpi::AsyncCommunicator asyncCommunicator;
-   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, 
-                                                  &asyncCommunicator) );
-
-   MallocerFreer::GetInstance()->Initialize<double>(overlapAOs,
-                                                    totalAONumber,
-                                                    totalAONumber);
+   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
 
-   // This loop A is parallelized with MPI
-   for(int A=totalAtomNumber-1; 0<=A; A--){
+   for(int A=0; A<totalAtomNumber; A++){
       const Atom& atomA = *molecule.GetAtom(A);
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int numValenceAOs = atomA.GetValenceSize();
       int calcRank = A%mpiSize;
       if(mpiRank == calcRank){
-         stringstream ompErrors;
 #pragma omp parallel 
          {
             double** diatomicOverlapAOs       = NULL;
@@ -3959,12 +3944,11 @@ void Cndo2::CalcOverlapAOs(double** overlapAOs, const Molecule& molecule) const{
                   this->CalcRotatingMatrix(rotatingMatrix, atomA, atomB);
                   this->RotateDiatmicOverlapAOsToSpaceFrame(diatomicOverlapAOs, rotatingMatrix, tmpDiatomicOverlapAOs, tmpOldDiatomicOverlapAOs, tmpMatrixBC, tmpVectorBC);
                   this->SetOverlapAOsElement(overlapAOs, diatomicOverlapAOs, atomA, atomB, symmetrize);
-               } // end of loop B parallelized with openMP
-
-            }  // end of try
+               } 
+            }
             catch(MolDSException ex){
 #pragma omp critical
-               ex.Serialize(ompErrors);
+               ex.Serialize(errorStream);
             }
             this->FreeDiatomicOverlapAOsAndRotatingMatrix(&diatomicOverlapAOs, &rotatingMatrix);
             MallocerFreer::GetInstance()->Free<double>(&tmpDiatomicOverlapAOs,
@@ -3977,28 +3961,27 @@ void Cndo2::CalcOverlapAOs(double** overlapAOs, const Molecule& molecule) const{
                                                        OrbitalType_end);
             MallocerFreer::GetInstance()->Free<double>(&tmpVectorBC,              
                                                        OrbitalType_end*OrbitalType_end);
-         }  // end of omp-parallelized region
-         // Exception throwing for omp-region
-         if(!ompErrors.str().empty()){
-            throw MolDSException::Deserialize(ompErrors);
          }
-      } // end of if(mpiRank == calcRnak)
-
-      // set data to gather in mpiHeadRank with asynchronous MPI 
-      int tag                      = A;
-      int source                   = calcRank;
-      int dest                     = mpiHeadRank;
-      double* buff                 = overlapAOs[firstAOIndexA];
-      MolDS_mpi::molds_mpi_int num = totalAONumber*numValenceAOs;
-      if(mpiRank == mpiHeadRank && mpiRank != calcRank){
-         asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
-      }
-      if(mpiRank != mpiHeadRank && mpiRank == calcRank){
-         asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+      } 
+      if(errorStream.str().empty()){
+         int tag                      = A;
+         int source                   = calcRank;
+         int dest                     = mpiHeadRank;
+         double* buff                 = overlapAOs[firstAOIndexA];
+         MolDS_mpi::molds_mpi_int num = totalAONumber*numValenceAOs;
+         if(mpiRank == mpiHeadRank && mpiRank != calcRank){
+            asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
+         }
+         if(mpiRank != mpiHeadRank && mpiRank == calcRank){
+            asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+         }
       }
-   }  // end of loop A parallelized with MPI
+   }
    asyncCommunicator.Finalize();
    communicationThread.join();
+   if(!errorStream.str().empty()){
+      throw MolDSException::Deserialize(errorStream);
+   }
    double* buff                 = &overlapAOs[0][0];
    MolDS_mpi::molds_mpi_int num = totalAONumber*totalAONumber;
    MolDS_mpi::MpiProcess::GetInstance()->Broadcast(buff, num, mpiHeadRank);

--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -3487,30 +3487,28 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
 
 void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore, 
                               const Molecule& molecule) const{
-   int totalNumberAtoms = molecule.GetNumberAtoms();
-
-   // MPI setting of each rank
-   int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
-   int mpiSize       = MolDS_mpi::MpiProcess::GetInstance()->GetSize();
-   int mpiHeadRank   = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
-   MolDS_mpi::AsyncCommunicator asyncCommunicator;
-   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, 
-                                                  &asyncCommunicator) );
 #ifdef MOLDS_DBG
    if(twoElecTwoCore == NULL){
       throw MolDSException(this->errorMessageCalcTwoElecTwoCoreNullMatrix);
    }
 #endif
+   int totalNumberAtoms = molecule.GetNumberAtoms();
    MallocerFreer::GetInstance()->Initialize<double>(twoElecTwoCore, 
                                                     totalNumberAtoms,
                                                     totalNumberAtoms,
                                                     dxy, dxy, dxy, dxy);
 
-   // this loop-a is MPI-parallelized
-   for(int a=totalNumberAtoms-1; 0<=a; a--){
+   // MPI setting of each rank
+   int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
+   int mpiSize       = MolDS_mpi::MpiProcess::GetInstance()->GetSize();
+   int mpiHeadRank   = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
+   stringstream errorStream;
+   MolDS_mpi::AsyncCommunicator asyncCommunicator;
+   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
+
+   for(int a=0; a<totalNumberAtoms; a++){
       int calcRank = a%mpiSize;
       if(mpiRank == calcRank){
-         stringstream ompErrors;
 #pragma omp parallel 
          {
             double**** diatomicTwoElecTwoCore    = NULL;
@@ -3533,14 +3531,12 @@ void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore,
                                                    tmpMatrixBC, 
                                                    tmpVectorBC, 
                                                    a, b);
-
                   int i=0;
                   for(int mu=0; mu<dxy; mu++){
                      for(int nu=mu; nu<dxy; nu++){
                         int j=0;
                         for(int lambda=0; lambda<dxy; lambda++){
                            for(int sigma=lambda; sigma<dxy; sigma++){
-                              //double value = diatomicTwoElecTwoCore[mu][nu][lambda][sigma];
                               this->twoElecTwoCoreMpiBuff[a][b][i][j] 
                                  = diatomicTwoElecTwoCore[mu][nu][lambda][sigma];
                               j++;
@@ -3549,36 +3545,34 @@ void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore,
                         i++;
                      }
                   }
-
-               }  // end of loop b parallelized with MPI
-
-            }  // end of try
+               }
+            }
             catch(MolDSException ex){
 #pragma omp critical
-               ex.Serialize(ompErrors);
+               ex.Serialize(errorStream);
             }
             MallocerFreer::GetInstance()->Free<double>(&diatomicTwoElecTwoCore,    dxy, dxy, dxy, dxy);
             MallocerFreer::GetInstance()->Free<double>(&tmpDiatomicTwoElecTwoCore, dxy*dxy*dxy*dxy);
             MallocerFreer::GetInstance()->Free<double>(&tmpRotMat,                 OrbitalType_end, OrbitalType_end);
             MallocerFreer::GetInstance()->Free<double>(&tmpMatrixBC,               dxy*dxy, dxy*dxy);
             MallocerFreer::GetInstance()->Free<double>(&tmpVectorBC,                     dxy*dxy*dxy*dxy);
-         }  // end of omp-parallelized region
-         // Exception throwing for omp-region
-         if(!ompErrors.str().empty()){
-            throw MolDSException::Deserialize(ompErrors);
          }
-      } // end of if(mpiRnak == calcRank)
-      // set data to gather in mpiHeadRank with asynchronous MPI 
-      if(a<totalNumberAtoms-1){
-         int b = a+1;
-         OrbitalType twoElecLimit = dxy;
-         int numBuff = (twoElecLimit+1)*twoElecLimit/2;
-         int num = (totalNumberAtoms-b)*numBuff*numBuff;
-         asyncCommunicator.SetBroadcastedMessage(&this->twoElecTwoCoreMpiBuff[a][b][0][0], num, calcRank);
       }
-   } // end of loop a parallelized with MPI
+      if(errorStream.str().empty()){
+         if(a<totalNumberAtoms-1){
+            int b = a+1;
+            OrbitalType twoElecLimit = dxy;
+            int numBuff = (twoElecLimit+1)*twoElecLimit/2;
+            int num = (totalNumberAtoms-b)*numBuff*numBuff;
+            asyncCommunicator.SetBroadcastedMessage(&this->twoElecTwoCoreMpiBuff[a][b][0][0], num, calcRank);
+         }
+      }
+   }
    asyncCommunicator.Finalize();
    communicationThread.join();
+   if(!errorStream.str().empty()){
+      throw MolDSException::Deserialize(errorStream);
+   }
 
 #pragma omp parallel for schedule(auto)
    for(int a=0; a<totalNumberAtoms; a++){

--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -2353,18 +2353,16 @@ void ZindoS::CalcCISMatrix(double** matrixCIS) const{
    int mpiRank     = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
    int mpiSize     = MolDS_mpi::MpiProcess::GetInstance()->GetSize();
    int mpiHeadRank = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
+   stringstream errorStream;
    MolDS_mpi::AsyncCommunicator asyncCommunicator;
-   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>,
-                                                  &asyncCommunicator) );
+   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
 
-   // this loop-a is MPI-parallelized
-   for(int k=this->matrixCISdimension-1; 0<=k; k--){
+   for(int k=0; k<this->matrixCISdimension; k++){
       int calcRank = k%mpiSize;
       if(calcRank == mpiRank){
          // single excitation from I-th (occupied)MO to A-th (virtual)MO
          int moI = this->GetActiveOccIndex(*this->molecule, k);
          int moA = this->GetActiveVirIndex(*this->molecule, k);
-         stringstream ompErrors;
 #pragma omp parallel for schedule(auto)
          for(int l=k; l<this->matrixCISdimension; l++){
             try{
@@ -2401,30 +2399,29 @@ void ZindoS::CalcCISMatrix(double** matrixCIS) const{
             }
             catch(MolDSException ex){
 #pragma omp critical
-               ex.Serialize(ompErrors);
+               ex.Serialize(errorStream);
             }
-         } // end of l-loop
-         // Exception throwing for omp-region
-         if(!ompErrors.str().empty()){
-            throw MolDSException::Deserialize(ompErrors);
          }
-      } // end of if(calcRank == mpiRank)
-      // Send data to head rank
-      int tag      = k;
-      int source   = calcRank;
-      int dest     = mpiHeadRank;
-      int num      = this->matrixCISdimension - k;
-      double* buff = &this->matrixCIS[k][k];
-      if(mpiRank == mpiHeadRank && mpiRank != calcRank){
-         asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
       }
-      if(mpiRank != mpiHeadRank && mpiRank == calcRank){
-         asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+      if(errorStream.str().empty()){
+         int tag      = k;
+         int source   = calcRank;
+         int dest     = mpiHeadRank;
+         int num      = this->matrixCISdimension - k;
+         double* buff = &this->matrixCIS[k][k];
+         if(mpiRank == mpiHeadRank && mpiRank != calcRank){
+            asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
+         }
+         if(mpiRank != mpiHeadRank && mpiRank == calcRank){
+            asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+         }
       }
-   } // end of k-loop which is MPI-parallelized
+   }
    asyncCommunicator.Finalize();
    communicationThread.join();
-   // Broadcast data to all rank
+   if(!errorStream.str().empty()){
+      throw MolDSException::Deserialize(errorStream);
+   }
    for(int k=0; k<this->matrixCISdimension; k++){
       int     num  = this->matrixCISdimension - k;
       double* buff = &this->matrixCIS[k][k];
@@ -3319,20 +3316,19 @@ void ZindoS::CalcAuxiliaryVector(double* y,
 // Note taht K_{NR} is not calculated.
 void ZindoS::CalcGammaNRMinusKNRMatrix(double** gammaNRMinusKNR, const vector<MoIndexPair>& nonRedundantQIndeces) const{
    int nonRedundantQIndecesSize = nonRedundantQIndeces.size();
-   //MPI setting of each rank
+   // MPI setting of each rank
    int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
    int mpiSize       = MolDS_mpi::MpiProcess::GetInstance()->GetSize();
    int mpiHeadRank   = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
+   stringstream errorStream;
    MolDS_mpi::AsyncCommunicator asyncCommunicator;
-   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>,
-                                                  &asyncCommunicator) );
-   // this loop-i is MPI-parallelized
-   for(int i=nonRedundantQIndecesSize-1; 0<=i; i--){
+   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
+
+   for(int i=0; i<nonRedundantQIndecesSize; i++){
       int calcRank = i%mpiSize;
       if(mpiRank == calcRank){
          int moI = nonRedundantQIndeces[i].moI;
          int moJ = nonRedundantQIndeces[i].moJ;
-         stringstream ompErrors;
 #pragma omp parallel for schedule(auto)
          for(int j=i; j<nonRedundantQIndecesSize; j++){
             try{
@@ -3343,30 +3339,29 @@ void ZindoS::CalcGammaNRMinusKNRMatrix(double** gammaNRMinusKNR, const vector<Mo
             } // end of try
             catch(MolDSException ex){
 #pragma omp critical
-               ex.Serialize(ompErrors);
+               ex.Serialize(errorStream);
             }
-         } //end of loop j parallelized with openMP
-         // Exception throwing for omp-region
-         if(!ompErrors.str().empty()){
-            throw MolDSException::Deserialize(ompErrors);
-         } 
-      } /// end of if(mpiRnak == calcRank)
-      // Send data to head rank
-      int tag    = i;
-      int source = calcRank;
-      int dest   = mpiHeadRank;
-      int num    = nonRedundantQIndecesSize - i;
-      double* buff = &gammaNRMinusKNR[i][i];
-      if(mpiRank == mpiHeadRank && mpiRank != calcRank){
-         asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
+         }
       }
-      if(mpiRank != mpiHeadRank && mpiRank == calcRank){
-         asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+      if(errorStream.str().empty()){
+         int tag      = i;
+         int source   = calcRank;
+         int dest     = mpiHeadRank;
+         int num      = nonRedundantQIndecesSize - i;
+         double* buff = &gammaNRMinusKNR[i][i];
+         if(mpiRank == mpiHeadRank && mpiRank != calcRank){
+            asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
+         }
+         if(mpiRank != mpiHeadRank && mpiRank == calcRank){
+            asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+         }
       }
-   } // end of loop-i parallelized with MPI
+   }
    asyncCommunicator.Finalize();
    communicationThread.join();
-   // broadcast data to all rank
+   if(!errorStream.str().empty()){
+      throw MolDSException::Deserialize(errorStream);
+   } 
    for(int i=0; i<nonRedundantQIndecesSize; i++){
       int     num  = nonRedundantQIndecesSize - i;
       double* buff = &gammaNRMinusKNR[i][i];
@@ -3382,20 +3377,19 @@ void ZindoS::CalcKRDagerGammaRInvMatrix(double** kRDagerGammaRInv,
                                       const vector<MoIndexPair>& redundantQIndeces) const{
    int nonRedundantQIndecesSize = nonRedundantQIndeces.size();
    int redundantQIndecesSize    = redundantQIndeces.size();
-   //MPI setting of each rank
+   // MPI setting of each rank
    int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
    int mpiSize       = MolDS_mpi::MpiProcess::GetInstance()->GetSize();
    int mpiHeadRank   = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
+   stringstream errorStream;
    MolDS_mpi::AsyncCommunicator asyncCommunicator;
-   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>,
-                                                  &asyncCommunicator) );
-   // this loop-i is MPI-parallelized
+   boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
+
    for(int i=0; i<nonRedundantQIndecesSize; i++){
       int calcRank = i%mpiSize;
       if(mpiRank == calcRank){
          int moI = nonRedundantQIndeces[i].moI;
          int moJ = nonRedundantQIndeces[i].moJ;
-         stringstream ompErrors;
 #pragma omp parallel for schedule(auto)
          for(int j=0; j<redundantQIndecesSize; j++){
             try{
@@ -3406,30 +3400,29 @@ void ZindoS::CalcKRDagerGammaRInvMatrix(double** kRDagerGammaRInv,
             } // end of try
             catch(MolDSException ex){
 #pragma omp critical
-               ex.Serialize(ompErrors);
+               ex.Serialize(errorStream);
             }
-         } // end of loop-j parallelized with openMP
-         // Exception throwing for omp-region
-         if(!ompErrors.str().empty()){
-            throw MolDSException::Deserialize(ompErrors);
          }
-      } // // end of if(mpiRnak == calcRank)
-      // Send data to head rank
-      int tag      = i;
-      int source   = calcRank;
-      int dest     = mpiHeadRank;
-      int num      = redundantQIndecesSize;
-      double* buff = &kRDagerGammaRInv[i][0];
-      if(mpiRank == mpiHeadRank && mpiRank != calcRank){
-         asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
       }
-      if(mpiRank != mpiHeadRank && mpiRank == calcRank){
-         asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+      if(errorStream.str().empty()){
+         int tag      = i;
+         int source   = calcRank;
+         int dest     = mpiHeadRank;
+         int num      = redundantQIndecesSize;
+         double* buff = &kRDagerGammaRInv[i][0];
+         if(mpiRank == mpiHeadRank && mpiRank != calcRank){
+            asyncCommunicator.SetRecvedMessage(buff, num, source, tag);
+         }
+         if(mpiRank != mpiHeadRank && mpiRank == calcRank){
+            asyncCommunicator.SetSentMessage(buff, num, dest, tag);
+         }
       }
-   } // end of loop-i parallelized with MPI
+   }
    asyncCommunicator.Finalize();
    communicationThread.join();
-   // broadcast data to all rank
+   if(!errorStream.str().empty()){
+      throw MolDSException::Deserialize(errorStream);
+   }
    for(int i=0; i<nonRedundantQIndecesSize; i++){
       int     num  = redundantQIndecesSize;
       double* buff = &kRDagerGammaRInv[i][0];