Commit: 204fdd58bab73315cca63fed089ff5430452f151 - MolDS (git)

Commit MetaInfo

• Revision: 204fdd58bab73315cca63fed089ff5430452f151 (tree)
• Time: 2013-02-05 21:27:01
• Author: Katsuhiko Nishimra <ktns.87@gmai...>
• Commiter: Katsuhiko Nishimra

Log Message

Merge trunk into branch automake. #28588

git-svn-id: https://svn.sourceforge.jp/svnroot/molds/branches/automake@1289 1136aad2-a195-0410-b898-f5ea1d11b9d8

Change Summary

• modified: doc/README.txt
• modified: src/Main.cpp
• modified: src/am1/Am1.cpp
• modified: src/am1/Am1D.cpp
• modified: src/base/Molecule.cpp
• modified: src/base/Molecule.h
• modified: src/base/Parameters.cpp
• modified: src/base/atoms/Atom.cpp
• modified: src/base/atoms/Atom.h
• modified: src/cndo/Cndo2.cpp
• modified: src/cndo/Cndo2.h
• modified: src/mndo/Mndo.cpp
• modified: src/pm3/Pm3.cpp
• modified: src/pm3/Pm3D.cpp
• modified: src/pm3/Pm3Pddg.cpp
• modified: src/wrappers/Blas.cpp
• modified: src/wrappers/Blas.h
• modified: src/zindo/ZindoS.cpp
• modified: src/zindo/ZindoS.h
• modified: test/c2h6_pm3_directCIS_singlet.dat
• modified: test/c2h6_pm3_directCIS_singlet.in

Incremental Difference

--- a/doc/README.txt
+++ b/doc/README.txt
@@ -1,6 +1,7 @@
 //************************************************************************//
 // Copyright (C) 2011-2012 Mikiya Fujii                                   // 
 // Copyright (C) 2012-2012 Katsuhiko Nishimra                             // 
+// Copyright (C) 2013-2013 Michihiro Okuyama                              // 
 //                                                                        // 
 // This file is part of MolDS.                                            // 
 //                                                                        // 
@@ -66,6 +67,9 @@ COMPILE(using GNUmake):
 
    For both case, the compile succeeded if you could fine "MolDS.out" in the "src" directory. 
    Type "$ make clean" when you wanna clean the compilation.
+   If you want to compile MolDS in debug-mode, -g and -DMOLDS_DBG 
+   should be added to CFLAGBASE in the Makefile,
+   i.e. -O0 -openmp -openmp-report2 -DMKL_INT=intptr_t -g -DMOLDS_DBG
 
 ==============================================================================
 CARRY OUT MolDS:
@@ -276,7 +280,7 @@ HOW TO WRITE INPUT:
       -options
        "davidson", "active_occ", "active_vir", "max_iter", "max_dim", "norm_tol", 
        "nstates", "exciton_energies", "all_transition_dipole_moments", 
-       "mulliken", and "num_print_coefficients" are prepared as options.
+       "mulliken", "unpaired_electron_population", and "num_print_coefficients" are prepared as options.
 
        "davidson" should be set as "yes" or "no". 
        The default value of the "davidson" is "yes".
@@ -318,9 +322,16 @@ HOW TO WRITE INPUT:
 
        "mulliken" is a option of mulliken popultaion analysis of the excited state.
        When "mulliken x" is included in CIS-directive, the mulliken popultaion of xth excited state is calculated.
-       Mulitiple indication of these mulliken options is possible. 
+       Multiple indication of these mulliken options is possible. 
        Note that "mulliken 0" is ignored because 0th excited state is the ground state.
-       Default settign of this "mulliken" option is nothing.
+       Default setting of this "mulliken" option is nothing.
+
+       "unpaired_electron_population" is a option of unpaired electron population(UEP) analysis of the excited state.
+       When "unpaired electron population yes" and "mulliken x" are included in CIS-directive, 
+       the UEP of xth excited state is calculated.
+       By multiple indication of these mulliken option, the UEP on multiple excited states are possible.
+       Note that the UEP on ground state is ignored. 
+       Default setting is "unpaired_electron_population" option is nothing.
 
        "num_print_coefficients" is a number of the coefficients of CIS-eigenvector shown in output.
        The default value of the "num_print_coefficients" is 1.
@@ -336,6 +347,7 @@ HOW TO WRITE INPUT:
             norm_tol 0.000001
             mulliken 1 
             mulliken 2
+	         unpaired_electron_population yes
          CIS_END
 
    <Hole Plot>

--- a/src/Main.cpp
+++ b/src/Main.cpp
@@ -24,6 +24,7 @@
 #include<boost/shared_ptr.hpp>
 #include<boost/format.hpp>
 #include"base/PrintController.h"
+#include"base/MolDSException.h"
 #include"base/Enums.h"
 #include"base/EularAngle.h"
 #include"base/atoms/Atom.h"

--- a/src/am1/Am1.cpp
+++ b/src/am1/Am1.cpp
@@ -25,6 +25,7 @@
 #include<vector>
 #include<boost/format.hpp>
 #include"../base/PrintController.h"
+#include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
 #include"../base/Enums.h"
 #include"../base/EularAngle.h"

--- a/src/am1/Am1D.cpp
+++ b/src/am1/Am1D.cpp
@@ -25,6 +25,7 @@
 #include<vector>
 #include<boost/format.hpp>
 #include"../base/PrintController.h"
+#include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
 #include"../base/Enums.h"
 #include"../base/EularAngle.h"

