#### 2. Molecule Class

```cpp
#include <vector>
#include <string>

class Molecule {
public:
    struct Atom {
        std::string element;
        double x, y, z;
    };

    void addAtom(const std::string& element, double x, double y, double z);
    void setCharge(int charge);
    void setMultiplicity(int multiplicity);

    const std::vector<Atom>& getAtoms() const;
    int getCharge() const;
    int getMultiplicity() const;

private:
    std::vector<Atom> atoms;
    int charge;
    int multiplicity;
};

void Molecule::addAtom(const std::string& element, double x, double y, double z) {
    atoms.push_back({ element, x, y, z });
}

void Molecule::setCharge(int charge) {
    this->charge = charge;
}

void Molecule::setMultiplicity(int multiplicity) {
    this->multiplicity = multiplicity;
}

const std::vector<Molecule::Atom>& Molecule::getAtoms() const {
    return atoms;
}

int Molecule::getCharge() const {
    return charge;
}

int Molecule::getMultiplicity() const {
    return multiplicity;
}
```

#### 3. Energy Class

No changes needed here.

```cpp
#include <vector>

class Energy {
public:
    void addEnergy(double energy);
    const std::vector<double>& getEnergies() const;

private:
    std::vector<double> energies;
};

void Energy::addEnergy(double energy) {
    energies.push_back(energy);
}

const std::vector<double>& Energy::getEnergies() const {
    return energies;
}
```

#### 4. VibrationalFrequency Class

No changes needed here.

```cpp
#include <vector>

class VibrationalFrequency {
public:
    void addFrequency(double frequency);
    const std::vector<double>& getFrequencies() const;

private:
    std::vector<double> frequencies;
};

void VibrationalFrequency::addFrequency(double frequency) {
    frequencies.push_back(frequency);
}

const std::vector<double>& VibrationalFrequency::getFrequencies() const {
    return frequencies;
}
```



#### 1. LogFile Class

This class will now handle parsing and storing additional types of data.

```cpp
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <regex>
#include "Molecule.h"
#include "Energy.h"
#include "VibrationalFrequency.h"
#include "Thermochemistry.h"
#include "SCFData.h"
#include "BasisSet.h"
#include "Method.h"

class LogFile {
public:
    LogFile(const std::string& filename) : filename(filename) {}

    bool readFile();
    void parseFile();

    const Molecule& getMolecule() const { return molecule; }
    const Energy& getEnergy() const { return energy; }
    const VibrationalFrequency& getVibrationalFrequency() const { return vibrationalFrequency; }
    const Thermochemistry& getThermochemistry() const { return thermochemistry; }
    const SCFData& getSCFData() const { return scfData; }
    const BasisSet& getBasisSet() const { return basisSet; }
    const Method& getMethod() const { return method; }

private:
    std::string filename;
    std::vector<std::string> lines;
    Molecule molecule;
    Energy energy;
    VibrationalFrequency vibrationalFrequency;
    Thermochemistry thermochemistry;
    SCFData scfData;
    BasisSet basisSet;
    Method method;

    void parseEnergy(const std::smatch& match);
    void parseMolecule(const std::smatch& match);
    void parseVibrationalFrequency(const std::smatch& match);
    void parseThermochemistry(const std::smatch& match);
    void parseSCFData(const std::smatch& match);
    void parseBasisSet(const std::smatch& match);
    void parseMethod(const std::smatch& match);
};

bool LogFile::readFile() {
    std::ifstream file(filename);
    if (!file.is_open()) {
        std::cerr << "Could not open file: " << filename << std::endl;
        return false;
    }

    std::string line;
    while (std::getline(file, line)) {
        lines.push_back(line);
    }
    file.close();
    return true;
}

void LogFile::parseFile() {
    std::regex energyRegex(R"(SCF Done:.*=\s*(-?\d+\.\d+))");
    std::regex atomRegex(R"(\s+(\d+)\s+(\d+)\s+(\d+)\s+([A-Z][a-z]?)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+))");
    std::regex freqRegex(R"(Frequencies --\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+))");
    std::regex thermoRegex(R"(Thermochemistry\s+.*\s+.*\s+.*)");
    std::regex scfRegex(R"(SCF Iteration\s+(\d+)\s+Energy\s+=\s+(-?\d+\.\d+))");
    std::regex basisSetRegex(R"(Basis set:([^\n]+))");
    std::regex methodRegex(R"(Method:([^\n]+))");

    for (const auto& line : lines) {
        std::smatch match;

        if (std::regex_search(line, match, energyRegex)) {
            parseEnergy(match);
        } else if (std::regex_search(line, match, atomRegex)) {
            parseMolecule(match);
        } else if (std::regex_search(line, match, freqRegex)) {
            parseVibrationalFrequency(match);
        } else if (std::regex_search(line, match, thermoRegex)) {
            parseThermochemistry(match);
        } else if (std::regex_search(line, match, scfRegex)) {
            parseSCFData(match);
        } else if (std::regex_search(line, match, basisSetRegex)) {
            parseBasisSet(match);
        } else if (std::regex_search(line, match, methodRegex)) {
            parseMethod(match);
        }
    }
}

void LogFile::parseEnergy(const std::smatch& match) {
    energy.addEnergy(std::stod(match[1].str()));
}

void LogFile::parseMolecule(const std::smatch& match) {
    molecule.addAtom(match[4].str(), std::stod(match[5].str()), std::stod(match[6].str()), std::stod(match[7].str()));
}

void LogFile::parseVibrationalFrequency(const std::smatch& match) {
    vibrationalFrequency.addFrequency(std::stod(match[1].str()));
    vibrationalFrequency.addFrequency(std::stod(match[2].str()));
    vibrationalFrequency.addFrequency(std::stod(match[3].str()));
}

void LogFile::parseThermochemistry(const std::smatch& match) {
    // Implement logic to parse thermochemistry data
    // Example: extract specific thermochemical properties
}

void LogFile::parseSCFData(const std::smatch& match) {
    int iteration = std::stoi(match[1].str());
    double energy = std::stod(match[2].str());
    scfData.addIterationData(iteration, energy);
}

void LogFile::parseBasisSet(const std::smatch& match) {
    basisSet.setBasisSet(match[1].str());
}

void LogFile::parseMethod(const std::smatch& match) {
    method.setMethod(match[1].str());
}
```

