Commit: 82f9cbd8ff2b6cfe0642b380a6348ce7caaf9454 - MolDS (git)

Commit MetaInfo

• Révision: 82f9cbd8ff2b6cfe0642b380a6348ce7caaf9454 (tree)
• l'heure: 2013-12-31 17:19:34
• Auteur: mikiya_fujii <mikiya_fujii@1136...>
• Commiter: mikiya_fujii

Message de Log

Refactring of access to Molecule::atomVect and Molecule::epcVect #32752

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

Change Summary

• modified: src/am1/Am1.cpp
• modified: src/base/InputParser.cpp
• modified: src/base/Molecule.cpp
• modified: src/base/Molecule.h
• modified: src/base/loggers/DensityLogger.cpp
• modified: src/base/loggers/MOLogger.cpp
• modified: src/cndo/Cndo2.cpp
• modified: src/indo/Indo.cpp
• modified: src/mc/MC.cpp
• modified: src/md/MD.cpp
• modified: src/mndo/Mndo.cpp
• modified: src/nasco/NASCO.cpp
• modified: src/optimization/BFGS.cpp
• modified: src/optimization/ConjugateGradient.cpp
• modified: src/optimization/GEDIIS.cpp
• modified: src/optimization/Optimizer.cpp
• modified: src/pm3/Pm3Pddg.cpp
• modified: src/rpmd/RPMD.cpp
• modified: src/zindo/ZindoS.cpp

Modification

--- a/src/am1/Am1.cpp
+++ b/src/am1/Am1.cpp
@@ -140,8 +140,8 @@ double Am1::GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const{
    double mndoTerm = Mndo::GetDiatomCoreRepulsionEnergy(indexAtomA, indexAtomB);
 
    // additional term, Eq. (4) in [S_1989].
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double distance   = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
    double ang2AU     = Parameters::GetInstance()->GetAngstrom2AU();
    double alphaA     = atomA.GetNddoAlpha(this->theory);
@@ -174,8 +174,8 @@ double Am1::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
 
    // additional term, first derivative of eq. (4) in [S_1989]
    double ang2AU     = Parameters::GetInstance()->GetAngstrom2AU();
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double alphaA     = atomA.GetNddoAlpha(this->theory);
    double alphaB     = atomB.GetNddoAlpha(this->theory);
    double distance   = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
@@ -215,8 +215,8 @@ double Am1::GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
 
    // additional term, first derivative of eq. (4) in [S_1989]
    double ang2AU     = Parameters::GetInstance()->GetAngstrom2AU();
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double alphaA     = atomA.GetNddoAlpha(this->theory);
    double alphaB     = atomB.GetNddoAlpha(this->theory);
    double distance   = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);

--- a/src/base/InputParser.cpp
+++ b/src/base/InputParser.cpp
@@ -487,7 +487,7 @@ int InputParser::ParseMolecularGeometry(Molecule* molecule, vector<string>* inpu
       else if((*inputTerms)[parseIndex] == "s"){
          atomType = S;
       }
-      int index = molecule->GetNumberAtoms();
+      int index = molecule->GetAtomVect().size();
       Atom* atom = AtomFactory::Create(atomType, index, x, y, z);
       molecule->AddAtom(atom);
       parseIndex += 4;
@@ -508,7 +508,7 @@ int InputParser::ParseEpcsConfiguration(Molecule* molecule, vector<string>* inpu
          parseIndex += 2;
       }
       AtomType atomType = EPC;
-      int index = molecule->GetNumberEpcs();
+      int index = molecule->GetEpcVect().size();
       Atom* atom = AtomFactory::Create(atomType, index, x, y, z, charge);
       molecule->AddEpc(atom);
    }