--- a/src/base/Molecule.cpp
+++ b/src/base/Molecule.cpp
@@ -95,6 +95,7 @@ void Molecule::Initialize(){
    this->wasCalculatedXyzCOC = false;
    this->xyzCOM = NULL;
    this->xyzCOC = NULL;
+   this->atomVect = NULL;
    try{
       this->atomVect = new vector<Atom*>;
       MallocerFreer::GetInstance()->Malloc<double>(&this->xyzCOM, CartesianType_end);
@@ -124,7 +125,9 @@ void Molecule::Finalize(vector<Atom*>** atomVect, double** xyzCOM, double**xyzCO
 }
 
 void Molecule::SetMessages(){
-   this->errorMessageGetAtomVectNull = "Error in base::Molecule::GetAtomVect: atomVect is NULL.\n";
+   this->errorMessageGetAtomNull = "Error in base::Molecule::GetAtom: atomVect is NULL.\n";
+   this->errorMessageAddAtomNull = "Error in base::Molecule::AddAtom: atomVect is NULL.\n";
+   this->errorMessageGetNumberAtomsNull = "Error in base::Molecule::GetNumberAtoms: atomVect is NULL.\n";
    this->errorMessageGetXyzCOCNull = "Error in base::Molecule::GetXyzCOC: xyzCOC is NULL.\n";
    this->errorMessageGetXyzCOMNull = "Error in base::Molecule::GetXyzCOM: xyzCOM is NULL.\n";
    this->messageTotalNumberAOs = "\tTotal number of valence AOs: ";
@@ -161,15 +164,17 @@ void Molecule::SetMessages(){
 }
 
 void Molecule::AddAtom(Atom* atom){
+#ifdef MOLDS_DBG
+   if(this->atomVect==NULL) throw MolDSException(this->errorMessageAddAtomNull);
+#endif
    this->atomVect->push_back(atom);
+   
 }
 
 double* Molecule::GetXyzCOM() const{
-   if(this->xyzCOM==NULL){
-      stringstream ss;
-      ss << this->errorMessageGetXyzCOMNull;
-      throw MolDSException(ss.str());
-   }
+#ifdef MOLDS_DBG
+   if(this->xyzCOM==NULL) throw MolDSException(this->errorMessageGetXyzCOMNull);
+#endif
    return this->xyzCOM;
 }
 
@@ -181,11 +186,9 @@ double* Molecule::GetXyzCOM(){
 }
 
 double* Molecule::GetXyzCOC() const{
-   if(this->xyzCOC==NULL){
-      stringstream ss;
-      ss << this->errorMessageGetXyzCOCNull;
-      throw MolDSException(ss.str());
-   }
+#ifdef MOLDS_DBG
+   if(this->xyzCOC==NULL) throw MolDSException(this->errorMessageGetXyzCOCNull);
+#endif
    return this->xyzCOC;
 }
 

--- a/src/base/Molecule.h
+++ b/src/base/Molecule.h
@@ -26,8 +26,18 @@ public:
    explicit Molecule(const Molecule& rhs);
    Molecule& operator=(const Molecule& rhs);
    ~Molecule();
-   int GetNumberAtoms() const{return this->atomVect->size();}
-   MolDS_base_atoms::Atom* GetAtom(int atomIndex) const{return (*this->atomVect)[atomIndex];}
+   inline int GetNumberAtoms() const{
+#ifdef MOLDS_DBG
+      if(this->atomVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetNumberAtomsNull);
+#endif
+      return this->atomVect->size();
+   }
+   inline MolDS_base_atoms::Atom* GetAtom(int atomIndex) const{
+#ifdef MOLDS_DBG
+      if(this->atomVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetAtomNull);
+#endif
+      return (*this->atomVect)[atomIndex];
+   }
    void AddAtom(MolDS_base_atoms::Atom* atom);
    double* GetXyzCOM() const;
    double* GetXyzCOM();
@@ -85,7 +95,9 @@ private:
                                  double rotatingAngle, 
                                  MolDS_base::EularAngle rotatingEularAngles)const;
    void OutputTranslatingConditions(double const* translatingDifference) const;
-   std::string errorMessageGetAtomVectNull;
+   std::string errorMessageGetAtomNull;
+   std::string errorMessageAddAtomNull;
+   std::string errorMessageGetNumberAtomsNull;
    std::string errorMessageGetXyzCOCNull;
    std::string errorMessageGetXyzCOMNull;
    std::string messageTotalNumberAOs;

--- a/src/base/Parameters.cpp
+++ b/src/base/Parameters.cpp
@@ -452,11 +452,9 @@ EularAngle Parameters::GetRotatingEularAngles() const{
 
 // methods for MOPlot
 vector<int>* Parameters::GetIndecesMOPlot() const{
-   if(this->indecesMOPlot==NULL){
-      stringstream ss;
-      ss << this->errorMessageGetIndecesMOPlotNull; 
-      throw MolDSException(ss.str());
-   }
+#ifdef MOLDS_DBG
+   if(this->indecesMOPlot==NULL) throw MolDSException(this->errorMessageGetIndecesMOPlotNull);
+#endif
    return this->indecesMOPlot;
 }
 
@@ -501,11 +499,9 @@ void Parameters::SetFrameLengthMOPlot(double lx, double ly, double lz){
 
 // methods for HolePlot
 vector<int>* Parameters::GetElecIndecesHolePlot() const{
-   if(this->elecIndecesHolePlot==NULL){
-      stringstream ss;
-      ss << this->errorMessageGetIndecesHolePlotNull; 
-      throw MolDSException(ss.str());
-   }
+#ifdef MOLDS_DBG
+   if(this->elecIndecesHolePlot==NULL) throw MolDSException(this->errorMessageGetIndecesHolePlotNull);
+#endif
    return this->elecIndecesHolePlot;
 }
 
@@ -550,11 +546,9 @@ void Parameters::SetFrameLengthHolePlot(double lx, double ly, double lz){
 
 // methods for ParticlePlot
 vector<int>* Parameters::GetElecIndecesParticlePlot() const{
-   if(this->elecIndecesParticlePlot==NULL){
-      stringstream ss;
-      ss << this->errorMessageGetIndecesParticlePlotNull; 
-      throw MolDSException(ss.str());
-   }
+#ifdef MOLDS_DBG
+   if(this->elecIndecesParticlePlot==NULL) throw MolDSException(this->errorMessageGetIndecesParticlePlotNull);
+#endif
    return this->elecIndecesParticlePlot;
 }
 
@@ -687,11 +681,9 @@ void Parameters::SetRequiresAllTransitionDipoleMomentsCIS(bool requiresAllTransi
 }
 
 vector<int>* Parameters::GetElectronicStateIndecesMullikenCIS() const{
-   if(this->electronicStateIndecesMullikenCIS==NULL){
-      stringstream ss;
-      ss << this->errorMessageGetElectronicStateIndecesMullikenCISNull; 
-      throw MolDSException(ss.str());
-   }
+#ifdef MOLDS_DBG
+   if(this->electronicStateIndecesMullikenCIS==NULL) throw MolDSException(this->errorMessageGetElectronicStateIndecesMullikenCISNull);
+#endif
    return this->electronicStateIndecesMullikenCIS;
 }
 

--- a/src/base/atoms/Atom.cpp
+++ b/src/base/atoms/Atom.cpp
@@ -130,7 +130,9 @@ void Atom::SetMessages(){
       = "Error in base_atoms::Atom::GetPm3PddgParameterDa: Bad index for parameter Da(daIndex). Only 0, and 1 are permitted.\n";
    this->errorMessageDaIndex  = "daIndex = ";
    this->errorMessageGetXyzCoordinatesNull = "Error in base_atoms::Atom::GetXyz: xyz is NULL\n";
+   this->errorMessageSetXyzCoordinatesNull = "Error in base_atoms::Atom::SetXyz: xyz is NULL\n";
    this->errorMessageGetPxyzMomentaNull  = "Error in base_atoms::Atom::GetPxyz: pxyz is NULL\n";
+   this->errorMessageSetPxyzMomentaNull  = "Error in base_atoms::Atom::SetPxyz: pxyz is NULL\n";
 }
 
 AtomType Atom::GetAtomType() const{
@@ -146,33 +148,31 @@ double Atom::GetCoreMass() const{
 }
 
 double* Atom::GetXyz() const{
-   if(this->xyz==NULL){
-      stringstream ss;
-      ss << this->errorMessageGetXyzCoordinatesNull;
-      throw MolDSException(ss.str());
-   }
+#ifdef MOLDS_DBG
+   if(this->xyz==NULL) throw MolDSException(this->errorMessageGetXyzCoordinatesNull);
+#endif
    return this->xyz;
 }
 
 double* Atom::GetPxyz() const{
-   if(this->pxyz==NULL){
-      stringstream ss;
-      ss << this->errorMessageGetPxyzMomentaNull;
-      throw MolDSException(ss.str());
-   }
+#ifdef MOLDS_DBG
+   if(this->pxyz==NULL) throw MolDSException(this->errorMessageGetPxyzMomentaNull);
+#endif
    return this->pxyz;
 }
 
 void Atom::SetXyz(double x, double y, double z) const{
-   xyz[0]= x;
-   xyz[1]= y;
-   xyz[2]= z;
+#ifdef MOLDS_DBG
+   if(this->xyz==NULL) throw MolDSException(this->errorMessageSetXyzCoordinatesNull);
+#endif
+   xyz[0]= x; xyz[1]= y; xyz[2]= z;
 }
 
 void Atom::SetPxyz(double px, double py, double pz) const{
-   pxyz[0]= px;
-   pxyz[1]= py;
-   pxyz[2]= pz;
+#ifdef MOLDS_DBG
+   if(this->pxyz==NULL) throw MolDSException(this->errorMessageSetPxyzMomentaNull);
+#endif
+   pxyz[0]= px; pxyz[1]= py; pxyz[2]= pz;
 }
 
 int Atom::GetValenceSize() const{

--- a/src/base/atoms/Atom.h
+++ b/src/base/atoms/Atom.h
@@ -273,7 +273,9 @@ private:
    std::string errorMessageGetNddoHspBadTheory;
    std::string errorMessageGetNddoHppBadTheory;
    std::string errorMessageGetXyzCoordinatesNull;
+   std::string errorMessageSetXyzCoordinatesNull;
    std::string errorMessageGetPxyzMomentaNull;
+   std::string errorMessageSetPxyzMomentaNull;
    void SetMessages();
    double GetRealAngularPartAO(double theta, 
                                double phi, 

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -115,9 +115,9 @@ Cndo2::~Cndo2(){
                                               this->molecule->GetTotalNumberAOs(),
                                               this->molecule->GetTotalNumberAOs());
    MallocerFreer::GetInstance()->Free<double>(&this->cartesianMatrix, 
+                                              CartesianType_end,
                                               this->molecule->GetTotalNumberAOs(),
-                                              this->molecule->GetTotalNumberAOs(),
-                                              CartesianType_end);
+                                              this->molecule->GetTotalNumberAOs());
    int electronicTransitionDipoleMomentsDim = 1;
    if(Parameters::GetInstance()->RequiresCIS()){
       electronicTransitionDipoleMomentsDim += Parameters::GetInstance()->GetNumberExcitedStatesCIS();
@@ -180,8 +180,8 @@ void Cndo2::SetMessages(){
       = "Error in cndo::Cndo2::RotateDiatmicOverlapAOsToSpaceFrame: rotatingMatrix is NULL.\n";
    this->errorMessageSetOverlapAOsElementNullDiaMatrix 
       = "Error in cndo::Cndo2::SetOverlapAOsElement: diatomicOverlapAOs is NULL.\n";
-   this->errorMessageGetElectronicTransitionDipoleMomentBadState
-      = "Error in cndo::Cndo2::GetElectronicTransitionDipoleMoment: Bad eigen state is set. In SCF module, the transition dipole moment of only between ground states can be calculated. Note taht state=0 means the ground state and other state = i means the i-th excited state in below.\n";
+   this->errorMessageCalcElectronicTransitionDipoleMomentBadState
+      = "Error in cndo::Cndo2::CalcElectronicTransitionDipoleMoment: Bad eigen state is set. In SCF module, the transition dipole moment of only between ground states can be calculated. Note taht state=0 means the ground state and other state = i means the i-th excited state in below.\n";
    this->errorMessageCalcFrequenciesNormalModesBadTheory
       = "Error in cndo::Cndo2::CalcFrequenciesNormalModesBadTheory: CNDO2 is not supported for frequency (normal mode) analysis.\n";
    this->errorMessageCalcOverlapAOsDifferentConfigurationsDiffAOs
@@ -279,9 +279,9 @@ void Cndo2::SetMolecule(Molecule* molecule){
                                                 this->molecule->GetTotalNumberAOs(), 
                                                 this->molecule->GetTotalNumberAOs());
    MallocerFreer::GetInstance()->Malloc<double>(&this->cartesianMatrix, 
+                                                CartesianType_end,
                                                 this->molecule->GetTotalNumberAOs(), 
-                                                this->molecule->GetTotalNumberAOs(),
-                                                CartesianType_end);
+                                                this->molecule->GetTotalNumberAOs());
    int electronicTransitionDipoleMomentsDim = 1;
    if(Parameters::GetInstance()->RequiresCIS()){
       electronicTransitionDipoleMomentsDim += Parameters::GetInstance()->GetNumberExcitedStatesCIS();
@@ -505,12 +505,11 @@ double Cndo2::GetDiatomVdWCorrection2ndDerivative(int indexAtomA,
 void Cndo2::DoSCF(bool requiresGuess){
    this->OutputLog(this->messageStartSCF);
    double ompStartTime = omp_get_wtime();
-
+#ifdef MOLDS_DBG
    if(this->molecule == NULL){
-      stringstream ss;
-      ss << this->errorMessageMoleculeNotSet;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageMoleculeNotSet);
    }
+#endif
 
    // temporary matrices for scf
    double**  oldOrbitalElectronPopulation = NULL;
@@ -697,11 +696,11 @@ double Cndo2::GetElectronicEnergy(int elecState) const{
       return this->elecSCFEnergy;
    }
    else{
+#ifdef MOLDS_DBG
       if(this->excitedEnergies == NULL){
-         stringstream ss;
-         ss << this->errorMessageGetElectronicEnergyNULLCISEnergy;
-         throw MolDSException(ss.str());
+         throw MolDSException(this->errorMessageGetElectronicEnergyNULLCISEnergy);
       }
+#endif
       int numberExcitedStates = Parameters::GetInstance()->GetNumberExcitedStatesCIS();
       if(numberExcitedStates < elecState){
          stringstream ss;
@@ -1655,61 +1654,56 @@ void Cndo2::CalcElectronicDipoleMomentGroundState(double*** electronicTransition
                                                   double const* const* orbitalElectronPopulation,
                                                   double const* const* overlapAOs) const{
    int groundState = 0;
-   for(int axis=0; axis<CartesianType_end; axis++){
-      electronicTransitionDipoleMoments[groundState][groundState][axis] = this->GetElectronicTransitionDipoleMoment(
-                                                                                groundState,
-                                                                                groundState,
-                                                                                static_cast<CartesianType>(axis),
-                                                                                NULL,
-                                                                                NULL,
-                                                                                cartesianMatrix,
-                                                                                molecule,
-                                                                                orbitalElectronPopulation,
-                                                                                overlapAOs,
-                                                                                NULL);
-   }
+   this->CalcElectronicTransitionDipoleMoment(electronicTransitionDipoleMoments[groundState][groundState],
+                                              groundState,
+                                              groundState,
+                                              NULL,
+                                              NULL,
+                                              cartesianMatrix,
+                                              molecule,
+                                              orbitalElectronPopulation,
+                                              overlapAOs,
+                                              NULL);
 }
 
-double Cndo2::GetElectronicTransitionDipoleMoment(int to, int from, CartesianType axis,
-                                                  double const* const* fockMatrix,
-                                                  double const* const* matrixCIS,
-                                                  double const* const* const* cartesianMatrix,
-                                                  const MolDS_base::Molecule& molecule, 
-                                                  double const* const* orbitalElectronPopulation,
-                                                  double const* const* overlapAOs,
-                                                  double const* groundStateDipole) const{
+void Cndo2::CalcElectronicTransitionDipoleMoment(double* transitionDipoleMoment,
+                                                 int to, int from,
+                                                 double const* const* fockMatrix,
+                                                 double const* const* matrixCIS,
+                                                 double const* const* const* cartesianMatrix,
+                                                 const MolDS_base::Molecule& molecule, 
+                                                 double const* const* orbitalElectronPopulation,
+                                                 double const* const* overlapAOs,
+                                                 double const* groundStateDipole) const{
    int groundState = 0;
    if(from == groundState && to == groundState){
-      double value = 0.0;
+      double const* centerOfDipole = molecule.GetXyzCOC();
       int totalAONumber = molecule.GetTotalNumberAOs();
-      stringstream ompErrors;
-#pragma omp parallel for reduction(+:value) schedule(auto) 
-      for(int mu=0; mu<totalAONumber; mu++){
-         try{
-            double threadValue = 0.0;
-            for(int nu=0; nu<totalAONumber; nu++){
-               threadValue -= orbitalElectronPopulation[mu][nu]
-                             *(cartesianMatrix[mu][nu][axis]-molecule.GetXyzCOC()[axis]*overlapAOs[mu][nu]);
-            }
-            value += threadValue;
-         }
-         catch(MolDSException ex){
-#pragma omp critical
-            ompErrors << ex.what() << endl ;
-         }
-      }
-      // Exception throwing for omp-region
-      if(!ompErrors.str().empty()){
-         throw MolDSException(ompErrors.str());
-      }
-      return value;
+      transitionDipoleMoment[XAxis] = 0.0;
+      transitionDipoleMoment[YAxis] = 0.0;
+      transitionDipoleMoment[ZAxis] = 0.0;
+      transitionDipoleMoment[XAxis] -= MolDS_wrappers::Blas::GetInstance()->Ddot(totalAONumber*totalAONumber,
+                                                                                 &orbitalElectronPopulation[0][0],
+                                                                                 &cartesianMatrix[XAxis][0][0]);
+      transitionDipoleMoment[YAxis] -= MolDS_wrappers::Blas::GetInstance()->Ddot(totalAONumber*totalAONumber,
+                                                                                 &orbitalElectronPopulation[0][0],
+                                                                                 &cartesianMatrix[YAxis][0][0]);
+      transitionDipoleMoment[ZAxis] -= MolDS_wrappers::Blas::GetInstance()->Ddot(totalAONumber*totalAONumber,
+                                                                                 &orbitalElectronPopulation[0][0],
+                                                                                 &cartesianMatrix[ZAxis][0][0]);
+      // set orign of dipole
+      double temp = MolDS_wrappers::Blas::GetInstance()->Ddot(totalAONumber*totalAONumber,
+                                                              &orbitalElectronPopulation[0][0],
+                                                              &overlapAOs[0][0]);
+      transitionDipoleMoment[XAxis] += centerOfDipole[XAxis]*temp;
+      transitionDipoleMoment[YAxis] += centerOfDipole[YAxis]*temp;
+      transitionDipoleMoment[ZAxis] += centerOfDipole[ZAxis]*temp;
    }
    else{
       stringstream ss;
-      ss << this->errorMessageGetElectronicTransitionDipoleMomentBadState;
+      ss << this->errorMessageCalcElectronicTransitionDipoleMomentBadState;
       ss << this->errorMessageFromState << from << endl;
       ss << this->errorMessageToState << to << endl;
-      ss << this->errorMessageCartesianType << CartesianTypeStr(axis) << endl;
       throw MolDSException(ss.str());
    }
 }
@@ -1735,15 +1729,10 @@ void Cndo2::CalcCartesianMatrixByGTOExpansion(double*** cartesianMatrix,
                for(int b=0; b<atomB.GetValenceSize(); b++){
                   int mu = firstAOIndexAtomA + a;      
                   int nu = firstAOIndexAtomB + b;      
-                  for(int i=0; i<CartesianType_end; i++){
-                     double value = this->GetCartesianMatrixElementByGTOExpansion(atomA, 
-                                                                                  a, 
-                                                                                  atomB, 
-                                                                                  b, 
-                                                                                  static_cast<CartesianType>(i),
-                                                                                  stonG);
-                     cartesianMatrix[mu][nu][i] = value;
-                  }
+                  this->CalcCartesianMatrixElementsByGTOExpansion(cartesianMatrix[XAxis][mu][nu], 
+                                                                  cartesianMatrix[YAxis][mu][nu],
+                                                                  cartesianMatrix[ZAxis][mu][nu],
+                                                                  atomA, a, atomB, b, stonG);
                }
             }
             
@@ -1758,28 +1747,19 @@ void Cndo2::CalcCartesianMatrixByGTOExpansion(double*** cartesianMatrix,
    if(!ompErrors.str().empty()){
       throw MolDSException(ompErrors.str());
    }
-   /* 
-   this->OutputLog("cartesian matrix\n"); 
-   for(int o=0; o<molecule.GetTotalNumberAOs(); o++){
-      for(int p=0; p<molecule.GetTotalNumberAOs(); p++){
-         for(int i=0; i<CartesianType_end; i++){
-            this->OutputLog(boost::format("%lf\t") % cartesianMatrix[o][p][i]);
-         }
-         this->OutputLog("\n");
-      }
-      this->OutputLog("\n");
-   }
-   this->OutputLog("\n");
-   */
 }
 
 // Calculate elements of Cartesian matrix between atomic orbitals. 
 // The analytic Cartesian matrix is calculated with Gaussian expansion technique written in [DY_1977]
-double Cndo2::GetCartesianMatrixElementByGTOExpansion(const Atom& atomA, int valenceIndexA, 
+void Cndo2::CalcCartesianMatrixElementsByGTOExpansion(double& xComponent,
+                                                      double& yComponent,
+                                                      double& zComponent,
+                                                      const Atom& atomA, int valenceIndexA, 
                                                       const Atom& atomB, int valenceIndexB,
-                                                      CartesianType axis,
                                                       STOnGType stonG) const{
-   double value = 0.0;
+   xComponent=0.0;
+   yComponent=0.0;
+   zComponent=0.0;
    ShellType shellTypeA = atomA.GetValenceShellType();
    ShellType shellTypeB = atomB.GetValenceShellType();
    OrbitalType valenceOrbitalA = atomA.GetValence(valenceIndexA);
@@ -1795,7 +1775,10 @@ double Cndo2::GetCartesianMatrixElementByGTOExpansion(const Atom& atomA, int val
    double Rab = sqrt( pow(atomA.GetXyz()[XAxis]-atomB.GetXyz()[XAxis], 2.0) 
                      +pow(atomA.GetXyz()[YAxis]-atomB.GetXyz()[YAxis], 2.0) 
                      +pow(atomA.GetXyz()[ZAxis]-atomB.GetXyz()[ZAxis], 2.0) );
-   double temp = 0.0;
+   double temp  = 0.0;
+   double tempX = 0.0;
+   double tempY = 0.0;
+   double tempZ = 0.0;
    for(int i=0; i<=stonG; i++){
       for(int j=0; j<=stonG; j++){
          temp = GTOExpansionSTO::GetInstance()->GetCoefficient(stonG, 
@@ -1816,20 +1799,23 @@ double Cndo2::GetCartesianMatrixElementByGTOExpansion(const Atom& atomA, int val
                                                                          shellTypeB, 
                                                                          valenceOrbitalB, 
                                                                          j);
-         temp *= this->GetGaussianCartesianMatrix(atomA.GetAtomType(), 
-                                                  valenceOrbitalA, 
-                                                  gaussianExponentA, 
-                                                  atomA.GetXyz(),
-                                                  atomB.GetAtomType(), 
-                                                  valenceOrbitalB, 
-                                                  gaussianExponentB,
-                                                  atomB.GetXyz(), 
+         tempX = this->GetGaussianCartesianMatrix(atomA.GetAtomType(), valenceOrbitalA, gaussianExponentA, atomA.GetXyz(),
+                                                  atomB.GetAtomType(), valenceOrbitalB, gaussianExponentB, atomB.GetXyz(), 
                                                   Rab,
-                                                  axis);
-         value += temp;
+                                                  XAxis);
+         tempY = this->GetGaussianCartesianMatrix(atomA.GetAtomType(), valenceOrbitalA, gaussianExponentA, atomA.GetXyz(),
+                                                  atomB.GetAtomType(), valenceOrbitalB, gaussianExponentB, atomB.GetXyz(), 
+                                                  Rab,
+                                                  YAxis);
+         tempZ = this->GetGaussianCartesianMatrix(atomA.GetAtomType(), valenceOrbitalA, gaussianExponentA, atomA.GetXyz(),
+                                                  atomB.GetAtomType(), valenceOrbitalB, gaussianExponentB, atomB.GetXyz(), 
+                                                  Rab,
+                                                  ZAxis);
+         xComponent += temp*tempX;
+         yComponent += temp*tempY;
+         zComponent += temp*tempZ;
       }
    }
-   return value;
 }
 
 // calculate gaussian Caretesian integrals. 
@@ -3511,11 +3497,11 @@ void Cndo2::CalcOverlapAOsWithAnotherConfiguration(double** overlapAOs,
       ss << this->errorMessageRhs << rhsMolecule->GetNumberAtoms() << endl;
       throw MolDSException(ss.str());
    }
+#ifdef MOLDS_DBG
    if(overlapAOs == NULL){
-      stringstream ss;
-      ss << this->errorMessageCalcOverlapAOsDifferentConfigurationsOverlapAOsNULL;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageCalcOverlapAOsDifferentConfigurationsOverlapAOsNULL);
    }
+#endif
    
    int totalAONumber = lhsMolecule.GetTotalNumberAOs();
    int totalAtomNumber = lhsMolecule.GetNumberAtoms();
@@ -5036,11 +5022,11 @@ double Cndo2::GetGaussianOverlapAOs1stDerivative(AtomType atomTypeA,
 void Cndo2::CalcRotatingMatrix(double** rotatingMatrix, 
                                const Atom& atomA, 
                                const Atom& atomB) const{
+#ifdef MOLDS_DBG
    if(rotatingMatrix==NULL){
-      stringstream ss;
-      ss << this->errorMessageCalcRotatingMatrixNullRotMatrix;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageCalcRotatingMatrixNullRotMatrix);
    }
+#endif
    MallocerFreer::GetInstance()->Initialize<double>(rotatingMatrix,  OrbitalType_end, OrbitalType_end);
 
    double x = atomB.GetXyz()[0] - atomA.GetXyz()[0];
@@ -5338,11 +5324,11 @@ void Cndo2::CalcRotatingMatrix2ndDerivatives(double**** rotMat2ndDerivatives,
 void Cndo2::CalcDiatomicOverlapAOsInDiatomicFrame(double** diatomicOverlapAOs, 
                                                   const Atom& atomA, 
                                                   const Atom& atomB) const{
+#ifdef MOLDS_DBG
    if(diatomicOverlapAOs==NULL){
-      stringstream ss;
-      ss << this->errorMessageCalDiaOverlapAOsDiaFrameNullMatrix;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageCalDiaOverlapAOsDiaFrameNullMatrix);
    }
+#endif
    int na = atomA.GetValenceShellType() + 1;
    int nb = atomB.GetValenceShellType() + 1;
    int m = 0;
@@ -5625,16 +5611,14 @@ void Cndo2::CalcDiatomicOverlapAOs2ndDerivativeInDiatomicFrame(double** diatomic
 // see (B.63) in Pople book.
 void Cndo2::RotateDiatmicOverlapAOsToSpaceFrame(double** diatomicOverlapAOs, 
                                                 double const* const* rotatingMatrix) const{
+#ifdef MOLDS_DBG
    if(diatomicOverlapAOs==NULL){
-      stringstream ss;
-      ss << this->errorMessageRotDiaOverlapAOsToSpaceFrameNullDiaMatrix;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageRotDiaOverlapAOsToSpaceFrameNullDiaMatrix);
    }
    if(rotatingMatrix==NULL){
-      stringstream ss;
-      ss << this->errorMessageRotDiaOverlapAOsToSpaceFrameNullRotMatrix;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageRotDiaOverlapAOsToSpaceFrameNullRotMatrix);
    }
+#endif
    double** oldDiatomicOverlapAOs = NULL;
    try{
       MallocerFreer::GetInstance()->Malloc<double>(&oldDiatomicOverlapAOs, OrbitalType_end, OrbitalType_end);
@@ -5679,11 +5663,11 @@ void Cndo2::SetOverlapAOsElement(double** overlapAOs,
                                  const Atom& atomA, 
                                  const Atom& atomB,
                                  bool isSymmetricOverlapAOs) const{
+#ifdef MOLDS_DBG
    if(diatomicOverlapAOs==NULL){
-      stringstream ss;
-      ss << this->errorMessageSetOverlapAOsElementNullDiaMatrix;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageSetOverlapAOsElementNullDiaMatrix);
    }
+#endif
 
    int firstAOIndexAtomA = atomA.GetFirstAOIndex();
    int firstAOIndexAtomB = atomB.GetFirstAOIndex();

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -76,7 +76,7 @@ protected:
    std::string errorMessageGetElectronicEnergyEnergyNotCalculated;
    std::string errorMessageGetElectronicEnergyNumberCISStates;
    std::string errorMessageGetElectronicEnergySetElecState;
-   std::string errorMessageGetElectronicTransitionDipoleMomentBadState;
+   std::string errorMessageCalcElectronicTransitionDipoleMomentBadState;
    std::string errorMessageCalcFrequenciesNormalModesBadTheory;
    std::string errorMessageFromState;
    std::string errorMessageToState;
@@ -119,14 +119,15 @@ protected:
    virtual void CalcSCFProperties();
    virtual void CalcCISProperties();
    virtual void CalcNormalModes(double** normalModes, double* normalForceConstants, const MolDS_base::Molecule& molecule) const;
-   virtual double GetElectronicTransitionDipoleMoment(int to, int from, MolDS_base::CartesianType axis,
-                                                      double const* const* fockMatrix,
-                                                      double const* const* matrixCIS,
-                                                      double const* const* const* cartesianMatrix,
-                                                      const MolDS_base::Molecule& molecule, 
-                                                      double const* const* orbitalElectronPopulation,
-                                                      double const* const* overlapAOs,
-                                                      double const* groundStateDipole) const;
+   virtual void CalcElectronicTransitionDipoleMoment(double* transitionDipoleMoment,
+                                                     int to, int from,
+                                                     double const* const* fockMatrix,
+                                                     double const* const* matrixCIS,
+                                                     double const* const* const* cartesianMatrix,
+                                                     const MolDS_base::Molecule& molecule, 
+                                                     double const* const* orbitalElectronPopulation,
+                                                     double const* const* overlapAOs,
+                                                     double const* groundStateDipole) const;
    double GetBondingAdjustParameterK(MolDS_base::ShellType shellA, 
                                      MolDS_base::ShellType shellB) const;
    virtual double GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const;
@@ -343,11 +344,13 @@ private:
    void CalcCartesianMatrixByGTOExpansion(double*** cartesianMatrix,
                                           const MolDS_base::Molecule& molecule, 
                                           MolDS_base::STOnGType stonG) const; 
-   double GetCartesianMatrixElementByGTOExpansion(const MolDS_base_atoms::Atom& atomA, 
+   void CalcCartesianMatrixElementsByGTOExpansion(double& xComponent,
+                                                  double& yComponent,
+                                                  double& zComponent,
+                                                  const MolDS_base_atoms::Atom& atomA, 
                                                   int valenceIndexA, 
                                                   const MolDS_base_atoms::Atom& atomB, 
                                                   int valenceIndexB,
-                                                  MolDS_base::CartesianType axis,
                                                   MolDS_base::STOnGType stonG) const;
    double GetGaussianCartesianMatrix(MolDS_base::AtomType atomTypeA, 
                                      MolDS_base::OrbitalType valenceOrbitalA, 

--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -748,7 +748,39 @@ void Mndo::CalcCISMatrix(double** matrixCIS) const{
                   if(A!=B){
                      for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
                         for(int nu=mu; nu<=lastAOIndexA; nu++){
+                           double tmpMuNu01 = 2.0*fockMatrix[moA][mu]
+                                                 *fockMatrix[moI][nu];
+                           double tmpMuNu02 = 2.0*fockMatrix[moJ][mu]
+                                                 *fockMatrix[moB][nu];
+                           double tmpMuNu03 = fockMatrix[moA][mu]
+                                             *fockMatrix[moB][nu];
+                           double tmpMuNu04 = fockMatrix[moI][mu]
+                                             *fockMatrix[moJ][nu];
+                           double tmpMuNu09 = 2.0*fockMatrix[moI][mu]
+                                                 *fockMatrix[moA][nu];
+                           double tmpMuNu10 = 2.0*fockMatrix[moB][mu]
+                                                 *fockMatrix[moJ][nu];
+                           double tmpMuNu11 = fockMatrix[moB][mu]
+                                             *fockMatrix[moA][nu];
+                           double tmpMuNu12 = fockMatrix[moJ][mu]
+                                             *fockMatrix[moI][nu];
                            for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
+                              double tmpMuNuLamda01 = tmpMuNu01*fockMatrix[moJ][lambda];
+                              double tmpMuNuLamda02 = tmpMuNu02*fockMatrix[moA][lambda];
+                              double tmpMuNuLamda03 = tmpMuNu03*fockMatrix[moI][lambda];
+                              double tmpMuNuLamda04 = tmpMuNu04*fockMatrix[moA][lambda];
+                              double tmpMuNuLamda05 = tmpMuNu01*fockMatrix[moB][lambda];
+                              double tmpMuNuLamda06 = tmpMuNu02*fockMatrix[moI][lambda];
+                              double tmpMuNuLamda07 = tmpMuNu03*fockMatrix[moJ][lambda];
+                              double tmpMuNuLamda08 = tmpMuNu04*fockMatrix[moB][lambda];
+                              double tmpMuNuLamda09 = tmpMuNu09*fockMatrix[moJ][lambda];
+                              double tmpMuNuLamda10 = tmpMuNu10*fockMatrix[moA][lambda];
+                              double tmpMuNuLamda11 = tmpMuNu11*fockMatrix[moI][lambda];
+                              double tmpMuNuLamda12 = tmpMuNu12*fockMatrix[moA][lambda];
+                              double tmpMuNuLamda13 = tmpMuNu09*fockMatrix[moB][lambda];
+                              double tmpMuNuLamda14 = tmpMuNu10*fockMatrix[moI][lambda];
+                              double tmpMuNuLamda15 = tmpMuNu11*fockMatrix[moJ][lambda];
+                              double tmpMuNuLamda16 = tmpMuNu12*fockMatrix[moB][lambda];
                               for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
                                  OrbitalType orbitalSigma = atomB.GetValence(sigma-firstAOIndexB);
                                  gamma = this->twoElecTwoCore[A]
@@ -758,75 +790,27 @@ void Mndo::CalcCISMatrix(double** matrixCIS) const{
                                                              [lambda-firstAOIndexB]
                                                              [sigma-firstAOIndexB];
    
-                                 value += 2.0*gamma*fockMatrix[moA][mu]
-                                                   *fockMatrix[moI][nu]
-                                                   *fockMatrix[moJ][lambda]
-                                                      *fockMatrix[moB][sigma];
-                                 value += 2.0*gamma*fockMatrix[moA][lambda]
-                                                   *fockMatrix[moI][sigma]
-                                                   *fockMatrix[moJ][mu]
-                                                   *fockMatrix[moB][nu];
-                                 value -= gamma*fockMatrix[moA][mu]
-                                               *fockMatrix[moB][nu]
-                                               *fockMatrix[moI][lambda]
-                                               *fockMatrix[moJ][sigma];
-                                 value -= gamma*fockMatrix[moA][lambda]
-                                               *fockMatrix[moB][sigma]
-                                               *fockMatrix[moI][mu]
-                                               *fockMatrix[moJ][nu];
+                                 value += gamma*tmpMuNuLamda01*fockMatrix[moB][sigma];
+                                 value += gamma*tmpMuNuLamda02*fockMatrix[moI][sigma];
+                                 value -= gamma*tmpMuNuLamda03*fockMatrix[moJ][sigma];
+                                 value -= gamma*tmpMuNuLamda04*fockMatrix[moB][sigma];
                                  if(lambda != sigma){
-                                    value += 2.0*gamma*fockMatrix[moA][mu]
-                                                      *fockMatrix[moI][nu]
-                                                      *fockMatrix[moJ][sigma]
-                                                      *fockMatrix[moB][lambda];
-                                    value += 2.0*gamma*fockMatrix[moA][sigma]
-                                                      *fockMatrix[moI][lambda]
-                                                      *fockMatrix[moJ][mu]
-                                                      *fockMatrix[moB][nu];
-                                    value -= gamma*fockMatrix[moA][mu]
-                                                  *fockMatrix[moB][nu]
-                                                  *fockMatrix[moI][sigma]
-                                                  *fockMatrix[moJ][lambda];
-                                    value -= gamma*fockMatrix[moA][sigma]
-                                                  *fockMatrix[moB][lambda]
-                                                  *fockMatrix[moI][mu]
-                                                  *fockMatrix[moJ][nu];
+                                    value += gamma*tmpMuNuLamda05*fockMatrix[moJ][sigma];
+                                    value += gamma*tmpMuNuLamda06*fockMatrix[moA][sigma];
+                                    value -= gamma*tmpMuNuLamda07*fockMatrix[moI][sigma];
+                                    value -= gamma*tmpMuNuLamda08*fockMatrix[moA][sigma];
                                  }
                                  if(mu != nu){
-                                    value += 2.0*gamma*fockMatrix[moA][nu]
-                                                      *fockMatrix[moI][mu]
-                                                      *fockMatrix[moJ][lambda]
-                                                      *fockMatrix[moB][sigma];
-                                    value += 2.0*gamma*fockMatrix[moA][lambda]
-                                                      *fockMatrix[moI][sigma]
-                                                      *fockMatrix[moJ][nu]
-                                                      *fockMatrix[moB][mu];
-                                    value -= gamma*fockMatrix[moA][nu]
-                                                  *fockMatrix[moB][mu]
-                                                  *fockMatrix[moI][lambda]
-                                                  *fockMatrix[moJ][sigma];
-                                    value -= gamma*fockMatrix[moA][lambda]
-                                                  *fockMatrix[moB][sigma]
-                                                  *fockMatrix[moI][nu]
-                                                  *fockMatrix[moJ][mu];
+                                    value += gamma*tmpMuNuLamda09*fockMatrix[moB][sigma];
+                                    value += gamma*tmpMuNuLamda10*fockMatrix[moI][sigma];
+                                    value -= gamma*tmpMuNuLamda11*fockMatrix[moJ][sigma];
+                                    value -= gamma*tmpMuNuLamda12*fockMatrix[moB][sigma];
                                  }
                                  if(mu != nu && lambda != sigma){
-                                    value += 2.0*gamma*fockMatrix[moA][nu]
-                                                      *fockMatrix[moI][mu]
-                                                      *fockMatrix[moJ][sigma]
-                                                      *fockMatrix[moB][lambda];
-                                    value += 2.0*gamma*fockMatrix[moA][sigma]
-                                                      *fockMatrix[moI][lambda]
-                                                      *fockMatrix[moJ][nu]
-                                                      *fockMatrix[moB][mu];
-                                    value -= gamma*fockMatrix[moA][nu]
-                                                  *fockMatrix[moB][mu]
-                                                  *fockMatrix[moI][sigma]
-                                                  *fockMatrix[moJ][lambda];
-                                    value -= gamma*fockMatrix[moA][sigma]
-                                                  *fockMatrix[moB][lambda]
-                                                  *fockMatrix[moI][nu]
-                                                  *fockMatrix[moJ][mu];
+                                    value += gamma*tmpMuNuLamda13*fockMatrix[moJ][sigma];
+                                    value += gamma*tmpMuNuLamda14*fockMatrix[moA][sigma];
+                                    value -= gamma*tmpMuNuLamda15*fockMatrix[moI][sigma];
+                                    value -= gamma*tmpMuNuLamda16*fockMatrix[moA][sigma];
                                  }
                               }
                            }
@@ -836,7 +820,13 @@ void Mndo::CalcCISMatrix(double** matrixCIS) const{
                   else{
                      for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
                         for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
+                           double tmpMuNu01 = 2.0*fockMatrix[moA][mu]
+                                                 *fockMatrix[moI][nu];
+                           double tmpMuNu02 = fockMatrix[moA][mu]
+                                             *fockMatrix[moB][nu];
                            for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
+                              double tmpMuNuLamda01 = tmpMuNu01*fockMatrix[moJ][lambda];
+                              double tmpMuNuLamda02 = tmpMuNu02*fockMatrix[moI][lambda];
                               for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
                                  if(mu==nu && lambda==sigma){
                                     OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
@@ -851,14 +841,8 @@ void Mndo::CalcCISMatrix(double** matrixCIS) const{
                                  else{
                                     gamma = 0.0;
                                  }
-                                 value += 2.0*gamma*fockMatrix[moA][mu]
-                                                   *fockMatrix[moI][nu]
-                                                   *fockMatrix[moJ][lambda]
-                                                   *fockMatrix[moB][sigma];
-                                 value -= gamma*fockMatrix[moA][mu]
-                                               *fockMatrix[moB][nu]
-                                               *fockMatrix[moI][lambda]
-                                               *fockMatrix[moJ][sigma];
+                                 value += gamma*tmpMuNuLamda01*fockMatrix[moB][sigma];
+                                 value -= gamma*tmpMuNuLamda02*fockMatrix[moJ][sigma];
                               }  
                            }
                         }
@@ -2947,11 +2931,11 @@ void Mndo::FreeTempMatricesSolveCPHF(double*** matrixCPHF,
 
 // see [PT_1996, PT_1997]
 void Mndo::CalcZMatrixForce(const vector<int>& elecStates){
+#ifdef MOLDS_DBG
    if(this->etaMatrixForce == NULL){
-      stringstream ss;
-      ss << this->errorMessageCalcZMatrixForceEtaNull;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageCalcZMatrixForceEtaNull);
    }
+#endif
    this->CheckZMatrixForce(elecStates); 
 
    // creat MO-index-pair for Q variables. 
@@ -3511,17 +3495,15 @@ void Mndo::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs, dou
 
 void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore, 
                               const Molecule& molecule) const{
+#ifdef MOLDS_DBG
    if(twoElecTwoCore == NULL){
-      stringstream ss;
-      ss << this->errorMessageCalcTwoElecTwoCoreNullMatrix;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageCalcTwoElecTwoCoreNullMatrix);
    }
-   else{
-      MallocerFreer::GetInstance()->Initialize<double>(twoElecTwoCore, 
-                                                       molecule.GetNumberAtoms(),
-                                                       molecule.GetNumberAtoms(),
-                                                       dxy, dxy, dxy, dxy);
-   } 
+#endif
+   MallocerFreer::GetInstance()->Initialize<double>(twoElecTwoCore, 
+                                                    molecule.GetNumberAtoms(),
+                                                    molecule.GetNumberAtoms(),
+                                                    dxy, dxy, dxy, dxy);
 
    stringstream ompErrors;
 #pragma omp parallel
@@ -3586,14 +3568,12 @@ void Mndo::CalcDiatomicTwoElecTwoCore(double**** matrix, int indexAtomA, int ind
       throw MolDSException(ss.str());
    }
 
+#ifdef MOLDS_DBG
    if(matrix == NULL){
-      stringstream ss;
-      ss << this->errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix);
    }
-   else{
-      MallocerFreer::GetInstance()->Initialize<double>(matrix, dxy, dxy, dxy, dxy);
-   } 
+#endif
+   MallocerFreer::GetInstance()->Initialize<double>(matrix, dxy, dxy, dxy, dxy);
 
    // calclation in diatomic frame
    for(int mu=0; mu<atomA.GetValenceSize(); mu++){
@@ -3668,19 +3648,17 @@ void Mndo::CalcDiatomicTwoElecTwoCore1stDerivatives(double***** matrix,
       throw MolDSException(ss.str());
    }
 
+#ifdef MOLDS_DBG
    if(matrix == NULL){
-      stringstream ss;
-      ss << this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix;
-      throw MolDSException(ss.str());
+      throw MolDSException(this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix);
    }
-   else{
-      MallocerFreer::GetInstance()->Initialize<double>(matrix, 
-                                                       dxy, 
-                                                       dxy, 
-                                                       dxy, 
-                                                       dxy, 
-                                                       CartesianType_end);
-   } 
+#endif
+   MallocerFreer::GetInstance()->Initialize<double>(matrix, 
+                                                    dxy, 
+                                                    dxy, 
+                                                    dxy, 
+                                                    dxy, 
+                                                    CartesianType_end);
 
    double**   rotatingMatrix = NULL;
    double***  rotMat1stDerivatives = NULL;
@@ -3766,20 +3744,18 @@ void Mndo::CalcDiatomicTwoElecTwoCore2ndDerivatives(double****** matrix,
       throw MolDSException(ss.str());
    }
 
+#ifdef MOLDS_DBG
    if(matrix == NULL){
-      stringstream ss;
-      ss << this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix;
-      throw MolDSException(ss.str());
-   }
-   else{
-      MallocerFreer::GetInstance()->Initialize<double>(matrix, 
-                                                       dxy, 
-                                                       dxy, 
-                                                       dxy, 
-                                                       dxy, 
-                                                       CartesianType_end,
-                                                       CartesianType_end);
-   } 
+      throw MolDSException(this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix);
+   }
+#endif
+   MallocerFreer::GetInstance()->Initialize<double>(matrix, 
+                                                    dxy, 
+                                                    dxy, 
+                                                    dxy, 
+                                                    dxy, 
+                                                    CartesianType_end,
+                                                    CartesianType_end);
 
    double** rotatingMatrix = NULL;
    double*** rotMat1stDerivatives = NULL;

--- a/src/pm3/Pm3.cpp
+++ b/src/pm3/Pm3.cpp
@@ -25,6 +25,7 @@
 #include<vector>
 #include<boost/format.hpp>
 #include"../base/PrintController.h"
+#include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
 #include"../base/Enums.h"
 #include"../base/EularAngle.h"

--- a/src/pm3/Pm3D.cpp
+++ b/src/pm3/Pm3D.cpp
@@ -25,6 +25,7 @@
 #include<vector>
 #include<boost/format.hpp>
 #include"../base/PrintController.h"
+#include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
 #include"../base/Enums.h"
 #include"../base/EularAngle.h"

--- a/src/pm3/Pm3Pddg.cpp
+++ b/src/pm3/Pm3Pddg.cpp
@@ -25,6 +25,7 @@
 #include<vector>
 #include<boost/format.hpp>
 #include"../base/PrintController.h"
+#include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
 #include"../base/Enums.h"
 #include"../base/EularAngle.h"

--- a/src/wrappers/Blas.cpp
+++ b/src/wrappers/Blas.cpp
@@ -39,6 +39,7 @@
 #include"../base/PrintController.h"
 #include"../base/MolDSException.h"
 #include"../base/Uncopyable.h"
+#include"../base/MallocerFreer.h"
 #include"Blas.h"
 using namespace std;
 using namespace MolDS_base;
@@ -310,6 +311,56 @@ void Blas::Dgemm(bool isColumnMajorMatrixA,
 #endif
 }
 
+// matrixD = matrixA*matrixB*matrixC
+//    matrixA: m*k-matrix (matrixA[m][k] in row-major (C/C++ style))
+//    matrixB: k*l-matrix (matrixB[k][n] in row-major (C/C++ style))
+//    matrixC: k*n-matrix (matrixC[k][n] in row-major (C/C++ style))
+//    matrixD: m*n-matrix (matrixC[m][n] in row-major (C/C++ style))
+void Blas::Dgemmm(molds_blas_int m, molds_blas_int n, molds_blas_int k, molds_blas_int l,
+                  double const* const* matrixA, 
+                  double const* const* matrixB, 
+                  double const* const* matrixC, 
+                  double**             matrixD) const{
+   bool isColumnMajorMatrixA = false; // because, in general, C/C++ style is row-major.
+   bool isColumnMajorMatrixB = false; // because, in general, C/C++ style is row-major.
+   bool isColumnMajorMatrixC = false; // because, in general, C/C++ style is row-major.
+   double alpha=1.0;
+   double beta =0.0;
+   this->Dgemmm(isColumnMajorMatrixA, isColumnMajorMatrixB, isColumnMajorMatrixC, m, n, k, l, alpha, matrixA, matrixB, matrixC, beta, matrixD);
+}
+
+// matrixD = alpha*matrixA*matrixB*matrixC + beta*matrixD
+//    matrixA: m*k-matrix 
+//    matrixB: k*l-matrix
+//    matrixC: l*n-matrix
+//    matrixD: m*n-matrix (matrixC[m][n] in row-major (C/C++ style))
+void Blas::Dgemmm(bool isColumnMajorMatrixA,
+                  bool isColumnMajorMatrixB, 
+                  bool isColumnMajorMatrixC, 
+                  molds_blas_int m, molds_blas_int n, molds_blas_int k, molds_blas_int l,
+                  double alpha,
+                  double const* const* matrixA,
+                  double const* const* matrixB,
+                  double const* const* matrixC,
+                  double beta,
+                  double** matrixD) const{
+   
+   double alphaBC = 1.0;
+   double betaBC  = 0.0;
+   bool isColumnMajorMatrixBC = false;
+   double** matrixBC = NULL;
+   try{
+      MallocerFreer::GetInstance()->Malloc<double>(&matrixBC, k, n); 
+      this->Dgemm(isColumnMajorMatrixB, isColumnMajorMatrixC,  k, n, l, alphaBC, matrixB, matrixC,  betaBC, matrixBC);
+      this->Dgemm(isColumnMajorMatrixA, isColumnMajorMatrixBC, m, n, k, alpha,   matrixA, matrixBC, beta,   matrixD );
+   }
+   catch(MolDSException ex){
+      MallocerFreer::GetInstance()->Free<double>(&matrixBC, k, n); 
+      throw ex;
+   }
+   MallocerFreer::GetInstance()->Free<double>(&matrixBC, k, n); 
+}
+
 // matrixC = matrixA*matrixA^T
 //    matrixA: n*k-matrix
 //    matrixC: n*n-matrix,symmetric (Use the upper triangular part, and copy it to the lower part.)

--- a/src/wrappers/Blas.h
+++ b/src/wrappers/Blas.h
@@ -102,15 +102,30 @@ public:
               double const* const* matrixB,
               double beta,
               double**             matrixC) const;
-              void Dsyrk(molds_blas_int n, molds_blas_int k,
-                         double const *const* matrixA,
-                         double**             matrixC)const;
-              void Dsyrk(molds_blas_int n, molds_blas_int k,
-                         bool isMatrixAColumnMajor,
-                         bool isMatrixATransposed,
-                         bool isLowerTriangularPartMatrixCUsed,
-                         double alpha, double const* const* matrixA,
-                         double beta,  double**             matrixC)const;
+   void Dgemmm(molds_blas_int m, molds_blas_int n, molds_blas_int k, molds_blas_int l,
+               double const* const* matrixA, 
+               double const* const* matrixB, 
+               double const* const* matrixC, 
+               double**             matrixD) const;
+   void Dgemmm(bool isColumnMajorMatrixA, 
+               bool isColumnMajorMatrixB, 
+               bool isColumnMajorMatrixC, 
+               molds_blas_int m, molds_blas_int n, molds_blas_int k, molds_blas_int l,
+               double alpha,
+               double const* const* matrixA,
+               double const* const* matrixB,
+               double const* const* matrixC,
+               double beta,
+               double**             matrixD) const;
+  void Dsyrk(molds_blas_int n, molds_blas_int k,
+             double const *const* matrixA,
+             double**             matrixC)const;
+  void Dsyrk(molds_blas_int n, molds_blas_int k,
+             bool isMatrixAColumnMajor,
+             bool isMatrixATransposed,
+             bool isLowerTriangularPartMatrixCUsed,
+             double alpha, double const* const* matrixA,
+             double beta,  double**             matrixC)const;
 private:
    Blas();
    ~Blas();

--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -129,8 +129,8 @@ void ZindoS::SetMessages(){
       = "Error in zindo::ZindoS::GetElectronicEnergy: Set electronic state is not calculated by CIS.\n";
    this->errorMessageGetElectronicEnergyNULLCISEnergy 
       = "Error in zindo::ZindoS::GetElectronicEnergy: excitedEnergies is NULL\n";
-   this->errorMessageGetElectronicTransitionDipoleMomentBadState
-      = "Error in zindo::ZindoS::GetElectronicTransitionDipoleMoment: Bad eigen state is set to calculate the transition dipole moment. Note taht state=0 means the ground state and other state = i means the i-th excited state in below.\n";
+   this->errorMessageCalcElectronicTransitionDipoleMomentBadState
+      = "Error in zindo::ZindoS::CalcElectronicTransitionDipoleMoment: Bad eigen state is set to calculate the transition dipole moment. Note taht state=0 means the ground state and other state = i means the i-th excited state in below.\n";
    this->errorMessageCalcFrequenciesNormalModesBadTheory
       = "Error in zindo::ZindoS::CalcFrequenciesNormalModesBadTheory: ZINDO/S is not supported for frequency (normal mode) analysis.\n";
    this->messageSCFMetConvergence = "\n\n\n\t\tZINDO/S-SCF met convergence criterion(^^b\n\n\n";
@@ -219,6 +219,7 @@ double ZindoS::GetFockDiagElement(const Atom& atomA,
             OrbitalType orbitalSigma;
             int sigma;
             int atomBNumberValence = atomB.GetValenceSize();
+            double rAB = molecule.GetDistanceAtoms(atomA, atomB);
             for(int i=0; i<atomBNumberValence; i++){
                sigma = i + atomB.GetFirstAOIndex();
                orbitalSigma = atomB.GetValence(i);
@@ -226,10 +227,11 @@ double ZindoS::GetFockDiagElement(const Atom& atomA,
                       *this->GetNishimotoMatagaTwoEleInt(atomA, 
                                                          orbitalMu, 
                                                          atomB, 
-                                                         orbitalSigma);
+                                                         orbitalSigma,
+                                                         rAB);
             }
             temp -= atomB.GetCoreCharge() 
-                   *this->GetNishimotoMatagaTwoEleInt(atomA, s, atomB, s);
+                   *this->GetNishimotoMatagaTwoEleInt(atomA, s, atomB, s, rAB);
          }
       }
       value += temp;
@@ -548,6 +550,13 @@ double ZindoS::GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, const
 double ZindoS::GetNishimotoMatagaTwoEleInt(const Atom& atomA, OrbitalType orbitalA, 
                                            const Atom& atomB, OrbitalType orbitalB) const{
    double r = this->molecule->GetDistanceAtoms(atomA, atomB);
+   return this->GetNishimotoMatagaTwoEleInt(atomA, orbitalA, atomB, orbitalB,r);
+}
+
+// ref. [MN_1957] and (5a) in [AEZ_1986]
+double ZindoS::GetNishimotoMatagaTwoEleInt(const Atom& atomA, OrbitalType orbitalA, 
+                                           const Atom& atomB, OrbitalType orbitalB,
+                                           const double rAB) const{
    double gammaAA;
    if(orbitalA == s || 
       orbitalA == px ||
@@ -596,7 +605,8 @@ double ZindoS::GetNishimotoMatagaTwoEleInt(const Atom& atomA, OrbitalType orbita
       throw MolDSException(ss.str());
    }  
 
-   return this->nishimotoMatagaParamA/( r+this->nishimotoMatagaParamB/(gammaAA+gammaBB) );
+   double gamma=gammaAA+gammaBB;
+   return this->nishimotoMatagaParamA/( rAB+this->nishimotoMatagaParamB/gamma );
 
 }
 
@@ -604,10 +614,23 @@ double ZindoS::GetNishimotoMatagaTwoEleInt(const Atom& atomA, OrbitalType orbita
 // For Nishimoto-Mataga, See ZindoS::GetNishimotoMatagaTwoEleInt 
 // or ref. [MN_1957] and (5a) in [AEZ_1986]
 double ZindoS::GetNishimotoMatagaTwoEleInt1stDerivative(const Atom& atomA, 
-                                                          OrbitalType orbitalA, 
-                                                          const Atom& atomB, 
-                                                          OrbitalType orbitalB,
-                                                          CartesianType axisA) const{
+                                                        OrbitalType orbitalA, 
+                                                        const Atom& atomB, 
+                                                        OrbitalType orbitalB,
+                                                        CartesianType axisA) const{
+   double r = this->molecule->GetDistanceAtoms(atomA, atomB);
+   return this->GetNishimotoMatagaTwoEleInt1stDerivative(atomA, orbitalA, atomB, orbitalB, r, axisA);
+}
+
+// First derivative of Nishimoto-Mataga related to the coordinate of atom A.
+// For Nishimoto-Mataga, See ZindoS::GetNishimotoMatagaTwoEleInt 
+// or ref. [MN_1957] and (5a) in [AEZ_1986]
+double ZindoS::GetNishimotoMatagaTwoEleInt1stDerivative(const Atom& atomA, 
+                                                        OrbitalType orbitalA, 
+                                                        const Atom& atomB, 
+                                                        OrbitalType orbitalB,
+                                                        const double rAB,
+                                                        CartesianType axisA) const{
    double gammaAA;
    if(orbitalA == s || 
       orbitalA == px ||
@@ -656,14 +679,51 @@ double ZindoS::GetNishimotoMatagaTwoEleInt1stDerivative(const Atom& atomA,
       throw MolDSException(ss.str());
    }  
 
-   double r = this->molecule->GetDistanceAtoms(atomA, atomB);
    double dCartesian = atomA.GetXyz()[axisA] - atomB.GetXyz()[axisA];
-   double value = -1.0*dCartesian/r;
+   double value = -1.0*dCartesian/rAB;
    value *= this->nishimotoMatagaParamA;
-   value *= pow( r+this->nishimotoMatagaParamB/(gammaAA+gammaBB) ,-2.0);
+   value *= pow( rAB+this->nishimotoMatagaParamB/(gammaAA+gammaBB) ,-2.0);
    return value;
 }
 
+void ZindoS::CalcNishimotoMatagaMatrix(double**** nishimotoMatagaMatrix, const Molecule& molecule) const{
+   int totalNumberAtoms = molecule.GetNumberAtoms();
+   stringstream ompErrors;
+#pragma omp parallel for schedule(auto)
+   for(int A=0; A<totalNumberAtoms; A++){
+      try{
+         const Atom& atomA = *molecule.GetAtom(A);
+         int firstAOIndexA = atomA.GetFirstAOIndex();
+         int lastAOIndexA  = atomA.GetLastAOIndex();
+         for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
+            OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
+            for(int B=A; B<totalNumberAtoms; B++){
+               const Atom& atomB = *molecule.GetAtom(B);
+               int firstAOIndexB = atomB.GetFirstAOIndex();
+               int lastAOIndexB  = atomB.GetLastAOIndex();
+               double rAB = molecule.GetDistanceAtoms(atomA, atomB);
+               for(int nu=firstAOIndexB; nu<=lastAOIndexB; nu++){
+                  OrbitalType orbitalNu = atomB.GetValence(nu-firstAOIndexB);
+                  nishimotoMatagaMatrix[A][orbitalMu][B][orbitalNu] = this->GetNishimotoMatagaTwoEleInt(atomA, 
+                                                                                                        orbitalMu, 
+                                                                                                        atomB, 
+                                                                                                        orbitalNu,
+                                                                                                        rAB);
+               }
+            }
+         }
+      }
+      catch(MolDSException ex){
+#pragma omp critical
+         ompErrors << ex.what() << endl ;
+      }
+   }
+   // Exception throwing for omp-region
+   if(!ompErrors.str().empty()){
+      throw MolDSException(ompErrors.str());
+   }
+}
+
 void ZindoS::CalcDiatomicOverlapAOsInDiatomicFrame(double** diatomicOverlapAOs, 
                                                    const Atom& atomA, 
                                                    const Atom& atomB) const{
@@ -839,7 +899,7 @@ void ZindoS::CalcOverlapESsWithAnotherElectronicStructure(double** overlapESs,
    MallocerFreer::GetInstance()->Free<double>(&tmpMatrix,            dimOverlapSingletSDs, dimOverlapESs);
 }
 
-// The order of mol, moJ, moK, moL is consistent with Eq. (9) in [RZ_1973]
+// The order of moI, moJ, moK, moL is consistent with Eq. (9) in [RZ_1973]
 double ZindoS::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL, 
                                            const Molecule& molecule, 
                                            double const* const* fockMatrix, 
@@ -984,24 +1044,146 @@ void ZindoS::DoCIS(){
 }
 
 void ZindoS::CalcCISProperties(){
+//calculate dipole moments and transitiondipolemoment
+{
+   double*** dipoleMOs = NULL;
+   double**  overlapMOs = NULL;
+   int totalNumberAOs = this->molecule->GetTotalNumberAOs();
+   int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
+   int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
+   int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
+   MallocerFreer::GetInstance()->Malloc<double>(&dipoleMOs,       CartesianType_end, totalNumberAOs, totalNumberAOs);
+   MallocerFreer::GetInstance()->Malloc<double>(&overlapMOs, totalNumberAOs, totalNumberAOs);
+   double alpha=1.0;
+   double beta =0.0;
+   //double ompStartTime = omp_get_wtime();
+   MolDS_wrappers::Blas::GetInstance()->Dgemmm(false, false, true, totalNumberAOs, totalNumberAOs, totalNumberAOs, totalNumberAOs,
+                                               alpha, 
+                                               this->fockMatrix,
+                                               this->cartesianMatrix[XAxis],
+                                               this->fockMatrix,
+                                               beta,
+                                               dipoleMOs[XAxis]);
+   MolDS_wrappers::Blas::GetInstance()->Dgemmm(false, false, true, totalNumberAOs, totalNumberAOs, totalNumberAOs, totalNumberAOs,
+                                               alpha, 
+                                               this->fockMatrix,
+                                               this->cartesianMatrix[YAxis],
+                                               this->fockMatrix,
+                                               beta,
+                                               dipoleMOs[YAxis]);
+   MolDS_wrappers::Blas::GetInstance()->Dgemmm(false, false, true, totalNumberAOs, totalNumberAOs, totalNumberAOs, totalNumberAOs,
+                                               alpha, 
+                                               this->fockMatrix,
+                                               this->cartesianMatrix[ZAxis],
+                                               this->fockMatrix,
+                                               beta,
+                                               dipoleMOs[ZAxis]);
+
+   double const* centerOfDipole = this->molecule->GetXyzCOC();
+   // set orign of dipole
+   MolDS_wrappers::Blas::GetInstance()->Dgemmm(false, false, true, totalNumberAOs, totalNumberAOs, totalNumberAOs, totalNumberAOs,
+                                               alpha, 
+                                               this->fockMatrix,
+                                               this->overlapAOs,
+                                               this->fockMatrix,
+                                               beta,
+                                               overlapMOs);
+   MolDS_wrappers::Blas::GetInstance()->Daxpy(totalNumberAOs*totalNumberAOs,
+                                              -centerOfDipole[XAxis], 
+                                              &overlapMOs[0][0],
+                                              &dipoleMOs[XAxis][0][0]);
+   MolDS_wrappers::Blas::GetInstance()->Daxpy(totalNumberAOs*totalNumberAOs,
+                                              -centerOfDipole[YAxis], 
+                                              &overlapMOs[0][0],
+                                              &dipoleMOs[YAxis][0][0]);
+   MolDS_wrappers::Blas::GetInstance()->Daxpy(totalNumberAOs*totalNumberAOs,
+                                              -centerOfDipole[ZAxis], 
+                                              &overlapMOs[0][0],
+                                              &dipoleMOs[ZAxis][0][0]);
+
 
    // dipole moments of excited states
-   this->CalcElectronicDipoleMomentsExcitedState(this->electronicTransitionDipoleMoments,
-                                                 this->fockMatrix,
-                                                 this->matrixCIS,
-                                                 this->cartesianMatrix,
-                                                 *this->molecule, 
-                                                 this->orbitalElectronPopulation,
-                                                 this->overlapAOs);
-   
+   //this->CalcElectronicDipoleMomentsExcitedStates(this->electronicTransitionDipoleMoments,
+   //                                               this->fockMatrix,
+   //                                               this->matrixCIS,
+   //                                               this->cartesianMatrix,
+   //                                               *this->molecule, 
+   //                                               this->orbitalElectronPopulation,
+   //                                               this->overlapAOs);
+   int groundState = 0;
+   for(int k=0; k<Parameters::GetInstance()->GetNumberExcitedStatesCIS(); k++){
+      int excitedState = k+1; // (k+1)-th excited state
+      this->electronicTransitionDipoleMoments[excitedState][excitedState][XAxis] = this->electronicTransitionDipoleMoments[groundState][groundState][XAxis];
+      this->electronicTransitionDipoleMoments[excitedState][excitedState][YAxis] = this->electronicTransitionDipoleMoments[groundState][groundState][YAxis];
+      this->electronicTransitionDipoleMoments[excitedState][excitedState][ZAxis] = this->electronicTransitionDipoleMoments[groundState][groundState][ZAxis];
+      for(int l=0; l<this->matrixCISdimension; l++){
+         // single excitation from I-th (occupied)MO to A-th (virtual)MO
+         int moI = this->GetActiveOccIndex(*this->molecule, l);
+         int moA = this->GetActiveVirIndex(*this->molecule, l);
+         double temp    = matrixCIS[k][l]*matrixCIS[k][l];
+         this->electronicTransitionDipoleMoments[excitedState][excitedState][XAxis] += temp*(-dipoleMOs[XAxis][moI][moI]+dipoleMOs[XAxis][moA][moA]);
+         this->electronicTransitionDipoleMoments[excitedState][excitedState][YAxis] += temp*(-dipoleMOs[YAxis][moI][moI]+dipoleMOs[YAxis][moA][moA]);
+         this->electronicTransitionDipoleMoments[excitedState][excitedState][ZAxis] += temp*(-dipoleMOs[ZAxis][moI][moI]+dipoleMOs[ZAxis][moA][moA]);
+      }
+   }
+
    // transition dipole moment
-   this->CalcElectronicTransitionDipoleMoments(this->electronicTransitionDipoleMoments,
-                                               this->fockMatrix,
-                                               this->matrixCIS,
-                                               this->cartesianMatrix,
-                                               *this->molecule, 
-                                               this->orbitalElectronPopulation,
-                                               this->overlapAOs);
+   //this->CalcElectronicTransitionDipoleMoments(this->electronicTransitionDipoleMoments,
+   //                                            this->fockMatrix,
+   //                                            this->matrixCIS,
+   //                                            this->cartesianMatrix,
+   //                                            *this->molecule, 
+   //                                            this->orbitalElectronPopulation,
+   //                                            this->overlapAOs);
+   for(int k=0; k<Parameters::GetInstance()->GetNumberExcitedStatesCIS(); k++){
+      int excitedState = k+1; // (k+1)-th excited state
+      this->electronicTransitionDipoleMoments[excitedState][groundState][XAxis] = 0.0;
+      this->electronicTransitionDipoleMoments[excitedState][groundState][YAxis] = 0.0;
+      this->electronicTransitionDipoleMoments[excitedState][groundState][ZAxis] = 0.0;
+      for(int l=0; l<this->matrixCISdimension; l++){
+         // single excitation from I-th (occupied)MO to A-th (virtual)MO
+         int moI = this->GetActiveOccIndex(*this->molecule, l);
+         int moA = this->GetActiveVirIndex(*this->molecule, l);
+         //double temp    = matrixCIS[k][l]*matrixCIS[k][l];
+         double temp    = this->matrixCIS[k][l]*sqrt(2.0);
+         this->electronicTransitionDipoleMoments[excitedState][groundState][XAxis] += temp*dipoleMOs[XAxis][moA][moI];
+         this->electronicTransitionDipoleMoments[excitedState][groundState][YAxis] += temp*dipoleMOs[YAxis][moA][moI];
+         this->electronicTransitionDipoleMoments[excitedState][groundState][ZAxis] += temp*dipoleMOs[ZAxis][moA][moI];
+      }
+   }
+   if(Parameters::GetInstance()->RequiresAllTransitionDipoleMomentsCIS()){
+      for(int k=0; k<Parameters::GetInstance()->GetNumberExcitedStatesCIS(); k++){
+         int departureExcitedState = k+1; // (k+1)-th excited state
+         for(int l=k+1; l<Parameters::GetInstance()->GetNumberExcitedStatesCIS(); l++){
+            int destinationExcitedState = l+1; // (l+1)-th excited state
+            this->electronicTransitionDipoleMoments[destinationExcitedState][departureExcitedState][XAxis] = 0.0;
+            this->electronicTransitionDipoleMoments[destinationExcitedState][departureExcitedState][YAxis] = 0.0;
+            this->electronicTransitionDipoleMoments[destinationExcitedState][departureExcitedState][ZAxis] = 0.0;
+            for(int l=0; l<this->matrixCISdimension; l++){
+               // single excitation from I-th (occupied)MO to A-th (virtual)MO
+               int moI = this->GetActiveOccIndex(*this->molecule, l);
+               int moA = this->GetActiveVirIndex(*this->molecule, l);
+               double temp    = matrixCIS[departureExcitedState-1][l]*matrixCIS[destinationExcitedState-1][l];
+               this->electronicTransitionDipoleMoments[destinationExcitedState][departureExcitedState][XAxis] += temp*(-dipoleMOs[XAxis][moI][moI]+dipoleMOs[XAxis][moA][moA]);
+               this->electronicTransitionDipoleMoments[destinationExcitedState][departureExcitedState][YAxis] += temp*(-dipoleMOs[YAxis][moI][moI]+dipoleMOs[YAxis][moA][moA]);
+               this->electronicTransitionDipoleMoments[destinationExcitedState][departureExcitedState][ZAxis] += temp*(-dipoleMOs[ZAxis][moI][moI]+dipoleMOs[ZAxis][moA][moA]);
+            }
+         }
+      }
+
+      for(int k=0; k<Parameters::GetInstance()->GetNumberExcitedStatesCIS()+1; k++){
+         for(int l=k+1; l<Parameters::GetInstance()->GetNumberExcitedStatesCIS()+1; l++){
+            for(int axis=0; axis<CartesianType_end; axis++){
+               this->electronicTransitionDipoleMoments[k][l][axis] 
+                  = this->electronicTransitionDipoleMoments[l][k][axis];
+            }
+         }
+      }
+   }
+   MallocerFreer::GetInstance()->Free<double>(&dipoleMOs, CartesianType_end, totalNumberAOs, totalNumberAOs);
+   MallocerFreer::GetInstance()->Free<double>(&overlapMOs, totalNumberAOs, totalNumberAOs);
+}// end of "calculate dipole moments and transitiondipolemoment"
+
 
 
    // free exciton energies
@@ -1028,30 +1210,27 @@ void ZindoS::CalcCISProperties(){
                                          *this->molecule);
 }
 
-void ZindoS::CalcElectronicDipoleMomentsExcitedState(double*** electronicTransitionDipoleMoments,
-                                                     double const* const* fockMatrix,
-                                                     double const* const* matrixCIS,
-                                                     double const* const* const* cartesianMatrix,
-                                                     const MolDS_base::Molecule& molecule, 
-                                                     double const* const* orbitalElectronPopulation,
-                                                     double const* const* overlapAOs) const{
+void ZindoS::CalcElectronicDipoleMomentsExcitedStates(double*** electronicTransitionDipoleMoments,
+                                                      double const* const* fockMatrix,
+                                                      double const* const* matrixCIS,
+                                                      double const* const* const* cartesianMatrix,
+                                                      const MolDS_base::Molecule& molecule, 
+                                                      double const* const* orbitalElectronPopulation,
+                                                      double const* const* overlapAOs) const{
    int groundState = 0;
    // dipole moment of excited states
    for(int k=0; k<Parameters::GetInstance()->GetNumberExcitedStatesCIS(); k++){
       int excitedState = k+1; // (k+1)-th excited state
-      for(int axis=0; axis<CartesianType_end; axis++){
-         electronicTransitionDipoleMoments[excitedState][excitedState][axis] 
-            = this->GetElectronicTransitionDipoleMoment(excitedState,
-                                                        excitedState,
-                                                        static_cast<CartesianType>(axis),
-                                                        fockMatrix,
-                                                        matrixCIS,
-                                                        cartesianMatrix,
-                                                        molecule,
-                                                        orbitalElectronPopulation,
-                                                        overlapAOs,
-                                                        electronicTransitionDipoleMoments[groundState][groundState]);
-      }
+      this->CalcElectronicTransitionDipoleMoment(electronicTransitionDipoleMoments[excitedState][excitedState],
+                                                 excitedState,
+                                                 excitedState,
+                                                 fockMatrix,
+                                                 matrixCIS,
+                                                 cartesianMatrix,
+                                                 molecule,
+                                                 orbitalElectronPopulation,
+                                                 overlapAOs,
+                                                 electronicTransitionDipoleMoments[groundState][groundState]);
    }
 }
 
@@ -1066,19 +1245,16 @@ void ZindoS::CalcElectronicTransitionDipoleMoments(double*** electronicTransitio
    // transition dipole moments from ground state to excited states
    for(int k=0; k<Parameters::GetInstance()->GetNumberExcitedStatesCIS(); k++){
       int excitedState = k+1; // (k+1)-th excited state
-      for(int axis=0; axis<CartesianType_end; axis++){
-         this->electronicTransitionDipoleMoments[excitedState][groundState][axis] 
-            = this->GetElectronicTransitionDipoleMoment(excitedState,
-                                                        groundState,
-                                                        static_cast<CartesianType>(axis),
-                                                        this->fockMatrix,
-                                                        this->matrixCIS,
-                                                        this->cartesianMatrix,
-                                                        *this->molecule,
-                                                        this->orbitalElectronPopulation,
-                                                        this->overlapAOs,
-                                                        this->electronicTransitionDipoleMoments[groundState][groundState]);
-      }
+      this->CalcElectronicTransitionDipoleMoment(electronicTransitionDipoleMoments[excitedState][groundState],
+                                                 excitedState,
+                                                 groundState,
+                                                 fockMatrix,
+                                                 matrixCIS,
+                                                 cartesianMatrix,
+                                                 molecule,
+                                                 orbitalElectronPopulation,
+                                                 overlapAOs,
+                                                 electronicTransitionDipoleMoments[groundState][groundState]);
    }
 
    if(Parameters::GetInstance()->RequiresAllTransitionDipoleMomentsCIS()){
@@ -1087,70 +1263,79 @@ void ZindoS::CalcElectronicTransitionDipoleMoments(double*** electronicTransitio
          int departureExcitedState = k+1; // (k+1)-th excited state
          for(int l=k+1; l<Parameters::GetInstance()->GetNumberExcitedStatesCIS(); l++){
             int destinationExcitedState = l+1; // (l+1)-th excited state
-            for(int axis=0; axis<CartesianType_end; axis++){
-               this->electronicTransitionDipoleMoments[destinationExcitedState][departureExcitedState][axis] 
-                  = this->GetElectronicTransitionDipoleMoment(destinationExcitedState,
-                                                              departureExcitedState,
-                                                              static_cast<CartesianType>(axis),
-                                                              this->fockMatrix,
-                                                              this->matrixCIS,
-                                                              this->cartesianMatrix,
-                                                              *this->molecule,
-                                                              this->orbitalElectronPopulation,
-                                                              this->overlapAOs,
-                                                              this->electronicTransitionDipoleMoments[groundState][groundState]);
-            }
+            this->CalcElectronicTransitionDipoleMoment(electronicTransitionDipoleMoments[destinationExcitedState][departureExcitedState],
+                                                       destinationExcitedState,
+                                                       departureExcitedState,
+                                                       fockMatrix,
+                                                       matrixCIS,
+                                                       cartesianMatrix,
+                                                       molecule,
+                                                       orbitalElectronPopulation,
+                                                       overlapAOs,
+                                                       electronicTransitionDipoleMoments[groundState][groundState]);
          }
       }
       for(int k=0; k<Parameters::GetInstance()->GetNumberExcitedStatesCIS()+1; k++){
          for(int l=k+1; l<Parameters::GetInstance()->GetNumberExcitedStatesCIS()+1; l++){
             for(int axis=0; axis<CartesianType_end; axis++){
-               this->electronicTransitionDipoleMoments[k][l][axis] 
-                  = this->electronicTransitionDipoleMoments[l][k][axis];
+               electronicTransitionDipoleMoments[k][l][axis] 
+                  = electronicTransitionDipoleMoments[l][k][axis];
             }
          }
       }
    }
 }
 
-double ZindoS::GetElectronicTransitionDipoleMoment(int to, int from, CartesianType axis,
-                                                   double const* const* fockMatrix,
-                                                   double const* const* matrixCIS,
-                                                   double const* const* const* cartesianMatrix,
-                                                   const MolDS_base::Molecule& molecule, 
-                                                   double const* const* orbitalElectronPopulation,
-                                                   double const* const* overlapAOs,
-                                                   double const* groundStateDipole) const{
-   double value = 0.0;
+void ZindoS::CalcElectronicTransitionDipoleMoment(double* transitionDipoleMoment,
+                                                  int to, int from,
+                                                  double const* const* fockMatrix,
+                                                  double const* const* matrixCIS,
+                                                  double const* const* const* cartesianMatrix,
+                                                  const MolDS_base::Molecule& molecule, 
+                                                  double const* const* orbitalElectronPopulation,
+                                                  double const* const* overlapAOs,
+                                                  double const* groundStateDipole) const{
+   double valueX = 0.0;
+   double valueY = 0.0;
+   double valueZ = 0.0;
+   double const* xyzCOC = molecule.GetXyzCOC();
    int groundState = 0;
+   int totalNumberAOs = molecule.GetTotalNumberAOs();
    stringstream ompErrors;
    if(Parameters::GetInstance()->GetNumberExcitedStatesCIS() < from ||
       Parameters::GetInstance()->GetNumberExcitedStatesCIS() < to ){
       stringstream ss;
-      ss << this->errorMessageGetElectronicTransitionDipoleMomentBadState;
+      ss << this->errorMessageCalcElectronicTransitionDipoleMomentBadState;
       ss << this->errorMessageFromState << from << endl;
       ss << this->errorMessageToState << to << endl;
-      ss << this->errorMessageCartesianType << CartesianTypeStr(axis) << endl;
       throw MolDSException(ss.str());
    }
 
    if(from != to){
       if(from != groundState && to != groundState){
          // transition dipole moment between different excited states
-#pragma omp parallel for reduction(+:value) schedule(auto)
+#pragma omp parallel for reduction(+:valueX,valueY,valueZ) schedule(auto)
          for(int l=0; l<this->matrixCISdimension; l++){
             try{
-               double temp = 0.0;
+               double temp  = 0.0;
+               double tempX = 0.0;
+               double tempY = 0.0;
+               double tempZ = 0.0;
                // single excitation from I-th (occupied)MO to A-th (virtual)MO
                int moI = this->GetActiveOccIndex(molecule, l);
                int moA = this->GetActiveVirIndex(molecule, l);
-               for(int mu=0; mu<molecule.GetTotalNumberAOs(); mu++){
-                  for(int nu=0; nu<molecule.GetTotalNumberAOs(); nu++){
-                     temp += (-1.0*fockMatrix[moI][mu]*fockMatrix[moI][nu] + fockMatrix[moA][mu]*fockMatrix[moA][nu])
-                             *(cartesianMatrix[mu][nu][axis] - molecule.GetXyzCOC()[axis]*overlapAOs[mu][nu]);
+               for(int mu=0; mu<totalNumberAOs; mu++){
+                  for(int nu=0; nu<totalNumberAOs; nu++){
+                     temp   = (-1.0*fockMatrix[moI][mu]*fockMatrix[moI][nu] + fockMatrix[moA][mu]*fockMatrix[moA][nu]);
+                     tempX += temp*(cartesianMatrix[XAxis][mu][nu] - xyzCOC[XAxis]*overlapAOs[mu][nu]);
+                     tempY += temp*(cartesianMatrix[YAxis][mu][nu] - xyzCOC[YAxis]*overlapAOs[mu][nu]);
+                     tempZ += temp*(cartesianMatrix[ZAxis][mu][nu] - xyzCOC[ZAxis]*overlapAOs[mu][nu]);
                   }
                }
-               value += matrixCIS[from-1][l]*matrixCIS[to-1][l]*temp;
+               temp    = matrixCIS[from-1][l]*matrixCIS[to-1][l];
+               valueX += temp*tempX;
+               valueY += temp*tempY;
+               valueZ += temp*tempZ;
             }
             catch(MolDSException ex){
 #pragma omp critical
@@ -1161,25 +1346,34 @@ double ZindoS::GetElectronicTransitionDipoleMoment(int to, int from, CartesianTy
          if(!ompErrors.str().empty()){
             throw MolDSException(ompErrors.str());
          }
-         return value;
+         transitionDipoleMoment[XAxis] = valueX;
+         transitionDipoleMoment[YAxis] = valueY;
+         transitionDipoleMoment[ZAxis] = valueZ;
       }
       else if(from == groundState && to != groundState){
          // transition dipole moment from the ground to excited states
-#pragma omp parallel for reduction(+:value) schedule(auto)
+#pragma omp parallel for reduction(+:valueX,valueY,valueZ) schedule(auto)
          for(int l=0; l<this->matrixCISdimension; l++){
             try{
-               double temp = 0.0;
+               double temp  = 0.0;
+               double tempX = 0.0;
+               double tempY = 0.0;
+               double tempZ = 0.0;
                // single excitation from I-th (occupied)MO to A-th (virtual)MO
                int moI = this->GetActiveOccIndex(molecule, l);
                int moA = this->GetActiveVirIndex(molecule, l);
-               for(int mu=0; mu<molecule.GetTotalNumberAOs(); mu++){
-                  for(int nu=0; nu<molecule.GetTotalNumberAOs(); nu++){
-                     temp += fockMatrix[moA][mu]
-                             *(cartesianMatrix[mu][nu][axis] - molecule.GetXyzCOC()[axis]*overlapAOs[mu][nu])
-                             *fockMatrix[moI][nu];
+               for(int mu=0; mu<totalNumberAOs; mu++){
+                  for(int nu=0; nu<totalNumberAOs; nu++){
+                     temp   = fockMatrix[moA][mu]*fockMatrix[moI][nu];
+                     tempX += temp*(cartesianMatrix[XAxis][mu][nu] - xyzCOC[XAxis]*overlapAOs[mu][nu]);
+                     tempY += temp*(cartesianMatrix[YAxis][mu][nu] - xyzCOC[YAxis]*overlapAOs[mu][nu]);
+                     tempZ += temp*(cartesianMatrix[ZAxis][mu][nu] - xyzCOC[ZAxis]*overlapAOs[mu][nu]);
                   }
                }
-               value += this->matrixCIS[to-1][l]*sqrt(2.0)*temp;
+               temp    = this->matrixCIS[to-1][l]*sqrt(2.0);
+               valueX += temp*tempX;
+               valueY += temp*tempY;
+               valueZ += temp*tempZ;
             }
             catch(MolDSException ex){
 #pragma omp critical
@@ -1190,25 +1384,34 @@ double ZindoS::GetElectronicTransitionDipoleMoment(int to, int from, CartesianTy
          if(!ompErrors.str().empty()){
             throw MolDSException(ompErrors.str());
          }
-         return value;
+         transitionDipoleMoment[XAxis] = valueX;
+         transitionDipoleMoment[YAxis] = valueY;
+         transitionDipoleMoment[ZAxis] = valueZ;
       }
       else if(from != groundState && to == groundState){
          // transition dipole moment from the excited to ground states
-#pragma omp parallel for reduction(+:value) schedule(auto)
+#pragma omp parallel for reduction(+:valueX,valueY,valueZ) schedule(auto)
          for(int l=0; l<this->matrixCISdimension; l++){
             try{
-               double temp = 0.0;
+               double temp  = 0.0;
+               double tempX = 0.0;
+               double tempY = 0.0;
+               double tempZ = 0.0;
                // single excitation from I-th (occupied)MO to A-th (virtual)MO
                int moI = this->GetActiveOccIndex(molecule, l);
                int moA = this->GetActiveVirIndex(molecule, l);
-               for(int mu=0; mu<molecule.GetTotalNumberAOs(); mu++){
-                  for(int nu=0; nu<molecule.GetTotalNumberAOs(); nu++){
-                     temp += fockMatrix[moI][mu]
-                             *(cartesianMatrix[mu][nu][axis] - molecule.GetXyzCOC()[axis]*overlapAOs[mu][nu])
-                             *fockMatrix[moA][nu];
+               for(int mu=0; mu<totalNumberAOs; mu++){
+                  for(int nu=0; nu<totalNumberAOs; nu++){
+                     temp   = fockMatrix[moI][mu]*fockMatrix[moA][nu];
+                     tempX += temp*(cartesianMatrix[XAxis][mu][nu] - xyzCOC[XAxis]*overlapAOs[mu][nu]);
+                     tempY += temp*(cartesianMatrix[YAxis][mu][nu] - xyzCOC[YAxis]*overlapAOs[mu][nu]);
+                     tempZ += temp*(cartesianMatrix[ZAxis][mu][nu] - xyzCOC[ZAxis]*overlapAOs[mu][nu]);
                   }
                }
-               value += matrixCIS[from-1][l]*sqrt(2.0)*temp;
+               temp    = matrixCIS[from-1][l]*sqrt(2.0);
+               valueX += temp*tempX;
+               valueY += temp*tempY;
+               valueZ += temp*tempZ;
             }
             catch(MolDSException ex){
 #pragma omp critical
@@ -1219,26 +1422,36 @@ double ZindoS::GetElectronicTransitionDipoleMoment(int to, int from, CartesianTy
          if(!ompErrors.str().empty()){
             throw MolDSException(ompErrors.str());
          }
-         return value;
+         transitionDipoleMoment[XAxis] = valueX;
+         transitionDipoleMoment[YAxis] = valueY;
+         transitionDipoleMoment[ZAxis] = valueZ;
       }
    }
    else{
       if(from != groundState){
          // dipole moment of the excited state. It is needed that the dipole of ground state has been already calculated!!
-#pragma omp parallel for reduction(+:value) schedule(auto)
+#pragma omp parallel for reduction(+:valueX,valueY,valueZ) schedule(auto)
          for(int l=0; l<this->matrixCISdimension; l++){
             try{
-               double temp = 0.0;
+               double temp  = 0.0;
+               double tempX = 0.0;
+               double tempY = 0.0;
+               double tempZ = 0.0;
                // single excitation from I-th (occupied)MO to A-th (virtual)MO
                int moI = this->GetActiveOccIndex(molecule, l);
                int moA = this->GetActiveVirIndex(molecule, l);
-               for(int mu=0; mu<molecule.GetTotalNumberAOs(); mu++){
-                  for(int nu=0; nu<molecule.GetTotalNumberAOs(); nu++){
-                     temp += (-1.0*fockMatrix[moI][mu]*fockMatrix[moI][nu] + fockMatrix[moA][mu]*fockMatrix[moA][nu])
-                             *(cartesianMatrix[mu][nu][axis] - molecule.GetXyzCOC()[axis]*overlapAOs[mu][nu]);
+               for(int mu=0; mu<totalNumberAOs; mu++){
+                  for(int nu=0; nu<totalNumberAOs; nu++){
+                     temp   = (-1.0*fockMatrix[moI][mu]*fockMatrix[moI][nu] + fockMatrix[moA][mu]*fockMatrix[moA][nu]);
+                     tempX += temp*(cartesianMatrix[XAxis][mu][nu] - xyzCOC[XAxis]*overlapAOs[mu][nu]);
+                     tempY += temp*(cartesianMatrix[YAxis][mu][nu] - xyzCOC[YAxis]*overlapAOs[mu][nu]);
+                     tempZ += temp*(cartesianMatrix[ZAxis][mu][nu] - xyzCOC[ZAxis]*overlapAOs[mu][nu]);
                   }
                }
-               value += matrixCIS[from-1][l]*matrixCIS[to-1][l]*temp;
+               temp    = matrixCIS[from-1][l]*matrixCIS[to-1][l];
+               valueX += temp*tempX;
+               valueY += temp*tempY;
+               valueZ += temp*tempZ;
             }
             catch(MolDSException ex){
 #pragma omp critical
@@ -1249,22 +1462,22 @@ double ZindoS::GetElectronicTransitionDipoleMoment(int to, int from, CartesianTy
          if(!ompErrors.str().empty()){
             throw MolDSException(ompErrors.str());
          }
-         value += groundStateDipole[axis];
-         return value;
+         transitionDipoleMoment[XAxis] = valueX + groundStateDipole[XAxis];
+         transitionDipoleMoment[YAxis] = valueY + groundStateDipole[YAxis];
+         transitionDipoleMoment[ZAxis] = valueZ + groundStateDipole[ZAxis];
       }
       else{
          // dipole moment of the ground state
-         value = Cndo2::GetElectronicTransitionDipoleMoment(to,
-                                                            from,
-                                                            axis,
-                                                            fockMatrix,
-                                                            matrixCIS,
-                                                            cartesianMatrix,
-                                                            molecule,
-                                                            orbitalElectronPopulation,
-                                                            overlapAOs,
-                                                            NULL);
-         return value;
+         Cndo2::CalcElectronicTransitionDipoleMoment(transitionDipoleMoment,
+                                                     to,
+                                                     from,
+                                                     fockMatrix,
+                                                     matrixCIS,
+                                                     cartesianMatrix,
+                                                     molecule,
+                                                     orbitalElectronPopulation,
+                                                     overlapAOs,
+                                                     NULL);
       }
    }
 }
@@ -2027,273 +2240,267 @@ void ZindoS::CalcCISMatrix(double** matrixCIS) const{
    this->OutputLog(this->messageStartCalcCISMatrix);
    double ompStartTime = omp_get_wtime();
 
-   stringstream ompErrors;
-#pragma omp parallel for schedule(auto)
-   for(int k=0; k<this->matrixCISdimension; k++){
-      try{
-         // 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);
+   int totalNumberAtoms = this->molecule->GetNumberAtoms();
+   double**** nishimotoMatagaMatrix=NULL;
+   try{
+      MallocerFreer::GetInstance()->Malloc<double>(&nishimotoMatagaMatrix, totalNumberAtoms, OrbitalType_end, totalNumberAtoms, OrbitalType_end);
+      this->CalcNishimotoMatagaMatrix(nishimotoMatagaMatrix, *this->molecule);
 
-         for(int l=k; l<this->matrixCISdimension; l++){
-            // single excitation from J-th (occupied)MO to B-th (virtual)MO
-            int moJ = this->GetActiveOccIndex(*this->molecule, l);
-            int moB = this->GetActiveVirIndex(*this->molecule, l);
-            double value=0.0;
-         
-            // Fast algorithm, but this is not easy to read. 
-            // Slow algorithm is also written below.
-            double gamma;
-            double exchange;
-            double coulomb;
-            // Off diagonal term (right upper)
-            if(k<l){
-               for(int A=0; A<molecule->GetNumberAtoms(); A++){
-                  const Atom& atomA = *molecule->GetAtom(A);
-                  int firstAOIndexA = atomA.GetFirstAOIndex();
-                  int lastAOIndexA  = atomA.GetLastAOIndex();
+      stringstream ompErrors;
+#pragma omp parallel for schedule(auto)
+      for(int k=0; k<this->matrixCISdimension; k++){
+         try{
+            // 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);
+
+            for(int l=k; l<this->matrixCISdimension; l++){
+               // single excitation from J-th (occupied)MO to B-th (virtual)MO
+               int moJ = this->GetActiveOccIndex(*this->molecule, l);
+               int moB = this->GetActiveVirIndex(*this->molecule, l);
+
+               // Fast algorithm, but this is not easy to read. Slow algorithm is also written below.
+               if(k<l){
+                  // Off diagonal term (right upper)
+                  matrixCIS[k][l] = this->GetCISOffDiagElement(nishimotoMatagaMatrix, *this->molecule, this->fockMatrix, moI, moA, moJ, moB);
+               }
+               else if(k==l){
+                  // Diagonal term
+                  matrixCIS[k][l] = this->GetCISDiagElement(energiesMO, nishimotoMatagaMatrix, *this->molecule, this->fockMatrix, moI, moA);
+               } 
+               // End of the fast algorith.
+
+               /*// Slow algorith, but this is easy to read. Fast altorithm is also written above.
+               double value=0.0;
+               value = 2.0*this->GetMolecularIntegralElement(moA, moI, moJ, moB, 
+                                                             *this->molecule, 
+                                                             this->fockMatrix, 
+                                                             NULL)
+                          -this->GetMolecularIntegralElement(moA, moB, moI, moJ, 
+                                                             *this->molecule, 
+                                                             this->fockMatrix, 
+                                                             NULL);
+               if(k==l){
+                  value += this->energiesMO[moA] - this->energiesMO[moI];
+               }
+               matrixCIS[k][l] = value;
+               // End of the slow algorith. */
+            }
+         }
+         catch(MolDSException ex){
+#pragma omp critical
+            ompErrors << ex.what() << endl ;
+         }
+      } // end of k-loop
+      // Exception throwing for omp-region
+      if(!ompErrors.str().empty()){
+         throw MolDSException(ompErrors.str());
+      }
+   }
+   catch(MolDSException ex){
+      MallocerFreer::GetInstance()->Free<double>(&nishimotoMatagaMatrix, totalNumberAtoms, OrbitalType_end, totalNumberAtoms, OrbitalType_end);
+      throw ex;
+   }
+   MallocerFreer::GetInstance()->Free<double>(&nishimotoMatagaMatrix, totalNumberAtoms, OrbitalType_end, totalNumberAtoms, OrbitalType_end);
+   double ompEndTime = omp_get_wtime();
+   this->OutputLog(boost::format("%s%lf%s\n%s") % this->messageOmpElapsedTimeCalcCISMarix.c_str()
+                                                % (ompEndTime - ompStartTime)
+                                                % this->messageUnitSec.c_str()
+                                                % this->messageDoneCalcCISMatrix.c_str() );
+}
 
-                  for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
-                     OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
+double ZindoS::GetCISDiagElement(double const* energiesMO,
+                                 double const* const* const* const* nishimotoMatagaMatrix,
+                                 const MolDS_base::Molecule& molecule,
+                                 double const* const* fockMatrix, 
+                                 int moI,
+                                 int moA) const{
+   double value    = energiesMO[moA] - energiesMO[moI];
+   double gamma    = 0.0;
+   double exchange = 0.0;
+   double coulomb  = 0.0;
+   int totalNumberAtoms = molecule.GetNumberAtoms();
+   for(int A=0; A<totalNumberAtoms; A++){
+      const Atom& atomA = *molecule.GetAtom(A);
+      int firstAOIndexA = atomA.GetFirstAOIndex();
+      int lastAOIndexA  = atomA.GetLastAOIndex();
 
-                     // CNDO term
-                     for(int B=A; B<molecule->GetNumberAtoms(); B++){
-                        const Atom& atomB = *molecule->GetAtom(B);
-                        int firstAOIndexB = atomB.GetFirstAOIndex();
-                        int lastAOIndexB  = atomB.GetLastAOIndex();
-
-                        for(int nu=firstAOIndexB; nu<=lastAOIndexB; nu++){
-                           OrbitalType orbitalNu = atomB.GetValence(nu-firstAOIndexB);
-
-                           if(A<B){
-                              gamma = this->GetNishimotoMatagaTwoEleInt(atomA, 
-                                                                        orbitalMu, 
-                                                                        atomB, 
-                                                                        orbitalNu);
-                              value += 2.0*gamma*fockMatrix[moA][mu]
-                                                *fockMatrix[moI][mu]
-                                                *fockMatrix[moJ][nu]
-                                                *fockMatrix[moB][nu];
-                              value -=     gamma*fockMatrix[moA][mu]
-                                                *fockMatrix[moB][mu]
-                                                *fockMatrix[moI][nu]
-                                                *fockMatrix[moJ][nu];
-                              value += 2.0*gamma*fockMatrix[moA][nu]
-                                                *fockMatrix[moI][nu]
-                                                *fockMatrix[moJ][mu]
-                                                *fockMatrix[moB][mu];
-                              value -=     gamma*fockMatrix[moA][nu]
-                                                *fockMatrix[moB][nu]
-                                                *fockMatrix[moI][mu]
-                                                *fockMatrix[moJ][mu];
-                           }
-                           else{
-                              gamma = atomA.GetZindoF0ss();
-                              value += 2.0*gamma*fockMatrix[moA][mu]
-                                                *fockMatrix[moI][mu]
-                                                *fockMatrix[moJ][nu]
-                                                *fockMatrix[moB][nu];
-                              value -=     gamma*fockMatrix[moA][mu]
-                                                *fockMatrix[moB][mu]
-                                                *fockMatrix[moI][nu]
-                                                *fockMatrix[moJ][nu];
-                           }  
-                        }
-                     }
+      for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
+         OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
+         double tmp1 = fockMatrix[moA][mu]*fockMatrix[moI][mu]; 
+         double tmp2 = fockMatrix[moI][mu]*fockMatrix[moI][mu]; 
+         double tmp3 = fockMatrix[moA][mu]*fockMatrix[moA][mu];
 
-                     // Aditional term for INDO or ZIND/S, see Eq. (10) in [RZ_1973]
-                     for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
-                        OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
-
-                        if(mu!=nu){
-                           exchange = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
-                           value += 2.0*exchange*fockMatrix[moA][mu]
-                                                *fockMatrix[moI][nu]
-                                                *fockMatrix[moJ][nu]
-                                                *fockMatrix[moB][mu];
-                           value += 2.0*exchange*fockMatrix[moA][mu]
-                                                *fockMatrix[moI][nu]
-                                                *fockMatrix[moJ][mu]
-                                                *fockMatrix[moB][nu];
-                           value -=     exchange*fockMatrix[moA][mu]
-                                                *fockMatrix[moB][nu]
-                                                *fockMatrix[moI][nu]
-                                                *fockMatrix[moJ][mu];
-                           value -=     exchange*fockMatrix[moA][mu]
-                                                *fockMatrix[moB][nu]
-                                                *fockMatrix[moI][mu]
-                                                *fockMatrix[moJ][nu];
-                        }
+         // CNDO term
+         for(int B=A; B<totalNumberAtoms; B++){
+            const Atom& atomB = *molecule.GetAtom(B);
+            int firstAOIndexB = atomB.GetFirstAOIndex();
+            int lastAOIndexB  = atomB.GetLastAOIndex();
 
-                        coulomb = this->GetCoulombInt(orbitalMu, orbitalNu, atomA);
+            for(int nu=firstAOIndexB; nu<=lastAOIndexB; nu++){
+               OrbitalType orbitalNu = atomB.GetValence(nu-firstAOIndexB);
 
-                        if( (orbitalMu == s || orbitalMu == px || orbitalMu == py || orbitalMu == pz) &&
-                            (orbitalNu == s || orbitalNu == px || orbitalNu == py || orbitalNu == pz) ){
-                              gamma = atomA.GetZindoF0ss();
-                        }
-                        else{
-                           stringstream ss;
-                           ss << this->errorMessageCalcCISMatrix;
-                           ss << this->errorMessageAtomType << AtomTypeStr(atomA.GetAtomType()) << "\n";
-                           ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalMu) << "\n";
-                           ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalNu) << "\n";
-#pragma omp critical
-                           ompErrors << ss.str() << endl ;
-                        }   
-
-                        value += 2.0*(coulomb-gamma)*fockMatrix[moA][mu]
-                                                    *fockMatrix[moI][mu]
-                                                    *fockMatrix[moJ][nu]
-                                                    *fockMatrix[moB][nu];
-                        value -=     (coulomb-gamma)*fockMatrix[moA][mu]
-                                                    *fockMatrix[moB][mu]
-                                                    *fockMatrix[moI][nu]
-                                                    *fockMatrix[moJ][nu];
-                     }
-                  }
+               if(A<B){
+                  gamma = nishimotoMatagaMatrix[A][orbitalMu][B][orbitalNu];
+                  value += 4.0*gamma*tmp1
+                                    *fockMatrix[moA][nu]
+                                    *fockMatrix[moI][nu];
+                  value -=     gamma*tmp2
+                                    *fockMatrix[moA][nu]
+                                    *fockMatrix[moA][nu];
+                  value -=     gamma*tmp3
+                                    *fockMatrix[moI][nu]
+                                    *fockMatrix[moI][nu];
                }
+               else{
+                  gamma = atomA.GetZindoF0ss();
+                  value += 2.0*gamma*tmp1
+                                    *fockMatrix[moA][nu]
+                                    *fockMatrix[moI][nu];
+                  value -=     gamma*tmp2
+                                    *fockMatrix[moA][nu]
+                                    *fockMatrix[moA][nu];
+               }     
             }
-            // Diagonal term
-            else if(k==l){
-               value = this->energiesMO[moA] - this->energiesMO[moI];
-               for(int A=0; A<molecule->GetNumberAtoms(); A++){
-                  const Atom& atomA = *molecule->GetAtom(A);
-                  int firstAOIndexA = atomA.GetFirstAOIndex();
-                  int lastAOIndexA  = atomA.GetLastAOIndex();
+         }
 
-                  for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
-                     OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
+         //Aditional term for INDO or ZIND/S, see Eq. (10) in [RZ_1973]
+         for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
+            OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
 
-                     // CNDO term
-                     for(int B=A; B<molecule->GetNumberAtoms(); B++){
-                        const Atom& atomB = *molecule->GetAtom(B);
-                        int firstAOIndexB = atomB.GetFirstAOIndex();
-                        int lastAOIndexB  = atomB.GetLastAOIndex();
-
-                        for(int nu=firstAOIndexB; nu<=lastAOIndexB; nu++){
-                           OrbitalType orbitalNu = atomB.GetValence(nu-firstAOIndexB);
-
-                           if(A<B){
-                              gamma = this->GetNishimotoMatagaTwoEleInt(atomA, 
-                                                                        orbitalMu, 
-                                                                        atomB, 
-                                                                        orbitalNu);
-                              value += 2.0*gamma*fockMatrix[moI][mu]
-                                                *fockMatrix[moA][mu]
-                                                *fockMatrix[moA][nu]
-                                                *fockMatrix[moI][nu];
-                              value -=     gamma*fockMatrix[moI][mu]
-                                                *fockMatrix[moI][mu]
-                                                *fockMatrix[moA][nu]
-                                                *fockMatrix[moA][nu];
-                              value += 2.0*gamma*fockMatrix[moI][nu]
-                                                *fockMatrix[moA][nu]
-                                                *fockMatrix[moA][mu]
-                                                *fockMatrix[moI][mu];
-                              value -=     gamma*fockMatrix[moI][nu]
-                                                *fockMatrix[moI][nu]
-                                                *fockMatrix[moA][mu]
-                                                *fockMatrix[moA][mu];
-                           }
-                           else{
-                              gamma = atomA.GetZindoF0ss();
-                              value += 2.0*gamma*fockMatrix[moI][mu]
-                                                *fockMatrix[moA][mu]
-                                                *fockMatrix[moA][nu]
-                                                *fockMatrix[moI][nu];
-                              value -=     gamma*fockMatrix[moI][mu]
-                                                *fockMatrix[moI][mu]
-                                                *fockMatrix[moA][nu]
-                                                *fockMatrix[moA][nu];
-                           }     
-                        }
-                     }
+            if(mu!=nu){
+               exchange = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
+               value += 2.0*exchange*tmp2
+                                    *fockMatrix[moA][nu]
+                                    *fockMatrix[moA][nu];
+            }
 
-                     // Aditional term for INDO or ZIND/S, see Eq. (10) in [RZ_1973]
-                     for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
-                        OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
-
-                        if(mu!=nu){
-                           exchange = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
-                           value += 2.0*exchange*fockMatrix[moI][mu]
-                                                *fockMatrix[moA][nu]
-                                                *fockMatrix[moA][nu]
-                                                *fockMatrix[moI][mu];
-                           value += 2.0*exchange*fockMatrix[moI][mu]
-                                                *fockMatrix[moA][nu]
-                                                *fockMatrix[moA][mu]
-                                                *fockMatrix[moI][nu];
-                           value -=     exchange*fockMatrix[moI][mu]
-                                                *fockMatrix[moI][nu]
-                                                *fockMatrix[moA][nu]
-                                                *fockMatrix[moA][mu];
-                           value -=     exchange*fockMatrix[moI][mu]
-                                                *fockMatrix[moI][nu]
-                                                *fockMatrix[moA][mu]
-                                                *fockMatrix[moA][nu];
-                        }
+            coulomb = this->GetCoulombInt(orbitalMu, orbitalNu, atomA);
+
+            if( (orbitalMu == s || orbitalMu == px || orbitalMu == py || orbitalMu == pz) &&
+                (orbitalNu == s || orbitalNu == px || orbitalNu == py || orbitalNu == pz) ){
+                  gamma = atomA.GetZindoF0ss();
+            }
+            else{
+               stringstream ss;
+               ss << this->errorMessageCalcCISMatrix;
+               ss << this->errorMessageAtomType << AtomTypeStr(atomA.GetAtomType()) << "\n";
+               ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalMu) << "\n";
+               ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalNu) << "\n";
+               throw MolDSException(ss.str());
+            }   
+            value += 2.0*(coulomb-gamma)*tmp1
+                                        *fockMatrix[moA][nu]
+                                        *fockMatrix[moI][nu];
+            value -=     (coulomb-gamma)*tmp2
+                                        *fockMatrix[moA][nu]
+                                        *fockMatrix[moA][nu];
+         }
+      }
+   }
+   return value;
+}
 
-                        coulomb = this->GetCoulombInt(orbitalMu, orbitalNu, atomA);
+double ZindoS::GetCISOffDiagElement(double const* const* const* const* nishimotoMatagaMatrix,
+                                    const MolDS_base::Molecule& molecule,
+                                    double const* const* fockMatrix, 
+                                    int moI,
+                                    int moA,
+                                    int moJ,
+                                    int moB) const{
+   double value    = 0.0;
+   double gamma    = 0.0;
+   double exchange = 0.0;
+   double coulomb  = 0.0;
+   int totalNumberAtoms = molecule.GetNumberAtoms();
+   for(int A=0; A<totalNumberAtoms; A++){
+      const Atom& atomA = *molecule.GetAtom(A);
+      int firstAOIndexA = atomA.GetFirstAOIndex();
+      int lastAOIndexA  = atomA.GetLastAOIndex();
 
-                        if( (orbitalMu == s || orbitalMu == px || orbitalMu == py || orbitalMu == pz) &&
-                            (orbitalNu == s || orbitalNu == px || orbitalNu == py || orbitalNu == pz) ){
-                              gamma = atomA.GetZindoF0ss();
-                        }
-                        else{
-                           stringstream ss;
-                           ss << this->errorMessageCalcCISMatrix;
-                           ss << this->errorMessageAtomType << AtomTypeStr(atomA.GetAtomType()) << "\n";
-                           ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalMu) << "\n";
-                           ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalNu) << "\n";
-#pragma omp critical
-                           ompErrors << ss.str() << endl ;
-                        }   
-
-                        value += 2.0*(coulomb-gamma)*fockMatrix[moI][mu]
-                                                    *fockMatrix[moA][mu]
-                                                    *fockMatrix[moA][nu]
-                                                    *fockMatrix[moI][nu];
-                        value -=     (coulomb-gamma)*fockMatrix[moI][mu]
-                                                    *fockMatrix[moI][mu]
-                                                    *fockMatrix[moA][nu]
-                                                    *fockMatrix[moA][nu];
-                     }
-                  }
+      for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){ 
+         OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA); 
+         double tmp1 = fockMatrix[moA][mu]*fockMatrix[moI][mu]; 
+         double tmp2 = fockMatrix[moA][mu]*fockMatrix[moB][mu]; 
+         double tmp3 = fockMatrix[moJ][mu]*fockMatrix[moB][mu];
+         double tmp4 = fockMatrix[moI][mu]*fockMatrix[moJ][mu];
+
+         // CNDO term
+         for(int B=A; B<totalNumberAtoms; B++){
+            const Atom& atomB = *molecule.GetAtom(B);
+            int firstAOIndexB = atomB.GetFirstAOIndex(); 
+            int lastAOIndexB  = atomB.GetLastAOIndex(); 
+            for(int nu=firstAOIndexB; nu<=lastAOIndexB; nu++){
+               OrbitalType orbitalNu = atomB.GetValence(nu-firstAOIndexB);
+               
+               if(A<B){
+                  gamma = nishimotoMatagaMatrix[A][orbitalMu][B][orbitalNu];
+                  value += 2.0*gamma*tmp1
+                                    *fockMatrix[moJ][nu]
+                                    *fockMatrix[moB][nu];
+                  value -=     gamma*tmp2
+                                    *fockMatrix[moI][nu]
+                                    *fockMatrix[moJ][nu];
+                  value += 2.0*gamma*tmp3
+                                    *fockMatrix[moA][nu]
+                                    *fockMatrix[moI][nu];
+                  value -=     gamma*tmp4
+                                    *fockMatrix[moA][nu]
+                                    *fockMatrix[moB][nu];
                }
+               else{
+                  gamma = atomA.GetZindoF0ss();
+                  value += 2.0*gamma*tmp1
+                                    *fockMatrix[moJ][nu]
+                                    *fockMatrix[moB][nu];
+                  value -=     gamma*tmp2
+                                    *fockMatrix[moI][nu]
+                                    *fockMatrix[moJ][nu];
+               }  
             }
-            // End of the fast algorith.
-            /* 
-            // Slow algorith, but this is easy to read. Fast altorithm is also written above.
-            value = 2.0*this->GetMolecularIntegralElement(moA, moI, moJ, moB, 
-                                                          this->molecule, 
-                                                          this->fockMatrix, 
-                                                          NULL)
-                       -this->GetMolecularIntegralElement(moA, moB, moI, moJ, 
-                                                          this->molecule, 
-                                                          this->fockMatrix, 
-                                                          NULL);
-            if(k==l){
-               value += this->energiesMO[moA] - this->energiesMO[moI];
+         }
+
+         // Aditional term for INDO or ZIND/S, see Eq. (10) in [RZ_1973]
+         double tmp5 = fockMatrix[moA][mu]*fockMatrix[moB][mu];
+         double tmp6 = fockMatrix[moA][mu]*fockMatrix[moJ][mu];
+         double tmp7 = fockMatrix[moA][mu]*fockMatrix[moI][mu];
+         for(int nu=firstAOIndexA; nu<=lastAOIndexA; nu++){
+            OrbitalType orbitalNu = atomA.GetValence(nu-firstAOIndexA);
+
+            if(mu!=nu){
+               exchange = this->GetExchangeInt(orbitalMu, orbitalNu, atomA);
+               value += 2.0*exchange
+                       *(tmp5*fockMatrix[moJ][nu] + tmp6*fockMatrix[moB][nu])
+                       *fockMatrix[moI][nu];
+               value -= exchange
+                       *(tmp6*fockMatrix[moI][nu] + tmp7*fockMatrix[moJ][nu])
+                       *fockMatrix[moB][nu];
+            }
+
+            coulomb = this->GetCoulombInt(orbitalMu, orbitalNu, atomA);
+
+            if( (orbitalMu == s || orbitalMu == px || orbitalMu == py || orbitalMu == pz) &&
+                (orbitalNu == s || orbitalNu == px || orbitalNu == py || orbitalNu == pz) ){
+                  gamma = atomA.GetZindoF0ss();
             }
-            // End of the slow algorith.
-            */
-            matrixCIS[k][l] = value;
+            else{
+               stringstream ss;
+               ss << this->errorMessageCalcCISMatrix;
+               ss << this->errorMessageAtomType << AtomTypeStr(atomA.GetAtomType()) << "\n";
+               ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalMu) << "\n";
+               ss << this->errorMessageOrbitalType << OrbitalTypeStr(orbitalNu) << "\n";
+               throw MolDSException(ss.str());
+            }   
+            value += 2.0*(coulomb-gamma)*tmp7
+                                        *fockMatrix[moJ][nu]
+                                        *fockMatrix[moB][nu];
+            value -=     (coulomb-gamma)*tmp5
+                                        *fockMatrix[moI][nu]
+                                        *fockMatrix[moJ][nu];
          }
       }
-      catch(MolDSException ex){
-#pragma omp critical
-         ompErrors << ex.what() << endl ;
-      }
-   }
-   // Exception throwing for omp-region
-   if(!ompErrors.str().empty()){
-      throw MolDSException(ompErrors.str());
    }
-   double ompEndTime = omp_get_wtime();
-   this->OutputLog(boost::format("%s%lf%s\n%s") % this->messageOmpElapsedTimeCalcCISMarix.c_str()
-                                                % (ompEndTime - ompStartTime)
-                                                % this->messageUnitSec.c_str()
-                                                % this->messageDoneCalcCISMatrix.c_str() );
+   return value;
 }
 
 void ZindoS::CheckMatrixForce(const vector<int>& elecStates){
@@ -2372,6 +2579,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
                   const Atom& atomB = *molecule->GetAtom(b);
                   int firstAOIndexB = atomB.GetFirstAOIndex();
                   int lastAOIndexB  = atomB.GetLastAOIndex();
+                  double rAB = this->molecule->GetDistanceAtoms(atomA, atomB);
 
                   // calc. first derivative of overlapAOs.
                   this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs, atomA, atomB);
@@ -2389,7 +2597,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
                                              +atomB.GetCoreCharge()
                                              *atomicElectronPopulation[a])
                                              *this->GetNishimotoMatagaTwoEleInt1stDerivative(
-                                                    atomA, s, atomB, s, (CartesianType)i);
+                                                    atomA, s, atomB, s, rAB, static_cast<CartesianType>(i));
                   }
                   for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
                      OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
@@ -2408,7 +2616,8 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
                                                   -0.5*pow(this->orbitalElectronPopulation[mu][nu],2.0))
                                                   *this->GetNishimotoMatagaTwoEleInt1stDerivative(
                                                          atomA, orbitalMu, atomB, orbitalNu,
-                                                         (CartesianType)i);
+                                                         rAB,
+                                                         static_cast<CartesianType>(i));
                         }
                      }
                   }

--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -48,14 +48,15 @@ protected:
    virtual void SetMessages();
    virtual void SetEnableAtomTypes();
    virtual void CalcCISProperties();
-   virtual double GetElectronicTransitionDipoleMoment(int to, int from, MolDS_base::CartesianType axis,
-                                                      double const* const* fockMatrix,
-                                                      double const* const* matrixCIS,
-                                                      double const* const* const* cartesianMatrix,
-                                                      const MolDS_base::Molecule& molecule, 
-                                                      double const* const* orbitalElectronPopulation,
-                                                      double const* const* overlapAOs,
-                                                      double const* groundStateDipole) const;
+   virtual void CalcElectronicTransitionDipoleMoment(double* transitionDipoleMoment,
+                                                     int to, int from,
+                                                     double const* const* fockMatrix,
+                                                     double const* const* matrixCIS,
+                                                     double const* const* const* cartesianMatrix,
+                                                     const MolDS_base::Molecule& molecule, 
+                                                     double const* const* orbitalElectronPopulation,
+                                                     double const* const* overlapAOs,
+                                                     double const* groundStateDipole) const;
    virtual void CalcGammaAB(double** gammaAB, const MolDS_base::Molecule& molecule) const;
    virtual double GetFockDiagElement(const MolDS_base_atoms::Atom& atomA, 
                                      int indexAtomA, 
@@ -100,6 +101,19 @@ protected:
                                               double const* const* fockMatrix, 
                                               double const* const* gammaAB) const;
    virtual void CalcCISMatrix(double** matrixCIS) const;
+   double GetCISDiagElement(double const* energiesMO,
+                            double const* const* const* const* nishimotoMatagaMatrix,
+                            const MolDS_base::Molecule& molecule,
+                            double const* const* fockMatrix, 
+                            int moI,
+                            int moA) const;
+   double GetCISOffDiagElement(double const* const* const* const* nishimotoMatagaMatrix,
+                               const MolDS_base::Molecule& molecule,
+                               double const* const* fockMatrix, 
+                               int moI,
+                               int moA,
+                               int moJ,
+                               int moB) const;
    virtual void CalcForce(const std::vector<int>& elecStates);
    int GetSlaterDeterminantIndex(int activeOccIndex, int activeVirIndex) const;
    int GetActiveOccIndex(const MolDS_base::Molecule& molecule, int matrixCISIndex) const;
@@ -157,13 +171,13 @@ private:
    void CalcAtomicUnpairedPopulationCIS(double*** atomicUnpairedPopulationCIS,
                                         double const* const* const* orbitalElectronPopulationCIS, 
                                         const MolDS_base::Molecule& molecule) const; 
-   void CalcElectronicDipoleMomentsExcitedState(double*** electronicTransitionDipoleMoments,
-                                                double const* const* fockMatrix,
-                                                double const* const* matrixCIS,
-                                                double const* const* const* cartesianMatrix,
-                                                const MolDS_base::Molecule& molecule, 
-                                                double const* const* orbitalElectronPopulation,
-                                                double const* const* overlapAOs) const;
+   void CalcElectronicDipoleMomentsExcitedStates(double*** electronicTransitionDipoleMoments,
+                                                 double const* const* fockMatrix,
+                                                 double const* const* matrixCIS,
+                                                 double const* const* const* cartesianMatrix,
+                                                 const MolDS_base::Molecule& molecule, 
+                                                 double const* const* orbitalElectronPopulation,
+                                                 double const* const* overlapAOs) const;
    void CalcElectronicTransitionDipoleMoments(double*** electronicTransitionDipoleMoments,
                                               double const* const* fockMatrix,
                                               double const* const* matrixCIS,
@@ -175,11 +189,24 @@ private:
                                       MolDS_base::OrbitalType orbitalA, 
                                       const MolDS_base_atoms::Atom& atomB, 
                                       MolDS_base::OrbitalType orbitalB) const; // ref. [MN_1957] and (5a) in [AEZ_1986]
+   double GetNishimotoMatagaTwoEleInt(const MolDS_base_atoms::Atom& atomA, 
+                                      MolDS_base::OrbitalType orbitalA, 
+                                      const MolDS_base_atoms::Atom& atomB, 
+                                      MolDS_base::OrbitalType orbitalB,
+                                      const double rAB) const; // ref. [MN_1957] and (5a) in [AEZ_1986]
+   double GetNishimotoMatagaTwoEleInt1stDerivative(const MolDS_base_atoms::Atom& atomA, 
+                                                   MolDS_base::OrbitalType orbitalA, 
+                                                   const MolDS_base_atoms::Atom& atomB, 
+                                                   MolDS_base::OrbitalType orbitalB,
+                                                   MolDS_base::CartesianType axisA) const;// ref. [MN_1957] and (5a) in [AEZ_1986]
    double GetNishimotoMatagaTwoEleInt1stDerivative(const MolDS_base_atoms::Atom& atomA, 
                                                    MolDS_base::OrbitalType orbitalA, 
                                                    const MolDS_base_atoms::Atom& atomB, 
                                                    MolDS_base::OrbitalType orbitalB,
+                                                   const double rAB,
                                                    MolDS_base::CartesianType axisA) const;// ref. [MN_1957] and (5a) in [AEZ_1986]
+   void CalcNishimotoMatagaMatrix(double**** nishimotoMatagaMatrix, 
+                                  const MolDS_base::Molecule& molecule) const;
    void CalcRitzVector(double* ritzVector, 
                        double const* const* expansionVectors, 
                        double const* const* interactionMatrix, 

--- a/test/c2h6_pm3_directCIS_singlet.dat
+++ b/test/c2h6_pm3_directCIS_singlet.dat
@@ -1,6 +1,6 @@
 
 
-     >>>>>  Welcome to the MolDS world at 2013/1/23(Wed.) 7:0:44  <<<<<
+     >>>>>  Welcome to the MolDS world at 2013/1/27(Sun.) 15:14:52  <<<<<
 
 
 **********  START: Parse input  **********
@@ -53,28 +53,29 @@ Input terms:
 theory | pm3 | theory_end | scf | max_iter | 50 | rms_density | 0.000001 | damping_thresh | 1.0 | 
 damping_weight | 0.0 | diis_num_error_vect | 5 | diis_start_error | 0.1 | diis_end_error | 0.00000002 | scf_end | cis | 
 davidson | no | active_occ | 7 | active_vir | 7 | nstates | 49 | mulliken | 3 | 
-mulliken | 100 | mulliken | 20 | cis_end | geometry | c | 0.0000 | 0.0200 | 0.0000 | 
-c | 1.4938 | -0.0150 | 0.0020 | h | -0.3500 | 1.0411 | 0.0010 | h | -0.3681 | 
--0.5205 | -0.9200 | h | -0.3700 | -0.5208 | 0.9016 | h | 1.8519 | 0.5200 | -0.9007 | 
-h | 1.8300 | 0.5240 | 0.9100 | h | 1.8600 | -1.0401 | 0.0000 | geometry_end | 
+mulliken | 100 | mulliken | 20 | unpaired_electron_population | yes | cis_end | geometry | c | 0.0000 | 
+0.0200 | 0.0000 | c | 1.4938 | -0.0150 | 0.0020 | h | -0.3500 | 1.0411 | 0.0010 | 
+h | -0.3681 | -0.5205 | -0.9200 | h | -0.3700 | -0.5208 | 0.9016 | h | 1.8519 | 
+0.5200 | -0.9007 | h | 1.8300 | 0.5240 | 0.9100 | h | 1.8600 | -1.0401 | 0.0000 | 
+geometry_end | 
 
 **********  DONE: Parse input  ***********
 
 
 **********  START: PM3-SCF  **********
-	SCF iter 0:	RMS density = 5.291502622129180	DIIS error = 0.000000e+00
-	SCF iter 1:	RMS density = 1.910456125626944	DIIS error = 0.000000e+00
-	SCF iter 2:	RMS density = 1.044203238495391	DIIS error = 3.872759e-01
-	SCF iter 3:	RMS density = 0.569275110921829	DIIS error = 2.820935e-01
+	SCF iter 0:	RMS density = 5.291502622129181	DIIS error = 0.000000e+00
+	SCF iter 1:	RMS density = 1.910456125626941	DIIS error = 0.000000e+00
+	SCF iter 2:	RMS density = 1.044203238495393	DIIS error = 3.872759e-01
+	SCF iter 3:	RMS density = 0.569275110921830	DIIS error = 2.820935e-01
 	SCF iter 4:	RMS density = 0.308641545774631	DIIS error = 1.728758e-01
-	SCF iter 5:	RMS density = 0.166813060251975	DIIS error = 9.584873e-02
-	SCF iter 6:	RMS density = 0.000657312351891	DIIS error = 5.177961e-02 (DIIS was applied)
-	SCF iter 7:	RMS density = 0.000224796079707	DIIS error = 2.374135e-04 (DIIS was applied)
-	SCF iter 8:	RMS density = 0.000083564976841	DIIS error = 9.666685e-05 (DIIS was applied)
-	SCF iter 9:	RMS density = 0.000026986732846	DIIS error = 3.256104e-05 (DIIS was applied)
-	SCF iter 10:	RMS density = 0.000004938275825	DIIS error = 8.689524e-06 (DIIS was applied)
-	SCF iter 11:	RMS density = 0.000002003884721	DIIS error = 1.920254e-06 (DIIS was applied)
-	SCF iter 12:	RMS density = 0.000000169506350	DIIS error = 7.444364e-07 (DIIS was applied)
+	SCF iter 5:	RMS density = 0.166813060251973	DIIS error = 9.584873e-02
+	SCF iter 6:	RMS density = 0.000657312351887	DIIS error = 5.177961e-02 (DIIS was applied)
+	SCF iter 7:	RMS density = 0.000224796079661	DIIS error = 2.374135e-04 (DIIS was applied)
+	SCF iter 8:	RMS density = 0.000083564977008	DIIS error = 9.666685e-05 (DIIS was applied)
+	SCF iter 9:	RMS density = 0.000026986732836	DIIS error = 3.256104e-05 (DIIS was applied)
+	SCF iter 10:	RMS density = 0.000004938275875	DIIS error = 8.689524e-06 (DIIS was applied)
+	SCF iter 11:	RMS density = 0.000002003884716	DIIS error = 1.920254e-06 (DIIS was applied)
+	SCF iter 12:	RMS density = 0.000000169506349	DIIS error = 7.444364e-07 (DIIS was applied)
 
 
 
@@ -123,13 +124,13 @@ h | 1.8300 | 0.5240 | 0.9100 | h | 1.8600 | -1.0401 | 0.0000 | geometry_end |
 	Mulliken charge(SCF):	0	6	H	1.000000e+00	3.744313e-02
 	Mulliken charge(SCF):	0	7	H	1.000000e+00	4.258006e-02
 
-	Elapsed time(omp) for the SCF = 0.063144[s].
+	Elapsed time(omp) for the SCF = 0.460603[s].
 **********  DONE: PM3-SCF  **********
 
 
 **********  START: PM3-CIS  **********
 ----------- START: Calculation of the CIS matrix -----------
-	Elapsed time(omp) for the calc. of the CIS matrix = 0.087306[s].
+	Elapsed time(omp) for the calc. of the CIS matrix = 0.080024[s].
 ----------- DONE: Calculation of the CIS matrix -----------
 
 	======  START: Direct-CIS  =====
@@ -362,20 +363,39 @@ h | 1.8300 | 0.5240 | 0.9100 | h | 1.8600 | -1.0401 | 0.0000 | geometry_end |
 	Mulliken charge:	20	6	H	1.000000e+00	6.650230e-02
 	Mulliken charge:	20	7	H	1.000000e+00	-7.729498e-02
 
-
-	Elapsed time(omp) for the CIS = 0.229337[s].
+				| k-th eigenstate | i-th atom | atom type | Unpaired electron population[a.u.]| 
+	Unpaired electron population:	3	0	C	1.004203e+00
+	Unpaired electron population:	3	1	C	1.265832e+00
+	Unpaired electron population:	3	2	H	2.440939e-01
+	Unpaired electron population:	3	3	H	3.361993e-01
+	Unpaired electron population:	3	4	H	2.099692e-01
+	Unpaired electron population:	3	5	H	3.263566e-01
+	Unpaired electron population:	3	6	H	4.841361e-01
+	Unpaired electron population:	3	7	H	2.469543e-01
+
+	Unpaired electron population:	20	0	C	7.356038e-01
+	Unpaired electron population:	20	1	C	5.615877e-01
+	Unpaired electron population:	20	2	H	2.741903e-01
+	Unpaired electron population:	20	3	H	1.371744e-01
+	Unpaired electron population:	20	4	H	1.535949e-01
+	Unpaired electron population:	20	5	H	1.006022e-01
+	Unpaired electron population:	20	6	H	1.013346e-01
+	Unpaired electron population:	20	7	H	2.323766e-01
+
+
+	Elapsed time(omp) for the CIS = 0.219919[s].
 **********  DONE: PM3-CIS  **********
 
 
 	Summary for memory usage:
 		Max Heap: 0.296152[MB].
-		Current Heap(Leaked): 0.000000[MB].
+		Current Heap(Leaked): 0.000144[MB].
 
 
      >>>>>  The MolDS finished normally!  <<<<<
-     >>>>>  CPU time: 0.3[s].  <<<<<
-     >>>>>  Elapsed time: 0[s].  <<<<<
-     >>>>>  Elapsed time(OMP): 0.297168[s].  <<<<<
+     >>>>>  CPU time: 0.26[s].  <<<<<
+     >>>>>  Elapsed time: 1[s].  <<<<<
+     >>>>>  Elapsed time(OMP): 0.685674[s].  <<<<<
      >>>>>  See you.  <<<<<
 
 

--- a/test/c2h6_pm3_directCIS_singlet.in
+++ b/test/c2h6_pm3_directCIS_singlet.in
@@ -20,6 +20,7 @@ CIS
    mulliken 3
    mulliken 100
    mulliken 20
+   unpaired_electron_population yes
 CIS_END
 
 GEOMETRY