#### 5. Thermochemistry Class

Extend the class to handle specific thermochemical properties.

```cpp
#include <vector>
#include <string>

class Thermochemistry {
public:
    struct Data {
        std::string property;
        double value;
    };

    void addData(const std::string& property, double value);
    const std::vector<Data>& getData() const;
    
    void setEnthalpy(double enthalpy);
    void setEntropy(double entropy);
    void setFreeEnergy(double freeEnergy);

    double getEnthalpy() const;
    double getEntropy() const;
    double getFreeEnergy() const;

private:
    std::vector<Data> data;
    double enthalpy;
    double entropy;
    double freeEnergy;
};

void Thermochemistry::addData(const std::string& property, double value) {
    data.push_back({ property, value });
}

const std::vector<Thermochemistry::Data>& Thermochemistry::getData() const {
    return data;
}

void Thermochemistry::setEnthalpy(double enthalpy) {
    this->enthalpy = enthalpy;
}

void Thermochemistry::setEntropy(double entropy) {
    this->entropy = entropy;
}

void Thermochemistry::setFreeEnergy(double freeEnergy) {
    this->freeEnergy = freeEnergy;
}

double Thermochemistry::getEnthalpy() const {
    return enthalpy;
}

double Thermochemistry::getEntropy() const {
    return entropy;
}

double Thermochemistry::getFreeEnergy() const {
    return freeEnergy;
}
```

#### 6. SCFData Class

Add a class to handle SCF (Self-Consistent Field) iteration data.

```cpp
#include <vector>

class SCFData {
public:
    struct IterationData {
        int iteration;
        double energy;
    };

    void addIterationData(int iteration, double energy);
    const std::vector<IterationData>& getIterationData() const;

private:
    std::vector<IterationData> iterationData;
};

void SCFData::addIterationData(int iteration, double energy) {
    iterationData.push_back({ iteration, energy });
}

const std::vector<SCFData::IterationData>& SCFData::getIterationData() const {
    return iterationData;
}
```

#### 7. BasisSet Class

Add a class to handle basis set information.

```cpp
#include <string>

class BasisSet {
public:
    void setBasisSet(const std::string& basisSet);
    const std::string& getBasisSet() const;

private:
    std::string basisSet;
};

void BasisSet::setBasisSet(const std::string& basisSet) {
    this->basisSet = basisSet;
}

const std::string& BasisSet::getBasisSet() const {
    return basisSet;
}
```

#### 8. Method Class

Add a class to handle computational method information.

```cpp
#include <string>

class Method {
public:
    void setMethod(const std::string& method);
    const std::string& getMethod() const;

private:
    std::string method;
}

void Method::setMethod(const std::string& method) {
    this->method = method;
}

const std::string& Method::getMethod() const {
    return method;
}
```

### Example Usage

Here's how you might use these classes in a main function:

```cpp
#include "LogFile.h

```cpp
#include "Molecule.h"
#include "Energy.h"
#include "VibrationalFrequency.h"
#include "Thermochemistry.h"
#include "SCFData.h"
#include "BasisSet.h"
#include "Method.h"