@@ -1360,7 +1360,7 @@ void InputParser::ValidateVdWConditions() const{
 }
 
 void InputParser::ValidateEpcConditions(const Molecule& molecule) const{
-   if(molecule.GetNumberEpcs()<=0){return;}
+   if(molecule.GetEpcVect().empty()){return;}
    TheoryType theory = Parameters::GetInstance()->GetCurrentTheory();
    // Validate theory
    if(theory == MNDO || 

--- a/src/base/Molecule.cpp
+++ b/src/base/Molecule.cpp
@@ -223,12 +223,10 @@ void Molecule::Finalize(vector<Atom*>** atomVect,
 }
 
 void Molecule::SetMessages(){
-   this->errorMessageGetAtomNull = "Error in base::Molecule::GetAtom: atomVect is NULL.\n";
-   this->errorMessageGetEPCNull  = "Error in base::Molecule::GetEnviromentalPointCharge: enviromentalPointChargeVect is NULL.\n";
+   this->errorMessageGetAtomVectNull  = "Error in base::Molecule::GetAtomVect: atomVect is NULL.\n";
+   this->errorMessageGetEPCVectNull  = "Error in base::Molecule::GetEpcVect: epcVect is NULL.\n";
    this->errorMessageAddAtomNull = "Error in base::Molecule::AddAtom: atomVect is NULL.\n";
-   this->errorMessageAddEPCNull  = "Error in base::Molecule::AddEnviromentalPointCharge: enviromentalPointChargeVect is NULL.\n";
-   this->errorMessageGetNumberAtomsNull = "Error in base::Molecule::GetNumberAtoms: atomVect is NULL.\n";
-   this->errorMessageGetNumberEPCsNull  = "Error in base::Molecule::GetNumberEnviromentalPointChargess: epcVect is NULL.\n";
+   this->errorMessageAddEPCNull  = "Error in base::Molecule::AddEnviromentalPointCharge: epcVect is NULL.\n";
    this->errorMessageGetXyzCOMNull = "Error in base::Molecule::GetXyzCOM: xyzCOM is NULL.\n";
    this->errorMessageGetXyzCOCNull = "Error in base::Molecule::GetXyzCOC: xyzCOC is NULL.\n";
    this->errorMessageCalcXyzCOMNull = "Error in base::Molecule::CalcXyzCOM: xyzCOM is NULL.\n";
@@ -832,9 +830,9 @@ void Molecule::OutputTranslatingConditions(double const* translatingDifference)
 }
 
 void Molecule::SynchronizeConfigurationTo(const Molecule& ref){
-   for(int a=0; a<this->GetNumberAtoms(); a++){
-      Atom& atom = *this->GetAtom(a);
-      const Atom& refAtom = *ref.GetAtom(a);
+   for(int a=0; a<this->GetAtomVect().size(); a++){
+      Atom& atom = *this->GetAtomVect()[a];
+      const Atom& refAtom = *ref.GetAtomVect()[a];
       for(int i=0; i<CartesianType_end; i++){
          atom.GetXyz()[i] = refAtom.GetXyz()[i];
       }
@@ -843,9 +841,9 @@ void Molecule::SynchronizeConfigurationTo(const Molecule& ref){
 }
 
 void Molecule::SynchronizeMomentaTo(const Molecule& ref){
-   for(int a=0; a<this->GetNumberAtoms(); a++){
-      Atom& atom = *this->GetAtom(a);
-      const Atom& refAtom = *ref.GetAtom(a);
+   for(int a=0; a<this->GetAtomVect().size(); a++){
+      Atom& atom = *this->GetAtomVect()[a];
+      const Atom& refAtom = *ref.GetAtomVect()[a];
       for(int i=0; i<CartesianType_end; i++){
          atom.GetPxyz()[i] = refAtom.GetPxyz()[i];
       }
@@ -853,9 +851,9 @@ void Molecule::SynchronizeMomentaTo(const Molecule& ref){
 }
 
 void Molecule::SynchronizePhaseSpacePointTo(const Molecule& ref){
-   for(int a=0; a<this->GetNumberAtoms(); a++){
-      Atom& atom = *this->GetAtom(a);
-      const Atom& refAtom = *ref.GetAtom(a);
+   for(int a=0; a<this->GetAtomVect().size(); a++){
+      Atom& atom = *this->GetAtomVect()[a];
+      const Atom& refAtom = *ref.GetAtomVect()[a];
       for(int i=0; i<CartesianType_end; i++){
          atom.GetXyz() [i] = refAtom.GetXyz() [i];
          atom.GetPxyz()[i] = refAtom.GetPxyz()[i];
@@ -865,19 +863,19 @@ void Molecule::SynchronizePhaseSpacePointTo(const Molecule& ref){
 }
 
 void Molecule::BroadcastConfigurationToAllProcesses(int root) const{
-   int numTransported = this->GetNumberAtoms()*CartesianType_end;
+   int numTransported = this->GetAtomVect().size()*CartesianType_end;
    double* tmp=NULL;
    try{
       MolDS_base::MallocerFreer::GetInstance()->Malloc<double>(&tmp, numTransported);
-      for(int a=0; a<this->GetNumberAtoms(); a++){
-         Atom& atom = *this->GetAtom(a);
+      for(int a=0; a<this->GetAtomVect().size(); a++){
+         Atom& atom = *this->GetAtomVect()[a];
          for(int i=0; i<CartesianType_end; i++){
             tmp[a*CartesianType_end+i] = atom.GetXyz()[i];
          }
       }
       MolDS_mpi::MpiProcess::GetInstance()->Broadcast(tmp, numTransported, root);
-      for(int a=0; a<this->GetNumberAtoms(); a++){
-         Atom& atom = *this->GetAtom(a);
+      for(int a=0; a<this->GetAtomVect().size(); a++){
+         Atom& atom = *this->GetAtomVect()[a];
          for(int i=0; i<CartesianType_end; i++){
             atom.GetXyz()[i] = tmp[a*CartesianType_end+i];
          }
@@ -891,19 +889,19 @@ void Molecule::BroadcastConfigurationToAllProcesses(int root) const{
 }
 
 void Molecule::BroadcastMomentaToAllProcesses(int root) const{
-   int numTransported = this->GetNumberAtoms()*CartesianType_end;
+   int numTransported = this->GetAtomVect().size()*CartesianType_end;
    double* tmp=NULL;
    try{
       MolDS_base::MallocerFreer::GetInstance()->Malloc<double>(&tmp, numTransported);
-      for(int a=0; a<this->GetNumberAtoms(); a++){
-         Atom& atom = *this->GetAtom(a);
+      for(int a=0; a<this->GetAtomVect().size(); a++){
+         Atom& atom = *this->GetAtomVect()[a];
          for(int i=0; i<CartesianType_end; i++){
             tmp[a*CartesianType_end+i] = atom.GetPxyz()[i];
          }
       }
       MolDS_mpi::MpiProcess::GetInstance()->Broadcast(tmp, numTransported, root);
-      for(int a=0; a<this->GetNumberAtoms(); a++){
-         Atom& atom = *this->GetAtom(a);
+      for(int a=0; a<this->GetAtomVect().size(); a++){
+         Atom& atom = *this->GetAtomVect()[a];
          for(int i=0; i<CartesianType_end; i++){
             atom.GetPxyz()[i] = tmp[a*CartesianType_end+i];
          }
@@ -917,12 +915,12 @@ void Molecule::BroadcastMomentaToAllProcesses(int root) const{
 }
 
 void Molecule::BroadcastPhaseSpacePointToAllProcesses(int root) const{
-   int numTransported = 2*this->GetNumberAtoms()*CartesianType_end;
+   int numTransported = 2*this->GetAtomVect().size()*CartesianType_end;
    double* tmp=NULL;
    try{
       MolDS_base::MallocerFreer::GetInstance()->Malloc<double>(&tmp, numTransported);
-      for(int a=0; a<this->GetNumberAtoms(); a++){
-         Atom& atom = *this->GetAtom(a);
+      for(int a=0; a<this->GetAtomVect().size(); a++){
+         Atom& atom = *this->GetAtomVect()[a];
          for(int i=0; i<CartesianType_end; i++){
             int k = a*CartesianType_end+i;
             tmp[2*k  ] = atom.GetXyz() [i];
@@ -930,8 +928,8 @@ void Molecule::BroadcastPhaseSpacePointToAllProcesses(int root) const{
          }
       }
       MolDS_mpi::MpiProcess::GetInstance()->Broadcast(tmp, numTransported, root);
-      for(int a=0; a<this->GetNumberAtoms(); a++){
-         Atom& atom = *this->GetAtom(a);
+      for(int a=0; a<this->GetAtomVect().size(); a++){
+         Atom& atom = *this->GetAtomVect()[a];
          for(int i=0; i<CartesianType_end; i++){
             int k = a*CartesianType_end+i;
             atom.GetXyz() [i] = tmp[2*k  ];

--- a/src/base/Molecule.h
+++ b/src/base/Molecule.h
@@ -26,29 +26,17 @@ public:
    explicit Molecule(const Molecule& rhs);
    Molecule& operator=(const Molecule& rhs);
    ~Molecule();
-   inline int GetNumberAtoms() const{
+   inline const std::vector<MolDS_base_atoms::Atom*>& GetAtomVect() const{
 #ifdef MOLDS_DBG
-      if(this->atomVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetNumberAtomsNull);
+      if(this->atomVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetAtomVectNull);
 #endif
-      return this->atomVect->size();
+      return *this->atomVect;
    }
-   inline MolDS_base_atoms::Atom* GetAtom(int atomIndex) const{
+   inline const std::vector<MolDS_base_atoms::Atom*>& GetEpcVect() const{
 #ifdef MOLDS_DBG
-      if(this->atomVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetAtomNull);
+      if(this->epcVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetEpcVectNull);
 #endif
-      return (*this->atomVect)[atomIndex];
-   }
-   inline int GetNumberEpcs() const{
-#ifdef MOLDS_DBG
-      if(this->epcVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetNumberEPCsNull);
-#endif
-      return this->epcVect->size();
-   }
-   inline MolDS_base_atoms::Atom* GetEpc(int epcIndex) const{
-#ifdef MOLDS_DBG
-      if(this->epcVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetEPCNull);
-#endif
-      return (*this->epcVect)[epcIndex];
+      return *this->epcVect;
    }
    void AddAtom(MolDS_base_atoms::Atom* atom);
    void AddEpc(MolDS_base_atoms::Atom* epc);
@@ -128,12 +116,10 @@ private:
                                  double rotatingAngle, 
                                  MolDS_base::EularAngle rotatingEularAngles)const;
    void OutputTranslatingConditions(double const* translatingDifference) const;
-   std::string errorMessageGetAtomNull;
-   std::string errorMessageGetEPCNull;
+   std::string errorMessageGetAtomVectNull;
+   std::string errorMessageGetEPCVectNull;
    std::string errorMessageAddAtomNull;
    std::string errorMessageAddEPCNull;
-   std::string errorMessageGetNumberAtomsNull;
-   std::string errorMessageGetNumberEPCsNull;
    std::string errorMessageGetXyzCOMNull;
    std::string errorMessageGetXyzCOCNull;
    std::string errorMessageCalcXyzCOMNull;

--- a/src/base/loggers/DensityLogger.cpp
+++ b/src/base/loggers/DensityLogger.cpp
@@ -257,8 +257,8 @@ double DensityLogger::GetMOValue(int moIndex,
                                  double const* const* forckMatrix,
                                  double x, double y, double z) const{
    double moValue = 0.0;
-   for(int a=0; a<this->molecule->GetNumberAtoms(); a++){
-      Atom* atomA = this->molecule->GetAtom(a);
+   for(int a=0; a<this->molecule->GetAtomVect().size(); a++){
+      Atom* atomA = this->molecule->GetAtomVect()[a];
       int firstAOIndexA = atomA->GetFirstAOIndex();
       int numberAOsA = atomA->GetValenceSize();
       for(int mu=firstAOIndexA; mu<firstAOIndexA+numberAOsA; mu++){
@@ -281,10 +281,10 @@ void DensityLogger::OutputHeaderToFile(ofstream& ofs, double const* origin, doub
    // output header to the cube file
    ofs << this->messageCubeHeaderComment1;
    ofs << this->messageCubeHeaderComment2;
-   sprintf(data,"\t%d\t%e\t%e\t%e\n", this->molecule->GetNumberAtoms(),
-                                      origin[XAxis], 
-                                      origin[YAxis], 
-                                      origin[ZAxis]);
+   sprintf(data,"\t%ld\t%e\t%e\t%e\n", this->molecule->GetAtomVect().size(),
+                                       origin[XAxis], 
+                                       origin[YAxis], 
+                                       origin[ZAxis]);
    ofs << string(data);
    memset(data,0,sizeof(data));
    sprintf(data,"\t%d\t%e\t%e\t%e\n", gridNumber[XAxis], dx, 0.0, 0.0);
@@ -300,8 +300,8 @@ void DensityLogger::OutputHeaderToFile(ofstream& ofs, double const* origin, doub
 void DensityLogger::OutputMoleculeToFile(ofstream& ofs, const Molecule& molecule) const{
    char data[1000] = "";
    // output molecule to the cube file
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      const Atom& atomA = *molecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      const Atom& atomA = *molecule.GetAtomVect()[a];
       memset(data,0,sizeof(data));
       sprintf(data,"\t%d\t%d\t%e\t%e\t%e\n", atomA.GetAtomType()+1, 
                                        atomA.GetNumberValenceElectrons(),

--- a/src/base/loggers/MOLogger.cpp
+++ b/src/base/loggers/MOLogger.cpp
@@ -165,10 +165,10 @@ void MOLogger::OutputHeaderToFile(ofstream& ofs, double const* origin, double dx
    // output header to the cube file
    ofs << this->messageCubeHeaderComment1;
    ofs << this->messageCubeHeaderComment2;
-   sprintf(data,"\t%d\t%e\t%e\t%e\n", this->molecule->GetNumberAtoms(),
-                                      origin[XAxis], 
-                                      origin[YAxis], 
-                                      origin[ZAxis]);
+   sprintf(data,"\t%ld\t%e\t%e\t%e\n", this->molecule->GetAtomVect().size(),
+                                       origin[XAxis], 
+                                       origin[YAxis], 
+                                       origin[ZAxis]);
    ofs << string(data);
    memset(data,0,sizeof(data));
    sprintf(data,"\t%d\t%e\t%e\t%e\n", gridNumber[XAxis], dx, 0.0, 0.0);
@@ -184,8 +184,8 @@ void MOLogger::OutputHeaderToFile(ofstream& ofs, double const* origin, double dx
 void MOLogger::OutputMoleculeToFile(ofstream& ofs, const Molecule& molecule) const{
    char data[1000] = "";
    // output molecule to the cube file
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      Atom* atomA = molecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      Atom* atomA = molecule.GetAtomVect()[a];
       memset(data,0,sizeof(data));
       sprintf(data,"\t%d\t%d\t%e\t%e\t%e\n", atomA->GetAtomType()+1, 
                                        atomA->GetNumberValenceElectrons(),
@@ -228,8 +228,8 @@ double MOLogger::GetMoValue(int moIndex,
                             double y, 
                             double z) const{
    double moValue = 0.0;
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      Atom* atomA = molecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      Atom* atomA = molecule.GetAtomVect()[a];
       int firstAOIndexA = atomA->GetFirstAOIndex();
       int numberAOsA = atomA->GetValenceSize();
       for(int mu=firstAOIndexA; mu<firstAOIndexA+numberAOsA; mu++){

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -117,7 +117,7 @@ Cndo2::~Cndo2(){
                                               this->molecule->GetTotalNumberAOs(),
                                               this->molecule->GetTotalNumberAOs());
    MallocerFreer::GetInstance()->Free<double>(&this->atomicElectronPopulation, 
-                                              this->molecule->GetNumberAtoms());
+                                              this->molecule->GetAtomVect().size());
    MallocerFreer::GetInstance()->Free<double>(&this->overlapAOs, 
                                               this->molecule->GetTotalNumberAOs(),
                                               this->molecule->GetTotalNumberAOs());
@@ -136,8 +136,8 @@ Cndo2::~Cndo2(){
    MallocerFreer::GetInstance()->Free<double>(&this->coreDipoleMoment, 
                                               CartesianType_end);
    MallocerFreer::GetInstance()->Free<double>(&this->gammaAB, 
-                                              this->molecule->GetNumberAtoms(),
-                                              this->molecule->GetNumberAtoms());
+                                              this->molecule->GetAtomVect().size(),
+                                              this->molecule->GetAtomVect().size());
    //this->OutputLog("cndo deleted\n");
 }
 
@@ -293,7 +293,7 @@ void Cndo2::SetMolecule(Molecule* molecule){
                                                 this->molecule->GetTotalNumberAOs(), 
                                                 this->molecule->GetTotalNumberAOs());
    MallocerFreer::GetInstance()->Malloc<double>(&this->atomicElectronPopulation,
-                                                this->molecule->GetNumberAtoms());
+                                                this->molecule->GetAtomVect().size());
    MallocerFreer::GetInstance()->Malloc<double>(&this->overlapAOs, 
                                                 this->molecule->GetTotalNumberAOs(), 
                                                 this->molecule->GetTotalNumberAOs());
@@ -313,8 +313,8 @@ void Cndo2::SetMolecule(Molecule* molecule){
                                                 CartesianType_end);
    if(this->theory == CNDO2 || this->theory == INDO){
       MallocerFreer::GetInstance()->Malloc<double>(&this->gammaAB,
-                                                   this->molecule->GetNumberAtoms(), 
-                                                   this->molecule->GetNumberAtoms());
+                                                   this->molecule->GetAtomVect().size(), 
+                                                   this->molecule->GetAtomVect().size());
    }
 }
 
@@ -327,8 +327,8 @@ void Cndo2::CheckNumberValenceElectrons(const Molecule& molecule) const{
 }
 
 void Cndo2::CheckEnableAtomType(const Molecule& molecule) const{
-   for(int i=0; i<molecule.GetNumberAtoms(); i++){
-      AtomType atomType = molecule.GetAtom(i)->GetAtomType();
+   for(int i=0; i<molecule.GetAtomVect().size(); i++){
+      AtomType atomType = molecule.GetAtomVect()[i]->GetAtomType();
       bool enable = false;
       for(int j=0; j<this->enableAtomTypes.size(); j++){
          if(atomType == this->enableAtomTypes[j]){
@@ -348,18 +348,18 @@ void Cndo2::CheckEnableAtomType(const Molecule& molecule) const{
 void Cndo2::CalcCoreRepulsionEnergy(){
    // interaction between atoms
    double energy = 0.0;
-   for(int i=0; i<this->molecule->GetNumberAtoms(); i++){
-      for(int j=i+1; j<this->molecule->GetNumberAtoms(); j++){
+   for(int i=0; i<this->molecule->GetAtomVect().size(); i++){
+      for(int j=i+1; j<this->molecule->GetAtomVect().size(); j++){
          energy += this->GetDiatomCoreRepulsionEnergy(i, j);
       }
    }
    this->coreRepulsionEnergy = energy;
 
    // interaction between atoms and epcs
-   if(this->molecule->GetNumberEpcs()<=0){return;}
+   if(this->molecule->GetEpcVect().empty()){return;}
    energy = 0.0;
-   for(int i=0; i<this->molecule->GetNumberAtoms(); i++){
-      for(int j=0; j<this->molecule->GetNumberEpcs(); j++){
+   for(int i=0; i<this->molecule->GetAtomVect().size(); i++){
+      for(int j=0; j<this->molecule->GetEpcVect().size(); j++){
          energy += this->GetAtomCoreEpcCoulombEnergy(i, j);
       }
    }
@@ -368,8 +368,8 @@ void Cndo2::CalcCoreRepulsionEnergy(){
 }
 
 double Cndo2::GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
    return atomA.GetCoreCharge()*atomB.GetCoreCharge()/distance; 
 }
@@ -384,8 +384,8 @@ double Cndo2::GetAtomCoreEpcCoulombEnergy(int indexAtom, int indexEpc) const{
 double Cndo2::GetDiatomCoreRepulsion1stDerivative(int indexAtomA, int indexAtomB, 
                                                   CartesianType axisA) const{
    double value=0.0;
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
    value = atomA.GetCoreCharge()*atomB.GetCoreCharge();
    value *= (atomA.GetXyz()[axisA] - atomB.GetXyz()[axisA])/distance;
@@ -407,8 +407,8 @@ double Cndo2::GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
 // See (2) in [G_2004] ((11) in [G_2006])
 void Cndo2::CalcVdWCorrectionEnergy(){
    double value = 0.0;
-   for(int i=0; i<this->molecule->GetNumberAtoms(); i++){
-      for(int j=i+1; j<this->molecule->GetNumberAtoms(); j++){
+   for(int i=0; i<this->molecule->GetAtomVect().size(); i++){
+      for(int j=i+1; j<this->molecule->GetAtomVect().size(); j++){
          value += this->GetDiatomVdWCorrectionEnergy(i, j);
       }
    }
@@ -442,8 +442,8 @@ double Cndo2::GetVdwDampingValue2ndDerivative(double vdWDistance, double distanc
 
 // See (2) in [G_2004] ((11) in [G_2006])
 double Cndo2::GetDiatomVdWCorrectionEnergy(int indexAtomA, int indexAtomB) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
    double vdWDistance = atomA.GetVdWRadii() + atomB.GetVdWRadii();
    double vdWCoefficients = 2.0*atomA.GetVdWCoefficient()*atomB.GetVdWCoefficient()
@@ -458,8 +458,8 @@ double Cndo2::GetDiatomVdWCorrectionEnergy(int indexAtomA, int indexAtomB) const
 // See (2) in [G_2004] ((11) in [G_2006]).
 double Cndo2::GetDiatomVdWCorrection1stDerivative(int indexAtomA, int indexAtomB, 
                                                   CartesianType axisA) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
    double vdWDistance = atomA.GetVdWRadii() + atomB.GetVdWRadii();
    double vdWCoefficients = 2.0*atomA.GetVdWCoefficient()*atomB.GetVdWCoefficient()
@@ -484,8 +484,8 @@ double Cndo2::GetDiatomVdWCorrection2ndDerivative(int indexAtomA,
                                                      int indexAtomB, 
                                                      CartesianType axisA1,
                                                      CartesianType axisA2) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
    double dCartesian1 = atomA.GetXyz()[axisA1] - atomB.GetXyz()[axisA1];
    double dCartesian2 = atomA.GetXyz()[axisA2] - atomB.GetXyz()[axisA2];
@@ -1013,13 +1013,13 @@ void Cndo2::OutputSCFEnergies() const{
    this->OutputLog(boost::format("%s\t%e\t%e\n") % this->messageElecEnergy
                                                  % this->elecSCFEnergy
                                                  % (this->elecSCFEnergy/eV2AU));
-   if(Parameters::GetInstance()->RequiresVdWSCF() && this->molecule->GetNumberEpcs()<=0){
+   if(Parameters::GetInstance()->RequiresVdWSCF() && this->molecule->GetEpcVect().empty()){
       this->OutputLog(this->messageNoteElecEnergyVdW);
    }
-   else if(Parameters::GetInstance()->RequiresVdWSCF() && 0<this->molecule->GetNumberEpcs()){
+   else if(Parameters::GetInstance()->RequiresVdWSCF() && !this->molecule->GetEpcVect().empty()){
       this->OutputLog(this->messageNoteElecEnergyEpcVdW);
    }
-   else if(!Parameters::GetInstance()->RequiresVdWSCF() && 0<this->molecule->GetNumberEpcs()){
+   else if(!Parameters::GetInstance()->RequiresVdWSCF() && !this->molecule->GetEpcVect().empty()){
       this->OutputLog(this->messageNoteElecEnergyEpc);
    }
    else{
@@ -1033,7 +1033,7 @@ void Cndo2::OutputSCFEnergies() const{
                                                    % (this->coreRepulsionEnergy/eV2AU));
 
    // output coulomb interaction between atoms and epcs
-   if(0<this->molecule->GetNumberEpcs()){
+   if(!this->molecule->GetEpcVect().empty()){
       this->OutputLog(this->messageCoreEpcCoulombTitle);
       this->OutputLog(boost::format("%s\t%e\t%e\n\n") % this->messageCoreEpcCoulomb
                                                       % this->coreEpcCoulombEnergy 
@@ -1113,8 +1113,8 @@ void Cndo2::OutputSCFDipole() const{
 void Cndo2::OutputSCFMulliken() const{
    int groundState = 0;
    this->OutputLog(this->messageMullikenAtomsTitle);
-   for(int a=0; a<this->molecule->GetNumberAtoms(); a++){
-      Atom* atom = this->molecule->GetAtom(a);
+   for(int a=0; a<this->molecule->GetAtomVect().size(); a++){
+      Atom* atom = this->molecule->GetAtomVect()[a];
       this->OutputLog(boost::format("%s\t%d\t%d\t%s\t%e\t%e\n") % this->messageMullikenAtomsSCF
                                                                 % groundState
                                                                 % a
@@ -1129,7 +1129,7 @@ void Cndo2::OutputNormalModes(double const* const* normalModes,
                               double const* normalForceConstants, 
                               const Molecule& molecule) const{
 
-   int hessianDim = CartesianType_end*molecule.GetNumberAtoms();
+   int hessianDim = CartesianType_end*molecule.GetAtomVect().size();
    double ang2AU = Parameters::GetInstance()->GetAngstrom2AU();
    double kayser2AU = Parameters::GetInstance()->GetKayser2AU();
 
@@ -1149,8 +1149,8 @@ void Cndo2::OutputNormalModes(double const* const* normalModes,
                                                                % sqrt(fabs(normalForceConstants[i]))
                                                                % (sqrt(fabs(normalForceConstants[i]))/kayser2AU));
       // normal modes
-      for(int a=0; a<molecule.GetNumberAtoms(); a++){
-         const double sqrtCoreMass = sqrt(molecule.GetAtom(a)->GetCoreMass());
+      for(int a=0; a<molecule.GetAtomVect().size(); a++){
+         const double sqrtCoreMass = sqrt(molecule.GetAtomVect()[a]->GetCoreMass());
          for(int j=XAxis; j<CartesianType_end; j++){
             int hessianIndex = CartesianType_end*a+j;
             this->OutputLog(boost::format("\t%e") % normalModes[i][hessianIndex]);
@@ -1178,8 +1178,8 @@ void Cndo2::OutputNormalModes(double const* const* normalModes,
                                                                % (sqrt(fabs(normalForceConstants[i]))/kayser2AU));
 
       double normSquare=0.0;
-      for(int a=0; a<molecule.GetNumberAtoms(); a++){
-         const double sqrtCoreMass = sqrt(molecule.GetAtom(a)->GetCoreMass());
+      for(int a=0; a<molecule.GetAtomVect().size(); a++){
+         const double sqrtCoreMass = sqrt(molecule.GetAtomVect()[a]->GetCoreMass());
          for(int j=XAxis; j<CartesianType_end; j++){
             int hessianIndex = CartesianType_end*a+j;
             normSquare += pow(normalModes[i][hessianIndex]/(sqrtCoreMass*ang2AU),2.0);
@@ -1188,8 +1188,8 @@ void Cndo2::OutputNormalModes(double const* const* normalModes,
       double norm = sqrt(normSquare);
 
       // normal modes
-      for(int a=0; a<molecule.GetNumberAtoms(); a++){
-         const double sqrtCoreMass = sqrt(molecule.GetAtom(a)->GetCoreMass());
+      for(int a=0; a<molecule.GetAtomVect().size(); a++){
+         const double sqrtCoreMass = sqrt(molecule.GetAtomVect()[a]->GetCoreMass());
          for(int j=XAxis; j<CartesianType_end; j++){
             int hessianIndex = CartesianType_end*a+j;
             this->OutputLog(boost::format("\t%e") % (normalModes[i][hessianIndex]/(sqrtCoreMass*ang2AU*norm)));
@@ -1248,7 +1248,7 @@ void Cndo2::CalcElecSCFEnergy(double* elecSCFEnergy,
                                                    totalNumberAOs, 
                                                    totalNumberAOs);
       MallocerFreer::GetInstance()->Malloc<double>(&dammyAtomicElectronPopulation,
-                                                   molecule.GetNumberAtoms());
+                                                   molecule.GetAtomVect().size());
       bool isGuess = false;
       this->CalcFockMatrix(fMatrix, 
                            molecule, 
@@ -1328,7 +1328,7 @@ void Cndo2::FreeElecEnergyMatrices(double*** fMatrix,
                                               totalNumberAOs,
                                               totalNumberAOs);
    MallocerFreer::GetInstance()->Free<double>(dammyAtomicElectronPopulation,
-                                              this->molecule->GetNumberAtoms());
+                                              this->molecule->GetAtomVect().size());
 }
 
 // The order of moI, moJ, moK, moL is consistent with Eq. (9) in [RZ_1973]
@@ -1337,13 +1337,13 @@ double Cndo2::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
                                           double const* const* fockMatrix, 
                                           double const* const* gammaAB) const{
    double value = 0.0;
-   for(int A=0; A<molecule.GetNumberAtoms(); A++){
-      const Atom& atomA = *molecule.GetAtom(A);
+   for(int A=0; A<molecule.GetAtomVect().size(); A++){
+      const Atom& atomA = *molecule.GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
 
-      for(int B=0; B<molecule.GetNumberAtoms(); B++){
-         const Atom& atomB = *molecule.GetAtom(B);
+      for(int B=0; B<molecule.GetAtomVect().size(); B++){
+         const Atom& atomB = *molecule.GetAtomVect()[B];
          int firstAOIndexB = atomB.GetFirstAOIndex();
          int lastAOIndexB  = atomB.GetLastAOIndex();
          double gamma = gammaAB[A][B];
@@ -1434,7 +1434,7 @@ void Cndo2::CalcFockMatrix(double** fockMatrix,
                            double const* const* const* const* const* const* twoElecsTwoAtomCores, 
                            bool isGuess) const{
    int totalNumberAOs   = molecule.GetTotalNumberAOs();
-   int totalNumberAtoms = molecule.GetNumberAtoms();
+   int totalNumberAtoms = molecule.GetAtomVect().size();
    MallocerFreer::GetInstance()->Initialize<double>(fockMatrix, totalNumberAOs, totalNumberAOs);
 
    // MPI setting of each rank
@@ -1446,7 +1446,7 @@ void Cndo2::CalcFockMatrix(double** fockMatrix,
    boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
 
    for(int A=0; A<totalNumberAtoms; A++){
-      const Atom& atomA = *molecule.GetAtom(A);
+      const Atom& atomA = *molecule.GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
       for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
@@ -1455,7 +1455,7 @@ void Cndo2::CalcFockMatrix(double** fockMatrix,
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
             for(int B=A; B<totalNumberAtoms; B++){
                try{
-                  const Atom& atomB = *molecule.GetAtom(B);
+                  const Atom& atomB = *molecule.GetAtomVect()[B];
                   int firstAOIndexB = atomB.GetFirstAOIndex();
                   int lastAOIndexB  = atomB.GetLastAOIndex();
                   for(int nu=firstAOIndexB; nu<=lastAOIndexB; nu++){
@@ -1553,9 +1553,9 @@ double Cndo2::GetFockDiagElement(const Atom& atomA,
       value += temp*gammaAB[indexAtomA][indexAtomA];
 
       temp = 0.0;
-      for(int BB=0; BB<molecule.GetNumberAtoms(); BB++){
+      for(int BB=0; BB<molecule.GetAtomVect().size(); BB++){
          if(BB != indexAtomA){
-            const Atom& atomBB = *molecule.GetAtom(BB);
+            const Atom& atomBB = *molecule.GetAtomVect()[BB];
             temp += ( atomicElectronPopulation[BB] - atomBB.GetCoreCharge()  )
                      *gammaAB[indexAtomA][BB];
          }
@@ -1631,12 +1631,12 @@ void Cndo2::CalcOrbitalElectronPopulation(double** orbitalElectronPopulation,
 void Cndo2::CalcAtomicElectronPopulation(double* atomicElectronPopulation,
                                          double const* const* orbitalElectronPopulation, 
                                          const Molecule& molecule) const{
-   int totalNumberAtoms = molecule.GetNumberAtoms();
+   int totalNumberAtoms = molecule.GetAtomVect().size();
    MallocerFreer::GetInstance()->Initialize<double>(atomicElectronPopulation, totalNumberAtoms);
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
    for(int A=0; A<totalNumberAtoms; A++){
-      int firstAOIndex = molecule.GetAtom(A)->GetFirstAOIndex();
-      int numberAOs = molecule.GetAtom(A)->GetValenceSize();
+      int firstAOIndex = molecule.GetAtomVect()[A]->GetFirstAOIndex();
+      int numberAOs    = molecule.GetAtomVect()[A]->GetValenceSize();
       for(int i=firstAOIndex; i<firstAOIndex+numberAOs; i++){
          atomicElectronPopulation[A] += orbitalElectronPopulation[i][i];
       }
@@ -1646,7 +1646,7 @@ void Cndo2::CalcAtomicElectronPopulation(double* atomicElectronPopulation,
 
 // calculate gammaAB matrix. (B.56) and (B.62) in J. A. Pople book.
 void Cndo2::CalcGammaAB(double** gammaAB, const Molecule& molecule) const{
-   int totalAtomNumber = molecule.GetNumberAtoms();
+   int totalAtomNumber = molecule.GetAtomVect().size();
 
    // MPI setting of each rank
    int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
@@ -1659,14 +1659,14 @@ void Cndo2::CalcGammaAB(double** gammaAB, const Molecule& molecule) const{
    for(int A=0; A<totalAtomNumber; A++){
       int calcRank = A%mpiSize;
       if(mpiRank == calcRank){
-         const Atom& atomA = *molecule.GetAtom(A);
+         const Atom& atomA = *molecule.GetAtomVect()[A];
          int na = atomA.GetValenceShellType() + 1;
          double orbitalExponentA = atomA.GetOrbitalExponent(
                                          atomA.GetValenceShellType(), s, this->theory);
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
          for(int B=A; B<totalAtomNumber; B++){
             try{
-               const Atom& atomB = *molecule.GetAtom(B);
+               const Atom& atomB = *molecule.GetAtomVect()[B];
                int nb = atomB.GetValenceShellType() + 1;
                double orbitalExponentB = atomB.GetOrbitalExponent(
                                                atomB.GetValenceShellType(), s, this->theory);
@@ -1768,9 +1768,9 @@ void Cndo2::CalcCoreDipoleMoment(double* coreDipoleMoment,
 
    for(int i=0; i<CartesianType_end; i++){
       coreDipoleMoment[i] = 0.0;
-      for(int A=0; A<molecule.GetNumberAtoms(); A++){
-         coreDipoleMoment[i] += molecule.GetAtom(A)->GetCoreCharge()
-                               *(molecule.GetAtom(A)->GetXyz()[i] - molecule.GetXyzCOC()[i]);
+      for(int A=0; A<molecule.GetAtomVect().size(); A++){
+         coreDipoleMoment[i] += molecule.GetAtomVect()[A]->GetCoreCharge()
+                               *(molecule.GetAtomVect()[A]->GetXyz()[i] - molecule.GetXyzCOC()[i]);
       }
    }
 }
@@ -1841,7 +1841,7 @@ void Cndo2::CalcCartesianMatrixByGTOExpansion(double*** cartesianMatrix,
                                               const Molecule& molecule, 
                                               STOnGType stonG) const{
    int totalAONumber   = molecule.GetTotalNumberAOs();
-   int totalAtomNumber = molecule.GetNumberAtoms();
+   int totalAtomNumber = molecule.GetAtomVect().size();
 
    // MPI setting of each rank
    int mpiRank       = MolDS_mpi::MpiProcess::GetInstance()->GetRank();
@@ -1852,7 +1852,7 @@ void Cndo2::CalcCartesianMatrixByGTOExpansion(double*** cartesianMatrix,
    boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
 
    for(int A=0; A<totalAtomNumber; A++){
-      const Atom& atomA  = *molecule.GetAtom(A);
+      const Atom& atomA  = *molecule.GetAtomVect()[A];
       int firstAOIndexA  = atomA.GetFirstAOIndex();
       int numValenceAOsA = atomA.GetValenceSize();
       int calcRank = A%mpiSize;
@@ -1862,7 +1862,7 @@ void Cndo2::CalcCartesianMatrixByGTOExpansion(double*** cartesianMatrix,
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
             for(int B=0; B<totalAtomNumber; B++){
                try{
-                  const Atom& atomB  = *molecule.GetAtom(B);
+                  const Atom& atomB  = *molecule.GetAtomVect()[B];
                   int firstAOIndexB  = atomB.GetFirstAOIndex();
                   int numValenceAOsB = atomB.GetValenceSize();
                   for(int b=0; b<numValenceAOsB; b++){
@@ -3804,11 +3804,11 @@ void Cndo2::CalcOverlapAOsWithAnotherConfiguration(double** overlapAOs,
       ss << this->errorMessageRhs << rhsMolecule->GetTotalNumberAOs() << endl;
       throw MolDSException(ss.str());
    }
-   if(lhsMolecule.GetNumberAtoms() != rhsMolecule->GetNumberAtoms()){
+   if(lhsMolecule.GetAtomVect().size() != rhsMolecule->GetAtomVect().size()){
       stringstream ss;
       ss << this->errorMessageCalcOverlapAOsDifferentConfigurationsDiffAtoms;
-      ss << this->errorMessageLhs << lhsMolecule.GetNumberAtoms() << endl;
-      ss << this->errorMessageRhs << rhsMolecule->GetNumberAtoms() << endl;
+      ss << this->errorMessageLhs << lhsMolecule.GetAtomVect().size() << endl;
+      ss << this->errorMessageRhs << rhsMolecule->GetAtomVect().size() << endl;
       throw MolDSException(ss.str());
    }
 #ifdef MOLDS_DBG
@@ -3818,7 +3818,7 @@ void Cndo2::CalcOverlapAOsWithAnotherConfiguration(double** overlapAOs,
 #endif
    
    int totalAONumber = lhsMolecule.GetTotalNumberAOs();
-   int totalAtomNumber = lhsMolecule.GetNumberAtoms();
+   int totalAtomNumber = lhsMolecule.GetAtomVect().size();
    MallocerFreer::GetInstance()->Initialize<double>(overlapAOs, totalAONumber, totalAONumber);
 
    stringstream ompErrors;
@@ -3841,8 +3841,8 @@ void Cndo2::CalcOverlapAOsWithAnotherConfiguration(double** overlapAOs,
          bool isSymmetricOverlapAOs = false;
 #pragma omp for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
          for(int A=0; A<totalAtomNumber; A++){
-            const Atom& lhsAtom = *lhsMolecule.GetAtom(A);
-            const Atom& rhsAtom = *rhsMolecule->GetAtom(A);
+            const Atom& lhsAtom = *lhsMolecule.GetAtomVect()[A];
+            const Atom& rhsAtom = *rhsMolecule->GetAtomVect()[A];
             int firstAOIndexLhsAtom = lhsAtom.GetFirstAOIndex();
             int firstAOIndexRhsAtom = rhsAtom.GetFirstAOIndex();
             for(int i=0; i<lhsAtom.GetValenceSize(); i++){
@@ -3939,7 +3939,7 @@ void Cndo2::CalcOverlapESsWithAnotherElectronicStructure(double** overlapESs,
 // calculate OverlapAOs matrix. E.g. S_{\mu\nu} in (3.74) in J. A. Pople book.
 void Cndo2::CalcOverlapAOs(double** overlapAOs, const Molecule& molecule) const{
    int totalAONumber = molecule.GetTotalNumberAOs();
-   int totalAtomNumber = molecule.GetNumberAtoms();
+   int totalAtomNumber = molecule.GetAtomVect().size();
    MallocerFreer::GetInstance()->Initialize<double>(overlapAOs, totalAONumber, totalAONumber);
 
    // MPI setting of each rank
@@ -3951,7 +3951,7 @@ void Cndo2::CalcOverlapAOs(double** overlapAOs, const Molecule& molecule) const{
    boost::thread communicationThread( boost::bind(&MolDS_mpi::AsyncCommunicator::Run<double>, &asyncCommunicator) );
 
    for(int A=0; A<totalAtomNumber; A++){
-      const Atom& atomA = *molecule.GetAtom(A);
+      const Atom& atomA = *molecule.GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int numValenceAOs = atomA.GetValenceSize();
       int calcRank = A%mpiSize;
@@ -3985,7 +3985,7 @@ void Cndo2::CalcOverlapAOs(double** overlapAOs, const Molecule& molecule) const{
                bool symmetrize = false;
 #pragma omp for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
                for(int B=A+1; B<totalAtomNumber; B++){
-                  const Atom& atomB = *molecule.GetAtom(B);
+                  const Atom& atomB = *molecule.GetAtomVect()[B];
                   this->CalcDiatomicOverlapAOsInDiatomicFrame(diatomicOverlapAOs, atomA, atomB);
                   this->CalcRotatingMatrix(rotatingMatrix, atomA, atomB);
                   this->RotateDiatmicOverlapAOsToSpaceFrame(diatomicOverlapAOs, rotatingMatrix, tmpDiatomicOverlapAOs, tmpOldDiatomicOverlapAOs, tmpMatrixBC, tmpVectorBC);
@@ -4189,8 +4189,8 @@ void Cndo2::CalcDiatomicOverlapAOs1stDerivatives(double*** diatomicOverlapAOs1st
                                               tmpRotatedDiatomicOverlapVec,
                                               tmpMatrixBC,
                                               tmpVectorBC,
-                                              *this->molecule->GetAtom(indexAtomA),
-                                              *this->molecule->GetAtom(indexAtomB));
+                                              *this->molecule->GetAtomVect()[indexAtomA],
+                                              *this->molecule->GetAtomVect()[indexAtomB]);
 }
 
 // Second derivative of diatomic overlapAOs integrals between AOs in space fixed flame.
@@ -4330,8 +4330,8 @@ void Cndo2::CalcDiatomicOverlapAOs2ndDerivatives(double**** diatomicOverlapAOs2n
                                               tmpRotMat,
                                               tmpRotMat1stDerivs,
                                               tmpRotMat2ndDerivs,
-                                              *this->molecule->GetAtom(indexAtomA),
-                                              *this->molecule->GetAtom(indexAtomB));
+                                              *this->molecule->GetAtomVect()[indexAtomA],
+                                              *this->molecule->GetAtomVect()[indexAtomB]);
 }
 
 double Cndo2::Get2ndDerivativeElementFromDistanceDerivatives(double firstDistanceDeri,
@@ -4358,8 +4358,8 @@ double Cndo2::Get2ndDerivativeElementFromDistanceDerivatives(double firstDistanc
 void Cndo2::CalcOverlapAOsByGTOExpansion(double** overlapAOs, 
                                          const Molecule& molecule, 
                                          STOnGType stonG) const{
-   int totalAONumber = molecule.GetTotalNumberAOs();
-   int totalAtomNumber = molecule.GetNumberAtoms();
+   int totalAONumber   = molecule.GetTotalNumberAOs();
+   int totalAtomNumber = molecule.GetAtomVect().size();
 
    // calculation overlapAOs matrix
    for(int mu=0; mu<totalAONumber; mu++){
@@ -4370,10 +4370,10 @@ void Cndo2::CalcOverlapAOsByGTOExpansion(double** overlapAOs,
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE) 
    for(int A=0; A<totalAtomNumber; A++){
       try{
-         const Atom& atomA = *molecule.GetAtom(A);
+         const Atom& atomA = *molecule.GetAtomVect()[A];
          int firstAOIndexAtomA = atomA.GetFirstAOIndex();
          for(int B=A+1; B<totalAtomNumber; B++){
-            const Atom& atomB = *molecule.GetAtom(B);
+            const Atom& atomB = *molecule.GetAtomVect()[B];
             int firstAOIndexAtomB = atomB.GetFirstAOIndex();
             for(int a=0; a<atomA.GetValenceSize(); a++){
                for(int b=0; b<atomB.GetValenceSize(); b++){

--- a/src/indo/Indo.cpp
+++ b/src/indo/Indo.cpp
@@ -137,9 +137,9 @@ double Indo::GetFockDiagElement(const Atom& atomA,
       value += temp;
    
       temp = 0.0;
-      for(int B=0; B<molecule.GetNumberAtoms(); B++){
+      for(int B=0; B<molecule.GetAtomVect().size(); B++){
          if(B != indexAtomA){
-            const Atom& atomB = *molecule.GetAtom(B);
+            const Atom& atomB = *molecule.GetAtomVect()[B];
             temp += ( atomicElectronPopulation[B] - atomB.GetCoreCharge()  )
                      *gammaAB[indexAtomA][B];
          }
@@ -207,8 +207,8 @@ double Indo::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
    value = Cndo2::GetMolecularIntegralElement(moI, moJ, moK, moL, molecule, fockMatrix, gammaAB);
 
    // Aditional terms for INDO, see Eq. (10) in [RZ_1973]
-   for(int A=0; A<molecule.GetNumberAtoms(); A++){
-      const Atom& atomA = *molecule.GetAtom(A);
+   for(int A=0; A<molecule.GetAtomVect().size(); A++){
+      const Atom& atomA = *molecule.GetAtomVect()[A];
       firstAOIndexA = atomA.GetFirstAOIndex();
       numberAOsA = atomA.GetValenceSize();
 

--- a/src/mc/MC.cpp
+++ b/src/mc/MC.cpp
@@ -161,9 +161,9 @@ void MC::CreateTrialConfiguration(Molecule* trial,
                                   > (*realRand),
                                   double stepWidth) const{
    // disturb an atom in trial molecule
-   int movedAtomIndex = static_cast<int>((*realRand)()*this->molecule->GetNumberAtoms());
-   const Atom& reffAtom = *current.GetAtom(movedAtomIndex);
-   Atom* trialAtom = trial->GetAtom(movedAtomIndex);
+   int movedAtomIndex = static_cast<int>((*realRand)()*this->molecule->GetAtomVect().size());
+   const Atom& reffAtom = *current.GetAtomVect()[movedAtomIndex];
+   Atom* trialAtom = trial->GetAtomVect()[movedAtomIndex];
    double dr[CartesianType_end] = {0.0, 0.0, 0.0};
    for(int i=0; i<CartesianType_end; i++){
       dr[i] = stepWidth*(2.0*(*realRand)() -1.0);
@@ -177,8 +177,8 @@ void MC::CreateTrialConfiguration(Molecule* trial,
    for(int i=0; i<CartesianType_end; i++){
       coreCenterShift[i] = dr[i]*trialAtomCoreMass/totalCoreMass;
    }
-   for(int a=0; a<current.GetNumberAtoms(); a++){
-      Atom* trialAtom = trial->GetAtom(a);
+   for(int a=0; a<current.GetAtomVect().size(); a++){
+      Atom* trialAtom = trial->GetAtomVect()[a];
       for(int i=0; i<CartesianType_end; i++){
          trialAtom->GetXyz()[i] -= coreCenterShift[i];
       }

--- a/src/md/MD.cpp
+++ b/src/md/MD.cpp
@@ -138,8 +138,8 @@ void MD::DoMD(){
 
 void MD::UpdateMomenta(const Molecule& molecule, double const* const* matrixForce, double dt) const{
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      Atom* atom = molecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      Atom* atom = molecule.GetAtomVect()[a];
       for(int i=0; i<CartesianType_end; i++){
          atom->GetPxyz()[i] += 0.5*dt*(matrixForce[a][i]);
       }
@@ -148,8 +148,8 @@ void MD::UpdateMomenta(const Molecule& molecule, double const* const* matrixForc
 
 void MD::UpdateCoordinates(Molecule& molecule, double dt) const{
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
-      for(int a=0; a<molecule.GetNumberAtoms(); a++){
-         Atom* atom = molecule.GetAtom(a);
+      for(int a=0; a<molecule.GetAtomVect().size(); a++){
+         Atom* atom = molecule.GetAtomVect()[a];
          double coreMass = atom->GetCoreMass();
          for(int i=0; i<CartesianType_end; i++){
             atom->GetXyz()[i] += dt*atom->GetPxyz()[i]/coreMass;
@@ -183,8 +183,8 @@ double MD::OutputEnergies(boost::shared_ptr<ElectronicStructure> electronicStruc
    int elecState = Parameters::GetInstance()->GetElectronicStateIndexMD();
    double eV2AU  = Parameters::GetInstance()->GetEV2AU();
    double coreKineticEnergy = 0.0;
-   for(int a=0; a<this->molecule->GetNumberAtoms(); a++){
-      Atom* atom = this->molecule->GetAtom(a);
+   for(int a=0; a<this->molecule->GetAtomVect().size(); a++){
+      Atom* atom = this->molecule->GetAtomVect()[a];
       double coreMass = atom->GetCoreMass();
       for(int i=0; i<CartesianType_end; i++){
          coreKineticEnergy += 0.5*pow(atom->GetPxyz()[i],2.0)/coreMass;

--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -78,70 +78,70 @@ Mndo::Mndo() : MolDS_zindo::ZindoS(){
 Mndo::~Mndo(){
    OrbitalType twoElecLimit = dxy;
    MallocerFreer::GetInstance()->Free<double>(&this->twoElecsTwoAtomCores, 
-                                              this->molecule->GetNumberAtoms(),
-                                              this->molecule->GetNumberAtoms(),
+                                              this->molecule->GetAtomVect().size(),
+                                              this->molecule->GetAtomVect().size(),
                                               twoElecLimit,
                                               twoElecLimit,
                                               twoElecLimit,
                                               twoElecLimit);
    MallocerFreer::GetInstance()->Free<double>(&this->twoElecsAtomEpcCores, 
-                                              this->molecule->GetNumberAtoms(),
-                                              this->molecule->GetNumberEpcs(),
+                                              this->molecule->GetAtomVect().size(),
+                                              this->molecule->GetEpcVect().size(),
                                               twoElecLimit,
                                               twoElecLimit,
                                               twoElecLimit,
                                               twoElecLimit);
    int numBuff = (twoElecLimit+1)*twoElecLimit/2;
    MallocerFreer::GetInstance()->Free<double>(&this->twoElecsTwoAtomCoresMpiBuff, 
-                                              this->molecule->GetNumberAtoms(),
-                                              this->molecule->GetNumberAtoms(),
+                                              this->molecule->GetAtomVect().size(),
+                                              this->molecule->GetAtomVect().size(),
                                               numBuff,
                                               numBuff);
    MallocerFreer::GetInstance()->Free<double>(&this->twoElecsAtomEpcCoresMpiBuff, 
-                                              this->molecule->GetNumberAtoms(),
-                                              this->molecule->GetNumberEpcs(),
+                                              this->molecule->GetAtomVect().size(),
+                                              this->molecule->GetEpcVect().size(),
                                               numBuff,
                                               numBuff);
    MallocerFreer::GetInstance()->Free<double>(&this->normalForceConstants,
-                                              CartesianType_end*molecule->GetNumberAtoms());
+                                              CartesianType_end*molecule->GetAtomVect().size());
    MallocerFreer::GetInstance()->Free<double>(&this->normalModes,
-                                              CartesianType_end*molecule->GetNumberAtoms(),
-                                              CartesianType_end*molecule->GetNumberAtoms());
+                                              CartesianType_end*molecule->GetAtomVect().size(),
+                                              CartesianType_end*molecule->GetAtomVect().size());
 }
 
 void Mndo::SetMolecule(Molecule* molecule){
    ZindoS::SetMolecule(molecule);
    OrbitalType twoElecLimit = dxy;
    MallocerFreer::GetInstance()->Malloc<double>(&this->twoElecsTwoAtomCores,
-                                                molecule->GetNumberAtoms(),
-                                                molecule->GetNumberAtoms(),
+                                                molecule->GetAtomVect().size(),
+                                                molecule->GetAtomVect().size(),
                                                 twoElecLimit,
                                                 twoElecLimit,
                                                 twoElecLimit,
                                                 twoElecLimit);
    MallocerFreer::GetInstance()->Malloc<double>(&this->twoElecsAtomEpcCores,
-                                                molecule->GetNumberAtoms(),
-                                                molecule->GetNumberEpcs(),
+                                                molecule->GetAtomVect().size(),
+                                                molecule->GetEpcVect().size(),
                                                 twoElecLimit,
                                                 twoElecLimit,
                                                 twoElecLimit,
                                                 twoElecLimit);
    int numBuff = (twoElecLimit+1)*twoElecLimit/2;
    MallocerFreer::GetInstance()->Malloc<double>(&this->twoElecsTwoAtomCoresMpiBuff, 
-                                                this->molecule->GetNumberAtoms(),
-                                                this->molecule->GetNumberAtoms(),
+                                                this->molecule->GetAtomVect().size(),
+                                                this->molecule->GetAtomVect().size(),
                                                 numBuff,
                                                 numBuff);
    MallocerFreer::GetInstance()->Malloc<double>(&this->twoElecsAtomEpcCoresMpiBuff, 
-                                                this->molecule->GetNumberAtoms(),
-                                                this->molecule->GetNumberEpcs(),
+                                                this->molecule->GetAtomVect().size(),
+                                                this->molecule->GetEpcVect().size(),
                                                 numBuff,
                                                 numBuff);
    MallocerFreer::GetInstance()->Malloc<double>(&this->normalForceConstants,
-                                                CartesianType_end*molecule->GetNumberAtoms());
+                                                CartesianType_end*molecule->GetAtomVect().size());
    MallocerFreer::GetInstance()->Malloc<double>(&this->normalModes,
-                                                CartesianType_end*molecule->GetNumberAtoms(),
-                                                CartesianType_end*molecule->GetNumberAtoms());
+                                                CartesianType_end*molecule->GetAtomVect().size(),
+                                                CartesianType_end*molecule->GetAtomVect().size());
 }
 void Mndo::SetMessages(){
    this->errorMessageSCFNotConverged 
@@ -315,8 +315,8 @@ double Mndo::GetAuxiliaryDiatomCoreRepulsionEnergy2ndDerivative(const Atom& atom
 }
 
 double Mndo::GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double tmp = this->GetAuxiliaryDiatomCoreRepulsionEnergy(atomA,
                                                             atomB,
                                                             this->molecule->GetDistanceAtoms(atomA, atomB));
@@ -327,8 +327,8 @@ double Mndo::GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const{
 }
 
 double Mndo::GetAtomCoreEpcCoulombEnergy(int indexAtom, int indexEpc) const{
-   const Atom& atom = *this->molecule->GetAtom(indexAtom);
-   const Atom& epc  = *this->molecule->GetAtom(indexEpc);
+   const Atom& atom = *this->molecule->GetAtomVect()[indexAtom];
+   const Atom& epc  = *this->molecule->GetAtomVect()[indexEpc];
    double distance = this->molecule->GetDistanceAtomEpc(indexAtom, indexEpc);
    return atom.GetCoreCharge()*epc.GetCoreCharge()/distance; 
 }
@@ -339,8 +339,8 @@ double Mndo::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
                                                  int indexAtomB, 
                                                  CartesianType axisA) const{
    double value =0.0;
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double distanceAB = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
    double twoElecInt = this->twoElecsTwoAtomCores[indexAtomA][indexAtomB][s][s][s][s];
    double twoElecInt1stDeriv = this->GetNddoRepulsionIntegral1stDerivative(
@@ -359,8 +359,8 @@ double Mndo::GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
                                                     CartesianType axisA1,
                                                     CartesianType axisA2) const{
    double value =0.0;
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double distanceAB = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
    double twoElecInt = this->twoElecsTwoAtomCores[indexAtomA][indexAtomB][s][s][s][s];
    double twoElecInt1stDeriv1 = this->GetNddoRepulsionIntegral1stDerivative(atomA, s, s, 
@@ -403,8 +403,8 @@ void Mndo::CalcHeatsFormation(double* heatsFormation,
                               const Molecule& molecule) const{
    int groundState = 0;
    *heatsFormation = this->GetElectronicEnergy(groundState);
-   for(int A=0; A<molecule.GetNumberAtoms(); A++){
-      const Atom& atom = *molecule.GetAtom(A);
+   for(int A=0; A<molecule.GetAtomVect().size(); A++){
+      const Atom& atom = *molecule.GetAtomVect()[A];
       *heatsFormation -= atom.GetMndoElecEnergyAtom();
       *heatsFormation += atom.GetMndoHeatsFormAtom();
    }
@@ -420,7 +420,7 @@ void Mndo::CalcNormalModes(double** normalModes, double* normalForceConstants, c
    bool isMassWeighted = true;
    this->CalcHessianSCF(normalModes, isMassWeighted);
    bool calcEigenVectors = true;
-   int hessianDim = CartesianType_end*molecule.GetNumberAtoms();
+   int hessianDim = CartesianType_end*molecule.GetAtomVect().size();
    MolDS_wrappers::Lapack::GetInstance()->Dsyevd(normalModes,
                                                  normalForceConstants,
                                                  hessianDim,
@@ -464,10 +464,10 @@ double Mndo::GetFockDiagElement(const Atom& atomA,
       value += temp;
 
       temp = 0.0;
-      int totalNumberAtoms=molecule.GetNumberAtoms();
+      int totalNumberAtoms=molecule.GetAtomVect().size();
       for(int B=0; B<totalNumberAtoms; B++){
          if(B != indexAtomA){
-            const Atom& atomB = *molecule.GetAtom(B);
+            const Atom& atomB = *molecule.GetAtomVect()[B];
             int firstAOIndexB = atomB.GetFirstAOIndex();
             int valenceSizeB = atomB.GetValenceSize();
             for(int lambda=0; lambda<valenceSizeB; lambda++){
@@ -492,11 +492,11 @@ double Mndo::GetFockDiagElement(const Atom& atomA,
       value += temp;
       
       // coulomb repulsion with point charge *
-      int numEpcs = molecule.GetNumberEpcs();
+      int numEpcs = molecule.GetEpcVect().size();
       if(0<numEpcs){
          double elecCharge = -1.0;
          for(int i=0; i<numEpcs; i++){
-            double epcCharge = molecule.GetEpc(i)->GetCoreCharge();
+            double epcCharge = molecule.GetEpcVect()[i]->GetCoreCharge();
             value += elecCharge*epcCharge*twoElecsAtomEpcCores[indexAtomA][i][mu][mu][s][s];
          }
       }
@@ -538,10 +538,10 @@ double Mndo::GetFockOffDiagElement(const Atom& atomA,
          exchange = this->GetExchangeInt(orbitalMu, orbitalNu, atomA); 
          temp = (1.5*exchange - 0.5*coulomb)
                *orbitalElectronPopulation[mu+firstAOIndexA][nu+firstAOIndexB];
-         int totalNumberAtoms = molecule.GetNumberAtoms();
+         int totalNumberAtoms = molecule.GetAtomVect().size();
          for(int BB=0; BB<totalNumberAtoms; BB++){
             if(BB != indexAtomA){
-               const Atom& atomBB = *molecule.GetAtom(BB);
+               const Atom& atomBB = *molecule.GetAtomVect()[BB];
                int firstAOIndexBB = atomBB.GetFirstAOIndex();
                int valenceSizeBB = atomBB.GetValenceSize();
                for(int lambda=0; lambda<valenceSizeBB; lambda++){
@@ -564,11 +564,11 @@ double Mndo::GetFockOffDiagElement(const Atom& atomA,
             }
          }
          // coulomb repulsion with point charge *
-         int numEpcs = molecule.GetNumberEpcs();
+         int numEpcs = molecule.GetEpcVect().size();
          if(0<numEpcs){
             double elecCharge = -1.0;
             for(int i=0; i<numEpcs; i++){
-               double epcCharge = molecule.GetEpc(i)->GetCoreCharge();
+               double epcCharge = molecule.GetEpcVect()[i]->GetCoreCharge();
                value += elecCharge*epcCharge*twoElecsAtomEpcCores[indexAtomA][i][mu][nu][s][s];
             }
          }
@@ -663,7 +663,7 @@ double Mndo::GetElectronCoreAttraction(int indexAtomA,
                                        int mu, 
                                        int nu, 
                                        double const* const* const* const* const* const* twoElecsTwoAtomCores) const{
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    return -1.0*atomB.GetCoreCharge()*twoElecsTwoAtomCores[indexAtomA][indexAtomB][mu][nu][s][s];
 }
 
@@ -677,7 +677,7 @@ double Mndo::GetElectronCoreAttraction1stDerivative(int indexAtomA,
                                                     int nu, 
                                                     double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivatives,
                                                     CartesianType axisA) const{
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double value = -1.0*atomB.GetCoreCharge()
                   *diatomicTwoElecsTwoCores1stDerivatives[mu][nu][s][s][axisA];
    return value;
@@ -711,13 +711,13 @@ double Mndo::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
                                          double const* const* fockMatrix, 
                                          double const* const* gammaAB) const{
    double value = 0.0;
-   for(int A=0; A<molecule.GetNumberAtoms(); A++){
-      const Atom& atomA = *molecule.GetAtom(A);
+   for(int A=0; A<molecule.GetAtomVect().size(); A++){
+      const Atom& atomA = *molecule.GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
 
-      for(int B=A; B<molecule.GetNumberAtoms(); B++){
-         const Atom& atomB = *molecule.GetAtom(B);
+      for(int B=A; B<molecule.GetAtomVect().size(); B++){
+         const Atom& atomB = *molecule.GetAtomVect()[B];
          int firstAOIndexB = atomB.GetFirstAOIndex();
          int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -834,13 +834,13 @@ void Mndo::CalcCISMatrix(double** matrixCIS) const{
              
             // Fast algorith, but this is not easy to read. 
             // Slow algorithm is alos written below.
-            for(int A=0; A<molecule->GetNumberAtoms(); A++){
-               const Atom& atomA = *molecule->GetAtom(A);
+            for(int A=0; A<molecule->GetAtomVect().size(); A++){
+               const Atom& atomA = *molecule->GetAtomVect()[A];
                int firstAOIndexA = atomA.GetFirstAOIndex();
                int lastAOIndexA  = atomA.GetLastAOIndex();
 
-               for(int B=A; B<molecule->GetNumberAtoms(); B++){
-                  const Atom& atomB = *molecule->GetAtom(B);
+               for(int B=A; B<molecule->GetAtomVect().size(); B++){
+                  const Atom& atomB = *molecule->GetAtomVect()[B];
                   int firstAOIndexB = atomB.GetFirstAOIndex();
                   int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -1073,25 +1073,25 @@ void Mndo::MallocTempMatricesEachThreadCalcHessianSCF(double*****    diatomicOve
                                                       double***      tmpMatrixBC,
                                                       double**       tmpVectorBC) const{
    MallocerFreer::GetInstance()->Malloc<double>(diatomicOverlapAOs1stDerivs,
-                                                this->molecule->GetNumberAtoms(),
+                                                this->molecule->GetAtomVect().size(),
                                                 OrbitalType_end,
                                                 OrbitalType_end,
                                                 CartesianType_end);
    MallocerFreer::GetInstance()->Malloc<double>(diatomicOverlapAOs2ndDerivs,
-                                                this->molecule->GetNumberAtoms(),
+                                                this->molecule->GetAtomVect().size(),
                                                 OrbitalType_end,
                                                 OrbitalType_end,
                                                 CartesianType_end,
                                                 CartesianType_end);
    MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecsTwoCores1stDerivs,
-                                                this->molecule->GetNumberAtoms(),
+                                                this->molecule->GetAtomVect().size(),
                                                 dxy,
                                                 dxy,
                                                 dxy,
                                                 dxy,
                                                 CartesianType_end);
    MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecsTwoCores2ndDerivs,
-                                                this->molecule->GetNumberAtoms(),
+                                                this->molecule->GetAtomVect().size(),
                                                 dxy,
                                                 dxy,
                                                 dxy,
@@ -1173,25 +1173,25 @@ void Mndo::FreeTempMatricesEachThreadCalcHessianSCF(double*****    diatomicOverl
                                                     double***      tmpMatrixBC,
                                                     double**       tmpVectorBC) const{
    MallocerFreer::GetInstance()->Free<double>(diatomicOverlapAOs1stDerivs,
-                                              this->molecule->GetNumberAtoms(),
+                                              this->molecule->GetAtomVect().size(),
                                               OrbitalType_end,
                                               OrbitalType_end,
                                               CartesianType_end);
    MallocerFreer::GetInstance()->Free<double>(diatomicOverlapAOs2ndDerivs,
-                                              this->molecule->GetNumberAtoms(),
+                                              this->molecule->GetAtomVect().size(),
                                               OrbitalType_end,
                                               OrbitalType_end,
                                               CartesianType_end,
                                               CartesianType_end);
    MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecsTwoCores1stDerivs,
-                                              this->molecule->GetNumberAtoms(),
+                                              this->molecule->GetAtomVect().size(),
                                               dxy,
                                               dxy,
                                               dxy,
                                               dxy,
                                               CartesianType_end);
    MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecsTwoCores2ndDerivs,
-                                              this->molecule->GetNumberAtoms(),
+                                              this->molecule->GetAtomVect().size(),
                                               dxy,
                                               dxy,
                                               dxy,
@@ -1265,8 +1265,8 @@ double Mndo::GetAuxiliaryHessianElement1(int mu,
                                          CartesianType axisA2,
                                          double const* const* orbitalElectronPopulation,
                                          double const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomC = *this->molecule->GetAtomVect()[indexAtomC];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    double value = orbitalElectronPopulation[mu]
                                            [nu]
@@ -1291,8 +1291,8 @@ double Mndo::GetAuxiliaryHessianElement2(int mu,
                                          CartesianType axisB,
                                          double const* const* const* const* orbitalElectronPopulation1stDerivs,
                                          double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomC = *this->molecule->GetAtomVect()[indexAtomC];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    double value = orbitalElectronPopulation1stDerivs[mu]
                                                     [nu]
@@ -1317,8 +1317,8 @@ double Mndo::GetAuxiliaryHessianElement3(int lambda,
                                          CartesianType axisA2,
                                          double const* const* orbitalElectronPopulation,
                                          double const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomC = *this->molecule->GetAtomVect()[indexAtomC];
    int firstAOIndexC = atomC.GetFirstAOIndex();
    double value = orbitalElectronPopulation[lambda]
                                            [sigma]
@@ -1343,8 +1343,8 @@ double Mndo::GetAuxiliaryHessianElement4(int lambda,
                                          CartesianType axisB,
                                          double const* const* const* const* orbitalElectronPopulation1stDerivs,
                                          double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomC = *this->molecule->GetAtomVect()[indexAtomC];
    int firstAOIndexC = atomC.GetFirstAOIndex();
    double value = orbitalElectronPopulation1stDerivs[lambda]
                                                     [sigma]
@@ -1369,8 +1369,8 @@ double Mndo::GetAuxiliaryHessianElement5(int mu,
                                          CartesianType axisA2,
                                          double const* const* orbitalElectronPopulation,
                                          double const* const* const* const* diatomicOverlapAOs2ndDerivs) const{
-   const Atom& atomA = *molecule->GetAtom(indexAtomA);
-   const Atom& atomC = *molecule->GetAtom(indexAtomC);
+   const Atom& atomA = *molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomC = *molecule->GetAtomVect()[indexAtomC];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexC = atomC.GetFirstAOIndex();
    double bondParameterA = atomA.GetBondingParameter(this->theory, atomA.GetValence(mu-firstAOIndexA));
@@ -1397,8 +1397,8 @@ double Mndo::GetAuxiliaryHessianElement6(int mu,
                                          CartesianType axisB,
                                          double const* const* const* const* orbitalElectronPopulation1stDerivs,
                                          double const* const* const* diatomicOverlapAOs1stDerivs) const{
-   const Atom& atomA = *molecule->GetAtom(indexAtomA);
-   const Atom& atomC = *molecule->GetAtom(indexAtomC);
+   const Atom& atomA = *molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomC = *molecule->GetAtomVect()[indexAtomC];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexC = atomC.GetFirstAOIndex();
    double bondParameterA = atomA.GetBondingParameter(this->theory, atomA.GetValence(mu-firstAOIndexA));
@@ -1428,8 +1428,8 @@ double Mndo::GetAuxiliaryHessianElement7(int mu,
                                          CartesianType axisA2,
                                          double const* const* orbitalElectronPopulation,
                                          double const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
-   const Atom& atomA = *molecule->GetAtom(indexAtomA);
-   const Atom& atomC = *molecule->GetAtom(indexAtomC);
+   const Atom& atomA = *molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomC = *molecule->GetAtomVect()[indexAtomC];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexC = atomC.GetFirstAOIndex();
    double temp1 = orbitalElectronPopulation[mu][nu]*orbitalElectronPopulation[lambda][sigma];
@@ -1459,8 +1459,8 @@ double Mndo::GetAuxiliaryHessianElement8(int mu,
                                          double const* const* orbitalElectronPopulation,
                                          double const* const* const* const* orbitalElectronPopulation1stDerivs,
                                          double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *molecule->GetAtom(indexAtomA);
-   const Atom& atomC = *molecule->GetAtom(indexAtomC);
+   const Atom& atomA = *molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomC = *molecule->GetAtomVect()[indexAtomC];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexC = atomC.GetFirstAOIndex();
    double temp1 = orbitalElectronPopulation1stDerivs[mu][nu]       [indexAtomB][axisB]
@@ -1494,12 +1494,12 @@ double Mndo::GetHessianElementSameAtomsSCF(int indexAtomA,
                                            double const* const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
    double value=0.0;
    int indexAtomB = indexAtomA;
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
-   for(int indexAtomC=0; indexAtomC<this->molecule->GetNumberAtoms(); indexAtomC++){
+   for(int indexAtomC=0; indexAtomC<this->molecule->GetAtomVect().size(); indexAtomC++){
       if(indexAtomA != indexAtomC){
-         const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
+         const Atom& atomC = *this->molecule->GetAtomVect()[indexAtomC];
          int firstAOIndexC = atomC.GetFirstAOIndex();
          int numberAOsC = atomC.GetValenceSize();
 
@@ -1630,8 +1630,8 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
                                                 double const* const* const* const* const* const*        diatomicTwoElecsTwoCores1stDerivs,
                                                 double const* const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
    double value=0.0;
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexB = atomB.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
@@ -1693,9 +1693,9 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
       }
    }
 
-   for(int indexAtomC=0; indexAtomC<this->molecule->GetNumberAtoms(); indexAtomC++){
+   for(int indexAtomC=0; indexAtomC<this->molecule->GetAtomVect().size(); indexAtomC++){
       if(indexAtomA != indexAtomC){
-         const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
+         const Atom& atomC = *this->molecule->GetAtomVect()[indexAtomC];
          int firstAOIndexC = atomC.GetFirstAOIndex();
          int numberAOsC = atomC.GetValenceSize();
 
@@ -1786,7 +1786,7 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
       MallocerFreer::GetInstance()->Malloc<double>(&orbitalElectronPopulation1stDerivs, 
                                                    totalNumberAOs,
                                                    totalNumberAOs,
-                                                   this->molecule->GetNumberAtoms(),
+                                                   this->molecule->GetAtomVect().size(),
                                                    CartesianType_end);
       this->CalcOrbitalElectronPopulation1stDerivatives(orbitalElectronPopulation1stDerivs);
 
@@ -1834,14 +1834,14 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
                                                              &tmpMatrixBC,
                                                              &tmpVectorBC);
 #pragma omp for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)                                                 
-            for(int indexAtomA=0; indexAtomA<this->molecule->GetNumberAtoms(); indexAtomA++){
-               const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
+            for(int indexAtomA=0; indexAtomA<this->molecule->GetAtomVect().size(); indexAtomA++){
+               const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
                int firstAOIndexA = atomA.GetFirstAOIndex();
                int lastAOIndexA  = atomA.GetLastAOIndex();
                for(int axisA = XAxis; axisA<CartesianType_end; axisA++){
             
                   // calculation of derivatives of the overlapAOss and two electron integrals
-                  for(int indexAtomB=0; indexAtomB<this->molecule->GetNumberAtoms(); indexAtomB++){
+                  for(int indexAtomB=0; indexAtomB<this->molecule->GetAtomVect().size(); indexAtomB++){
                      if(indexAtomA != indexAtomB){
                         this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs[indexAtomB], 
                                                                    tmpDiaOverlapAOsInDiaFrame,        
@@ -1885,10 +1885,10 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
             
                   // calculation of each hessian element
                   int k = indexAtomA*CartesianType_end + axisA; // hessian index, i.e. hessian[k][l]
-                  for(int indexAtomB=indexAtomA; indexAtomB<this->molecule->GetNumberAtoms(); indexAtomB++){
+                  for(int indexAtomB=indexAtomA; indexAtomB<this->molecule->GetAtomVect().size(); indexAtomB++){
                      // hessian element (atomA != atomB)
                      if(indexAtomA!=indexAtomB){
-                        const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+                        const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
                         for(int axisB = XAxis; axisB<CartesianType_end; axisB++){
                            int l = indexAtomB*CartesianType_end + axisB; // hessian index, i.e. hessian[k][l]
                            hessianSCF[k][l] = this->GetHessianElementDifferentAtomsSCF(indexAtomA, 
@@ -1963,23 +1963,23 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
       MallocerFreer::GetInstance()->Free<double>(&orbitalElectronPopulation1stDerivs, 
                                                  totalNumberAOs,
                                                  totalNumberAOs,
-                                                 this->molecule->GetNumberAtoms(),
+                                                 this->molecule->GetAtomVect().size(),
                                                  CartesianType_end);
       throw ex;
    }
    MallocerFreer::GetInstance()->Free<double>(&orbitalElectronPopulation1stDerivs, 
                                               totalNumberAOs,
                                               totalNumberAOs,
-                                              this->molecule->GetNumberAtoms(),
+                                              this->molecule->GetAtomVect().size(),
                                               CartesianType_end);
-   int hessianDim = this->molecule->GetNumberAtoms()*CartesianType_end;
+   int hessianDim = this->molecule->GetAtomVect().size()*CartesianType_end;
    for(int k=0; k<hessianDim; k++){
       for(int l=k; l<hessianDim; l++){
          hessianSCF[l][k] = hessianSCF[k][l];
       }
    }
    /*
-   int hessianDim = this->molecule->GetNumberAtoms()*CartesianType_end;
+   int hessianDim = this->molecule->GetAtomVect().size()*CartesianType_end;
    for(int i=0; i<hessianDim; i++){
       for(int j=0; j<hessianDim; j++){
          printf("hess elem: %d %d %e\n",i,j,hessianSCF[i][j]);
@@ -1999,7 +1999,7 @@ void Mndo::CalcOrbitalElectronPopulation1stDerivatives(double**** orbitalElectro
    vector<MoIndexPair> redundantQIndeces;
    this->CalcActiveSetVariablesQ(&nonRedundantQIndeces, &redundantQIndeces, numberOcc, numberVir);
    int dimensionCPHF = nonRedundantQIndeces.size() + redundantQIndeces.size();
-   int numberCPHFs = this->molecule->GetNumberAtoms()*CartesianType_end;
+   int numberCPHFs = this->molecule->GetAtomVect().size()*CartesianType_end;
    double** solutionsCPHF = NULL; // solutions of CPHF
    double** transposedFockMatrix = NULL; // transposed Fock matrix
    try{
@@ -2012,7 +2012,7 @@ void Mndo::CalcOrbitalElectronPopulation1stDerivatives(double**** orbitalElectro
       for(int mu=0; mu<totalNumberAOs; mu++){
          try{
             for(int nu=0; nu<totalNumberAOs; nu++){
-               for(int indexAtomA=0; indexAtomA<this->molecule->GetNumberAtoms(); indexAtomA++){
+               for(int indexAtomA=0; indexAtomA<this->molecule->GetAtomVect().size(); indexAtomA++){
                   for(int axis=XAxis; axis<CartesianType_end; axis++){
          
                      int moI, moJ;
@@ -2047,7 +2047,7 @@ void Mndo::CalcOrbitalElectronPopulation1stDerivatives(double**** orbitalElectro
 
       /*
       // check the CPHF's solutions 
-      for(int indexAtomA=0; indexAtomA<this->molecule->GetNumberAtoms(); indexAtomA++){
+      for(int indexAtomA=0; indexAtomA<this->molecule->GetAtomVect().size(); indexAtomA++){
          for(int axis=XAxis; axis<CartesianType_end; axis++){
             double temp=0.0;
             printf("cphf: atom:%d axis:%s start\n ",indexAtomA,CartesianTypeStr(axis));
@@ -2077,7 +2077,7 @@ void Mndo::SolveCPHF(double** solutionsCPHF,
                      const vector<MoIndexPair>& nonRedundantQIndeces,
                      const vector<MoIndexPair>& redundantQIndeces) const{
    int dimensionCPHF = nonRedundantQIndeces.size() + redundantQIndeces.size();
-   int numberCPHFs = this->molecule->GetNumberAtoms()*CartesianType_end;
+   int numberCPHFs = this->molecule->GetAtomVect().size()*CartesianType_end;
    double** matrixCPHF = NULL; // (Gmamma - K matrix)N, see (40) - (46) to slove (34) in [PT_1996].
    try{
       this->MallocTempMatricesSolveCPHF(&matrixCPHF, dimensionCPHF);
@@ -2100,7 +2100,7 @@ void Mndo::CalcStaticFirstOrderFocks(double** staticFirstOrderFocks,
                                      const vector<MoIndexPair>& redundantQIndeces) const{
    stringstream ompErrors;
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
-   for(int indexAtomA=0; indexAtomA<this->molecule->GetNumberAtoms(); indexAtomA++){
+   for(int indexAtomA=0; indexAtomA<this->molecule->GetAtomVect().size(); indexAtomA++){
       try{
          for(int axisA=XAxis; axisA<CartesianType_end; axisA++){
             int k=indexAtomA*CartesianType_end + axisA;
@@ -2158,13 +2158,13 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
       MallocerFreer::GetInstance()->Malloc<double>(&tmpRotatedDiatomicOverlapVec,       OrbitalType_end*OrbitalType_end);
       MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrixBC,                        OrbitalType_end, OrbitalType_end);
       MallocerFreer::GetInstance()->Malloc<double>(&tmpVectorBC,                        OrbitalType_end*OrbitalType_end);
-      const Atom& atomA = *molecule->GetAtom(indexAtomA);
+      const Atom& atomA = *molecule->GetAtomVect()[indexAtomA];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
       int coreChargeA   = atomA.GetCoreCharge();
-      for(int indexAtomB=0; indexAtomB<this->molecule->GetNumberAtoms(); indexAtomB++){
+      for(int indexAtomB=0; indexAtomB<this->molecule->GetAtomVect().size(); indexAtomB++){
          if(indexAtomA != indexAtomB){
-            const Atom& atomB = *molecule->GetAtom(indexAtomB);
+            const Atom& atomB = *molecule->GetAtomVect()[indexAtomB];
             int firstAOIndexB = atomB.GetFirstAOIndex();
             int lastAOIndexB  = atomB.GetLastAOIndex();
             int coreChargeB   = atomB.GetCoreCharge();
@@ -2468,7 +2468,7 @@ void Mndo::CalcForceSCFElecCoreAttractionPart(double* force,
                                              int indexAtomA, 
                                              int indexAtomB,
                                              double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
    for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
@@ -2491,8 +2491,8 @@ void Mndo::CalcForceSCFOverlapAOsPart(double* force,
                                      int indexAtomA, 
                                      int indexAtomB,
                                      double const* const* const* diatomicOverlapAOs1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexB = atomB.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
@@ -2520,8 +2520,8 @@ void Mndo::CalcForceSCFTwoElecPart(double* force,
                                   int indexAtomA, 
                                   int indexAtomB,
                                   double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexB = atomB.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
@@ -2559,8 +2559,8 @@ void Mndo::CalcForceExcitedStaticPart(double* force,
                                       int indexAtomA, 
                                       int indexAtomB,
                                       double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexB = atomB.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
@@ -2592,7 +2592,7 @@ void Mndo::CalcForceExcitedElecCoreAttractionPart(double* force,
                                                   int indexAtomA, 
                                                   int indexAtomB,
                                                   double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
    for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
@@ -2616,8 +2616,8 @@ void Mndo::CalcForceExcitedTwoElecPart(double* force,
                                        int indexAtomA, 
                                        int indexAtomB,
                                        double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexB = atomB.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
@@ -2661,9 +2661,9 @@ void Mndo::CalcForce(const vector<int>& elecStates){
    }
 
    // this loop is MPI-parallelized
-   for(int a=0; a<this->molecule->GetNumberAtoms(); a++){
+   for(int a=0; a<this->molecule->GetAtomVect().size(); a++){
       if(a%mpiSize != mpiRank){continue;}
-      const Atom& atomA = *molecule->GetAtom(a);
+      const Atom& atomA = *molecule->GetAtomVect()[a];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
       stringstream ompErrors;
@@ -2697,9 +2697,9 @@ void Mndo::CalcForce(const vector<int>& elecStates){
                                               &tmpDiatomicTwoElecsTwoCores);
 
 #pragma omp for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
-            for(int b=0; b<this->molecule->GetNumberAtoms(); b++){
+            for(int b=0; b<this->molecule->GetAtomVect().size(); b++){
                if(a == b){continue;}
-               const Atom& atomB = *molecule->GetAtom(b);
+               const Atom& atomB = *molecule->GetAtomVect()[b];
                int firstAOIndexB = atomB.GetFirstAOIndex();
                int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -2848,7 +2848,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
    }// end of for(int a) with MPI parallelization
 
    // communication to reduce thsi->matrixForce on all node (namely, all_reduce)
-   int numTransported = elecStates.size()*this->molecule->GetNumberAtoms()*CartesianType_end;
+   int numTransported = elecStates.size()*this->molecule->GetAtomVect().size()*CartesianType_end;
    MolDS_mpi::MpiProcess::GetInstance()->AllReduce(&this->matrixForce[0][0][0], numTransported, std::plus<double>());
 }
 
@@ -2972,13 +2972,13 @@ double Mndo::GetSmallQElement(int moI,
    int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
    bool isMoPOcc = moP<numberOcc ? true : false;
    
-   for(int A=0; A<molecule->GetNumberAtoms(); A++){
-      const Atom& atomA = *molecule->GetAtom(A);
+   for(int A=0; A<molecule->GetAtomVect().size(); A++){
+      const Atom& atomA = *molecule->GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
 
-      for(int B=A; B<molecule->GetNumberAtoms(); B++){
-         const Atom& atomB = *molecule->GetAtom(B);
+      for(int B=A; B<molecule->GetAtomVect().size(); B++){
+         const Atom& atomB = *molecule->GetAtomVect()[B];
          int firstAOIndexB = atomB.GetFirstAOIndex();
          int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -3153,8 +3153,8 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
 
    // Fast algorith, but this is not easy to read. 
    // Slow algorithm is alos written below.
-   for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
-      const Atom& atomA = *this->molecule->GetAtom(A);
+   for(int A=0; A<this->molecule->GetAtomVect().size(); A++){
+      const Atom& atomA = *this->molecule->GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
 
@@ -3197,8 +3197,8 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
                         *this->fockMatrix[moK][mu];
             }
 
-            for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
-               const Atom& atomB = *this->molecule->GetAtom(B);
+            for(int B=A; B<this->molecule->GetAtomVect().size(); B++){
+               const Atom& atomB = *this->molecule->GetAtomVect()[B];
                int firstAOIndexB = atomB.GetFirstAOIndex();
                int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -3301,16 +3301,16 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
    double*  twiceMoJK = NULL;
    double*  tmpVector = NULL;
    int numAOs = this->molecule->GetTotalNumberAOs();
-   MallocerFreer::GetInstance()->Malloc<double>(&twoElec,   this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
-   for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
-      const Atom& atomA = *this->molecule->GetAtom(A);
+   MallocerFreer::GetInstance()->Malloc<double>(&twoElec,   this->molecule->GetAtomVect().size()*dxy*dxy, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoKL, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJL, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoJK, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&tmpVector, this->molecule->GetAtomVect().size()*dxy*dxy);
+   for(int A=0; A<this->molecule->GetAtomVect().size(); A++){
+      const Atom& atomA = *this->molecule->GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
       for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
@@ -3322,8 +3322,8 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
       }
    }
 
-   for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
-      const Atom& atomB = *this->molecule->GetAtom(B);
+   for(int B=0; B<this->molecule->GetAtomVect().size(); B++){
+      const Atom& atomB = *this->molecule->GetAtomVect()[B];
       int firstAOIndexB = atomB.GetFirstAOIndex();
       int lastAOIndexB  = atomB.GetLastAOIndex();
       for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
@@ -3335,12 +3335,12 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
       }
    }
 
-   for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
-      const Atom& atomA = *this->molecule->GetAtom(A);
+   for(int A=0; A<this->molecule->GetAtomVect().size(); A++){
+      const Atom& atomA = *this->molecule->GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
-      for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
-         const Atom& atomB = *this->molecule->GetAtom(B);
+      for(int B=A; B<this->molecule->GetAtomVect().size(); B++){
+         const Atom& atomB = *this->molecule->GetAtomVect()[B];
          int firstAOIndexB = atomB.GetFirstAOIndex();
          int lastAOIndexB  = atomB.GetLastAOIndex();
          double gamma = 0.0;
@@ -3389,29 +3389,29 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
          }
       }
    }
-   MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy, 
+   MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetAtomVect().size()*dxy*dxy, 
                                               twoElec, 
                                               twiceMoKL,
                                               tmpVector);
-   value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, tmpVector);
-   MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy, 
+   value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetAtomVect().size()*dxy*dxy,twiceMoIJ, tmpVector);
+   MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetAtomVect().size()*dxy*dxy, 
                                               twoElec, 
                                               twiceMoJL,
                                               tmpVector);
-   value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, tmpVector);
-   MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetNumberAtoms()*dxy*dxy, 
+   value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetAtomVect().size()*dxy*dxy,twiceMoIK, tmpVector);
+   MolDS_wrappers::Blas::GetInstance()->Dsymv(this->molecule->GetAtomVect().size()*dxy*dxy, 
                                               twoElec, 
                                               twiceMoJK,
                                               tmpVector);
-   value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, tmpVector);
-   MallocerFreer::GetInstance()->Free<double>(&twoElec,   this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoKL, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoJL, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoJK, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&tmpVector, this->molecule->GetNumberAtoms()*dxy*dxy);
+   value -= MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetAtomVect().size()*dxy*dxy,twiceMoIL, tmpVector);
+   MallocerFreer::GetInstance()->Free<double>(&twoElec,   this->molecule->GetAtomVect().size()*dxy*dxy, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoKL, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoJL, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoJK, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&tmpVector, this->molecule->GetAtomVect().size()*dxy*dxy);
    // End of algorithm using blas
    */
 
@@ -3425,14 +3425,14 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
    double** twiceMoB  = NULL;
    double** tmpMatrix = NULL;
    int numAOs = this->molecule->GetTotalNumberAOs();
-   MallocerFreer::GetInstance()->Malloc<double>(&twoElec,   this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
-   for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
-      const Atom& atomA = *this->molecule->GetAtom(A);
+   MallocerFreer::GetInstance()->Malloc<double>(&twoElec,   this->molecule->GetAtomVect().size()*dxy*dxy, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIJ, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIK, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoIL, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&twiceMoB, 3, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrix,3, this->molecule->GetAtomVect().size()*dxy*dxy);
+   for(int A=0; A<this->molecule->GetAtomVect().size(); A++){
+      const Atom& atomA = *this->molecule->GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
       for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
@@ -3444,8 +3444,8 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
       }
    }
 
-   for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
-      const Atom& atomB = *this->molecule->GetAtom(B);
+   for(int B=0; B<this->molecule->GetAtomVect().size(); B++){
+      const Atom& atomB = *this->molecule->GetAtomVect()[B];
       int firstAOIndexB = atomB.GetFirstAOIndex();
       int lastAOIndexB  = atomB.GetLastAOIndex();
       for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
@@ -3457,12 +3457,12 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
       }
    }
 
-   for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
-      const Atom& atomA = *this->molecule->GetAtom(A);
+   for(int A=0; A<this->molecule->GetAtomVect().size(); A++){
+      const Atom& atomA = *this->molecule->GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
-      for(int B=0; B<this->molecule->GetNumberAtoms(); B++){
-         const Atom& atomB = *this->molecule->GetAtom(B);
+      for(int B=0; B<this->molecule->GetAtomVect().size(); B++){
+         const Atom& atomB = *this->molecule->GetAtomVect()[B];
          int firstAOIndexB = atomB.GetFirstAOIndex();
          int lastAOIndexB  = atomB.GetLastAOIndex();
          double gamma = 0.0;
@@ -3513,23 +3513,23 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
    }
 
    MolDS_wrappers::Blas::GetInstance()->Dgemm(false, true, true,
-                                              this->molecule->GetNumberAtoms()*dxy*dxy,
+                                              this->molecule->GetAtomVect().size()*dxy*dxy,
                                               3,
-                                              this->molecule->GetNumberAtoms()*dxy*dxy,
+                                              this->molecule->GetAtomVect().size()*dxy*dxy,
                                               1.0,
                                               twoElec,
                                               twiceMoB,
                                               0.0,
                                               tmpMatrix);
-   value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIJ, &tmpMatrix[0][0]);
-   value -=    MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIK, &tmpMatrix[1][0]);
-   value -=    MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetNumberAtoms()*dxy*dxy,twiceMoIL, &tmpMatrix[2][0]);
-   MallocerFreer::GetInstance()->Free<double>(&twoElec,   this->molecule->GetNumberAtoms()*dxy*dxy, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&twiceMoB, 3, this->molecule->GetNumberAtoms()*dxy*dxy);
-   MallocerFreer::GetInstance()->Free<double>(&tmpMatrix,3, this->molecule->GetNumberAtoms()*dxy*dxy);
+   value = 4.0*MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetAtomVect().size()*dxy*dxy,twiceMoIJ, &tmpMatrix[0][0]);
+   value -=    MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetAtomVect().size()*dxy*dxy,twiceMoIK, &tmpMatrix[1][0]);
+   value -=    MolDS_wrappers::Blas::GetInstance()->Ddot(this->molecule->GetAtomVect().size()*dxy*dxy,twiceMoIL, &tmpMatrix[2][0]);
+   MallocerFreer::GetInstance()->Free<double>(&twoElec,   this->molecule->GetAtomVect().size()*dxy*dxy, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoIJ, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoIK, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoIL, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&twiceMoB, 3, this->molecule->GetAtomVect().size()*dxy*dxy);
+   MallocerFreer::GetInstance()->Free<double>(&tmpMatrix,3, this->molecule->GetAtomVect().size()*dxy*dxy);
    // End of second algorithm using blas
    */
 
@@ -3562,7 +3562,7 @@ void Mndo::CalcTwoElecsTwoAtomCores(double****** twoElecsTwoAtomCores,
       throw MolDSException(this->errorMessageCalcTwoElecsTwoAtomCoresNullMatrix);
    }
 #endif
-   int totalNumberAtoms = molecule.GetNumberAtoms();
+   int totalNumberAtoms = molecule.GetAtomVect().size();
    MallocerFreer::GetInstance()->Initialize<double>(twoElecsTwoAtomCores, 
                                                     totalNumberAtoms,
                                                     totalNumberAtoms,
@@ -3674,14 +3674,14 @@ void Mndo::CalcTwoElecsTwoAtomCores(double****** twoElecsTwoAtomCores,
 
 void Mndo::CalcTwoElecsAtomEpcCores(double****** twoElecsAtomEpcCores, 
                                     const Molecule& molecule) const{
-   if(molecule.GetNumberEpcs()<=0){return;}
+   if(molecule.GetEpcVect().empty()){return;}
 #ifdef MOLDS_DBG
    if(twoElecsAtomEpcCores == NULL){
       throw MolDSException(this->errorMessageCalcTwoElecsAtomEpcCoresNullMatrix);
    }
 #endif
-   int totalNumberAtoms = molecule.GetNumberAtoms();
-   int totalNumberEpcs  = molecule.GetNumberEpcs();
+   int totalNumberAtoms = molecule.GetAtomVect().size();
+   int totalNumberEpcs  = molecule.GetEpcVect().size();
    MallocerFreer::GetInstance()->Initialize<double>(twoElecsAtomEpcCores, 
                                                     totalNumberAtoms,
                                                     totalNumberEpcs,
@@ -3705,7 +3705,7 @@ void Mndo::CalcTwoElecsAtomEpcCores(double****** twoElecsAtomEpcCores,
             double**   tmpRotMat                   = NULL;
             double**   tmpMatrixBC                 = NULL;
             double*    tmpVectorBC                 = NULL;
-            const Atom& atom = *molecule.GetAtom(a);
+            const Atom& atom = *molecule.GetAtomVect()[a];
             try{
                MallocerFreer::GetInstance()->Malloc<double>(&diatomicTwoElecsTwoCores,    dxy, dxy, dxy, dxy);
                MallocerFreer::GetInstance()->Malloc<double>(&tmpDiatomicTwoElecsTwoCores, dxy*dxy*dxy*dxy);
@@ -3715,7 +3715,7 @@ void Mndo::CalcTwoElecsAtomEpcCores(double****** twoElecsAtomEpcCores,
                // note that terms with condition a==b are not needed to calculate. 
 //#pragma omp for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
                for(int b=0; b<totalNumberEpcs; b++){
-                  const Atom& epc  = *molecule.GetEpc(b);
+                  const Atom& epc  = *molecule.GetEpcVect()[b];
                   this->CalcDiatomicTwoElecsTwoCores(diatomicTwoElecsTwoCores, 
                                                      tmpDiatomicTwoElecsTwoCores,
                                                      tmpRotMat, 
@@ -3883,8 +3883,8 @@ void Mndo::CalcDiatomicTwoElecsTwoCores(double**** matrix,
                                         double*    tmpVectorBC,
                                         int indexAtomA, 
                                         int indexAtomB) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    this->CalcDiatomicTwoElecsTwoCores(matrix, 
                                       tmpVec,
                                       tmpRotMat, 
@@ -3908,8 +3908,8 @@ void Mndo::CalcDiatomicTwoElecsTwoCores1stDerivatives(double***** matrix,
                                                       double****  tmpDiatomicTwoElecsTwoCores,
                                                       int indexAtomA, 
                                                       int indexAtomB) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    if(indexAtomA == indexAtomB){
       stringstream ss;
       ss << this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesSameAtoms;
@@ -3992,8 +3992,8 @@ void Mndo::CalcDiatomicTwoElecsTwoCores2ndDerivatives(double****** matrix,
                                                       double*****  tmpDiatomicTwoElecsTwoCores1stDerivs,
                                                       int indexAtomA, 
                                                       int indexAtomB) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    if(indexAtomA == indexAtomB){
       stringstream ss;
       ss << this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesSameAtoms;

--- a/src/nasco/NASCO.cpp
+++ b/src/nasco/NASCO.cpp
@@ -177,8 +177,8 @@ void NASCO::DoNASCO(Molecule& molecule){
 void NASCO::UpdateMomenta(Molecule& molecule, 
                           double const* const* matrixForce, 
                           const double dt) const{
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      Atom* atom = molecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      Atom* atom = molecule.GetAtomVect()[a];
       for(int i=0; i<CartesianType_end; i++){
          atom->GetPxyz()[i] += 0.5*dt*(matrixForce[a][i]);
       }
@@ -188,9 +188,9 @@ void NASCO::UpdateMomenta(Molecule& molecule,
 void NASCO::UpdateCoordinates(Molecule& tmpMolecule, 
                               const Molecule& molecule,
                               const double dt) const{
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      Atom* atom    = molecule.GetAtom(a);
-      Atom* tmpAtom = tmpMolecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      Atom* atom    = molecule.GetAtomVect()[a];
+      Atom* tmpAtom = tmpMolecule.GetAtomVect()[a];
       double coreMass = tmpAtom->GetCoreMass();
       for(int i=0; i<CartesianType_end; i++){
          tmpAtom->GetXyz()[i] =  atom->GetXyz()[i] + dt*atom->GetPxyz()[i]/coreMass;
@@ -253,8 +253,8 @@ double NASCO::OutputEnergies(const ElectronicStructure& electronicStructure,
                              int elecState) const{
    double eV2AU = Parameters::GetInstance()->GetEV2AU();
    double coreKineticEnergy = 0.0;
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      Atom* atom = molecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      Atom* atom = molecule.GetAtomVect()[a];
       double coreMass = atom->GetCoreMass();
       for(int i=0; i<CartesianType_end; i++){
          coreKineticEnergy += 0.5*pow(atom->GetPxyz()[i],2.0)/coreMass;

--- a/src/optimization/BFGS.cpp
+++ b/src/optimization/BFGS.cpp
@@ -114,7 +114,7 @@ void BFGS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStruct
    double lineSearchInitialEnergy   = 0.0;
    double const* const* matrixForce = NULL;
    double const* vectorForce        = NULL;
-   const int dimension = molecule.GetNumberAtoms()*CartesianType_end;
+   const int dimension = molecule.GetAtomVect().size()*CartesianType_end;
    double** matrixHessian        = NULL;
    double*  vectorOldForce       = NULL;
    double*  vectorStep           = NULL;
@@ -156,7 +156,7 @@ void BFGS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStruct
          trustRadius=min(trustRadius,maxNormStep);
 
          //Calculate RFO step
-         MallocerFreer::GetInstance()->Malloc(&matrixStep, molecule.GetNumberAtoms(), CartesianType_end);
+         MallocerFreer::GetInstance()->Malloc(&matrixStep, molecule.GetAtomVect().size(), CartesianType_end);
          vectorStep = &matrixStep[0][0];
          this->CalcRFOStep(vectorStep, matrixHessian, vectorForce, trustRadius, dimension);
 
@@ -215,16 +215,16 @@ void BFGS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStruct
    catch(MolDSException ex){
       MallocerFreer::GetInstance()->Free(&matrixHessian, dimension, dimension);
       MallocerFreer::GetInstance()->Free(&vectorOldForce, dimension);
-      MallocerFreer::GetInstance()->Free(&matrixStep, molecule.GetNumberAtoms(), CartesianType_end);
-      MallocerFreer::GetInstance()->Free(&matrixDisplacement, molecule.GetNumberAtoms(), CartesianType_end);
-      MallocerFreer::GetInstance()->Free(&matrixOldCoordinates, molecule.GetNumberAtoms(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free(&matrixStep, molecule.GetAtomVect().size(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free(&matrixDisplacement, molecule.GetAtomVect().size(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free(&matrixOldCoordinates, molecule.GetAtomVect().size(), CartesianType_end);
       throw ex;
    }
    MallocerFreer::GetInstance()->Free(&matrixHessian, dimension, dimension);
    MallocerFreer::GetInstance()->Free(&vectorOldForce, dimension);
-   MallocerFreer::GetInstance()->Free(&matrixStep, molecule.GetNumberAtoms(), CartesianType_end);
-   MallocerFreer::GetInstance()->Free(&matrixDisplacement, molecule.GetNumberAtoms(), CartesianType_end);
-   MallocerFreer::GetInstance()->Free(&matrixOldCoordinates, molecule.GetNumberAtoms(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&matrixStep, molecule.GetAtomVect().size(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&matrixDisplacement, molecule.GetAtomVect().size(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&matrixOldCoordinates, molecule.GetAtomVect().size(), CartesianType_end);
 }
 
 void BFGS::CalcRFOStep(double* vectorStep,
@@ -356,7 +356,7 @@ void BFGS::ShiftHessianRedundantMode(double** matrixHessian,
                                      const Molecule& molecule) const{
    const double one                      = 1;
    const double largeEigenvalue          = 1.0e3;
-   const int    numAtoms                 = molecule.GetNumberAtoms();
+   const int    numAtoms                 = molecule.GetAtomVect().size();
    const int    dimension                = numAtoms *CartesianType_end;
    const int    numTranslationalModes    = 3;
    const int    numRotationalModes       = 3;
@@ -397,7 +397,7 @@ void BFGS::ShiftHessianRedundantMode(double** matrixHessian,
       for(int c=0; c<numRotationalModes;c++){
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
          for(int n=0;n<numAtoms;n++){
-            const double* xyz = molecule.GetAtom(n)->GetXyz();
+            const double* xyz = molecule.GetAtomVect()[n]->GetXyz();
             for(int d=0;d<CartesianType_end;d++){
                matrixesRedundantModes[c+numTranslationalModes][n][d] = 0.0;
                for(int e=0;e<CartesianType_end;e++){
@@ -548,8 +548,8 @@ void BFGS::RollbackMolecularGeometry(MolDS_base::Molecule& molecule,
    bool tempCanOutputLogs = molecule.CanOutputLogs();
    bool rollbackCanOutputLogs = true;
    molecule.SetCanOutputLogs(rollbackCanOutputLogs);
-   for(int i=0;i<molecule.GetNumberAtoms();i++){
-      const Atom* atom = molecule.GetAtom(i);
+   for(int i=0;i<molecule.GetAtomVect().size();i++){
+      const Atom* atom = molecule.GetAtomVect()[i];
       double*     xyz  = atom->GetXyz();
       for(int j=0;j<CartesianType_end;j++){
          xyz[j] = matrixOldCoordinates[i][j];
@@ -562,9 +562,9 @@ void BFGS::CalcDisplacement(double      *      *& matrixDisplacement,
                             double const* const*  matrixOldCoordinates,
                             const MolDS_base::Molecule& molecule)const{
    //Calculate displacement (K_k at Eq. (15) in [SJTO_1983])
-   MallocerFreer::GetInstance()->Malloc(&matrixDisplacement, molecule.GetNumberAtoms(), CartesianType_end);
-   for(int i=0;i<molecule.GetNumberAtoms();i++){
-      const Atom*   atom = molecule.GetAtom(i);
+   MallocerFreer::GetInstance()->Malloc(&matrixDisplacement, molecule.GetAtomVect().size(), CartesianType_end);
+   for(int i=0;i<molecule.GetAtomVect().size();i++){
+      const Atom*   atom = molecule.GetAtomVect()[i];
       const double* xyz  = atom->GetXyz();
       for(int j=0;j<CartesianType_end;j++){
          matrixDisplacement[i][j] = xyz[j] - matrixOldCoordinates[i][j];
@@ -575,9 +575,9 @@ void BFGS::CalcDisplacement(double      *      *& matrixDisplacement,
 void BFGS::StoreMolecularGeometry(double **& matrixCoordinates, 
                                   const MolDS_base::Molecule& molecule)const{
    //Store old coordinates
-   MallocerFreer::GetInstance()->Malloc(&matrixCoordinates, molecule.GetNumberAtoms(), CartesianType_end);
-   for(int i=0;i<molecule.GetNumberAtoms();i++){
-      const Atom*   atom = molecule.GetAtom(i);
+   MallocerFreer::GetInstance()->Malloc(&matrixCoordinates, molecule.GetAtomVect().size(), CartesianType_end);
+   for(int i=0;i<molecule.GetAtomVect().size();i++){
+      const Atom*   atom = molecule.GetAtomVect()[i];
       const double* xyz  = atom->GetXyz();
       for(int j=0;j<CartesianType_end;j++){
          matrixCoordinates[i][j] = xyz[j];

--- a/src/optimization/ConjugateGradient.cpp
+++ b/src/optimization/ConjugateGradient.cpp
@@ -88,9 +88,9 @@ void ConjugateGradient::SearchMinimum(boost::shared_ptr<ElectronicStructure> ele
    requireGuess = false;
    matrixForce = electronicStructure->GetForce(elecState);
    try{
-      MallocerFreer::GetInstance()->Malloc<double>(&oldMatrixForce, molecule.GetNumberAtoms(), CartesianType_end);
-      MallocerFreer::GetInstance()->Malloc<double>(&matrixSearchDirection, molecule.GetNumberAtoms(), CartesianType_end);
-      for(int a=0;a<molecule.GetNumberAtoms();a++){
+      MallocerFreer::GetInstance()->Malloc<double>(&oldMatrixForce, molecule.GetAtomVect().size(), CartesianType_end);
+      MallocerFreer::GetInstance()->Malloc<double>(&matrixSearchDirection, molecule.GetAtomVect().size(), CartesianType_end);
+      for(int a=0;a<molecule.GetAtomVect().size();a++){
          for(int i=0; i<CartesianType_end; i++){
             matrixSearchDirection[a][i] = matrixForce[a][i];
          }
@@ -120,12 +120,12 @@ void ConjugateGradient::SearchMinimum(boost::shared_ptr<ElectronicStructure> ele
       }
    }
    catch(MolDSException ex){
-      MallocerFreer::GetInstance()->Free<double>(&oldMatrixForce, molecule.GetNumberAtoms(), CartesianType_end);
-      MallocerFreer::GetInstance()->Free<double>(&matrixSearchDirection, molecule.GetNumberAtoms(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free<double>(&oldMatrixForce, molecule.GetAtomVect().size(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free<double>(&matrixSearchDirection, molecule.GetAtomVect().size(), CartesianType_end);
       throw ex;
    }
-   MallocerFreer::GetInstance()->Free<double>(&oldMatrixForce, molecule.GetNumberAtoms(), CartesianType_end);
-   MallocerFreer::GetInstance()->Free<double>(&matrixSearchDirection, molecule.GetNumberAtoms(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free<double>(&oldMatrixForce, molecule.GetAtomVect().size(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free<double>(&matrixSearchDirection, molecule.GetAtomVect().size(), CartesianType_end);
    *lineSearchedEnergy = lineSearchCurrentEnergy;
 }
 
@@ -135,7 +135,7 @@ void ConjugateGradient::UpdateSearchDirection(double const* const** matrixForce,
                                               boost::shared_ptr<ElectronicStructure> electronicStructure, 
                                               const MolDS_base::Molecule& molecule,
                                               int elecState) const{
-   for(int a=0;a<molecule.GetNumberAtoms();a++){
+   for(int a=0;a<molecule.GetAtomVect().size();a++){
       for(int i=0; i<CartesianType_end; i++){
          oldMatrixForce[a][i] = (*matrixForce)[a][i];
       }
@@ -143,14 +143,14 @@ void ConjugateGradient::UpdateSearchDirection(double const* const** matrixForce,
    *matrixForce = electronicStructure->GetForce(elecState);
    double beta=0.0;
    double temp=0.0;
-   for(int a=0;a<molecule.GetNumberAtoms();a++){
+   for(int a=0;a<molecule.GetAtomVect().size();a++){
       for(int i=0; i<CartesianType_end; i++){
          temp += pow(oldMatrixForce[a][i],2.0);
          beta += ((*matrixForce)[a][i] - oldMatrixForce[a][i])*(*matrixForce)[a][i];
       }
    }
    beta /= temp;
-   for(int a=0;a<molecule.GetNumberAtoms();a++){
+   for(int a=0;a<molecule.GetAtomVect().size();a++){
       for(int i=0; i<CartesianType_end; i++){
          matrixSearchDirection[a][i] *= beta;
          matrixSearchDirection[a][i] += (*matrixForce)[a][i];

--- a/src/optimization/GEDIIS.cpp
+++ b/src/optimization/GEDIIS.cpp
@@ -89,7 +89,7 @@ void GEDIIS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStru
    double lineSearchInitialEnergy   = 0.0;
    double const* const* matrixForce = NULL;
    double const* vectorForce        = NULL;
-   const int dimension = molecule.GetNumberAtoms()*CartesianType_end;
+   const int dimension = molecule.GetAtomVect().size()*CartesianType_end;
    double** matrixHessian           = NULL;
    double*  vectorOldForce          = NULL;
    double*  vectorStep              = NULL;
@@ -137,8 +137,8 @@ void GEDIIS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStru
          lineSearchInitialEnergy = lineSearchCurrentEnergy;
          double preRFOEnergy = lineSearchInitialEnergy;
 
-         MallocerFreer::GetInstance()->Malloc(&matrixGEDIISCoordinates, molecule.GetNumberAtoms(), CartesianType_end);
-         MallocerFreer::GetInstance()->Malloc(&matrixGEDIISForce,       molecule.GetNumberAtoms(), CartesianType_end);
+         MallocerFreer::GetInstance()->Malloc(&matrixGEDIISCoordinates, molecule.GetAtomVect().size(), CartesianType_end);
+         MallocerFreer::GetInstance()->Malloc(&matrixGEDIISForce,       molecule.GetAtomVect().size(), CartesianType_end);
          try{
             history.SolveGEDIISEquation(&preRFOEnergy, matrixGEDIISCoordinates, matrixGEDIISForce);
 
@@ -176,7 +176,7 @@ void GEDIIS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStru
          this->ShiftHessianRedundantMode(matrixHessian, molecule);
 
          //Calculate RFO step
-         MallocerFreer::GetInstance()->Malloc(&matrixStep, molecule.GetNumberAtoms(), CartesianType_end);
+         MallocerFreer::GetInstance()->Malloc(&matrixStep, molecule.GetAtomVect().size(), CartesianType_end);
          vectorStep = &matrixStep[0][0];
          this->CalcRFOStep(vectorStep, matrixHessian, vectorForce, trustRadius, dimension);
 
@@ -237,20 +237,20 @@ void GEDIIS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStru
    catch(MolDSException ex){
       MallocerFreer::GetInstance()->Free(&matrixHessian, dimension, dimension);
       MallocerFreer::GetInstance()->Free(&vectorOldForce, dimension);
-      MallocerFreer::GetInstance()->Free(&matrixStep             , molecule.GetNumberAtoms(), CartesianType_end);
-      MallocerFreer::GetInstance()->Free(&matrixDisplacement     , molecule.GetNumberAtoms(), CartesianType_end);
-      MallocerFreer::GetInstance()->Free(&matrixOldCoordinates   , molecule.GetNumberAtoms(), CartesianType_end);
-      MallocerFreer::GetInstance()->Free(&matrixGEDIISCoordinates, molecule.GetNumberAtoms(), CartesianType_end);
-      MallocerFreer::GetInstance()->Free(&matrixGEDIISForce      , molecule.GetNumberAtoms(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free(&matrixStep             , molecule.GetAtomVect().size(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free(&matrixDisplacement     , molecule.GetAtomVect().size(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free(&matrixOldCoordinates   , molecule.GetAtomVect().size(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free(&matrixGEDIISCoordinates, molecule.GetAtomVect().size(), CartesianType_end);
+      MallocerFreer::GetInstance()->Free(&matrixGEDIISForce      , molecule.GetAtomVect().size(), CartesianType_end);
       throw ex;
    }
    MallocerFreer::GetInstance()->Free(&matrixHessian, dimension, dimension);
    MallocerFreer::GetInstance()->Free(&vectorOldForce, dimension);
-   MallocerFreer::GetInstance()->Free(&matrixStep             , molecule.GetNumberAtoms(), CartesianType_end);
-   MallocerFreer::GetInstance()->Free(&matrixDisplacement     , molecule.GetNumberAtoms(), CartesianType_end);
-   MallocerFreer::GetInstance()->Free(&matrixOldCoordinates   , molecule.GetNumberAtoms(), CartesianType_end);
-   MallocerFreer::GetInstance()->Free(&matrixGEDIISCoordinates, molecule.GetNumberAtoms(), CartesianType_end);
-   MallocerFreer::GetInstance()->Free(&matrixGEDIISForce      , molecule.GetNumberAtoms(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&matrixStep             , molecule.GetAtomVect().size(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&matrixDisplacement     , molecule.GetAtomVect().size(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&matrixOldCoordinates   , molecule.GetAtomVect().size(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&matrixGEDIISCoordinates, molecule.GetAtomVect().size(), CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&matrixGEDIISForce      , molecule.GetAtomVect().size(), CartesianType_end);
 }
 
 GEDIIS::GEDIISHistory::GEDIISHistory():maxEntryCount(5){
@@ -287,12 +287,12 @@ void GEDIIS::GEDIISHistory::DiscardEntries(){
 GEDIIS::GEDIISHistory::Entry::Entry(double energy,
                                     const MolDS_base::Molecule& molecule,
                                     double const* const* matrixForce):
-   energy(energy),numAtoms(molecule.GetNumberAtoms()),matrixCoordinate(NULL),matrixForce(NULL) {
+   energy(energy),numAtoms(molecule.GetAtomVect().size()),matrixCoordinate(NULL),matrixForce(NULL) {
    MallocerFreer::GetInstance()->Malloc(&this->matrixCoordinate, this->numAtoms, CartesianType_end);
    MallocerFreer::GetInstance()->Malloc(&this->matrixForce,      this->numAtoms, CartesianType_end);
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
    for(int i = 0; i < this->numAtoms; i++){
-      const Atom*   atom = molecule.GetAtom(i);
+      const Atom*   atom = molecule.GetAtomVect()[i];
       const double* xyz  = atom->GetXyz();
       for(int j = 0; j < CartesianType_end; j++){
          this->matrixCoordinate[i][j] = xyz[j];

--- a/src/optimization/Optimizer.cpp
+++ b/src/optimization/Optimizer.cpp
@@ -131,16 +131,16 @@ void Optimizer::CheckEnableTheoryType(TheoryType theoryType) const{
 
 void Optimizer::ClearMolecularMomenta(Molecule& molecule) const{
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE) 
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      const Atom* atom = molecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      const Atom* atom = molecule.GetAtomVect()[a];
       atom->SetPxyz(0.0, 0.0, 0.0);
    }
 }
 
 void Optimizer::UpdateMolecularCoordinates(Molecule& molecule, double const* const* matrixForce, double dt) const{
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE) 
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      const Atom* atom = molecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      const Atom* atom = molecule.GetAtomVect()[a];
       double coreMass = atom->GetCoreMass();
       for(int i=0; i<CartesianType_end; i++){
          atom->GetXyz()[i] += dt*matrixForce[a][i]/coreMass;
@@ -151,8 +151,8 @@ void Optimizer::UpdateMolecularCoordinates(Molecule& molecule, double const* con
 
 void Optimizer::UpdateMolecularCoordinates(Molecule& molecule, double const* const* matrixForce) const{
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
-      const Atom* atom = molecule.GetAtom(a);
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
+      const Atom* atom = molecule.GetAtomVect()[a];
       for(int i=0; i<CartesianType_end; i++){
          atom->GetXyz()[i] += matrixForce[a][i];
       }
@@ -233,7 +233,7 @@ bool Optimizer::SatisfiesConvergenceCriterion(double const* const* matrixForce,
    double maxGradient = 0.0;
    double sumSqureGradient = 0.0;
    double energyDifference = currentEnergy - oldEnergy;
-   for(int a=0; a<molecule.GetNumberAtoms(); a++){
+   for(int a=0; a<molecule.GetAtomVect().size(); a++){
       for(int i=0; i<CartesianType_end; i++){
          if(maxGradient<fabs(matrixForce[a][i])){
             maxGradient = fabs(matrixForce[a][i]);
@@ -241,7 +241,7 @@ bool Optimizer::SatisfiesConvergenceCriterion(double const* const* matrixForce,
          sumSqureGradient += pow(matrixForce[a][i],2.0);
       }
    }
-   sumSqureGradient /= static_cast<double>(molecule.GetNumberAtoms()*CartesianType_end);
+   sumSqureGradient /= static_cast<double>(molecule.GetAtomVect().size()*CartesianType_end);
    double rmsGradient = sqrt(sumSqureGradient);
 
    // output logs

--- a/src/pm3/Pm3Pddg.cpp
+++ b/src/pm3/Pm3Pddg.cpp
@@ -143,8 +143,8 @@ double Pm3Pddg::GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) con
    double pm3Term = Pm3::GetDiatomCoreRepulsionEnergy(indexAtomA, indexAtomB);
 
    // pddg additional term, eq. (4) in [RCJ_2002]
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int na = atomA.GetNumberValenceElectrons();
    int nb = atomB.GetNumberValenceElectrons();
    double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
@@ -172,8 +172,8 @@ double Pm3Pddg::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
    double pm3Term = Pm3::GetDiatomCoreRepulsion1stDerivative(indexAtomA, indexAtomB, axisA);
 
    // pddg additional term, first derivative of eq. (4) in [RCJ_2002]
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int na = atomA.GetNumberValenceElectrons();
    int nb = atomB.GetNumberValenceElectrons();
    double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
@@ -203,8 +203,8 @@ double Pm3Pddg::GetDiatomCoreRepulsion2ndDerivative(int indexAtomA,
    double pm3Term = Pm3::GetDiatomCoreRepulsion2ndDerivative(indexAtomA, indexAtomB, axisA1, axisA2);
 
    // pddg additional term, first derivative of eq. (4) in [RCJ_2002]
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    double distance = this->molecule->GetDistanceAtoms(indexAtomA, indexAtomB);
    double dCartesian1 = (atomA.GetXyz()[axisA1] - atomB.GetXyz()[axisA1]);
    double dCartesian2 = (atomA.GetXyz()[axisA2] - atomB.GetXyz()[axisA2]);

--- a/src/rpmd/RPMD.cpp
+++ b/src/rpmd/RPMD.cpp
@@ -94,16 +94,16 @@ void RPMD::UpdateMomenta(const vector<boost::shared_ptr<Molecule> >& molecularBe
                          double temperature){
    double kB = Parameters::GetInstance()->GetBoltzmann();
    int numBeads = molecularBeads.size();
-   int numAtom  = molecularBeads[0]->GetNumberAtoms();
+   int numAtom  = molecularBeads[0]->GetAtomVect().size();
    for(int b=0; b<numBeads; b++){
       int preB  = b==0 ? numBeads-1 : b-1;
       int postB = b==numBeads-1 ? 0 : b+1;
       double const* const* electronicForceMatrix 
          = electronicStructureBeads[b]->GetForce(elecState);;
       for(int a=0; a<numAtom; a++){
-         Atom* atom      = molecularBeads[b]->GetAtom(a);
-         Atom* preAtom   = molecularBeads[preB]->GetAtom(a);
-         Atom* postAtom  = molecularBeads[postB]->GetAtom(a);
+         Atom* atom      = molecularBeads[b]->GetAtomVect()[a];
+         Atom* preAtom   = molecularBeads[preB]->GetAtomVect()[a];
+         Atom* postAtom  = molecularBeads[postB]->GetAtomVect()[a];
          double coreMass = atom->GetCoreMass();
          for(int i=0; i<CartesianType_end; i++){
             double beadsForce = -1.0*coreMass*pow(kB*temperature*static_cast<double>(numBeads),2.0)
@@ -118,10 +118,10 @@ void RPMD::UpdateMomenta(const vector<boost::shared_ptr<Molecule> >& molecularBe
 void RPMD::UpdateCoordinates(const vector<boost::shared_ptr<Molecule> >& molecularBeads,
                              double dt){
    int numBeads = molecularBeads.size();
-   int numAtom = molecularBeads[0]->GetNumberAtoms();
+   int numAtom = molecularBeads[0]->GetAtomVect().size();
    for(int b=0; b<numBeads; b++){
       for(int a=0; a<numAtom; a++){
-         Atom* atom = molecularBeads[b]->GetAtom(a);
+         Atom* atom = molecularBeads[b]->GetAtomVect()[a];
          double coreMass = atom->GetCoreMass();
          for(int i=0; i<CartesianType_end; i++){
             atom->GetXyz()[i] += dt*atom->GetPxyz()[i]/coreMass;
@@ -150,7 +150,7 @@ void RPMD::FluctuateBeads(const vector<boost::shared_ptr<Molecule> >& molecularB
       molecule->SetCanOutputLogs(false);
       mc->SetMolecule(molecule);
       mc->SetCanOutputLogs(false);
-      mc->DoMC(molecule->GetNumberAtoms(), elecState, temperature, stepWidth, seed+b);
+      mc->DoMC(molecule->GetAtomVect().size(), elecState, temperature, stepWidth, seed+b);
    }
 }
 
@@ -168,7 +168,7 @@ void RPMD::DoRPMD(const Molecule& refferenceMolecule){
    double dt          = Parameters::GetInstance()->GetTimeWidthRPMD();
    double kB          = Parameters::GetInstance()->GetBoltzmann();
    int numBeads       = Parameters::GetInstance()->GetNumberBeadsRPMD();
-   int numAtom        = refferenceMolecule.GetNumberAtoms();
+   int numAtom        = refferenceMolecule.GetAtomVect().size();
 
    // create Beads
    vector<boost::shared_ptr<Molecule> > molecularBeads;
@@ -249,12 +249,12 @@ double RPMD::OutputEnergies(const vector<boost::shared_ptr<Molecule> >& molecula
                             int elecState,
                             double temperature){
    int numBeads = molecularBeads.size();
-   int numAtom = molecularBeads[0]->GetNumberAtoms();
+   int numAtom = molecularBeads[0]->GetAtomVect().size();
    double beadsKineticEnergy = 0.0;;
    for(int b=0; b<numBeads; b++){
       double coreKineticEnergy = 0.0;
       for(int a=0; a<numAtom; a++){
-         Atom* atom = molecularBeads[b]->GetAtom(a);
+         Atom* atom = molecularBeads[b]->GetAtomVect()[a];
          double coreMass = atom->GetCoreMass();
          for(int i=0; i<CartesianType_end; i++){
             coreKineticEnergy += 0.5*pow(atom->GetPxyz()[i],2.0)/coreMass;
@@ -269,8 +269,8 @@ double RPMD::OutputEnergies(const vector<boost::shared_ptr<Molecule> >& molecula
       double harmonicEnergy = 0.0;
       int preB = b==0 ? numBeads-1 : b-1;
       for(int a=0; a<numAtom; a++){
-         Atom* atom    = molecularBeads[b]->GetAtom(a);
-         Atom* preAtom = molecularBeads[preB]->GetAtom(a);
+         Atom* atom    = molecularBeads[b]->GetAtomVect()[a];
+         Atom* preAtom = molecularBeads[preB]->GetAtomVect()[a];
          double coreMass = atom->GetCoreMass();
          double dx = atom->GetXyz()[XAxis] - preAtom->GetXyz()[XAxis];
          double dy = atom->GetXyz()[YAxis] - preAtom->GetXyz()[YAxis];

--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -91,9 +91,9 @@ ZindoS::ZindoS() : MolDS_cndo::Cndo2(){
 ZindoS::~ZindoS(){
    if(this->theory==ZINDOS){
       MallocerFreer::GetInstance()->Free<double>(&this->nishimotoMatagaMatrix, 
-                                                 this->molecule->GetNumberAtoms(), 
+                                                 this->molecule->GetAtomVect().size(), 
                                                  OrbitalType_end, 
-                                                 this->molecule->GetNumberAtoms(), 
+                                                 this->molecule->GetAtomVect().size(), 
                                                  OrbitalType_end);
    }
    MallocerFreer::GetInstance()->Free<double>(&this->matrixCIS, 
@@ -105,7 +105,7 @@ ZindoS::~ZindoS(){
                                               this->matrixCISdimension);
    MallocerFreer::GetInstance()->Free<double>(&this->matrixForce, 
                                               this->matrixForceElecStatesNum,
-                                              this->molecule->GetNumberAtoms(),
+                                              this->molecule->GetAtomVect().size(),
                                               CartesianType_end);
    MallocerFreer::GetInstance()->Free<double>(&this->zMatrixForce, 
                                               this->zMatrixForceElecStatesNum,
@@ -123,11 +123,11 @@ ZindoS::~ZindoS(){
                                                  this->molecule->GetTotalNumberAOs());
       MallocerFreer::GetInstance()->Free<double>(&this->atomicElectronPopulationCIS, 
                                                  elecStates->size(),
-                                                 this->molecule->GetNumberAtoms());
+                                                 this->molecule->GetAtomVect().size());
       if(Parameters::GetInstance()->RequiresUnpairedPopCIS()){
          MallocerFreer::GetInstance()->Free<double>(&this->atomicUnpairedPopulationCIS, 
                                                     elecStates->size(),
-                                                    this->molecule->GetNumberAtoms());
+                                                    this->molecule->GetAtomVect().size());
       }
    }
    //this->OutputLog("ZindoS deleted\n");
@@ -137,9 +137,9 @@ void ZindoS::SetMolecule(Molecule* molecule){
    Cndo2::SetMolecule(molecule);
    if(this->theory==ZINDOS){
       MallocerFreer::GetInstance()->Malloc<double>(&this->nishimotoMatagaMatrix, 
-                                                   this->molecule->GetNumberAtoms(), 
+                                                   this->molecule->GetAtomVect().size(), 
                                                    OrbitalType_end, 
-                                                   this->molecule->GetNumberAtoms(), 
+                                                   this->molecule->GetAtomVect().size(), 
                                                    OrbitalType_end);
    }
 }
@@ -255,10 +255,10 @@ double ZindoS::GetFockDiagElement(const Atom& atomA,
       value += temp;
    
       temp = 0.0;
-      int totalNumberAtoms = molecule.GetNumberAtoms();
+      int totalNumberAtoms = molecule.GetAtomVect().size();
       for(int B=0; B<totalNumberAtoms; B++){
          if(B != indexAtomA){
-            const Atom& atomB = *molecule.GetAtom(B);
+            const Atom& atomB = *molecule.GetAtomVect()[B];
             OrbitalType orbitalSigma;
             int sigma;
             int atomBNumberValence = atomB.GetValenceSize();
@@ -732,18 +732,18 @@ double ZindoS::GetNishimotoMatagaTwoEleInt1stDerivative(const Atom& atomA,
 }
 
 void ZindoS::CalcNishimotoMatagaMatrix(double**** nishimotoMatagaMatrix, const Molecule& molecule) const{
-   int totalNumberAtoms = molecule.GetNumberAtoms();
+   int totalNumberAtoms = molecule.GetAtomVect().size();
    stringstream ompErrors;
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
    for(int A=0; A<totalNumberAtoms; A++){
       try{
-         const Atom& atomA = *molecule.GetAtom(A);
+         const Atom& atomA = *molecule.GetAtomVect()[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);
+               const Atom& atomB = *molecule.GetAtomVect()[B];
                int firstAOIndexB = atomB.GetFirstAOIndex();
                int lastAOIndexB  = atomB.GetLastAOIndex();
                double rAB = molecule.GetDistanceAtoms(atomA, atomB);
@@ -956,8 +956,8 @@ double ZindoS::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
    double gamma;
    double exchange;
    double coulomb;
-   for(int A=0; A<molecule.GetNumberAtoms(); A++){
-      const Atom& atomA = *molecule.GetAtom(A);
+   for(int A=0; A<molecule.GetAtomVect().size(); A++){
+      const Atom& atomA = *molecule.GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
 
@@ -965,8 +965,8 @@ double ZindoS::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
          OrbitalType orbitalMu = atomA.GetValence(mu-firstAOIndexA);
 
          // CNDO term
-         for(int B=A; B<molecule.GetNumberAtoms(); B++){
-            const Atom& atomB = *molecule.GetAtom(B);
+         for(int B=A; B<molecule.GetAtomVect().size(); B++){
+            const Atom& atomB = *molecule.GetAtomVect()[B];
             int firstAOIndexB = atomB.GetFirstAOIndex();
             int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -1700,7 +1700,7 @@ void ZindoS::CalcAtomicElectronPopulationCIS(double*** atomicElectronPopulationC
    if(!Parameters::GetInstance()->RequiresMullikenCIS()){
       return;
    }
-   int totalNumberAtoms = molecule.GetNumberAtoms();
+   int totalNumberAtoms = molecule.GetAtomVect().size();
    vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS();
    // malloc or initialize free exciton energies
    if(*atomicElectronPopulationCIS == NULL){
@@ -1719,8 +1719,8 @@ void ZindoS::CalcAtomicElectronPopulationCIS(double*** atomicElectronPopulationC
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
       for(int a=0; a<totalNumberAtoms; a++){
          try{
-            int firstAOIndex = molecule.GetAtom(a)->GetFirstAOIndex();
-            int numberAOs = molecule.GetAtom(a)->GetValenceSize();
+            int firstAOIndex = molecule.GetAtomVect()[a]->GetFirstAOIndex();
+            int numberAOs = molecule.GetAtomVect()[a]->GetValenceSize();
             for(int i=firstAOIndex; i<firstAOIndex+numberAOs; i++){
                (*atomicElectronPopulationCIS)[k][a] += orbitalElectronPopulationCIS[k][i][i];
             }
@@ -1746,7 +1746,7 @@ void ZindoS::CalcAtomicUnpairedPopulationCIS(double*** atomicUnpairedPopulationC
    if(!Parameters::GetInstance()->RequiresUnpairedPopCIS()){
       return;
    }
-   int totalNumberAtoms = molecule.GetNumberAtoms();
+   int totalNumberAtoms = molecule.GetAtomVect().size();
    vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS();
    // malloc or initialize free exciton energies
    if(*atomicUnpairedPopulationCIS == NULL){
@@ -1765,8 +1765,8 @@ void ZindoS::CalcAtomicUnpairedPopulationCIS(double*** atomicUnpairedPopulationC
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
       for(int a=0; a<totalNumberAtoms; a++){
          try{
-            int firstAOIndex = molecule.GetAtom(a)->GetFirstAOIndex();
-            int numberAOs = molecule.GetAtom(a)->GetValenceSize();
+            int firstAOIndex = molecule.GetAtomVect()[a]->GetFirstAOIndex();
+            int numberAOs = molecule.GetAtomVect()[a]->GetValenceSize();
             (*atomicUnpairedPopulationCIS)[k][a] = 0.0; 
             for(int i=firstAOIndex; i<firstAOIndex+numberAOs; i++){
                double orbitalSquarePopulation = 0.0; 
@@ -1957,12 +1957,12 @@ void ZindoS::OutputCISMulliken() const{
    if(!Parameters::GetInstance()->RequiresMullikenCIS()){
       return;
    }
-   int totalNumberAtoms = this->molecule->GetNumberAtoms();
+   int totalNumberAtoms = this->molecule->GetAtomVect().size();
    this->OutputLog(this->messageMullikenAtomsTitle);
    vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS();
    for(int k=0; k<elecStates->size(); k++){
       for(int a=0; a<totalNumberAtoms; a++){
-         const Atom& atom = *this->molecule->GetAtom(a);
+         const Atom& atom = *this->molecule->GetAtomVect()[a];
          this->OutputLog(boost::format("%s\t%d\t%d\t%s\t%e\t%e\n") % this->messageMullikenAtoms
                                                                    % (*elecStates)[k]
                                                                    % a
@@ -1981,12 +1981,12 @@ void ZindoS::OutputCISUnpairedPop() const{
    if(!Parameters::GetInstance()->RequiresUnpairedPopCIS()){
       return;
    }
-   int totalNumberAtoms = this->molecule->GetNumberAtoms();
+   int totalNumberAtoms = this->molecule->GetAtomVect().size();
    this->OutputLog(this->messageUnpairedAtomsTitle);
    vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS();
    for(int k=0; k<elecStates->size(); k++){
       for(int a=0; a<totalNumberAtoms; a++){
-         const Atom& atom = *this->molecule->GetAtom(a);
+         const Atom& atom = *this->molecule->GetAtomVect()[a];
          this->OutputLog(boost::format("%s\t%d\t%d\t%s\t%e\n") % this->messageUnpairedAtoms
                                                                % (*elecStates)[k]
                                                                % a
@@ -2447,9 +2447,9 @@ double ZindoS::GetCISDiagElement(double const* energiesMO,
    double gamma    = 0.0;
    double exchange = 0.0;
    double coulomb  = 0.0;
-   int totalNumberAtoms = molecule.GetNumberAtoms();
+   int totalNumberAtoms = molecule.GetAtomVect().size();
    for(int A=0; A<totalNumberAtoms; A++){
-      const Atom& atomA = *molecule.GetAtom(A);
+      const Atom& atomA = *molecule.GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
 
@@ -2461,7 +2461,7 @@ double ZindoS::GetCISDiagElement(double const* energiesMO,
 
          // CNDO term
          for(int B=A; B<totalNumberAtoms; B++){
-            const Atom& atomB = *molecule.GetAtom(B);
+            const Atom& atomB = *molecule.GetAtomVect()[B];
             int firstAOIndexB = atomB.GetFirstAOIndex();
             int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -2540,9 +2540,9 @@ double ZindoS::GetCISOffDiagElement(double const* const* const* const* nishimoto
    double gamma    = 0.0;
    double exchange = 0.0;
    double coulomb  = 0.0;
-   int totalNumberAtoms = molecule.GetNumberAtoms();
+   int totalNumberAtoms = molecule.GetAtomVect().size();
    for(int A=0; A<totalNumberAtoms; A++){
-      const Atom& atomA = *molecule.GetAtom(A);
+      const Atom& atomA = *molecule.GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
 
@@ -2555,7 +2555,7 @@ double ZindoS::GetCISOffDiagElement(double const* const* const* const* nishimoto
 
          // CNDO term
          for(int B=A; B<totalNumberAtoms; B++){
-            const Atom& atomB = *molecule.GetAtom(B);
+            const Atom& atomB = *molecule.GetAtomVect()[B];
             int firstAOIndexB = atomB.GetFirstAOIndex(); 
             int lastAOIndexB  = atomB.GetLastAOIndex(); 
             for(int nu=firstAOIndexB; nu<=lastAOIndexB; nu++){
@@ -2636,7 +2636,7 @@ void ZindoS::CheckMatrixForce(const vector<int>& elecStates){
    if(this->matrixForce == NULL){
       MallocerFreer::GetInstance()->Malloc<double>(&this->matrixForce, 
                                                    elecStates.size(),
-                                                   this->molecule->GetNumberAtoms(), 
+                                                   this->molecule->GetAtomVect().size(), 
                                                    CartesianType_end);
       this->matrixForceElecStatesNum = elecStates.size();
    }
@@ -2644,7 +2644,7 @@ void ZindoS::CheckMatrixForce(const vector<int>& elecStates){
       MallocerFreer::GetInstance()->
       Initialize<double>(this->matrixForce,
                          elecStates.size(),
-                         this->molecule->GetNumberAtoms(),
+                         this->molecule->GetAtomVect().size(),
                          CartesianType_end);
    }
 }
@@ -3121,13 +3121,13 @@ double ZindoS::GetSmallQElement(int moI,
    int numberOcc = this->molecule->GetTotalNumberValenceElectrons()/2;
    bool isMoPOcc = moP<numberOcc ? true : false;
    
-   for(int A=0; A<molecule->GetNumberAtoms(); A++){
-      const Atom& atomA = *molecule->GetAtom(A);
+   for(int A=0; A<molecule->GetAtomVect().size(); A++){
+      const Atom& atomA = *molecule->GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
 
-      for(int B=A; B<molecule->GetNumberAtoms(); B++){
-         const Atom& atomB = *molecule->GetAtom(B);
+      for(int B=A; B<molecule->GetAtomVect().size(); B++){
+         const Atom& atomB = *molecule->GetAtomVect()[B];
          int firstAOIndexB = atomB.GetFirstAOIndex();
          int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -3518,8 +3518,8 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
 
    // Fast algorith, but this is not easy to read. 
    // Slow algorithm is alos written below.
-   for(int A=0; A<this->molecule->GetNumberAtoms(); A++){
-      const Atom& atomA = *this->molecule->GetAtom(A);
+   for(int A=0; A<this->molecule->GetAtomVect().size(); A++){
+      const Atom& atomA = *this->molecule->GetAtomVect()[A];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
 
@@ -3566,8 +3566,8 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
          }
 
          //  (A==B || A!=B) && (mu==nu && lambda==sigma for (mu nu|lamba sigma))
-         for(int B=A; B<this->molecule->GetNumberAtoms(); B++){
-            const Atom& atomB = *this->molecule->GetAtom(B);
+         for(int B=A; B<this->molecule->GetAtomVect().size(); B++){
+            const Atom& atomB = *this->molecule->GetAtomVect()[B];
             int firstAOIndexB = atomB.GetFirstAOIndex();
             int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -3624,10 +3624,10 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
 void ZindoS::CalcDiatomicTwoElecsTwoCores1stDerivatives(double*** matrix,  
                                                         int indexAtomA, 
                                                         int indexAtomB) const{
-   const Atom& atomA = *molecule->GetAtom(indexAtomA);
+   const Atom& atomA = *molecule->GetAtomVect()[indexAtomA];
    const int firstAOIndexA = atomA.GetFirstAOIndex();
    const int lastAOIndexA  = atomA.GetLastAOIndex();
-   const Atom& atomB = *molecule->GetAtom(indexAtomB);
+   const Atom& atomB = *molecule->GetAtomVect()[indexAtomB];
    const int firstAOIndexB = atomB.GetFirstAOIndex();
    const int lastAOIndexB  = atomB.GetLastAOIndex();
    const double rAB = this->molecule->GetDistanceAtoms(atomA, atomB);
@@ -3657,9 +3657,9 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
    }
   
    // this loop is MPI-parallelized
-   for(int a=0; a<this->molecule->GetNumberAtoms(); a++){
+   for(int a=0; a<this->molecule->GetAtomVect().size(); a++){
       if(a%mpiSize != mpiRank){continue;}
-      const Atom& atomA = *molecule->GetAtom(a);
+      const Atom& atomA = *molecule->GetAtomVect()[a];
       int firstAOIndexA = atomA.GetFirstAOIndex();
       int lastAOIndexA  = atomA.GetLastAOIndex();
       stringstream ompErrors;
@@ -3689,9 +3689,9 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
                                         &tmpMatrixBC,
                                         &tmpVectorBC);
 #pragma omp for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
-            for(int b=0; b<this->molecule->GetNumberAtoms(); b++){
+            for(int b=0; b<this->molecule->GetAtomVect().size(); b++){
                if(a == b){continue;}
-               const Atom& atomB = *molecule->GetAtom(b);
+               const Atom& atomB = *molecule->GetAtomVect()[b];
                int firstAOIndexB = atomB.GetFirstAOIndex();
                int lastAOIndexB  = atomB.GetLastAOIndex();
                double rAB = this->molecule->GetDistanceAtoms(atomA, atomB);
@@ -3839,7 +3839,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
    }    // end of for(int a) with MPI parallelization
 
    // communication to reduce thsi->matrixForce on all node (namely, all_reduce)
-   int numTransported = elecStates.size()*this->molecule->GetNumberAtoms()*CartesianType_end;
+   int numTransported = elecStates.size()*this->molecule->GetAtomVect().size()*CartesianType_end;
    MolDS_mpi::MpiProcess::GetInstance()->AllReduce(&this->matrixForce[0][0][0], numTransported, std::plus<double>());
 
    /*
@@ -3848,9 +3848,9 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
    // calculated with GTO expansion technique.
    stringstream ompErrors;
 #pragma omp parallel for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
-   for(int a=0; a<this->molecule->GetNumberAtoms(); a++){
+   for(int a=0; a<this->molecule->GetAtomVect().size(); a++){
       try{
-         const Atom& atomA = *molecule->GetAtom(a);
+         const Atom& atomA = *molecule->GetAtomVect()[a];
          int firstAOIndexA = atomA.GetFirstAOIndex();
          int lastAOIndexA  = atomA.GetLastAOIndex();
          for(int i=0; i<CartesianType_end; i++){
@@ -3859,9 +3859,9 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
             double forceElecCoreAttPart = 0.0;
             double forceOverlapAOsPart = 0.0;
             double forceTwoElecPart = 0.0;
-            for(int b=0; b<this->molecule->GetNumberAtoms(); b++){
+            for(int b=0; b<this->molecule->GetAtomVect().size(); b++){
                if(a != b){
-                  const Atom& atomB = *molecule->GetAtom(b);
+                  const Atom& atomB = *molecule->GetAtomVect()[b];
                   int firstAOIndexB = atomB.GetFirstAOIndex();
                   int lastAOIndexB  = atomB.GetLastAOIndex();
 
@@ -3978,8 +3978,8 @@ void ZindoS::CalcForceExcitedStaticPart(double* force,
                                       int indexAtomA, 
                                       int indexAtomB,
                                       double const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexB = atomB.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
@@ -4005,8 +4005,8 @@ void ZindoS::CalcForceExcitedElecCoreAttractionPart(double* force,
                                                   int indexAtomA, 
                                                   int indexAtomB,
                                                   double const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
    for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
@@ -4023,8 +4023,8 @@ void ZindoS::CalcForceExcitedOverlapAOsPart(double* force,
                                             int indexAtomA, 
                                             int indexAtomB,
                                             double const* const* const* diatomicOverlapAOs1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexB = atomB.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();
@@ -4053,8 +4053,8 @@ void ZindoS::CalcForceExcitedTwoElecPart(double* force,
                                          int indexAtomA, 
                                          int indexAtomB,
                                          double const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
-   const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
-   const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
+   const Atom& atomA = *this->molecule->GetAtomVect()[indexAtomA];
+   const Atom& atomB = *this->molecule->GetAtomVect()[indexAtomB];
    int firstAOIndexA = atomA.GetFirstAOIndex();
    int firstAOIndexB = atomB.GetFirstAOIndex();
    int lastAOIndexA  = atomA.GetLastAOIndex();