int main() {
    LogFile logFile("path/to/gaussian.log");

    if (!logFile.readFile()) {
        return 1;
    }

    logFile.parseFile();

    const Molecule& molecule = logFile.getMolecule();
    const Energy& energy = logFile.getEnergy();
    const VibrationalFrequency& vibrationalFrequency = logFile.getVibrationalFrequency();
    const Thermochemistry& thermochemistry = logFile.getThermochemistry();
    const SCFData& scfData = logFile.getSCFData();
    const BasisSet& basisSet = logFile.getBasisSet();
    const Method& method = logFile.getMethod();

    // Output parsed data for the molecule
    std::cout << "Molecule Data:" << std::endl;
    for (const auto& atom : molecule.getAtoms()) {
        std::cout << atom.element << ": (" << atom.x << ", " << atom.y << ", " << atom.z << ")\n";
    }
    std::cout << "Charge: " << molecule.getCharge() << std::endl;
    std::cout << "Multiplicity: " << molecule.getMultiplicity() << std::endl;

    // Output parsed energy data
    std::cout << "\nEnergy Data:" << std::endl;
    for (const auto& e : energy.getEnergies()) {
        std::cout << "Energy: " << e << "\n";
    }

    // Output parsed vibrational frequency data
    std::cout << "\nVibrational Frequencies:" << std::endl;
    for (const auto& freq : vibrationalFrequency.getFrequencies()) {
        std::cout << "Frequency: " << freq << "\n";
    }

    // Output parsed thermochemistry data
    std::cout << "\nThermochemistry Data:" << std::endl;
    std::cout << "Enthalpy: " << thermochemistry.getEnthalpy() << std::endl;
    std::cout << "Entropy: " << thermochemistry.getEntropy() << std::endl;
    std::cout << "Free Energy: " << thermochemistry.getFreeEnergy() << std::endl;
    for (const auto& data : thermochemistry.getData()) {
        std::cout << data.property << ": " << data.value << "\n";
    }

    // Output parsed SCF data
    std::cout << "\nSCF Data:" << std::endl;
    for (const auto& iteration : scfData.getIterationData()) {
        std::cout << "Iteration: " << iteration.iteration << " Energy: " << iteration.energy << "\n";
    }

    // Output basis set information
    std::cout << "\nBasis Set: " << basisSet.getBasisSet() << std::endl;

    // Output method information
    std::cout << "\nMethod: " << method.getMethod() << std::endl;

    return 0;
}
```

### Updated Parsing Logic for Thermochemistry

Let's update the `parseThermochemistry` method to extract specific thermochemical properties like enthalpy, entropy, and free energy.

```cpp
void LogFile::parseThermochemistry(const std::smatch& match) {
    std::regex enthalpyRegex(R"(Enthalpy\s+=\s+(-?\d+\.\d+))");
    std::regex entropyRegex(R"(Entropy\s+=\s+(-?\d+\.\d+))");
    std::regex freeEnergyRegex(R"(Free Energy\s+=\s+(-?\d+\.\d+))");

    for (const auto& line : lines) {
        std::smatch match;

        if (std::regex_search(line, match, enthalpyRegex)) {
            thermochemistry.setEnthalpy(std::stod(match[1].str()));
        } else if (std::regex_search(line, match, entropyRegex)) {
            thermochemistry.setEntropy(std::stod(match[1].str()));
        } else if (std::regex_search(line, match, freeEnergyRegex)) {
            thermochemistry.setFreeEnergy(std::stod(match[1].str()));
        }
    }
}
```

fix parse functions to parse gaussian log file, do not use regex, here reference file
use following code as reference
bool GaussianOutputFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {
    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if (pmol == nullptr)
      return false;

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    const char* title = pConv->GetTitle();

    char buffer[BUFF_SIZE];
    string str,str1,str2,thermo_method;
    double x,y,z;
    OBAtom *atom;
    vector<string> vs,vs2;
    int total_charge = 0;
    unsigned int spin_multiplicity = 1;
    bool hasPartialCharges = false;
    string chargeModel; // descriptor for charges (e.g. "Mulliken")

    // Variable for G2/G3/G4 etc. calculations
    double ezpe,Hcorr,Gcorr,E0,CV;
    bool ezpe_set=false,Hcorr_set=false,Gcorr_set=false,E0_set=false,CV_set=false;
    double temperature = 0; /* Kelvin */
    std::vector<double> Scomponents;
    // Electrostatic potential. ESP is calculated
    // once unless the Opt and Pop jobs are combined.
    // In this case, ESP is calculated once before
    // the geometry optmization and once after. If this
    // happens, the second ESP must be added to OBMol.
    OBFreeGrid *esp   = nullptr;
    int NumEsp        = 1;
    int NumEspCounter = 0;
    bool ESPisAdded   = false;

    // coordinates of all steps
    // Set conformers to all coordinates we adopted
    std::vector<double*> vconf; // index of all frames/conformers
    std::vector<double> coordinates; // coordinates in each frame
    int natoms = 0; // number of atoms -- ensure we don't go to a new job with a different molecule

    // OBConformerData stores information about multiple steps
    // we can change attribute later if needed (e.g., IRC)
    OBConformerData *confData = new OBConformerData();
    confData->SetOrigin(fileformatInput);
    std::vector<unsigned short> confDimensions = confData->GetDimension(); // to be fair, set these all to 3D
    std::vector<double>         confEnergies   = confData->GetEnergies();
    std::vector< std::vector< vector3 > > confForces = confData->GetForces();

    //Vibrational data
    std::vector< std::vector< vector3 > > Lx;
    std::vector<double> Frequencies, Intensities, RamanActivities;
    //Rotational data
    std::vector<double> RotConsts(3);
    int RotSymNum=1;
    OBRotationData::RType RotorType = OBRotationData::UNKNOWN;

    // Translation vectors (if present)
    vector3 translationVectors[3];
    int numTranslationVectors = 0;

    //Electronic Excitation data
    std::vector<double> Forces, Wavelengths, EDipole,
      RotatoryStrengthsVelocity, RotatoryStrengthsLength;

    // Orbital data
    std::vector<double> orbitals;
    std::vector<std::string> symmetries;
    int aHOMO, bHOMO, betaStart;
    aHOMO = bHOMO = betaStart = -1;

    int i=0;
    bool no_symmetry=false;
    char coords_type[25];

    //Prescan file to find second instance of "orientation:"
    //This will be the kind of coords used in the chk/fchk file
    //Unless the "nosym" keyword has been requested
    while (ifs.getline(buffer,BUFF_SIZE))
      {
        if (strstr(buffer, "Symmetry turned off by external request.") != nullptr)
          {
            // The "nosym" keyword has been requested
            no_symmetry = true;
          }
        if (strstr(buffer, "orientation:") != nullptr && (strstr(buffer, "Dipole") == nullptr))
          {
            i++;
            tokenize (vs, buffer);
            // gotta check what types of orientation are present
            strncpy (coords_type, vs[0].c_str(), 24);
            strcat (coords_type, " orientation:");
          }
        if ((no_symmetry && i==1) || i==2)
           break;
      }
    // Reset end-of-file pointers etc.
    ifs.clear();
    ifs.seekg(0);  //rewind

    mol.BeginModify();
    while (ifs.getline(buffer,BUFF_SIZE))
      {
        if(strstr(buffer, "Entering Gaussian") != nullptr)
        {
          //Put some metadata into OBCommentData
          string comment("Gaussian ");

          if (nullptr != strchr(buffer, '='))
            {
            comment += strchr(buffer,'=')+2;
            comment += "";
            for(unsigned i=0; i<115 && ifs; ++i)
            {
              ifs.getline(buffer,BUFF_SIZE);
              if (strstr(buffer, "Revision") != nullptr)
                {
                  if (buffer[strlen(buffer)-1] == ',')
                    {
                      buffer[strlen(buffer)-1] = '\0';
                    }
                  add_unique_pairdata_to_mol(&mol,"program",buffer,0);
                }
              else if(buffer[1]=='#')
              {
                //the line describing the method
                if (strstr(buffer, "Opt") != nullptr)
                {
                    // It is expected to have two sets of ESP in
                    // the log file if Opt is combined with Pop.
                    NumEsp = 2;
                }
                comment += buffer;
                OBCommentData *cd = new OBCommentData;
                cd->SetData(comment);
                cd->SetOrigin(fileformatInput);
                mol.SetData(cd);

                tokenize(vs,buffer);
                if (vs.size() > 1)
                  {
                    char *str = strdup(vs[1].c_str());
                    char *ptr = strchr(str,'/');

                    if (nullptr != ptr)
                      {
                        *ptr = ' ';
                        add_unique_pairdata_to_mol(&mol,"basis",ptr,0);
                        *ptr = '\0';
                        add_unique_pairdata_to_mol(&mol,"method",str,0);
                      }
                  }

                break;
              }
            }
          }
        }

        else if (strstr(buffer, "Multiplicity") != nullptr)
          {
            tokenize(vs, buffer, " \t\n");
            if (vs.size() == 6)
              {
                total_charge = atoi(vs[2].c_str());
                spin_multiplicity = atoi(vs[5].c_str());
              }

            ifs.getline(buffer,BUFF_SIZE);
          }
        else if (strstr(buffer, coords_type) != nullptr)
          {
            numTranslationVectors = 0; // ignore old translationVectors
            ifs.getline(buffer,BUFF_SIZE);      // ---------------
            ifs.getline(buffer,BUFF_SIZE);      // column headings
            ifs.getline(buffer,BUFF_SIZE);	// column headings
            ifs.getline(buffer,BUFF_SIZE);	// ---------------
            ifs.getline(buffer,BUFF_SIZE);
            tokenize(vs,buffer);
            while (vs.size()>4)
              {
                int corr = vs.size()==5 ? -1 : 0; //g94; later versions have an extra column
                x = atof((char*)vs[3+corr].c_str());
                y = atof((char*)vs[4+corr].c_str());
                z = atof((char*)vs[5+corr].c_str());
                int atomicNum = atoi((char*)vs[1].c_str());

                if (atomicNum > 0) // translation vectors are "-2"
                  {
                    if (natoms == 0) { // first time reading the molecule, create each atom
                      atom = mol.NewAtom();
                      atom->SetAtomicNum(atoi((char*)vs[1].c_str()));
                    }
                    coordinates.push_back(x);
                    coordinates.push_back(y);
                    coordinates.push_back(z);
                  }
                else {
                  translationVectors[numTranslationVectors++].Set(x, y, z);
                }

                if (!ifs.getline(buffer,BUFF_SIZE)) {
                  break;
                }
                tokenize(vs,buffer);
              }
            // done with reading atoms
            natoms = mol.NumAtoms();
            if(natoms==0)
              return false;
            // malloc / memcpy
            double *tmpCoords = new double [(natoms)*3];
            memcpy(tmpCoords, &coordinates[0], sizeof(double)*natoms*3);
            vconf.push_back(tmpCoords);
            coordinates.clear();
            confDimensions.push_back(3); // always 3D -- OBConformerData allows mixing 2D and 3D structures
          }
        else if(strstr(buffer, "Dipole moment") != nullptr)
            {
              ifs.getline(buffer,BUFF_SIZE); // actual components   X ###  Y #### Z ###
              tokenize(vs,buffer);
              if (vs.size() >= 6)
                {
                  OBVectorData *dipoleMoment = new OBVectorData;
                  dipoleMoment->SetAttribute("Dipole Moment");
                  double x, y, z;
                  x = atof(vs[1].c_str());
                  y = atof(vs[3].c_str());
                  z = atof(vs[5].c_str());
                  dipoleMoment->SetData(x, y, z);
                  dipoleMoment->SetOrigin(fileformatInput);
                  mol.SetData(dipoleMoment);
                }
              if (!ifs.getline(buffer,BUFF_SIZE)) break;
            }
        else if (strstr(buffer, "Traceless Quadrupole moment") != nullptr)
            {
              ifs.getline(buffer,BUFF_SIZE); // actual components XX ### YY #### ZZ ###
              tokenize(vs,buffer);
              ifs.getline(buffer,BUFF_SIZE); // actual components XY ### XZ #### YZ ###
              tokenize(vs2,buffer);
              if ((vs.size() >= 6) && (vs2.size() >= 6))
                {
                  double Q[3][3];
                  OpenBabel::OBMatrixData *quadrupoleMoment = new OpenBabel::OBMatrixData;

                  Q[0][0] = atof(vs[1].c_str());
                  Q[1][1] = atof(vs[3].c_str());
                  Q[2][2] = atof(vs[5].c_str());
                  Q[1][0] = Q[0][1] = atof(vs2[1].c_str());
                  Q[2][0] = Q[0][2] = atof(vs2[3].c_str());
                  Q[2][1] = Q[1][2] = atof(vs2[5].c_str());
                  matrix3x3 quad(Q);

                  quadrupoleMoment->SetAttribute("Traceless Quadrupole Moment");
                  quadrupoleMoment->SetData(quad);
                  quadrupoleMoment->SetOrigin(fileformatInput);
                  mol.SetData(quadrupoleMoment);
                }
              if (!ifs.getline(buffer,BUFF_SIZE)) break;
            }
        else if (strstr(buffer, "Exact polarizability") != nullptr)
            {
              // actual components XX, YX, YY, XZ, YZ, ZZ
              double xx, xy, yy, xz, yz, zz;
              const char *ptr = buffer+strlen("Exact polarizability:   ");
              if (ptr &&
                  6 == sscanf(ptr, "%8lf%8lf%8lf%8lf%8lf%8lf",
                              &xx, &xy, &yy, &xz, &yz, &zz))
              {
                  double Q[3][3];
                  OpenBabel::OBMatrixData *pol_tensor = new OpenBabel::OBMatrixData;

                  Q[0][0] = xx;
                  Q[1][1] = yy;
                  Q[2][2] = zz;
                  Q[1][0] = Q[0][1] = xy;
                  Q[2][0] = Q[0][2] = xz;
                  Q[2][1] = Q[1][2] = yz;
                  matrix3x3 pol(Q);

                  if (mol.HasData("Exact polarizability"))
                    {
                      mol.DeleteData("Exact polarizability"); // Delete the old one to add the new one
                    }
                  pol_tensor->SetAttribute("Exact polarizability");
                  pol_tensor->SetData(pol);
                  pol_tensor->SetOrigin(fileformatInput);
                  mol.SetData(pol_tensor);
                }
              if (!ifs.getline(buffer,BUFF_SIZE)) break;
            }
        else if(strstr(buffer, "Total atomic charges") != nullptr ||
                strstr(buffer, "Mulliken atomic charges") != nullptr ||
                strstr(buffer, "Mulliken charges:") != nullptr)
          {
            hasPartialCharges = true;
            chargeModel = "Mulliken";
            /*
              Gaussian usually calculates the electronic
              properties more than once, before and after
              geometry optimization. The second one is what
              we should be interested in. Thus, here, we
              delete the previously added Data to store the
              new one.
             */
            if (mol.HasData("Mulliken charges"))
              {
                mol.DeleteData("Mulliken charges");
              }
            OBPcharge *Mulliken = new OpenBabel::OBPcharge();
            std::vector<double> MPA_q;

            ifs.getline(buffer,BUFF_SIZE);	// column headings
            ifs.getline(buffer,BUFF_SIZE);
            tokenize(vs,buffer);
            while (vs.size() >= 3 &&
                   strstr(buffer, "Sum of ") == nullptr)
              {
                atom = mol.GetAtom(atoi(vs[0].c_str()));
                if (!atom)
                  break;
                atom->SetPartialCharge(atof(vs[2].c_str()));
                MPA_q.push_back(atof(vs[2].c_str()));
                if (!ifs.getline(buffer,BUFF_SIZE)) break;
                tokenize(vs,buffer);

              }
            if (MPA_q.size() == mol.NumAtoms())
            {
                Mulliken->AddPartialCharge(MPA_q);
                Mulliken->SetAttribute("Mulliken charges");
                Mulliken->SetOrigin(fileformatInput);
                mol.SetData(Mulliken);
            }
            else
            {
                cout << "Read " << MPA_q.size() << " Mulliken charges for " << mol.NumAtoms() << " atoms\n";
            }
          }
        else if (strstr(buffer, "Hirshfeld charges") != nullptr &&
                 strstr(buffer, "CM5 charges") != nullptr)
          {
            /*
              Hirshfeld and CM5 charges are printed in the
              same block in the Gaussian log file.
             */
            hasPartialCharges = true;
            chargeModel = "Hirshfeld";
            if (mol.HasData("Hirshfeld charges"))
              {
                mol.DeleteData("Hirshfeld charges");
              }
            if (mol.HasData("CM5 charges"))
              {
                mol.DeleteData("CM5 charges");
              }
            OBPcharge *Hirshfeld = new OpenBabel::OBPcharge();
            OBPcharge *CM5       = new OpenBabel::OBPcharge();
            std::vector<double> HPA_q;
            std::vector<double> CM5_q;
            ifs.getline(buffer,BUFF_SIZE);	// column headings
            ifs.getline(buffer,BUFF_SIZE);
            tokenize(vs,buffer);
            while (vs.size() >= 8 &&
                   strstr(buffer, "Tot ") == nullptr)
              {
                atom = mol.GetAtom(atoi(vs[0].c_str()));
                if (!atom)
                  break;
                atom->SetPartialCharge(atof(vs[2].c_str()));
                HPA_q.push_back(atof(vs[2].c_str()));
                CM5_q.push_back(atof(vs[7].c_str()));
                if (!ifs.getline(buffer,BUFF_SIZE)) break;
                tokenize(vs,buffer);

              }
            if (CM5_q.size() == mol.NumAtoms() and
                HPA_q.size() == mol.NumAtoms())
            {
                Hirshfeld->AddPartialCharge(HPA_q);
                Hirshfeld->SetAttribute("Hirshfeld charges");
                Hirshfeld->SetOrigin(fileformatInput);
                CM5->AddPartialCharge(CM5_q);
                CM5->SetAttribute("CM5 charges");
                CM5->SetOrigin(fileformatInput);
                mol.SetData(CM5);
                mol.SetData(Hirshfeld);
            }
            else
            {
                cout << "Read " << HPA_q.size() << " Hirshfeld charges for " << mol.NumAtoms() << " atoms\n";
            }
          }
        else if (strstr(buffer, "Electrostatic Properties Using The SCF Density") != nullptr)
          {
              NumEspCounter++;
          }
        else if (strstr(buffer, "Atomic Center") != nullptr && !ESPisAdded)
          {
            // Data points for ESP calculation
            tokenize(vs,buffer);
            if (nullptr == esp)
              esp = new OpenBabel::OBFreeGrid();
            if (vs.size() == 8)
              {
                esp->AddPoint(atof(vs[5].c_str()),atof(vs[6].c_str()),
                              atof(vs[7].c_str()),0);
              }
            else if (vs.size() > 5)
              {
                double x,y,z;
                if (3 == sscanf(buffer+32,"%10lf%10lf%10lf",&x,&y,&z))
                  {
                    esp->AddPoint(x,y,z,0);
                  }
              }
          }
        else if (strstr(buffer, "ESP Fit Center") != nullptr && !ESPisAdded)
          {
            // Data points for ESP calculation
            tokenize(vs,buffer);
            if (nullptr == esp)
              esp = new OpenBabel::OBFreeGrid();
            if (vs.size() == 9)
              {
                esp->AddPoint(atof(vs[6].c_str()),atof(vs[7].c_str()),
                              atof(vs[8].c_str()),0);
              }
            else if (vs.size() > 6)
              {
                double x,y,z;
                if (3 == sscanf(buffer+32,"%10lf%10lf%10lf",&x,&y,&z))
                  {
                    esp->AddPoint(x,y,z,0);
                  }
              }
          }
        else if (strstr(buffer, "Electrostatic Properties (Atomic Units)") != nullptr && !ESPisAdded)
          {
            int i,np;
            OpenBabel::OBFreeGridPoint *fgp;
            OpenBabel::OBFreeGridPointIterator fgpi;
            for(i=0; (i<5); i++)
              {
                ifs.getline(buffer,BUFF_SIZE);	// skip line
              }
            // Assume file is correct and that potentials are present
            // where they should.
            np = esp->NumPoints();
            fgpi = esp->BeginPoints();
            i = 0;
            for(fgp = esp->BeginPoint(fgpi); nullptr != fgp; fgp = esp->NextPoint(fgpi))
              {
                ifs.getline(buffer,BUFF_SIZE);
                tokenize(vs,buffer);
                if (vs.size() >= 2)
                  {
                    fgp->SetV(atof(vs[2].c_str()));
                    i++;
                  }
              }
            if (NumEsp == NumEspCounter)
              {
                if (i == np)
                  {
                    esp->SetAttribute("Electrostatic Potential");
                    esp->SetOrigin(fileformatInput);
                    mol.SetData(esp);
                    ESPisAdded = true;
                  }
                else
                  {
                    cout << "Read " << esp->NumPoints() << " ESP points i = " << i << "\n";
                  }
              }
            else if (!ESPisAdded)
              {
                esp->Clear();
              }
          }
        else if (strstr(buffer, "Charges from ESP fit") != nullptr)
          {
            hasPartialCharges = true;
            chargeModel = "ESP";
            if (mol.HasData("ESP charges"))
              {
                mol.DeleteData("ESP charges");
              }
            OBPcharge *ESP = new OpenBabel::OBPcharge();
            std::vector<double> ESP_q;
            ifs.getline(buffer,BUFF_SIZE);	// Charge / dipole line
            ifs.getline(buffer,BUFF_SIZE); // column header
            ifs.getline(buffer,BUFF_SIZE); // real charges
            tokenize(vs,buffer);
            while (vs.size() >= 3 &&
                   strstr(buffer, "-----") == nullptr)
              {
                atom = mol.GetAtom(atoi(vs[0].c_str()));
                if (!atom)
                  break;
                atom->SetPartialCharge(atof(vs[2].c_str()));
                ESP_q.push_back(atof(vs[2].c_str()));
                if (!ifs.getline(buffer,BUFF_SIZE)) break;
                tokenize(vs,buffer);
              }
            if (ESP_q.size() == mol.NumAtoms())
            {
                ESP->AddPartialCharge(ESP_q);
                ESP->SetAttribute("ESP charges");
                ESP->SetOrigin(fileformatInput);
                mol.SetData(ESP);
            }
            else
            {
                cout << "Read " << ESP_q.size() << " ESP charges for " << mol.NumAtoms() << " atoms\n";
            }
          }
        else if (strstr(buffer, "Natural Population") != nullptr)
          {
            hasPartialCharges = true;
            chargeModel = "NBO";
            ifs.getline(buffer,BUFF_SIZE);	// column headings
            ifs.getline(buffer,BUFF_SIZE);  // again
            ifs.getline(buffer,BUFF_SIZE);  // again (-----)
            ifs.getline(buffer,BUFF_SIZE); // real data
            tokenize(vs,buffer);
            while (vs.size() >= 3 &&
                   strstr(buffer, "=====") == nullptr)
              {
                atom = mol.GetAtom(atoi(vs[1].c_str()));
                if (!atom)
                  break;
                atom->SetPartialCharge(atof(vs[2].c_str()));

                if (!ifs.getline(buffer,BUFF_SIZE)) break;
                tokenize(vs,buffer);
              }
          }
        else if(strstr(buffer, " Frequencies -- ")) //vibrational frequencies
        {
          //The info should appear only once as several blocks starting with this line
          tokenize(vs, buffer);
          for(unsigned int i=2; i<vs.size(); ++i)
            Frequencies.push_back(atof(vs[i].c_str()));
          ifs.getline(buffer,BUFF_SIZE); //Red. masses
          ifs.getline(buffer,BUFF_SIZE); //Frc consts
          ifs.getline(buffer,BUFF_SIZE); //IR Inten
          tokenize(vs, buffer);
          for(unsigned int i=3; i<vs.size(); ++i)
            Intensities.push_back(atof(vs[i].c_str()));

          ifs.getline(buffer, BUFF_SIZE); // column labels or Raman intensity
          if(strstr(buffer, "Raman Activ")) {
            tokenize(vs, buffer);
            for(unsigned int i=3; i<vs.size(); ++i)
              RamanActivities.push_back(atof(vs[i].c_str()));
            ifs.getline(buffer, BUFF_SIZE); // Depolar (P)

            while (strstr(buffer, "Atom") == nullptr)
              ifs.getline(buffer, BUFF_SIZE); // eventually column labels
          }
          ifs.getline(buffer, BUFF_SIZE); // actual displacement data
          tokenize(vs, buffer);
          vector<vector3> vib1, vib2, vib3;
          double x, y, z;
          while(vs.size() >= 5) {
            for (unsigned int i = 2; i < vs.size()-2; i += 3) {
              x = atof(vs[i].c_str());
              y = atof(vs[i+1].c_str());
              z = atof(vs[i+2].c_str());

              if (i == 2)
                vib1.push_back(vector3(x, y, z));
              else if (i == 5)
                vib2.push_back(vector3(x, y, z));
              else if (i == 8)
                vib3.push_back(vector3(x, y, z));
            }

            if (!ifs.getline(buffer, BUFF_SIZE))
              break;
            tokenize(vs,buffer);
          }
          Lx.push_back(vib1);
          if (vib2.size())
            Lx.push_back(vib2);
          if (vib3.size())
            Lx.push_back(vib3);
        }

        else if(strstr(buffer, " This molecule is "))//rotational data
        {
          if(strstr(buffer, "asymmetric"))
            RotorType = OBRotationData::ASYMMETRIC;
          else if(strstr(buffer, "symmetric"))
            RotorType = OBRotationData::SYMMETRIC;
          else if(strstr(buffer, "linear"))
            RotorType = OBRotationData::LINEAR;
          else
             RotorType = OBRotationData::UNKNOWN;
          ifs.getline(buffer,BUFF_SIZE); //symmetry number
          tokenize(vs, buffer);
          RotSymNum = atoi(vs[3].c_str());
        }

        else if(strstr(buffer, "Rotational constant"))
        {
          tokenize(vs, buffer);
          RotConsts.clear();
          for (unsigned int i=3; i<vs.size(); ++i)
            RotConsts.push_back(atof(vs[i].c_str()));
        }

        else if(strstr(buffer, "alpha electrons")) // # of electrons / orbital
        {
          tokenize(vs, buffer);
          if (vs.size() == 6) {
            // # alpha electrons # beta electrons
            aHOMO = atoi(vs[0].c_str());
            bHOMO = atoi(vs[3].c_str());
          }
        }
        else if(strstr(buffer, "rbital symmetries")) // orbital symmetries
          {
            symmetries.clear();
            std::string label; // used as a temporary to remove "(" and ")" from labels
            int iii,offset = 0;
            bool bDoneSymm;

            // Extract both Alpha and Beta symmetries
            ifs.getline(buffer, BUFF_SIZE); // skip the current line
            for(iii=0; (iii<2); iii++) {
              if (strstr(buffer, "electronic state"))
                break; // We've gone too far!
              while (!ifs.eof() &&
                     (nullptr != strstr(buffer,"Alpha") ||
                      nullptr != strstr(buffer,"Beta"))) {
                // skip the Alpha: and Beta: title lines
                ifs.getline(buffer, BUFF_SIZE);
              }
              do {
                bDoneSymm = nullptr == strstr(buffer, "(");
                if (!bDoneSymm) {
                  tokenize(vs, buffer);

                  if (nullptr != strstr(buffer, "Occupied") || nullptr != strstr(buffer, "Virtual")) {
                    offset = 1; // skip first token
                  } else {
                    offset = 0;
                  }
                  for (unsigned int i = offset; i < vs.size(); ++i) {
                    label = vs[i].substr(1, vs[i].length() - 2);
                    symmetries.push_back(label);
                  }
                  ifs.getline(buffer, BUFF_SIZE); // get a new line if we've been reading symmetries
                }
                // don't read a new line if we're done with symmetries
              } while (!ifs.eof() && !bDoneSymm);
            } // end alpha/beta section
          }
        else if (strstr(buffer, "Alpha") && strstr(buffer, ". eigenvalues --")) {
          orbitals.clear();
          betaStart = 0;
          while (strstr(buffer, ". eigenvalues --")) {
            tokenize(vs, buffer);
            if (vs.size() < 4)
              break;
            if (vs[0].find("Beta") !=string::npos && betaStart == 0) // mark where we switch from alpha to beta
              betaStart = orbitals.size();
            for (unsigned int i = 4; i < vs.size(); ++i) {
              orbitals.push_back(atof(vs[i].c_str()));
            }
            ifs.getline(buffer, BUFF_SIZE);
          }
        }
        else if(strstr(buffer, " Excited State")) // Force and wavelength data
        {
          // The above line appears for each state, so just append the info to the vectors
          tokenize(vs, buffer);
          if (vs.size() >= 9) {
            double wavelength = atof(vs[6].c_str());
            double force = atof(vs[8].substr(2).c_str()); // remove the "f=" part
            Forces.push_back(force);
            Wavelengths.push_back(wavelength);
          }
        }
        else if(strstr(buffer, " Ground to excited state Transition electric dipole moments (Au):"))
          // Electronic dipole moments
        {
          ifs.getline(buffer, BUFF_SIZE); // Headings
          ifs.getline(buffer, BUFF_SIZE); // First entry
          tokenize(vs, buffer);
          while (vs.size() == 5) {
            double s = atof(vs[4].c_str());
            EDipole.push_back(s);
            ifs.getline(buffer, BUFF_SIZE);
            tokenize(vs, buffer);
          }
        }
        else if(strstr(buffer, "       state          X           Y           Z     R(velocity)")) {
          // Rotatory Strengths
          ifs.getline(buffer, BUFF_SIZE); // First entry
          tokenize(vs, buffer);
          while (vs.size() == 5) {
            double s = atof(vs[4].c_str());
            RotatoryStrengthsVelocity.push_back(s);
            ifs.getline(buffer, BUFF_SIZE);
            tokenize(vs, buffer);
          }
        }
        else if(strstr(buffer, "       state          X           Y           Z     R(length)")) {
          // Rotatory Strengths
          ifs.getline(buffer, BUFF_SIZE); // First entry
          tokenize(vs, buffer);
          while (vs.size() == 5) {
            double s = atof(vs[4].c_str());
            RotatoryStrengthsLength.push_back(s);
            ifs.getline(buffer, BUFF_SIZE);
            tokenize(vs, buffer);
          }
        }

        else if (strstr(buffer, "Forces (Hartrees/Bohr)"))
          {
            ifs.getline(buffer, BUFF_SIZE); // column headers
            ifs.getline(buffer, BUFF_SIZE); // ------
            ifs.getline(buffer, BUFF_SIZE); // real data
          }

        else if (strstr(buffer, "Isotropic = ")) // NMR shifts
          {
            tokenize(vs, buffer);
            if (vs.size() >= 4)
              {
                atom = mol.GetAtom(atoi(vs[0].c_str()));
                OBPairData *nmrShift = new OBPairData();
                nmrShift->SetAttribute("NMR Isotropic Shift");

                string shift = vs[4].c_str();
                nmrShift->SetValue(shift);

                atom->SetData(nmrShift);
              }
          }
        else if (strstr(buffer, "SCF Done:") != nullptr)
          {
            tokenize(vs,buffer);
            mol.SetEnergy(atof(vs[4].c_str()) * HARTEE_TO_KCALPERMOL);
            confEnergies.push_back(mol.GetEnergy());
          }
/* Temporarily commented out until the handling of energy in OBMol is sorted out
        // MP2 energies also use a different syntax

        // PM3 energies use a different syntax
        else if(strstr(buffer,"E (Thermal)") != nullptr)
          {
            ifs.getline(buffer,BUFF_SIZE); //Headers
            ifs.getline(buffer,BUFF_SIZE); //Total energy; what we want
            tokenize(vs,buffer);
            mol.SetEnergy(atof(vs[1].c_str()));
            confEnergies.push_back(mol.GetEnergy());
            }
*/
        else if (strstr(buffer, "Standard basis:") != nullptr)
          {
            add_unique_pairdata_to_mol(&mol,"basis",buffer,2);
          }
        else if (strstr(buffer, "Zero-point correction=") != nullptr)
          {
            tokenize(vs,buffer);
            ezpe = atof(vs[2].c_str());
            ezpe_set = true;
          }
        else if (strstr(buffer, "Thermal correction to Enthalpy=") != nullptr)
          {
            tokenize(vs,buffer);
            Hcorr = atof(vs[4].c_str());
            Hcorr_set = true;
          }
        else if (strstr(buffer, "Thermal correction to Gibbs Free Energy=") != nullptr)
          {
            tokenize(vs,buffer);
            Gcorr = atof(vs[6].c_str());
            Gcorr_set = true;
          }
        else if (strstr(buffer, "CV") != nullptr)
          {
              ifs.getline(buffer,BUFF_SIZE); //Headers
              ifs.getline(buffer,BUFF_SIZE); //Total heat capacity
              tokenize(vs,buffer);
              if (vs.size() == 4)
              {
                  if (vs[0].compare("Total") == 0)
                  {
                      CV = atof(vs[2].c_str());
                      CV_set = true;
                  }
              }
              ifs.getline(buffer,BUFF_SIZE); //Electronic
              ifs.getline(buffer,BUFF_SIZE); //Translational
              tokenize(vs,buffer);
              if ((vs.size() == 4) && (vs[0].compare("Translational") == 0) )
              {
                  Scomponents.push_back(atof(vs[3].c_str()));
              }
              ifs.getline(buffer,BUFF_SIZE); //Rotational
              tokenize(vs,buffer);
              if ((vs.size() == 4) && (vs[0].compare("Rotational") == 0))
              {
                  Scomponents.push_back(atof(vs[3].c_str()));
              }
              ifs.getline(buffer,BUFF_SIZE); //Vibrational
              tokenize(vs,buffer);
              if ((vs.size() == 4) && (vs[0].compare("Vibrational") == 0))
              {
                  Scomponents.push_back(atof(vs[3].c_str()));
              }
          }
        else if (strstr(buffer,"Temperature=") != nullptr &&
                 strstr(buffer,"Pressure=") != nullptr)
          {
              tokenize(vs,buffer);
              temperature = atof(vs[1].c_str());
          }
        else if (strstr(buffer, "(0 K)") != nullptr)
          {
            /* This must be the last else */
            int i,nsearch;
            const char *search[] = { "CBS-QB3 (0 K)", "G2(0 K)", "G3(0 K)", "G4(0 K)", "W1BD (0 K)", "W1U  (0 K)" };
            const char *mymeth[] = { "CBS-QB3", "G2", "G3", "G4", "W1BD", "W1U" };
            const int myindex[] = { 3, 2, 2, 2, 3, 3 };

            nsearch = sizeof(search)/sizeof(search[0]);
            for(i=0; (i<nsearch); i++)
              {
                if (strstr(buffer, search[i]) != nullptr)
                  {
                    tokenize(vs,buffer);
                    E0 = atof(vs[myindex[i]].c_str());
                    E0_set = 1;
                    thermo_method = mymeth[i];
                    break;
                  }
              }
          }
      } // end while

    if (mol.NumAtoms() == 0) { // e.g., if we're at the end of a file PR#1737209
      mol.EndModify();
      return false;
    }

    mol.EndModify();

    // Check whether we have data to extract heat of formation.
    if (ezpe_set && Hcorr_set && Gcorr_set && E0_set &&
        CV_set && (thermo_method.size() > 0))
      {
          extract_thermo(&mol,thermo_method,temperature,ezpe,
                         Hcorr,Gcorr,E0,CV,RotSymNum,Scomponents);
      }

    // Attach orbital data, if there is any
    if (orbitals.size() > 0)
      {
        OBOrbitalData *od = new OBOrbitalData;
        if (aHOMO == bHOMO) {
          od->LoadClosedShellOrbitals(orbitals, symmetries, aHOMO);
        } else {
          // we have to separate the alpha and beta vectors
          std::vector<double>      betaOrbitals;
          std::vector<std::string> betaSymmetries;
          unsigned int initialSize = orbitals.size();
          unsigned int symmSize = symmetries.size();
          if (initialSize != symmSize || betaStart == -1)
            {
              cerr << "Inconsistency: orbitals have " << initialSize << " elements while symmetries have " << symmSize << endl;
            }
          else
            {
              for (unsigned int i = betaStart; i < initialSize; ++i) {
                betaOrbitals.push_back(orbitals[i]);
                if (symmetries.size() > 0)
                  betaSymmetries.push_back(symmetries[i]);
              }
              // ok, now erase the end elements of orbitals and symmetries
              for (unsigned int i = betaStart; i < initialSize; ++i) {
                orbitals.pop_back();
                if (symmetries.size() > 0)
                  symmetries.pop_back();
              }
              // and load the alphas and betas
              od->LoadAlphaOrbitals(orbitals, symmetries, aHOMO);
              od->LoadBetaOrbitals(betaOrbitals, betaSymmetries, bHOMO);
            }
        }
        od->SetOrigin(fileformatInput);
        mol.SetData(od);
      }

    //Attach vibrational data, if there is any, to molecule
    if(Frequencies.size()>0)
    {
      OBVibrationData* vd = new OBVibrationData;
      if (RamanActivities.size() != 0) {
        // check to see if they're all zero
        bool allZero = true;
        for (auto &i : RamanActivities) {
          if (i != 0.0) {
            allZero = false;
            break;
          }
        }
        if (!allZero) {
          vd->SetData(Lx, Frequencies, Intensities, RamanActivities);
        } else { // zero Raman
          vd->SetData(Lx, Frequencies, Intensities);
        }
      } else // no Raman
        vd->SetData(Lx, Frequencies, Intensities);
      vd->SetOrigin(fileformatInput);
      mol.SetData(vd);
    }
    //Attach rotational data, if there is any, to molecule
    if(RotConsts[0]!=0.0)
    {
      OBRotationData* rd = new OBRotationData;
      rd->SetData(RotorType, RotConsts, RotSymNum);
      rd->SetOrigin(fileformatInput);
      mol.SetData(rd);
    }
    // Attach unit cell translation vectors if found
    if (numTranslationVectors > 0) {
      OBUnitCell* uc = new OBUnitCell;
      uc->SetData(translationVectors[0], translationVectors[1], translationVectors[2]);
      uc->SetOrigin(fileformatInput);
      mol.SetData(uc);
    }
    //Attach electronic transition data, if there is any, to molecule
    if(Forces.size() > 0 && Forces.size() == Wavelengths.size())
    {
      OBElectronicTransitionData* etd = new OBElectronicTransitionData;
      etd->SetData(Wavelengths, Forces);
      if (EDipole.size() == Forces.size())
        etd->SetEDipole(EDipole);
      if (RotatoryStrengthsLength.size() == Forces.size())
        etd->SetRotatoryStrengthsLength(RotatoryStrengthsLength);
      if (RotatoryStrengthsVelocity.size() == Forces.size())
        etd->SetRotatoryStrengthsVelocity(RotatoryStrengthsVelocity);
      etd->SetOrigin(fileformatInput);
      mol.SetData(etd);
    }

    // set some default coordinates
    // ConnectTheDots will remove conformers, so we add those later
    mol.SetCoordinates(vconf[vconf.size() - 1]);

    if (!pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.ConnectTheDots();
    if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.PerceiveBondOrders();

    // Set conformers to all coordinates we adopted
    // but remove last geometry -- it's a duplicate
    if (vconf.size() > 1)
      vconf.pop_back();

    mol.SetConformers(vconf);
    mol.SetConformer(mol.NumConformers() - 1);
    // Copy the conformer data too
    confData->SetDimension(confDimensions);
    confData->SetEnergies(confEnergies);
    confData->SetForces(confForces);
    mol.SetData(confData);

    if (hasPartialCharges) {
      mol.SetPartialChargesPerceived();

      // Annotate that partial charges come from Mulliken
      OBPairData *dp = new OBPairData;
      dp->SetAttribute("PartialCharges");
      dp->SetValue(chargeModel); // Mulliken, ESP, etc.
      dp->SetOrigin(fileformatInput);
      mol.SetData(dp);
    }
    mol.AssignTotalChargeToAtoms(total_charge);
    mol.SetTotalCharge(total_charge);
    mol.SetTotalSpinMultiplicity(spin_multiplicity);

    mol.SetTitle(title);
    return(true);
  }